Thursday, March 12, 2009

GAMS: how to get a list of parameters in a GDX file

> How can I get a list of the parameters in the farm.gdx? just the parameters, only the names

$call gdxdump farm.gdx symbols noheader | awk "{ if ($4==\"Par\") print $2 }" > farm_parameters.txt 

This is using a Unix tool awk, but this also works under windows as GAMS provides a copy of a number of such tools.

Wednesday, March 11, 2009

How to linearize this function

I'm trying to solve an optimization problem via linear or mixed integer programming model. But, the objective function includes the following expression f(x)=1/x + 1/x^2 Please let me know how to formulate this function as linear or mixed integer programming model

There are several possibilities:

  1. Solve using a piecewise linear formulation (e.g. using SOS2 variables).
  2. Solve using as SLP (Successive LP). In this case for simplicity you may want to try just linearizing (Taylor approximation) the objective each cycle.
  3. Use an NLP solver (e.g. MINOS, SNOPT or CONOPT). You could solve first as LP (eg by dropping this term) and use that solution as starting point.

Monday, March 9, 2009

Sparse data in MSF OML

The modeling language of the Microsoft Solver Foundation does not support sparse data. I.e. if you have a parameter a[i,j] you need to have values for all a[i,j]’s and you cannot skip the zero’s. In some cases this is a little bit of problem. Adding all zero records may be a lot of work and may take a large amount of space in your spreadsheet. Sometimes it is possible to work around this problem by adding extra sets and columns. Probably the easiest way to explain this is with a small example. A network optimization problem is a good example where a sparse data structure can help. In this case we have a max flow model with a small but sparse network. In GAMS we can write:

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

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;


The parameter capacity defines the network: where there is no capacity, we assume there no arc. We can calculate two-dimensional set arcs(i,j) to indicate which arcs exist. This can be used to simplify the flow-balance equations. Note that in GAMS only the variables that appear in the equations are actually part of the model. I.e. although we declare x(i,j), in the equations we only have x(arcs). This means we only have 22 variables x(i,j).



In OML we cannot use a sparse table for the capacities: we would need a complete table with 11x11=121 records. That would also mean that we have 121 variables. The trick is to have an explicit set arcs in the data table. As OML does not have 2-dimensional sets, I used parameters From[.], To[.] to describe the network. The complete data is below:

maxflowdata 
The OML model can now look like:


maxflowmodel 
The node balance equation becomes a little bit unwieldy and less readable with this approach. This is a small, artificial example but note that quite a few models have some network component. It illustrates that leaving out sparse data handling and multidimensional sets to simplify a modeling language, although not a show stopper, comes with a cost in the form of additional complexity for the modeler. Note that in this model we still have a possibly large node table to maintain (the GAMS model handles this sparse).

AMPL vs. GAMS on radiotherapy treatment planning

 

hi all,
i would just like to solicit your advice regarding whether to get an
academic license for AMPL or GAMS. i'm currently doing a master's
thesis on cancer radiotherapy treatment planning (RTP), and i've yet
to decide between AMPL or GAMS for prototyping the models. some of my
decision criteria include language syntax, user interface, software
cost and available literature on RTP.
i think that the syntax for AMPL is better compared to GAMS because
the former more closely resembles the algebraic model formulations.
however, there seem to be more literature on RTP in GAMS and there's
also an integrated development environment for it. thanks for your
help.

 

Both systems are very good in what they are doing and widely used, so you cannot really go wrong with either choice. I would probably suggest to add extra points for the modeling system that is used by your colleagues and collaborators. That makes exchanging models and data easier and also is easier when discussing problems, tricks, issues etc.

Bob Fourer (AMPL) answered:

It's hard to find someone who can give equally expert advice on two competing
systems, as once you become familiar with one of them you don't usually have
much incentive to keep learning about the other.  But here are a few comments
from my hardly unbiased view.

AMPL was designed with the idea of being much closer to mathematical notation
and generally much more natural to use than GAMS, and it's superior on that
score.  A GAMS model typically relies on more special conventions and
reformulations than its AMPL counterpart; a case in point is the often extensive
use of the GAMS $ operator to impose various conditions.  Also, the IDE
notwithstanding, GAMS is fundamentally more of a batch system whereas AMPL
offers a more flexible option of interactively exploring models and results.
Finally while in certain areas GAMS is established through long use, still I see
modelers in these areas choosing AMPL, particularly when they are undertaking
new projects that do not depend on existing GAMS models.

