Monday, January 16, 2017

MIP Starts sometimes disappoint

Modern MIP solvers allow to use a good integer solution to “warm start” the branch and bound process. This looks like a very good approach to cut down on total solution times, and it sometimes is. But in many cases it turns out, finding good solutions and passing them on to the MIP solver is not worth the while.

Some reasons can be:

  • Solvers on their own are very good in finding good solutions quickly: they contain batteries of heuristics to find good integer solutions and improve them.
  • Time you spend on calculating a good initial solution can be used by the solver. If you spent n seconds to construct an initial solution, that is time taken away from the solver.
  • The solver spends most time in the tail: getting the gap down to zero (or close).

This is probably best illustrated as follows:


The blue line represents the currently best integer solution. I.e. we are minimizing here. When we provide a good integer solution we basically eat away some time from the left of the picture. But that is exactly the part where MIP solvers are very effective: the blue line decreases there very fast (the gradient with respect to time is large). In the picture above if we would have provided a starting solution with an objective value of 1000, we would essentially cut out the shaded area. The long tail at the right is where we spend most cycles.

The red line is the lower bound: the best possible objective. When we skip the shaded box we also miss out on some improvement in this bound.

It is clear that in this case, if we cannot find good solutions very cheaply ourselves, we may as well leave this to the solver.

Of course there are cases where starting with good solutions may make good sense, e.g.:

  • We can obtain good solutions quickly.
  • We may even have a feasible operating point by nature of the problem (e.g. when optimizing a certain process we may observe the physical state and use that operating point as a starting point).
  • Good solutions found outside the solver may improve reliability of the overall process: if the MIP solver fails completely we still have a good solution.

As suggested in the comments stopping on a gap percentage may be more effective in reducing the solution time:


Note that this is a depiction of quite a typical behavior of MIP lower and upper bounds: almost all non-trivial models show this type of narrowing of the gap. Often things are even more pronounced than in this smaller example. Of course stopping on gap does not likely give you the optimal solution, but in practical cases a gap of say a few percent or so may be quite acceptable. Setting a gap tolerance is a very standard approach when solving large MIP models.

Friday, January 6, 2017

New Book: R for Data Science; Visualize, Model, Transform, Tidy and Import Data

Recommended introductory text. It discusses some packages I was not really familiar with (like tibble).

Cover image


Part I. Explore

1. Data Visualization with ggplot2
2. Workflow: Basics
3. Data Transformation with dplyr
4. Workflow: Scripts
5. Exploratory Data Analysis
6. Workflow: Projects

Part II. Wrangle

7. Tibbles with tibble
8. Data Import with readr
9. Tidy with tidyr
10. Relational Data with dplyr
11. Strings with stringr
12. Factors with forcats
13. Dates and Times with lubridate

Part III. Program

14. Pipes with magrittr
15. Functions
16. Vectors
17. Iteration with purrr

Part IV. Model

18. Model Basics with modelr
19. Model Building
20. Many Models with purrr and broom

Part V. Communicate

21. R Markdown
22. Graphics for Communication with ggplot2
23. R Markdown Formats
24. R Markdown Workflow


  1. Hadley Wickham, Garrett Grolemund, R for Data Science: Import, Tidy, Transform, Visualize, and Model Data, O'Reilly Media; 1 edition (January 5, 2017), 522 pages, price: about $0.08 per page (glossy, with many color pictures) 
  2. Website:

Tuesday, January 3, 2017

Solving LP by LU???

Yes, according to this paper [link]:


Abstract This paper presents a new approach for the solution of Linear Programming Problems with the help of LU Factorization Method of matrices. This method is based on the fact that a square matrix can be factorized into the product of unit lower triangular matrix and upper triangular matrix. In this method, we get direct solution without iteration. We also show that this method is better than simplex method.

Basically the paper claims that we can solve an LP

\[\begin{align}\max\>& c^Tx\\
&Ax \le b \\
& x\ge 0\end{align}\]

as follows:

  1. Make sure \(A\) is square (by adding rows or column if needed using vacuous constraints or extra slacks).
  2. Solve \(Ax=b\) using LU decomposition
  3. Done

Wow. This is baloney of course.  But here is someone actually believing this stuff.


Of course we can use this nonsensical method and apply it to a different problem:


