$title Generate Sample
$ontext
    generateSampleGrid.gms
    GAMS code for implementing the bilevel estimation program of chapter 2,
    plus some alternative solution algorithms. Note that storage costs are
    also included in contrast to paper, but with number of periods equal to
    one (the set t contains one element), equivalence is obtained.
    Progress feedback works only when run from DOS prompt. Minimize DOS
    window to improve speed considerably, as screen updating takes time.

    Usage:
    This program can be started in several instances parallel.
    First it should be run in gridmode==0 to generate the models to
    estimate, which are then stored in "models.gdx". Then, the program
    should be started again with gridmode set to 1 and with up to five
    command line arguments
    --gridmode=     "0" to generate true models, "1" to estimate models
                    already generated with gridmode=0
    <below: only with gridmode=0>
    --firstmodel=   number of the first model to estimate
    --lastmodel=    number of the last model to estimate
    --pnr=          number of the process (unique, 1 to 4)
    --o=            name of list file (unique string)

    For example, if you have access to a computer with four processors,
    you would like to run the following commands from four different DOS-
    prompts (called DOS1, DOS2, DOS3, and DOS4):
    DOS1:\>gams generateSampleGrid --gridmode=0
    <wait until process finishes, then run the following four processes from one window each>
    DOS1:\>gams generateSampleGrid --gridmode=1 --firstmodel=1  --lastmodel=25  --pnr=1 o=proc1.lst
    DOS2:\>gams generateSampleGrid --gridmode=1 --firstmodel=26 --lastmodel=50  --pnr=2 o=proc2.lst
    DOS3:\>gams generateSampleGrid --gridmode=1 --firstmodel=51 --lastmodel=75  --pnr=3 o=proc3.lst
    DOS4:\>gams generateSampleGrid --gridmode=1 --firstmodel=76 --lastmodel=100 --pnr=4 o=proc4.lst


    Torbjoern Jansson
    LEI, The Hague, NL
$offtext
$offlisting

* Temporary directory: set to a fast local drive
$setglobal tempdir C:\TEMP

* Enable multiple processes working on partitions of model set
*$setglobal gridmode 1
$ifi %gridmode%==0 $setglobal pnr
$ifi %gridmode%==0 $setglobal firstmodel 1
$ifi %gridmode%==0 $setglobal lastmodel 0

* Create new set of random numbers for each partition
*$ifi %gridmode%==1  execseed=1+gsecond(jnow);
$ifi %gridmode%==1  execseed=%pnr%;

* Write option file for conopt
$onecho > conopt.opt
rtmaxj = 1.00e+9
rtmaxv = 1.00e+9
lfilos =  100
lfilog =  100
$offecho

* Write log file to remember what is being executed
$onecho > process_%pnr%.txt
Process nr:     %pnr%
First model:    %firstmodel%
Last model:     %lastmodel%
Grid mode:      %gridmode%
$offecho


scalar saveGDX "Set to 1 to save results for analyseSample.gms" /1/;

set m Number of models to generate  /m1*m100/;
set n Number of estimation attempts /n1*n1000/;
set mrun(m) Models to analyse;

set a Solution approaches implemented/
  a0  Not a solution method just the raw data saved
  a1  Direct solution bound to fail
  a2  Perpendicular starting point working well in symmetric case
  a3  Global optimisation with Baron as NLP
  a4  Global optimisation with Baron as MINLP
  a5  Simpler but not better own method
  a6  Facchinei et al or NLPEC solver smooth approximation
  a7  Yet another enumeration method not working
  a8  Traditional method not using price observations
  a9  Penalty function iterative approximation as in NLPEC solver
  a10 Best of a6 and a9
  /;

set aon(a) Solution approaches to test /a6,a8,a9,a10/;
alias(aon,aon2);

* Declarations for the economic model (the inner problem)

set i Regions /i1*i10/;
alias(i,j,k,ii,jj);

set t Periods /t1*t1/;
alias(t,t1,t2);

set tnext(t,t1);
  tnext(t,t1) $  (ord(t1) = (ord(t)+1))                = yes;
  tnext(t,t1) $ ((ord(t1) = 1) and (ord(t) = card(t))) = yes;

set im(i) Regions with market balance (one is dropped);
    im(i) = yes$(not sameas(i,"i1"));