In my opinion AMPL and GAMS are closer in practice than suggested here (e.g. where you use $ in GAMS, one would use : in AMPL). I actually slightly prefer the GAMS syntax when doing real work, as it is a little bit more compact and it is obvious where a summation ends (in AMPL this is based on operator priority, in GAMS a sum is visually bracketed by parentheses).

Friday, March 6, 2009

A conference scheduling problem (2)

In this posting I described the solution for a conference scheduling problem. Now the question is posed: can we add a constraint on the room capacity: up to 15 people per room. The first thing I did was to inspect my solution. Indeed we are exceeding this limit here and there:
schedule2 
In row total some of the occupancy counts are larger than 15. Note that this was a solution with objective=1491. We can add the room capacity constraint as follows:

roomcap 
Using this additional constraint, the solution looks like:
schedule3
This reduces the optimal objective from 1491 to 1487. I.e. we pay for this with a small decrease w.r.t. the preferences. The solution time was 2 minutes with GAMS/Cplex.

I use here the fact that all rooms are identical. If there are five rooms with different capacity, we need to know how talks are assigned to individual rooms. That requires extra binary variables (either by adding a room index to variable x or by introducing a new variable that handles the room assignment).

A conference scheduling problem

 

From: GagoX <elgrang...@gmail.com>
Date: Mar 3, 11:43 am
Subject: scheduling problem optimizing preferences
To: sci.math

Hi everyone,

this is a real life problem and I have no idea even how to start!

I need to organize a small conference:

Np=50 participants (P1 to P50) and Nt=13 talks (T1 to T13).
There are 5 time spots each with 5 rooms available.
T1 will be given 3 times, T2 to T7 will be given 2 times each and T8
to T13 will be given only once. Consequently there are 22 talks. The
constraint is that the repeated talks should not be scheduled in the
same time slot.

Each participant has given his preference to which talk to attend,
from 1st preference to 8th preference. And each participant should
attend 5 talks, one in each time slot.

Problem: What talks  should be scheduled in each time slot so to
optimize the preferences?

The optimum solution would be to fulfill that each participand attend
their 5 first preferences.

How is this problem called and is it possible to solve...what should
be the strategy.Or even better, does someone know the solution?

I guess if I knew how to ennumerate  each possibility (I do not know),
I could then quantify the fulfilling of the preferences and then take
the best one (brute force) (of course programming the whole
thing)...but I am afraid that there are too many possibilities!

Thaks a lot in advance!!!!

This is a task that is sometimes very well suited to be implemented as mixed integer programming problem. (If they become too big, scheduling problems can be better implemented using different techniques.)

Lets try to do this quickly in GAMS. We start with the sets, and the basic data:

set
p 'participants' /p1*p50/
t 'presenter' /t1*t13/
n 'talk number' /n1*n3/
s 'timeslot' /slot1*slot5/

;

scalar rooms 'number of parallel sessions' /5/

set talk(t,n) 'each presenter can give several talks' /
t1.(n1,n2,n3)
(t2*t7).(n1,n2)
(t8*t13).n1
/;

scalar numtalk 'total number of talks';
numtalk=card(talk);
display talk,numtalk;


** actually the count is 21
** abort$(numtalk<>22) "oops, check your sets";

parameter ntalk(t) 'numbers of times a talk is given';
ntalk(t) = sum(talk(t,n),1);
display ntalk;


We see that the total number of talks is actually 21.



The preferences can be modeled as a parameter preferences(p,t). We assume eight of each row has numbers 1,..,8 and the rest is zero (no preference). We assume further 8 is highest preference so we will maximize. To be able to run a model we need to simulate some data. We can use random data, but need to consider that talks t1 and t2 through t7 are a little bit more popular (they are repeated a few times). The algorithm is a little bit complicated, but here it is:



parameter italk(t,n) 'number all talks (all 21 of them)';
scalar counter/0/;
loop(talk(t,n),
counter=counter+1;
italk(talk)=counter;
);
display italk;

