Monday, July 29, 2013

Microeconomic Principles

Another model illustrating tables 3.3 and 3.4 from “Complementarity Modeling in Energy Markets”.

$ontext

  
Example from chapter 3: "Some Microeconomic Principles"

  
Reference:
     
Steven A. Gabriel, Antonio J. Conejo, J. David Fuller,
     
Benjamin F. Hobbs, Carlos Ruiz
     
Complementarity Modeling in Energy Markets
     
Springer 2012

$offtext


SETS
 i 
'firms' /i1,i2,i3/
 b 
'production process (mines)' /b1,b2/
;

TABLE c(i,b) 'variable cost'
      
b1    b2
 
i1  0.55  0.81
 
i2  0.62  1.25
 
i3  0.78  1.35
;

TABLE K(i,b) 'process capacities'
      
b1     b2
 
i1  21000  16000
 
i2  17000  22000
 
i3  18000  14000
;

parameters
   alpha
'inverse demand intercept'  / 2.5 /
   beta 
'inverse demand slope'      / 0.000016666666667 /
;


*-------------------------------------------------------------------------------
* Reporting macro
*-------------------------------------------------------------------------------

parameters Revenue(i), VarCost(i), ProdSurp, TotProdSurp, SocWel;
parameters
Price, ConsPay, ConsSurp, TotProdSurp;

Parameter
PResults(*,*,*,*);
option
PResults:3:2:2;
Parameter
CResults(*,*);
option
CResults:3:1:1;


$macro report(name)  \
  Price = alpha-beta*q.l;                                           \
  ConsPay = Price*q.l;                                              \
  ConsSurp = alpha*q.l -0.5*beta*q.l**2 - Price*q.l;                \
  Revenue(i) = Price*
sum
(b, x.l(i,b));                              \
  VarCost(i)=
sum
(b, C(i,b)*x.l(i,b));                              \
  ProdSurp(i) = Revenue(i)-VarCost(i);                              \
  TotProdSurp =
sum
(i, ProdSurp(i));                                \
  SocWel = ConsSurp + TotProdSurp;                                  \
 
display
x.l, Revenue, VarCost, ProdSurp;                          \
 
display
q.l, Price, ConsPay, ConsSurp, TotprodSurp, SocWel;       \
                                                                    \
  Presults(name,
'x'
,i,b) = x.l(i,b);                                \
  Presults(name,
'revenue',i,'b1'
) = Revenue(i);                     \
  Presults(name,
'var.cost',i,'b1'
) = VarCost(i);                    \
  Presults(name,
'surplus',i,'b1'
) = ProdSurp(i);                    \
 
display
PResults;                                                 \
                                                                    \
  Cresults(name,
'q'
) = q.l;                                         \
  Cresults(name,
'p'
) = Price;                                       \
  Cresults(name,
'ConsPay'
) = ConsPay;                               \
  Cresults(name,
'ConsSurpl'
) = ConsSurp;                            \
  Cresults(name,
'ProdSurpl'
) = TotProdSurp;                         \
  Cresults(name,
'SocWelfare'
) = SocWel;                             \
 
display
CResults;                                                 \


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

* Maximize Social Welfare
* NLP Model
*-------------------------------------------------------------------------------

POSITIVE VARIABLES
    x(i,b)
'production by i from process b'
    q     
'demand quantity'
;
VARIABLES
    W     
'social welfare'
;

EQUATIONS
   SocialWelDef  
'social welfare definition'
   Capacity(i,b) 
'upper limit of generating output'
   MarketClear   
'market clearing'
;

SocialWelDef..   W =e= alpha*q-0.5*beta*q**2 -
sum((i,b), C(i,b)*x(i,b));
Capacity(i,b)..  x(i,b) =l= K(i,b);
MarketClear..    q =e=
sum
((i,b), x(i,b));

MODEL SocWelMax /SocialWelDef,Capacity,MarketClear/
;
SOLVE
SocWelMax maximizing W using nlp;

report(
'PriceTakers(NLP)'
)


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

* Maximize Social Welfare
* MCP Model
*-------------------------------------------------------------------------------