Abstract. In this paper the linear fractional programming problem is converted to a linear programming problem and for solving it, LU Factorization method is introduced. Finally, to illustrate the proposed method, it is verified by means of the numerical example and the results conclude that this method gives more reliable solution or the equal solution compared to the existing method.

My goodness.


  1. S.M.Chincole, A.P.Bhadane, LU Factorization Method To Solve Linear Programming Problem, International Journal of Emerging Technology and Advanced Engineering ,Volume 4, Issue 4, 176-180, Warning: this is nonsensical. Note that the example has all constraints binding in the optimal solution.
  2. R. Sophia Porchelvi, A. Anna Sheela, Solving Linear Fractional Programming Problem Using LU Factorization Method, Asian Journal of Research in Social Sciences and Humanities Vol. 6, No. 10, October 2016, pp. 1929-1938, Warning squared: astonishingly this is a paper based on (1).

Sunday, January 1, 2017

Solving Takuzu puzzles as a MIP using xor

The Takuzu puzzle (1) can be solved using a MIP model. Interestingly we can use an xor condition in the model.

Assuming an  \(n \times n\) board (with \(n\) even), we need to fill the cells \((i,j)\), with values 0 or 1, with the following restrictions:

  1. The number of ones in each row and in each column is equal to \(\displaystyle\frac{n}{2}\).
  2. Then we have also: The number of zeroes in each row and in each column is equal to \(\displaystyle\frac{n}{2}\).
  3. There can be no duplicate rows and no duplicate columns.
  4. In a row or column there cannot be more than two consecutive zeroes or ones.

A small example is given in (1):


Of course our central variable will be \(x_{i,j}\in\{0,1\}\). Condition 1 can be easily modeled using two constraints:

  &\sum_j x_{i,j} = \frac{n}{2}\>\> \forall i\\
  &\sum_i x_{i,j} = \frac{n}{2}\>\> \forall j

Condition 2 we don’t need to model explicitly: this follows from the above constraints.

Condition 3 (no duplicate rows or columns) is not very easy. One way to handle this, is to calculate a value for each row and column as follows:

&v_i = \sum_j 2^{j-1} x_{i,j}\\
&w_j = \sum_i 2^{i-1} x_{i,j}

and then make sure there are no duplicate \(v_i\)’s or duplicate \(w_j\)’s. This can be done using some all-different formulation (3). A big-M approach to implement \(v_i \le v_{i’}-1 \textbf{ or } v_i \ge v_{i’}+1, \forall i\ne i’\) can look like:

&v_i \le v_{i’}-1 + M \delta_{i,i’}\\
&v_i \ge v_{i’}+1 – M (1-\delta_{i,i’})\\
&\delta_{i,i’} \in \{0,1\}

However I want to attack this in a different way, using an xor condition (4). The constraints involved can be conceptualized as:

&\sum_j |x_{i,j} – x_{i’,j}| \ge 1 \>\> \forall i\ne i’\\
&\sum_i |x_{i,j} – x_{i,j’}| \ge 1 \>\> \forall j\ne j’

The value \(z_{i,j,i’,j’}=|x_{i,j}-x_{i’,j’}|\) can be interpreted as \(z_{i,j,i’,j’}=x_{i,j} \textbf{ xor } x_{i’,j’}\): zero if \(x_{i,j}\) and \(x_{i’,j’}\) are identical and one if they are different. This in turn can be written as a set of inequalities (4):

&z_{i,j,i’,j’}\le x_{i,j} + x_{i’,j’}\\
&z_{i,j,i’,j’}\ge x_{i,j} - x_{i’,j’}\\
&z_{i,j,i’,j’}\ge x_{i’,j’} - x_{i,j}\\
&z_{i,j,i’,j’}\le 2 - x_{i,j} - x_{i’,j’}\\
&z_{i,j,i’,j’} \in \{0,1\}

After this we just need:

&\sum_j z_{i,j,i’,j} \ge 1 \>\> \forall i\ne i’\\
&\sum_i z_{i,j,i,j’} \ge 1 \>\> \forall j\ne j’

Due to the nature of these conditions on \(z\) we can limit the xor constraints to

&z_{i,j,i’,j’}\le x_{i,j} + x_{i’,j’}\\
&z_{i,j,i’,j’}\le 2 - x_{i,j} - x_{i’,j’}\\
&z_{i,j,i’,j’} \in \{0,1\}