parameter preferences(p,t) 'preference = 0..8, 8 is highest preference';
preferences(p,t)=0;

scalar kk,kt;
alias(t,tt);
set tcurrent(t);
loop((p,tt)$(ord(tt)<=8),

repeat
* draw from 1..numtalk
kk = UniformInt(1,numtalk);
kt = sum(talk(t,n)$(italk(t,n)=kk),ord(t));
* preferences(p,'element with ord kt') = ord(t)
tcurrent(t)$(ord(t)=kt) = yes;
if (sum(tcurrent,preferences(p,tcurrent))=0,
preferences(p,tcurrent) = ord(tt);
else
kk = 0
);
tcurrent(tcurrent) = no;
until kk<>0;

);

display "after",preferences


We renumber talks in italks as follows:



----     84 PARAMETER italk  number all talks (all 21 of them)

n1 n2 n3

t1 1.000 2.000 3.000
t2 4.000 5.000
t3 6.000 7.000
t4 8.000 9.000
t5 10.000 11.000
t6 12.000 13.000
t7 14.000 15.000
t8 16.000
t9 17.000
t10 18.000
t11 19.000
t12 20.000
t13 21.000


This allows us to draw from unique talks 1..21, which gives a better distribution than if we would have drawn from 1..13 (e.g. we would generate too few attendants for the popular talk t1 which is given three times). Using this we can generate some data for preferences:



----    110 PARAMETER preferences  preference = 0..8, 8 is highest preference

t1 t2 t3 t4 t5 t6 t7 t8 t9

p1 6.000 1.000 4.000 5.000 7.000 3.000
p2 4.000 6.000 7.000 8.000 1.000 5.000 3.000
p3 2.000 3.000 8.000 1.000 4.000 6.000 7.000
p4 1.000 3.000 5.000 2.000 6.000 8.000 7.000
p5 3.000 5.000 4.000 2.000 1.000 7.000 6.000 8.000
p6 4.000 7.000 1.000 8.000 5.000 2.000 3.000 6.000
p7 3.000 7.000 2.000 1.000 8.000 6.000 5.000
p8 7.000 1.000 2.000 4.000 8.000 5.000
p9 1.000 6.000 7.000 3.000 4.000 2.000 5.000
p10 2.000 5.000 3.000 1.000 7.000 6.000 8.000
p11 7.000 6.000 3.000 5.000 2.000 4.000
p12 7.000 1.000 5.000 8.000 2.000 6.000 3.000
p13 8.000 1.000 4.000 6.000 2.000 3.000
p14 2.000 1.000 5.000 3.000 4.000 6.000
p15 2.000 3.000 6.000 5.000 8.000 1.000 7.000
p16 5.000 7.000 1.000 3.000 8.000 2.000
p17 4.000 1.000 8.000 6.000 5.000
p18 3.000 1.000 2.000 5.000 7.000 6.000 4.000
p19 6.000 8.000 7.000 5.000 1.000 4.000
p20 5.000 1.000 6.000 8.000 4.000
p21 1.000 2.000 5.000 3.000 7.000
p22 4.000 1.000 8.000 5.000 3.000 6.000 7.000
p23 2.000 5.000 6.000 8.000 3.000
p24 1.000 4.000 5.000 7.000 6.000 2.000 3.000 8.000
p25 2.000 4.000 6.000 5.000 1.000 7.000
p26 3.000 5.000 7.000 8.000 4.000 6.000
p27 2.000 6.000 8.000 1.000 4.000 3.000 7.000
p28 4.000 1.000 3.000 6.000 7.000 8.000
p29 2.000 3.000 4.000 8.000 1.000 5.000 6.000
p30 4.000 1.000 3.000 8.000 2.000 6.000 5.000
p31 5.000 6.000 2.000 7.000 1.000 4.000
p32 6.000 3.000 4.000 5.000 7.000 8.000
p33 2.000 4.000 6.000 5.000 1.000 8.000
p34 7.000 4.000 3.000 6.000 5.000
p35 1.000 7.000 8.000 2.000 6.000 4.000 3.000
p36 2.000 5.000 3.000 6.000 1.000 7.000
p37 1.000 3.000 5.000 4.000 6.000 2.000
p38 7.000 3.000 8.000 6.000 2.000 1.000
p39 2.000 1.000 4.000 8.000 7.000 5.000
p40 5.000 1.000 4.000 8.000 7.000 2.000 3.000
p41 4.000 7.000 2.000 3.000 5.000 6.000
p42 2.000 8.000 7.000 5.000 4.000 6.000
p43 5.000 2.000 7.000 6.000 1.000 3.000
p44 4.000 7.000 1.000 3.000 5.000
p45 3.000 7.000 1.000 4.000 5.000 8.000 2.000
p46 2.000 5.000 6.000 8.000 4.000 7.000
p47 7.000 5.000 3.000 8.000 2.000 1.000
p48 5.000 7.000 8.000 6.000 3.000
p49 2.000 5.000 1.000 3.000 8.000
p50 2.000 5.000 1.000 3.000 4.000 7.000

