$title Test dual solution of SOCP

$onText
Essentially maximizes x+y such that sqr(x) + sqr(y) <= 8.
Thus, x = y = 2.

However, to test conic constraints, we write the constraint
as sqr(x) + sqr(y) <= sqr(z), 0 <= z <= sqrt(8).

Contributor: Stefan Vigerske
$offText

$if not set TESTTOL $set TESTTOL 1e-6

scalar mchecks / 0 /;
$if not %QPMCHECKS% == 0 $if not %QCPMCHECKS% == 0 mchecks = 1;

Variables x, y, z, objvar;
Equation e1, e2, obj;

e1..  sqr(x) + sqr(y) =L=  sqr(z);
e2.. -sqr(x) - sqr(y) =G= -sqr(z);
z.lo = 0;
z.up = sqrt(8);

obj.. objvar =E= x + y;

Model m1 / e1, obj /;
Model m2 / e2, obj /;

* don't less GAMS mess around with the solution
m1.tolproj = 0;
m2.tolproj = 0;

* from global solvers, we want optimal solutions
m1.optcr = 0;
m2.optcr = 0;

* for local solvers, the problem might look nonconvex, so lets start close to optimum
x.l = 1.9;
y.l = 1.9;
z.l = sqrt(sqr(x.l) + sqr(y.l));

Solve m1 max objvar using QCP;

abort$(m1.solvestat <> %solveStat.normalCompletion%)  'solve status is not normal';
abort$(m1.modelstat <> %modelStat.locallyOptimal% and m1.modelstat <> %modelStat.optimal%) 'model status does not indicate (local) optimality';

abort$(abs(x.l - 2) > sqrt(%TESTTOL%))       'wrong primal solution for x, expected 2';
abort$(abs(y.l - 2) > sqrt(%TESTTOL%))       'wrong primal solution for y, expected 2';
abort$(abs(z.l - sqrt(8)) > sqrt(%TESTTOL%)) 'wrong primal value for z, expected sqrt(8)';
abort$(abs(objvar.l - 4) > %TESTTOL%)  'wrong optimal value, expected 4';
abort$(abs(e1.l) > %TESTTOL%)          'wrong activity for e1, expected 0';

if( m1.marginals and mchecks,
  abort$(abs(x.m) > %TESTTOL%)           'wrong dual solution for x, expected 0';
  abort$(abs(y.m) > %TESTTOL%)           'wrong dual solution for y, expected 0';
  abort$(abs(z.m - sqrt(2)) > %TESTTOL%) 'wrong dual solution for z, expected sqrt(2)';
  abort$(abs(e1.m - 0.25) > %TESTTOL%)   'wrong dual solution for e1, expected 0.25';
);

Solve m2 max objvar using QCP;

abort$(m2.solvestat <> %solveStat.normalCompletion%)  'solve status is not normal';
abort$(m2.modelstat <> %modelStat.locallyOptimal% and m2.modelstat <> %modelStat.optimal%) 'model status does not indicate (local) optimality';

abort$(abs(x.l - 2) > sqrt(%TESTTOL%))       'wrong primal solution for x, expected 2';
abort$(abs(y.l - 2) > sqrt(%TESTTOL%))       'wrong primal solution for y, expected 2';
abort$(abs(z.l - sqrt(8)) > sqrt(%TESTTOL%)) 'wrong primal value for z, expected sqrt(8)';
abort$(abs(objvar.l - 4) > %TESTTOL%)  'wrong optimal value, expected 4';
abort$(abs(e2.l) > %TESTTOL%)          'wrong activity for e2, expected 0';

if( m2.marginals and mchecks,
  abort$(abs(x.m) > %TESTTOL%)           'wrong dual solution for x, expected 0';
  abort$(abs(y.m) > %TESTTOL%)           'wrong dual solution for y, expected 0';
  abort$(abs(z.m - sqrt(2)) > %TESTTOL%) 'wrong dual solution for z, expected sqrt(2)';
  abort$(abs(e2.m + 0.25) > %TESTTOL%)   'wrong dual solution for e2, expected -0.25';
);

* do another solve of m1 where the SOCP constraint e1 must not be active
* solution should be x = y = 1 and some z >= 2, thus sqr(x) + sqr(y) = 2 < 4 <= sqr(z)
x.up = 1;
y.up = 1;
z.lo = 2;
Solve m1 max objvar using QCP;

abort$(m1.solvestat <> %solveStat.normalCompletion%)  'solve status is not normal';
abort$(m1.modelstat <> %modelStat.locallyOptimal% and m1.modelstat <> %modelStat.optimal%) 'model status does not indicate (local) optimality';

abort$(abs(x.l - 1) > sqrt(%TESTTOL%))       'wrong primal solution for x, expected 1';
abort$(abs(y.l - 1) > sqrt(%TESTTOL%))       'wrong primal solution for y, expected 1';
abort$(abs(objvar.l - 2) > %TESTTOL%)  'wrong optimal value, expected 2';
abort$(e1.l > -2 + %TESTTOL%)          'wrong activity for e1, expected at most -2';

if( m1.marginals and mchecks,
  abort$(abs(x.m - 1) > %TESTTOL%)  'wrong dual solution for x, expected 1';
  abort$(abs(y.m - 1) > %TESTTOL%)  'wrong dual solution for y, expected 1';
  abort$(abs(z.m) > %TESTTOL%)      'wrong dual solution for z, expected 0';
  abort$(abs(e1.m) > %TESTTOL%)     'wrong dual solution for e1, expected 0';
);