Network formulation

I was looking at a model containing underlying traffic network. Such a network can be quite elegantly formulated in GAMS, but apparently many users still have problems with this. The mathematical formulation is:

The GAMS formulation can look like:

$ontext   max flow network example   Data from example in      Mitsuo Gen, Runwei Cheng, Lin Lin      Network Models and Optimization: Multiobjective Genetic Algorithm Approach      Springer, 2008   Erwin Kalvelagen,   Amsterdam Optimization,   May 2008 $offtext sets   i 'nodes' /node1*node11/   source(i) /node1/   sink(i) /node11/ ; alias(i,j); abort$(card(source)<>1) "We want one source node"; abort$(card(sink)<>1) "We want one sink node"; parameter capacity(i,j) /    node1.node2 60,   node1.node3 60,   node1.node4 60,   node2.node3 30    node2.node5 40,   node2.node6 30,   node3.node4 30,   node3.node6 50    node3.node7 30,   node4.node7 40,   node5.node8 60,   node6.node5 20    node6.node8 30,   node6.node9 40,   node6.node10 30,  node7.node6 20    node7.node10 40,  node8.node9 30,   node8.node11 60,  node9.node10 30    node9.node11 50,  node10.node11 50 /; set arcs(i,j); arcs(i,j)$capacity(i,j) = yes; display arcs; parameter rhs(i); rhs(source) = -1; rhs(sink) = 1; variables   x(i,j) 'flow along arcs'   f      'total flow' ; positive variables x; x.up(i,j) = capacity(i,j); equations    flowbal(i) 'flow balance' ; flowbal(i)..    sum(arcs(j,i), x(j,i)) - sum(arcs(i,j), x(i,j)) =e= f*rhs(i); model m/flowbal/; solve m maximizing f using lp;

MIP Gap not decreasing

Here is a small fragment of a Cplex log:

Nodes                                         Cuts/    Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap    3186  2520     -666.6217  4096      956.6330     -667.2010  1313338  169.74%    3226  2560     -666.6205  4097      956.6330     -667.2010  1323797  169.74%    3266  2600     -666.6201  4095      956.6330     -667.2010  1335602  169.74% Elapsed real time = 2801.61 sec. (tree size = 77.54 MB, solutions = 2) *  3324+ 2656                         -125.5775     -667.2010  1363079  431.31%    3334  2668     -666.5811  4052     -125.5775     -667.2010  1370748  431.31%    3380  2714     -666.5799  4017     -125.5775     -667.2010  1388391  431.31%    3422  2756     -666.5791  4011     -125.5775     -667.2010  1403440  431.31%

Usually we expect the gap to decrease while the MIP solver is working its way. Here we see a different behavior. Because we go through zero the gap% is actually increasing here.

I believe the gap calculation used by Cplex is:

Obviously this is sensitive w.r.t. to the sign of BestInteger. Of course in practice we don’t see this behavior very often (in many models BestInteger and BestNode have the same sign).

Writing solutions back into a database