+ t10 t11 t12 t13

p1 2.000 8.000
p2 2.000
p3 5.000
p4 4.000
p7 4.000
p8 6.000 3.000
p9 8.000
p10 4.000
p11 1.000 8.000
p12 4.000
p13 5.000 7.000
p14 8.000 7.000
p15 4.000
p16 4.000 6.000
p17 7.000 3.000 2.000
p18 8.000
p19 3.000 2.000
p20 2.000 7.000 3.000
p21 4.000 6.000 8.000
p22 2.000
p23 7.000 1.000 4.000
p25 3.000 8.000
p26 2.000 1.000
p27 5.000
p28 5.000 2.000
p29 7.000
p30 7.000
p31 3.000 8.000
p32 2.000 1.000
p33 3.000 7.000
p34 8.000 2.000 1.000
p35 5.000
p36 8.000 4.000
p37 8.000 7.000
p38 4.000 5.000
p39 3.000 6.000
p40 6.000
p41 8.000 1.000
p42 1.000 3.000
p43 8.000 4.000
p44 6.000 8.000 2.000
p45 6.000
p46 1.000 3.000
p47 6.000 4.000
p48 2.000 1.000 4.000
p49 6.000 4.000 7.000
p50 6.000 8.000


Of course in practice we have his data available and we don't have to worry about this step (actually the most complicated part of the whole model). The rest of the model was implemented as follows:



variable z 'objective variable';
binary variables
x(t,p,s) 'assign talk/person/slot'
xts(t,s) 'assign talk/slot'
xtp(t,p) 'assign talk/person'
;

* we can relax xts and xtp
xts.prior(t,s) = INF;
xtp.prior(t,p) = INF;


equations
defxts1(t,s) 'calculate xts: all x=0 => xts = 0'
defxts2(t,p,s) 'calculate xts: any x=1 => xts = 1'

defxtp(t,p) 'calculate xtp'

talkcount(t) 'talk can be repeated ntalk times'
slotcount(s) 'up to 5 talks each period'

listen(p,s) 'p can visit only one talk in each time period'
zdef 'objective: maximize preferences'
;

* x(t,p,s)=0 for all p => xts(t,s) = 0
defxts1(t,s).. xts(t,s) =L= sum(p, x(t,p,s));
* any x(t,p,s)=1 => xts(t,s) = 1
defxts2(t,p,s).. xts(t,s) =G= x(t,p,s);

* p visits talk t in any s
defxtp(t,p).. xtp(t,p) =E= sum(s, x(t,p,s));

* each talk happens exactly ntalk times
talkcount(t).. sum(s, xts(t,s)) =E= ntalk(t);

* timeslot can only hold up to 5 talks
slotcount(s).. sum(t, xts(t,s)) =L= rooms;

* p can only visit one talk in each time period
listen(p,s).. sum(t,x(t,p,s)) =E= 1;

* we don't allow here (t,p) combinations without preference
xtp.fx(t,p)$(preferences(p,t)=0) = 0;

* objective: maximize
zdef.. z =e= sum((t,p), preferences(p,t) * xtp(t,p));


model scheduling /all/;


option mip=cplex;
scheduling.optfile=1;
scheduling.iterlim=100000;
scheduling.optcr=0;
solve scheduling using mip maximizing z;


