## Tuesday, March 31, 2009

### MTSP/VRP example

> You have an example of a VRP (Vehicle Routing Problem)?

This is an example of a Multiple traveling salesman problem or uncapacitated VRP. A very simple formulation though.

### AMPL, OML version of GAMS Scheduling Model

In the good tradition of letting blogs link to each other, here is a blog posting with an AMPL and OML version of this GAMS model. Those are nicely formatted listings.

## Sunday, March 29, 2009

### Efficient Portfolio Frontier in MSF (2)

Here is an interesting but obscure reformulation that sometimes performs better than the traditional x’Qx term in a mean-variance portfolio model. We notice that the covariances can be calculated as:

where μi is the mean return and di,t are the historical data. We now can write:

where

This is especially a good formulation if there are more instruments than time periods (in large practical problems this is often the case as older historical data is not considered as having much predictive value). The corresponding OML model, where we solve K sub-problems at the same time, can look like (see post for more info):

Model[
Parameters[Sets,I,K,T],
Parameters[Integers,Selected[I]],
Parameters[Reals,Return[I]],
Parameters[Reals,Lambda[K]],
Parameters[Integers,TT[]],

Decisions[Reals[0,Infinity],Alloc[K,I]],
Decisions[Reals,w[K,T]],

Constraints[
Foreach[{k,K},FilteredSum[{i,I},Selected[i]>0,Alloc[k,i]]==1],
],

Goals[Minimize["Overall"->
Sum[{k,K},
Sum[{t,T},w[k,t]^2]/TT[]
- Lambda[k]*FilteredSum[{i,I},Selected[i]>0,Return[i]*Alloc[k,i]]
]
]]
]

As an aside: this formulation also shows how this can be linearized: replace ∑w2 by ∑|w|.

## Saturday, March 28, 2009

### Yet another big-M problem, don’t blame Cplex

Hi,

I posted this message at the ILOG CPLEX Forum but I guess it could be
of interest here. I know that some people watch both forums, so I
apologize for the cross-post.

Any idea is appreciated. Could it be a problem on our cplex
installation ? We are using Ubuntu (one of the 2008 versions, I can
check exactly which one if anyone finds that information important).
---

Dear all,

I am really puzzled with the following:

I have an IP model. ( …) . If I solve this model (I am using Cplex 11.2, interative
optimizer, under Ubuntu ), I obtain an "optimal" integer solution with
cost 23 (…
).

HOWEVER, I know of a solution with cost 21 ( … ). Indeed, if I
pass this solution to CPLEX as an mst file, Cplex starts using it as
the "Best Integer found" ( …
).

Am I missing anything ? Any help is appreciated,

It was suggested by another poster to use some clever but perfectly sensible options:

