$title  The Illustrative PIES Model -- Hogan, William (1975)

*	"Energy policy models for Project Independence", 
*	Computers and Operations Research 2

set     j       Consumption regoins /j1, j2/
        i       Coal regions        /i1, i2/
        k       Oil regions         /k1, k2/
        r       Refineries          /r1, r2/

        c       Supply increments for coal regions /L, M, H/,
        o       Supply increments for oil regoins  /L, H/,

        p       Energy products         /Coal, Light, Heavy/,
        g(p)    Grades of refined oil   / Light, Heavy/,

	res	Resources		/steel, newcap/;

alias (p,pp,ip,jp);

*	------------------------------------------ end of set definitions

table table1(i,c,*) Resource requirements for production levels (p 257)

*       cap     Production capacity (tons per day)
*       c0      Minimum price / ton ($)
*       newcap  New capital / ton
*       steel   Steel / ton

        cap     c0      newcap  steel
i1.L    300     5       1       1
i1.M    300     6       5       2
i1.H    400     8       10      3

i2.L    200     4       1       1
i2.M    300     5       5       4
i2.H    600     7       6       5;

table table2(i,j)  Transport costs ($ per ton)

                j1      j2
        i1      1.00    2.50
        i2      0.75    2.75;


table table3(k,o,*)  Oil resource requirements

                cap     c0      newcap  steel
        k1.L    1100    1       0       0
        k1.H    1200    1.5     10      4
        k2.L    1300    1.25    0       0
        k2.H    1100    1.50    15      2;

table table4(k,r)  Oil transport costs ($ per barrel)
                r1      r2
        k1      2       3
        k2      4       2;

table table5(*,r)       Refinery yields and cost
                        r1      r2
        Light           0.6     0.5
        Heavy           0.4     0.5
        cost            6.5     5.0;

table table6(r,j)  Transport costs for refined products ($ per barrel)
                j1      j2
        r1      1       1.2
        r2      1       1.5;

table  table7(p,*)  Elasticities of final demand

                RefP    RefQ    Light   Heavy   Coal
        Light   16      1200   -0.5     0.2     0.1
        Heavy   12      1000    0.1    -0.5     0.2
        Coal    12      1000    0.1     0.2    -0.75;

table  table8(p,j) Demand without K or S constraints


          j1        j2
light    1252      1266
heavy    1041      1055
coal     1102       998   ;

table  table9(p,j) Demand with  K or S constraints

          j1        j2
light    1205      1229
heavy     996      1020
coal      996       910   ;

*	------------------------------------------ end of tabular data
*	The model assumes the same demand function in all regions:

parameter	pref(p)		Reference price, 
		qref(p)		Reference demand;

pref(p) = table7(p,"refp");
qref(p) = table7(p,"refq");


parameter	sigma(p,pp)	Asymmetric slutsky matrix used for the AS model
		slutsky(p,pp)	Symmetric slutsky matrix used for CP model;

*	Asymmetric demand system:

sigma(p,pp) = table7(p,pp)*qref(p)/pref(pp);

*	Related symmetric demand system:

slutsky(p,pp) = 0.5*(sigma(p,pp)+sigma(pp,p));

*	Calculate the inverse demand coefficients:

variable	SLUTSKYINV(p,p)	Inverse demand curve;

equation	invdef;

invdef(ip,jp)..	sum(p, SLUTSKYINV(ip,p)*slutsky(p,jp)) =e= 1$sameas(ip,jp);

model invert /invdef/;
solve invert using mcp;

*	------------------------------------------ end of demand system calculation

parameter	delta(p,j)	Demand adjustment (off-diagonal terms),
		rs(res)		Resource supply,
		resutil		Resource utilization,
		report		Pivot report data,
		op		Flag for the diagonal own-price demand model (QP or MCP) /1/
		cp		Flag for symmetric cross-price demand model (QP or MCP) /0/
		as		Flag for the asymmetric cross-price demand model (MCP) /0/;


