emp03.gms : Check correctness of NLP->MCP reform in EMP

Description

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


Small Model of Type : GAMS


Category : GAMS Test library


Main file : emp03.gms

$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