Note that we do not need to compute all \(z\) for all combinations \((i,j)\), \((i’,j’)\) (e.g. we never compare \(x_{2,3}\) with \(x_{3,2}\)). Furthermore we can relax \(z\) to be continuous between zero and one. Also note that if we compared row \(i\) and \(i’\) then we don’t need to compare row \(i’\) and row \(i\). Similarly for the columns. So we can do:

&\sum_j z_{i,j,i’,j} \ge 1 \>\> \forall i< i’\\
&\sum_i z_{i,j,i,j’} \ge 1 \>\> \forall j< j’

Finally we need to handle restriction 4. This can be implemented as:

&1\le x_{i,j}+x_{i,j+1}+x_{i,j+2}\le 2\\
&1\le x_{i,j}+x_{i+1,j}+x_{i+2,j}\le 2

Usually we need to write this as four inequalities (unless we can specify ranged equations). If the summation was longer I would have introduced extra variables with bounds one and two to reduce the number of constraints.

Note that the lower bound of one can also be inferred from \((1-x_{i,j})+(1-x_{i,j+1})+(1-x_{i,j+2})\le 2 \Rightarrow x_{i,j}+x_{i+1,j}+x_{i+2,j} \ge 1\).

Complete model

The variables are:


The dummy variable is related to a dummy objective: we look for a feasible integer solution only.

The different \(z\) values we need to compare is encoded in a set compare:


Finally the model equations are:


The result looks like:

----     77 VARIABLE x.L 

            i1          i2          i3          i4

i1                       1           1
i2           1                                   1
i3                                   1           1
i4           1           1

A large problem

The site (2) contains some “very hard” 14 x 14 puzzles. An example is:


It is interesting to see how this model behaves on this problem. It turns out the problem can be solved completely in preprocessing. E.g. Cplex and Gurobi report that zero Simplex iterations were needed and solve this problem blazingly fast:

Tried aggregator 3 times.
MIP Presolve eliminated 5491 rows and 2301 columns.
MIP Presolve modified 18 coefficients.
Aggregator did 396 substitutions.
Reduced MIP has 92 rows, 48 columns, and 298 nonzeros.
Reduced MIP has 48 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.03 sec. (7.10 ticks)
Found incumbent of value 0.000000 after 0.05 sec. (8.15 ticks)

Root node processing (before b&c):
  Real time             =    0.05 sec. (8.24 ticks)
Sequential b&c:
  Real time             =    0.00 sec. (0.00 ticks)
Total (root+branch&cut) =    0.05 sec. (8.24 ticks)
MIP status(101): integer optimal solution
Cplex Time: 0.05sec (det. 8.25 ticks)

Verifying Uniqueness

A well-formed puzzle has only one solution. To verify whether a solution is unique we can use a similar strategy as in (6):

  1. First solve the model to obtain a solution \(x^*_{i,j}\).
  2. Formulate a constraint to forbid this solution.
  3. Add this constraint to the model and resolve. The solver should not find a new solution but rather report that this model is infeasible.

Here the cut can be formulated as:

\[\sum_{i,j} x^*_{i,j} x_{i,j} \le \frac{n^2}{2}-1\]
Data Entry

Entering the problem data:


is an interesting issue. In GAMS “zero” and “does not exist” is identical, so we can not just use:

table init(i,j)
i1 i2 i3 i4

i1     1     0
i2        0
i3     0
i4  1  1     0


I used the following, which is not very clean:

table init(i,j)
i1 i2 i3 i4

i1     1     2
i2        2
i3     2
i4  1  1     2

x.fx(i,j)$(init(i,j)=1) = 1;
x.fx(i,j)$(init(i,j)=2) = 0;

Not sure if this is better:

parameter init(i,j) / #i.#j NA /;
i1 i2 i3 i4

i1     1     0
i2        0
i3     0
i4  1  1     0
) = init(i,j);

Note that AMPL has a "default" facility, so we can use in a data file:

param init default -1 : i1 i2 i3 i4 :=
                   i1    .  1  .  0
                   i2    .  .  0  .
                   i3    .  0  .  .
                   i4    1  1  .  0

It makes sense to use Excel for data entry. The initial configuration in Excel can look like:


To read this into GAMS and retain the zeros we can use the following code:

parameter init(i,j);

$call 'gdxxrw i=takuzu.xlsx sq=n par=init rng=b3 rdim=1 cdim=1'

