Tuesday, July 18, 2017

Ubuntu Bash on Windows

Windows 10 has a beta feature: Bash.exe with an Ubuntu Linux subsystem (see (1) and (2)), also known as WSL (Windows Subsystem for Linux). This will allow you to run Linux command line tools from within Windows. This is not the same as running Linux in a Virtual Machine (VM) using virtualization software like VirtualBox or VMware. It is also different from Cygwin or MINGW which is a windows port of GNU tools. In the case of  WSL, real, unmodified Ubuntu Linux binaries are executed.

learnbash

LXSS diagram

There is a single Linux instance, so when invoking bash.exe several times, you talk to the same instance:


wsl


The main issue I noticed is that some networking commands (like ping) requires bash.exe to be started as administrator.

References

  1. Jack Hammons, Bash on Ubuntu on Windows, https://msdn.microsoft.com/commandline/wsl
  2. Learning about Bash on Windows Subsystem for Linux, https://blogs.msdn.microsoft.com/commandline/learn-about-bash-on-windows-subsystem-for-linux/
  3. VirtualBox, https://www.virtualbox.org/
  4. VMware, https://www.vmware.com/
  5. Cygwin, https://www.cygwin.com/
  6. MINGW, Minimalist GNU for Windows, http://www.mingw.org/

Monday, July 17, 2017

Modeling flood damages

While working on a investment planning model to combat damages resulting from flooding, I received the results from a rainfall model that calculates damages as a result of excessive rain and flooding. The picture the engineers produced looks like:

image 

This picture has the scenario on the x-axis (sorted by damage) and the damages on the y-axis. This picture is very much like a load-duration curve in power generation.

For a more “statistical” picture we can use standard histogram (after binning the data):

image

Gamma distribution

We can use standard techniques to fit a distribution. When considering a Gamma distribution (1), a simple approach is the method of moments. The mean and the variance of the Gamma distribution with parameters \(k\) and \(\theta\) are given by:

\[\begin{align} &E[X] = k \cdot \theta\\ & Var[X] = k \cdot \theta^2 \end{align} \]

Using sample mean \(\mu\) and standard deviation \(\sigma\), we can solve:

\[\begin{align} & k \cdot \theta = \mu\\ & \sqrt{k} \cdot \theta = \sigma \end{align} \]

This can be easily solved numerically and it actually seems to work:

image

Weibull distribution

An alternative distribution that is sometimes suggested is the Weibull distribution (2).  The method of moments estimator for the Weibull distribution with parameters \(\lambda\) and \(k\) can be found by solving the system:

\[\begin{align} & \lambda \Gamma(1+1/k) = \mu\\ &\lambda \sqrt{\Gamma(1+2/k)+\left(\Gamma(1+1/k)\right)^2} = \sigma \end{align}\]
I was unable to get a solution from this: solvers failed miserably on this.

An alternative approach would be to use an MLE (Maximum Likelihood Estimation) technique. This yields a system of equations:

\[\begin{align} & \lambda^k = \frac{1}{n}\sum_{i=1}^n x_i^k\\ &k^{-1} = \frac{\sum_{i=1}^n x_i^k \ln x_i}{\sum_{i=1}^n x_i^k} – \frac{1}{n}\sum_{i=1}^n \ln x_i \end{align}\]
Note that we can solve the second equation first to solve for \(k\) and then calculate \(\lambda\) using the first equation. A solver like CONOPT will do this automatically by recognizing the triangular structure of this system. Note that our data set contains many \(x_i=0\). These are replaced by \(x_i=0.001\) so we can take the logarithm.

This gives is a very similar picture:

image

References
  1. Gamma distribution, https://en.wikipedia.org/wiki/Gamma_distribution
  2. Weibull distribution, https://en.wikipedia.org/wiki/Weibull_distribution

Tuesday, July 11, 2017

Rectangles: no-overlap constraints

The question of how we can formulate constraints that enforce rectangles not to overlap comes up regularly (1). There are basically four cases to consider when considering two rectangles \(i\) and \(j\):

image

