$title Test modsolstat & solution correctness - multiple QCons & binaries (MIQCP03,SEQ=603)

$onText
Modified from QCP04 to include binary variables

Contributor: Toni Lastusilta
$offText

$if not set DEMOSIZE      $set DEMOSIZE   0
$if not set GLOBALSIZE    $set GLOBALSIZE 0
$if not %DEMOSIZE% == 0   $set DEMOSIZE   1
$if not %GLOBALSIZE% == 0 $set GLOBALSIZE 1

* set default N based on license and type of solver
$if %DEMOSIZE%%GLOBALSIZE% == 00 $set NN  25
$if %DEMOSIZE%%GLOBALSIZE% == 01 $set NN   6
$if %DEMOSIZE%%GLOBALSIZE% == 10 $set NN  14
$if %DEMOSIZE%%GLOBALSIZE% == 11 $set NN   3

$if not set N       $set N %NN%
$if not set MTYPE   $set MTYPE miqcp
$if not set TESTTOL $set TESTTOL 1e-5
scalar mchecks / 0 /;
$if not set QCPMCHECKS $goTo qpmchecks
$if not %QCPMCHECKS% == 0 mchecks = 1;
$goTo donemcheck
$label qpmchecks
$if not  %QPMCHECKS% == 0 mchecks = 1;
$label donemcheck

$eolCom //

set i /1*%N%/;
variables x(i), y(i), z;
binary variable b(i);
scalar bigM /100/;
parameter xp(i), yp(i);

xp(i) = uniform(0,1);
yp(i) = uniform(0,1);
* set some inactive bounds to keep the global solvers happy
x.lo(i) = -1.1;  x.up(i) = 1.1;
y.lo(i) = -1.1;  y.up(i) = 1.1;

equation defc, defz, afew;

defc(i).. sqr(x(i)) + sqr(y(i)) =L= 1 + bigM*b(i);
defz..    z =E= sum(i,xp(i)*x(i) + yp(i)*y(i));
afew..    sum(i, b(i)) =l= max(2,card(i)*0.01);

model m /all/;
m.limrow=0; m.limcol=0;
option optcr = 0;
solve m min z using %MTYPE%;

* disable mchecks if no marginals available
if(not m.marginals, mchecks = 0);

scalars
    vtol    /    1e-8 /,
    tol     /    %TESTTOL% /,
    objval;
parameters
    dLdx(i)  'dLangrangian/dx',
    dLdy(i)  'dLangrangian/dy';

* capability problems is an OK return
if {(m.solvestat = %solveStat.capabilityProblems%),
  abort$(m.modelstat <> %modelStat.noSolutionReturned%)           'wrong modelstat for capability error';
  display 'Solver capability error: further tests suppressed';
else
  abort$(m.modelstat <> %modelStat.optimal% and m.modelstat <> %modelStat.locallyOptimal% and m.modelstat <> %modelStat.feasibleSolution% and m.modelstat <> %modelStat.integerSolution%)  'do not have feasible solution';

* check primal solution feasibility
  objval = sum(i,xp(i)*x.l(i) + yp(i)*y.l(i));
  abort$(abs(z.l-objval) > tol)                       'bad z.l';
  abort$(abs(defz.l) > tol)                           'bad defz.l';
  abort$(smax{i,defc.l(i) - 1} > tol)                 'bad defc.l';

  if {mchecks,
* check dual solution feasibility
    abort$(abs(z.m) > tol)                            'bad z.m';
    abort$(abs(defz.m-1) > tol)                       'bad defz.m';
    abort$(smax{i,defc.m(i)} > tol)                   'bad defc.m';

* check first order optimality condition (assume b fixed)
    dLdx(i) = (xp(i) - defc.m(i) * 2 * x.l(i))$(b.l(i)<0.5);
    dLdy(i) = (yp(i) - defc.m(i) * 2 * y.l(i))$(b.l(i)<0.5);
    display dLdx, dLdy;
    abort$(smax{i,abs(dLdx(i))} > tol)                'bad dLdx';
    abort$(smax{i,abs(dLdy(i))} > tol)                'bad dLdy';
  };
  display 'All tests passed';
};