$gdxin takuzu.gdx
$loaddc init


The Squeeze=no option in the call to GDXXRW will tell it to retain the zeros. The $OnEps command tells GAMS to read the zero as an EPS special value. EPS is very useful: it is physically stored in the sparse data structures, but as soon as you do a numerical operation on it, it becomes zero. This will show:

----     12 PARAMETER init 

            i1          i2          i3          i4

i1                   1.000                     EPS
i2                                 EPS
i3                     EPS
i4       1.000       1.000                     EPS

This will now allow us to fix variables as:

binary variable x(i,j);
x.fx(i,j)$init(i,j) = init(i,j);

This model turned out to have quite a few interesting angles.

  1. Takuzu,
  2. Online Takuzu puzzles,
  3. All-different and Mixed Integer Programming,
  4. xor as Inequalities,
  5. Similar puzzles are Sudoku and KenKen:
  6. Unique Solutions in KenKen,

Saturday, December 31, 2016

Smart Slime Mold on page 3 of the Washington Post



  1. Washington Post article:
  2. Slime mold and linear programming:

MIP Modeling

In an investment planning model I needed to model the following:


We can build a facility during any time period. So we have:

\[\begin{align}&build_t \in \{0,1\} \\
                   &\sum_t build_t \le 1\end{align}\]

The variable \(open_t\in \{0,1\}\) indicates we have an open facility, ready to do business. A facility is considered open after it is built (we don’t close facilities during the planning horizon). To be more precise: it becomes available one period after it has been built.

There are actually a few ways to model this:

alternative 1

Look if the facility was built in any period \(t’<t\):

\[open_t = \sum_{t’|t’<t} build_{t’}\]
alternative 2

Use an ‘or’ formulation:

\[open_t = build_{t-1} \textbf{ or } open_{t-1}\]

This can be linearized as:

\[\begin{align}&open_t \ge build_{t-1}\\
                    &open_t \ge open_{t-1}\\
                    &open_t \le build_{t-1}+open_{t-1}\end{align}\]

alternative 3

We can simplify the recursive construct:

\[open_t = build_{t-1}+open_{t-1}\]

With this, we now explicitly forbid \(build_{t-1}+open_{t-1}=2\) which was not allowed (implicitly) anyway.

The last formulation seems preferable in terms of simplicity and the number of non-zero elements needed.

As usual with lags we need to be careful what happens near the borders, in this case when \(t=t1\). We just can fix \(open_{t1}=0\). When using GAMS you can even skip that, as indexing with lags outside the domain returns zero. Here that means \(open_{t1} = build_{t1-1} + open_{t1-1} = 0 + 0\). We can also verify this in the equation listing:

---- eqopen2  =E= 

eqopen2(t1)..  open(t1) =E= 0 ; (LHS = 0)
eqopen2(t2)..  - build(t1) - open(t1) + open(t2) =E= 0 ; (LHS = 0)
eqopen2(t3)..  - build(t2) - open(t2) + open(t3) =E= 0 ; (LHS = 0)

Wednesday, December 28, 2016

Unique Solutions in KenKen

In (1) a MIP model is proposed so solve the KenKen puzzle. During a discussion, the question came up if I could prove the uniqueness of a solution. In the Mixed Integer Programming model I used a standard formulation for a solution: 

\[x_{i,j,k} = \begin{cases}1 & \text{if cell $(i,j)$ has the value $k$}\\
                                    0 & \text{otherwise}\end{cases}\]

A general approach could be to use the technique described in (2): add a cut to forbid the current solution and solve again. If this second solve is infeasible we have established that the original solution was unique.

In this case we can use a more specialized cut that is simpler:

\[\sum_{i,j,k} x^*_{i,j,k} x_{i,j,k} \le n^2-1\]

where \(x^*\) is the previous solution and \(n \times n\) is the size of the puzzle.

To test this with the model and problem data shown in (1) I used:


Note that \(\displaystyle\sum_{i,j,k|x^*_{i,j,k}=1} x_{i,j,k}\) is identical to \(\displaystyle\sum_{i,j,k} x^*_{i,j,k} x_{i,j,k}\). To make sure things work correctly with solution values like 0.9999, I actually used a somewhat generous tolerance: \(\displaystyle\sum_{i,j,k|x^*_{i,j,k}>0.5} x_{i,j,k}\).