rs(res) = 0;
delta(p,j) = 0;

*	------------------------------------------ end of assigned parameters

*	Set up a linear programming model to check the formulation:

NONNEGATIVE
VARIABLES
	QC(i,c)		Quantity of coal extracted by region i at cost level c,
	QO(k,o)		Quantity of oil resources extracted -- region k at cost level o
	QR(r)		Quantity of oil refined -- refinery r,

	D(p,j)		Demand -- energy product p in consumption region j,

	XC(i,j)		Quantity of coal transported from region i to market j
	XO(k,r)		Quantity of oil resources shipped from region k to refinery r
	XR(r,g,j)	Quantity of oil grade g transported from refinery r to market j;

variable	COST		Total cost;

equations	oilresource, crudeoil, refinedoil, coalsupply, demand, costdef, resource;

*	Oil supply from region (k) by cost increment (o) equals oil shipments to refineries (r):

oilresource(k).. sum(o, QO(k,o)) =g= sum(r, XO(k,r));

*	Oil supply from regions (k) equals refinery (r) output:

crudeoil(r)..	sum(k, XO(k,r)) =g= QR(r);

*	Refined oil supply equals refined oil shipments:

refinedoil(r,g).. QR(r)*table5(g,r) =g= sum(j, XR(r,g,j));

*	Coal supply equals coal demand:

coalsupply(i)..	sum(c, QC(i,c)) =g= sum(j, XC(i,j));


demand(j,p)..	sum((r,g(p)), XR(r,g,j)) + sum(i, XC(i,j))$sameas(p,"coal") =g= D(p,j) + delta(p,j);


resource(res)$rs(res)..	rs(res) =g= sum((i,c), table1(i,c,res)*QC(i,c))  + 
				    sum((k,o), table3(k,o,res)*QO(k,o));

costdef..	COST =e=  sum((r,g,j), XR(r,g,j)*table6(r,j)) + 
			  sum((k,r),   XO(k,r)*table4(k,r)) +
			  sum((i,j),   XC(i,j)*table2(i,j)) +

			  sum((k,o),   QO(k,o)*table3(k,o,"c0")) +
			  sum((i,c),   QC(i,c)*table1(i,c,"c0")) +

			  sum(r,       QR(r)*table5("cost",r));

*	Upper bounds on coal and oil production:

QC.UP(i,c) = table1(i,c,"cap");
QO.UP(k,o) = table3(k,o,"cap");

*	Solve the model as a linear program with fixed product demand
*	in each market and no resource constraints:

D.FX(p,j) = table8(p,j);

model pies_lp /all/;
solve pies_lp using LP minimizing COST;

report("price",p,j,"LP") = demand.m(j,p);
report("quantity",p,j,"LP") = D.l(p,j);
resutil(res,"LP") = sum((i,c), table1(i,c,res)*QC.L(i,c))  + 
			sum((k,o), table3(k,o,res)*QO.L(k,o));

report("price",p,j,"LP") = demand.m(j,p);
report("quantity",p,j,"LP") = D.l(p,j);

*	------------------------------------------ end of LP cross check

*	Define social surplus models (quadratic programs) with own-price
*	and cross-price demand systems:

*	Consumer surplus -- integral of the inverse demand curve, defined for
*	own-price (_op) demand and for cross-price (_cp) demand:

$macro	CS_op(p,j)	(D(p,j)*(pref(p) + 1/sigma(p,p) * (D(p,j)/2-qref(p))))

$macro	CS_cp(p,j)	(D(p,j)*(pref(p) + sum(pp, slutskyinv.L(p,pp) * (D(pp,j)/2-qref(pp)))))

variable	NSS	Negative of social surplus;

equation	nssdef;

nssdef..	NSS =e= COST - sum((p,j), CS_op(p,j)$op + CS_cp(p,j)$cp);

model pies_qcp /all/;

*	Use CPLEX as the quadratic programming solver:

