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

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

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

Equity in optimization models

In optimization models, we often use an aggregate measure in the objective function, such as total profit, the sum of tardiness of jobs, and countrywide GDP. This can lead to particularly bad results for some individuals or groups.

Here is an example I have used on several occasions. 

Problem Statement

We have \(P\) persons. They must be assigned to \(M\) groups or teams. For simplicity, we can assume \(n\) is a multiple of \(m\), and the group size is \[\frac{N}{M}\] Each person \(p_1\) specifies some preferences to be placed in the same group as a person \(p_2\). A negative preference can be used to indicate that I prefer to be in a different group. Find an optimal assignment taking into account these preferences. 

Modeling surprises

Here is an example where the PuLP modeling tool goes berserk.

In standard linear programming, only \(\ge\), \(=\) and \(\le\) constraints are supported. Some tools also allow \(\ne\), which for MIP models needs to be reformulated into a disjunctive constraint. Here is an attempt to do this in PuLP [1]. PuLP does not support this relational operator in its constraints, so we would expect a meaningful error message.

Rounding inside an optimization model

In [1], the question was asked: how can I round to two decimal places inside an optimization model? I.e., \[\color{darkred}y_{i,j} = \mathbf{round}(\color{darkred}x_{i,j},2)\] To get this off my chest first: I have never encountered a situation like this. Rounding to two decimal places is more for reporting than something we want inside model equations. Given that, let me look into this modeling problem a bit more as an exercise. 

Instead of integers use binaries

In [1], a small (fragment of a) model is proposed:

High-Level Model
\[\begin{align} \min\> & \sum_i | \color{darkblue}a_i\cdot \color{darkred}x_i| \\ & \max_i |\color{darkred}x_i| = 1 \\ & \color{darkred}x_i \in \{-1,0,1\} \end{align}\]

Can we formulate this as a straight MIP? 

Math vs Programming

 A programmer writes about this blog:

(It is old, but I just came across this).

In my previous post, I just argued the other way around. To make sure: I don't hate programmers.

BTW, in quite a few programming languages for loops are very slow, and need to be replaced by something like sum(). Examples: Python, R, SQL. 

Small non-convex MINLP: Pyomo vs GAMS

 In [1], the following Pyomo model (Python fragment) is presented:

model.x = Var(name="Number of batches", domain=NonNegativeIntegers, initialize=10)                    
model.a = Var(name="Batch Size", domain=NonNegativeIntegers, bounds=(5,20))

# Objective function
def total_production(model):
    return model.x * model.a
model.total_production = Objective(rule=total_production, sense=minimize)

# Constraints
# Minimum production of the two output products
def first_material_constraint_rule(model):
    return sum(0.2 * model.a * i for i in range(1, value(model.x)+1)) >= 70
model.first_material_constraint = Constraint(rule=first_material_constraint_rule)

def second_material_constraint_rule(model):
    return sum(0.8 * model.a * i for i in range(1, value(model.x)+1)) >= 90
model.second_material_constraint = Constraint(rule=second_material_constraint_rule)

# At least one production run
def min_production_rule(model):
    return model.x >= 1
model.min_production = Constraint(rule=min_production_rule)

String Art

 

In [1], a greyscale picture is approximated by strings (lines) between points around the image. Here, I will try something similar with a formal optimization model.

Scheduling Team Meetings

This is a simple problem from [1]:

I'm rusty on constraint optimization and am looking for help in this particular

use case. There are individuals who are each member to several teams. This is a

fixed many-to-many relationship and is determined a-priori. There are 3 time

slots where the teams can be scheduled to conduct a business meeting, but if an

individual is a member of more than one team which are both meeting at a given

time slot, they'll only be able to attend one. The objective is to schedule the

teams into the time slots, minimizing the number of overlaps of individuals.

For beginners, it is often a good idea to split the task in two:

1. Formulate a mathematical model (on a piece of paper)
2. Implement the model in code

Linear Programming Nonsense?

1. Inventory balance constraints

I came accross this text [1]:

Inventory Balance Constraint

Finding common patterns

In [1], the following problem is stated:

Given a boolean matrix, with \(m\) rows and \(n\) columns, find the largest pattern of ones that is found in at least \(\color{darkblue}K\) rows. We can ignore cells where the pattern has a zero value: they don't count.

A small example [1] is given: 

row 1:[0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1]
row 2:[0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1]
row 3:[0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1]
row 4:[1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1]
row 5:[1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1]
row 6:[1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1]

 With \(K=3\), we can form a pattern with 10 nonzero elements:

 

row 1:  [0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1]
row 2:  [0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1]
row 3:  [0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1]
row 4:  [1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1]
row 5:  [1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1]
row 6:  [1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1] 
pattern:[1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1]

The pattern is shared by rows 4, 5, and 6.

Can we formulate a MIP model for this? My first attempt is as follows.

Assigning jobs to machines without overlap

 Here we consider the following problem from [1]:

• We have jobs with a given start time and completion time
• Jobs can be repeated on given days (e.g. job 1 needs to run on Monday, Wednesday, and Friday)
• We want to assign jobs to machines in such a way that there is no overlap
• The objective is to minimize the number of machines needed to execute all jobs

The planning horizon is a week, and the problems has 50 jobs.

Supplier selection: an easy MIP

This is a fairly standard supplier selection model. It was posted in [1]. It is interesting as it is a good fit for a MIP model: it is simple to model, it solves fast, and it is not obvious how to solve without an optimization model.

The problem is:

• We want to order items in different quantities from suppliers.
• Suppliers have an available inventory for these items. This can be zero.
• We can split the ordering over different suppliers.
• The cost structure is as follows:
• Shipping cost is a fixed cost per supplier.
• Item cost is a variable per-unit cost. 

High Level MIP Modeling

Most LP/MIP modeling tools stay close to what is the underpinning of a Mixed-Integer Programming problem: a system of linear equations (equalities and inequalities) plus a linear objective. Examples are Pulp and GAMS. In Constraint Programming (CP), modeling tools need to provide access to higher-level global constraints. Without these global constraints (such as the famous all-different constraint), CP solvers would not perform very well. 

A strange series

In [1] a question is posed about a somewhat strange series. Let \(i \in \{ 0,\dots,n\}\). Then we require:

• \(\color{darkred}x_i \in \{0,1,2,\dots\}\),
• \(\color{darkred}x_i\) is equal to the number of times the value \(i\) occurs in the series. In mathematical terms, one could say: \[\color{darkred}x_i = |\{j:\color{darkred}x_j=i\}|\] (here \(|S|\) is the cardinality of the set \(S\)).

An example for \(n=4\) is: 

index  0 1 2 3 4 

value  2 1 2 0 0