In the picture above, we indicate by \((x_i,y_i)\) the left-lower corner of a rectangle (these are typically decision variables). The height and the width are \(h_i\) and \(w_i\) (these are typically constants). From the above, we can formulate the “no-overlap”  constraints as:

\[
\bbox[lightcyan,10px,border:3px solid darkblue] {
\begin{align}
&x_i+w_i \le x_j  & \text{or} \\
&x_j+w_j \le x_i  & \text{or} \\
&y_i+h_i \le y_j  & \text{or} \\
&y_j+h_j \le x_i
\end{align}
}\]

The “or” condition can be modeled with the help of binary variables \(\delta\) and big-M constraints:

\[
\bbox[lightcyan,10px,border:3px solid darkblue] {
\begin{align}
&x_i+w_i \le x_j + M_1 \delta_{i,j,1}\\
&x_j+w_j \le x_i + M_2 \delta_{i,j,2}\\
&y_i+h_i \le y_j + M_3 \delta_{i,j,3}\\
&y_j+h_j \le x_i+ M_4 \delta_{i,j,4}\\
&\sum_k \delta_{i,j,k} \le 3\\
& \delta_{i,j,k} \in \{0,1\}
\end{align}
}\]

The binary variables make sure at least one of the constraints is active (not relaxed).

If we have multiple rectangles we need to compare all combinations, i.e. we generate these constraints for  \(\forall i,j\). Except: we can not compare a rectangle with itself. A first approach would be to generate the above constraints for all \(i\ne j\). This is correct, but we can apply an optimization. Note that we actually compare things twice: if rectangles \(i\) and \(j\) are non-overlapping then we don’t need to check the pair of rectangles \(j\) and \(i\). That means we only need constraints and variables \(\delta_{i,j,k}\) for \(i<j\).

It is important to make constants \(M_k\) as small as possible. In general we have a good idea about good values for these \(M_{k}\)’s. E.g. consider the case where we need to pack boxes in a container (2):

image

If the container has length and width \(H\) and \(W\) then we can easily see: \(M_1=M_2=W\) and \(M_3=M_4=H\).

References
  1. How to write this constraint in LPSolve?, https://stackoverflow.com/questions/44858353/how-to-write-this-constraint-in-lp-solver-logical-and-and-doesnot-exist
  2. Filling up a container with boxes, http://yetanothermathprogrammingconsultant.blogspot.com/2016/06/filling-up-container-with-boxes.html

Monday, July 10, 2017

GLPK was used in the proof of the 300 year old Kepler conjecture

image

A lot of authors.

Apparently 43,078 LPs had to be solved.

References
  1. Thomas Hales e.a., A Formal Proof of the Kepler Conjecture, Forum of mathematics, Pi (2017), vol 5, e2, https://www.cambridge.org/core/services/aop-cambridge-core/content/view/78FBD5E1A3D1BCCB8E0D5B0C463C9FBC/S2050508617000014a.pdf/formal_proof_of_the_kepler_conjecture.pdf
  2. Original message appeared in the Help-GLPK mailing list,  https://lists.gnu.org/archive/html/help-glpk/2017-07/msg00001.html

Tuesday, June 27, 2017

Who needs Cplex or Gurobi: solving LPs using LU

The strange notion that linear programming problems can be solved simply by applying an LU decomposition (not as part of the Simplex method but just LU) seems to take hold. The original idea is from (1) but now we also can use LU as a solver for “Fuzzy Linear Fractional Programming Problems” (2).


image 

I agree with the abstract: the proposed approach is simple.

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 method does not work except in some very special cases (e.g. all constraints binding in the optimal solution).
  2. S. Muruganandam and P. Ambika,  Solving Fuzzy Linear Fractional Programming Problem using LU Decomposition Method, Annals of Pure and Applied Mathematics, Vol. 13, No 1., 2017, pp 89-97.
  3. Solving LP by LU??, http://yetanothermathprogrammingconsultant.blogspot.com/2017/01/solving-lp-by-lu.html

Wednesday, June 21, 2017