$onecho > cplex.opt
threads 4
$offecho


Note that variables xts and xtp are automatically integer, so we relax them by assigning a priority of infinity to them.



Cplex solves this quite quickly:



--- Job conference2.gms Start 03/05/09 17:58:00 WEX-VIS 23.0.2 x86/MS Windows       
GAMS Rev 230 Copyright (C) 1987-2009 GAMS Development. All rights reserved
Licensee: Erwin Kalvelagen G090213/0001CV-WIN
Amsterdam Optimization Modeling Group DC4572
--- Starting compilation
--- conference2.gms(177) 3 Mb
--- Starting execution: elapsed 0:00:00.003
--- conference2.gms(169) 4 Mb
--- Generating MIP model scheduling
--- conference2.gms(170) 5 Mb
--- 4,234 rows 3,966 columns 17,496 non-zeroes
--- 3,250 discrete-columns
*** 465 relaxed-columns WARNING
--- Executing CPLEX: elapsed 0:00:00.034

ILOG CPLEX Feb 14, 2009 23.0.2 WIN 9101.9411 VIS x86/MS Windows
Cplex 11.2.1, GAMS Link 34
Cplex licensed for 1 use of lp, qp, mip and barrier, with 4 parallel threads.

Reading parameter(s) from "C:\projects\test\cplex.opt"
>> threads 4
Finished reading from "C:\projects\test\cplex.opt"
Reading data...
Starting Cplex...
Tried aggregator 1 time.
MIP Presolve eliminated 1501 rows and 1501 columns.
Reduced MIP has 2733 rows, 2465 columns, and 10595 nonzeros.
Reduced MIP has 2000 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec.
Clique table members: 650.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: opportunistic, using up to 4 threads.
Parallel mode: opportunistic, using up to 4 threads for concurrent optimization.
Tried aggregator 1 time.
No LP presolve or aggregator reductions.
Using devex.
Using devex.

Total real time on 4 threads = 0.17 sec.
Root relaxation solution time = 0.17 sec.

Nodes Cuts/
Node Left Objective IInf Best Integer Best Node ItCnt Gap

0 0 1500.0000 906 1500.0000 5
* 0+ 0 1272.0000 1500.0000 5 17.92%
0 0 1500.0000 779 1272.0000 Fract: 1 3652 17.92%
* 0+ 0 1308.0000 1500.0000 3652 14.68%
* 0+ 0 1353.0000 1500.0000 3652 10.86%
* 0+ 0 1446.0000 1500.0000 3652 3.73%
0 2 1500.0000 779 1446.0000 1500.0000 3652 3.73%
* 4+ 0 1460.0000 1500.0000 10828 2.74%
* 4+ 0 1471.0000 1500.0000 10683 1.97%
* 40+ 33 1472.0000 1498.0000 37861 1.77%
100 95 1492.8125 568 1472.0000 1498.0000 54702 1.77%
* 127+ 120 1478.0000 1498.0000 58492 1.35%
200 194 1487.6000 458 1478.0000 1498.0000 74738 1.35%
300 291 1492.6364 499 1478.0000 1498.0000 88945 1.35%
* 346+ 323 1483.0000 1498.0000 101352 1.01%
* 374+ 323 1488.0000 1498.0000 108821 0.67%
400 253 1493.0667 555 1488.0000 1498.0000 115116 0.67%
500 153 1495.0000 654 1488.0000 1498.0000 148604 0.67%
600 159 1491.4444 448 1488.0000 1497.3333 208207 0.63%
700 205 1492.8333 526 1488.0000 1496.6667 247331 0.58%
800 243 1494.6667 438 1488.0000 1496.6667 263583 0.58%
900 312 cutoff 1488.0000 1496.5000 282888 0.57%
1000 385 1494.9375 465 1488.0000 1496.0000 309440 0.54%
Elapsed real time = 32.45 sec. (tree size = 0.49 MB, solutions = 10)
1100 457 1492.1875 579 1488.0000 1495.7500 326655 0.52%
1200 524 cutoff 1488.0000 1495.5000 347873 0.50%
1300 588 1490.5000 386 1488.0000 1495.3333 371387 0.49%
1400 655 1489.4286 429 1488.0000 1495.0000 394081 0.47%
1500 706 cutoff 1488.0000 1495.0000 418373 0.47%
1600 760 cutoff 1488.0000 1494.9167 445114 0.46%
1700 810 1489.5000 355 1488.0000 1494.7500 466940 0.45%
1800 865 cutoff 1488.0000 1494.6667 487877 0.45%
1900 924 1490.3333 385 1488.0000 1494.5000 510279 0.44%
2000 972 1492.5000 523 1488.0000 1494.3333 536279 0.43%
Elapsed real time = 43.77 sec. (tree size = 1.27 MB, solutions = 10)
2100 1022 cutoff 1488.0000 1494.1429 561289 0.41%
* 2127 391 integral 0 1491.0000 1494.0714 567476 0.21%
2200 380 cutoff 1491.0000 1494.0000 590101 0.20%
2300 363 1492.0000 530 1491.0000 1493.6667 627287 0.18%
2400 358 cutoff 1491.0000 1493.5000 652433 0.17%
2500 326 cutoff 1491.0000 1493.2000 684303 0.15%
2600 272 1492.5625 605 1491.0000 1493.0000 722077 0.13%
2700 245 cutoff 1491.0000 1493.0000 749425 0.13%
2800 191 1492.0000 334 1491.0000 1492.5000 771850 0.10%
2900 107 cutoff 1491.0000 1492.2609 793601 0.08%
3000 15 cutoff 1491.0000 1492.0000 808231 0.07%
Elapsed real time = 57.28 sec. (tree size = 0.05 MB, solutions = 11)