option qcp=cplex; 

*	Solve the market equilibrium with diagonal demand 
*	as a quadratic program:

op = yes;
cp = no;
as = no;

delta(p,j) = 0;
D.UP(p,j) = +inf;
D.LO(p,j) = 0;
rs(res) = 0;
solve pies_qcp using QCP minimizing NSS;

*	Report the equilibrium values with own-price demand function, ignoring
*	cross-price effects:

resutil(res,"op") = sum((i,c), table1(i,c,res)*QC.L(i,c))  + 
				sum((k,o), table3(k,o,res)*QO.L(k,o));
report("price",p,j,"op") = demand.m(j,p);
report("quantity",p,j,"op") = D.l(p,j);

*	Solve the market equilibrium with cross-price demand 
*	as a quadratic program:

op = no;
cp = yes;
as = no;

delta(p,j) = 0;
D.UP(p,j) = +inf;
D.LO(p,j) = 0;
rs(res) = 0;
solve pies_qcp using QCP minimizing NSS;

resutil(res,"cp") = sum((i,c), table1(i,c,res)*QC.L(i,c))  + 
				sum((k,o), table3(k,o,res)*QO.L(k,o));
report("price",p,j,"cp") = demand.m(j,p);
report("quantity",p,j,"cp") = D.l(p,j);

*	------------------------------------------ end of QP setup and diagonal solution

*	Define the complementarity model.

NONNEGATIVE VARIABLES
	P_O(k)		Supply price of oil,
	P_XO(r)		Delivered price of oil,
	P_R(r,g)	Price of refinery outputs,

	P_C(i)		Supply price of coal,
	P_D(j,p)	Demand price of all products (oil and coal)

	PR(res)		Resource price;

equations
	prf_QC(i,c)		Quantity of coal extracted by region i at cost level c,
	prf_QO(k,o)		Quantity of oil resources extracted -- region k at cost level o
	prf_QR(r)		Quantity of oil refined -- refinery r,
	def_D(p,j)		Demand -- energy product p in consumption region j,
	prf_XR(r,g,j)		Quantity of oil transported from refinery r to market j
	prf_XC(i,j)		Quantity of coal transported from region i to market j
	prf_XO(k,r)		Quantity of oil resources shipped from region k to refinery r
	coalprice(j,p)		Coal price constraint;

prf_QC(i,c)..		table1(i,c,"c0") + sum(res,PR(res)*table1(i,c,res)) =e= P_C(i);

prf_QO(k,o)..		table3(k,o,"c0") + sum(res,PR(res)*table3(k,o,res)) =e= P_O(k) ;

prf_QR(r)..		P_XO(r) + table5("cost",r) =g= sum(g,P_R(r,g)*table5(g,r));

prf_XR(r,g,j)..		P_R(r,g) + table6(r,j)	=g= P_D(j,g);

prf_XC(i,j)..		P_C(i) + table2(i,j) =G= P_D(j,"coal");

prf_XO(k,r)..		P_O(k) + table4(k,r) =G= P_XO(r);

*	Define the own-price demand function:

$macro	D_op(p,j,price)  (qref(p) + slutsky(p,p)*(price(j,p)-pref(p)))

*	Define the symmetric cross-price demand function:

$macro	D_cp(p,j,price)  (qref(p) + sum(pp, slutsky(p,pp)*(price(j,pp)-pref(pp))))

*	Define the asymmetric cross-price demand function:

$macro	D_as(p,j,price)  (qref(p) + sum(pp, sigma(p,pp)*(price(j,pp)-pref(pp))))

def_D(p,j)..	D(p,j) =e= D_op(p,j,P_D)$op + D_cp(p,j,P_D)$cp + D_as(p,j,P_D)$as; 

model pies_mcp /
	oilresource.P_O, crudeoil.P_XO, refinedoil.P_R, coalsupply.P_C, demand.P_D,
	prf_QC.QC, prf_QO.QO, prf_QR.QR, def_D.D, prf_XR.XR, prf_XC.XC, prf_XO.XO,
	resource.PR /;


