$title Linearization Techniques for extremal-Nash Equilibria (LINEARNE,SEQ=382)

$onText
This model proposes a new approach to the computation
of extremal-Nash equilibria in a wholesale power market with
transmission constraints. The approach uses linearization techniques
to formulate the extremal-Nash equilibrium problem as
a single-stage mixed-integer linear programming problem which
can be solved with standard software. Through the introduced
concept of extremal-Nash equilibria, the derived structure can
efficiently locate all Nash equilibria of the game.

This example is loosely based on the NSW market. There are 7 portfolios
each controlling 2 units.

Reference:
Hesamzadeh, M R, and Biggar, D R, Computation of extremal-nash equilibria
in a wholesale power market using a single-stage MILP.
IEEE Transactions on Power Systems 27, 3 (2012), 1706-1707.

This code has been developed by:

Dr. Mohammad R. Hesamzadeh,
School of Electrical Engineering,
KTH Royal Institute of Technology, Sweden
(email:mrhesamzadeh@ee.kth.se)

Dr. Darryl R. Biggar,
Regulatory Development Branch,
Australian Competition and Consumer Commission, Australia
(email:darryl.biggar@accc.gov.au)

Keywords: extended mathematical programming, energy economics, Nash equilibria,
          linearization techniques, power markets
$offText

Set
   p 'portfolios owning strategic generating units'
     / delta, eraringp, macgen, sithe, uranq, redbank, snowy /
   u 'generating units'
     / colongra,    munmorah,   eraring,    shoalhaven, bayswater
       huntvgt,     sithe1,     sithe2,     uranq1,     uranq2
       redbank1,    redbank2,   snowy1,     snowy2,     delta-con
       eraring-con, macgen-con, tallawarra, VOLL                  /
   owns(p,u) 'ownership relation between portfolios and units'
             / delta.(colongra, munmorah),  eraringp.(eraring, shoalhaven)
               macgen.(bayswater, huntvgt), sithe.(sithe1, sithe2)
               uranq.(uranq1, uranq2),      redbank.(redbank1, redbank2)
               snowy.(snowy1, snowy2)                                      /
   us(u) 'list of strategic generators'
   un(u) 'list of non-strategic generators';

option us < owns;
un(u) = not us(u);

Table unitData(u,*) 'unit data'
                     cost   cap
   colongra         87.19   560
   munmorah         19.87  1080
   eraring          17.42   660
   shoalhaven      300.00   240
   bayswater        12.78  1200
   huntvgt         383.47    44
   sithe1           37.16    80
   sithe2           37.16    80
   uranq1           78.32   320
   uranq2           78.32   320
   redbank1         12.49   100
   redbank2         12.49    45
   snowy1          300.00  1250
   snowy2          300.00   666
   delta-con        19.87  3240
   eraring-con      17.42  1980
   macgen-con       12.78  3600
   tallawarra       27.54   422
   VOLL          10000.00 90000;

Set ds 'range of demands to be considered' / d1*d35 /;

Parameter
   cost(u)    'variable cost of generating unit u dollars per MWh'
   cap(u)     'maximum production capacity for unit u in MW'
   Demand(ds) 'demand in each run of the model'
   D          'total demand';

cost(u)    = unitData(u,'cost');
cap(u)     = unitData(u,'cap');
Demand(ds) = 8000 + (ord(ds) - 1)*200;

* The following variables will need to be changed according to the
* number of actions per unit
Set
   a           'possible actions for each unit'               / q1*q4  /
   s           'enumeration of strategies for each portfolio' / s1*s16 /
   k           'index of binary digits - b1 LSB'              / b1*b2  /
   sl(p,s,u,a) 'list of allowed strategies for each portfolio'
   ps(p,s)     'list of strategies defined for each portfolio';

abort$(round(card(a)**smax(p, sum(owns(p,u),1))) <> card(s))
   "scenarios and max number units in portfolio don't match";

