Hi All, I'm trying to do the following constrained optimization example. Maximize x1*(1-x1) + x2*(1-x2) + x3*(1-x3) s.t. x1 + x2 + x3 = 1 x1 >= 0 and x1 <= 1 x2 >= 0 and x2 <= 1 x3 >= 0 and x3 <= 1 which are the constraints. I'm expecting the answer x1=x2=x3 = 1/3. I tried the "constrOptim" function in R and I'm running into some issues. I first start off by setting my equalities to inequalities x1+x2+x3 >= 1 x1+x2+x3 <= 1.001 However, I get a "Error in constrOptim(c(1.00004, 0, 0), fr, grr, ui = t(A), ci = b) : initial value not feasible Execution halted" error. The values 1.00004,0,0 are well within the constraints domain. fr = function(x){ x1 = x[1] x2 = x[2] x3 = x[3] x1*(1-x1) + x2*(1-x2) + x3*(1-x3) } grr = function(x){ x1 = x[1] x2 = x[2] x3 = x[3] c(1-2*x1,1-2*x2,1-2*x3) } A = matrix(c(1,1,1,-1,-1,-1,1,0,0,0,1,0,0,0,1,1,0,0,0,1,0,0,0,1),3,8,byrow=T) b = c(1,-1.001,0,0,0,1,1,1) y = constrOptim(c(1.00004,0,0),fr,grr,ui = t(A),ci = b) Any help/pointers greatly appreciated. Thanks [[alternative HTML version deleted]]
Your initial value is indeed "infeasible", as the error message says. Your x[1] is 1.00004, which is not in the interval [0,1]. To incorporate both equalities and inequalities (linear or nonlinear), you could try my function `constrOptim.nl'. Contact me off the list if you are interested. Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvaradhan at jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml ---------------------------------------------------------------------------- -------- -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of s t Sent: Wednesday, November 04, 2009 5:48 PM To: r-help at r-project.org Subject: [R] Constrained Optimization Hi All, I'm trying to do the following constrained optimization example. Maximize x1*(1-x1) + x2*(1-x2) + x3*(1-x3) s.t. x1 + x2 + x3 = 1 x1 >= 0 and x1 <= 1 x2 >= 0 and x2 <= 1 x3 >= 0 and x3 <= 1 which are the constraints. I'm expecting the answer x1=x2=x3 = 1/3. I tried the "constrOptim" function in R and I'm running into some issues. I first start off by setting my equalities to inequalities x1+x2+x3 >= 1 x1+x2+x3 <= 1.001 However, I get a "Error in constrOptim(c(1.00004, 0, 0), fr, grr, ui = t(A), ci = b) : initial value not feasible Execution halted" error. The values 1.00004,0,0 are well within the constraints domain. fr = function(x){ x1 = x[1] x2 = x[2] x3 = x[3] x1*(1-x1) + x2*(1-x2) + x3*(1-x3) } grr = function(x){ x1 = x[1] x2 = x[2] x3 = x[3] c(1-2*x1,1-2*x2,1-2*x3) } A matrix(c(1,1,1,-1,-1,-1,1,0,0,0,1,0,0,0,1,1,0,0,0,1,0,0,0,1),3,8,byrow=T) b = c(1,-1.001,0,0,0,1,1,1) y = constrOptim(c(1.00004,0,0),fr,grr,ui = t(A),ci = b) Any help/pointers greatly appreciated. Thanks [[alternative HTML version deleted]]
Hi, This probably does not answer your question, which I presume is about the workings of constrOptim function, but I have a couple of comments and different solutions of your problem. 1. In general, in the problem of this type, one can incorporate the equality constraint(s) into the objective function by adding a penalty term for each constraint. A pretty crude solution would be to redefine your objective function as, say, function(x){ x1 = x[1] x2 = x[2] x3 = x[3] -(x1*(1-x1) + x2*(1-x2) + x3*(1-x3)) + 1000*(x1+x2+x3-1)^2 } and solve a minimization problem with simple 0/1 box constraints. Notice the minus in front of your objective function to turn the problem into the minimization one. 2. Second possibility would be to formulate the Lagrangian using just the equality constraint. Then we can analytically differentiate the Lagrangian WRT x1, x2, x3, and the multiplier and equate the derivatives to zero. This will form a set of, generally nonlinear, equations to solve. Again, one can use a minimization with box constraints to minimize the sum of squared left-hand-sides of the equations. This is probably not the most efficient method of solving a system of equations but it often works and is easy to program. Here is an example of the objective function in this case function(x){ x1 = x[1] x2 = x[2] x3 = x[3] x4 = x[4] (1-2*x1+x4)^2 + (1-2*x2+x4)^2 + (1-2*x3+x4)^2 + (x1+x2+x3-1)^2 } 3. This is just a follow-up to 2. In your case, the partial derivatives of the Lagrangian and the equality constraint are linear. So if this was not a toy problem but a real one, one can just solve a system of linear equations and check that the solution is inside the box constraints. In your case> aa[,1] [,2] [,3] [,4] [1,] -2 0 0 1 [2,] 0 -2 0 1 [3,] 0 0 -2 1 [4,] 1 1 1 0> bb[1] -1 -1 -1 1> solve(aa, bb)[1] 0.3333333 0.3333333 0.3333333 -0.3333333 Hope this helps, Andy __________________________________ Andy Jaworski 518-1-01 Process Laboratory 3M Corporate Research Laboratory ----- E-mail: apjaworski@mmm.com Tel: (651) 733-6092 Fax: (651) 736-3122 From: s t <thampy77@yahoo.com> To: r-help@r-project.org Date: 11/04/2009 04:48 PM Subject: [R] Constrained Optimization Sent by: <r-help-bounces@r-project.org> Hi All, I'm trying to do the following constrained optimization example. Maximize x1*(1-x1) + x2*(1-x2) + x3*(1-x3) s.t. x1 + x2 + x3 = 1 x1 >= 0 and x1 <= 1 x2 >= 0 and x2 <= 1 x3 >= 0 and x3 <= 1 which are the constraints. I'm expecting the answer x1=x2=x3 = 1/3. I tried the "constrOptim" function in R and I'm running into some issues. I first start off by setting my equalities to inequalities x1+x2+x3 >= 1 x1+x2+x3 <= 1.001 However, I get a "Error in constrOptim(c(1.00004, 0, 0), fr, grr, ui = t(A), ci = b) : initial value not feasible Execution halted" error. The values 1.00004,0,0 are well within the constraints domain. fr = function(x){ x1 = x[1] x2 = x[2] x3 = x[3] x1*(1-x1) + x2*(1-x2) + x3*(1-x3) } grr = function(x){ x1 = x[1] x2 = x[2] x3 = x[3] c(1-2*x1,1-2*x2,1-2*x3) } A = matrix(c(1,1,1,-1,-1,-1,1,0,0,0,1,0,0,0,1,1,0,0,0,1,0,0,0,1),3,8,byrow=T) b = c(1,-1.001,0,0,0,1,1,1) y = constrOptim(c(1.00004,0,0),fr,grr,ui = t(A),ci = b) Any help/pointers greatly appreciated. Thanks [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]]
On Wed, 4 Nov 2009 14:48:08 -0800 (PST) s t <thampy77 at yahoo.com> wrote:> I'm trying to do the following constrained optimization example. > Maximize x1*(1-x1) + x2*(1-x2) + x3*(1-x3) > s.t. x1 + x2 + x3 = 1 > x1 >= 0 and x1 <= 1 > x2 >= 0 and x2 <= 1 > x3 >= 0 and x3 <= 1 > which are the constraints. > I'm expecting the answer x1=x2=x3 = 1/3.This is a quadratic programming problem (mininimising/maximising a quadratic function under linear equality and inequality constraints). There are several R packages that solve such problems, e.g.: R> library(quadprog) R> Dmat <- diag(3) R> dvec <- rep(1,3)/3 R> Amat <- cbind(rep(1,3), diag(3), -diag(3)) R> bvec <- c(1, rep(0,3), rep(-1,3)) R> solve.QP(Dmat, dvec, Amat, bvec, meq=1) $solution [1] 0.3333333 0.3333333 0.3333333 $value [1] -0.1666667 $unconstrainted.solution [1] 0.3333333 0.3333333 0.3333333 $iterations [1] 2 0 $iact [1] 1 You may want to consult: http://cran.r-project.org/web/views/Optimization.html HTH. Cheers, Berwin ========================== Full address ===========================Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Maths and Stats (M019) +61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009 e-mail: berwin at maths.uwa.edu.au Australia http://www.maths.uwa.edu.au/~berwin