Root node processing (before b&c):
Real time = 2.64
Parallel b&c, 4 threads:
Real time = 54.99
Sync time (average) = 2.47
Wait time (average) = 0.00
-------
Total (root+branch&cut) = 57.63 sec.
MIP status(101): integer optimal solution
Fixing integer variables, and solving final LP...
Parallel mode: opportunistic, using up to 4 threads for concurrent optimization.
Tried aggregator 1 time.
LP Presolve eliminated 4234 rows and 3966 columns.
All rows and columns eliminated.

Total real time on 4 threads = 0.00 sec.
Fixed MIP status(1): optimal

Proven optimal solution.

MIP Solution: 1491.000000 (814080 iterations, 3023 nodes)
Final Solve: 1491.000000 (0 iterations)

Best possible: 1491.000000
Absolute gap: 0.000000
Relative gap: 0.000000

--- Restarting execution
--- conference2.gms(170) 0 Mb
--- Reading solution for model scheduling
--- conference2.gms(170) 3 Mb
--- GDX File C:\projects\test\x.gdx
*** Status: Normal completion
--- Job conference2.gms Stop 03/05/09 17:58:58 elapsed 0:00:57.817

Note that we found a proven global optimum in less than a minute here. The complete schedule looks like:




schedule 
It is often the simplest to start with a single large multidimensional “cube” (our variable x(t,p,s)), and then derive auxiliary variables from that cube (our variables xts, xtp). If the problem becomes very large, this may not be possible (x just becomes too large). However, if this structure fits, in many cases the performance is very good. In this case we have 3250 binary variables. This is a large number. Is is quite a feat that we can solve this quickly to optimality.



See also [follow-up].

Wednesday, March 4, 2009

Excel COM problem: Call was rejected by callee

I was asked to help with a problem in GAMS/GDXXRW. It threw an exception with the error "Call was rejected by callee". After some experimenting, it turned out that the culprit was a pivot table (and graph). This pivot table was using data being updated by GDXXRW and after removing the pivot table the problem went away. A command line option is added to GDXXRW to give Excel some extra time after the COM Automation Server (i.e. Excel) is started. This will allow the pivot table to initialize without interference from GDXXRW.

Tuesday, March 3, 2009

big-M and log-normal

The terms big-M (in MIP modeling) and log-normal (the statistical distribution) are examples of concepts where the name has unintended, misleading connotations.

In MIP models a big-M constant should really be chosen as small as possible (for an example see here). I have seen many models where innocent students write something like M=10000000000 causing all kind of problems. Of course if this constant would have been called small-m, we probably would not have seen this as often.

