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:

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

Duopoly Example

This example is taken from the book “Complementarity Modeling in Energy Markets.”

We have two firms i producing q(i). The total production Qtot is then defined by adding up these two quantities. Because of say transmission capacity or available drilling riggs, there is a limit on the total production Qtot. The cost function of each firm is linear: Cost(i) = γ(i)*q(i). The demand price P is determined by an inverse demand function of the form P=α-β*Qtot.

The two firms each face an optimization problem: maximize profit = revenue – cost. To combine the optimization problem of both agents we form the first-order conditions (KKT conditions) and combine these in one model. For details see the book.

A “scalar” model is given in the appendix of the book. The model does not use GAMS indexing so it is merely using GAMS as a calculator. Here is a GAMSified version:

$ontext    duopoly example    Reference:       Steven A. Gabriel, Antonio J. Conejo, J. David Fuller,       Benjamin F. Hobbs, Carlos Ruiz       Complementarity Modeling in Energy Markets       Springer 2012 $offtext set    i /firm1,firm2/; parameters    gamma(i)  'unit cost' /firm1 1, firm2 2/    beta      'coefficient inverse demand function' /5/    alpha     'coefficient inverse demand function' /10/    Qmax      'upper limit on total quantity produced' / 1.06 / ; positive variables    q(i)        'production at each firm'    lambda_MAX  'dual of MaxQ constraint'    P           'price' ; free variable    Qtot  'total production: automatic positive -- q(i) is positive' ; equations    FOCfirm(i)  'first-order conditions: derivative of lagrangian wrt q(i)'    MaxQ        'max total production restriction'    QDef        'definition of Qtot (total production)'    Demand      'inverse demand function' ; FOCFirm(i)..   0 =g= P - beta*q(i) - gamma(i) - lambda_MAX; MaxQ..         Qmax =g= Qtot; QDef..         Qtot =e= sum(i, q(i)); Demand..       P =g= alpha - beta*Qtot; model problem2 /        FOCfirm.q        MaxQ.lambda_MAX        QDef.QTot        Demand.P /; solve problem2 using mcp; parameter results(*,*,*); $macro report(a) \ results(a, 'unit cost', i)       = gamma(i);      \ results(a, 'q',         i)       = q.l(i);        \ results(a, 'q',         'total') = Qtot.l;        \ results(a, 'q',         'max')   = Qmax;          \ results(a, 'p',         'total') = P.l;           \ results(a, 'dual',      'max')   = lambda_MAX.l;  \ display results; report('base case') * experiments: * make firm 2 very expensive gamma('firm2') = 6; solve problem2 using mcp; report('firm2 expensive') * make firm 2 very cheap gamma('firm2') = .5; solve problem2 using mcp; report('firm2 cheap')

The results look like:

----     75 PARAMETER results                                  firm1       firm2       total         max base case      .unit cost       1.000       2.000 base case      .q               0.630       0.430       1.060       1.060 base case      .p                                       4.700 base case      .dual                                                0.550 firm2 expensive.unit cost       1.000       6.000 firm2 expensive.q               0.900                   0.900       1.060 firm2 expensive.p                                       5.500 firm2 cheap    .unit cost       1.000       0.500 firm2 cheap    .q               0.480       0.580       1.060       1.060 firm2 cheap    .p                                       4.700 firm2 cheap    .dual                                                1.300

Notes:

• In the base case the capacity constraints is binding and we have a dual.
• In the second case the capacity constraint is no longer binding and the dual is zero. Furthermore, firm 2 is too expensive to  make it worthwhile to produce anything (the unit cost is larger than the price).
• If we make firm 2 very cheap we see that firm 1 stays in business. The capacity constraint effectively puts a floor on the price, and at this price firm 1 is still able to produce.
• In the model, both firms share a single dual λ in the FOC equation.
• This model is easily extended to more firms (no equations need to be changed), in which case we have a oligopoly.
• Qtot is a free variable and we match it to a =e= constraint. We could also make it a positive variable. In that case we could make QDef a =g= equation.