*   im(i) = yes;

set AD(i,j) Set of admissible transport flows;
    AD(i,j) $ (not sameas(i,j)) = yes;

parameter c(i,j) Observed (true) transport costs;
parameter so(i)  Observed (true) storage costs;
parameter p(t,i) Observed prices;
parameter e(t,i) Excess demand;

scalar fp Fix price that shifts (the observed) price system /120/;
scalar mu Smoothing parameter of relaxation                 /0/;

scalar w1 Weight for transport cost in estimation /1/;
scalar w2 Weight for storage cost in estimation   /1/;
scalar w3 Weight for prices in estimation         /1/;

* *** Declarations belonging to the models ***

free variables
    pe(t,i)   Estimated price of region i
    z         Free objective variable
    zz        Another free objective variable
    zpen      Penalty function value;

positive variables
    x(t,i,j)  Transport stream from i to j
    st(t,i)   Storage
    ce(i,j)   Estimated transport cost
    se(i)     Estimated storage   cost
    pi(t,i,j) Dual value of lower bound on x
    ro(t,i)   Dual value of lower bound on st;

binary variables
    bt(t,i,j) Is there a trade flow
    bs(t,i)   Is there storage;

equations
    F         Objective function definition
    h(t,i)    Market balance i
    dx(t,i,j) First order conditions for optimal transportation
    ds(t,i)   First order conditions for optimal storage
    cx        Complementarity restriction
    cs
    pen       Penalty function for complementary slackness
    FPen      Objective function penalty term

*     Binary tree for Baron
    bdx(t,i,j)    Alternative foc tp
    bds(t,i)      Alternative foc st
    bcx(t,i,j)    Alternateve cmp tp
    bcs(t,i)      Alternateve cmp st

    DUM           Dummy
    TPC     Transport cost;

*   Weighted least squares deviation

F       ..   z =E=

*                  Squared deviations of costs from observations.
*                  Symmetry ==> only upper triangle less diagonal

            w1*sum(AD(i,j) $ (ord(i) lt ord(j)), sqr(ce(i,j) - c(i,j)))

*                  Squared deviations of prices to observations.

          + w3*sum((t,i), sqr(pe(t,i) - p(t,i)))

*                  Squared deviations of storage costs from observations

          + w2*sum(i, sqr(se(i) - so(i)));

*   Penalty function approach

Pen ..  zpen =e= mu*sum((t,i,j) $ AD(i,j), pi(t,i,j) * x(t,i,j))
               + mu*sum((t,i)            , ro(t,i) * st(t,i));

FPen .. zz =e= z + zpen;

*   Market balance

h(t,im)   ..   sum(AD(im,j), x(t,j,im)-x(t,im,j)) - st(t,im)
             + sum(tnext(t1,t), st(t1,im)) =E= e(t,im);

*   First order condition for transport problem

dx(t,i,j) $ AD(i,j) ..   (ce(i,j) $ (ord(i) lt ord(j))
                         +ce(j,i) $ (ord(i) gt ord(j)))
                       + pe(t,i) - pe(t,j) =E= pi(t,i,j);

ds(t,i) ..  se(i) - sum(tnext(t,t1), pe(t1,i)) + pe(t,i) =E= ro(t,i);

*   Complementary slackness condition for transport problem

cx (t,i,j) $ AD(i,j)  ..  pi(t,i,j) * x(t,i,j) =L= mu;

cs (t,i) ..   ro(t,i) * st(t,i) =L= mu;

*   BARONS binary first order condition for transport problem

bdx(t,i,j) $ AD(i,j)  ..  pi(t,i,j) =l= pi.up(t,i,j) * (1-bt(t,i,j));
bds(t,i)              ..  ro(t,i)   =l= ro.up(t,i)   * (1-bs(t,i));

bcx (t,i,j) $ AD(i,j) ..  x(t,i,j)  =l= x.up(t,i,j)  * bt(t,i,j);
bcs (t,i)             ..  st(t,i)   =l= st.up(t,i)   * bs(t,i);

DUM .. z =l= 10;

*   Transport cost minimisation