Minimizing the k-th largest x

Minimizing the largest \(x_i\) is an easy exercise in LP modeling:
\[\bbox[lightcyan,10px,border:3px solid darkblue] { \begin{align} \min\>& z\\ & z \ge x_i \end{align}} \]
This is sometimes called MINIMAX.

What about minimizing the \(k^{\text{th}}\) largest value? Here is one possible MIP formulation:

\[\bbox[lightcyan,10px,border:3px solid darkblue] { \begin{align} \min\>& z\\ & z \ge x_i - \delta_i M\\ & \sum_i \delta_i = k-1\\ & \delta_i \in \{0,1\} \end{align}} \]

I.e., we have \(k-1\) exceptions on the constraint \(z \ge x_i\). The objective will make sure the exceptions are the largest \(x_i\).  We should think long and hard about making the big-M’s as small as possible. If you have no clue about proper bounds on \(x_i\) one could use indicator constraints.

Interestingly, minimizing the sum of the \(k\) largest can be modeled as a pure LP (1) :

\[\bbox[lightcyan,10px,border:3px solid darkblue] { \begin{align} \min\>& \sum_i v_i + k\cdot q\\ & v_i \ge x_i - q\\ &v_i \ge 0 \\ &q \text{ free } \end{align}} \]

You almost would think there is an LP formulation for minimizing the \(k^{\text{th}}\) largest value. I don't see it however.
References
  1. Comment by Michael Grant in http://orinanobworld.blogspot.se/2015/08/optimizingpartoftheobjectivefunction-ii.html.

Wednesday, June 14, 2017

Modeling production runs with length exactly three

In (1) the question was posed how to model production runs that have an exact length. Although the question was asked in terms of integer variables, this is much easier to deal with when we have binary variables. Let’s introduce a binary variable:

\[x_{t} = \begin{cases}1 & \text{when the unit $x$ is turned on in period $t$}\\0 &\text{otherwise}\end{cases}\]

The question was formulated as an implication:

\[x_t=1 \implies x_{t+1}=x_{t+2}=1\]

This implication is not 100% correct: it would turn on \(x\) forever. But we understand what the poster meant, if a unit is turned on, it should stay on for exactly three periods.

image

This condition can be more correctly stated as the implication:

\[x_{t-1}=0 \text{ and } x_{t}=1 \implies x_{t+1}=x_{t+2}=1 \text{ and } x_{t+3}=0\]

We can make linear inequalities out of this as follows:

\[\bbox[lightcyan,10px,border:3px solid darkblue]{\begin{align}&x_{t+1}\ge x_{t}-x_{t-1}\\&x_{t+2}\ge x_{t}-x_{t-1}\\&1-x_{t+3}\ge x_{t}-x_{t-1}\end{align}}\]

If we want to model the condition “a unit should stay on for at least three periods”, we can drop the last constraint and just keep:

\[\begin{align}&x_{t+1}\ge x_{t}-x_{t-1}\\&x_{t+2}\ge x_{t}-x_{t-1}\end{align}\]


References
  1. https://stackoverflow.com/questions/44496473/block-of-consecutive-variables-to-have-same-value-in-mixed-integer-linear-progra

Friday, June 9, 2017

A Staffing Problem

In (1) a "simple problem" is stated (problems are rarely as simple as they seem):

For the next 18 weeks there is some data about demand for staffing resources:

image

Actually we don’t need the data for the individual projects: just the totals. It is noted that the maximum number of “bodies” we need is 52 (in week 3).

We start with 48 temps available in the first period. We can let a temp staffer go and we can hire new ones. However when hiring a staffer it will cost 10 days (2 weeks) to train this person. During training a staffer is not productive.  We can also keep a staffer idling for a few periods.

To make a difference between idling and training, as shown here:

image

I will assume that hiring + training is slightly more expensive than keeping someone idle. I think it makes sense to assign some cost to the hiring and training process. This property was not in the original post but I think I have convinced myself that this is actually a reasonable assumption.

To model this I introduce two sets of binary variables:

\[\begin{align}
& r_{i,t} = \begin{cases} 1 & \text{if a staffer $i$ is available for training or work during period $t$}\\
0 & \text{otherwise} \end{cases} \\
& h_{i,t} = \begin{cases} 1 & \text{if a staffer $i$ is hired at the beginning of period $t$}\\
0 & \text{otherwise} \end{cases}
\end{align}\]

We can link \(h_{i,t}\) and \(r_{i,t}\) as follows:

\[h_{i,t}\ge r_{i,t}-r_{i,t-1}\]

This implements the implication:

\[r_{i,t-1}=0 \text{ and } r_{i,t}=1 \implies h_{i,t}=1\]

that is: if we change \(r\) from 0 to 1, we have a new hire. We will add a (small) cost to a hire to the objective so we don’t need to add constraint to enforce the other way around:

\[ r_{i,t-1}=1 \text{ or } r_{i,t}=0 \implies h_{i,t}=0\]

There is one wrinkle here: we need to make sure that the 48 staffers already hired before can work in the first period without being considered to be hired. We can explicitly model this:

\[h_{i,t}\ge \begin{cases} r_{i,t}-r_{i,t-1} &\text{if $t>1$}\\
r_{i,t}&\text{if $t=1$ and $i>48$}\end{cases}\]

Actually in the GAMS model below I approached this slightly differently, but with the same net result.

The main equation in the model is to make sure we have enough staffing in each period \(t\). This can be modeled as:

\[\sum_i \left( r_{i,t} – h_{i,t} - h_{i,t-1}\right) \ge \mathit{demand}_t \]

The complete model can look like:

image

Notes:

  • \(\mathit{rinit}_{i,t}\) is used to indicate initial staffing when we start period 1. It is sparse: it assumes the value 1 only if \(i=1\) and \(t \le 48\).
  • This allows us to write equation EHire as one single, clean equation.
  • Note that GAMS will assume a variable is zero outside its domain:  \(r_{i,t-1}\) is zero when \(t=1\).
  • I used the maximum demand (52 in week 3) to dimension the the problem: set \(i\) has 52 members.
  • The model does not decide which workers are idle. We can see in each period how many workers we have and how many are in training. The remaining ones are working or sitting idle.
  • We reuse resource numbers \(i\). In the results below worker \(i=21\) is likely to be two different persons.

The results look like:

image

 image

Here a green cell with code=2 means hired and working (or idle, but not in training). A red cell with code=1 means in training.  The row “total” is the number of staffers (either working, idling or training). The row “Work” indicates the number of workers available (not in training). We are overshooting demand a bit: some workers are idling.

A picture is always a good idea, so here we see demand vs. staffing available for work.

image

We see two place where we increase the workforce: we hire in week 1 to handle demand in week 3, and we hire again in week 7 to deal with demand in week 9.

References

Wednesday, June 7, 2017

Minimum down- and up-time

In machine scheduling models we sometimes want to impose minimum up-time and minimum down-time restrictions. E.g., from (1):

My question is how do I add in additional constraints that if the factory switches off then its needs to stay off for 3 months, and if it switches back on then it needs to stay on for 4 months?

One possible solution is the following. Let us define our binary decision variable by

\[x_{i,t} = \begin{cases}1 & \text{if factory $i$ is operating in month $t$} \\
  0&\text{otherwise} \end{cases}\]
Method 1

We really want to forbid short run patterns such as 010 (i.e. off-on-off), 0110, 01110 and 101, 1001. Forbidding patterns 010, 0110, 01110 will ensure a factory is up at least 4 consecutive periods. By not allowing 101,1001 we really make sure that a down time period is at least three months. We can model these restrictions in a linear fashion as follows (2):

Forbid 010 \(-x_{i,t}+x_{i,t+1}-x_{i,t+2}\le 0\)
Forbid 0110 \(-x_{i,t}+x_{i,t+1}+x_{i,t+2}-x_{i,t+3}\le 1\)
Forbid 01110 \(-x_{i,t}+x_{i,t+1}+x_{i,t+2}+x_{i,t+3}-x_{i,t+4}\le 2\)
Forbid 101 \(x_{i,t}-x_{i,t+1}+x_{i,t+2}\le 1\)
Forbid 1001 \(x_{i,t}-x_{i,t+1}-x_{i,t+2}+x_{i,t+3}\le 1\)
Method 2

