*==>mypowell.gms $TITLE Powell's nonlinear program $offsymxref offsymlist offuellist offuelxref $inlinecom /* */ /* ****************************************************** Reference: M.J.D. Powell, "A Method for Nonlinear Constraints in Minimization Problems", Optimization (London) (R. Fletcher, ed.) Academic Press, 1969, pp. 238-298. Powell gives the following problem: min f(x) := exp (x1 x2 x3 x4 x5) h1(x) := - 10 = 0 s.t. h2(x) := x2 x3 - 5 x4 x5 = 0 h3(x) := x1**3 + x2**3 + 1 = 0 Powell's reported solution: x = (-1.7172, 1.5957, 1.8272, -.7636, -.7637) (correct to 3 sig figs) By relaxing the first constraint to >= 10, and rewriting the last two as pairs of inequalities, and writing xi as xi+ - xi- (i.e. zi - yi), and taking the KKT conditions, we get the NCP solved by Harker and Xiao. However, the above relaxation of the first constraint does not yield an optimal solution to the original problem: the objective function decreases to zero, as the norm of x increases without bound for a set of otherwise feasible points. We use the set SGN = (pos, neg) to control which sense(s) of the inequality are being used. The variable ui("pos") corresponds to the constraint hi(x) <= 0; E.G. h1("pos") corresponds to the constraint - 10 <= 0; fixing u1("pos") to 0 turns off this constraint. ****************************************************** */ sets run / run1 * run6 /, I / 1 * 5 /, SGN / pos, neg /; alias (J,I); parameters sense(SGN) / pos 1 neg -1 /, includerun (run) / run1 01 run2 01 run3 01 run4 01 run5 01 run6 01 /; table ipt(I,run) /* initial points */ run1 run2 run3 run4 run5 run6 1 -2 -3 -10 -5 -.85 -1.7171 2 2 3 10 5 -.72 1.596 3 2 3 10 5 1.8272 4 -1 -1 -1 -1 -.76364 5 -1 -1 -1 -1 -.76364 ; scalar xnorm / 10 /, objval ; positive variable z(I), /* x = z - y */ y(I), u1(SGN), u2(SGN), u3(SGN); equations delz(I), /* del of Lagrangian wrt z */ dely(I), /* del of Lagrangian wrt y */ negh1(SGN), negh2(SGN), negh3(SGN); delz(I) .. exp(prod(J, z(J)-y(J))) * prod(J$(ord(I) ne ord(J)), z(J)-y(J)) /* delf */ /* u1'delh1 */ + 2 * (u1("pos")-u1("neg")) * (z(I) - y(I)) +( 3*power(z("1")-y("1"),2) * u3("pos")-u3("neg") ) $(ord(I) = 1) +( (z("3")-y("3"))*(u2("pos")-u2("neg")) +3*power(z("2")-y("2"),2)*(u3("pos")-u3("neg")) ) $(ord(I) = 2) +( (z("2")-y("2"))*(u2("pos")-u2("neg")) ) $(ord(I) = 3) +( -5*(z("5")-y("5"))*(u2("pos")-u2("neg")) ) $(ord(I) = 4) +( -5*(z("4")-y("4"))*(u2("pos")-u2("neg")) ) $(ord(I) = 5) =g= 0; dely(I) .. -exp(prod(J, z(J)-y(J))) * prod(J$(ord(I) ne ord(J)), z(J)-y(J)) /* -delf */ /* u1'delh1 */ - 2 * (u1("pos")-u1("neg")) * (z(I) - y(I)) -( 3*power(z("1")-y("1"),2) * u3("pos")-u3("neg") ) $(ord(I) = 1) -( (z("3")-y("3"))*(u2("pos")-u2("neg")) +3*power(z("2")-y("2"),2)*(u3("pos")-u3("neg")) ) $(ord(I) = 2) -( (z("2")-y("2"))*(u2("pos")-u2("neg")) ) $(ord(I) = 3) -( -5*(z("5")-y("5"))*(u2("pos")-u2("neg")) ) $(ord(I) = 4) -( -5*(z("4")-y("4"))*(u2("pos")-u2("neg")) ) $(ord(I) = 5) =g= 0; negh1(SGN) .. 0 =g= sense(SGN) * ( sum (J, power(z(J)-y(J), 2)) - xnorm ); negh2(SGN) .. 0 =g= sense(SGN) * ( (z("2")-y("2"))*(z("3")-y("3")) - 5*(z("4")-y("4"))*(z("5")-y("5")) ); negh3(SGN) .. 0 =g= sense(SGN) * ( power(z("1")-y("1"), 3) + power(z("2")-y("2"), 3) + 1 ); model powell / delz.z, dely.y, negh1.u1, negh2.u2, negh3.u3 /; option limrow = 0; option limcol = 0; option iterlim = 10000; option reslim = 600; option domlim = 100; file out /powell.out/; put out; loop (run, if (includerun(run), z.l(I) = max (0, ipt(I,run)); y.l(I) = max (0, -ipt(I,run)); u1.l(SGN) = 0; u2.l(SGN) = 0; u3.l(SGN) = 0; /* u1.fx("pos") = 0; */ solve powell using mcp; put "Run ", ord(run):2:0 /; put " x"/; put " --------------"/; loop (I, put (z.l(I) - y.l(I)):15:8 /; ); put " u"/; put " --------------"/; put (u1.l("pos") - u1.l("neg")):15:8 /; put (u2.l("pos") - u2.l("neg")):15:8 /; put (u3.l("pos") - u3.l("neg")):15:8 /; objval = exp( prod(I, z.l(I)-y.l(I)) ); put " Objective value: ", objval:15:8 //; ));