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!

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

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:

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:
  
  This is due to the way we had to declare the sets for use with the projection option.

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).

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.

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:

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

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));

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.

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:

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

Of 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 x₁,x₂ and x₃ 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.

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.