TPC .. z =E= sum((t,i,j)$AD(i,j), x(t,i,j)*(ce.l(i,j)$ (ord(i) lt ord(j))
                                           +ce.l(j,i)$ (ord(i) gt ord(j))))
           + sum((t,i), st(t,i)*se.l(i));

model EstimNLP  NLP formulation of MPEC     /F,h,dx,cx,ds,cs/;
model EstimPEN  Penalty formulation of MPEC /F,h,dx,   ds,pen,FPen/;
model TPmin     Transportation model        /TPC,h  /;
model EstimPre  Relaxed version of MPEC     /F,h,dx,ds/;
model BaronNLP  Binary formulation of MPEC  /F,h,dx,bdx,bcx,ds,bds,bcs/;
model DUMMY /DUM/;

EstimNLP.limcol    = 0;
EstimPEN.limcol    = 0;
TPmin.limcol       = 0;
TPmin.limrow       = 0;
TPmin.optfile      = 1;

TPmin.solprint     = 2;
EstimNLP.solprint  = 2;
EstimPEN.solprint  = 2;
EstimPRE.solprint  = 2;

EstimPRE.solvelink = 2;
EstimNLP.solvelink = 2;
EstimPEN.solvelink = 2;
TPmin.solvelink    = 2;

EstimNLP.workspace = 100;
EstimNLP.optfile   = 1;
EstimPEN.workspace = 100;
EstimPEN.optfile   = 1;

DUMMY.solprint = 2;

*   Declarations of items used to save program output

parameter ptru(m,t,i)     True price;
parameter ctru(m,i,j)     True transportation cost;
parameter stru(m,i)       True storage cost;
parameter etru(m,t,i)     True excess demand;

parameter pest(m,n,t,i,a) Estimated price;
parameter cest(m,n,i,j,a) Estimated transportation cost;
parameter sest(m,n,i,a)   Estimated storage cost;

parameter objes(*,m,n);

set enum /e1*e1000/;
set cenum/s,x,pi,rho/;
parameter zenum(m,n,cenum,enum);
scalar nenum /1/;

*   For the log window

scalar logcount 'Counter for iteration log' /0/;
scalar iflog    'Iteration log frequency'   /0/;
scalar progress 'Progress fraction'         /0/;
file batch /%tempdir%\titlebatch%pnr%.bat/;
batch.lw = 0;
iflog = (card(m)*card(n))/1000;

*   Bound prices to help Baron

pe.lo(t,i) = 0;
pe.up(t,i) = fp*2;

$onecho > init.gms
*   Automatically generated include file that initializes problem at a
*   feasible point, using observations.
*   1) Reset all variables with bounds and equations of importance

    option kill=pe;
    option kill=ce;
    option kill=se;
    option kill=x;
    option kill=F;
    option kill=h;
    option kill=dx;
    option kill=ds;
    option kill=cx;
    option kill=cs;
    option kill=ro;
    option kill=pi;
    option kill=z;

*   2) Restore bounds on transport costs

    ce.up(i,j) = 200 $ (AD(i,j) and (ord(i) lt ord(j)));

*   3) Starting point for costs is observation

    ce.l(i,j) $ (AD(i,j) and (ord(i) lt ord(j))) = c(i,j);
    se.l(i) = so(i);

*   4) Find a corresponding feasible price vector by solving TP problem

    solve TPmin using nlp minimising z;

*   pe.l(t,i) = fp + h.m(t,i) $ im(i);
    pe.l(t,i) = p("t1","i1") + h.m(t,i) $ im(i);

*   5) Initialise the first order conditions using duals of TP problem

    pi.l(t,i,j) $ AD(i,j)
        = ce.l(i,j) $ (ord(i) lt ord(j))
        + ce.l(j,i) $ (ord(i) gt ord(j)) + pe.l(t,i) - pe.l(t,j);

    ro.l(t,i) = se.l(i) - sum(tnext(t,t1), pe.l(t1,i)) + pe.l(t,i);
$offecho

$onechov > savesol.gms
* Automatically generated batinclude that saves the outcome of a solution
objes("%1",m,n) = z.l;
pest(m,n,t,i,"%1") = pe.l(t,i);
cest(m,n,i,j,"%1") = ce.l(i,j);
sest(m,n,i,"%1")   = se.l(i);
$offecho

