GAMS has strict domain checking (type checking). This has some interesting consequences. Basically all good: much better protection against errors compared to simply using integer indices \(1,\dots,n\). Consider the square matrix \(\color{darkblue}A_{i,j}\):
|
set |
> How to linearize sum{i in I} binvar[i]*intvar[i] ≤ C where binvar is a binary variable and intvar is an integer variable.
In general we can linearize y[i] = binvar[i]*intvar[i] as:
0 ≤ y[i] ≤ intvar.up[i]
y[i] ≤ binvar[i]*intvar.up[i]
y[i] ≤ intvar[i]
y[i] ≥ intvar[i] - intvar.up[i]*(1-binvar[i])
In this case we can drop the ≤ inequalities, so we have:
0 ≤ y[i] ≤ intvar.up[i]
y[i] ≥ intvar[i] - intvar.up[i]*(1-binvar[i])
sum{i in I} y[i] ≤ C
> How to formulate c[i]=min(a[i],b[i]) for use with Cplex. A and b are integer.
There is not much we can exploit knowing that a, b are integer. In this formulation I use an extra binary variable bv[i] that implements the OR-functionality. The AMPL/Cplex formulation can look like:
var bv{I} binary;
e1{i in I}: c[i] <= a[i];
e2{i in I}: c[i] <= b[i];
e3{i in I}: bv[i]=0 ==> c[i]=a[i] else c[i]=b[i];
e3 implements c[i] = a[i] or c[i] = b[i] and e1,e2 make sure the minimum of the two is chosen.
> How to formulate c=min(a,b) where a,b are binary in a MIP.
This is like c = a×b. See: http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-binary-variables.html.
I.e.
c ≤ aSimilarly c = max(a,b) can be modeled via the constraints:
c ≤ b
c ≥ a + b − 1
c in [0,1] (continuous)
a,b in {0,1} (binary)
c ≥ a
c ≥ b
c ≤ a + b
c in [0,1] (continuous)
a,b in {0,1} (binary)
> How to linearize z = x × b with x integer and b binary?
Use the formulation from http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-continuous-and-binary.html.
> You had a note on how to model oligopolistic producer behavior in GAMS?
Uploaded to http://www.amsterdamoptimization.com/pdf/oligopoly.pdf.
The paper shows that an MCP formulation is natural for this. It also compares prices and quantities with the situation where the firms are merged and a large monopolist is formed. This behavior is easily modeled using a simple NLP model. A numerical experiment like this can be quickly implemented using a modeling language.
Here is an example how to do probit estimation with GAMS. We use the dataset Table F21.1 from Greene. We estimate by forming a likelihood function which we can maximize using a standard NLP solver. First we do OLS to get a good starting point for the NLP. The OLS estimation is done through the specialized LS solver.
The reason to use a specialized least squares solver is twofold: (1) it is not easy to do least squares reliably using standard (linear) optimization tools. The normal equations allow for using an LP but in forming the normal equations you may create numerical difficulties. Of course one could solve a QP. (2) Some of the statistics useful in reporting OLS results are not trivial to compute in GAMS. A specialized solver can solve LS problems quickly and reliably with linear technology and in additional produce all kind of statistics useful in assessing the quality of the fit.
The max likelihood estimation problem is to maximize:
Here y is the dependent (binary) variable, and x are the independent variables (note: the variables are data in the optimization problem). β is the vector of coefficients to estimate (these are the decision variables in the optimization problem). Φ is the CDF of the standard normal distribution. This expression is reformulated into:
$ontext
Probit Estimation
We use OLS to get a good starting point
Erwin Kalvelagen, Amsterdam Optimization, 2009
Data:
http://pages.stern.nyu.edu/~wgreene/Text/tables/TableF21-1.txt
$offtext
set i /1*32/;
table data(i,*)
GPA TUCE PSI GRADE
1 2.66 20 0 0
2 2.89 22 0 0
3 3.28 24 0 0
4 2.92 12 0 0
5 4.00 21 0 1
6 2.86 17 0 0
7 2.76 17 0 0
8 2.87 21 0 0
9 3.03 25 0 0
10 3.92 29 0 1
11 2.63 20 0 0
12 3.32 23 0 0
13 3.57 23 0 0
14 3.26 25 0 1
15 3.53 26 0 0
16 2.74 19 0 0
17 2.75 25 0 0
18 2.83 19 0 0
19 3.12 23 1 0
20 3.16 25 1 1
21 2.06 22 1 0
22 3.62 28 1 1
23 2.89 14 1 0
24 3.51 26 1 0
25 3.54 24 1 1
26 2.83 27 1 1
27 3.39 17 1 1
28 2.67 24 1 0
29 3.65 21 1 1
30 4.00 23 1 1
31 3.10 21 1 0
32 2.39 19 1 1
;
set k 'independent variables' /constant,gpa,tuce,psi/;
parameters
y(i) 'grade'
x(k,i) 'independent variables'
;
y(i) = data(i,'grade');
x('constant',i) = 1;
x(k,i)$(not sameas(k,'constant')) = data(i,k);
parameter estimate(k,*);
*-----------------------------------------------------------
* O L S
*-----------------------------------------------------------
variable sse,coeff(k);
equation obj,fit(i);
obj.. sse =n= 0;
fit(i).. y(i) =e= sum(k, coeff(k)*x(k,i));
model ols /obj,fit/;
option lp=ls;
solve ols using lp minimizing sse;
estimate(k,'OLS') = coeff.l(k);
*-----------------------------------------------------------
* P R O B I T
*-----------------------------------------------------------
variable logl;
equation like;
like.. logl =e= sum(i$(y(i)=1), log(errorf(sum(k,coeff(k)*x(k,i)))))
+sum(i$(y(i)=0), log(1-errorf(sum(k,coeff(k)*x(k,i)))));
model mle /like/;
solve mle using nlp maximizing logl;
estimate(k,'Probit') = coeff.l(k);
display estimate;
The results are identical to the numbers published in Greene.
---- 99 PARAMETER estimate
OLS Probit
GPA 0.464 1.626
TUCE 0.010 0.052
PSI 0.379 1.426
constant -1.498 -7.452
Updated the LS solver docs to include this example.
Hi,
For a project I'm working on, I need to solve the following problem. Suppose
we have a squared-grid (size=n) with n^2 cells with a 4-connexity
neighborhood. Each cell must contains exactly one building. We have several
building colors (blue, red and green) with an increasing amount of points.
The game rules are the following:
* No constraints on blue buildings.
* Red buildings must have at least one blue building in its neighborhood.
* Green buildings must have at least one red building and at least one blue
building in its neighborhood.
The goal is to maximize the number of points by having the biggest buildings
(green > red > blue). I have started to write the LP but I have some
difficulties to express the constraint on red buildings because it directly
depends on the values of the variables. So, I would like to transform the
following constraint
param n, integer, > 0, default 4;
var x{1..n, 1..n}, integer, >=0, <=2; /*blue=0, red=1 and green=2 */
s.t. r{i in 1..n, j in 1..n:x[i,j]=1}: sum{a in i-1..i+1, b in j-1..j+1:a>=1
and a<=n and b>=1 and b<=n and i!=a and b!=j and x[a,b]=0} x[i,j] >= 1;
into a valid one. I think there's a way to express it by using binary
variables but I don't see how. Does anybody can help me ? Thanks in advance
for your help.
Indeed you can not use a variable to drive sets (basically the sets are handled by the modeling system and the variables by the solver, so there is a phase difference; another way to look at it is that the solver expects a system of purely linear inequalities). It is better to use a binary variable with an extra index for the color. That makes implementing the logical condition much easier. The model is simple once you decided on how the variables are organized. Here is the model:
set C; # colors
param n;
set I := 1..n;
param points{C};set neighbor{i in I, j in I} := setof{i1 in max(i-1,1)..min(i+1,n),j1 in max(j-1,1)..min(j+1,n)}(i1,j1);
var x{I,I,C} binary;
maximize obj:
sum{i in I, j in I, c in C} points[c]*x[i,j,c];
OneColor{i in I, j in I}:
sum{c in C} x[i,j,c] = 1;RedRequirement{i in I, j in I}:
sum{(i1,j1) in neighbor[i,j]} x[i1,j1,'blue'] >= x[i,j,'red'];
GreenRequirement1{i in I, j in I}:
sum{(i1,j1) in neighbor[i,j]} x[i1,j1,'red'] >= x[i,j,'green'];
GreenRequirement2{i in I, j in I}:
sum{(i1,j1) in neighbor[i,j]} x[i1,j1,'blue'] >= x[i,j,'green'];
solve;for{i in I}
{
for{j in I}
{
printf if x[i,j,'red']>0.5 then 'R '
else if x[i,j,'green']>0.5 then 'G '
else if x[i,j,'blue']>0.5 then 'B ';
}
printf "\n";
}data;
set C := red green blue;
param n := 5;
param points := green 10, red 5, blue 2;
end;
This gives:
Reading model section from colors.mod...
Reading data section from colors.mod...
48 lines were read
Generating obj...
Generating OneColor...
Generating RedRequirement...
Generating GreenRequirement1...
Generating GreenRequirement2...
Model has been successfully generated
ipp_basic_tech: 1 row(s) and 0 column(s) removed
ipp_reduce_bnds: 1 pass(es) made, 0 bound(s) reduced
ipp_basic_tech: 0 row(s) and 0 column(s) removed
ipp_reduce_coef: 1 pass(es) made, 0 coefficient(s) reduced
glp_intopt: presolved MIP has 100 rows, 75 columns, 657 non-zeros
glp_intopt: 75 integer columns, all of which are binary
Scaling...
A: min|aij| = 1.000e+000 max|aij| = 1.000e+000 ratio = 1.000e+000
Problem data seem to be well scaled
Crashing...
Size of triangular part = 100
Solving LP relaxation...
0: obj = 2.500000000e+002 infeas = 5.000e+001 (0)
* 54: obj = 2.060000000e+002 infeas = 0.000e+000 (0)
* 78: obj = 2.110000000e+002 infeas = 2.029e-015 (0)
OPTIMAL SOLUTION FOUND
Integer optimization begins...
Gomory's cuts enabled
MIR cuts enabled
Cover cuts enabled
Clique cuts enabled
Creating the conflict graph...
The conflict graph has 2*75 vertices and 150 edges
+ 78: mip = not found yet <= +inf (1; 0)
+ 352: >>>>> 1.930000000e+002 <= 2.070000000e+002 7.3% (13; 0)
+ 909: >>>>> 1.980000000e+002 <= 2.060000000e+002 4.0% (30; 8)
+ 6635: mip = 1.980000000e+002 <= tree is empty 0.0% (0; 211)
INTEGER OPTIMAL SOLUTION FOUND
Time used: 1.0 secs
Memory used: 0.8 Mb (855246 bytes)
G G G G G
B R G B R
G G G G G
G B G G G
R G G R B
Model has been successfully processed
Someone else suggested to use:
# binaries indicating color;
var is_blue{i in N, j in N}, binary;
var is_green{i in N, j in N}, binary;
var is_red{i in N, j in N}, binary;
# color counts
var blues;
var greens;
var reds;
I don't think this results in a very clean model (invalid link:thread has been removed). If you have similar variables that behave in the same way, in general it is better to add an index position so the variables can be folded into one larger structure. We can even take this approach one step further: instead of three similar constraints we can use a single constraint block using an additional set describing the implications. Note that the suggested model here (invalid link:thread has been removed) gives a different solution than my model (even after changing the objective to make the points identical) because the border is handled differently from what the poster proposed.
Update: the set neighbor should be rewritten as:
set neighbor{i in I, j in I} := setof{i1 in max(i-1,1)..min(i+1,n),j1 in max(j-1,1)..min(j+1,n): i1!=i and j1!=j}(i1,j1);
display neighbor;
to reflect what the poster indicated (see the comments).
The other post that caught my eye was this one:
Hi Everyone,
I am looking for some help with the following problem:
Let's say I have a number of print jobs (tax statements, bills etc.)
Each print job has 1 basestock and up to 6 brochures.
e.g. A tax statement prints on blue basestock and has 2 brochures (My Tax and Tax Breaks for Pensioners)
The machines that the print jobs run on have 6 hoppers for brochures and 2 trays for basestock.
I can combine the print jobs into sets, thus minimising the set up time for the machine.
e.g. If I combine 3 jobs into 1 batch, then when I run this batch on the machine I only need to to one setup.
The problem seems to be, to find the minimum number of collections that the sets can be grouped into.
e.g. I have the following print jobs:
Job 1:
Basestock A, Brochures 1, 2, 3, 4Job 2:
Basestock B, Brochures 2,3Job 3:
Basestock A, Brochures 5,6,7I can combine Jobs 2 and 3, becuase this gives me a batch with 2 basestocks and 6 brochures which is within the capabilities of the machine.
So I am left with two batches:
Batch 1: Job 1Batch 2: Jobs 2 and 3
I am looking for an algorithm that will provide the smallest number of
collectionsAny assistance would be much appreciated
I don’t think the description or my understanding of it is completely correct. E.g. I think job 2 + job 3 leads to a batch of 5 brochures instead of 6. If we combine jobs 1 and 2 we use brochures 2,3 in for both jobs. I assume this is ok, and that it actually saves a spot for brochures. The advice was given to use a MIP model (this is actually not a bad idea at all), but the suggested model is not that good:
Perhaps you can make a ILP model to solve it. Suppose you have n jobs,you can define binary variable x[i,j] (1 <= i <= n,1 <= j <= n) to show if the job i and the job j are in the same batch, As jobs are given,you can get a parameter c[i,j] to show if the job i and the job j can be in the same batch (for your example, c[2,3]=1 ,c[1,2]=0).Then
you set up such constraints:x[i,j]=1 (i!=j) means that job i and job j are in the same batch
x[i,j]=1 (i==j)means that job i are put in the batch solelyfor any i: sum( x[i,j])=1 1 <= j <= n
for any i,j :x[i,j]=x[j,i]
for any i: if c[i,j]==0 then x[i,j]=0 1 <= j <= nLast,you set up the object:
min = sum( x[i,j]) 1 <= i <= n, i <= j <= n
This does not look right. First you need probably something like a binary variable x[i,k] indicating whether job i goes into batch k. Then the rest of the model should follow. Assuming I interpret the problem correctly, the model could look like
set BROCHURES;
set BASESTOCK;
set JOBS;
set BATCHES;set Brochures{JOBS};
set BaseStock{JOBS};param numHoppers;
param numTrays;var x{JOBS,BATCHES} binary;
var batchUse{BATCHES} >= 0, <= 1;
var batchBrochure{BATCHES,BROCHURES} >= 0, <= 1;
var batchBaseStock{BATCHES,BASESTOCK} >= 0, <= 1;minimize obj: sum{k in BATCHES} batchUse[k];
batchUsage{i in JOBS, k in BATCHES}: x[i,k] <= batchUse[k];
allJobs{i in JOBS}: sum{k in BATCHES} x[i,k] = 1;
calcBatchBrochures{k in BATCHES,b in BROCHURES,i in JOBS:b in Brochures[i]}:
batchBrochure[k,b] >= x[i,k];
brochureCapacity{k in BATCHES}: sum{b in BROCHURES} batchBrochure[k,b] <= numHoppers;
calcBatchBaseStock{k in BATCHES,b in BASESTOCK,i in JOBS:b in BaseStock[i]}:
batchBaseStock[k,b] >= x[i,k];
baseStockCapacity{k in BATCHES}: sum{b in BASESTOCK} batchBaseStock[k,b] <= numTrays;solve;
display x;data;
set BROCHURES := 1 2 3 4 5 6 7;
set BASESTOCK := A B;
set JOBS := job1 job2 job3;
set BATCHES := batch1 batch2 batch3;
set Brochures['job1'] := 1 2 3 4;
set Brochures['job2'] := 2 3;
set Brochures['job3'] := 5 6 7;
set BaseStock['job1'] := A;
set BaseStock['job2'] := B;
set BaseStock['job3'] := A;
param numHoppers := 6;
param numTrays := 2;
end;
I am a little bit careful here in exploiting that some variable are automatically integer when they become binding. So we only need to keep x as binary. These relaxed variables (batchUse, batchBrochure, batchBaseStock) can be considered as bounds, and they are free to move when not binding. In the case of batchUse the objective will drive the variable down to zero. In the case of batchBrochure, batchBaseStock we don’t really care about their precise value, as long as they fulfill their role in making sure the capacity is not exceeded. Putting it differently, the calcBatchBrochures and calcBatchBaseStock inequalities implement the logical condition:
if x[i,k]=1 then batchBrochure[k,b]=1 else leave batchBrochure[k,b] unrestricted
if x[i,k]=1 then batchBaseStock[k,b]=1 else leave batchBaseStock[k,b] unrestricted
When batchBrochure[k,b], batchBaseStock[k,b] are left floating, the capacity constraints brochureCapacity and baseStockCapacity may or may not force them to zero. So at the end of the solve, some of the values of batchBrochure[k,b], batchBaseStock[k,b] may be not equal to zero while there is no usage. The only thing you can rely on, is that if there is usage of brochures or base stock these value will be one (the reverse will not hold in general).
It is not known how large the BATCH set should be in advance. In this case we could have used just 2 elements. An upper bound on the number of elements needed can be established by allowing as many batches as there are jobs. The result from this model will look like:
Reading model section from print.mod...
Reading data section from print.mod...
44 lines were read
Generating obj...
Generating batchUsage...
Generating allJobs...
Generating calcBatchBrochures...
Generating brochureCapacity...
Generating calcBatchBaseStock...
Generating baseStockCapacity...
Model has been successfully generated
ipp_basic_tech: 4 row(s) and 0 column(s) removed
ipp_reduce_bnds: 1 pass(es) made, 0 bound(s) reduced
ipp_basic_tech: 0 row(s) and 0 column(s) removed
ipp_reduce_coef: 1 pass(es) made, 0 coefficient(s) reduced
glp_intopt: presolved MIP has 51 rows, 39 columns, 120 non-zeros
glp_intopt: 9 integer columns, all of which are binary
Scaling...
A: min|aij| = 1.000e+000 max|aij| = 1.000e+000 ratio = 1.000e+000
Problem data seem to be well scaled
Crashing...
Size of triangular part = 51
Solving LP relaxation...
0: obj = 0.000000000e+000 infeas = 1.400e+001 (0)
* 25: obj = 1.000000000e+000 infeas = 0.000e+000 (0)
* 30: obj = 1.000000000e+000 infeas = 0.000e+000 (0)
OPTIMAL SOLUTION FOUND
Integer optimization begins...
Gomory's cuts enabled
MIR cuts enabled
Cover cuts enabled
Clique cuts enabled
Creating the conflict graph...
The conflict graph has 2*9 vertices and 18 edges
+ 30: mip = not found yet >= -inf (1; 0)
+ 105: >>>>> 2.000000000e+000 >= 1.333333333e+000 33.3% (4; 0)
+ 120: mip = 2.000000000e+000 >= tree is empty 0.0% (0; 7)
INTEGER OPTIMAL SOLUTION FOUND
Time used: 0.0 secs
Memory used: 0.2 Mb (239814 bytes)
Display statement at line 29
x[job1,batch1] = 1
x[job1,batch2] = 0
x[job1,batch3] = 0
x[job2,batch1] = 0
x[job2,batch2] = 1
x[job2,batch3] = 0
x[job3,batch1] = 0
x[job3,batch2] = 1
x[job3,batch3] = 0
Model has been successfully processed
These post illustrate that good modeling is difficult. It requires paying attention to details as well as keeping an eye on the model as a whole. In that sense it is often more difficult than programming where often once break down difficult tasks into more manageable pieces.
Hi Dear all,
I am new to GAMS and I am working on a MILP. I had been my code in AMPL.
Because of my solution capacity for the n > 8 had was too large, and I have a student version of AMPL then I want to change it in the GAMS representation.
I have two loop using FOR in my code, as:
for{i in 1..N} {
for{j in 1..N} {
for{p in 1..N} {
let a[i,j,p,i,j,p] := f[j,j] * d [ i,i] * f[p,p] * d[i,i] ;
for{k in 1..N} {
if(k <> i) then
for{n in 1..N} {
if(n <> j) then
for{q in 1..N} {
if(q <> p) then
let a[i,j,p,k,n,q]:= f[j,n] * d[i,k] * f[p,q] * d[i,k];
}
}
}
}
}
};In the loop, f and d are two-dimensional matrix and i,j,p,k,n,q mention to elements of the matrices and N is a scalar as follow.
I have written the same loop in GAMS as:
* ----------------------------------------------
sets i /1*5/,
j /1*5/,
p /1*5/,
k /1*5/,
n /1*5/,
q /1*5/ ;scalar M size of problem /5/;
parameter c(i,j,p,k,n,q);
Scalar N /5/ ;
* ----------------------------------------------
Table f(i,i)
1 2 3 4 5
1 1 1 2 3
2 1 2 1 2
3 1 2 1 2
4 2 1 1 1
5 3 2 2 1 ;* ------------------------------------------------------
Table d(k,k)
1 2 3 4 5
1 5 2 4 1
2 5 3 2
3 2 3
4 4 5
5 1 2 5 ;* ------------------------------------------------------
for(i=1 to M,
for(j=1 to M,
for(p=1 to M,
a(i,j,p,i,j,p) = f(j,j) * d(i,i) * f(p,p) * d(i,i) ;
for(k=1 to M,
if(k ne i,
for(n=1 to M,
if(n ne j,
for(q=1 to M,
if(q ne p,
a(i,j,p,k,n,q) = f(j,n) * d(i,k) * f(p,q) * d(i,k) ;
);
);
);
);
);
);
);
);
);* ----------------------------------------------
But there is an error regard to FOR. I know that it should be a scalar instead of i,j,p,k,n,q. but then who can I define the matrices?
Thanks in advance
This does not make much sense. Wow, this is very ugly and convoluted GAMS. At least from a syntactical point of view better would be:
* ----------------------------------------------
sets i /1*5/;
alias(i,j,p,k,n,q);* ----------------------------------------------
parameter a(i,j,p,k,n,q);
Table f(i,i)
1 2 3 4 5
1 1 1 2 3
2 1 2 1 2
3 1 2 1 2
4 2 1 1 1
5 3 2 2 1 ;
* ------------------------------------------------------
Table d(k,k)
1 2 3 4 5
1 5 2 4 1
2 5 3 2
3 2 3
4 4 5
5 1 2 5 ;
* ------------------------------------------------------
a(i,j,p,k,n,q) = f(j,n) * f(p,q) * sqr(d(i,k));
When you use GAMS in most cases you don’t need explicit loops: an assignment is an implicit loop. Also there is no advantage to split this into two assignments as suggested by the poster. Essentially he does:
a(i,j,p,i,j,p) = f(j,j) * f(p,p) * sqr(d(i,i));
a(i,j,p,k,n,q)$(ord(k)<>ord(i) and ord(n)<>ord(j) and ord(q)<>ord(p)) = f(j,n) * f(p,q) * sqr(d(i,k));
These two statements have the same effect as the single statement shown before.
There seems to be a big misunderstanding here how to write GAMS and actually how to write loops (the splitting of the assignment is an indication of that).
> How can I solve a Minimum Spanning Tree problem with GAMS?
It is possible to implement an algorithm in GAMS (Prim, Kruskal, see model library model mst.gms). It is also possible to formulate a MIP. There are several formulations (see Magnanti, Wolsey, “Optimal Trees”, in Network Models, volume 7, Handbooks of Operations Research, 1995). If the graph is planar, there is a formulation described in Justin C. Williams, “A linear-size zero - one programming model for the minimum spanning tree problem in planar graphs”, Networks, 39(1), 2002.
Here is one of the flow formulations mst.gms.
| $ontext |
As the MIP is automatically integer (total unimodularity) we can solve as an RMIP. This is a big LP however: for a 42 city problem with all edges i → j allowed, I get:
MODEL STATISTICS
BLOCKS OF EQUATIONS 3 SINGLE EQUATIONS 72,325
BLOCKS OF VARIABLES 3 SINGLE VARIABLES 72,325
NON ZERO ELEMENTS 284,131 DISCRETE VARIABLES 72,324
PS: I believe the model mst.gms from the model library is incorrect. The assignment mst(a,u,ip)$(dmin eq darc(a,u)) = ord(s); should be changed so that if multiple arcs (a,u) have length dmin, only one is chosen. A corrected version for the above instance is mst2.gms. The fix was needed for this data set: without it a suboptimal tree was returned: obj=662 instead of obj=591 (as demonstrated by msterr.gms). This is another good reason to implement this as an LP: to debug algorithms.