pies.gms : PIES Energy Equilibrium

Description

A linear program with a variable rhs in the constraint system
is expressed as a complementarity problem.

LP:     min     <c,x>
                Ax = q(p),
        s.t.    Bx = b,
                x >= 0

where the prices p are the duals to the first constraint.

MCP:           A'p + B'v + c >= 0, x >= 0, comp.
        -Ax + q(p)           =  0, p free, comp.
        -Bx              + b =  0, v free, comp.

Of course, the variables x and v are going to be
split up further in the GAMS model.

References:
William W. Hogan, Energy Policy Models for Project Independence,
Computers \& Operations Research (2), 1975.

N. Josephy, A Newton Method for the PIES energy model,
Tech Report, Mathematics Research Center, UW-Madison, 1979.

S.P. Dirkse and M.C. Ferris, "MCPLIB: A Collection of Nonlinear Mixed
Complementarity Problems", Optimization Methods in Software 5 (1995), 319-345.

Contributor: Michael Ferris & Steven Dirkse, October 2010


Small Model of Type : EQUIL


Category : GAMS EMP library


Main file : pies.gms

$TITLE PIES Energy Equilibrium (PIES,SEQ=56)

$ontext
A linear program with a variable rhs in the constraint system
is expressed as a complementarity problem.

LP:     min     <c,x>
                Ax = q(p),
        s.t.    Bx = b,
                x >= 0

where the prices p are the duals to the first constraint.

MCP:           A'p + B'v + c >= 0, x >= 0, comp.
        -Ax + q(p)           =  0, p free, comp.
        -Bx              + b =  0, v free, comp.

Of course, the variables x and v are going to be
split up further in the GAMS model.

References:
William W. Hogan, Energy Policy Models for Project Independence,
Computers \& Operations Research (2), 1975.

N. Josephy, A Newton Method for the PIES energy model,
Tech Report, Mathematics Research Center, UW-Madison, 1979.

S.P. Dirkse and M.C. Ferris, "MCPLIB: A Collection of Nonlinear Mixed
Complementarity Problems", Optimization Methods in Software 5 (1995), 319-345.

Contributor: Michael Ferris & Steven Dirkse, October 2010
$offtext


SETS
  comod /
    C  'coal'
    L  'light oil'
    H  'heavy oil'
  /
  R  'resources' /
    C  'capital'
    S  'steel'
  /
  creg   'coal producing regions'         / 1 * 2 /
  oreg   'crude oil producing regions'    / 1 * 2 /
  ctyp   'increments of coal production'  / 1 * 3 /
  otyp   'increments of oil production'   / 1 * 2 /
  refin  'refineries'                     / 1 * 2 /
  users  'consumption regions'            / 1 * 2 /
  ;
alias (comod,cc);

PARAMETERS
  rmax(R) 'maximum resource usage' /
    C       35000
    S       12000
  /
  cmax(creg,ctyp) 'coal prod. limits' /
    1.1     300
    1.2     300
    1.3     400
    2.1     200
    2.2     300
    2.3     600
  /
  omax(oreg,otyp) 'oil prod. limits' /
    1.1     1100
    1.2     1200
    2.1     1300
    2.2     1100
  /
  rcost(refin) 'refining cost' /
    1       6.5
    2       5
  /
  q0(comod) 'base demand for commodities' /
    C       1000
    L       1200
    H       1000
  /
  p0(comod) 'base prices for commodities' /
    C       12
    L       16
    H       12
  /
  demand(comod,users) 'computed at optimality'
  output(refin,*)  '% output of light/heavy oil' /
    1.L     .6
    1.H     .4
    2.L     .5
    2.H     .5
  /
  ;
TABLE esub(comod,cc)    'cross-elasticities of substitution'
        C       L       H
C      -.75     .1      .2
L       .1     -.5      .2
H       .2      .1     -.5 ;

TABLE cruse(R,creg,ctyp) 'resource use in coal prod'
        1       2       3
C.1     1       5      10
C.2     1       5       6
S.1     1       2       3
S.2     1       4       5 ;

table oruse(R,oreg,otyp) 'resource use in oil prod'
        1       2
C.1     0      10
C.2     0      15
S.1     0       4
S.2     0       2 ;

table ccost(creg,ctyp)  'coal prod. cost'
        1       2       3
1       5       6       8
2       4       5       7 ;

table ocost(oreg,otyp)  'oil prod. cost'
        1       2
1       1       1.5
2       1.25    1.5 ;

table ctcost(creg,users)
        1       2
1       1       2.5
2       .75     2.75 ;

table otcost(oreg,refin)
        1       2
1       2       3
2       4       2 ;

