# The facility location problem

It was originally contributed by Mathieu Tanneau (@mtanneau) and Alexis Montoison (@amontoison).

```
using JuMP
import HiGHS
import LinearAlgebra
import Plots
import Random
```

## Uncapacitated facility location

### Problem description

We are given

- A set $M=\{1, \dots, m\}$ of clients
- A set $N=\{ 1, \dots, n\}$ of sites where a facility can be built

**Decision variables** Decision variables are split into two categories:

- Binary variable $y_{j}$ indicates whether facility $j$ is built or not
- Binary variable $x_{i, j}$ indicates whether client $i$ is assigned to facility $j$

**Objective** The objective is to minimize the total cost of serving all clients. This costs breaks down into two components:

- Fixed cost of building a facility.

In this example, this cost is $f_{j} = 1, \ \forall j$.

- Cost of serving clients from the assigned facility.

In this example, the cost $c_{i, j}$ of serving client $i$ from facility $j$ is the Euclidean distance between the two.

**Constraints**

- Each customer must be served by exactly one facility
- A facility cannot serve any client unless it is open

### MILP formulation

The problem can be formulated as the following MILP:

\[\begin{aligned} \min_{x, y} \ \ \ & \sum_{i, j} c_{i, j} x_{i, j} + \sum_{j} f_{j} y_{j} \\ s.t. & \sum_{j} x_{i, j} = 1, && \forall i \in M \\ & x_{i, j} \leq y_{j}, && \forall i \in M, j \in N \\ & x_{i, j}, y_{j} \in \{0, 1\}, && \forall i \in M, j \in N \end{aligned}\]

where the first set of constraints ensures that each client is served exactly once, and the second set of constraints ensures that no client is served from an unopened facility.

### Problem data

```
Random.seed!(314)
# number of clients
m = 12
# number of facility locations
n = 5
# Clients' locations
Xc = rand(m)
Yc = rand(m)
# Facilities' potential locations
Xf = rand(n)
Yf = rand(n)
# Fixed costs
f = ones(n);
# Distance
c = zeros(m, n)
for i in 1:m
for j in 1:n
c[i, j] = LinearAlgebra.norm([Xc[i] - Xf[j], Yc[i] - Yf[j]], 2)
end
end
```

Display the data

```
Plots.scatter(
Xc,
Yc,
label = "Clients",
markershape = :circle,
markercolor = :blue,
)
Plots.scatter!(
Xf,
Yf,
label = "Facility",
markershape = :square,
markercolor = :white,
markersize = 6,
markerstrokecolor = :red,
markerstrokewidth = 2,
)
```

### JuMP implementation

Create a JuMP model

`ufl = Model(HiGHS.Optimizer)`

```
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: HiGHS
```

Variables

```
@variable(ufl, y[1:n], Bin);
@variable(ufl, x[1:m, 1:n], Bin);
```

Each client is served exactly once

`@constraint(ufl, client_service[i in 1:m], sum(x[i, j] for j in 1:n) == 1);`

A facility must be open to serve a client

`@constraint(ufl, open_facility[i in 1:m, j in 1:n], x[i, j] <= y[j]);`

Objective

`@objective(ufl, Min, f'y + sum(c .* x));`

Solve the uncapacitated facility location problem with HiGHS

`optimize!(ufl)`

```
Presolving model
72 rows, 65 cols, 180 nonzeros
72 rows, 65 cols, 180 nonzeros
Solving MIP model with:
72 rows
65 cols (65 binary, 0 integer, 0 implied int., 0 continuous)
180 nonzeros
( 0.0s) Starting symmetry detection
( 0.0s) No symmetry present
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 0 inf inf 0 0 0 0 0.0s
T 0 0 0 0.00% 0 5.22641797 100.00% 0 0 0 32 0.0s
Solving report
Status Optimal
Primal bound 5.22641797047
Dual bound 5.22641797047
Gap 0% (tolerance: 0.01%)
Solution status feasible
5.22641797047 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.00 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 32 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
```

`println("Optimal value: ", objective_value(ufl))`

`Optimal value: 5.226417970467934`

### Visualizing the solution

The threshold 1e-5 ensure that edges between clients and facilities are drawn when x[i, j] ≈ 1.

```
x_ = value.(x) .> 1 - 1e-5
y_ = value.(y) .> 1 - 1e-5
```

```
5-element BitVector:
0
0
1
1
0
```

Display clients

```
p = Plots.scatter(
Xc,
Yc,
markershape = :circle,
markercolor = :blue,
label = nothing,
)
```

Show open facility

```
mc = [(y_[j] ? :red : :white) for j in 1:n]
Plots.scatter!(
Xf,
Yf,
markershape = :square,
markercolor = mc,
markersize = 6,
markerstrokecolor = :red,
markerstrokewidth = 2,
label = nothing,
)
```

Show client-facility assignment

