$title Equilibrium model with and without containing bilevel (EQUIL_BILEVEL,SEQ=99)

$onText

Model to demonstrate that the agents in an equilibrium model are
specified in the EMPinfo file and translated by JAMS in the same way
whether they are:
  A. part of an equilibrium model (no outer optimizing agent) or
  B. part of an inner model (outer optimizing agent with bilevel)

To show this we first create an equilibrium model with two VI agents,
and then create a bilevel model whose inner problem reduces to the
equilibrium model when the outer variables are fixed.  The outer
optimization is created so that solutions to the inner problem and
equilibrium problem are the same, but in such a way that this depends
on the implicit relationship between the inner and outer variables
enforced by the VI agents.

A casual inspection of the scalar models produced by JAMS is enough to
demonstrate the similarity between the equilibrium and inner models.

Contributor: Steven Dirkse, August 2013

$offText

sets
  i  / i1 * i2 /
  ii / 1  *  2 /
  j  / j1 * j4 /
  k  / k1 * k5 /
  t  / t1 * t2 /
  ;
alias(j,jj);

parameters
  c(j) /
    j1  4
    j2  -9
    j3  6
    j4  2
  /
  b(i) /
    i1    2
    i2    1
  /
  bb(ii) /
    1    -1
    2     1
  /;
table Q(jj,jj)
    j1  j2  j3  j4
j1   4  -2  -1  -1
j2  -2   4  -2  -1
j3  -1  -2   4  -2
j4  -1  -1  -2   4
;
table A(i,j)
    j1  j2  j3  j4
i1   1           1
i2       1   1   1
;
table AA(ii,k)
    k1  k2  k3  k4  k5
1    1           1   1
2        1   1   1
;
table E(t,k)
    k1  k2  k3  k4  k5
t1   2  -1   2   2
t2   -1  3       3   3
;

scalars
  out1bar 'target value for out1 at outer obj solution'    / 1 /
  v1bar   "target value for v('j1') at out obj solution"
  d1      'dv1/dout1' / -0.5 /
  out2bar 'target value for out1 at outer obj solution'    / 2 /
  w2bar   "target value for w('k2') at out obj solution"
  d2      'dw2/dout2' / [5/6] /
  ;
positive variable v(j);
positive variable w(k);
free variables
  out1   'var for outer min model'
  out2   'var for outer min model'
  s(t)
  ;
equation F1(j) 'VI func perp to v';
equation g(i);
equation F2(t) 'VI func perp to s';
equation sdef(t);
equation h(ii);

F1(j) .. sum{jj, Q(j,jj)*v(jj)} + c(j) =N= 0;
g(i)  .. sum {j, A(i,j)*v(j)} =G= b(i) - out1;

F2(t)   .. 2 * s(t) =N= 0;
sdef(t) .. sum{k, E(t,k)*w(k)} =E= s(t);
h(ii)   .. sum {k, AA(ii,k)*w(k)} =G= bb(ii) + out2;

v.up(j) = 20;
v.up('j2') = 2;
out1.fx = out1bar;
out2.fx = out2bar;

model equil / F1, g, F2, sdef, h /;

file opt1 / 'jams.opt' /;
putclose opt1 'fileName scalar1.gms'
            / 'dict dict1.txt'
            /;
file empinfo  / '%emp.info%' /;
putclose empinfo  'equilibrium'
                / '  vi     F1 v  g' /
                / '  vi  w  F2 s  sdef h' / ;

equil.holdfixed = 1;
equil.optfile = 1;
option limrow=99, limcol=99;
solve equil using emp;
abort$[equil.solvestat <> %solveStat.normalCompletion%] 'equil not solved', equil.solvestat;

v1bar = v.l('j1');
w2bar = w.l('k2');

scalars delta, o0, o1, odiff, marg;

* this ontext/offtext block checks correctness of d1 := dv1/dout1
$onText
delta = 0.01;
o0 = v.l('j1');

out1.fx = out1bar + delta;
solve equil using emp;
abort$[equil.solvestat <> %solveStat.normalCompletion%] 'equil not solved', equil.solvestat;

o1 = v.l('j1');

odiff = o1 - o0;
marg = odiff / delta;

display delta, o0, o1, odiff, marg, d1;
$exit
$offText

* this ontext/offtext block checks correctness of d1 := dv1/dout1
$onText
delta = 0.01;
o0 = w.l('k2');

out2.fx = out2bar + delta;
solve equil using emp;
abort$[equil.solvestat <> %solveStat.normalCompletion%] 'equil not solved', equil.solvestat;

o1 = w.l('k2');

odiff = o1 - o0;
marg = odiff / delta;

display delta, o0, o1, odiff, marg, d2;
$exit
$offText

delta = 0.01;
scalars ztarget, c_out1, c_v1, c_out2, c_w2;

c_out1 = out1bar + delta;
c_v1   = v1bar - (1/d1) * delta;
c_out2 = out2bar + delta;
c_w2   = w2bar - (1/d2) * delta;
ztarget = sqr(delta) + sqr(delta/d1)
        + sqr(delta) + sqr(delta/d2);

variable z 'outer objective';
equation zdef;
zdef .. sqr(v('j1') - c_v1) + sqr(out1 - c_out1)
      + sqr(w('k2') - c_w2) + sqr(out2 - c_out2)
        =e= z;

model bilev / F1, g, F2, sdef, h
              zdef
            /;

file opt2 / 'jams.op2' /;
putclose opt2 'fileName scalar2.gms'
            / 'dict dict2.txt'
           /;
putclose empinfo  'bilevel out1 out2'
                / '  vi     F1 v  g' /
                / '  vi  w  F2 s  sdef h' / ;

out1.lo = out1bar - 1;
out1.up = out1bar + 1;
out1.lo = -INF;
out1.up =  INF;
out2.lo = out2bar - 1;
out2.up = out2bar + 1;
out2.lo = -INF;
out2.up =  INF;

bilev.optfile = 2;
solve bilev using emp min z;
abort$[bilev.solvestat <> %solveStat.normalCompletion%] 'bilev not solved', bilev.solvestat;

scalars
 v1, w2
 tol / 1e-5 /;
v1 = v.l('j1');
w2 = w.l('k2');

display z.l, ztarget, out1.l, out1bar, v1, v1bar,
        out2.l, out2bar, w2, w2bar;
abort$[abs(z.l - ztarget) > tol] 'bilevel objective not as expected', z.l, ztarget;
abort$[abs(out1.l - out1bar) > tol] 'variable out1 not as expected', out1.l, out1bar;
abort$[abs(out2.l - out2bar) > tol] 'variable out2 not as expected', out2.l, out2bar;

* uncomment these to keep all the intermediate files around
execute 'rm -f scalar1.??? scalar2.??? dict1.txt dict2.txt'
execute 'rm -f scalar?pf.pf jams.opt jams.op2'