table ltcost(refin,users)  'light oil trans costs'
        1       2
1       1       1.2
2       1       1.5 ;

table htcost(refin,users)  'heavy oil trans costs'
        1       2
1       1       1.2
2       1       1.5 ;


positive variables
  c(creg,ctyp)           'coal production'
  o(oreg,otyp)           'oil production'
  ct(creg,users)         'coal transportation levels'
  ot(oreg,refin)         'crude oil transportation levels'
  lt(refin,users)        'light transportation levels'
  ht(refin,users)        'heavy transportation levels'
  p(comod,users)         'commodity prices'
  ;
c.up(creg,ctyp) = cmax(creg,ctyp);
o.up(oreg,otyp) = omax(oreg,otyp);

equations
  dembal(comod,users)  'excess supply of product'
  cmbal(creg)          'coal material balance'
  ombal(oreg)          'oil material balance'
  lmbal(refin)         'light material balance'
  hmbal(refin)         'heavy material balance'
  ruse(R)              'resource use constraints'
  ;

dembal(comod,users) ..
  (sum(creg,ct(creg,users)))$[sameas(comod,'C')]
  + (sum(refin,lt(refin,users)))$[sameas(comod,'L')]
  + (sum(refin,ht(refin,users)))$[sameas(comod,'H')]
  =g=
  q0(comod) * prod(cc, (p(cc,users)/p0(cc))**esub(comod,cc));

cmbal(creg) ..
  sum(ctyp,c(creg,ctyp))  =e=  sum(users,ct(creg,users));

ombal(oreg) ..
  sum(otyp,o(oreg,otyp))  =e=  sum(refin,ot(oreg,refin));

lmbal(refin) ..
  sum(oreg, ot(oreg,refin)) * output(refin,"L") =e=
  sum(users,lt(refin,users));

hmbal(refin) ..
  sum(oreg, ot(oreg,refin)) * output(refin,"H") =e=
  sum(users,ht(refin,users));

ruse(R) ..
  rmax(R) =g=
  sum(creg, sum(ctyp, c(creg,ctyp)*cruse(R,creg,ctyp)))
  + sum(oreg, sum(otyp, o(oreg,otyp)*oruse(R,oreg,otyp)));

option limrow = 0;
option limcol = 0;
option iterlim = 1000;
option reslim = 120;

table i_c(creg,ctyp)
        1       2       3
1       300     300     400
2       200     300     600 ;

table i_o(oreg,otyp)
        1       2
1       1100    1000
2       1300    1000 ;

table i_ct(creg,users)  'initial trans'
        1       2
1       0       828
2       1016    84 ;

table i_ot(oreg,refin)  'initial trans'
        1       2
1       2075    0
2       0       2358 ;

table i_lt(refin,users) 'initial trans'
        1       2
1       22      1223
2       1179    0 ;

table i_ht(refin,users) 'initial trans'
        1       2
1       0       830
2       998     180 ;

table iprice(comod,users) 'initial price estimate'
        1       2
C       11.7    13.7
L       15.8    16.0
H       11.9    12.4 ;

c.l(creg,ctyp) = i_c(creg,ctyp);
o.l(oreg,otyp) = i_o(oreg,otyp);
ct.l(creg,users) = i_ct(creg,users);
ot.l(oreg,refin) = i_ot(oreg,refin);
lt.l(refin,users) = i_lt(refin,users);
ht.l(refin,users) = i_ht(refin,users);
p.lo(comod,users) = .1;
* p.fx(comod,users) = iprice(comod,users);
p.l(comod,users) = iprice(comod,users);

* initialize the multiplies.....
cmbal.m(creg) = 1;
ombal.m(oreg) = 1;
lmbal.m(refin) = 1;
hmbal.m(refin) = 1;
ruse.m(R) = 1;

variables obj;
equation defobj;

defobj.. obj =e= sum((creg,ctyp), ccost(creg,ctyp)*c(creg,ctyp))
         + sum((oreg,otyp), ocost(oreg,otyp)*o(oreg,otyp))
         + sum((creg,users), ctcost(creg,users)*ct(creg,users))
         + sum((oreg,refin), (otcost(oreg,refin) + rcost(refin))*ot(oreg,refin))
         + sum((refin,users), ltcost(refin,users)*lt(refin,users))
         + sum((refin,users), htcost(refin,users)*ht(refin,users));

model piesemp / defobj, dembal, cmbal, ombal, lmbal, hmbal, ruse /;

file myinfo /'%emp.info%'/;
putclose myinfo 'dualVar p dembal';

solve piesemp using emp min obj;

* the following sets up the MCP generated explicitly

