Gurobi Standalone version available

The Gurobi standalone solver is now available. More information is available on the website www.gurobi.com. There is reference material online to peruse. I have some experience with using Gurobi, and although I have not done any formal benchmarks, on my models it has shown to be very competitive with the other leading LP/MIP solvers and remarkably stable. Gurobi seems especially attractive when you have a machine with many cores/cpu’s available (and soon we all have cell-phones with 16 cores).

Interestingly the shell is based on Python. I am probably one of the very few ones who remember ABC, a small language developed at the CWI that was in many respects the origin of Python.

CSP: indexing by a variable

A table lookup construct y=value[x] where x is an integer variable can be implemented in a MIP model as follows:

Parameters[Sets[Integers],I],
Parameters[Reals,value[I]],
Decisions[Integers[0,1],b[I]],
Constraints[
    Sum[{i,I},b[i]]==1,
    x == Sum[{i,I},i*b[i]],
    y == Sum[{i,I},value[i]*b[i]],
    ....

Here x and y are variables (they can be relaxed to be continuous as they are automatically integer – that is depending on value[] for y). Some experienced MIP users may recognize a SOS1 structure here. In a CSP model this can sometimes be expressed more succinctly as:

y == value[x]

This is also (and more coherently) discussed in the AMPL extensions document: http://www.ampl.com/NEW/FUTURE/logic.html#varsub. More info here: http://users.iems.northwestern.edu/~4er/WRITINGS/extmodcp.pdf. These references contain some motivating examples.

The Microsoft Solver Foundation group will add this to their OML language: http://code.msdn.microsoft.com/solverfoundation/Thread/View.aspx?ThreadId=1739. This will allow for very natural and intuitive formulations for some models.

One-liner optimization problems

As long as your problem fits on one line, Wolfram’s Alpha engine may be able to solve it:  http://www.wolframalpha.com/input/?i=optimization+examples. Not sure if we can solve any larger problems with this.

Probit estimation with GAMS

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

This model looks much simpler than what is discussed here.

Updated the LS solver docs to include this example.

MS OML model using C# and data from Access

Looking at http://code.msdn.microsoft.com/solverfoundation/Thread/List.aspx it is noticeable that many users have problems coding simple problems in C#. One possible reason is that they use C# to assemble the model. Although not extremely difficult, it is gives very unwieldy models: the signal-to-noise ratio in the code is small as you need lots of syntactic clutter just to specify all variables and constraints compared to a specialized modeling language. Large models can easily have dozens of blocks of variables and equations. In this post I want to emphasize an alternative that is somewhat underrated: it is possible to use OML directly in your C# application. That will immediately make the model much more compact and readable.

The second issue is that the data-binding is often not completely straightforward. Many questions are related to data binding. Below is the simplest solution I could come up with for the following architecture: the math programming model is a simple transportation model and all data is stored in an Access database. The goal is to provide a skeleton example that is both readable and simple while being useful as a starting point for larger, more complex applications. The advanced features of LINQ as used throughout the Solver Foundation documentation are largely geared towards SQL Server. Therefore I wanted to explore how a simpler database like Access could be handled. If Access is working, there is good reason to believe that any other major database will also work, as we are working with the lowest common denominator in some respects. Many databases are accessible through OleDb. In practice it may be a problem that all data has to come from the database: OML has no facilities for data manipulation. Large models often have large amount of data, which may need some form of processing (aggregation etc.). Even if your real database is say Oracle, it may be useful to use Access as front-end tool to perform these data manipulation steps.

The model is the trnsport.gms model from the GAMS model library. It is small and has a few small parameters. For more info see http://www.amsterdamoptimization.com/models/msf/oml.pdf. In OML the model looks like:

Model[
  Parameters[Sets,Plants,Markets],
  Parameters[Reals,Capacity[Plants],Demand[Markets],Cost[Plants,Markets]],

  Decisions[Reals[0,Infinity],x[Plants,Markets],TotalCost],

  Constraints[
     TotalCost == Sum[{i,Plants},{j,Markets},Cost[i,j]*x[i,j]],
     Foreach[{i,Plants}, Sum[{j,Markets},x[i,j]]<=Capacity[i]],
     Foreach[{j,Markets}, Sum[{i,Plants},x[i,j]]>=Demand[j]]
  ],

  Goals[Minimize[TotalCost]]
]

In C# we can do:

/// <summary>
/// Holds the OML model
/// </summary>
string strModel = @"Model[
      Parameters[Sets,Plants,Markets],
      Parameters[Reals,Capacity[Plants],Demand[Markets],Cost[Plants,Markets]],

      Decisions[Reals[0,Infinity],x[Plants,Markets],TotalCost],

      Constraints[
         TotalCost == Sum[{i,Plants},{j,Markets},Cost[i,j]*x[i,j]],
         Foreach[{i,Plants}, Sum[{j,Markets},x[i,j]]<=Capacity[i]],
         Foreach[{j,Markets}, Sum[{i,Plants},x[i,j]]>=Demand[j]]
      ],

      Goals[Minimize[TotalCost]]
   ]";

