In [1]:
using JuMP, GLPK

Location of Warehouses

This is based on Example 5.1 (Location of Warehouses) from Applied Linear Programming

A firm has 5 distribution centers and we want to dertermine which subset of these should serve as a site for a warehouse. The goal is the build a minimum number of warehouses that can cover all distribution centers so that every warehouse is within 10 miles of each distribution center.

Per the problem statement, $m$ is the number of distribution centers.

We are given a table of distances between distribution centers, $D$:

In [2]:
m = 5
max_miles = 10

D = [0 10 15 20 18;
     10 0 20 15 10;
     15 20 0 8 17;
     20 15 8 0 5;
     18 10 17 5 0
    ]
Out[2]:
5×5 Array{Int64,2}:
  0  10  15  20  18
 10   0  20  15  10
 15  20   0   8  17
 20  15   8   0   5
 18  10  17   5   0

For example, it is 18 miles between distribution centers 1 (column 1) and 5 (row 5).

To convert this to a binary coverage vector $A$, we convert each distance into a binary variable indicating whether the distribution centers are 10 or fewer miles from one another:

In [3]:
A = [Int(D[i, j] <= max_miles) for i=1:m, j=1:m]
Out[3]:
5×5 Array{Int64,2}:
 1  1  0  0  0
 1  1  0  0  1
 0  0  1  1  0
 0  0  1  1  1
 0  1  0  1  1

Now we can model this problem using the JuMP package and the (open source) Cbc solver:

In [9]:
model = Model(with_optimizer(GLPK.GLPKOptimizerMIP))

# decision variable (binary): whether to build warehouse near distribution center i
@variable(model, y[1:m], Bin)

# Objective: minimize number of warehouses
@objective(model, Min, sum(y))

# Constraint: has to cover all warehouses
# (.>= is the element-wise dot comparison operator)
@constraint(model, A*y .>= 1)
Out[9]:
5-element Array{JuMP.ConstraintRef{JuMP.Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}},1}:
 y[1] + y[2] ≥ 1.0       
 y[1] + y[2] + y[5] ≥ 1.0
 y[3] + y[4] ≥ 1.0       
 y[3] + y[4] + y[5] ≥ 1.0
 y[2] + y[4] + y[5] ≥ 1.0

We have an additional constraint that at least 1 warehouse should be within 10 miles of distribution center 1, but our activity matrix $A$ already covers that, so technically we do not need this explicit constraint.

In [10]:
@constraint(model, y[1] + y[2] >= 1)
Out[10]:
$ y_{1} + y_{2} \geq 1.0 $
In [11]:
# Solve problem using MIP solver
JuMP.optimize(model)
In [16]:
println("Total # of warehouses: ", JuMP.objectivevalue(model))

println("Build warehouses at distribution center(s):")

for i=1:m 
    if JuMP.resultvalue(y[i]) == 1 
        println("Warehouse $i")
    end
end
Total # of warehouses: 2.0
Build warehouses at distribution center(s):
Warehouse 2
Warehouse 3