I verified this with CPLEX 11.0 (also on Ubuntu).  There are two ways (that I found) to get CPLEX to produce a correct solution (objective = 21) without the hot start:  turn off the presolver; or crank up the simplex feasibility tolerance (1e-9 worked for me).  Turning off the presolver seems to be the faster option.  (Disclaimer:  In both instances, I aborted the run after it found an integer-feasible solution with value 21 but before it exhausted the search tree, so it's possible one or the other option might have led to problems down the road.)

Cplex support in this forum adds more suggestions for even more exotic options (again they make perfect sense). However, I would approach this a little bit differently. First this problem is not purely a Cplex problem. Other solvers show the same problem. The reason lies actually in the formulation, which causes problems for relaxation-based branch & bound methods. It is always difficult to recognize the model structure from an LP file (please, use a proper modeling language) but there are big-M constants in the model with value 100,000. After I replaced all these instances by 10,000 Cplex shows the expected behavior: it does not stop with obj=23 but goes further to obj=21 on Ubuntu. So rather than tinkering with tolerances I would work on the model formulation and provide big-M values that are as small as possible. These values are often related to some upper bound. Experienced modelers (you can sometimes recognize them by their grey haires and thick glasses) will try to prevent using just a big arbitrary value like 1e5 but rather work on deducing a tight bound. Once you have proper formulation you will see that the model will behave better for any MIP solver. I would argue that a solver independent fix is better than using a solver specific option that works around the problem. (If you absolutely cannot provide a good value for your big-M’s you may want to look into Cplex’s indicator constraints).

### Efficient Portfolio Frontier in MSF

I received a few questions about how to calculate an efficient portfolio frontier in MS Solver Foundation. The simplest algorithm is to solve a series of QP’s. (A better way may be a parametric QP analysis). As we don’t have looping constructs in OML, it is not immediately obvious how to do this with the Excel plug-in. With the API’s we could do the looping in C# and the Quadratic Programming Problem itself in OML. Here is an alternative approach: solve a bigger problem that has all the individual QP’s embedded. Conceptually this would look like:

 Obj: x1’Qx1-λ1 μ'x1 +x2’Qx2-λ2 μ'x2 +x3’Qx3-λ3 μ'x3 +x4’Qx4-λ4 μ'x4 S.t. ∑ix1,i=1 ∑ix2,i=1 ∑ix3,i=1 ∑ix4,i=1

Note that xk is a vector rather than a scalar, and that I ignored the non-negativity conditions xk≥0. A similar approach can be used for any problem where the individual problems are relatively small. A good example is DEA where we achieve better performance by batching small problems into bigger ones.

The OML formulation of this combined portfolio model can look like:

Model[
Parameters[Sets,I,K],
Parameters[Integers,Selected[I]],
Parameters[Reals,Return[I]],
Parameters[Reals,Lambda[K]],
Parameters[Reals,Covar[I,I]],

Decisions[Reals[0,Infinity],Alloc[K,I]],
Constraints[
Foreach[{k,K},FilteredSum[{i,I},Selected[i]>0,Alloc[k,i]]==1]
],

Goals[Minimize["Overall"->
Sum[{k,K},
FilteredSum[{i,I},Selected[i]>0,
FilteredSum[{j,I},Selected[j]>0,Alloc[k,i]*Covar[i,j]*Alloc[k,j]]]
- Lambda[k]*FilteredSum[{i,I},Selected[i]>0,Return[i]*Alloc[k,i]]
]
]]
]

Here the parameter Selected indicates whether an instrument can be part of the portfolio.

## Thursday, March 26, 2009

### Drawing TSP Tour in GAMS

Can I draw a TSP tour in GAMS?

Yes, use 2D vectors. Here are a few examples: a 29 city problem and a 51 city problem solved as a MIP.

## Wednesday, March 25, 2009

### A machine scheduling problem

Hello

I'm new to AMPL.  I have  a scheduling problem that is giving me
trouble. I have a set of unrelated machines.  I have a set of jobs.
Each job has a set of operations to perform in order.  I want to
minimize the time of completion of the last job/operation pair
processed. Below is my .mod file. In the data section I have two test
cases - one is commented and one is not.

The commented test has no problems. The solver finds the solution. The
uncommented test case gives the following:

CPLEX Error  1217: No solution exists.
CPLEX Error  1217: No solution exists.
CPLEX Error  1262: No basis exists.
CPLEX 7.1.0: integer infeasible.
1894 MIP simplex iterations
39 branch-and-bound nodes; no basis.

Like I said I am new to this type of programming. Does anybody smarter
than me see a problem with any of my constraints? Or a comment on what
usually causes error 1217 and 1262?

#################################### .mod file
set job;
set oper;
set mach;
param first > 0;
param last  > first;
set T := first .. last;

var max_val;
var asgn { job, oper, mach, T } binary;
var sched  { mach, T } >= 0;
var delay  { mach, T } >= 0;

param time { job, oper, mach } >= 0;

minimize z: max_val;

# the last scheduled task for each machine is less than or equal to
our objective
subject to max_val_constraint { m in mach }:  sched[m,last] <=
max_val;

# sched[m,t] is the completion time of the task machine m is doing at
time t
# sched[m,t] = sched[m,t] + delay[m,t] + { time to process assigned
subject to mach_sched_0 { m in mach, t in T: t>1}:
sched[m,t] = sched[m,t-1] + delay[m,t] +
sum{ i in job, j in oper } asgn[i,j,m,t] * time[i,j,m];

subject to mach_sched_1 { m in mach }:
sched[m,1] = delay[m,1] +
sum{ i in job, j in oper } asgn[i,j,m,1] * time[i,j,m];

# for each job/operation pair, only assign it to be processed once
subject to jobtask_assign_cnt { i in job, j in oper }:
sum{ m in mach, t in T } asgn[i,j,m,t] = 1;

# for each machine at time t, assign only one task to process
subject to machtime_assign_cnt { m in mach, t in T }:
sum{ i in job, j in oper } asgn[i,j,m,t] <= 1;

# for job i oper j, make sure job i, oper j-1 has been previously assigned
subject to jobtask_ordering { i in job, j in oper, t in T: j>1 }:
sum{ m in mach } asgn[i,j,m,t] >= sum{ m in mach } asgn[i,j-1,m,t];

subject to delay_constraint { m in mach, t in T }:
delay[m,t] >= 0;

data;

set job  := 1 2 3 ;
set oper := 1 2 3 4;
set mach := 1 2 ;

param first := 1;
param last  := 12;

param time :=
[1,*,*]: 1 2 :=
1 1 4
2 2 2
3 3 3
4 4 4

[2,*,*]: 1 2 :=
1 2 4
2 2 1
3 1 6
4 1 2

[3,*,*]: 1 2  :=
1 3 4
2 2 3
3 3 1
4 4 3;

# test that works
#set job  := 1 2;
#set oper := 1 2;
#set mach := 1 2 3;
#set T := 1 2 3 4;
#
#param time :=
#[1,*,*]: 1 2 3 :=
#1 1 2 3
#2 1 2 3
#
#[2,*,*]: 1 2 3 :=
#1 1 2 3
#2 1 2 3;

The constraint

# for job i oper j, make sure job i, oper j-1 has been previously assigned
subject to jobtask_ordering { i in job, j in oper, t in T: j>1 }:
sum{ m in mach } asgn[i,j,m,t] >= sum{ m in mach } asgn[i,j-1,m,t];

does not seem to be formulated correctly. It says that all operations have to be executed in the same time period t. For this data set you have four operations and only two machines available. The smaller data set had two operations which can be allocated to the available two machines. If you really want (i,j-1) to be finished before (i,j) you probably need to keep track when each task (i,j) finishes. I believe your current formulation mixes up time slots t and completion times modeled via sched[m,t]. If your execution times are small integer, then map everything to time slots. If your execution times are floating point numbers, you should look into a continuous time formulation (e.g. through x[i,j,ii,jj,m] = 1 iff (i,j) precedes (ii,jj) on machine m).

### MS Solver Foundation job

It was brought to my attention that already there are (short term) Microsoft Foundation jobs available: http://www.careerbuilder.com/JobSeeker/Jobs/JobDetails.aspx?job_did=J3G7ML643DJNGKK4315. That means people start to do serious work with this!

### rminlp’s

Hi all,
I am working in pharma and here I need you help for a problem related to a pharmacokinetic equation.
I would like to have the (X,Y,S, T, U, V, W) possible values that will lead to a value of Z=5.7.
Z= 1/X*1/(Y/(1+(U×V)/(W×(U+S)))+(1-Y)) )
knwing that:
0.57<X<0.79
0.00016<W<0.00032
0.89<Y<0.94
0.0936<U<0.967
V =0.05 ou 0.072
S = 29 ou 5.3

