Thursday, December 5, 2013

Populating a spreadsheet through VBA

Populating a large spreadsheet table can be slow if we push data cell-by-cell. Here is a small example.

Naïve approach

We start with creating a 1000x1000 table by directly accessing individual cells.

Private Sub run1()
  
Dim r As
Range
  
Set r = Range("Sheet1!B2"
)

   r.Resize(m +
1, n + 1
).ClearContents

  
Dim t0 As
Single
   t0 = Timer

  
Dim i As
Integer
  
Dim j As
Integer

  
For i = 1 To
m
      r(
1 + i, 1) = "row"
& i
  
Next
i
  
For j = 1 To
n
      r(
1, 1 + j) = "col"
& j
  
Next
j

  
For i = 1 To
m
     
For j = 1 To
n
         r(
1 + i, 1 + j) = 3.14

     
Next j
  
Next
i

  
Dim t1 As
Single
   t1 = Timer
   r(
0, 1) = "Time : " & Format(t1 - t0, "fixed"
)
End
Sub

This approach takes about 240 seconds on my machine.

A refinement is to turn off screen updating: Application.ScreenUpdating = False. This did not make a difference for me. Excel seems to turn off screen updating automatically after a little while.

Fast approach

A better approach is to create a VBA array and to write this array in one swoop to the spreadsheet. We use here a Variant array as it contains both strings and numbers. The code looks like:

Private Sub run3()
  
Dim r As
Range
  
Set r = Range("Sheet3!B2"
)

   r.Resize(m +
1, n + 1
).ClearContents

  
Dim t0 As
Single
   t0 = Timer

  
Dim i As
Integer
  
Dim j As
Integer

  
Dim a(1 To m + 1, 1 To n + 1
)
  
For i = 1 To
m
      a(
1 + i, 1) = "row"
& i
  
Next
i
  
For j = 1 To
n
      a(
1, 1 + j) = "col"
& j
  
Next
j

  
For i = 1 To
m
     
For j = 1 To
n
         a(
1 + i, 1 + j) = 3.14

     
Next j
  
Next
i

   r.Resize(m +
1, n + 1
).Value = a

  
Dim t1 As
Single
   t1 = Timer
   r(
0, 1) = "Time : " & Format(t1 - t0, "fixed"
)

End Sub

This code takes 0.46 seconds!

image

This strategy is also useful when calling Excel from other languages such as Delphi or C#.

Friday, November 15, 2013

Converting USDA PSD data to GAMS GDX

USDA has some good agricultural production and trade data available at http://www.fas.usda.gov/psdonline/psdHome.aspx. Here we try to read the CSV file psd_alldata.csv and convert it to a usable GAMS GDX file.

$ontext

  
Convert psd_alldata.csv to a GDX file

  
Source:
  
http://www.fas.usda.gov/psdonline/psdDownload.aspx

  
Sorry for the obscure syntax, I am sure this will confuse users

  
Erwin Kalvelagen, 2013

$offtext


$set csv psd_alldata.csv
$set gdx psd_alldata.gdx

$version 241

$call csv2gdx %csv% id=alldata useheader=y index=2,4,5,9,11 value=12


set year 'superset' /1900*2050/;
alias
(commodity,country,attribute,unit,*);
parameter
alldata(commodity,country,year,attribute,unit);

$gdxin %gdx%
$load alldata


sets

  commodity2(commodity)
  country2(country)
  year2(year)   
'market year'
  attribute2(attribute)
  unit2(unit)
;



option
  commodity2 < alldata
  country2 < alldata
  year2 < alldata
  attribute2 < alldata
  unit2 < alldata
;


execute_unload '%gdx%', alldata, commodity2=commodity, country2=country,
                        year2=year, attribute2=attribute, unit2=unit;

The final GDX file looks like:

image

This file can be read as:

$ontext

  
Example: Read PSD data

$offtext


$set gdx psd_alldata.gdx


sets
  commodity
  country
  year       
'market year'
  attribute
  unit
;


parameter alldata(commodity,country,year,attribute,unit);

$gdxin  %gdx%
$load   commodity country year attribute unit
$loaddc alldata