The log-normal distribution is not formed by taking the logarithm of a normally distributed random variable (e.g. this posting). I suspect a better name like exp-normal distribution would have prevented this misconception.

OML: MS vs Ketron

OML is the name of the modeling language (Optimization Modeling Language) of Microsoft Solver Foundation. I just remember seeing that TLA (Three Letter Acronym, itself a TLA) before. Ketron has OML defined as Optimization & Modeling Library. Ketron Management Science was famous for its MPSIII/Whizard mainframe based LP solver (competitor of IBM's MPSX). I believe Whizard used super-sparsity (this technology has largely disappeared together with the punch-card).

Sunday, March 1, 2009

TSP in OML (MS Solver Foundation)

In theory it is easy to formulate a TSP when the all-different constraint is present. Let x[i] be the ith city visited, an OML-like formulation could look like:

Model[
  Parameters[Integers,N=14],
  Parameters[Sets[Integers],city],
  Parameters[Reals,dist[city,city]],
  Decisions[Integers[0,N-1],x[city]],
  Constraints[
     Unequal[Foreach[{i,city},x[i]]]
  ],
  Goals[
    Minimize[Sum[{i,city},dist[x[i],x[i++1]]]]
  ]
]

This does not work completely. First, bounds can not contain a symbolic constant. Yuck: we need to write [0,13] as bounds:  Decisions[Integers[0,13],x[city]]. The objective has  a few more difficulties. There is no circular lead/lag operator in OML (like the -- or ++ operator in GAMS).  This can be solved by introducing a parameter next[city]. This can be calculated in Excel and imported in the OML model. Furthermore, we cannot use a variable as an index, so dist[x[i],x[next[i]]] is not a valid construct (some languages geared towards solving CSP problems allow this; it adds concise expressiveness to the language that may help the modeler). The error message indicates this is not related to OML per se but rather to solver support. Also it is not allowed to use something like FilteredSum[{j,city},j==x[i],...] (a condition in a FilteredSum cannot depend on a decision variable). The only thing I could think of was:

Minimize[Sum[{i,city},{j,city},{k,city},AsInt[j==x[i]]*AsInt[k==x[next[i]]]*dist[j,k]]]

Note that I used integers as set elements. I started with using data-binding through tables, using set elements {'city0','city1',...,'city13'}. This made the model somewhat more complicated, as x[i] is an integer. So the CSP formulation uses the simpler set {0,1,...,13} allowing us to operate on indices. (Some consider this bad modeling: in GAMS arithmetic on sets is actually discouraged, and needs a function ord()to convert the element — always a string in GAMS — to an integer). Unfortunately I was not able to solve the small 14 city example with the CSP solver using this formulation. (Even after adding City[0]==0 as constraint). To check the input, I also tried a MIP formulation with Gurobi. That solved this small instance just fine.

In the MIP formulation we needed to formulate:

FilteredSum[{i,city},i != 'city0', x['city0',i]]==1

From what I can see this is not possible: a set element not being an integer can not be used in OML. As a workaround I created sets by binding to some dummy parameters:

Sum[{i0,city0},{i,city2},x[i0,i]]==1

where city0 is a set containing only a single element 'city0', and city2 is a set {'city1',...,'city13'}.  The constraint

FilteredSum[{i,city},{j,city},i != j,x[i,j]] <= N

did not work as i and j are not numeric. A workaround could be to have a parameter num[city] which can be calculated in Excel and then imported. Then the condition can read: num[i]!=num[j]. I just used dist[i,j]>0 as condition, as that also can be used to exclude the diagonal (the Excel spreadsheet makes sure all diagonal distances dist[i,i] are zero). A more general approach would be to be able to introduce sparse 2-d sets (like one could do in AMPL and GAMS), but OML has only one dimensional sets as far as I know.

This exercise was really meant to explore the expressiveness of OML. The little model actually shows some interesting issues with OML and some possible workarounds. I do not want to suggest this is a good way to solve TSP's. Obviously these approaches are not suited for large problems. A cutting plane algorithm often is quite effective for slightly larger problems (see this 42 city problem). The excel plugin does not allow us to to implement such an algorithm: an OML model can only contain a single model and has no looping facilities. For the real large problems, you need to look elsewhere, such as Concorde.