This GAMS model calculates some solutions:

variables z,x,y,u,v,w,s;
equation e;

e.. z =e= (1/X)*1/(Y/(1+(U*V)/(W*(U+S)))+(1-Y));

x.lo = 0.57;    x.up = 0.79;
w.lo = 0.00016; w.up = 0.00032;
y.lo = 0.89;    y.up = 0.94;
u.lo = 0.0936;  u.up = 0.0967;

binary variable b1,b2;
equation e1,e2;
e1.. v =e= 0.05*b1 + 0.072*(1-b1);
e2.. s =e= 29*b2 + 5.3*(1-b2);

z.fx = 5.7;

model m /all/;

b1.fx=0; b2.fx=0; solve m minimizing z using rminlp;
b1.fx=0; b2.fx=1; solve m minimizing z using rminlp;
b1.fx=1; b2.fx=0; solve m minimizing z using rminlp;
b1.fx=1; b2.fx=1; solve m minimizing z using rminlp;

I could only find feasible solutions for b2=0 (i.e. s=5.3). This was confirmed by the global solver Baron. Here they are:

 LOWER          LEVEL          UPPER         MARGINAL ---- VAR z                  5.7000         5.7000         5.7000         1.0000      ---- VAR x                  0.5700         0.6019         0.7900          .          ---- VAR y                  0.8900         0.8900         0.9400         EPS         ---- VAR u                  0.0936         0.0936         0.0967         EPS         ---- VAR v                 -INF            0.0720        +INF             .          ---- VAR w                  0.0002         0.0003         0.0003         EPS         ---- VAR s                 -INF            5.3000        +INF             .          ---- VAR b1                  .              .              .             EPS         ---- VAR b2                  .              .              .             EPS LOWER          LEVEL          UPPER         MARGINAL ---- VAR z                  5.7000         5.7000         5.7000         1.0000      ---- VAR x                  0.5700         0.5700         0.7900         EPS         ---- VAR y                  0.8900         0.9400         0.9400         EPS         ---- VAR u                  0.0936         0.0965         0.0967          .          ---- VAR v                 -INF            0.0500        +INF             .          ---- VAR w                  0.0002         0.0003         0.0003         EPS         ---- VAR s                 -INF            5.3000        +INF             .          ---- VAR b1                 1.0000         1.0000         1.0000         EPS         ---- VAR b2                  .              .              .             EPS