POSITIVE VARIABLES
    gamma(i,b) 
'dual of capacity constraint'
;
FREE VARIABLES
    lambda     
'price'
;

EQUATIONS
   PriceDef       
'price definition as inverse demand'
   Mx(i,b)        
'KKT condition for derivative w.r.t. x'
   Capacity2(i,b) 
'=g= version of capacity constraint'
;

PriceDef..        lambda - (alpha-beta*q) =g= 0;
MX(i,b)..         C(i,b) - lambda + gamma(i,b) =g= 0;
Capacity2(i,b)..  - x(i,b) + K(i,b) =g= 0;


MODEL CompetComplem /PriceDef.q, Mx.x, Capacity2.gamma, MarketClear.lambda/;
solve
CompetComplem using MCP;

report(
'PriceTakers(MCP)'
)


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

* Monopoly
* NLP Model
*-------------------------------------------------------------------------------

VARIABLES
    TotProfit    
'total profit'
;

EQUATIONS
   TotalProfit   
'Total profit calculation'
;

TotalProfit..    TotProfit =e= (alpha-beta*q)*q -
sum((i,b), C(i,b)*x(i,b));

MODEL Monopoly /TotalProfit,Capacity,MarketClear/
;
SOLVE
Monopoly maximizing TotProfit using nlp;

report(
'Monopoly'
)

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

* Oligopoly (Nash-Cournot)
* MCP Model
*-------------------------------------------------------------------------------

alias (b,bb);

EQUATIONS

   Mx2(i,b)       
'KKT condition for derivative w.r.t. x'
;

MX2(i,b)..         C(i,b) - lambda + beta*
sum(bb,x(i,bb)) + gamma(i,b) =g= 0;

MODEL Oligopoly /PriceDef.q, Mx2.x, Capacity2.gamma, MarketClear.lambda/
;
solve
Oligopoly using MCP;

report(
'Cournot(MCP)'
)

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

* Oligopoly (Nash-Cournot)
* NLP Model
*-------------------------------------------------------------------------------

VARIABLES
    CournotObj     
'Obj for Cournot model as NLP'
;

EQUATIONS
   DefCournotObj   
'Obj for Cournot model as NLP'
;

DefCournotObj..    CournotObj =e= alpha*q-0.5*beta*q**2 -
sum((i,b), C(i,b)*x(i,b))
                                  - 0.5*
sum(i, beta*sqr(sum
(b,x(i,b))));

MODEL CournotNLP /DefCournotObj,Capacity,MarketClear/
;
SOLVE
CournotNLP maximizing CournotObj using nlp;

report(
'Cournot(NLP)'
)

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

* Cournot (firm1) + Price Takers (firm2+3)
* NLP Model
*-------------------------------------------------------------------------------


set
  iCournot(i)
/i1/
;

VARIABLES
   Obj            
'Obj for Cournot/Pricetaking model as NLP'
;

EQUATIONS
   DefObj         
'Obj for Cournot/Pricetaking model as NLP'
;

DefObj..    Obj =e= alpha*q-0.5*beta*q**2 -
sum((i,b), C(i,b)*x(i,b))
                                  - 0.5*
sum(iCournot,beta*sqr(sum
(b,x(iCournot,b))));

MODEL CombNLP /DefObj,Capacity,MarketClear/
;
SOLVE
CombNLP maximizing Obj using nlp;

report(
'CombA(NLP)'
)

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

* Cournot (firm1+firm2) + Price Takers (firm3)
* NLP Model
*-------------------------------------------------------------------------------

iCournot(
'i2') = yes;

SOLVE
CombNLP maximizing Obj using nlp;

report(
'CombB(NLP)'
)

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

* Cournot (firm1) + Price Takers (firm2+3)
* MCP Model
*-------------------------------------------------------------------------------

iCournot(i) =
no;
iCournot(
'i1') = yes
;

equation

   MXComb(i,b)
'First order condition';

MXComb(i,b)..     C(i,b) - lambda + (beta*
sum
(bb,x(i,bb)))$Icournot(i) + gamma(i,b) =g= 0;