```
for i in 1:m
for j in 1:n
if x_[i, j] == 1
Plots.plot!(
[Xc[i], Xf[j]],
[Yc[i], Yf[j]],
color = :black,
label = nothing,
)
end
end
end
p
```

## Capacitated facility location

### Problem formulation

The capacitated variant introduces a capacity constraint on each facility, i.e., clients have a certain level of demand to be served, while each facility only has finite capacity which cannot be exceeded.

Specifically,

- The demand of client $i$ is denoted by $a_{i} \geq 0$
- The capacity of facility $j$ is denoted by $q_{j} \geq 0$

The capacity constraints then write

\[\begin{aligned} \sum_{i} a_{i} x_{i, j} &\leq q_{j} y_{j} && \forall j \in N \end{aligned}\]

Note that, if $y_{j}$ is set to $0$, the capacity constraint above automatically forces $x_{i, j}$ to $0$.

Thus, the capacitated facility location can be formulated as follows

\[\begin{aligned} \min_{x, y} \ \ \ & \sum_{i, j} c_{i, j} x_{i, j} + \sum_{j} f_{j} y_{j} \\ s.t. & \sum_{j} x_{i, j} = 1, && \forall i \in M \\ & \sum_{i} a_{i} x_{i, j} \leq q_{j} y_{j}, && \forall j \in N \\ & x_{i, j}, y_{j} \in \{0, 1\}, && \forall i \in M, j \in N \end{aligned}\]

For simplicity, we will assume that there is enough capacity to serve the demand, i.e., there exists at least one feasible solution.

Demands

`a = rand(1:3, m);`

Capacities

`q = rand(5:10, n);`

Display the data

```
Plots.scatter(
Xc,
Yc,
label = nothing,
markershape = :circle,
markercolor = :blue,
markersize = 2 .* (2 .+ a),
)
Plots.scatter!(
Xf,
Yf,
label = nothing,
markershape = :rect,
markercolor = :white,
markersize = q,
markerstrokecolor = :red,
markerstrokewidth = 2,
)
```

### JuMP implementation

Create a JuMP model

`cfl = Model(HiGHS.Optimizer)`

```
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: HiGHS
```

Variables

```
@variable(cfl, y[1:n], Bin);
@variable(cfl, x[1:m, 1:n], Bin);
```

Each client is served exactly once

`@constraint(cfl, client_service[i in 1:m], sum(x[i, :]) == 1);`

Capacity constraint

`@constraint(cfl, capacity, x'a .<= (q .* y));`

Objective

`@objective(cfl, Min, f'y + sum(c .* x));`

Solve the problem

`optimize!(cfl)`

```
Presolving model
17 rows, 65 cols, 125 nonzeros
17 rows, 65 cols, 125 nonzeros
Solving MIP model with:
17 rows
65 cols (65 binary, 0 integer, 0 implied int., 0 continuous)
125 nonzeros
( 0.0s) Starting symmetry detection
( 0.0s) No symmetry present
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 0 inf inf 0 0 0 0 0.0s
S 0 0 0 0.00% 0 7.54947525 100.00% 0 0 0 0 0.0s
0 0 0 0.00% 5.636380012 7.54947525 25.34% 0 0 3 21 0.0s
S 0 0 0 0.00% 5.636380012 7.373597157 23.56% 12 1 3 21 0.0s
C 0 0 0 0.00% 6.078302661 7.204177775 15.63% 96 14 6 37 0.0s
L 0 0 0 0.00% 6.173715063 6.173715063 0.00% 276 23 6 50 0.0s
Solving report
Status Optimal
Primal bound 6.17371506253
Dual bound 6.17371506253
Gap 0% (tolerance: 0.01%)
Solution status feasible
6.17371506253 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.02 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 69 (total)
0 (strong br.)
29 (separation)
19 (heuristics)
```

`println("Optimal value: ", objective_value(cfl))`

`Optimal value: 6.17371506253207`

### Visualizing the solution

The threshold 1e-5 ensure that edges between clients and facilities are drawn when x[i, j] ≈ 1.

```
x_ = value.(x) .> 1 - 1e-5;
y_ = value.(y) .> 1 - 1e-5;
```

Display the solution

```
p = Plots.scatter(
Xc,
Yc,
label = nothing,
markershape = :circle,
markercolor = :blue,
markersize = 2 .* (2 .+ a),
)
mc = [(y_[j] ? :red : :white) for j in 1:n]
Plots.scatter!(
Xf,
Yf,
label = nothing,
markershape = :rect,
markercolor = mc,
markersize = q,
markerstrokecolor = :red,
markerstrokewidth = 2,
)
```

Show client-facility assignment

```
for i in 1:m
for j in 1:n
if x_[i, j] == 1
Plots.plot!(
[Xc[i], Xf[j]],
[Yc[i], Yf[j]],
color = :black,
label = nothing,
)
break
end
end
end
p
```

This tutorial was generated using Literate.jl. View the source `.jl`

file on GitHub.