## Tuesday, March 24, 2009

### Endogenous if-then-else

> I am trying to define a nonlinear function where the variable in the
> objective is dependent on a parameter (x) which varies across
> i=1,2,...,N cases and two variables, a cut-off value (cutoff) and a
> coefficient (coeff) which transforms x above the value of cutoff. The
> logic, in pseudocode is:
>
> if x[i] > cutoff {
>     x2[i] = x[i] * coeff
>     }
> else {
>     x2[i] = x[i]
>     }
>
> The objective function would be something like:
>
> minimize Objective: sum {i in N} [complicated function goes here] * x2[i] ;
>
> How, if at all, can this be implemented in AMPL?
>

You could try:
e{i in I}: x2[i] = if (x[i]>cutoff) then x[i]*coeff else x[i];
However this is a discontinuous function so beware that many NLP solvers assume smooth functions and may get into trouble near the jump. A smooth approximation would be better in that case (there are several possible functional forms for such a smooth approximation; one would need to invest some time to try a few of them; an example is shown here). If the model was otherwise linear, one could look into a piecewise linear formulation (AMPL has very good support for this; this approach is less interesting if the model is otherwise nonlinear as it would become an MINLP).

## Monday, March 23, 2009

### Indexing by a variable

> I've:
>
> param N_MAX = 1000;
> var N >= 1, <= N_MAX;
>
> var x {0 .. N_MAX} >= 0;
>
> I want to do *something like*:
>
> subject to {i in  {(N/2)-1, (N/2)+1} }: #calculation of indices
> simplified for sake of discussion
>     x[i] = 1000; #RHS simplified for sake of discussion
>
> But of course I can't use variables for indexing. Neither can I do it
> explicitly as there are too many possible values that N can take. So,
> I'm looking for suggestions.
>
> Sincerely,

param N_MAX := 10;

set I := {1..N_MAX};
set J := {0..N_MAX};
set IJ{i in I} := setof{j in J: ((j >= i/2-1) and (j <= i/2+1))}(j);
display IJ;

var N >= 1, <= N_MAX;
param xup := 10000;
var x {J} >= 0, <= xup;
var b {I} binary;

minimize obj: sum{j in J} x[j];
s.t. e1: sum{i in I} i*b[i] = N;
s.t. e2: sum{i in I} b[i] = 1;
s.t. e3{i in I, j in IJ[i]}:  x[j] <= 1000*b[i] + xup*(1-b[i]);
s.t. e4{i in I, j in IJ[i]}:  x[j] >= 1000*b[i];

solve;
display N,b,x;
end;

Usually it is best to introduce a binary variable b[i] with b[i]=1 <=> i=N. Then you can introduce constraints based on b[i]. xup is here an upper bound on x[j]. Choose this number as tight as possible: it is a big-M constant. As usual it is a good idea to move as much as possible complexity away from the constraints. In this case we use an intermediate set IJ to handle a more difficult concept. Such a set can be tested and debugged before a complete model is ready to run.

### MSF Examples

There were some requests to make my MS Solver Foundation examples available. I have placed them here. The examples include:

1. TSP in OML (post)
2. Tiling squares in OML (post1 and post2)
3. Sudoku (post)
4. Maxflow (post)
5. Langford’s problem (post)
6. Magic Squares (post)

## Thursday, March 19, 2009

### MIP model

Hi,
I am dealing with an optimization task as follows:
1) Tom is given a large box which contains 1000 bags of marbles.
2) Inside each bag, there are between 1 and 50 marbles.
3) Within a given bag, there are no duplicate colors of marbles, however there are duplicate colors of marbles across the bags.
4) The bags are labeled and Tom has a list of what color marbles each bag contains.
Objective:  Tom wants to collect as many UNIQUE colors of marbles as possible, but is only allowed to pick out 100 bags from the box.  Tom derives absolutely no value from getting duplicate colors of marbles.  Tell Tom which 100 bags to pick.
My question:
First of all, is this even a problem that's appropriate for linear programming? I was able to solve this quite simply on a small scale in Excel's Solver, but when I moved over to GLPK I couldn't figure out how to avoid making my objective statement non-linear.  Initially I wrote out the problem like below.  However, glpk complained about having a variable in my if/else statement in the objective function.  I'm still rather new to LP so if someone could tell me whether or not I should even be using LP for this, and if so, I'd appreciate some suggestions on how to model the problem such that GLPK doesn't complain. Thanks.
/* sets */
set COLORS;
set BAGS;
/* parameters */
param colors_to_bags {i in COLORS, j in BAGS}
/* decision variables */
var selected_bags {j in BAGS} binary >= 0;
/* objective function */
maximize uniques: sum{i in COLORS} if (sum{j in BAGS} colors_to_bags[i, j] * selected_bags[j]) > 1 then 1 else 0;
/* constraints */
s.t. max_bags: sum{j in BAGS} selected_bags{j} <= 100;
data;
set BAGS := bag1 bag2 .... bag1000;
set COLORS := green yellow blue ...etc; (assume there are many many colors)
/* 1 means the bag contains that color marble, 0 means it does not */
praram colors_to_bags: bag1 bag2 .... bag1000:=
green 0 1 0 0 ...
yellow 0 1 0 1 ...
blue 1 0 0 1 ...
etc. 0 1 1 0 ...;
end;

