example 2D experimental designs using AlgDesign library

basis_exp_dx.R

#
# experimental designs using different basis
#
library(AlgDesign)

# generate a full factorial combination of factor levels 
candidates = gen.factorial(levels=c(51,51), nVars=2, center=TRUE, varNames=c("X1","X2"))
cand.2 = gen.factorial(levels=c(7,7), nVars=2, center=TRUE, varNames=c("X1","X2"))
# add columns for the orthogonal basis 
candidates = transform(candidates, X1.p = poly(candidates$X1, degree=6), X2.p = poly(candidates$X2, degree=6))

# design using monomial basis 
dx.1 = optFederov(~(X1+I(X1^2)+I(X1^3)+I(X1^4)+I(X1^5)+I(X1^6))*(X2+I(X2^2)+I(X2^3)+I(X2^4)+I(X2^5)+I(X2^6)), data=candidates, nTrials=49)
dx.4 = optFederov(~(X1+I(X1^2)+I(X1^3)+I(X1^4)+I(X1^5)+I(X1^6))*(X2+I(X2^2)+I(X2^3)+I(X2^4)+I(X2^5)+I(X2^6)), data=cand.2, nTrials=49)
# design using orthogonal basis from poly()
dx.2 = optFederov(~(X1.p.1+X1.p.2+X1.p.3+X1.p.4+X1.p.5+X1.p.6)*(X2.p.1+X2.p.2+X2.p.3+X2.p.4+X2.p.5+X2.p.6), data=candidates, nTrials=49)
# the monte carlo method is more suitable for high dimension spaces 
dice = data.frame(var=c("X1","X2"), low=c(-25,-25), high=c(25,25), center=c(0,0), nLevels=c(101,101), round=4, factor=FALSE)
dx.3 = optMonteCarlo(~(X1+I(X1^2)+I(X1^3)+I(X1^4)+I(X1^5)+I(X1^6))*(X2+I(X2^2)+I(X2^3)+I(X2^4)+I(X2^5)+I(X2^6)),data=dice,nTrials=49)
# constraint on the design space
dxconstraint = function(xvec){
 X1=xvec[1]
 X2=xvec[2]
 if(sqrt(X1**2 + X2**2) <= 22){# keep the points within a 22 unit radius 
   return(TRUE)
 } else {
   return(FALSE)
 }
}
dx.5 = optMonteCarlo(~(X1+I(X1^2)+I(X1^3)+I(X1^4)+I(X1^5)+I(X1^6))*(X2+I(X2^2)+I(X2^3)+I(X2^4)+I(X2^5)+I(X2^6)),data=dice,nTrials=3*49,constraints=dxconstraint,RandomStart=TRUE,nCand=40*49)

png("2D_dx_points.png")
plot(candidates[dx.1$rows,]$X1, candidates[dx.1$rows,]$X2, xlab="X1",ylab="X2",col="blue")
points(candidates[dx.2$rows,]$X1, candidates[dx.2$rows,]$X2, pch=2, col="red")
points(dx.3$design$X1, dx.3$design$X2, pch=3, col="green")
points(8.33333333*dx.4$design$X1, 8.33333333*dx.4$design$X2, pch=4, col="darkblue")
legend(x=-20, y=-10, legend=c("monomial, optFederov","ortho poly, optFederov","monomial, optMonteCarlo","monomial, equi-distant"), col=c("blue","red","green","darkblue"), pch=c(1,2,3,4))
dev.off()

png("2D_dx_constraint.png")
plot(seq(-25,25), seq(-25,25), type="n", xlab="X1",ylab="X2")
points(dx.5$design$X1, dx.5$design$X2, pch=1, col="blue")
dev.off()