*	Verify that the equilibrium solved with the social surplus formlation
*	solves the MCP model with a symmetric cross-price demand curve:

op = yes;
cp = no;
as = no;

delta(p,j) = 0;
solve pies_qcp using QCP minimizing NSS;

PR.FX(res) = 0;
P_O.L(k) = oilresource.M(k);
P_XO.L(r) = crudeoil.M(r);
P_R.L(r,g) = refinedoil.M(r,g);
P_C.L(i) = coalsupply.m(i);
P_D.L(j,p) = demand.M(j,p);

pies_mcp.iterlim = 0;
solve pies_mcp using mcp;
abort$round(pies_mcp.objval) "MCP fails to replicate the QP model (own-price demand)";

op = no;
cp = yes;
as = no;

delta(p,j) = 0;
solve pies_qcp using QCP minimizing NSS;

PR.FX(res) = 0;
P_O.L(k) = oilresource.M(k);
P_XO.L(r) = crudeoil.M(r);
P_R.L(r,g) = refinedoil.M(r,g);
P_C.L(i) = coalsupply.m(i);
P_D.L(j,p) = demand.M(j,p);

pies_mcp.iterlim = 0;
solve pies_mcp using mcp;
abort$round(pies_mcp.objval) "MCP fails to replicate the QP model (cross-price demand)";

*	------------------------------------------ end of MCP setup and diagonal solution

*	Solve the model iteratively with both diagonal and cross-price demand systems.

set	iter	Iterations for demand adjustments /iter0*iter10/;

parameter	iterlog		Iteration log for iterative demand adjustments,
		qd(p,j)		Exact demand at current prices,
		dev		Maximum deviation /1/;

op = yes;
cp = no;
as = no;

dev = 1;
loop(iter$round(dev,4),
	iterlog(iter,p) = sum(j, sqr(delta(p,j) - (D_as(p,j,demand.M) - D.L(p,j))));
	iterlog(iter,j) = sum(p, sqr(delta(p,j) - (D_as(p,j,demand.M) - D.L(p,j))));
	delta(p,j) = D_as(p,j,demand.M) - D.L(p,j);
	solve pies_qcp using QCP minimizing NSS;
	dev = sum(p,iterlog(iter,p)) + sum(j,iterlog(iter,j));
);
display "Iteration log with diagonal (own-price) demand system:", iterlog;

resutil(res,"iter_op") = sum((i,c), table1(i,c,res)*QC.L(i,c))  + 
				sum((k,o), table3(k,o,res)*QO.L(k,o));
report("price",p,j,"iter_op") = demand.m(j,p);
report("quantity",p,j,"iter_op") = D.l(p,j) + delta(p,j);
report("delta%",p,j,"iter_op") = 100 * delta(p,j) /(D.l(p,j) + delta(p,j));


iterlog(iter,p) = 0;
iterlog(iter,j) = 0;
delta(p,j) = 0;

op = no;
cp = yes;
as = no;

dev = 1;
loop(iter$round(dev,4),
	iterlog(iter,p) = sum(j, sqr(delta(p,j) - (D_as(p,j,demand.M) - D.L(p,j))));
	iterlog(iter,j) = sum(p, sqr(delta(p,j) - (D_as(p,j,demand.M) - D.L(p,j))));
	delta(p,j) = D_as(p,j,demand.M) - D.L(p,j);
	solve pies_qcp using QCP minimizing NSS;
	dev = sum(p,iterlog(iter,p)) + sum(j,iterlog(iter,j));
);
display "Iteration log with symmetric cross-price demand system:", iterlog;


resutil(res,"iter_cp") = sum((i,c), table1(i,c,res)*QC.L(i,c))  + 
				sum((k,o), table3(k,o,res)*QO.L(k,o));
