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:
Look if the facility was built in any period \(t’<t\):
| \[open_t = \sum_{t’|t’<t} build_{t’}\] |
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}\] |
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) |
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 **** SOLVER STATUS 1 Normal Completion RESOURCE USAGE, LIMIT 0.031 1000.000 |
This approach can also be applied to the Sudoku MIP model.
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.
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.
 | 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:
|
For a project I was looking at a simple two stage stochastic integer programming problem (1).
I find it always very rewarding to reproduce the results in the paper. Re-implementing the model forces me to not skip the details. Compared to superficial reading it requires to pay attention to all aspects. Of course in many cases journal papers do not allow to reproduce results because crucial details or data are left out. Presenting a stylized model illustrating and illuminating the concepts is only done sporadically, may be because it is actually hard to create such a “simple” model, and it may look less “scientific”.
The model is a fairly standard 2 stage SLP. In stage 1 we can implement so-called permanent floodplain management actions (raising structures, implement warning systems, etc.) The realizations are floods of different sizes. In stage 2 we can apply certain emergency measures (evacuations, issue sandbags etc.).
After re-implementing the model I saw the same objective function values, but the stage 2 variables differed a bit.
The paper shows a table with the stage 2 results:
In GAMS we do not see any of the levee monitoring being implemented:
| ---- 145 VARIABLE xe.L stage 2 decisions: emergency options peak5-6 peak6-8 peak8-10 peak10+ evacuate 1.000 1.000 1.000 |
GAMS will report marginals (duals and reduced costs) even for a MIP. They fix the integer variables and resolve as an LP to get these. Here are the relevant marginals:
| ---- 145 VARIABLE xe.M stage 2 decisions: emergency options peak0-5 peak5-6 peak6-8 peak8-10 peak10+ evacuate 80000.000 7900.000 -12000.000 -8000.000 -9000.000 |
This EPS means this variable is non-basic with a zero reduced cost. This can indicate the presence of multiple optima. Indeed when using a different MIP solver I see:
| ---- 145 VARIABLE xe.L stage 2 decisions: emergency options peak5-6 peak6-8 peak8-10 peak10+ evacuate 1.000 1.000 1.000 |
which is the same as the results in the paper.
Digging a little bit further, we can see that the damage reduction by implementing emergency measure “monitor” is the same as its cost. This means we are indifferent whether to apply this option in this scenario.
| ---- 154 PARAMETER red2 damage reduction b(j,s)*xe.l(j,s) peak5-6 peak6-8 peak8-10 peak10+ evacuate 300000.000 500000.000 1000000.000
peak5-6 peak6-8 peak8-10 peak10+ evacuate 100000.000 100000.000 100000.000 |
I would guess the chance this will happen with real data is smaller than with these nicely rounded invented data.
The model itself is implemented as a fairly standard deterministic equivalent formulation:
Project homepage: http://pymprog.sourceforge.net/
Uses GLPK as solver. Here is a small example:
from pymprog import * c = (10, 6, 4) A = [ ( 1, 1, 1), ( 9, 4, 5), ( 2, 2, 6) ] b = (10, 60, 30) begin('basic') # begin modelling verbose(True) # be verbose x = var('x', 3) #create 3 variables maximize(sum(c[i]*x[i] for i in range(3))) for i in range(3): sum(A[i][j]*x[j] for j in range(3)) <= b[i] solve() # solve the model sensitivity() # sensitivity report end() #Good habit: do away with the model |
There is quite some competition in this space. Some of the major competitors:
This does not happen a lot, but in some cases we want to sort decision variables. This turns out not such an easy exercise.
If we deal with just the \(\min()\) or \(\max()\) function we can do things better than sorting. An important trick is to exploit any so-called convexity so we don’t need to add extra binary variables. Even if the problem is not convex, there is no need to sort things.
Interestingly, minimizing the sum of the K largest values also does not need sorting. See (4).
At first sight this is not very useful exercise: we can sort a parameter – constants in the model – outside of the model. However the problem demonstrates some of the concepts I will use later on. Assume \(a_j\) is our given parameter. A permutation \(y\) of \(a\) can be written as \(y=Pa\) where \(P\) is a permutation matrix (1). A permutation matrix is an identity matrix \(I\) with rows and columns interchanged. Such a matrix has exactly a 1 in each row and in each column. A different way to look at this is as an assignment problem with binary variables \(p_{i,j} \in \{0,1\}\).
The linear constraints to sort \(a_j\) (descending) can look like:
| \[\boxed{\begin{align} &\sum_i p_{i,j} = 1 \>\> \forall j\\ &\sum_j p_{i,j} = 1 \>\> \forall i\\ &y_i = \sum_j p_{i,j} a_j \\ &y_i \ge y_{i+1} \\ &p_{i,j} \in \{0,1\} \end{align}}\] |
Here indices \(i\) and \(j\) are from sets with the same cardinality (\(n\)), i.e. \(p_{i,j}\) is a square matrix. Note that we use \(n^2\) binary variables here, so don’t expect this to work for very large vectors \(a_j\) (2). The same structure can also be used in the context of all-different constraints (3).
| ---- 32 PARAMETER a j1 2, j2 4, j3 1, j4 3
j1 j2 j3 j4 i1 1
i1 4, i2 3, i3 2, i4 1 |
Instead of sorting a parameter \(a_j\), consider the problem of sorting a variable \(x_j\). We cannot just replace \(a_j\) by \(x_j\) as
| \[y_i = \sum_j p_{i,j} x_j \] |
is non-linear. There is a way to linearize the product \(q_{i,j} = p_{i,j} x_j\) as this is a multiplication of a binary variable with a continuous variable. Assuming \(x_j \in [0,U_j]\), we have:
| \[\begin{align} &0\le q_{i,j} \le U_j p_{i,j}\\ &x_j – U_j (1-p_{i,j}) \le q_{i,j} \le x_j \end{align}\] |
| ---- 41 VARIABLE x.L j1 0.172, j2 0.843, j3 0.550, j4 0.301
j1 j2 j3 j4 i1 1
j1 j2 j3 j4 i1 0.843
i1 0.843, i2 0.550, i3 0.301, i4 0.172 |
A different way to formulate a sorting scheme is to look at the problem as a special transportation problem. Transportation equations look like:
| \[\boxed{\begin{align} &\sum_i z_{i,j} = x_j \>\> \forall j\\ &\sum_j z_{i,j} = y_i \>\> \forall i \end{align}}\] |
As can be seen this is linear in \(x\). We still need to enforce that one supply node is linked to exactly one demand node. This can be done using an additional assignment block:
| \[\boxed{\begin{align} &\sum_i p_{i,j} = 1 \>\> \forall j\\ &\sum_j p_{i,j} = 1 \>\> \forall i\\ &\sum_i z_{i,j} = x_j \>\> \forall j\\ &\sum_j z_{i,j} = y_i \>\> \forall i\\ &y_i \ge y_{i+1} \\ &0\le z_{i,j} \le p_{i,j} U_j\\ &p_{i,j} \in \{0,1\} \end{align}}\] |
| ---- 43 VARIABLE x.L j1 0.172, j2 0.843, j3 0.550, j4 0.301
j1 j2 j3 j4 i1 1
j1 j2 j3 j4 i1 0.843
i1 0.843, i2 0.550, i3 0.301, i4 0.172 |
From this post:
I'm currently stuck with a MIP program where the interest rate, i, is based on the number of units produced for Housing Plan A. If the number of plan A houses sold is the highest among all four types then i=1. If the number of plan A houses sold is the second highest, then i=2 and so on up to i=4. The interest rate is basically 2i%. Not really sure how to add constraints that will represent the position of plan A houses and implement the correct interest rate in the objective function. The objective function maximizes the total profit (e.g 50,000A + 40,000B + 70,000C + 80,000D). Any ideas on how to use binary variables to represent position?
Some results can look like:
| ---- 47 VARIABLE x.L A 100.000, B 90.000, C 120.000, D 150.000
1 2 3 4 A 1.000
1 2 3 4 A 100.000
1 150.000, 2 120.000, 3 100.000, 4 90.000
A 0.060, B 0.080, C 0.040, D 0.020 |
There are a number of very fast CSV file readers available in R and Python. Lets have a quick test to see how they compare.
I generated a very simple, but large CSV file with 100 million records using a GAMS script as follows:
| set |
The generated CSV file looks like:
| D:\tmp\csv>head d.csv D:\tmp\csv>dir d.* Directory of D:\tmp\csv 12/08/2016 03:42 PM 3,656,869,678 d.csv D:\tmp\csv> |
We also see the CSV file is much larger than the intermediate (compressed) GAMS GDX file.
This is the default CSV reader in R.
> system.time(d<-read.csv("d.csv")) user system elapsed 1361.61 50.56 1434.39 |
read_csv is from the readr package, and it is much faster for large CSV files:
|
Would it help to read a compressed CSV file?
> system.time(d<-read_csv("d2.csv.gz")) Error in .Call("readr_read_connection_", PACKAGE = "readr", con, chunk_size) : negative length vectors are not allowed Timing stopped at: 57.53 4.43 62.29 |
Bummer. I have no idea what went wrong here. May be we hit some size limit (note the CSV file is larger than 2 gb; other compression formats gave the same result).
As mentioned in the comments, the package data.table has a function fread.
> system.time(dt<-fread("d.csv")) Read 100000000 rows and 5 (of 5) columns from 3.406 GB file in 00:04:40 user system elapsed 275.19 4.33 281.05 |
Quite fast:
t0=pc() df=pd.read_csv("d.csv") print(pc()-t0) |
158.2270488541103 |
The paratext library should be even faster.
Fuel is cheap. Why are we still paying to check bags?
Interesting video: [Video] A Huge Debate: R vs. Python for Data Science (presentation by Eduardo Ariño de la Rubia)
Packages discussed:
The new version of Cplex (12.7) has some interesting new facilities to make it easier to use a Benders Decomposition approach (1,3). Paul Rubin shares some experiences with this in (2).
The original paper by Jacques Benders can be found in (4). These old papers are always an interesting read. The style and notation is often very different from how papers look like these days, although in this case the notation is quite modern.
From this post we have a problem stated as:
My wife teaches AM and PM kindergarten classes. AM has 14 students and PM 11. At the beginning of each month, she puts out a new seating chart where she rotates students in such a way that they (ideally) sit at a different table and with different students for that month.
There are 3 students per table, but if numbers force the issue, the last may have more or less. We realize that, by the end of the year, there will be some unavoidable situations where students end up sitting with each other or at the same tables again, and this is okay.
Something like this I imagine:
A somewhat related problem is discussed in (1). As usual we have a set of binary variables that indicate the assignment of students to tables:
| \[x_{s,t,m} = \begin{cases}1&\text{if student $s$ is seated at table $t$ in month $m$}\\ 0&\text{otherwise}\end{cases}\] |
The first set of equations are somewhat straightforward. Each student sits at exactly one table:
| \[\sum_t x_{s,t,m} = 1 \>\>\forall s,m\] |
We cannot exceed the capacity of a table:
| \[\sum_s x_{s,t,m} \le cap_t \>\> \forall t,m\] |
In the complete model below I used an equality constraint for this: there is no slack capacity. Secondly we want to keep track of how many times students sit at the same table. We introduce binary variables \(meet_{s1,s2,t,m}\in\{0,1\}\) indicating if two students \(s1\) and \(s2\) sit at the same table. The obvious nonlinear equation would be:
| \[meet_{s1,s2,t,m} = x_{s1,t,m} \cdot x_{s2,t,m}\>\> \forall s1<s2,t,m\] |
Note that we can exploit symmetry here: if we already compared \(s1,s2\) we do not need to bother about \(s2,s1\). The standard linearization is
| \[\begin{align} & meet_{s1,s2,t,m} \le x_{s1,t,m}\\ & meet_{s1,s2,t,m} \le x_{s2,t,m}\\ & meet_{s1,s2,t,m} \ge x_{s1,t,m}+x_{s1,t,m}-1\end{align}\] |
As we push variable \(meet\) down, we can drop the first two inequalities and just keep the last. Note that variable \(meet\) is now a bound instead of the product \(x_{s1,t,m} \cdot x_{s2,t,m}\). When reporting we should use the underlying variables \(x\).
A count of the number of times students \(s1,s2\) sit at the same table is a simple aggregation:
| \[meetcount_{s1,s2} = \sum_{t,m} meet_{s1,s2,t,m}\] |
An obvious objective is to minimize the expression \(\max meetcount_{s1,s2}\). This is easily linearized:
| \[\begin{align}\min\>&maxmeetcount\\ &maxmeetcount \ge meetcount_{s1,s2} \>\>\forall s1<s2\end{align}\] |
To test this we use 14 students, 11 months (I assume one month vacation period). The capacities of the 5 tables is 3 except for one table with 2 students (this is to make the capacities sum up to 14).
The equations are exactly the same as discussed above:
To help the solver we fix the table assignment for the first month. In addition I assume student 1 is always sitting at table 1. This makes the model easier to solve. Note that I only fix some variables to 1. As a result a lot of other variables should be zero. E.g. if \(x_{s1,t1,m}=1\) we know that \(x_{s1,t2,m}=0\) and similar for other tables than \(t1\). Instead of fixing all these these myself to zero, I leave it to the MIP presolver to do that for me. This fixing step helps a lot.
We achieve a quick solution with an objective of 2. I.e. two students will sit at most twice at the same table.
In the above model we minimized the maximum number of times students sit at the same table, The optimal value is 2. However this means that we allow the number of times two students meet (the meet count, “meets” in the tables below) to float between 0 and 2 without a further incentive to decrease the average meet count within this band width. This leads to a distribution as follows:
On average two students sit at the same table with another student 1.57 times. Detail: note that we need to base the calculations on the optimal values of \(x\) instead of \(meetcount\) as \(meetcount\) is using variable \(meet\) which is just a bound on the product \(x_{s1,t,m} \cdot x_{s2,t,m}\).
We can try to bring the average down by minimizing the sum of the meet counts. The results show, this will not actually bring the average down. But the distribution is certainly different. This objective will produce a few student pairs that have a high meet count:
We can try to minimize the number of student pairs meeting twice or more while still pushing down the maximum by using a composite objective:
| \[\min maxmeetcount + \frac{count2}{100}\] |
where \(count2\) is the number of pairs meeting twice or more. This gives a slightly different distribution:
A depiction of this distribution:
Another interesting objectives:
When we make a picture of when student pairs sit at the same table we get something like this:
We color coded the cells whether we have 0, 1, 2 or 3 or more months in between (the cells with a ‘3’ indicate that meets are separated by 3 or more months). Note we did not display the pairs that only sit at the same table once.
To forbid the cases 0 and 1 in the above solution we can add the constraints:
| \[\begin{align} &meetm_{i1,i2,m} = \sum_t meet_{i1,i2,t,m}\\ &meetm_{i1,i2,m}+meetm_{i1,i2,m+1}+meetm_{i1,i2,m+2} \le 1\end{align} \] |
Here \(meetm\) is a simple aggregation of the variable \(meet\). The second constraint only allows one time the pair \((i1,i2)\) sit at the same table in each three month period.
Note that now we have the variable \(meetm\) we can optimize constraint that calculates the meet count:
| \[meetcount_{s1,s2} = \sum_m meetm_{s1,s2,m}\] |
This will reduce the number of nonzero elements in this equation, so this is a good idea.
The new picture
shows only 2s and 3s (where again 3 means at least 3 months in between). I.e. we successfully removed all cases where student pairs sit at the same table with zero or one month in between.