You need to introduce an extra variable have_color[i]:

var have_color{i in COLORS} binary;

/* all colors_to_bags[i,j]*selected_bags[j] = 0 =>have_color[i]=0 */
eq1{i in COLORS}: have_color[i] <= sum{j in BAGS} colors_to_bags[i,j]*selected_bags[j];

/* any colors_to_bags[i,j]*selected_bags[j] = 1 =>have_color[i]=1 */
eq2{i in COLORS, j in BAGS}: have_color[i] >= colors_to_bags[i,j]*selected_bags[j];

Then the obj becomes:

maximize uniques: sum{i in COLORS} have_color[i];

You can improve performance by relaxing have_color to be a continuous variable between 0 and 1 (it is automatically integer). In addition you may drop eq2 as you are maximizing this anyway. We are lucky: eq2 is adding many equations.

### MSF Solver Foundation API Question

Hi,
I am trying to port an Excel sheet that uses Solver the add-in from Frontline Systems to a C# program that uses the Solver Fondation API.
The sheet uses solver to find the combination of stocks in a portfolio for a given variance.
I have followed the Markowitz sample that comes with the express installation but one of the constraints that have to be added is that the portfolio's variance should be a certain value. The variance of the portfolio is calculated the following way:
SQRT(MMULT(MMULT(D38:L38;\$B\$10:\$J\$18);TRANSPOSE(D38:L38)))*SQRT(252)
where D38:L38 is the allocation vector (the result) and \$B\$10:\$J\$18 is the covariance matrix of the available stocks.
thank,

I believe this ends up as a nonlinear constraint. In OML-like notation (I invent the Sqrt function here):

`  Sqrt[Sum[{i,I},{j,I},x[i]*covar[i,j]*x[j]]] * Sqrt[252] == SomeValue`

This shows an obvious advantage of using MSF/OML compared to Excel formulas: it is much easier to write readable constraints. Excel's Solver and its bigger brother from Frontline Systems can handle such nonlinear constraints through its GRG2 based NLP solver, but MSF does not support this. MSF supports convex QP's, but I don't think this can be reformulated into a linearly constrained QP.

Notes:

1. the lack of Sqrt is not a show-stopping problem, as we could write:

Sum[{i,I},{j,I},x[i]*covar[i,j]*x[j]] == SomeOtherValue

2. For large problems this construct can be improved somewhat, first by using symmetry of Covar, and depending on the modeling system/solver by reformulating the quadratic form. As in your case the number of instruments is small, there is no need for that.

3. For very large problems sometimes a linear approximation is used (easy if minimizing variance, actually not so easy in this case).

## Tuesday, March 17, 2009

### MSF CSP: Tiling Squares (2)

In this post I displayed a CSP model for Microsoft Solver Foundation that solves the Tiling Squares problem. It was observed that a symmetry breaking constraint was really helpful. Basically it forced similar squares to be number in the x-axis direction. Of course a small extension would be to say: if same squares have the same x-coordinate, number them in the y direction. The symmetry breaking constraints can look like:

// this is to reduce some symmetry
Foreach[{i,T},{j,i+1,T},Implies[b[i]&b[j]&(Side[i]==Side[j]),x[j]>=x[i]]],
Foreach[{i,T},{j,i+1,T},Implies[
b[i]&b[j]&(Side[i]==Side[j])&(x[j]==x[i]),
y[j]>=y[i]]]

This reduces the calculation time further:

A bigger improvement can be achieved by observing that we can force lower numbered squares to be added first.

E.g. squares 9,10 and 11 have size 4. If we turn on square 10 we want to require that square 9 is also used. This can be implemented with the constraint:

FilteredForeach[{i,1,T},Side[i]==Side[i-1],Implies[b[i],b[i-1]]]