followed by:

SolverContext context;
context.LoadModel(FileFormat.OML, new StringReader(strModel));
Solution solution = context.Solve();
Console.Write("{0}", solution.GetReport());

This was easy and as short as can be. Now we need to get the data. The database is organized as:

We will use the tables Capacity and Demand and the Query Cost.  The data looks like:

To bind the data we use the following code:

/// <summary>
/// Solve the problem
/// </summary>
public void Solve()
{
    context.LoadModel(FileFormat.OML, new StringReader(strModel));

    foreach (Parameter p in context.CurrentModel.Parameters)
    {
        switch (p.Name)
        {
            case "Capacity":
                setBinding(p,"select plant,capacity from capacity",
                    "capacity", new string[]{"plant"});
                break;
            case "Demand":
                setBinding(p,"select market,demand from demand",
                    "demand", new string[]{"market"});
                break;
            case "Cost":
                setBinding(p,"select plant,market,cost from cost",
                    "cost", new string[]{"plant", "market"});
                break;
        }

    }

    Solution solution = context.Solve();
    Console.Write("{0}", solution.GetReport());

}

In each binding operation we specify:

1. The SFS parameter, which we retrieve from the CurrentModel
2. The query to be used against the database
3. The name of the data column
4. The names of the index columns (passed on as an array of strings)

The complete model looks like:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Data;
using System.Data.OleDb;
using System.Data.Linq;
using System.Text;
using Microsoft.SolverFoundation.Services;
using System.IO;

namespace OML1
{
    class Trnsport
    {
        /// <summary>
        /// Called by the OS
        /// </summary>
        /// <param name="args"></param>
        static void Main(string[] args)
        {
            Trnsport t = new Trnsport();
            t.Solve();
        }

        /// <summary>
        /// Holds the OML model
        /// </summary>
        string strModel = @"Model[
              Parameters[Sets,Plants,Markets],
              Parameters[Reals,Capacity[Plants],Demand[Markets],Cost[Plants,Markets]],

              Decisions[Reals[0,Infinity],x[Plants,Markets],TotalCost],

              Constraints[
                 TotalCost == Sum[{i,Plants},{j,Markets},Cost[i,j]*x[i,j]],
                 Foreach[{i,Plants}, Sum[{j,Markets},x[i,j]]<=Capacity[i]],
                 Foreach[{j,Markets}, Sum[{i,Plants},x[i,j]]>=Demand[j]]
              ],

              Goals[Minimize[TotalCost]]
           ]";

        /// <summary>
        /// Connection string for MS Access 
        /// Use x86 architecture!
        /// </summary>
        string connection = @"Provider=Microsoft.ACE.OLEDB.12.0;Data Source=C:\projects\ms\OML1\OML1\trnsport.accdb;Persist Security Info=False;";

        /// <summary>
        ///  SFS
        /// </summary>
        SolverContext context;

        /// <summary>
        ///  Constructor
        /// </summary>
        public Trnsport()
        {
            context = SolverContext.GetContext();
        }

        /// <summary>
        /// get query result as DataSet
        /// </summary>
        /// <param name="connection">connection string</param>
        /// <param name="query">query as string</param>
        /// <returns></returns>
        private DataSet SelectOleDbSrvRows(string connection, string query)
        {
            DataSet ds = new DataSet();
            OleDbConnection conn = new OleDbConnection(connection);
            OleDbDataAdapter adapter = new OleDbDataAdapter();
            adapter.SelectCommand = new OleDbCommand(query, conn);
            adapter.Fill(ds);
            return ds;
        }

