Description

Test of correctness of the levels and marginals returned.
We have QP terms in the objective only - all QCP solvers accept this.
All cases are considered, e.g.
  1) =L=, =G=, =E= constraints (should we add =N=?)
  2) variables at lower and upper bound, and not at bound
  3) min [convex obj] or max [concave obj]
  4) special attention paid to the form of the obj. constraint, i.e.
     cz * z = xQx + cx + b where cz and b take different values

Small Model of Type : QCP

Category : GAMS Test library

Main file : qcp01.gms

$title Test of correctness for levels & marginals of QCP (QCP01,SEQ=74)
$onText
Test of correctness of the levels and marginals returned.
We have QP terms in the objective only - all QCP solvers accept this.
All cases are considered, e.g.
  1) =L=, =G=, =E= constraints (should we add =N=?)
  2) variables at lower and upper bound, and not at bound
  3) min [convex obj] or max [concave obj]
  4) special attention paid to the form of the obj. constraint, i.e.
     cz * z = xQx + cx + b where cz and b take different values
$offText

$if not set MTYPE   $set MTYPE qcp
$if not set TESTTOL $set TESTTOL 1e-6
$if not set SLOWOK $set SLOWOK 0
scalar slowOK 'slow solves are OK: just abort.noerror in this case' / %SLOWOK% /;
scalar mchecks / 0 /;
$if not %QPMCHECKS% == 0 mchecks = 1;

scalar failed / 0 /;
$escape =
$echo if{%=1, display '%=1 failed', '%=2'; failed=1}; > gtest.gms
$escape %

set J / j1 * j7 /;
set GT / gt1 * gt2 /;
set LT / lt1 * lt2 /;
set EQ / eq1 /;

scalars
    cz,
    cb;

parameters
    c(J) / j1   0
           j2   0
           j3   1
           j4   2
           j5   -1
           j6   -1
           j7   3 /
    bgt(GT) /
           gt1  2
           gt2  4 /,
    blt(LT) /
           lt1  -17
           lt2  10 /,
    beq(EQ) /
           eq1  9 /;

table Agt(GT,J)
        j1      j2
gt1     1       1
gt2     4       1 ;

table Alt(LT,J)
        j3      j4
lt1     -6      1
lt2     1       2 ;

parameter Aeq(EQ,J);
Aeq('eq1',j)$(ord(j) le 7) = 1;


variable
    x(J)
    z;

* these bounds should never be active - but they help the global solvers
x.lo(J) = -2;
x.up(J) =  5;

* these are active and tested to be so
x.lo('j5') = 1;
x.up('j6') = 3;

equations
    gte(GT),
    lte(LT),
    eqe(EQ),
    obj;

obj..   cz*z =E= cb + sum{J,c(J)*x(J)}
         + sqr(x('j1')) + sqr(x('j2')) - x('j1')*x('j2')
         + sqr(x('j3')-x('j4'))
         + sqr(x('j5')+1) + sqr(x('j6')-4)
         + sqr(x('j7')-1);
gte(GT)..  sum{J, Agt(GT,J)*x(J)} =G= bgt(GT);
lte(LT)..  sum{J, Alt(LT,J)*x(J)} =L= blt(LT);
eqe(EQ)..  sum{J, Aeq(EQ,J)*x(J)} =E= beq(EQ);

model m / all /;

set czvals 'obj multipliers' /
'cz=1',
'cz=0.5'
'cz=3'
'cz=-1'
'cz=-0.5'
'cz=-3'
/;
parameter czv(czvals) /
'cz=1'          1
'cz=0.5'        0.5
'cz=3'          3
'cz=-1'         -1
'cz=-0.5'       -0.5
'cz=-3'         -3
/;

set cbvals 'obj constants' /
'cb=0'
'cb=2'
'cb=-2'
/;
parameter cbv(cbvals) /
'cb=0'          0
'cb=2'          2
'cb=-2'         -2
/;

scalars
    isMin   /    0 /,
    tol     /    %TESTTOL% /,
    obj_l   /    0.0 /
    obj_m   /    1.0 /
    objval  /   12.0 /;

parameters
   x_l(J)    / j1       1.0
               j2       1.0
               j3       3.0
               j4       1.0
               j5       1.0
               j6       3.0
               j7       -1.0 /
   x_m(J)    / j1       0.0
               j2       0.0
               j3       0.0
               j4       0.0
               j5       4.0
               j6       -2.0
               j7       0.0 /
   gte_l(GT) / gt1      2.0
               gt2      5.0 /
   gte_m(GT) / gt1      2.0
               gt2      0.0 /
   lte_l(LT) / lt1      -17.0
               lt2      5.0 /
   lte_m(LT) / lt1      -1.0
               lt2      0.0 /
   eqe_l(EQ) / eq1      9.0 /
   eqe_m(EQ) / eq1      -1.0 /;

loop {czvals$[ord(czvals) <= INF],
 cz = czv(czvals);
 loop {cbvals$[ord(cbvals) <= INF],
  cb = cbv(cbvals);
  if {(cz > 0),
    isMin = 1;
    solve m using %MTYPE% min z;
  else
    isMin = -1;
    solve m using %MTYPE% max z;
  };
  display m.solvestat, m.modelstat;
  abort.noError$[slowOK and %solveStat.resourceInterrupt% = m.solvestat] 'Solve too slow';
$ batInclude gtest "( m.solvestat <> %solveStat.normalCompletion% or (m.modelstat > %modelStat.locallyOptimal% and m.modelstat <> %modelStat.feasibleSolution%))" "wrong status codes"

$ batInclude gtest "(abs(cz*z.l-cb-objval) > tol)"             "bad z.l"
$ batInclude gtest "(abs(obj.l-cb) > tol)"                     "bad obj.l"
$ batInclude gtest "(smax(J, abs(x.l(j)-x_l(j))) > tol)"       "bad x.l"
$ batInclude gtest "(smax(GT,abs(gte.l(GT)-gte_l(GT))) > tol)" "bad gte.l"
$ batInclude gtest "(smax(LT,abs(lte.l(LT)-lte_l(LT))) > tol)" "bad lte.l"
$ batInclude gtest "(smax(EQ,abs(eqe.l(EQ)-eqe_l(EQ))) > tol)" "bad eqe.l"

  if {mchecks,
$   batInclude gtest "(abs(z.m) > tol)"                             "bad z.m"
$   batInclude gtest "(abs(cz*obj.m-obj_m) > tol)"                  "bad obj.m"
$   batInclude gtest "(smax(J, abs(cz*x.m(j)-x_m(j))) > tol)"       "bad x.m"
$   batInclude gtest "(smax(GT,abs(cz*gte.m(GT)-gte_m(GT))) > tol)" "bad gte.m"
$   batInclude gtest "(smax(LT,abs(cz*lte.m(LT)-lte_m(LT))) > tol)" "bad lte.m"
$   batInclude gtest "(smax(EQ,abs(cz*eqe.m(EQ)-eqe_m(EQ))) > tol)" "bad eqe.m"
  };

$if set doscalar
  if (failed, option %MTYPE%=convert; if {(cz > 0), solve m using %MTYPE% min z; else solve m using %MTYPE% max z;})
  abort$failed 'test failed', cz, cb;

 };
};