This has much effect on the computation time:

If we add numbers to the squares we see how the numbering follows these rules:

The reference sequence A036444 shows a(19)=2. I.e. we can use just a set with a maximum of two identical tiles to cover an area of (19 × 19). Indeed we can verify this with this model. The data and results are:

## Monday, March 16, 2009

### MSF CSP: Tiling Squares

The tiling problem is not very amenable for solving by MIP solvers (see tiling.pdf). In this problem we try to cover an (n × n) square with smaller squares. In a MIP we need big-M formulations to express the non-overlap condition. With CSP we can write this down more directly and conveniently. Instead of using something like:

we can write using Microsoft Solver Foundation:

Foreach[{i,T},{j,i+1,T},
x[i]>=x[j]+Side[j] |
x[i]+Side[i]<=x[j] |
y[i]>=y[j]+Side[j] |
y[i]+Side[i]<=y[j]
]

Here the vertical bar means OR. A small (7 × 7) problem can be implemented quickly as follows:

Here we already used knowledge about which tiles to use. We only needed to determine where to place them. The real problem is slightly more complicated: we don’t know which tiles to use in advance. Below, we are looking for a configuration for the (9 x 9) problem where no tile can be used more than three times. This is number h(9) in Integer Square Tilings and a(9) in sequence A036444. To formulate this problem we add boolean variables b(i) indicating if tile i is being used.

Also we added an extra constraint to reduce symmetry (this caused significant speed-up). The solution times with and without the symmetry breaking constraint are:

 Without symmetry breaking constraint With symmetry breaking constraint

I did not make an exhaustive examination of all solver options, so there may be better performance possible. The complete model looks like:

``Model[    Parameters[Integers,N=9],   // size of square to fill    Parameters[Integers,T=24],   // number of tiles      Parameters[Sets,Tiles],       Parameters[Integers,Side[Tiles]],           Decisions[Booleans,b[Tiles]],    Decisions[Integers[0,8],x[Tiles],y[Tiles]],       Constraints[           Foreach[{i,Tiles},Implies[b[i],x[i]<=N-Side[i] & y[i]<=N-Side[i]]],           // Note: Use T instead of Tiles         Foreach[{i,T},{j,i+1,T},Implies[b[i]&b[j],             x[i]>=x[j]+Side[j] |             x[i]+Side[i]<=x[j] |             y[i]>=y[j]+Side[j] |             y[i]+Side[i]<=y[j]         ]],          Sum[{i,Tiles},Side[i]^2 * AsInt[b[i]]] == N^2,                  // this is to reduce some symmetry        Foreach[{i,T},{j,i+1,T},Implies[b[i]&b[j]&(Side[i]==Side[j]),x[j]>=x[i]]]    ]  ]``

We use the Implies[] construct here to implement if-then. Note that we used T instead of Tiles when we needed to operate on the upper-triangle of the (i,j) pairs. In this case we could not use Tiles as OML then assumes the indices are non-numeric (well, they are actually integers).

Update: more symmetry breaking constraints and their effect on solution times are discussed in this follow-up.

## Sunday, March 15, 2009

### Fixed LP after MIP

Hi everybody,

After I solve a problem as an IP, I fix some of the variables
(complicating ones) to integer values of optimal solution to the IP
and solve the problem as a LP. The objective function only depends on
the fixed variables in both problems and there are constraints which
will be binding also those that will not. The problem is highly
degenerated.

The problem now is that I can not get any non-zero values for dual
variables using getDual or getDuals.

Anyone likes to put a comment on?

Cheers,

GAMS does this with all its MIP solvers to produce marginals for the rows and columns.

At the end of your fixed LP you should have a complete basis with m basics and n−m non-basics. The (non-binding) rows are likely to be basic, hence your zero's. The fixed integer variables are likely to be non-basic at bound, so you should see some nonzero marginals there.

I am sure your solver can produce a solution listing or a basis file so you can have a look where the non-basics end up and their marginal value (the dual value or reduced cost).

## Saturday, March 14, 2009

### AMPL/Cplex: No solver needed

Hello,

Anyone can help?
I have run the following model using AMPL code with CPLEX 11.2 solver
in window xp. The program run well and take a couple seconds for the
small size of problem. But once the problem size is more than 200 (of
j and k), the solve time increases dramatically. When the problem size
of set j and set k is 900 ,  it took so long ( longer than 10 hrs) to
find the solution. I wonder, since the model looks simple with only 2
variables per constraint. I'm not sure did I write the AMPL code wrong
solution from CPLEX after it done. The solution is correct). Thank you
very much.

Best regards,
Alex

