$title Check on MIP solution value (MIP02,SEQ=497)

$onText
In this model, we test that the solution returned by the MIP solver is
the one expected.  The test includes a check on the marginal values
returned.  The model used is the one from magic in the GAMS model
library, adjusted to remove degeneracy.

Contributor: Steve Dirkse, July 2010
$offText

 Sets  t  demand blocks / 12pm-6am, 6am-9am, 9am-3pm, 3pm-6pm, 6pm-12pm /
       g  generators    / type-1, type-2, type-3 /

 Parameters dem(t)  demand (1000mw)   / 12pm-6am  15, 6am-9am   30, 9am-3pm   25, 3pm-6pm  40, 6pm-12pm   27 /
            dur(t)  duration (hours)  / 12pm-6am   6, 6am-9am    3, 9am-3pm    6, 3pm-6pm    3, 6pm-12pm   6 /

 Table data(g,*)  generation data

         min-pow  max-pow  cost-min  cost-inc  start    number
*        (1000mw) (1000mw)  (l/h)    (l/h/mw)   (l)     (units)

 type-1    .85      2.0      1000       2.0     2000      12
 type-2   1.25      1.75     2600       1.3     1000      10
 type-3   1.5       4.0      3000       3.0      500       5

 Parameters peak     peak power (1000mw)
            ener(t)  energy demand in load block (1000mwh)
            tener    total energy demanded (1000mwh)
            lf       load factor ;

  peak = smax(t, dem(t));  ener(t) = dur(t)*dem(t);  tener = sum(t, ener(t));  lf = tener/(peak*24);
  display peak, tener, lf, ener;

$eject
 Variables  x(g,t)  generator output (1000mw)
            n(g,t)  number of generators in use
            s(g,t)  number of generators started up
            cost    total operating cost (l)

 Integer Variables n; Positive Variable s;

 Equations pow(t)    demand for power (1000mw)
           res(t)    spinning reserve requirements (1000mw)
           st(g,t)   start-up definition
           minu(g,t) minimum generation level (1000mw)
           maxu(g,t) maximum generation level (1000mw)
           cdef      cost definition (l);

 pow(t)..  sum(g, x(g,t)) =g= dem(t);

 res(t)..  sum(g, data(g,"max-pow")*n(g,t)) =g= 1.15*dem(t);

* we add a small RHS to equation st() to get unique marginals
* without it both s.lo and st.lo are binding when n(t) = n(t--1)
 st(g,t).. s(g,t) =g= n(g,t) - n(g,t--1) + .0001;

 minu(g,t)..  x(g,t) =g= data(g,"min-pow")*n(g,t);

* we add a small RHS to equation st maxu() to get unique marginals
 maxu(g,t)..  x(g,t) =l= data(g,"max-pow")*n(g,t) + 1e-4;

 cdef.. cost =e= sum((g,t), dur(t)*data(g,"cost-min")*n(g,t) + data(g,"start")*s(g,t)
               + 1000*dur(t)*data(g,"cost-inc")*(x(g,t)-data(g,"min-pow")*n(g,t)) );

  n.up(g,t) = data(g,"number");

 Model william / all /;

 william.optcr = 0;
 william.savepoint = 1;
 Solve william minimizing cost using mip;

scalar rc;
abort$[william.solvestat <> %solveStat.normalCompletion% or
       william.modelstat <> %modelStat.optimal%] 'did not even solve';

* compare duals only if the MIP solver reported them
if( not william.marginals,
  execute 'gdxdiff william_p.gdx mip02_solu.gdx eps=1e-5 releps=1e-4 Field=L > %system.nullfile%';
else
  execute 'gdxdiff william_p.gdx mip02_solu.gdx eps=1e-5 releps=1e-4 > %system.nullfile%';
);
rc=errorlevel;
abort$[rc <> 0] 'Solution not as expected';