abort$(round(2**card(k)) <> card(a)) "action and binary digits don't match";

* Enumeration all strategies
Parameter pow(u), pu(p), cont;

pu(p) = sum(owns(p,u), 1);
loop(p,
   pow(u) = 0;
   loop(s$(ord(s) <= round(card(a)**pu(p))),
      loop((owns(p,u),a)$sameas('q1',a), sl(p,s,u,a+pow(u)) = yes);
      cont = 1;
      loop(owns(p,u)$cont,
         if(pow(u) < card(a) - 1,
            pow(u) = pow(u)  + 1;
            cont   = 0;
         else
            pow(u) = 0;
         );
      );
   );
);
display sl;
option ps < sl;

Alias (p,cp), (s,cs);

Parameter
   b0   'additive offset for constant'
   b(k) 'conversion factor from binary to production';

b0   = power(2,- card(k));
b(k) = power(2,-(card(k) - ord(k) + 1));

Set y(a,k) 'mapping from unit action to share of output';
y(a,k) = mod(trunc(ord(a)/power(2,ord(k) - 1)),2);

Variable
   social_cost 'objective variable'
   pi0(p)      'profit of portfolio p'
   pi(cp,cs,p) 'ditto in case (cp cs)'
   pr0         'dual of sum(u,g0(u)) = D'
   pr(cp,cs)   'ditto in case (cp cs)';

Positive Variable
   g0(u)        'production of unit'
   g(cp,cs,u)   'ditto in case (cp cs)'
   gh0(u)       'volume of capacity offered to the market of stragegic unit'
   gh(cp,cs,u)  'ditto in case (cp cs)'
   v0(u)        'dual of g0(u) >= 0'
   v(cp,cs,u)   'ditto in case (cp cs)'
   w0(u)        'dual of g0(u) <= gh0'
   w(cp,cs,u)   'ditto in case (cp cs)'
   z0(u,k)      'z0(u,k) = w0(u)*x0(u,k)'
   z(cp,cs,u,k) 'ditto in case (cp cs)';

Binary Variable
   x0(u,k)      'select digits for volume offered'
   x(cp,cs,u,k) 'ditto in case (cp cs)';

Equation
   obj                'objective worst social cost'
   kkt0(u)            'KKT system'
   kkt(cp,cs,u)       'ditto in case (cp cs)'
   demandsat0         'demand satisfaction'
   demandsat(cp,cs)   'ditto in case (cp cs)'
   caplimit0(u)       'capacity limit'
   caplimit(cp,cs,u)  'ditto in case (cp cs)'
   defcap0(u)         'capacity definition for stratgic units'
   defcap(cp,cs,u)    'ditto in case (cp cs)'
   defprofit0(p)      'profit definition'
   defprofit(cp,cs,p) 'ditto in case (cp cs)'
   defscen(cp,cs,u,k) 'scenario capacity fixing'
   defNE(cp,cs)       'Nash equilibrium';

obj ..          social_cost =e= sum(u, cost(u)*g0(u));

* Optimal dispatch for the candidate WNE
kkt0(u)..      -cost(u) + pr0 + v0(u) - w0(u) =e= 0;

demandsat0..    sum(u, g0(u)) =e= D;

caplimit0(u)..  gh0(u)  =g= g0(u);

defcap0(us)..   gh0(us) =e= (b0 + sum(k, b(k)*x0(us,k)))*cap(us);

defprofit0(p).. pi0(p)  =e= sum(owns(p,us), (w0(us)*b0 + sum(k, b(k)*z0(us,k)))*cap(us));

gh0.fx(un) = cap(un);

* Complementary conditions
* disjunction: v0() = 0 or g0() = 0 hold, w0() = 0 or g0() = gh0()
Equation v0is0(u), g0is0(u), w0is0(u), g0isgh0(u);

v0is0(u)..   v0(u) =e= 0;

g0is0(u)..   g0(u) =e= 0;

