Sunday, July 23, 2017

The CDF of the Gamma distribution in GAMS

For a model (1) I needed to calculate the cumulative distribution function \(F(x)\) of the Gamma distribution. There is no uniform consensus:  Gamma distribution has different parametrizations. I use the \(k, \theta\) parametrization from (2). This would yield (3):

\[F(x) =  \gamma( \frac{x}{\theta},k)\]

where \(\gamma(x,a)\) is the incomplete Gamma function defined by:

\[\gamma(x,a) = \frac{1}{\Gamma(a)} \int_0^x t^{a-1}e^{-t} dt\]

The function \(\Gamma(x)\) is the (complete) Gamma function:

\[\Gamma(x) = \int_0^{\infty} t^{x-1} e^{-t} dt \]

Unfortunately there a few alternative definitions of the incomplete Gamma function. In (2) the following definition is used:

\[\gamma(s,x) = \int_0^x t^{s-1} e^{-t} dt\]

So watch out: different definitions are being used.

As  a result we can implement the CDF of the Gamma distribution with parameter \(k\) and \(\theta\) as follows (2):


There is also a more obscure way to calculate this using an extrinsic function. The syntax is unfortunately needlessly complicated and the documentation is really, really bad (4).

/ 0.363 /
/ 27.863 /
/ 1.5 /

p = gammareg(x/theta,k);
display p;

$FuncLibIn stodclib stodclib
function cdfGamma / stodclib.cdfGamma /;
p = cdfGamma(x,k,theta);


  1. Modeling Flood Damages,
  2. Gamma Distribution,
  3. New Special Functions in GAMS,
  4. Stochastic Library,

no loops please

In (1) a piece of matlab code is presented:


There may be better way to formulate this in Matlab.  In many languages, including Matlab, R and GAMS loops are considered bad if there are alternatives. To be honest, I don’t see immediately how this particular fragment can be vectorized in Matlab.

The suggested  GAMS version is:


I strongly disagree. I believe this can be formulated without loops as:

 c(i,j,k)$(ord(k)<=plan(j,i)) = sub(j,i); 

Also notice that MATLAB will store this in a fully allocated or dense matrix (unless explicitly using a sparse matrix) while GAMS always stores things in a sparse format.


  1. Generate Data using for loop in GAMS.

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.


LXSS diagram

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


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


  1. Jack Hammons, Bash on Ubuntu on Windows,
  2. Learning about Bash on Windows Subsystem for Linux,
  3. VirtualBox,
  4. VMware,
  5. Cygwin,
  6. MINGW, Minimalist GNU for Windows,

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:


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


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:


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:


  1. Gamma distribution,
  2. 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\):


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] {
&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

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

\bbox[lightcyan,10px,border:3px solid darkblue] {
&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\}

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


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

  1. How to write this constraint in LPSolve?,
  2. Filling up a container with boxes,

Monday, July 10, 2017

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


A lot of authors.

Apparently 43,078 LPs had to be solved.

  1. Thomas Hales e.a., A Formal Proof of the Kepler Conjecture, Forum of mathematics, Pi (2017), vol 5, e2,
  2. Original message appeared in the Help-GLPK mailing list,

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


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

  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, 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??,

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.
  1. Comment by Michael Grant in

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.


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}\]


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:


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:


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:

& 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}

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:



  • \(\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:



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.


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.