Good afternoon,
In the code below, I have a set of functions (m1,m2,m3,s1,s2, and s3) which
represent response surface designs for the mean and variance for three response
variables, followed by an objective function that uses the "Big M"
method to minimize variance (that is, push s1, s2, and s3 as close to 0 as
possible) and hit targets for each of the three means (which are 0, 10, and 100,
respectively). The results are OK but s1 and s3 are negative. I want to
incorporate a constraint or set of constraints that requires s1, s2, and s3 to
be >= 0. I apologize if this is a "dumb" question and the answer
may be staring me in the face, but after several hours of tinkering with min and
max functions within the objective function and performing google searches for
"adding constraints to optimization functions" or the like, I am at a
loss. I am also sure there is a much more elegant way to code what I have done
and so apologize for any crudeness there.
Thank you for any assistance. Code follows:
#Define the Response Surface Designs
m1<-function(x) {
x1<-x[1]
x2<-x[2]
2.1754-0.2219*x1-0.1493*x2-0.1656*x1^2-0.2911*x2^2-0.0862*x1*x2}
m2<-function(x) {
x1<-x[1]
x2<-x[2]
10.0005+0.0465*x1+0.0492*x2-0.0139*x1^2-0.0050*x2^2-0.0325*x1*x2}
m3<-function(x) {
x1<-x[1]
x2<-x[2]
95.1074+0.5288*x1+0.6521*x2-0.1746*x1^2-0.1357*x2^2-0.0712*x1*x2}
s1<-function(x) {
x1<-x[1]
x2<-x[2]
0.0311+0.0000*x1+0.00032*x2-0.01226*x1^2-0.01209*x2^2-0.00075*x1*x2}
s2<-function(x) {
x1<-x[1]
x2<-x[2]
0.003588-0.00022*x1-0.001967*x2+0.001482*x1^2+0.000245*x2^2+0.001375*x1*x2}
s3<-function(x) {
x1<-x[1]
x2<-x[2]
0.17789+0.00683*x1+0.006478*x2-0.07143*x1^2-0.06860*x2^2+0.01338*x1*x2}
# Define the "Big M"
M <- 100000
#Defining the Objective Function
objective1<-function(x) {
x1<-x[1]
x2<-x[2]
M*(s1(c(x1,x2)))+M*(s2(c(x1,x2))) + M*(s3(c(x1,x2))) +
(1/3)*m1(c(x1,x2)) + (1/3)*abs(m2(c(x1,x2))-10) + (1/3)*(100-m3(c(x1,x2)))}
#Optimization
result1 <- nlminb(start=c(-0.3976,1.5541), objective1, gradient = NULL,
hessian = NULL, lower = c(-1.682,-1.682), upper = c(1.682,1.682))
result1$objective
m1(c(result1$par))
m2(c(result1$par))
m3(c(result1$par))
s1(c(result1$par))
s2(c(result1$par))
s3(c(result1$par))
Thanks for any help,
Greg