Notes:

  • We use the tool csv2gdx to read the csv and produce a raw gdx file (the gdx file will be overwritten with the final version after adding the supporting sets).
  • The projection option is used to extract the sets from the parameter alldata. This is non-intuitive (i.e. ugly) syntax but it is fast.
  • The set year is explicitly declared to make sure the years are ordered. In the csv file the (marketing) years are not ordered. Refinement: in order to make sure the superset of years is not chosen too small, we probably should have used a $loaddc alldata instead of $load alldata.
  • The declaration of the sets look a bit funny in the gdx viewer:
    image
    This is due to the way we had to declare the sets for use with the projection option.

Monday, October 28, 2013

Turbo Pascal Compiler in JavaScript: http://www.teamten.com/lawrence/projects/turbo_pascal_compiler/.


There must be more of these interesting efforts. I remember a Basic interpreter in TeX (http://texcatalogue.sarovar.org/entries/basix.html).

Wednesday, October 23, 2013

Paul Krugmans review of the new book by Bill Nordhaus:  www.nybooks.com

(More info about book see: www.amazon.com)



Of course it is useful to keep in mind this criticism: pindyck.html.

Sunday, October 13, 2013

gdx2sqlite

SQLITE (http://www.sqlite.org/) seems to be a reasonable alternative to MS Access as a file based database system. Basically we want to use a database rather than a spreadsheet to disseminate model results (the results are too big to comfortably fit in a spreadsheet). One of the advantages of Sqlite is that it can handle files that exceed 2 gb. In the GAMS tool GDX2ACCESS I tried to be smart and do the bulk inserts as  fast as I could (often that means using a csv file in between). In Sqlite there are no special facilities to do bulk inserts (but of course we can be judicious with using transactions). Here are comparisons with some test data:

image

Even though I use just INSERT statements, SQLITE is pretty fast, and produces smaller files than MS Access.

image

Sunday, September 29, 2013

Ill-posed problem

A client was solving a (large, non-linear) problem and got as message back:

ill-posed problem, solution status is unknown

The manual at http://www.gams.com/dd/docs/solvers/mosek.pdf did not give any help. This message seems to indicate the modeler (me) did something wrong. However I am really sure the model is just fine: smooth and differentiable everywhere and even convex.

Solver messages are really important (they’ll help the user when they are in trouble), and I think many of these messages are not very good; they often seem to be made up after a long day of coding without too much thought. A good message would probably try to answer the obvious question: “And what now?”. If at all possible it should help the user fix the problem (may be even adding some information about which equations and variables are involved).

Instead of fixing the so-called ill-posed problem we used a newer version of the solver, and now we obtained an optimal solution.

Sometimes it can help to have different solvers available. In GAMS you can then automate:

option nlp=minos;
solve
mod1 minimizing z using nlp;
if
(mod1.modelstat > 2,
*   not optimal or locally optimal

  
option nlp=conopt;
  
solve
mod1 minimizing z using nlp;
);

Complex loop vs complex sum in GAMS

My original formulation was a sum. Ugly and complex but it works; unfortunately it took a long time 1616 seconds. After a few “grande dark roasts” I rewrote it as a loop. Much better: 0 seconds! I don’t understand exactly why.

* fast: 0 seconds
loop((TechYieldname,fpumap(fpuall,fpu),cropmap(cropAll,crop),cropjmap(crop,j),LandMap(LandAll,Land),tech,pt(p,t))$UserTechYieldDriver(TechYieldname,FpuAll,CropAll,LandAll,Tech,p),
   YldTechGRX00a(TechYieldname,j,fpu,land,tech,t) =
        UserTechYieldDriver(TechYieldname,FpuAll,CropAll,LandAll,Tech,p);
);


* slow: 1616 seconds

YldTechGRX00b(TechYieldname,j,fpu,land,tech,t) =
sum((fpumap(fpuall,fpu),cropmap(cropAll,crop),cropjmap(crop,j),LandMap(LandAll,Land),pt(p,t)),UserTechYieldDriver(TechYieldname,FpuAll,CropAll,LandAll,Tech,p));

Thursday, August 29, 2013

Pindyck

Robert Pindyck in “CLIMATE CHANGE POLICY: WHAT DO THE MODELS TELL US?”:

ABSTRACT: Very little. ….

One of the problems with these long range models is the use a discount rate. This rate can have large effects on the model results. This is true in many models with long planning horizons. One should be aware of this, e.g. by running different scenarios with different discount rates to see what the impact is on the solutions.

Listen also to: http://www.econtalk.org/archives/2013/08/pindyck_on_clim.html.

Tuesday, August 27, 2013

Cplex solution pool

I am working on a fairly complex but small MIP model, where we want to find a large number of non-dominated, feasible solutions (Pareto solutions). In one of the models we have developed for this, we use Cplex to quickly give us all solutions in a certain neighborhood using the solution pool facilities (we filter out dominated solutions later on). As the models are fairly small, it was not immediately obvious if more cores or cpus would help. It turns out it really does:

image

2D Interpolation With SOS2 variables

Using SOS2 variables to implement a 1D interpolation scheme is fairly easy (see: http://yetanothermathprogrammingconsultant.blogspot.com/2009/06/gams-piecewise-linear-functions-with.html). However, a 2D problem is already much more difficult. Here is an example taken from the Lindo web site:

$ontext

 
2D Interpolation with SOS2 variables

 
See: http://www.lindo.com/cgi-bin/modelf.cgi?Piecelin2Dsos3.txt;LINGO

$offtext


$set  n1 3
$set  n2 3
$eval n3 (%n1%+%n2%-1)


sets
   i
/i1*i%n1%/
   j
/j1*j%n2%/
   k
/k1*k%n3%/
;

table data(i,j,*)
          
x   y   f

i1.j1    195 1800  20
i1.j2    217 1900  26
i1.j3    240 2000  30
i2.j1    195 3500  52
i2.j2    217 3600  61
i2.j3    240 4100  78
i3.j1    195 5100  69
i3.j2    217 5200  80
i3.j3    240 5600  93
;

parameters
   xv(i,j)
   yv(i,j)
   fv(i,j)
;
xv(i,j) = data(i,j,
'x');
yv(i,j) = data(i,j,
'y'
);
fv(i,j) = data(i,j,
'f'
);

sos2 variables

  wx(i)
  wy(j)
  wd(k)
;



positive variables
  WGT(i,j)
  xa
  ya
  fa
;


variables z;

equations

   xconvex
   yconvex
   dconvex
   ewx
   ewy
   ewd

   compx
   compy
   compfv

   obj
;

xconvex..
sum(i, wx(i)) =e= 1;
yconvex..
sum
(j, wy(j)) =e= 1;
dconvex..
sum
(k, wd(k)) =e= 1;

ewx(i)..  wx(i) =e=
sum
(j, wgt(i,j));
ewy(j)..  wy(j) =e=
sum
(i, wgt(i,j));
ewd(k)..  wd(k) =e=
sum((i,j)$(ord(i)+ord(j)-1=ord
(k)), wgt(i,j));

compx.. xa =e=
sum
((i,j), xv(i,j)*wgt(i,j));
compy.. ya =e=
sum
((i,j), yv(i,j)*wgt(i,j));
compfv.. fa =e=
sum
((i,j), fv(i,j)*wgt(i,j));


obj..  z =e= YA + 15*XA;

fa.lo = 67;
xa.lo = 227;
xa.up = 229;



model m /all/
;
solve
m minimizing z using mip;

4D Interpolation

For a engineering design application we were offered a 4D lookup table that predicts

imageOf course rather than “lookup” it is often better to do interpolation, so values in between are accepted. This still looks pretty bad. Our model was a MIP and using SOS variables we can easily do 1D and with some effort do 2D. For more dimensions, things are really difficult.

Luckily it turned out that we know the values of x1,x2 and x3 before we start the optimization model. I.e. we can do a 3D Interpolation in GAMS, followed by a 1D interpolation with SOS2 variables in Cplex.

For the 3D interpolation we used a simple IDW (Inverse Distance Weighting) scheme (explained in: http://en.wikipedia.org/wiki/Inverse_distance_weighting). We do this a number of times so we can estimate a good set of control points for the remaining 1D interpolation. This was done by a piecewise linear function, as explained here: http://yetanothermathprogrammingconsultant.blogspot.com/2009/06/gams-piecewise-linear-functions-with.html.

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.