MODEL Comb1MCP /PriceDef.q, MxComb.x, Capacity2.gamma, MarketClear.lambda/
;   ;
SOLVE
Comb1MCP using mcp;


report(
'CombA(MCP)'
)

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

* Cournot (firm1+firm2) + Price Takers (firm3)
* MCP Model
*-------------------------------------------------------------------------------

iCournot(
'i2') = yes;

SOLVE
Comb1MCP using mcp;

report(
'CombB(MCP)'
)

The interesting aspect here is that we model here a market with different combinations of players: some are price takers (perfect competition) and some follow a Nash-Cournot oligopolistic behavior. Actually in some energy market models (not discussed in the book) we have coefficient (parameter) δ indicating “Market Power” where δ=0 means perfect competitive behavior and δ=1 means Cournot behavior. The values 0<δ<1 indicate some behavior in between.

Sunday, July 28, 2013

GAMSification

The book “Complementarity Modeling in Energy Markets” has a number of interesting GAMS models in the back. Some of these models are presented in scalar format, i.e. bypassing the indexing facilities. Here my attempt to express these models in a more generic format. Often that means making the models more general, which may require some extra effort and thinking upfront. Besides making the models more extensible I also hope this makes the models more structure-revealing. Besides delivering results a model should also convey the concepts on which is built.

In this model we consider a network of markets with suppliers in each market:

energynetwork

We have two type of producers: exporters and producers that only serve their own market. They optimize their profits and their optimization problem is actually very much alike.

Another player is the TCO: the Transportation System Operator. This agent also maximizes its profit.

Combining these systems of First Order Conditions (KKT conditions) with a Market Clearing condition leads to a small MCP model. Here follows my version:

 

$ontext

  
network of energy suppliers

  
Reference:
     
Steven A. Gabriel, Antonio J. Conejo, J. David Fuller,
     
Benjamin F. Hobbs, Carlos Ruiz
     
Complementarity Modeling in Energy Markets
     
Springer 2012

$offtext


Sets
   s
'producers (supplier)' /A,B,C,D/
   n
'nodes in network' /node1,node2/
   sn(s,n)
'network topology' /
     
(A,B).node1
     
(C,D).node2
  
/
   exports(s,n,n)
'export links' /(A,B).node1.node2/
   exporter(s)
'exporters are subset of suppliers'
   link(n,n)
'export flow'
;
alias(n,n1,n2);

* interesting: the next rhs are the same nut the calculate something

* different
exporter(s) =
sum(exports(s,n1,n2),1);
link(n1,n2) =
sum
(exports(s,n1,n2),1);
display
sn,exports,exporter,link;

parameters

   tau_reg(n,n) 
'regulated tariff for export (exogenous)'  /node1.node2 0.5/
   gamma(s)  
'unit production cost for each supplier'
               
/ A 10, B 12, C 15, D 18 /
   a(n)       
'coefficient for demand function'
               
/node1 20, node2 40/
   b(n)        
'coefficient for demand function'
               
/node1 1, node2 2/
   q_max(s)   
'upper bound on production (capacity)'
               
/A 10, B 10, C 5, D 5 /
   g_max(n,n)  
'max flow over link (capacity)' /node1.node2 5/
   gamma_TSO   
'unit cost for TSO' /1/
;

positive variables
   sell(s,n)   
'quantity sold'
   q(s,n)      
'quantity produced'
   f(s,n,n)    
'quantity exported'
   lambda(s,n) 
'dual on upper bound production'
   g(n,n)      
'flows'
   veps(n,n)   
'dual on transportation capacity condition'
;

free variables
   delta(s,n) 
'dual on sell=production-export'
   pi(n)      
'prices'
   tau(n,n)   
'congestion tariff (endogenous)'
;


equations
   FOC_supplier_sell(s,n)   
'derivative of Lagrangian wrt variable sell'
   FOC_supplier_q(s,n)      
'derivative of Lagrangian wrt variable q'
   FOC_supplier_exp(s,n,n)  
'derivative of Lagrangian wrt variable f'
   supplier_max_prod(s,n)   
