[R] Creating Optimization Constraints

[boylangc at earthlink.net]

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