The Mathematical Model is:

Objective: Min Delta
S.t. (1) Xj *(Bj + Sum(i) M(ICij + TC1ij + TC2ijk)) <= Delta,     For all j, k

(2) Sum(i) Xj = 1

(3) Xj = {0, 1}

The AMPL code is:

set stag ;
set sup;
set w;

param tc1{i in sup,j in stag}>=0;
param b{j in stag}>=0;
param ic {i in sup,j in stag}>=0;
param M integer;
param tc2{i in sup,j in stag,k in w}>=0;

var x{j in stag} binary;
var Delta;

minimize cost: Delta;
subject to delt{k in w, j in stag}:  x[j]*(b[j]+ sum{i in sup} (M*(ic[i,j]+ tc1[i,j]+ tc2[i,j,k]) )) <= Delta;
subject to c: sum{j in stag} x[j] = 1;

Ten hours is a long time for 900 binary variables in this rather simple model. The problem has many rows as a result of the way the max is formulated. This construct where we minimize the max is sometimes difficult for MIP solvers. Note that some MIP models are just very difficult to solve (example). However in this case, I believe this problem can be solved extremely fast without even using Cplex. First we calculate a[j,k] such that we can simplify the problem to:

minimize cost: Delta;
subject to delt{k in w, j in stag}:  x[j]*a[j,k] <= Delta;
subject to c: sum{j in stag} x[j] = 1;

The first constraint can be interpreted as: if (x[j]=1) then Delta = maxk a[j,k]. The complete model can now be restated as:

param a{j in stag, k in w} := b[j] + sum{i in sup} (M*(ic[i,j]+ tc1[i,j]+ tc2[i,j,k]));
param d{j in stag} := max{k in w} a[j,k];
param dmin := min{j in stag} d[j];
var x{j in stag} binary := if (d[j]=dmin) then 1 else 0;
var Delta := dmin;
display x,Delta;

(I left as an exercise to make sure there is only one x[j]=1 in case multiple d[j]’s assume the value dmin). If you really want to call Cplex (especially after shelling out big bucks for it), the following should be very fast:

var x{j in stag} binary;
var Delta;

minimize cost: Delta;
subject to c: sum{j in stag} x[j] = 1;
subject to c2: sum{j in stag} x[j]*d[j] = dmin;
subject to c3: Delta = dmin;

option solver cplex;
solve;

This formulation automatically deals with the problem of multiple d[j]’s being equal to the minimum.

## Friday, March 13, 2009

### MSF QP problem

There is some issue with MS Solver Foundation’s QP. First I thought it was related to a reported problem about not reporting the constant term in objective values. But this must be something else:

gives:

The MSF people reported they know about the problem and it will be fixed for the next release. A workaround is to simplify the objective to obj->x*y-x-y+1. In that case the objective value is not reporting the constant term, but the levels of the variables x and y are correct.

## 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:

The OML model can now look like:

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

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:

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:

