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 /;
``````