report("price",p,j,"iter_cp") = demand.m(j,p);
report("quantity",p,j,"iter_cp") = D.l(p,j) + delta(p,j);
report("delta%",p,j,"iter_cp") = 100 * delta(p,j) /(D.l(p,j) + delta(p,j));

*	Verify that the equilibrium is consistent with the MCP model:

op = no;
cp = no;
as = yes;

D.L(p,j) = D_as(p,j,demand.M);

delta(p,j) = 0;

PR.FX(res) = 0;
P_O.L(k) = oilresource.M(k);
P_XO.L(r) = crudeoil.M(r);
P_R.L(r,g) = refinedoil.M(r,g);
P_C.L(i) = coalsupply.m(i);
P_D.L(j,p) = demand.M(j,p);

pies_mcp.iterlim = 0;
solve pies_mcp using mcp;
abort$round(pies_mcp.objval) "MCP fails to replicate asymmetric equilibrium.";

*	-----------------------------------------------------------------------

*	Install constraints on steel and new capital and then solve the
*	model with diagonal demand system (ignoring cross-price elasticities):

rs("steel")  = 12000;
rs("newcap") = 35000;

op = yes;
cp = no;
as = no;

delta(p,j)    = 0;

solve pies_qcp using QCP minimizing NSS;

resutil(res,"con_op") = sum((i,c), table1(i,c,res)*QC.L(i,c))  + 
			sum((k,o), table3(k,o,res)*QO.L(k,o));
report("price",p,j,"con_op") = demand.m(j,p);
report("quantity",p,j,"con_op") = D.l(p,j);

op = no;
cp = yes;
as = no;

delta(p,j)    = 0;

solve pies_qcp using QCP minimizing NSS;

resutil(res,"con_cp") = sum((i,c), table1(i,c,res)*QC.L(i,c))  + 
			sum((k,o), table3(k,o,res)*QO.L(k,o));
report("price",p,j,"con_cp") = demand.m(j,p);
report("quantity",p,j,"con_cp") = D.l(p,j);

*	Solve the QP model iteratively with cross-price elasticities of demand:

op = no;
cp = yes;
as = no;
dev = 1;
iterlog(iter,p) = 0;
iterlog(iter,j) = 0;

loop(iter$round(dev,4),
	iterlog(iter,p) = sum(j, sqr(delta(p,j) - (D_as(p,j,demand.M) - D.L(p,j))));
	iterlog(iter,j) = sum(p, sqr(delta(p,j) - (D_as(p,j,demand.M) - D.L(p,j))));
	delta(p,j) = D_as(p,j,demand.M) - D.L(p,j);
	solve pies_qcp using QCP minimizing NSS;
	dev = sum(p,iterlog(iter,p)) + sum(j,iterlog(iter,j));
);
display "Iteration log with resource constraints:", iterlog;

resutil(res,"iter_con") = sum((i,c), table1(i,c,res)*QC.L(i,c))  + 
				sum((k,o), table3(k,o,res)*QO.L(k,o));
report("price",p,j,"iter_con") = demand.m(j,p);
report("quantity",p,j,"iter_con") = D.l(p,j) + delta(p,j);
report("delta%",p,j,"iter_con") = 100 * delta(p,j) /(D.l(p,j) + delta(p,j));

*	Verify consistency with the MCP model:

op = no;
cp = no;
as = yes;

D.L(p,j) = D_as(p,j,demand.M);
delta(p,j) = 0;

PR.UP(res) = +inf;
PR.L(res) = resource.M(res);
P_O.L(k) = oilresource.M(k);
P_XO.L(r) = crudeoil.M(r);
P_R.L(r,g) = refinedoil.M(r,g);
P_C.L(i) = coalsupply.m(i);
P_D.L(j,p) = demand.M(j,p);

pies_mcp.iterlim = 0;
solve pies_mcp using mcp;
abort$round(pies_mcp.objval) "MCP fails to replicate constrained equilibrium.";


option resutil:1;
display resutil;

option report:2:2:1;
display report;