Speeding up SMAX

The following GAMS code was very slow:

ACCRLRLY(B,BA,Y,U)$SUM((H,T,UR), XSB7(B,BA,"NP",H,Y,T,U,UR))                              = ACRLRLY(BA,Y, U) * SMAX((H,T,UR), XSB7(B,BA,"NP",H,Y,T,U,UR));

Note that the $ condition on the left is correct: it makes sure we don’t take the SMAX of an empty set. This is really needed as SMAX over an empty set is unreliable (see:http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/gams-smin-bug.html).

However the SMAX will not skip zeros like a SUM, making this operation very slow. A possible fix is:

ACCRLRLY(B,BA,Y,U)$SUM((H,T,UR), XSB7(B,BA,"NP",H,Y,T,U,UR))     = ACRLRLY(BA,Y, U) * SMAX((H,T,UR)$XSB7(B,BA,"NP",H,Y,T,U,UR), XSB7(B,BA,"NP",H,Y,T,U,UR));

The difference in timing is very large:

Original	75 seconds
New version	0.9 seconds

This only works if it is safe to skip zeros. I.e. if there are positive elements in each smax. We can check for this e.g. by using an abort$(cond) statement.

QP is convex or not?

I have a QP problem that is rejected (not convex) or solved depending on the solver or method:

Solver	Objective
Gurobi	1.01022764e+02 (23 iterations, 0.16 seconds)
Cplex	*** CPLEX Error  5002: Q in %s is not positive semi-definite.
MOSEK (QP)	Return code - 1295 [MSK_RES_ERR_OBJ_Q_NOT_PSD]: The quadratic coefficient matrix in the objective is not PSD.
MOSEK (NLP)	1.0102276353e+002 (24 iterations, 0.30 seconds)
XPRESS	1.0102297e+002 (12 iterations, 1 second)
KNITRO	1.01022803552155e+002 (21 iterations, 0.5 seconds)
IPOPT (MUMPS)	Out of memory in MUMPS.
IPOPT (MA27)	1.0102276352077803e+002 (40 iterations, 0.7 seconds)

Try this to explain in a coherent fashion to a client…

The IPOPT issues may be unrelated to the convexity issue. Replacing the linear solver MUMPS by Harwell’s MA27 will often speed things up but from this we also see a more reliable behavior.

More background on this behavior: the Q matrix is diagonal with one negative element, but that variable is constraint to be zero. The left model in the table below illustrates this. The right model is the same, but now the variable is fixed to zero. 

GAMS model	variables z,x1,x2; equations obj,e; obj.. z =e= sqr(x1-1) - sqr(x2-2); e..   x2 =e= 0; model m/all/; solve m minimizing z using qcp;	variables z,x1,x2; equations obj; obj.. z =e= sqr(x1-1) - sqr(x2-2); x2.fx=0; model m/all/; solve m minimizing z using qcp;
Cplex	*** CPLEX Error  5002: Q in %s is not positive semi-definite.	*** CPLEX Error  5002: Q in %s is not positive semi-definite.
Gurobi	OK	*** Objective Q not PSD (negative diagonal entry)
Xpress	OK	?899 Warning: The quadratic objective is not convex
MOSEK	Return code - 1295 [MSK_RES_ERR_OBJ_Q_NOT_PSD]: The quadratic coefficient matrix in the objective is not PSD.	Return code - 1295 [MSK_RES_ERR_OBJ_Q_NOT_PSD]: The quadratic coefficient matrix in the objective is not PSD.
Cplex LP file	\ENCODING=ISO-8859-1 \Problem name: gamsmodel Minimize _obj: z - 2 x1 + 4 x2 + [ 2 x1 ^2 - 2 x2 ^2 ] / 2 Subject To obj#0: z  = -3 _e#1:  x2  = 0 Bounds       z Free       x1 Free       x2 Free End	\ENCODING=ISO-8859-1 \Problem name: gamsmodel Minimize _obj: z - 2 x1 + 4 x2 + [ 2 x1 ^2 - 2 x2 ^2 ] / 2 Subject To obj: z  = -3 Bounds       z Free       x1 Free       x2 = 0 End
Gurobi LP file	Minimize   - 2 x1 + 4 x2 - 3 Constant + [ 2 x1 ^ 2 - 2 x2 ^ 2 ] / 2 Subject To e: x2 = 0 Bounds x1 free x2 free Constant = 1 End	Minimize   - 2 x1 + 4 x2 - 3 Constant + [ 2 x1 ^ 2 - 2 x2 ^ 2 ] / 2 Subject To Bounds x1 free x2 = 0 Constant = 1 End

