CLP is not doing so poorly on this problem (note the log/log scale however).
Monday, June 6, 2016
COIN-OR CLP vs Cplex on assignment problem
http://orinanobworld.blogspot.com/2016/06/using-clp-with-java.html
CLP is not doing so poorly on this problem (note the log/log scale however).
CLP is not doing so poorly on this problem (note the log/log scale however).
For rent: PC with 2TB of RAM
https://aws.amazon.com/ec2/instance-types/x1/#
Just 4 bucks an hour. I guess not really meant to run MS Word.
But what is the difference between a GiB and a GB? I had to look this up here.
Just 4 bucks an hour. I guess not really meant to run MS Word.
Friday, June 3, 2016
R Markdown
I am starting to get into the habit to write R code directly in R markdown. I like it. Here is an excellent cheat sheet. I noticed I need to slightly change the point sizes in my ggplot plots:
Monday, May 30, 2016
All-different and mixed integer programming
There are quite a few models that use the ‘all-different’ constraint, i.e. a set of integer variables \(x_i\), \(i=\{1,...,n\}\) is feasible only if they are all different, i.e. \(x_i\ne x_j\) for \(i\ne j\). We can probably say that most of these models are of the “educational” type and are not practical production models.
Constraint programming solvers have typically a built-in “all-different” global constraint. This makes it easy to write down such a constraint and the solver will have knowledge about the constraint which it can exploit. That means we can expect better performance than for instance a bunch of pairwise not-equal constraints. So the first observation is: if you have a model that is largely built around an all-different constraint, consider to implement the model using a constraint programming solver.
Approach 1
One way to implement this construct for use in a MIP model is to use pairwise comparison: \(x_i\ne x_j\) for \(i\ne j\). In a MIP we need binary variables for this:
| (I) | \[\boxed{\begin{align}&x_i \le x_j - 1 + M^{(1)}_{i,j}\delta_{i,j}\\ &x_i \ge x_j + 1 – M^{(2)}_{i,j}(1-\delta_{i,j})\\ &\delta_{i,j} \in \{0,1\} \end{align}}\] | for all \(i\lt j\) |
Note that we only need to compare \(x_i\) and \(x_j\) if \(i<j\). This means we need \(\frac{n(n-1)}{2}\) binary variables. It is important to choose good values for \(M\) so let’s work on that for a minute (see here for an example where things go wrong if we don’t pay attention to this). Note that instead of using a single \(M\) we use different values \(M^{(1)}_{i,j}\) and \(M^{(2)}_{i,j}\).
The value of \(M^{(1)}_{i,j}\) should be chosen as small as possible subject to \(x_j-1+M^{(1)}_{i,j}\ge x^{up}_i\). This means \(M^{(1)}_{i,j}\ge x^{up}_i +1 – x_j\). This yields the following optimal value: \(M^{(1)}_{i,j} = x^{up}_i +1 – x^{lo}_j\). Similarly we want to choose the smallest value \(M^{(2)}_{i,j}\) such that \(x_j+1-M^{(2)}_{i,j} \le x^{lo}_i\). This gives us: \(M^{(2)}_{i,j} = x^{up}_j +1 - x^{lo}_i\).
Approach 2
Here we consider the special case where each \(x_i \in \{1,…,n\}\) where \(i=\{1,..,n\}\). Now we can write:
| (II) | \[\boxed{\begin{align}&\sum_j \delta_{i,j} = 1 \> \forall i\\ &\sum_i \delta_{i,j} = 1 \> \forall j \\ &x_i = \sum_j j \delta_{i,j} \\ &\delta_{i,j} \in \{0,1\} \end{align}}\] |
The matrix \(\delta\) can be looked at as a permutation matrix. This permutation matrix \(\delta\) is a row- and column-permuted identity matrix, i.e. it has a single entry equal to one in each row and in each column. Another way of looking at the first two equations is as an assignment problem. Note that we have \(n^2\) binary variables here, but no big-M’s.
We can make one extra simplification: we can drop the first assignment constraint. The resulting equations are:
| (III) | \[\boxed{\begin{align}&\sum_j j \delta_{i,j} = x_i\> \forall i\\ &\sum_i \delta_{i,j} = 1 \> \forall j \\ &\delta_{i,j} \in \{0,1\} \end{align}}\] |
This simplification in only possible in this special case. We see below that some minor extensions make this simplification fail.
Extension 1
The second approach was tailored to a very special case. Instead of \(x_i \in \{1,…,n\}\) where \(i=\{1,..,n\}\) now consider a slightly more general problem, where \(x_i \in \{a_1,…,a_n\}\) with \(a_{i+1} \gt a_i\) (that is: no duplicates in the set). This problem is easily handled by a simple extension to model II:
| (IIa) | \[\boxed{\begin{align}&\sum_j \delta_{i,j} = 1 \> \forall i\\ &\sum_i \delta_{i,j} = 1 \> \forall j \\ &x_i = \sum_j a_j \delta_{i,j} \\ &\delta_{i,j} \in \{0,1\} \end{align}}\] |
Can we use the same trick as in model III? The answer is no. For example take \(x_i \in \{0,1,3,5\}\). Then a model using:
| (IIIa) | \[\boxed{\begin{align}&\sum_j a_j \delta_{i,j} = x_i\> \forall i\\ &\sum_i \delta_{i,j} = 1 \> \forall j \\ &\delta_{i,j} \in \{0,1\} \end{align}}\] | Incorrect! |
Extension 2
We can allow explicit duplicates by saying: \(x_i \in \{a_1,…,a_n\}\) with \(a_{i+1} \ge a_i\). E.g. the data array \(a=\{1,2,2,3,3,3\}\) would allow 1 one, 2 twos, and 3 threes. This can be handled directly by model IIa.
Extension 3
Allow a subset. I.e. consider \(x_j \in \{a_1,…,a_n\}\) where \(j=\{1,..,m\}\) and \(m<n\). For example choose two different values \(x_j\) from the set \( \{1,2,3\}\). To model this let’s use \(i\) as the index of \(a_i\). So we have \(j\) being a subset of \(i=\{1,…,n\}\). We need to revise model II as follows:
| (IIb) | \[\boxed{\begin{align}&\sum_j \delta_{i,j} \le 1 \> \forall i\\ &\sum_i \delta_{i,j} = 1 \> \forall j \\ &x_j = \sum_i a_i \delta_{i,j} \\ &\delta_{i,j} \in \{0,1\} \end{align}}\] |
If the range of integer variables \(x_j\) is somewhat small (i.e. we don't have \(x^{up}_j \gg x^{lo}_j\)), we can solve the same problem as handled by approach 1.
References
- H.P. Williams, Hong Yan, "Representations of the all-different Predicate of Constraint Satisfaction in Integer Programming," INFORMS Journal on Computing, Vol. 13 (2001) 96-103
- W.J. van Hoeve, “The alldifferent Constraint: A Survey,” 6th Annual workshop of the ERCIM Working Group on Constraints, 2001, [link]
- http://yetanothermathprogrammingconsultant.blogspot.com/2016/10/mip-modeling-from-sudoku-to-kenken.html
Saturday, May 28, 2016
MIP Modeling
From this post:
Imagine a village with people trading goods. Each person has his own offer in this format: Giving amount
aof goodbfor amountcof goodd.I have made a simple table to point out what I am trying to explain:
In this case there are three different goods: Wood, Sand and Gras.
I also live in the village and noticed that the prices of the traders vary greatly. I have 1 wood and want to increase it by simply trading between the five dealers. But: I must not visit a dealer more than once.
There would be different routes, for example I could do
Dave - Adam, which would result in +1 wood for me. A better route would beEarl - Berta - Adam, because it would mean +2 wood.
It is possible to model this with a linear Mixed Integer Programming (MIP) model. I tried this quickly. The optimal solution is:
| ---- 65 VARIABLE x.L trade t1 t2 t3 t4 Adam 1
t1 t2 t3 t4 Adam 4
t1 t2 t3 t4 t5 wood 4 4 |
The complete MIP model is below. It is built around an assignment problem structure, which makes sure we visit at most one dealer each period, and that we do not revisit dealers.
Extensions:
- The data for a trades can be fractional (1.3 grass for 2.4 wood).
- The multiples of a trade are in the model restricted to integers. We can easily allow any trade size > 0 by making y a continuous variable. We may need to change equation xy2 to \(y_{p,t} \ge 0.01 x_{p,t}\).
- We can easily apply limits on each trade, in the form of a capacity constraint. We already applied a fixed upper bound of 100 on each trade. This can be made part of the input instead (and can be made different for each trade).
Reading a spreadsheet
I want to import the data for a three-dimensional parameter p(i,j,k) that is stored in in k excel sheets but GAMS does not let me use dollar control statements in loops. Is there any way to do that using loops or other flow control statements like 'for' or 'while'?
Let’s make some data as follows:
 |  |  |
We can read this as follows:
| $set xls d:\tmp\test2.xlsx |
The output will look like:
| ---- 23 PARAMETER a all data Sheet1 Sheet2 Sheet3 i1.j1 1.000 2.000 3.000 |
It is quite slow however as we do a call to gdxxrw for each sheet (we usually prefer to read all data into a single GDX file using just one invocation of gdxxrw).
Friday, May 27, 2016
More nature-inspired heuristics
In the book
a large number of meta-heuristics based loosely on Nature are described. A list is: Simulated Annealing, Genetic Algorithms, Differential Evolution, Particle Swam Optimization, Firefly Algorithms, Cuckoo Search, Bat Algorithms, Flower Pollination Algorithms, But there are actually more to be found:
(Updated to reflect comments)
a large number of meta-heuristics based loosely on Nature are described. A list is: Simulated Annealing, Genetic Algorithms, Differential Evolution, Particle Swam Optimization, Firefly Algorithms, Cuckoo Search, Bat Algorithms, Flower Pollination Algorithms, But there are actually more to be found:
- Moth-flame optimization
- Whale Optimization
- Sine Cosine Algorithm
- Dragonfly Algorithm
- Multi-Verse Optimization
- Ant Lion Optimizer
- Grey Wolfe Optimizer
- Slime Mold Solver (apparently the favorite of Paul Rubin)
- and many more
I am note sure if the development of all these different methods is a good sign. How many should one consider or try out on a given problem?
A small collection of more interestingly named heuristics is here: http://yetanothermathprogrammingconsultant.blogspot.com/2016/08/egyptian-vultures-leaping-frogs-and.html.
A small collection of more interestingly named heuristics is here: http://yetanothermathprogrammingconsultant.blogspot.com/2016/08/egyptian-vultures-leaping-frogs-and.html.
References
A useful paper with a similar point of view:Kenneth Sörensen, “Metaheuristics—the metaphor exposed,” International Transactions in Operational Research, Volume 22, Issue 1, January 2015, Pages 3–18
(Updated to reflect comments)
Monday, May 23, 2016
Big-M to the extreme
A major issue in MIP modeling is choosing good values for big-M constants. The poorly chosen name 'big-M' indicates we should use a really big value which is exactly the opposite of what we should do, Here we see an extreme example:
It is my conjecture that just because of this name 'big-M' we have a lot of models ill-behaving as a result of bad numerics. If textbooks just would call this 'small-m' instead, new modelers would not have the urge to use these ridiculously large numbers.
It is my conjecture that just because of this name 'big-M' we have a lot of models ill-behaving as a result of bad numerics. If textbooks just would call this 'small-m' instead, new modelers would not have the urge to use these ridiculously large numbers.
As indicated in the comments some solvers support indicator constraints (Cplex, as well as Scip and Xpress I believe) which allow you to formulate implications without big-M constants.
Sunday, May 22, 2016
Mixed Integer Programming Class Library (MIPCL)
Open source MIP solver (GNU LGPL license)
Just came across this solver. Looks interesting,- Short paper: https://www.researchgate.net/publication/301802844_Mixed_Integer_Programming_Class_Library_MIPCL
- Web site: http://www.mipcl-cpp.appspot.com/
- Documentation: http://www.mipcl-cpp.appspot.com/static/docs/mipcl/html/index.html
- Could not find any sources, so may be not open source after all. I was under the impression that LGPL implied making the source code available for modification.
Thursday, May 19, 2016
Strange scheduling problem
A somewhat strange scheduling model was presented to me. We operate a machine in one of several operating modes \(i\). We have time periods \(t\) and the operating cost \(c^{op}_{i,t}\) changes over time. Obviously, the best schedule would be to pick in each period \(t\) the cheapest operating mode. Now we add a changeover cost: when we change from mode \(i\) to mode \(j\) we pay some cost \(c^{ch}_{i,j}\). This would require a real optimization model to find the optimal operating sequence.
Here is a simple model to handle this:
We have used some random data here, which looks like:
| ---- 12 PARAMETER opcost operating cost period1 period2 period3 period4 period5 period6 period7 mode1 0.172 0.843 0.550 0.301 0.292 0.224 0.350
mode1 mode2 mode3 mode4 mode5 mode1 0.572 1.188 1.445 1.256 |
The solution looks like:
| ---- 34 VARIABLE x.L operating schedule period1 period2 period3 period4 period5 period6 period7 mode1 1 1 1
period4 mode3.mode1 1
|
We see one changeover between periods 4 and 5. Notice that \(y\) is defined in the model above as:
\[y_{i,j,t} = \begin{cases} 1 \> \text{if $x_{i,t}=1$ and $x_{i,t+1}=1$} \\
0 \> \text{otherwise}\end{cases}\]
This can also be interpreted as a nonlinear constraint \(y_{i,j,t} = x_{i,t}x_{i,t+1}\). In the model we have linearized this. A refinement of the model would have us look at the operating mode just before we start scheduling (i.e. the operating mode in period 0). To handle that easily it is helpful to change the definition of \(y\) to:
\[y_{i,j,t} = \begin{cases} 1 \> \text{if $x_{i,t-1}=1$ and $x_{i,t}=1$} \\
0 \> \text{otherwise}\end{cases}\]
Only equation order needs to change:
Note that \(x0\) is an extremely sparse matrix: it has only one element corresponding to the operating mode for the machine just before period 1. (We used here a GAMS feature: addressing lags outside the domain cause the symbol to be dropped; so \(x_{i,t-1}\) for \(t=period1\) will disappear; in that special case the parameter \(x0\) will kick in). The final results with this initial changeover cost is:
| ---- 37 VARIABLE x.L operating schedule period1 period2 period3 period4 period5 period6 period7 mode1 1 1 1
period4 period5 mode3.mode1 1
|
Note that the initial mode in period 0 was 4, and now we keep operating using that mode for a little while.
Subscribe to:
Posts (Atom)