        /// <summary>
        /// Perform some magic to make sure the query output arrives in OML model.
        /// </summary>
        /// <param name="p">OML/SFS parameter</param>
        /// <param name="query">database query</param>
        /// <param name="valueColumn">column with values</param>
        /// <param name="IndexColumns">columns with indices</param>
        private void setBinding(Parameter p, string query, string valueColumn, string[] IndexColumns)
        {
             DataSet ds = SelectOleDbSrvRows(connection, query);
             DataTable dt = ds.Tables[0];
             p.SetBinding(dt.AsEnumerable(), valueColumn, IndexColumns);
        }

        /// <summary>
        /// Solve the problem
        /// </summary>
        public void Solve()
        {
            context.LoadModel(FileFormat.OML, new StringReader(strModel));

            foreach (Parameter p in context.CurrentModel.Parameters)
            {
                switch (p.Name)
                {
                    case "Capacity": 
                        setBinding(p,"select plant,capacity from capacity",
                            "capacity",new string[]{"plant"});
                        break;
                    case "Demand":
                        setBinding(p,"select market,demand from demand",
                            "demand", new string[]{"market"});
                        break;
                    case "Cost":
                        setBinding(p,"select plant,market,cost from cost",
                            "cost", new string[]{"plant", "market"});
                        break;
  
                }

            }

            Solution solution = context.Solve();
            Console.Write("{0}", solution.GetReport());

        }

    }
}

Some notes:

• This should work with almost any database. Just change the connection string accordingly.


• MS Access drivers are 32 bit so make sure you compile as 32 bit application. When targeting a 64 bit environment I got an exception about not being able to find an appropriate driver.


• The contents of the sets are derived from the parameter bindings: set elements are the union of the set elements used in the parameter binding.


• It may be useful to test the model beforehand using the Excel plug-in.


• Should use only one connection: make OleDbConnection conn a field of the object.


• How to write results back? This is often more complicated.

MS .NET Chart: how to turn point labels on/off

There is no property to turn on or off the visibility of data point labels (i.e. something called enabled or so). The following trick works: make the color transparent if you don’t want to see the labels. Here is some code that works for me:

        private void FlipLabels(Series s)
        {
            if (s.LabelForeColor == System.Drawing.Color.Transparent)
            {
                s.LabelForeColor = System.Drawing.Color.Black;
                s.SmartLabelStyle.CalloutLineColor = System.Drawing.Color.Black;
            }
            else
            {
                s.LabelForeColor = System.Drawing.Color.Transparent;
                s.SmartLabelStyle.CalloutLineColor = System.Drawing.Color.Transparent;
            }
        }

See also: http://social.msdn.microsoft.com/Forums/en-US/MSWinWebChart/threads/

Update: too bad I can not change the line style to arrows in a line chart. See: http://social.msdn.microsoft.com/Forums/en-US/MSWinWebChart/thread/01b61a40-203d-45dd-8834-7e8f49b974c7/.

Modeling posts today

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, 4

Job 2:
Basestock B, Brochures 2,3

Job 3:
Basestock A, Brochures 5,6,7

I 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 1

Batch 2: Jobs 2 and 3

I am looking for an algorithm that will provide the smallest number of
collections

Any 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 solely

for 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 <= n

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

Plot of Max Likelihood surface

> How can I plot the surface of a log likelihood function in GAMS.

Create a parameter that evaluates the log-likelihood function around the optimum. Lets use this example. We add the code:

set j /j1*j50/;
alias (j,jj);
scalar clo; clo = 0.8*c.l;
scalar cup; cup = 1.2*c.l;
scalar cstep; cstep = (cup-clo)/(card(j)-1);
scalar slo; slo = 0.95*sigma.l;
scalar sup; sup = 1.05*sigma.l;
scalar sstep; sstep = (sup-slo)/(card(j)-1);

scalar pc,psigma;

parameter surface(j,jj);
loop((j,jj),
   pc = clo + cstep*(ord(j)-1);
   psigma = slo + sstep*(ord(jj)-1);
   surface(j,jj) = m*log(pc) - m*pc*log(psigma)
                 + (pc-1)*sum(k,log(x(k)-theta))
                 - sum(i,((x(i)-theta)/psigma)**pc);
);
execute_unload "weibull",surface;

Now we can use Chart|Two dimensions|3D Charts|Surface.

loop madness

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

* ----------------------------------------------
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       ;
* ------------------------------------------------------

parameter a(i,j,p,k,n,q);
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).

Database Connection Strings

Somehow this is always one of the biggest hurdles to get database connectivity to work. Here is a good quick overview: http://www.connectionstrings.com/.