* -------------------------------------------------------------------------
*   Generate m models
* -------------------------------------------------------------------------

scalar gridmode /%gridmode%/;
loop(m $ (not gridmode),

*   START: Generate true model with prices consistent with transport costs

    c(i,j) $ (ord(i) lt ord(j)) = uniform(20,100);
    c(i,j) $ (ord(i) gt ord(j)) = c(j,i);

    so(i) = uniform(3,10);

    e(t,i)    = uniform(-10,10);
*   Make sure solution exists by adjusting excess demand on one point
    e("t1","i1") = -sum((t,i) $ (not (    sameas(t,"t1")
                                      and sameas(i,"i1"))), e(t,i));

    ce.l(i,j) =  c(i,j);
    se.l(i)   = so(i);

    solve TPmin using nlp minimising z;

    p(t,i) = fp + h.m(t,i) $ im(i);

    ptru(m,t,i) =  p(t,i);
    ctru(m,i,j) =  c(i,j);
    stru(m,i)   = so(i);
    etru(m,t,i) =  e(t,i);
);
$ifi %gridmode%==0 execute_unload "models.gdx" m,t,i,j,ptru,ctru,stru,etru;
$ifi %gridmode%==0 $exit

* 1: Process number
* 2: First model to estimate
* 3: Last model to estimate
if(gridmode,
    execute_load "models.gdx" ptru,ctru,stru,etru;
    mrun(m) = yes $ ((ord(m) ge %firstmodel%) and (ord(m) le %lastmodel%));
else
    mrun(m) = yes;
);
display "Running for firstmodel=%firstmodel%, lastmodel=%lastmodel%";
display "Process number: %pnr%", mrun;
*$exit

* -------------------------------------------------------------------------
*       Add error terms to the model n times and estimate m
* -------------------------------------------------------------------------

loop(m $ mrun(m),

    e(t,i)             = etru(m,t,i);

    loop(n,

*       Add error terms to costs and prices, truncate to avoid negatives

        c(i,j) $ (ord(i) lt ord(j)) = ctru(m,i,j)
                                    + min(19, max(-19,normal(0,6)));
        c(i,j) $ (ord(i) gt ord(j)) = c(j,i);

        p(t,i)           = ptru(m,t,i) + min(19, max(-19,normal(0,6)));
        so(i)            = max(stru(m,i)   + normal(0,(3/2)), 0.2);

        pest(m,n,t,i,"a0") =  p(t,i);
        cest(m,n,i,j,"a0") =  c(i,j);
        sest(m,n,i,"a0")   = so(i);

*       Show progress in title bar of DOS window

        logcount=logcount-1;
        if(logcount le 0,
            logcount=iflog;
            progress=100-100*((ord(m)-%firstmodel%)*card(n)+ord(n)-1)/card(n)/card(mrun);
            putclose batch "title M: ",m.tl,", ",progress:0:1,"%% left";
            execute "%tempdir%\titlebatch%pnr%";
        );

*       -------------------- approach 1 -----------------------------------
        if(aon("a1"),

$include "init.gms"

*           a) Solve using NLP from the starting point just obtained

            solve EstimNLP using NLP minimizing z;

$batinclude "savesol.gms" a1
        );
*       -------------------- approach 2 -----------------------------------
        if(aon("a2"),

$include "init.gms"

*        a) Solve the relaxed MPEC

            solve EstimPre using NLP minimizing z;

*        b) Solve TP-problem again

            solve TPmin using nlp minimizing z;

            pe.l(t,i) = fp + h.m(t,i) $ im(i);

            pi.l(t,i,j) $ AD(i,j)
                = ce.l(i,j) $ (ord(i) lt ord(j))
                + ce.l(j,i) $ (ord(i) gt ord(j)) + pe.l(t,i) - pe.l(t,j);

            ro.l(t,i) = se.l(i) - sum(tnext(t,t1), pe.l(t1,i)) + pe.l(t,i);

*        c) Solve the full MPEC from this new point using NLP

            mu = 0;
            solve EstimNLP using NLP minimising z;

$batinclude "savesol.gms" a2
        );