A different approach is as follows. First define binary variables:

\[ \begin{align} &\delta^{\mathit{on}}_{i,t} = \begin{cases}1 & \text{if factory $i$ is turned on in month $t$} \\
  0&\text{otherwise} \end{cases}\\ &\delta^{\mathit{off}}_{i,t} = \begin{cases}1 & \text{if factory $i$ is turned off in month $t$} \\
  0&\text{otherwise} \end{cases} \end{align} \]

Mathematically we write this as:

\[ \begin{align} &\delta^{\mathit{on}}_{i,t} = (1-x_{i,t-1})\cdot x_{i,t}\\ &\delta^{\mathit{off}}_{i,t} = x_{i,t-1}\cdot (1-x_{i,t})\\ \end{align} \]

We can linearize these non-linear equations by:

\[ \begin{align} &\delta^{\mathit{on}}_{i,t} \le 1-x_{i,t-1}\\ &\delta^{\mathit{on}}_{i,t} \le x_{i,t}\\ &\delta^{\mathit{on}}_{i,t} \ge x_{i,t}-x_{i,t-1}\\ &\delta^{\mathit{off}}_{i,t} \le x_{i,t-1}\\ &\delta^{\mathit{off}}_{i,t} \le 1-x_{i,t}\\ &\delta^{\mathit{off}}_{i,t} \ge x_{i,t-1}-x_{i,t}\\ \end{align} \]

With these variables we can implement the implications:

\[ \begin{align} &\delta^{\mathit{on}}_{i,t} = 1 \implies x_{i,t} + x_{i,t+1} + x_{i,t+2} + x_{i,t+3} = 4 \\ &\delta^{\mathit{off}}_{i,t} = 1 \implies x_{i,t} + x_{i,t+1} + x_{i,t+2} = 0 \end{align} \]

This can be linearized as:

\[ \begin{align} & x_{i,t+1} + x_{i,t+2} + x_{i,t+3} \ge 3 \delta^{\mathit{on}}_{i,t} \\ & x_{i,t+1} + x_{i,t+2} \le 2 (1-\delta^{\mathit{off}}_{i,t}) \end{align} \]

Note that I dropped \(x_{i,t}\) in both inequalities. These are already known from the definition of \(\delta^{\mathit{on}}_{i,t}\) and \(\delta^{\mathit{off}}_{i,t}\).

In the comments it is mentioned by Rob Pratt that we can strengthen this a bit (the math does not look very good in a comment, so I repeat it here):

\[\begin{align} & x_{i,t+1} \ge \delta^{\mathit{on}}_{i,t}\\
& x_{i,t+2} \ge \delta^{\mathit{on}}_{i,t}\\
& x_{i,t+3} \ge \delta^{\mathit{on}}_{i,t}\\
& x_{i,t+1} \le 1-\delta^{\mathit{off}}_{i,t} \\
& x_{i,t+2} \le 1-\delta^{\mathit{off}}_{i,t} \end{align}\]

References
  1. How do I add a constraint to keep a factory switched on or off for a certain period of time in PuLP? https://stackoverflow.com/questions/44281389/how-do-i-add-a-constraint-to-keep-a-factory-switched-on-or-off-for-a-certain-per/44293592
  2. Integer cuts, http://yetanothermathprogrammingconsultant.blogspot.com/2011/10/integer-cuts.html

Friday, May 26, 2017

Working on an R package….

The C++ code (parsing Excel formulas so we can for instance execute them) is not working 100% correctly as of now…

Capture

Time to make it work under a debugger (amazingly I did not need a debugger until now).

Update: found it (without debugger). Dereferencing a nil-pointer. Is it unreasonable to expect a (descriptive) exception to be raised when this happens?