The LP files are reproduced here to show what Gurobi and Cplex think they are solving (these files were generated with a solver option). I am not sure why some solvers accept the version with an explicit equation x2=0 while reject the version with x2=0 as bounds. It is all a mystery to me.

The client solved this model as an NLP using MINOS. It is interesting to see how many issues we face just by solving it as a QP!

Aggregation: database vs GAMS

The following aggregation code in GAMS aggregates PRODUCTION and AREA HARVESTED from a relatively large data set from FAO (2.2 million records). Afterwards we recalculate YIELDS, as aggregating yields can not be done in the same way as PRODUCTION and AREA (this also famously occurs when aggregating prices: first do volume and value and then recalculate prices).

*------------------------------------------- * Aggregation *------------------------------------------- set    y(year) /1990*2000/    aggr_type(type) 'these quantities are aggregated directly' /           'Area Harvested'             Ha           'Production Quantity'        Tonnes     / ; parameter AggrCropData(cropgroup,region,y,type); AggrCropData(cropgroup,region,y,aggr_type) =      sum((cropmap(crop,cropgroup),regionmap(country, region),unit,flag),              cropdata(country,crop,aggr_type,y,unit,flag)); *------------------------------------------- * Recalculate Yield (Hg/Ha) *------------------------------------------- AggrCropData(cropgroup,region,y,'Yield')$AggrCropData(cropgroup,region,y,'Area Harvested') =     1e5*AggrCropData(cropgroup,region,y,'Production Quantity') /       AggrCropData(cropgroup,region,y,'Area Harvested')

See also: http://yetanothermathprogrammingconsultant.blogspot.com/2013/05/large-scale-aggregation-example.html.

Of course we can also do this in SQL:

--- --- Aggregation --- use aggregation;   --- --- if target table exists, drop it --- IF EXISTS (SELECT * FROM sys.objects  WHERE object_id = OBJECT_ID('[dbo].[AggregatedData]') AND type='U') DROP TABLE AggregatedData;   --- --- Step 1: --- Aggregate Area and Production ---   SELECT C.CropGroup, R.Region, A.Element, A.[Year], Sum(A.Value) AS [Value] INTO AggregatedData FROM Production_Crops_E_All_Data AS A,      cropmap AS C,      regionmap AS R WHERE      A.Element In ('Area Harvested','Production Quantity') And      A.[Year]>='1990' And A.[Year]<='2000' And      A.Country=R.Country And      A.Item=C.Crop GROUP BY C.CropGroup, R.Region, A.Element, A.[Year];     --- --- Step 2: --- Recalculate Yield --- INSERT INTO AggregatedData ( CropGroup, Region, Element, [Year], [Value] ) SELECT A.CropGroup, A.Region, 'Yield', A.[Year], 1.0e5*A.[Value]/B.[Value] FROM AggregatedData AS A,      AggregatedData AS B WHERE      A.CropGroup = B.CropGroup AND      A.Region = B.Region AND      A.[Year] = B.[Year] AND      A.Element = 'Production Quantity' AND      B.Element = 'Area Harvested' AND      B.Value > 0;

Interesting to see how well SQL server is doing compared to GAMS. Here are timings in seconds:

	Step1	Step2
GAMS	1.4	0.08
MS Access	5.4	0.04
SQL Server	1.4	0.06