'capacity constraint'
   supplier_supply(s,n)     
'supply is locally sold or exported'

   MC(n)    
'market clearing'

   FOC_TSO_g(n,n)           
'derivative of Lagragian wrt variable g'
   TSO_max(n,n)             
'capacity constraint'
   TSO_sum(n,n)             
'summation of flows'
;


FOC_supplier_sell(sn(s,n))..
     -pi(n)+delta(s,n) =g= 0;

FOC_supplier_q(sn(s,n))..
     gamma(s)+lambda(s,n)-delta(s,n) =g= 0;

FOC_supplier_exp(exports(s,n1,n2))..
     -pi(n2)+(tau_reg(n1,n2)+tau(n1,n2))+delta(s,n1) =g= 0;

supplier_max_prod(sn(s,n))..
     q_max(s)-q(s,n) =g= 0;

supplier_supply(sn(s,n))..
     sell(s,n)-q(s,n)+
sum(exports(s,n,n2),f(s,n,n2))=e=0;

MC(n)..
    
sum(sn(s,n),sell(s,n)) + sum
(exports(s,n1,n),f(s,n1,n)) -(a(n)-b(n)*pi(n)) =e= 0;

FOC_TSO_g(link(n1,n2))..
     -tau_reg(n1,n2) - tau(n1,n2) + gamma_TSO + veps(n1,n2) =g= 0;

TSO_max(link(n1,n2))..
     g_max(n1,n2) - g(n1,n2) =g= 0;

TSO_sum(link(n1,n2))..
     g(n1,n2) -
sum
(exports(s,n1,n2),f(s,n1,n2)) =e= 0;



MODEL compl1
/

  
FOC_supplier_sell.sell
  
FOC_supplier_q.q
  
FOC_supplier_exp.f
  
supplier_max_prod.lambda
  
supplier_supply.delta

  
MC.pi

  
FOC_TSO_g.g
  
TSO_max.veps
  
TSO_sum.tau

/;

solve
compl1 using mcp;


parameter
results(*,*,*);

$macro report(a)  \
    results(a,
'unit cost'
, s) = gamma(s);                \
    results(a,
'sell',  s) = sum
(n, sell.l(s,n));         \
    results(a,
'q',     s) = sum
(n, q.l(s,n));            \
    results(a,
'qmax'
,  s) = q_max(s);                    \
    results(a,
'export',s) = sum
((n1,n2),f.l(s,n1,n2));   \
    results(a,
'flow',  n) = sum
(n2, g.l(n,n2));          \
    results(a,
'maxflow',n) = sum
(n2, g_max(n,n2));       \
    results(a,
'price'
, n) = pi.l(n);                     \
    results(a,
'exo.tariff',n) = sum
(n2,tau_reg(n,n2));   \
    results(a,
'end.tariff',n) = sum
(n2,tau.l(n,n2));     \
   
display
results;


report(
'base case'
)


* make A,B more expensive

gamma(
'A') = 15;
gamma(
'B'
) = 17;

solve
compl1 using mcp;

report(
'A,B expensive'
)

 

Some of the results:

----    151 PARAMETER results 

                                   A           B           C           D       node1       node2

base case    .unit cost       10.000      12.000      15.000      18.000
base case    .sell             7.439       0.561       5.000
base case    .q               10.000       3.000       5.000
base case    .qmax            10.000      10.000       5.000       5.000
base case    .export           2.561       2.439
base case    .flow                                                             5.000
base case    .maxflow                                                          5.000
base case    .price                                                           12.000      15.000
base case    .exo.tariff                                                       0.500
base case    .end.tariff                                                       2.500
A,B expensive.unit cost       15.000      17.000      15.000      18.000
A,B expensive.sell             5.000                   5.000
A,B expensive.q                8.000                   5.000
A,B expensive.qmax            10.000      10.000       5.000       5.000
A,B expensive.export           3.000
A,B expensive.flow                                                             3.000
A,B expensive.maxflow                                                          5.000
A,B expensive.price                                                           15.000      16.000
A,B expensive.exo.tariff                                                       0.500
A,B expensive.end.tariff                                                       0.500