positive variables
  cv(creg)        'dual to cmbal'
  ov(oreg)
  lv(refin)
  hv(refin)
  mu(R)           'dual to ruse cons., i.e. marginal utility'
  ;
equations
  delc(creg, ctyp)
  delo(oreg, otyp)
  delct(creg,users)
  delot(oreg,refin)
  dellt(refin,users)
  delht(refin,users)
  ;

delc(creg,ctyp) ..
  ccost(creg,ctyp) + sum(R, cruse(R,creg,ctyp)*mu(R))
  =g=   cv(creg);

delo(oreg,otyp) ..
  ocost(oreg,otyp) + sum(R, oruse(R,oreg,otyp)*mu(R))
  =g=   ov(oreg);

delct(creg,users) ..
  ctcost(creg,users) + cv(creg)  =g=  p("C",users);

delot(oreg,refin) ..
  otcost(oreg,refin) + rcost(refin) + ov(oreg)   =g=
  output(refin,"L") * lv(refin) + output(refin,"H") * hv(refin);

dellt(refin,users) ..
  ltcost(refin,users) + lv(refin)
  =g=  p("L",users);

delht(refin,users) ..
  htcost(refin,users) + hv(refin)
  =g=  p("H",users);

model pies / delc.c, delo.o, delct.ct, delot.ot, dellt.lt, delht.ht,
        dembal.p, cmbal.cv, ombal.ov, lmbal.lv, hmbal.hv, ruse.mu /;

cv.l(creg) = cmbal.m(creg);
ov.l(oreg) = ombal.m(oreg);
lv.l(refin) = lmbal.m(refin);
hv.l(refin) = hmbal.m(refin);
mu.l(R) = ruse.m(R);

pies.iterlim = 0;
solve pies using mcp;
abort$[pies.solveStat <> %SOLVESTAT.NORMAL COMPLETION%] 'EMP solution should be optimal for PIES MCP';
abort$[pies.modelStat <> %MODELSTAT.OPTIMAL%] 'EMP solution should be optimal for PIES MCP';

file  out /pies.out/;
put out;

put "Coal Prod:      type 1     type 2     type 3" /;
put "region 1   ", c.l("1","1"):11:3, c.l("1","2"):11:3, c.l("1","3"):11:3 /;
put "region 2   ", c.l("2","1"):11:3, c.l("2","2"):11:3, c.l("2","3"):11:3 /;
put /;
put "Oil Prod:       type 1     type 2" /;
put "region 1   ", o.l("1","1"):11:3, o.l("1","2"):11:3 /;
put "region 2   ", o.l("2","1"):11:3, o.l("2","2"):11:3 /;
put /;
put "Coal Trans:     user 1     user 2" /;
put "region 1   ", ct.l("1","1"):11:3, ct.l("1","2"):11:3 /;
put "region 2   ", ct.l("2","1"):11:3, ct.l("2","2"):11:3 /;
put /;
put "Oil Trans:     refin 1    refin 2" /;
put "region 1   ", ot.l("1","1"):11:3, ot.l("1","2"):11:3 /;
put "region 2   ", ot.l("2","1"):11:3, ot.l("2","2"):11:3 /;
put /;
put "Light Trans:    user 1     user 2" /;
put "refin 1    ", lt.l("1","1"):11:3, lt.l("1","2"):11:3 /;
put "refin 2    ", lt.l("2","1"):11:3, lt.l("2","2"):11:3 /;
put /;
put "Heavy Trans:    user 1     user 2" /;
put "refin 1    ", ht.l("1","1"):11:3, ht.l("1","2"):11:3 /;
put "refin 2    ", ht.l("2","1"):11:3, ht.l("2","2"):11:3 /;
put /;
put "Prices:         user 1     user 2" /;
put "Coal       ", p.l("C","1"):11:3, p.l("C","2"):11:3 /;
put "Light oil  ", p.l("L","1"):11:3, p.l("L","2"):11:3 /;
put "Heavy oil  ", p.l("H","1"):11:3, p.l("H","2"):11:3 /;
put /;
demand(comod,users) = q0(comod) *
         prod(cc, (p.l(cc,users)/p0(cc))**esub(comod,cc));
put "Demand:         user 1     user 2" /;
put "Coal       ", demand("C","1"):11:3, demand("C","2"):11:3 /;
put "Light oil  ", demand("L","1"):11:3, demand("L","2"):11:3 /;
put "Heavy oil  ", demand("H","1"):11:3, demand("H","2"):11:3 /;
put /;
put "Capital usage: ", (0-ruse.l("C")):10:2, "    dual price: ", mu.l("C"):9:3 /;
put "  Steel usage: ", (0-ruse.l("S")):10:2, "    dual price: ", mu.l("S"):9:3 /;