In http://rpubs.com/KGFeMan/PM_IP an example of a cutting stock problem in solved in R using LpSolve. In this example all possible cutting patterns are generated in advance using some facilities of the iterpc package. A different way is to generate interesting patterns on the fly using a Column Generation scheme (see (1)).
I was present at a GAMS demonstration yesterday, and I saw a strange thing. The model is a very small transportation model (model 1 in the GAMS model library):
| \[\boxed{\begin{align}\min\>&\sum_{i,j} c_{i,j} x_{i,j}& \\ &\sum_j x_{i,j} \le a_i & \forall i \\ &\sum_i x_{i,j} \ge b_j & \forall j\\ &x_{i,j}\ge 0 \end{align}}\] |
The twist is to make \(x_{i,j}\) semi-continuous: \(x_{i,j}=0\) or between \([100,1000]\). This can be modeled as:
| \[\boxed{\begin{align} & 100 y_{i,j}\le x_{i,j} \le 1000 y_{i,j} \\ & y_{i,j} \in \{0,1\} \end{align}}\] |
When looking at the log we saw:
| Solution satisfies tolerances. MIP Solution: 159.525000 (3 iterations, 0 nodes) Best possible: 153.675000 | I expect the highlighted numbers to be the same |
Some background : by default GAMS solves a MIP with relative gap tolerance of 10%, so the MIP solution of 159.525 is not globally optimal but rather within 10% of the best possible solution. As GAMS wants to report marginals (duals and reduced costs) for a MIP, all integer variables are fixed and then the problem is resolved as an LP. This is the “Final Solve”. One would expect the same objective as the MIP solution, but here this is not the case. What is happening? Is this a bug?
Let’s compare the solutions.
| MIP Solution (MIP solution found within gap tolerance) ---- 78 VARIABLE y.L new-york chicago topeka seattle 1.000
new-york chicago topeka seattle 325.000
---- 78 VARIABLE y.L new-york chicago topeka seattle 1.000
new-york chicago topeka seattle 325.000
|
We see the discrete variables (\(y\)) are the same but one of the continuous variables (\(x\)) has changed.
Hypothesis: my guess is that the solution with \(z=159.525\) was found using a heuristic instead of by solving this node using the standard dual simplex method. As a result this node may not be LP optimal. The log gives some hints:
| Iteration Dual Objective In Variable Out Variable Iteration Dual Objective In Variable Out Variable Nodes Cuts/ * 0+ 0 159.5250 12.6000 92.10% Nodes Cuts/ 0 0 153.6750 1 159.5250 153.6750 3 3.67% |
Indeed it seems it took zero LP iterations to find the solution of 159.525 possibly by applying some heuristics on the root node solution. Somehow Cplex never solved that node to optimality after applying the heuristic, so we don’t see the 156.375 solution.
Now is this a bug? At first sight it is not. However an integer solution that is not LP optimal is actually often nonsensical. In the above example, the MIP solutions just ships too much to Topeka:
| ---- 86 PARAMETER sol new-york chicago topeka capacity shipped seattle 325.000 350.000 325.000 |
Obviously this solution does not pass the laugh test. I can see presenting such a solution to your boss may get you into trouble. As a result, I think reported MIP solutions should be LP optimal to prevent such meaningless solutions.
So this GAMS strategy of doing this final solve actually has a good side effect: cleaning up integer solutions. If you don’t use GAMS and use a gap tolerance: recheck your model solutions, because there may be issues with the continuous variables being non-optimal.
Note: this experiment was done with GAMS 24.7.3 and Cplex 12.6.3.0.
Conclusion: when using a gap tolerance Cplex can return solutions that are non-optimal with respect to the continuous variables. This is probably a bug. Although these solutions are (integer) feasible they can have undesirable characteristics. The way to repair these solutions is to fix the integer variables and resolve the model as an LP.
Many problems solved by dynamic programming can also be formulated as a Mixed Integer Programming problem. Here is a very simple problem. The participant is asked to provide a Dynamic Programming algorithm that solves this problem to optimality. A MIP approach would be even simpler. Given \(n\) gold mines located along a river. Each mine has a location \(x_i\) (measured along the river) and has \(w_i\) tons of gold. The goal is to select a subset of of the mines as pick-up points. The number of pick-up points is given: \(k\). The cost to move gold from mine \(i\) to pick-up point \(j\) is \(c_{i,j}=w_i |x_i-x_j|\).
This actually looks very much like an assignment problem:
| \[\boxed{\begin{align} \min\>&\sum_{i,j} c_{i,j} x_{i,j}& \\ & \sum_j x_{i,j} = 1 &\forall i\\& \sum_i x_{i,j} = 1 &\forall j\\ & x_{i,j}\in \{0,1\} \end{align}}\] |
We need to allow a destination \(j\) to serve multiple sources \(i\), so first we drop the second constraint. The second issue is that we need to count the number of destinations that are used. We can do this by introducing a binary variable \(y_j\) with \(\sum_j y_j=k\). If \(y_j=0\) we need to ensure \(x_{i,j}=0\). This can simply be modeled as \(x_{i,j}\le y_j\). The complete model can now look like:
| \[\boxed{\begin{align} \min\>&\sum_{i,j} c_{i,j} x_{i,j} & \\ & \sum_j x_{i,j} = 1 & \forall i\\& x_{i,j}\le y_j & \forall i,j\\ &\sum_j y_j = k \\& x_{i,j},y_j \in \{0,1\} \end{align}}\] |
The solution for the example looks like
| ---- 42 PARAMETER c cost of moving i to j mine1 mine2 mine3 mine1 10.000 20.000
mine2 mine1 1.000
mine2 1.000 |
Notes:
Some more exotic meta-heuristics with intriguing names:
See also: http://yetanothermathprogrammingconsultant.blogspot.com/2016/05/more-nature-inspired-heuristics.html
In this post a somewhat complex class of functions was incorrectly assumed to be convex:
(From: http://www.sfu.ca/~ssurjano/boha.html, the plot is \(f_1(x)\)).
Indeed the picture gives the impression this is a simple bowl shaped thing and should be easy to minimize. Most functions with a cosine term are not convex, and indeed we can see more details here:
(From: http://al-roomi.org/benchmarks/unconstrained/2-dimensions/11-bohachevsky-s-function-no-2; this plot is \(f_2(x)\)).
All the functions have a global optimal solution \(x^*_1=0\), \(x^*_2=0\), \(f^*=0\).
It is not so difficult to solve this with a good local solver: use a simple multi-start strategy. I.e. using the NLP solver CONOPT and 20 trials with a random starting point \(x^0_1,x^0_2 \in [-100,100]\) gives:
| ---- 27 PARAMETER xsol x.1 x.2 f.obj trial1 0.6186 1.35041E-22 0.4129 |
(blank means zero in this GAMS output). There was no problem in finding the global optimum \(x^*_1=0\) , \(x^*_2=0\) , \(f^*=0\) here.
Not all solvers support the \(cos(x)\) function:
It would be nice if we had a tool that could detect convexity. A modeling system like AMPL and GAMS has knowledge about the analytical form of the optimization problem, so in theory it could try to find out if the problem is convex. This is not a trivial task (see (1)). A global solver like Baron does in essence a similar thing before solving the problem. A more heuristic approach is used in (2).
The Huffington Post is now even reporting on optimization models:
 | In this post I described a simple way to export a large GAMS data set to a CSV file which can be used as an external data source by Pivot Tables. Basically this was needed as Excel has a limit of a million rows, so we cannot store the data in Excel itself as a table. However newer versions of Excel have a facility called PowerPivot that allows us to store a very large data (or even multiple tables) in what is called a “data model”. This is an internal memory based database (“in-memory analytics engine”), and it does not have a strict row limit. Subsequently, we can create a Pivot Table based on this data. Here is how to set this up interactively. |
The limits are discussed here.
The .xlsx file containing the data can become large. After reading a the above cube.csv file (63 MB) the resulting .xlsx file is only 3 MB. So the compression may actually be doing very well.
The row limit of a table in the data model seems to be: 1,999,999,997, well beyond the row limit of Excel sheets.
In this post I showed how a small non-convex QP could be reformulated as a MIP using KKT conditions plus a special (linear) objective. Will this work in practice? How does this compare performance-wise with global solvers? Let’s see if we can get a feeling for this.
One easy way to make the small test problem bigger is to duplicate the problem:
| Objective | \(c'x_1 + x'_1 Q x_1\) | \(+c'x_2 + x'_2 Q x_2\) | \(\cdots\) | \(+c'x_K + x'_K Q x_K\) |
| Constraints | \(A x_1 \le b\) \(\ell \le x_1 \le u\) | |||
| \(A x_2 \le b\) \(\ell \le x_2 \le u\) | ||||
| \(\ddots\) | ||||
| \(A x_K \le b\) \(\ell \le x_K \le u\) |
Basically we need to add an index K to the variables \(x_i\), and to the constraints. Finally we add a summation to the objective. The new model looks like:
For the KKT model we can do something similar.
Here are some results with K=5 and K=10. I stopped after 1000 seconds (if not proven optimal by then, I report the gap).
Note: The global QP solver GloMiQo behaves the same as Antigone.
At least for this model the MIP formulation is doing very well. The Cplex MIP solver beats the Cplex global QP solver by a large margin!
Yalmip has reported similar results produced by its KKTQP method.
The Cplex global QP solver has some problems providing marginals (duals and reduced costs) for a non-convex QP. Here is a small model taken from (1):
To solve this problem to optimality in GAMS we add ‘option optcr=0’ and we provide a Cplex option ‘OptimalityTarget = 3’.
The solution looks like:
| LOWER LEVEL UPPER ---- EQU Obj . . . Obj objective function ---- VAR x LOWER LEVEL UPPER i1 . 1.0000 1.0000 LOWER LEVEL UPPER ---- VAR f -INF -17.0000 +INF |
The marginals are missing. Indeed the log contains message:
| *** No duals possible for models with quadratic objective when solved globally. |
The simplest way to get the marginals of this problem is to use a global solver that provides duals:
These solvers do not provide (correct) duals:
One way to get marginals with a global solver that does not provide them is resolving the problem with a local solver. As in:
For this model it takes Cplex 28 iterations to find the global solution, and then Conopt needs 4 iterations to recover the duals. The results look like:
| LOWER LEVEL UPPER MARGINAL ---- EQU Obj . . . 1.0000 Obj objective function ---- VAR x LOWER LEVEL UPPER MARGINAL i1 . 1.0000 1.0000 -58.0000 LOWER LEVEL UPPER MARGINAL ---- VAR f -INF -17.0000 +INF . |
An alternative is to solve the non-convex QP as a MIP using a model with KKT conditions and a special objective (see here).
(1) Floudas e.a., Handbook of Test Problems in Local and Global Optimization, Kluwer, 1999.
