This is an example of a Multiple traveling salesman problem or uncapacitated VRP. A very simple formulation though.
| $ontext |
Tuesday, March 31, 2009
MTSP/VRP example
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
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:
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[Reals,MeanAdjRet[I,T]],
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],
Foreach[{k,K},{t,T},w[k,t]==FilteredSum[{i,I},Selected[i]>0,MeanAdjRet[i,t]*Alloc[k,i]]]
],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
| 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:
Here the parameter Selected indicates whether an instrument can be part of the portfolio.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]]
]
]]
]
Thursday, March 26, 2009
Drawing TSP Tour in GAMS
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
task}
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 |
| LOWER LEVEL UPPER MARGINAL ---- VAR z 5.7000 5.7000 5.7000 1.0000 |
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:
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.