Revisiting a continuous facility location problem

I am revisiting here a problem from [1]:

We have \(n\) demand points and their locations. How many facilities do we need to service these customers? And where do we place them? We have a restriction: there is a maximum distance between customer and facility.   

Data

We randomly generated 75 demand points inside the \([0,1]\times[0,1]\) square. The maximum allowed distance between a demand point and a facility is 0.25. 

Random Demand Points

Nonconvex problem: local vs multistart vs global

In [1] a somewhat abstract non-convex problem is given:

\[\begin{align}\min_x & - x_1^2 - x_2^2 - x_3^2 - x_4^2\\ & Ax \le b \end{align}\]

This is a nonconvex objective with some linear constraints. Of course, the easiest way to gain some insight into how to solve this is to perform some quick and dirty tests with a global QP solver designed for these types of models. As the constraints are linear, a solver like Cplex or Gurobi may be a good starting point (pun intended).

If you have access to a local NLP solver, it may not always be easy to find a good starting point. One possible approach is to use a multistart algorithm (some NLP solvers, like Baron and Knitro, have this built-in). This will not guarantee a global solution (and most likely, it doesn't give you one), but at least we can prevent really bad solutions.  

Towers of Hanoi: inventory and network formulation

The towers of Hanoi problem [1] is a famous puzzle demonstrating recursion. The task is to move a stack of disks from one peg to another. We can only move one disk at a time, and we need to obey the rule that larger disks can never be on top of a smaller disk. The moves for the 3 disk problem are:

Wolf, goat and cabbage problem: MIP and network model

The Wolf, goat, and cabbage problem can be stated as [1]:

A farmer with a wolf, a goat, and a cabbage must cross a river by boat. The boat can carry only the farmer and a single item. If left unattended together, the wolf would eat the goat, or the goat would eat the cabbage. How can they cross the river without anything being eaten?

A long transcontinental flight was a good opportunity to try to attack this problem. Here are some approaches to model this. Not at all a very useful or practical model, but still interesting, I think (although I may be in a small minority here). 

Inventory model

In this model, we keep track of the inventory of items (wolf, goat, cabbage). We assume the starting inventory is: all items are on the left bank of the river. The final inventory should be: all items are on the right bank. In this model, I assume the following numbering scheme:

----     31 SET trip  trips

trip1 ,    trip2 ,    trip3 ,    trip4 ,    trip5 ,    trip6 ,    trip7 ,    trip8 ,    trip9 ,    trip10


----     31 SET dir  direction of trip

L->R,    R->L


----     31 SET tripDir  trip direction combos

              L->R        R->L

trip1          YES
trip2                      YES
trip3          YES
trip4                      YES
trip5          YES
trip6                      YES
trip7          YES
trip8                      YES
trip9          YES
trip10                     YES

Small MIP, proving optimality is difficult

This is a simple, smallish MIP model that is difficult to solve to proven optimality. The problem is stated as [1]:

I have grid of dimensions H and W of boolean variables. The only constraint is that if a variable is false then at least one of the adjacent variables (top, right, left, bottom, diagonals don't count) must be true. The goal is to minimize the number of true values in the grid.

A high-level mathematical representation of this can be:

High-level model
\[\begin{align}\min&\sum_{i,j}\color{darkred}x_{i,j}\\ & \color{darkred}x_{i,j}=0 \implies \color{darkred}x_{i-1,j} + \color{darkred}x_{i+1,j} + \color{darkred}x_{i,j-1} + \color{darkred}x_{i,j+1} \ge 1 \\ & \color{darkred}x_{i,j}\in \{0,1\} \end{align}\]

Programming vs Modeling

I found another interesting GAMS model. Here, I wanted to demonstrate the difference between programming (implementing algorithms in a programming language) and modeling (building models and implementing them in a modeling tool). The problem is solving a shortest path problem given a sparse network with random costs. We compare three strategies:

1. Implement Dijkstra's shortest path algorithm in GAMS,
2. implement Dijkstra's shortest path algorithm in Python and
3. use a simple LP model.  

One of the practical problems is how to represent and display the optimal shortest path in a meaningful and compact way. In some cases, this is more difficult than implementing the pure shortest path algorithm or model. 

Network models with node capacities

Min-cost flow network model

A standard min-cost flow network model, formulated as a linear programming model, can look like:

Model 1: standard min-cost flow network model
\[\begin{align}\min& \sum_{(i,j)\in A}\color{darkblue}c_{i,j}\cdot\color{darkred}f_{i,j}\\ & \sum_{j|(j,i)\in A}\color{darkred}f_{j,i} + \color{darkblue}{\mathit{supply}}_i = \sum_{j|(i,j)\in A}\color{darkred}f_{i,j} + \color{darkblue}{\mathit{demand}}_i && \forall i \\ & \color{darkred}f_{i,j} \in [0,\color{darkblue}{\mathit{capacity}}_{i,j}] \end{align}\]

Sorting using a MIP model

This is not, per se, very useful, but sorting a parameter inside a MIP model is not very easy for MIP solvers. Obviously, if you need a sorted parameter in the model, it is better to use a sorting algorithm. But useless models can still be interesting.

I use: 

    Input: a 1-dimensional parameter with values.
    Output: a 1-dimensional variable with the values sorted in ascending order.

We can implement this with a permutation matrix \(\color{darkred}X\), which is a permuted identity matrix. In a MIP context, this becomes a binary variable with some assignment constraints.

MIP Model for sorting \(p_i\)
\[\begin{align}\min\>&\color{darkred}z=0\\& \sum_i \color{darkred}x_{i,j} = 1&&\forall j\\ & \sum_j \color{darkred}x_{i,j} = 1&&\forall i\\ & \color{darkred}y_i = \sum_j \color{darkred}x_{i,j}\cdot\color{darkblue}p_j\\& \color{darkred}y_i \ge \color{darkred}y_{i-1}\\ & \color{darkred}x_{i,j} \in \{0,1\}\end{align}\]

PuLP surprises

Formulating optimization models inside traditional programming languages such as Python is very popular. The main tool the developers use to make this possible is operator overloading. There are cases, where we can write code that looks somewhat reasonable, is accepted and processed without any warning or error messages, but is total nonsense. It is rather difficult with this approach to make things airtight. Especially error handling. In [1], we see a good example. I have created a small fragment here that illustrates the problem.

Non-convex Quadratic Integer Programming

Here, I want to revisit a particular model from [1]:

Model 3: Quadratic Preemptive Model

\[\begin{align}\max\>&\color{darkred}z_{model3}=\sum_p \color{darkred}z_p \\ & \color{darkred}z_p = \sum_{p',g} \color{darkblue}{\mathit{pref}}_{p,p'}\cdot \color{darkred}{\mathit{assign}}_{p,g}\cdot\color{darkred}{\mathit{assign}}_{p',g} \\ & \color{darkblue}z_{model2}^* \le \color{darkred}z_p & \forall p\\ & \sum_g \color{darkred}{\mathit{assign}}_{p,g} = 1 & \forall p \\ & \sum_p \color{darkred}{\mathit{assign}}_{p,g} = \color{darkblue}{\mathit{groupSize}} & \forall g \\& \color{darkred}{\mathit{assign}}_{p,g} \in \{0,1\}\end{align}\]