powercone1.gms : test conic problem using power cone

Description

Test of correctness of the levels and marginals returned.
Although this test should work for all NLP solvers,
it is designed specifically to test a conic solvers, e.g., Mosek.
It also tests whether slight variations in the algebraic formulation
of the power cone work.

Small Model of Type : NLP

Category : GAMS Test library

Main file : powercone1.gms

$title Test of correct handling of power cone in DNLP (POWERCONE1,SEQ=798)

$onText
Test of correctness of the levels and marginals returned.
Although this test should work for all NLP solvers,
it is designed specifically to test a conic solvers, e.g., Mosek.
It also tests whether slight variations in the algebraic formulation
of the power cone work.


The model is taken from (6.8) in
  https://docs.mosek.com/9.0/capi/tutorial-pow-shared.html

Contributor: Stefan
$offText

$if not set TESTTOL $set TESTTOL 1e-3
$if not set MCHECKS $set MCHECKS 1

Nonnegative Variables x0, x1, x2;
Variable x3, x4, x5, z;
Equations e1, e2, e3a, e3b, e3c, e3d, e3e, e3f, objdef;
Scalar x4obj / 1 /;

objdef.. x3 + x4obj * x4 - x0 =E= z;
e1.. x0 + x1 + 0.5*x2 =E= 2;
e2.. x0**0.2 * x1**0.8 =G= abs(x3);
e3a.. x2**0.4 * x5**0.6 =G= abs(x4);
e3b.. -x2**0.4 * x5**0.6 =L= -abs(x4);
e3c.. x2**0.4 * x5**0.6 =G= sqrt(x4**2);
e3d.. -x2**0.4 * x5**0.6 =L= -sqrt(sqr(x4));
e3e.. x2**0.4 * x5**0.6 =G= x4;
e3f.. x5**0.6 * (-x2**0.4) =L= x4;
x5.fx = 1;

Model ma / objdef, e1, e2, e3a /;
Model mb / objdef, e1, e2, e3b /;
Model mc / objdef, e1, e2, e3c /;
Model md / objdef, e1, e2, e3d /;
Model me / objdef, e1, e2, e3e /;
Model mf / objdef, e1, e2, e3f /;

option bratio = 1;

Scalar modelnr;
Scalar solvestat, modelstat, e3l;
for( modelnr = 1 to 6,
  x0.l = 0.05;
  x1.l = 1;
  x2.l = 3;

  if( modelnr = 5, x4.lo = 0 )
  if( modelnr = 6, x4.lo = -inf; x4.up = 0; x4obj = -1 )

  if( modelnr = 1, Solve ma max z using DNLP; solvestat = ma.solvestat; modelstat = ma.modelstat; e3l = e3a.l )
  if( modelnr = 2, Solve mb max z using DNLP; solvestat = mb.solvestat; modelstat = mb.modelstat; e3l = e3b.l )
  if( modelnr = 3, Solve mc max z using DNLP; solvestat = mc.solvestat; modelstat = mc.modelstat; e3l = e3c.l )
  if( modelnr = 4, Solve md max z using DNLP; solvestat = md.solvestat; modelstat = md.modelstat; e3l = e3d.l )
  if( modelnr = 5, Solve me max z using DNLP; solvestat = me.solvestat; modelstat = me.modelstat; e3l = e3e.l )
  if( modelnr = 6, Solve mf max z using DNLP; solvestat = mf.solvestat; modelstat = mf.modelstat; e3l = e3f.l )

  abort$(solvestat <> %solveStat.normalCompletion%) "wrong solver status, expected normal completion", solvestat ;
  abort$(modelstat > %modelStat.locallyOptimal% and modelstat <> %modelStat.feasibleSolution%) "wrong model status, expected at least feasibility", modelstat ;

  abort$(abs(x0.l - 0.06389795) > %TESTTOL%) "wrong x0.l, expected 0.06389795", x0.l ;
  abort$(abs(x1.l - 0.78336975) > %TESTTOL%) "wrong x1.l, expected 0.78336975", x1.l ;
  abort$(abs(x2.l - 2.30546462) > %TESTTOL%) "wrong x2.l, expected 2.30546462", x2.l ;
  abort$(abs(x3.l - 0.47453779) > %TESTTOL%) "wrong x3.l, expected 0.47453779", x3.l ;
  abort$(abs(x4.l - x4obj*1.39670081) > %TESTTOL%) "wrong x4.l, expected x4obj*1.39670081", x4.l ;
  abort$(abs(e1.l - 2) > %TESTTOL%) "wrong e1.l, expected 2", e1.l ;
  abort$(abs(e2.l) > %TESTTOL%) "wrong e2.l, expected zero", e2.l ;
  abort$(abs(e3l) > %TESTTOL%)  "wrong e3?.l, expected zero", e3l ;
  abort$(abs(z.l - 1.8073406571) > %TESTTOL%) "wrong z.l, expected 1.8073406571", z.l ;

  if( %MCHECKS% <> 0,
    abort$(abs(x0.m) > %TESTTOL%) "wrong x0.m, expected zero", x0.m ;
    abort$(abs(x1.m) > %TESTTOL%) "wrong x1.m, expected zero", x1.m ;
    abort$(abs(x2.m) > %TESTTOL%) "wrong x2.m, expected zero", x2.m ;
    abort$(abs(x3.m) > %TESTTOL%) "wrong x3.m, expected zero", x3.m ;
    abort$(abs(x4.m) > %TESTTOL%) "wrong x4.m, expected zero", x4.m ;
    abort$(abs(x5.m - 0.83801326) > %TESTTOL%) "wrong x5.m, expected 0.83801326", x5.m ;
    abort$(abs(e1.m - 0.48466370) > %TESTTOL%) "wrong e1.m, expected 0.48466370", e1.m ;
    abort$(abs(e2.m + 1) > %TESTTOL%) "wrong e2.m, expected -1", e2.m ;
    abort$(modelnr = 1 and abs(e3a.m + 1) > %TESTTOL%) "wrong e3a.m, expected  1", e3a.m ;
    abort$(modelnr = 2 and abs(e3b.m - 1) > %TESTTOL%) "wrong e3b.m, expected -1", e3b.m ;
    abort$(modelnr = 3 and abs(e3c.m + 1) > %TESTTOL%) "wrong e3c.m, expected  1", e3c.m ;
    abort$(modelnr = 4 and abs(e3d.m - 1) > %TESTTOL%) "wrong e3d.m, expected -1", e3d.m ;
    abort$(modelnr = 5 and abs(e3e.m + 1) > %TESTTOL%) "wrong e3e.m, expected  1", e3e.m ;
    abort$(modelnr = 6 and abs(e3f.m - 1) > %TESTTOL%) "wrong e3f.m, expected -1", e3f.m ;
  )
)