*       -------------------- approach 3 -----------------------------------
        if(aon("a3"),

*           Comment out to let Baron start from previous solution
$include "init.gms"

*           Extra: bounds on all variables, just for the Baron..

            se.up(i)     = 100*smax(t,pe.up(t,i)-pe.lo(t,i));
            pi.up(t,i,j) = ce.up(i,j) + pe.up(t,i)-pe.lo(t,i);
            ro.up(t,i)   = se.up(i)   + pe.up(t,i)-pe.lo(t,i);
            x.up(t,i,j)  = sum((t1,k), ABS(e(t1,k)));
            st.up(t,i)   = sum((t1,j), ABS(e(t1,j)));

*        a) Solve using Baron from the starting point just generated

            option NLP=Baron;
            solve EstimNLP using NLP minimizing z;
            option NLP=Conopt;

$batinclude "savesol.gms" a3
        );

*       -------------------- approach 4 -----------------------------------
        if(aon("a4"),

*           Comment out to let Baron start from previous solution
$include "init.gms"

*           Extra: bounds on all variables, just for the Baron..

            se.up(i)     = 100*smax(t,pe.up(t,i)-pe.lo(t,i));
            pi.up(t,i,j) = ce.up(i,j) + pe.up(t,i)-pe.lo(t,i);
            ro.up(t,i)   = se.up(i)   + pe.up(t,i)-pe.lo(t,i);
            x.up(t,i,j)  = sum((t1,k), ABS(e(t1,k)));
            st.up(t,i)   = sum((t1,j), ABS(e(t1,j)));

*           Initialise binary variables to feasible point

            bt.l(t,i,j) = 1 $  x.l(t,i,j);
            bs.l(t,i)   = 1 $ st.l(t,i);

*        a) Solve using Baron MINLP from the starting point just generated

            option MINLP=Baron;
            BaronNLP.solprint = 2;
            solve BaronNLP using MINLP minimizing z;
            BaronNLP.solprint = 2;

$batinclude "savesol.gms" a4
        );

*       -------------------- approach 5 -----------------------------------
        if(aon("a5"),

$include "init.gms"

*        a) Solve the MPEC using NLP from the starting point just generated

            solve EstimNLP using NLP minimizing z;

$batinclude "savesol.gms" a5
        );
*       -------------------- approach 6 -----------------------------------
        if(aon("a6"),

$include "init.gms"

*        a) Solve problem without complementarity constraints

            solve EstimPre using NLP minimizing z;

*        b) Iterate over increasingly better approximations

            mu = 1;

            while(mu,

*               b1) Make feasible

                solve TPmin using nlp minimizing z;

                pe.l(t,i) = fp + h.m(t,i) $ im(i);

                pi.l(t,i,j) $ AD(i,j)
                    = ce.l(i,j) $ (ord(i) lt ord(j))
                    + ce.l(j,i) $ (ord(i) gt ord(j)) + pe.l(t,i)-pe.l(t,j);

                ro.l(t,i) = se.l(i) -sum(tnext(t,t1),pe.l(t1,i))+pe.l(t,i);

*               b2) Solve approximation

                mu $ (mu<0.001) = 0;

                solve EstimNLP using NLP minimising z;

                mu = mu/2;
            );

$batinclude "savesol.gms" a6
        );
*       -------------------- approach 7 -----------------------------------
        if(aon("a7"),

$include "init.gms"

*           Force one column after another into the basis

            nenum=1;

*        a) Spatial price equilibrium

            loop((t2,ii,jj) $ AD(ii,jj),
                pi.up(t2,ii,jj) = 0;
                solve EstimPRE using nlp minimizing z;
                solve TPmin using nlp minimizing z;
                pe.l(t,i) = fp + h.m(t,i) $ im(i);

*               Initialise the first order condition as well

                pi.l(t,i,j) $ AD(i,j)
                    = ce.l(i,j) $ (ord(i) lt ord(j))
                    + ce.l(j,i) $ (ord(i) gt ord(j)) + pe.l(t,i)-pe.l(t,j);

                ro.l(t,i) = se.l(i) -sum(tnext(t,t1),pe.l(t1,i))+pe.l(t,i);
                solve EstimNLP using nlp minimizing z;
                pi.up(t2,ii,jj) = inf;
                zenum(m,n,"pi",enum) $ (ord(enum)=nenum) = z.l;
                nenum=nenum+1;
            );

