Thursday, August 10, 2017

Element renaming

In GAMS set elements (used for indexing) are strings. This means it is not necessary to use non-descriptive numbers as in:

set i /1*50/;

Actually, even with numbered elements, it is good practice to prefix such elements, as in:

set i /i1*i50/;

This will help when inspecting large, multi-dimensional data-structures.

As an example consider the code fragment:

set
i /1*8/
j /1*4/
k /1*6/
;

parameter p(i,j,k);
p(i,j,k) = uniform(0,1);

I see these type of numbered sets quite often. Here is why this may be a bad idea: especially after pivoting a bit, it is really difficult to see what is going on:

If we prefix the set elements with a name, we achieve much more clarity:

When I receive a model with numbered set elements the first thing I try to do is to prefix the labels.

Of course, in practical models, we often can use meaningful names, but sometimes we just have numbers as IDs. I have also seen cases where the ID is a number or a very long description. In that case, the number can be a better candidate for set element.

Exporting set elements

In GAMS a parameter is really a sparse (multi-dimensional) matrix. This corresponds quite nicely to a table in a RDBMS (relational database) or to a dataframe in R or Python. A matrix in R or Python is dense, so those data structures do not work as well for large, sparse parameters, although R has nice facilities to use matrices with row and column names.

Example

Let’s generate some data in GAMS and export to an SQLite database:

set
i 'observations' /i1*i50/
s 'scenarios' /scen1*scen3/
;

parameter results(s,i);
results('scen1',i) = normal(0,1);
results('scen2',i) = normal(-1+4*ord(i)/card(i),0.4);
results('scen3',i) = uniform(4,5);

execute "gdx2sqlite -i results.gdx -o results.db";

We can try to plot this data in R as follows:

library("RSQLite")
library("ggplot2")

# connect to the sqlite file
sqlite = dbDriver("SQLite")
db = dbConnect(sqlite, "results.db")

# retrieve data
df = dbGetQuery(db,"select * from results")

ggplot(data=df,aes(x=i,y=value,color=s)) + geom_line()

dbDisconnect(db)


s  i      value
1 scen1 i1 -0.3133429
2 scen1 i2  0.3276748
3 scen1 i3  0.4635588
4 scen1 i4 -1.8299478
5 scen1 i5 -0.7316124
6 scen1 i6 -0.9715983
Unfortunately, things don’t work exactly as expected. We see an error message:
geom_path: Each group consists of only one observation. Do you need to adjust
the group aesthetic?

And the plot does not look correct:

The reason is the variable df$i.If we convert this back to integers we are ok. To do this conversion we actually have to do two things: 1. Drop the first character ‘i’ from df$i.
2. Convert df$i from string to numeric (type casting). Interestingly we can fix this in several stages: Partial GAMS solution We can not really rename set elements in GAMS. However we can create a new parameter, a mapping set and use a sum statement to map between sets: set i 'observations' /i1*i50/ s 'scenarios' /scen1*scen3/ i2 'observations as number' /1*50/ map(i2,i) '1-1 mapping set' ; map(i2,i)$(ord(i2)=ord(i)) = yes;

parameter results(s,i);
results('scen1',i) = normal(0,1);
results('scen2',i) = normal(-1+4*ord(i)/card(i),0.4);
results('scen3',i) = uniform(4,5);

parameter results2(s,i2);
results2(s,i2) = sum(map(i2,i),results(s,i));

Now we can export the parameter results2 instead of results. Note that index i2 will still be a string when it arrives in R. The type casting to a numeric type has to be done in R or in SQL (see next sections how this can be done).

Conversion in R

The conversion to a numeric type can be done in R in just one line

df$i <- as.numeric(substring(df$i,2))

Now things look better:

Conversion in SQL

We can also perform the conversion in SQL while reading the data from the SQLite database file:

select s,
cast(substr(i,2) as numeric) as i,
value
from results

Scripting

This whole thing can be easily scripted in GAMS.

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

GammaReg(x/theta,k)

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

 scalars  k / 0.363 /  theta / 27.863 /  x / 1.5 /  p;p = gammareg(x/theta,k);display p;$FuncLibIn stodclib stodclibfunction cdfGamma / stodclib.cdfGamma /;p = cdfGamma(x,k,theta);display p; References 1. Modeling Flood Damages, http://yetanothermathprogrammingconsultant.blogspot.com/2017/07/modeling-flood-damages.html 2. Gamma Distribution, https://en.wikipedia.org/wiki/Gamma_distribution 3. New Special Functions in GAMS, http://www.amsterdamoptimization.com/pdf/specfun.pdf 4. Stochastic Library, https://www.gams.com/latest/docs/userguides/userguide/_u_g__extrinsic_functions.html#UG_ExtrinsicFunctions_StochasticLibrary 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.

References

1. Generate Data using for loop in GAMS. https://stackoverflow.com/questions/45167648/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.

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.

References

1. Jack Hammons, Bash on Ubuntu on Windows, https://msdn.microsoft.com/commandline/wsl
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:

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:

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

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

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

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

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.

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