The issue of writing optimization results back into a database is often a little bit more complex than reading input data (see http://yetanothermathprogrammingconsultant.blogspot.com/2011/04/live-databases.html for some thoughts on reading data).

First with writing more things can go wrong: you need additional permissions, you can overwrite data, and conceptual there are some issues. Here are some questions and possibilities that come to mind:

1. Do we have a say in how the tables look like or do we need to follow some existing data model?
2. Update or Insert. Do we use an SQL INSERT or SQL UPDATE to populate solution tables?
3. Some databases have SQL extensions that implement an UPSERT. E.g. MySQL has INSERT INTO … ON DUPLICATE KEY UPDATE, or REPLACE INTO, Oracle and SQL Server have MERGE.
4. It may be better to erase the solution tables first.
5. Or always INSERT with a new key indicating the run id. A new solution gets a new run id.
6. For very large data some RDBMS systems have a BULK INSERT facility (SQL Server: BCP or BULK INSERT, SQL*Loader for Oracle).
7. We can write to temporary tables and deal with the problem of merging results into the real tables further down the road.
8. Just write a well-designed temp database (e.g. in Access) and let the IT people deal with getting the data into their systems.

In general I prefer 8, let others take responsibility for plugging things into the database. Note that this step can be expensive. A few weeks ago I was discussing a complex model where we needed to get the results into the RDBMS. To achieve a reasonable turnaround time, I was given a time limit of about an hour to come up with solutions for the model, while the system guy allocated 2 hours for himself to put the solution back into the database (much time is spent on checks).

Remote Desktop adventure

Suddenly my remote desktop window went black. Even after reconnecting, I got again just a black screen. Took me a little while to figure out that with Ctrl-Alt-End I could send a Ctrl-Alt-Del to the remote session, so I could log out. That solved the problem: after logging on again everything was fine.

GAMS/Gurobi and Ctrl-C

I am seeing some strange results when I interrupt a GAMS/Gurobi run with Ctrl-c (in the IDE using the Interrupt button). The GAMS link should give in that case the best integer solution found so far, basically identical to an iteration or time limit. However the solution seems wrong: many of the binary variables x(i,j) in a large model are off (they are reported as zero instead of one). I suspect this is a problem in the GAMS link (i.e. not in Gurobi itself). From what I can see the link it passing back levels of integer variables incorrectly in these circumstances (I have seen this before: integer variables should have basis status SUPERBASIC to prevent only a basis status being sent back).

Btw, for this (huge) model the Gurobi heuristics do a fantastic job finding good solutions:

Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time      0     0 179791.544    0 6744 1511067.07 179791.544  88.1%     -   91s H    0     0                    247836.12285 179791.544  27.5%     -   93s H    0     0                    232489.21120 179791.544  22.7%     -   96s      0     0 181272.276    0 7366 232489.211 181272.276  22.0%     -  397s H    0     0                    226736.43968 181272.276  20.1%     -  399s H    0     0                    218612.71830 181272.276  17.1%     -  404s

Just need to remember not to use the Interrupt button.

PS. Another minor issue with the link:

MIP  Solution:        7259.007800    (-2147483648 iterations, 0 nodes) Best possible:        6916.876514 Absolute gap:          342.131286 Relative gap:            0.047132

it has troubles with the iteration count (which actually is 120570, a number that should not overflow).

Optimal Spread (2)

Is turns out we also have to distribute x(t)’s such that we follow a non-uniform parameter. In the previous section http://yetanothermathprogrammingconsultant.blogspot.com/2011/05/optimal-spread.html we used a function v(t) = 1 with cumulative function w(t) = t.

In this case we consider a generic parameter v(t), with its cumulative cousin w(t)=w(t-1)+v(t). To trace this function as closely as possible we can use:

E.g. when we use v(t)=t, then we see:

----     77 PARAMETER result                v           x i1        1.000 i2        2.000 i3        3.000 i4        4.000 i5        5.000       1.000 i6        6.000 i7        7.000 i8        8.000 i9        9.000       1.000 i10      10.000 i11      11.000 i12      12.000       1.000 i13      13.000 i14      14.000       1.000 i15      15.000 i16      16.000       1.000 i17      17.000 i18      18.000       1.000 i19      19.000 i20      20.000       1.000

Optimal spread

In a multi-objective MIP model I am working I have a binary variable x(t). This variable is mostly zero but now and then it is one. One of the goals is to achieve an even spread of the ones. I.e.

The first row has the ones clustered in the middle. The second one looks pretty good.

The following idea is based on how Q-Q plots are used to find distributions in statistical data (http://en.wikipedia.org/wiki/Q-Q_plot).

First we form a running sum:

 

The “predicted’ running sum for an evenly distributed x(t) can be expressed as tn/T where n is the total number of ones, t is the current index and T is the last t:

 

Now we try to keep the values s(t) as close as possible to the predicted values tn/T. I.e.:

Note that n is a variable in my model: I don’t know in advance how many x(t)’s are turned on. Indeed when I run this (with n=3) I see:

----     42 VARIABLE x.L  t3  1.000,    t7  1.000,    t10 1.000

Here are the results when we compare the two solutions (click to enlarge):

Indeed this seems to confirm this approach could work, and we can add d(max)-d(min) as an extra objective to minimize.

Excel: Trace Dependents

I often use the “Trace Dependents” tool in Excel to quickly see which formulas use a given cell. Here is an example where I was a bit overwhelmed:

Large (almost) triangular spreadsheet?

A spreadsheet model without circular references can be called “triangular”. Excel can form a dependency graph and ordering such that recalculation can be done in one iteration. If there are circular references the spreadsheet needs to perform more iterations before the results converge.

I was given a very large spreadsheet implementing an agricultural country trade model. The question was how near triangular is this? With the Excel tool to find circular references, we only find a single case:

Luckily I have a tool that takes an Excel spreadsheet and parses Excel formulas and produces a GAMS representation of this. Essentially each cell forms an equation:

The whole spreadsheet can be viewed as a fixed point expression:

This is of course just a system of nonlinear equations:

When we solve this spreadsheet model as an NLP with Conopt we get some statistics:

---   67,618 rows  68,707 columns  232,364 non-zeroes ---   776,017 nl-code  103,528 nl-non-zeroes    Iter Phase Ninf   Infeasibility   RGmax    NSB   Step InItr MX OK       0   0        7.3221660778E-07 (Input point)                                 Pre-triangular equations:    22507                                 Post-triangular equations:   13897       1   0        1.2425686171E-07 (After pre-processing)       2   0        3.4325009848E-11 (After scaling)

Graphically, after reordering rows and columns, we have:

This certainly gives an indication we are not close to a triangular model and we need a simultaneous equation solver to handle this type of model.

Initial Point for NLP models

This example was introduced to illustrate the need for good initial points in NLP models (see http://amsterdamoptimization.com/pdf/snopt.pdf).

The problem is just a three variable model: find the smallest enclosing circle around n points. GAMS has a default initial point of zero which is often very bad. Only IPOPT can find an optimal solution quickly. But if we add a good starting point, all solvers find the optimal very quickly. Even the performance of IPOPT improves.

For a presentation I cleaned up the example a little bit:

 

Note: actually the model is easier to solve using a slightly different formulation:

(i.e. minimize the square of the radius r).