*        b) Intertemporal price equilibrium

            loop(ii,
                ro.up(t2,ii) = 0;
                solve EstimPRE using nlp minimizing z;
                solve TPmin using nlp minimizing z;
                pe.l(t,i) = fp + h.m(t,i) $ im(i);

*               Initialise the first order condition as well

                pi.l(t,i,j) $ AD(i,j)
                    = ce.l(i,j) $ (ord(i) lt ord(j))
                    + ce.l(j,i) $ (ord(i) gt ord(j)) + pe.l(t,i)-pe.l(t,j);

                ro.l(t,i) = se.l(i) -sum(tnext(t,t1),pe.l(t1,i))+pe.l(t,i);
                solve EstimNLP using nlp minimizing z;
                ro.up(t2,ii) = inf;
                zenum(m,n,"rho",enum) $ (ord(enum)=nenum) = z.l;
                nenum=nenum+1;
            );

            z.l = smin((cenum,enum)$(    zenum(m,n,cenum,enum)
                                                 and (ord(enum) lt nenum)),
                                    zenum(m,n,cenum,enum));

$batinclude "savesol.gms" a7
        );
*       -------------------- approach 8 -----------------------------------
        if(aon("a8"),

$include "init.gms"

            z.l =
               sum(AD(i,j) $ (ord(i) lt ord(j)), sqr(ce.l(j,i)-c(i,j)))*w1
             + sum((t,i),                        sqr(pe.l(t,i)-p(t,i)))*w3
             + sum(i,                            sqr(se.l(i)  -so(i))) *w2
               ;

$batinclude "savesol.gms" a8
        );

*       -------------------- approach 9 -------------------------------------------
*       A penalty function approach,
*       with ordinary complementary slackness times mu in objective

        if(aon("a9"),

$include "init.gms"

*        b) Solve problem without complementarity constraints

            solve EstimPre using NLP minimizing z;

*        c) Iterate over increasingly better approximations

            mu = 0.1;
            zpen.l = 100;

            while(mu<1001,

*           d) Solve approximation. Stop if comp. gap < 0.01

                mu $ ((mu > 1000) or (zpen.l/mu lt 0.01)) = 100000;

                solve EstimPEN using NLP minimising zz;

                mu=mu*2;
            );

$batinclude "savesol.gms" a9
        );

*       -------------------- approach 10 ----------------------------------
*       a6 most often produces the best solution, but when it is NOT best,
*       it is much worse than a9 (penalty function). Thus, check which one
*       was the best and use that one.

        if(aon("a10"),

            if((aon("a6") and aon("a9")),
                if((objes("a9",m,n) lt objes("a6",m,n)),
                    objes("a10",m,n) = objes("a9",m,n);
                    pest(m,n,t,i,"a10") = pest(m,n,t,i,"a9");
                    cest(m,n,i,j,"a10") = cest(m,n,i,j,"a9");
                    sest(m,n,i,"a10")   = sest(m,n,i,"a9");
                else
                    objes("a10",m,n) = objes("a6",m,n);
                    pest(m,n,t,i,"a10") = pest(m,n,t,i,"a6");
                    cest(m,n,i,j,"a10") = cest(m,n,i,j,"a6");
                    sest(m,n,i,"a10")   = sest(m,n,i,"a6");
                );
            );
        );
    );
);

parameter objes1;
set objitems /pre2,dds,nds,ddx,ndx/;
objes1(m,n,a) = objes(a,m,n);
objes1(m,n,objitems) = objes(objitems,m,n);
display "All objective values transposed for readability:", objes1;
if(aon("a7"), display zenum);

* Rate solution approaches. Find the best solution, and give one point to
* the approach that found it

scalar bestz /0/;
parameter rating(*);
rating(aon) $ (not sameas(aon,"a0"))
= sum((m,n) $ (objes1(m,n,aon) lt (smin(aon2 $ (not sameas(aon2,"a0")),
    objes1(m,n,aon2)) + 0.0001)),1);

rating("total") = card(n)*card(m);
display rating;

option kill=rating;
option kill=objes;
option kill=objes1;

if(saveGDX,
execute_unload "smp%pnr%.gdx" m mrun n i t a aon ptru ctru stru etru pest cest sest);