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

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

#### Contents

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

#### References

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: http://r4ds.had.co.nz/ ## 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. ##### Update 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. #### References 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, http://www.ijetae.com/files/Volume4Issue4/IJETAE_0414_31.pdf. 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, https://www.aijsh.com/shop/articlepdf/2016/10/1475416103154.pdf. 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:  \begin{align} &\sum_j x_{i,j} = \frac{n}{2}\>\> \forall i\\ &\sum_i x_{i,j} = \frac{n}{2}\>\> \forall j\end{align} 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:  \begin{align}&v_i = \sum_j 2^{j-1} x_{i,j}\\&w_j = \sum_i 2^{i-1} x_{i,j}\end{align} 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:  \begin{align}&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\}\end{align} However I want to attack this in a different way, using an xor condition (4). The constraints involved can be conceptualized as:  \begin{align}&\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’\end{align} 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):  \begin{align}&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\}\end{align} After this we just need:  \begin{align}&\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’\end{align} Due to the nature of these conditions on $$z$$ we can limit the xor constraints to  \begin{align}&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\}\end{align} 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:  \begin{align}&\sum_j z_{i,j,i’,j} \ge 1 \>\> \forall i< i’\\&\sum_i z_{i,j,i,j’} \ge 1 \>\> \forall j< j’\end{align} Finally we need to handle restriction 4. This can be implemented as:  \begin{align}&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\end{align} 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 1i2 1 1i3 1 1i4 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 solutionCplex 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 /;$onmulti
table
init(i,j)

i1 i2 i3 i4

i1     1     0

i2        0

i3     0

i4  1  1     0
;
variable
x(i,j);
x.fx(i,j)$(init(i,j)<> NA ) = 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'

$oneps$gdxin takuzu.gdx
$loaddc init display 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                     EPSi2                                 EPSi3                     EPSi4       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. ##### References 1. Takuzu, https://en.wikipedia.org/wiki/Takuzu 2. Online Takuzu puzzles, http://www.binarypuzzle.com/ 3. All-different and Mixed Integer Programming, http://yetanothermathprogrammingconsultant.blogspot.com/2016/05/all-different-and-mixed-integer.html 4. xor as Inequalities, http://yetanothermathprogrammingconsultant.blogspot.com/2016/02/xor-as-linear-inequalities.html 5. Similar puzzles are Sudoku and KenKen: http://yetanothermathprogrammingconsultant.blogspot.com/2016/10/mip-modeling-from-sudoku-to-kenken.html 6. Unique Solutions in KenKen, http://yetanothermathprogrammingconsultant.blogspot.com/2016/12/unique-solutions-in-kenken.html ## Saturday, December 31, 2016 ### Smart Slime Mold on page 3 of the Washington Post #### References ### 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’ ##### 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) REMAINING 3 ENTRIES SKIPPED ## 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.000ITERATION COUNT, LIMIT 0 2000000000 This approach can also be applied to the Sudoku MIP model. ##### References 1. KenKen puzzle solved using a MIP model: http://yetanothermathprogrammingconsultant.blogspot.com/2016/10/mip-modeling-from-sudoku-to-kenken.html 2. Forbid a given 0-1 solution: http://yetanothermathprogrammingconsultant.blogspot.com/2011/10/integer-cuts.html ## 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. ##### References 1. Imre Pólik, How good is the MILP solver in SAS/OR?, December 2016, http://blogs.sas.com/content/operations/2016/12/19/good-milp-solver-in-sasor/ 2. Philipp M. Christophel e.a., The SAS MILP Solver: Current Status and Future Development, INFORMS 2016 meeting, http://blogs.sas.com/content/operations/files/2016/12/Christophel_INFORMS2016.pdf 3. H. Mittelmann, The solvable MIPLIB instances (MIPLIB2010), http://plato.asu.edu/ftp/solvable.html ## 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. ##### MonetDBLite 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. ##### References 1. MonetDB The column-store pioneer, https://www.monetdb.org/Home 2. Daniel Adabi, Distinguishing Two Major Types of Column-Stores, March 2010, http://dbmsmusings.blogspot.com/2010/03/distinguishing-two-major-types-of_29.html 3. https://www.monetdb.org/blog/monetdblite-r 4. Anthony Damico, MonetDBLite because fast, June 2016, https://www.r-bloggers.com/monetdblite-because-fast/ 5. Using MonetDB[Lite] with real-world CSV files, https://www.r-bloggers.com/using-monetdblite-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). Chapters: Introduction A short tour of combinatorial optimization and computational complexity Solution construction and greedy algorithms Local search GRASP: The basic heuristic Runtime distributions Extended construction heuristics Path-relinking GRASP with path-relinking Parallel GRASP heuristics GRASP for continuous optimization Case studies ##### References 1. Mauricio G. C. Resende and Celso C. Ribiero, Optimization by GRASP, Springer 2016. Price: about0.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.