$title 'Check correctness of NLP->MCP reform in EMP' (EMP03,SEQ=388) $ontext Test conventions for EMP rewriting as an MCP (aka NLPD). model 1 is a max model, model 2 a min. They are identical except for 1 is max z=f(x), the other is min z=-f(x) and of course the equation duals have opposite sign Key observation: the NLPD option dualVar x f essentially equates f.m with variable x, so we don't compute derivatives wrt. x as we would if x were a primal variable. Since this is so, we assume x has the same sign as f.m would in the NLP model, and we need to do a reformulation that preserves the sign of f.m=x. The model shows how we can reformulate a max model and its equivalent min so that we preserve the sign of the equation duals in the MCP version. Contributor: Steve Dirkse, April 2008 $offtext variable z; positive variables x, y; x.lo = 1e-2; equations f1 'objective definition for max model' f2 'objective definition for min model' e 'equality' up 'upper bound' lo 'lower bound' ; f1.. z =E= x + y; f2.. z =E= -x - y; e.. y =E= sqrt(x); up.. exp(x) =L= 20; lo.. y =G= .1 * x; model nlp1 'max of z=f(x)' / f1, e, up, lo /; model nlp2 'min of z=-f(x)' / f2, e, up, lo /; free variables v1 'perp to e in nlp1' v2 'perp to e in nlp2'; positive variables uUp1 'perp to up in nlp1' uLo2 'perp to lo in nlp2'; negative variable uLo1 'perp to lo in nlp1' uUp2 'perp to up in nlp2'; equations dLdx1, dLdy1 dLdx2, dLdy2 eNeg upNeg loNeg ; dLdx1.. -1 - v1*(0.5/sqrt(x)) + uUp1*exp(x) - uLo1*.1 =N= 0; dLdy1.. -1 + v1 + uUp1*0 + uLo1 =N= 0; dLdx2.. -1 + v2*(0.5/sqrt(x)) - uUp2*exp(x) + uLo2*.1 =N= 0; dLdy2.. -1 - v2 - uUp2*0 - uLo2 =N= 0; eNeg.. sqrt(x) =E= y; upNeg.. 20 =G= exp(x); loNeg.. .1 * x =L= y; model mcp1 / dLdx1.x, dLdy1.y, eNeg.v1, upNeg.uUp1, loNeg.uLo1 /; model mcp2 / dLdx2.x, dLdy2.y, e.v2, up.uUp2, lo.uLo2 /; solve nlp1 using nlp max z; * now set all the duals: * the duals for mcp1 use the same sign as nlp1: THAT IS KEY! v1.l = e.m; uUp1.l = up.m; uLo1.l = lo.m; * the duals for mcp2 flip sign v2.l = -e.m; uUp2.l = -up.m; uLo2.l = -lo.m; * and of course the duals for nlp2 flip sign also z.l = -z.l; f2.m = -f1.m; e.m = -e.m; up.m = -up.m; lo.m = -lo.m; * these next three solves should not take any iterations option nlp = minos; solve nlp2 using nlp min z; abort$[nlp2.iterusd > 0] 'nlp2 took some iterations'; option mcp = path; mcp1.iterlim = 0; solve mcp1 using mcp; abort$[mcp1.objval > 1e-4] 'mcp1 was not given a solution'; mcp2.iterlim = 0; solve mcp2 using mcp; abort$[mcp2.objval > 1e-4] 'mcp2 was not given a solution'; nlp1.optfile = 99; nlp2.optfile = 99; $onecho > jams.o99 dict nlpDict.txt fileName genMCP.gms $offecho $echo modeltype mcp > "%emp.info%" solve nlp2 using emp min z; abort$[nlp2.iterusd > 0] 'nlp2 took some iterations'; z.l = -z.l; f1.m = -f2.m; e.m = -e.m; up.m = -up.m; lo.m = -lo.m; solve nlp1 using emp max z; abort$[nlp1.iterusd > 0] 'nlp1 took some iterations'; $if not set DUMPMCP $exit * now we make it really easy to debug this in case the nlp models * do not solve as expected with EMP. The MCP models below were known * to work $onecho > emp03mcp1.gms * written by GAMS/EMP at 04/12/08 08:20:11 * Variables x2,x3,u2,u3,u4; Negative Variables u4; Positive Variables x3; Positive Variables u3; Equations e2,e3,e4,dL_dx2,dL_dx3; e2.. 0 =E= - sqrt(x2) + x3; e3.. 20 =G= exp(x2); e4.. 0 =L= - 0.1*x2 + x3; dL_dx2.. -1 + ( - 0.5/sqrt(x2))*u2 + (exp(x2))*u3 - 0.1*u4 =N= 0; dL_dx3.. -1 + u2 + u4 =N= 0; * set non-default bounds x2.lo = 0.01; * set non-default levels x2.l = 2.99573227355399; x3.l = 1.73081838260229; u2.l = 1; u3.l = 0.0644440342506719; Model m / e2.u2,e3.u3,e4.u4,dL_dx2.x2,dL_dx3.x3 /; m.limrow=0; m.limcol=0; Solve m using MCP; $offecho $onecho > emp03mcp2.gms * written by GAMS/EMP at 04/12/08 08:21:27 * Variables x2,x3,u2,u3,u4; Negative Variables u3; Positive Variables x3; Positive Variables u4; Equations e2,e3,e4,dL_dx2,dL_dx3; e2.. - sqrt(x2) + x3 =E= 0; e3.. exp(x2) =L= 20; e4.. - 0.1*x2 + x3 =G= 0; dL_dx2.. -1 - ( - 0.5/sqrt(x2))*u2 - (exp(x2))*u3 + 0.1*u4 =N= 0; dL_dx3.. -1 - u2 - u4 =N= 0; * set non-default bounds x2.lo = 0.01; * set non-default levels x2.l = 2.99573227355399; x3.l = 1.73081838260229; u2.l = -1; u3.l = -0.0644440342506719; Model m / e2.u2,e3.u3,e4.u4,dL_dx2.x2,dL_dx3.x3 /; m.limrow=0; m.limcol=0; Solve m using MCP; $offecho