w0is0(u)..   w0(u) =e= 0;

g0isgh0(u).. g0(u) =e= gh0(u);

* z0(us,k) = w0(us)*x0(us,k)
* Disjunction if x0() = 1 then z0isw0() holds else z0is0 holds
Equation z0is0(u,k), z0isw0(u,k);

z0is0(us,k)..  z0(us,k) =e= 0;

z0isw0(us,k).. z0(us,k) =e= w0(us);

* Optimal dispatch for the alternative case cp ss
kkt(ps,u)..       -cost(u) + pr(ps) + v(ps,u) - w(ps,u) =e= 0;

demandsat(ps)..    sum(u, g(ps,u)) =e= D;

caplimit(ps,u)..   gh(ps,u)  =g= g(ps,u);

defcap(ps,us)..    gh(ps,us) =e= (b0 + sum(k, b(k)*x(ps,us,k)))*cap(us);

defprofit(ps,p)..  pi(ps,p)  =e= sum(owns(p,us), (w(ps,us)*b0
                              +  sum(k, b(k)*z(ps,us,k)))*cap(us));

defscen(ps(cp,cs),us,k)..
   x(ps,us,k) =e= x0(us,k)$(not owns(cp,us)) + sum((sl(ps,us,a), y(a,k)),1)$owns(cp,us);

gh.fx(ps,un) = cap(un);

* Complementary conditions
* disjunction: v() = 0 or g() = 0 hold, w() = 0 or g() = gh()
Equation vis0(cp,cs,u), gis0(cp,cs,u), wis0(cp,cs,u), gisgh(cp,cs,u);

vis0(ps(cp,cs),u)..  v(ps,u) =e= 0;

gis0(ps(cp,cs),u)..  g(ps,u) =e= 0;

wis0(ps(cp,cs),u)..  w(ps,u) =e= 0;

gisgh(ps(cp,cs),u).. g(ps,u) =e= gh(ps,u);

* z(ps,us,k) = w(ps,us)*x(ps,us,k)
* Disjunction if x() = 1 then zisw() holds else zis0 holds
Equation zis0(cp,cs,u,k), zisw(cp,cs,u,k);

zis0(ps,us,k).. z(ps,us,k) =e= 0;

zisw(ps,us,k).. z(ps,us,k) =e= w(ps,us);

* These equations ensure that the selected strategy combination is a NE
defNE(ps(cp,cs)).. pi0(cp) =g= pi(ps,cp);

Model LinearNE / all /;

File emp / '%emp.info%' /;
put  emp '* problem %gams.i%' / 'default bigM 1e5';
loop(u,
   put / 'disjunction *' v0is0(u) 'else' g0is0(u);
   put / 'disjunction *' w0is0(u) 'else' g0isgh0(u);
);
loop((us,k),
   put / 'disjunction ' x0(us,k) z0isw0(us,k) 'else' z0is0(us,k);
);
loop((ps,u),
   put / 'disjunction *' vis0(ps,u) 'else' gis0(ps,u);
   put / 'disjunction *' wis0(ps,u) 'else' gisgh(ps,u);
);
loop((ps,us,k),
   put / 'disjunction ' x(ps,us,k) zisw(ps,us,k) 'else' zis0(ps,us,k);
);
putClose;

option optCr = 0, limRow = 0, limCol = 0, solPrint = off;

Parameter results(ds,*);
results(ds,'D')   = Demand(ds);
results(ds,'obj') = na;
loop(ds,
   D = Demand(ds);
   solve LinearNE using emp minimizing social_cost;
   results(ds,'obj')$[LinearNE.modelStat = %modelStat.optimal%] = LinearNE.objVal;
);

display results;
abort$[smax{ds, abs(mapval(results(ds,'obj')))} > 0]
   'Some cases not solved', results;
abort$[smin{ds$[ord(ds) > 1], results(ds,'obj') - results(ds,'obj')} < 0]
   'Objectives computed are not increasing', results;