Indeed the solution from the first solve was unique. The second solve yielded:

               S O L V E      S U M M A R Y

     MODEL   kenken2             OBJECTIVE  z
     TYPE    MIP                 DIRECTION  MINIMIZE
     SOLVER  CPLEX               FROM LINE  115

**** SOLVER STATUS     1 Normal Completion        
**** MODEL STATUS      10 Integer Infeasible      
**** OBJECTIVE VALUE               NA

RESOURCE USAGE, LIMIT          0.031      1000.000
ITERATION COUNT, LIMIT         0    2000000000

This approach can also be applied to the Sudoku MIP model.

  1. KenKen puzzle solved using a MIP model:
  2. Forbid a given 0-1 solution:

Tuesday, December 20, 2016

SAS/OR MIP solver

In (1) the performance of the SAS/OR MIP solver is discussed. It uses the MIPLIB problem set. Timings from different studies are always difficult to compare, but if we take the 214 or 218 feasible problems, and measure how many of those problems could be solved within in 2 hours we see:

Solver Number solved within 2 hours (12 threads) Source
CBC 119 of 218 (3)
SAS/OR 170 of 214 (1)
Cplex 204 of 218 (3)
Gurobi 210 of 218 (3)

I think it is fair to guess that the SAS solver is faster than CBC but a little bit slower than the market leaders Cplex and Gurobi. Again: I don’t know the difference in speed between the machines these benchmarks were executed on, so take this with a grain of salt.

I have quite a few clients with large SAS licenses, but they typically do not use SAS/OR. As a result I have never used the SAS/OR solvers. I usually use SAS as a database. My typical workflow is: export the SAS data, combine with other data, develop and solve the models with other tools and report solutions in Excel.

  1. Imre PĆ³lik, How good is the MILP solver in SAS/OR?, December 2016,
  2. Philipp M. Christophel e.a., The SAS MILP Solver: Current Status and Future Development, INFORMS 2016 meeting, 
  3. H. Mittelmann, The solvable MIPLIB instances (MIPLIB2010),

Sunday, December 18, 2016

MonetDB and R


MonetDB (1) is a so-called column-oriented database or column-store (2). It is server based, SQL compliant and open-source. The column orientation can lead to better performance for certain types of tasks, especially OLAP and ETL (i.e. analytical work). Traditional row-wise databases are said to be more appropriate for OLTP workloads.  


There exists a CRAN package to let R talk to a MonetDB server (MonetDB.R). There is also a different package called MonetDBLite. This contains an in-process version of MonetDB.  This means MonetDBLite is more or less an alternative for RSQLite (4). A picture comparing the performance of MonetDBLite is from (4):

Basically the more towards the left-lower corner the better.

In (3) there are some timings comparing MonetDBLite to SQLite. E.g. converting a (large) table to a data frame:

There is lots of data traffic from the database to R and the server based MonetDB.R does not finish the larger tables (traffic has to go through a TCP/IP socket). But the embedded database is very fast compared to SQLite.

This looks like an interesting database for R related work.

  1. MonetDB The column-store pioneer,
  2. Daniel Adabi, Distinguishing Two Major Types of Column-Stores, March 2010,
  4. Anthony Damico, MonetDBLite because fast, June 2016,
  5. Using MonetDB[Lite] with real-world CSV files,

Saturday, December 17, 2016

New book: Optimization by GRASP; Greedy Randomized Adaptive Search Procedures

This is a well known meta-heuristic (with a somewhat less flashy name than some of the new heuristics that get invented every week). Apparently this is the first book fully dedicated to GRASP. GRASP was introduced in 1989 by Feo and Resende (2).

  1. Introduction
  2. A short tour of combinatorial optimization and computational complexity
  3. Solution construction and greedy algorithms
  4. Local search
  5. GRASP: The basic heuristic
  6. Runtime distributions
  7. Extended construction heuristics
  8. Path-relinking
  9. GRASP with path-relinking
  10. Parallel GRASP heuristics
  11. GRASP for continuous optimization
  12. Case studies

  1. Mauricio G. C. Resende and Celso C. Ribiero, Optimization by GRASP,  Springer 2016. Price: about $0.25 per page (some of them with color pictures).
  2. T.A. Feo and M.G.C. Resende, A probabilistic heuristic for a computationally difficult set covering problem, Operations Research Letters, 8:67–71, 1989.