Using this additional constraint, the solution looks like:

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!

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.000t2        4.000       5.000t3        6.000       7.000t4        8.000       9.000t5       10.000      11.000t6       12.000      13.000t7       14.000      15.000t8       16.000t9       17.000t10      18.000t11      19.000t12      20.000t13      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.000p2        4.000       6.000       7.000                   8.000       1.000       5.000                   3.000p3        2.000       3.000       8.000       1.000                   4.000       6.000                   7.000p4        1.000       3.000       5.000                   2.000       6.000       8.000       7.000p5        3.000       5.000       4.000       2.000       1.000       7.000       6.000                   8.000p6        4.000       7.000       1.000       8.000                   5.000       2.000       3.000       6.000p7        3.000       7.000       2.000       1.000       8.000                               6.000       5.000p8        7.000       1.000       2.000       4.000                                           8.000       5.000p9        1.000       6.000       7.000       3.000       4.000       2.000       5.000p10       2.000       5.000       3.000       1.000       7.000       6.000                               8.000p11       7.000       6.000                   3.000       5.000       2.000       4.000p12       7.000       1.000       5.000                   8.000                   2.000       6.000       3.000p13       8.000       1.000                   4.000                   6.000       2.000                   3.000p14       2.000       1.000                               5.000       3.000                   4.000       6.000p15       2.000       3.000       6.000       5.000       8.000       1.000       7.000p16       5.000       7.000                   1.000                   3.000       8.000                   2.000p17       4.000                               1.000       8.000       6.000       5.000p18       3.000                   1.000                   2.000       5.000       7.000       6.000       4.000p19       6.000       8.000                               7.000       5.000       1.000                   4.000p20       5.000       1.000       6.000                   8.000                               4.000p21       1.000                   2.000       5.000       3.000                   7.000p22       4.000       1.000       8.000       5.000       3.000                               6.000       7.000p23       2.000                               5.000       6.000       8.000       3.000p24       1.000       4.000       5.000       7.000       6.000       2.000       3.000       8.000p25       2.000       4.000                   6.000                   5.000       1.000       7.000p26       3.000       5.000       7.000                   8.000       4.000                               6.000p27       2.000       6.000       8.000       1.000       4.000       3.000       7.000p28       4.000       1.000       3.000       6.000       7.000                               8.000p29       2.000       3.000       4.000                   8.000       1.000       5.000       6.000p30       4.000                   1.000       3.000       8.000       2.000                   6.000       5.000p31       5.000                   6.000                   2.000                   7.000       1.000       4.000p32                   6.000       3.000                   4.000       5.000       7.000                   8.000p33       2.000                   4.000       6.000       5.000       1.000                               8.000p34       7.000       4.000       3.000                               6.000                   5.000p35       1.000       7.000       8.000       2.000       6.000       4.000       3.000p36       2.000       5.000                   3.000       6.000       1.000       7.000p37       1.000                   3.000       5.000       4.000       6.000                   2.000p38                               7.000       3.000       8.000       6.000       2.000                   1.000p39       2.000                               1.000       4.000       8.000       7.000       5.000p40       5.000       1.000       4.000       8.000       7.000       2.000                   3.000p41                   4.000                   7.000       2.000       3.000       5.000                   6.000p42       2.000       8.000       7.000       5.000                   4.000       6.000p43       5.000       2.000       7.000       6.000       1.000       3.000p44                                                       4.000       7.000       1.000       3.000       5.000p45       3.000       7.000                   1.000       4.000       5.000       8.000                   2.000p46       2.000       5.000       6.000                   8.000       4.000                   7.000p47       7.000       5.000       3.000       8.000       2.000                   1.000p48       5.000       7.000       8.000       6.000       3.000p49       2.000       5.000       1.000       3.000                               8.000p50       2.000       5.000       1.000       3.000       4.000                                           7.000  +         t10         t11         t12         t13 p1        2.000                               8.000p2                                            2.000p3        5.000p4                    4.000p7                                4.000p8                                6.000       3.000p9                                8.000p10                               4.000p11                               1.000       8.000p12                               4.000p13       5.000                               7.000p14       8.000                               7.000p15                               4.000p16                               4.000       6.000p17       7.000       3.000       2.000p18                                           8.000p19       3.000                               2.000p20       2.000       7.000                   3.000p21       4.000       6.000       8.000p22                               2.000p23       7.000       1.000       4.000p25       3.000                   8.000p26       2.000                               1.000p27                               5.000p28       5.000                   2.000p29       7.000p30                               7.000p31                   3.000                   8.000p32                   2.000                   1.000p33       3.000                               7.000p34                   8.000       2.000       1.000p35                               5.000p36       8.000                               4.000p37       8.000                   7.000p38                   4.000                   5.000p39       3.000                   6.000p40                   6.000p41       8.000                               1.000p42       1.000                   3.000p43       8.000       4.000p44       6.000                   8.000       2.000p45                   6.000p46                   1.000                   3.000p47                   6.000       4.000p48       2.000       1.000       4.000p49       6.000                   4.000       7.000p50                   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 xtpxts.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) = 0defxts1(t,s).. xts(t,s) =L= sum(p, x(t,p,s));* any x(t,p,s)=1 => xts(t,s) = 1defxts2(t,p,s).. xts(t,s) =G= x(t,p,s); * p visits talk t in any sdefxtp(t,p).. xtp(t,p) =E= sum(s, x(t,p,s)); * each talk happens exactly ntalk timestalkcount(t).. sum(s, xts(t,s)) =E= ntalk(t); * timeslot can only hold up to 5 talksslotcount(s).. sum(t, xts(t,s)) =L= rooms; * p can only visit one talk in each time periodlisten(p,s)..  sum(t,x(t,p,s)) =E= 1; * we don't allow here (t,p) combinations without preferencextp.fx(t,p)\$(preferences(p,t)=0) = 0; * objective: maximizezdef.. 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.optthreads 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 reservedLicensee: 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 WindowsCplex 11.2.1, GAMS Link 34Cplex licensed for 1 use of lp, qp, mip and barrier, with 4 parallel threads. Reading parameter(s) from "C:\projects\test\cplex.opt">>  threads 4Finished 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.64Parallel 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 solutionFixing 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.000000Absolute gap:            0.000000Relative 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:

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.

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