4. Networks 1: Assignment and Transportation Problems¶

ISE 754, Fall 2024

Package Used: No new packages used.

Many problems can be modeled as networks, where the nodes of the network represent entities and the arcs represent either flows or relationships, for example. All of the network models that we will be covering can be solved as linear programs. In the past, modeling a problem as a network allowed network-specific procedures to be used that were more efficient than solving a general linear programming formulation of the problem, but now even large-scale linear programs can be solved efficiently, and the use of network-specific procedures are not necessary for many problems.

1. Linear Assignment Problem (LAP)¶

The linear assignment problem can be modeled as a complete bipartite network, which is a network consisting of two sets of nodes with all arcs from one set to the other. The first set of "supply" nodes represent an entity that can be assigned to one of the second set of "demand" nodes. Each arc represents the cost of making the assignment, and the problem is solved by determining the least-cost assignment. Although the Hungarian method is a specialized procedure for the LAP, modeling the LAP as a network flow problem and then solving it as an LP works is also effective and easier to implement.

Ex: Unter's Passenger-to-Car Assignment¶

An example of the linear assignment problem would be its use by a ride-hailing service to assign available cars to passengers wanting a ride. One way of doing the assignment is to make the assignment based on which car is closest to any of the passengers and then continue until all the assignments are complete. Unfortunately, this simple greedy procedure does not, in general, provide an optimal assignment; instead, an optimal assignment can be determined modeled as a network flow problem, where one unit of supply provided by a car can be used to satisfy one unit of demand required for each passenger, where the arc connecting each card to each passenger represents the cost (distance) of the car from the passenger.

image.png

In [1]:
C = [4 8 5; 2 1 9; 6 3 7]     # Car-to-passenger arc cost (distance) matrix
Out[1]:
3×3 Matrix{Int64}:
 4  8  5
 2  1  9
 6  3  7
In [2]:
idx = argmin(C)               # Find closest passenger
Out[2]:
CartesianIndex(2, 2)
In [3]:
C[idx]
Out[3]:
1
In [4]:
C[idx[1],:] .= maximum(C)+1   # Assign large value so that assignment not made again
C[:,idx[2]] .= maximum(C)+1
C
Out[4]:
3×3 Matrix{Int64}:
  4  11   5
 10  11  10
  6  11   7
In [5]:
idx = argmin(C)               # Find next closest passenger
Out[5]:
CartesianIndex(1, 1)
In [6]:
C[idx]
Out[6]:
4
In [7]:
C[idx[1],:] .= maximum(C)+1
C[:,idx[2]] .= maximum(C)+1
C
Out[7]:
3×3 Matrix{Int64}:
 13  12  12
 13  11  10
 13  11   7
In [8]:
idx = argmin(C)
Out[8]:
CartesianIndex(3, 3)
In [9]:
C[idx]
Out[9]:
7
In [10]:
TC = 1 + 4 + 7                # Final cost of assignment
Out[10]:
12

This procedure can be made into a function to allow it to be easily reused:

In [11]:
C = [4 8 5; 2 1 9; 6 3 7]            # Need to re-define matrix since it was changed above
function greedyassignment(C)
    C = copy(C)                      # Copy-by-value instead of default copy-by-reference
    X = zeros(size(C))
    BigM = maximum(C) + 1               # Large value so that assignment not made again
    for i = 1:size(C,2)
        idx = argmin(C)
        X[idx] = 1
        C[idx[1],:] .= BigM
        C[:,idx[2]] .= BigM
    end
    return X
end
X = greedyassignment(C)
Out[11]:
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
In [12]:
TC = sum(C.*X)
Out[12]:
12.0

What if the number of cars does not equal the number of passengers?

In [13]:
greedyassignment(C[1:2,:])   # Not feasible to have one car serve multiple passengers
Out[13]:
2×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
In [14]:
greedyassignment(C[:,1:2])   # Feasible to have extra car not assigned
Out[14]:
3×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0
 0.0  0.0

The LAP can be modeled as a network flow problem and then solving as an LP:

Linear Assignment Problem (LAP):

$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{i \in \text{Cars}} \sum_{j \in \text{Pass}} c_{ij} x_{i,j} \\ \quad \mbox{subject to} \quad \sum_{j \in \text{Pass}} x_{ij} &\le 1\mbox{,} &\quad i \in \text{Cars}\mbox{,}\quad\mbox{: each car is assigned up to one passenger}\\ \sum_{i \in \text{Cars}} x_{ij} &= 1\mbox{,} &\quad j \in \text{Pass}\mbox{,}\quad\mbox{: each passenger is assigned a car}\\ x_{i,j} &\ge 0 \end{eqnarray*} $

In [15]:
using JuMP, Cbc

Cars, Pass = 1:size(C,1), 1:size(C,2)
model = Model(Cbc.Optimizer)
@variable(model, 0 <= x[Cars,Pass])
@objective(model, Min, sum(C[i,j]*x[i,j] for i ∈ Cars, j ∈ Pass))
@constraint(model, [i ∈ Cars], sum(x[i,j] for j ∈ Pass) <= 1 )
@constraint(model, [j ∈ Pass], sum(x[i,j] for i ∈ Cars) == 1 )
print(model)
optimize!(model)
TCᵒ, Xᵒ = objective_value(model), Array(value.(x))
$$ \begin{aligned} \min\quad & 4 x_{1,1} + 8 x_{1,2} + 5 x_{1,3} + 2 x_{2,1} + x_{2,2} + 9 x_{2,3} + 6 x_{3,1} + 3 x_{3,2} + 7 x_{3,3}\\ \text{Subject to} \quad & x_{1,1} + x_{2,1} + x_{3,1} = 1\\ & x_{1,2} + x_{2,2} + x_{3,2} = 1\\ & x_{1,3} + x_{2,3} + x_{3,3} = 1\\ & x_{1,1} + x_{1,2} + x_{1,3} \leq 1\\ & x_{2,1} + x_{2,2} + x_{2,3} \leq 1\\ & x_{3,1} + x_{3,2} + x_{3,3} \leq 1\\ & x_{1,1} \geq 0\\ & x_{2,1} \geq 0\\ & x_{3,1} \geq 0\\ & x_{1,2} \geq 0\\ & x_{2,2} \geq 0\\ & x_{3,2} \geq 0\\ & x_{1,3} \geq 0\\ & x_{2,3} \geq 0\\ & x_{3,3} \geq 0\\ \end{aligned} $$ 
Presolve 6 (0) rows, 9 (0) columns and 18 (0) elements
0  Obj 0 Primal inf 2.9999997 (3)
4  Obj 10
Optimal - objective value 10
Optimal objective 10 - 4 iterations time 0.002
Out[15]:
(10.0, [0.0 0.0 1.0; 1.0 0.0 0.0; 0.0 1.0 0.0])
In [16]:
Xᵒ
Out[16]:
3×3 Matrix{Float64}:
 0.0  0.0  1.0
 1.0  0.0  0.0
 0.0  1.0  0.0

Ex: Distance as Cost for Ex 1¶

Given the (x, y) coordinates of each car and passenger, determine the assignment based on minimizing the total Euclidean distance of the cars from their passengers.

In [17]:
using Random

d2(x1, x2) = length(x1) == length(x2) ? sqrt(sum((x1 .- x2).^2)) : 
    error("Inputs not same length.")

D2(X₁, X₂) = [d2(i, j) for i in eachrow(X₁), j in eachrow(X₂)]

function randX(P, n)
    (xmn, xmx), (ymn, ymx) = extrema(P, dims=1)
    return [xmn ymn] .+ rand(n, 2) .* [(xmx - xmn) (ymx - ymn)]
end

Random.seed!(87345)

n = 3
Xc = randX([0 0; 100 100], n)                     # Car cooridinates
Xp = randX([0 0; 100 100], n)                     # Passenger cooridinates
D = D2(Xc, Xp)
Out[17]:
3×3 Matrix{Float64}:
 102.045    100.488   50.0907
  78.7744    92.9669  36.1811
   2.13914   60.4517  53.4778
In [18]:
X = greedyassignment(D)
Out[18]:
3×3 Matrix{Float64}:
 0.0  1.0  0.0
 0.0  0.0  1.0
 1.0  0.0  0.0
In [19]:
sum(D.*X)
Out[19]:
138.80817459089374
In [20]:
using JuMP, Cbc

function LAP(c)
    N₁, N₂ = 1:size(c,1), 1:size(c,2)
    model = Model(Cbc.Optimizer)
    @variable(model, 0 <= x[N₁,N₂])
    @objective(model, Min, sum(c[i,j] * x[i, j] for i ∈ N₁, j ∈ N₂))
    @constraint(model, [i ∈ N₁], sum(x[i,j] for j ∈ N₂) <= 1)
    @constraint(model, [j ∈ N₂], sum(x[i,j] for i ∈ N₁) == 1)
    optimize!(model)
    if termination_status(model) == JuMP.OPTIMAL
        return objective_value(model), Array(value.(x))
    end
    println(termination_status(model))
end

TCᵒ, Xᵒ = LAP(D)

Presolve 6 (0) rows, 9 (0) columns and 18 (0) elements
0  Obj 0 Primal inf 2.9999997 (3)
5  Obj 138.80817
Optimal - objective value 138.80817
Optimal objective 138.8081746 - 5 iterations time 0.002
Out[20]:
(138.80817459089374, [0.0 1.0 0.0; 0.0 0.0 1.0; 1.0 0.0 0.0])

What if the number of cars does not equal the number of passengers?

In [21]:
LAP(D[1:2,:])   # Not feasible to have one car serve multiple passengers

INFEASIBLE_OR_UNBOUNDED
Presolve 0 (-5) rows, 0 (-6) columns and 0 (-12) elements
Optimal - objective value 222.96446
After Postsolve, objective 222.96446, infeasibilities - dual 0 (0), primal 1.9999998 (2)
Presolved model was optimal, full model needs cleaning up
0  Obj 222.96446 Primal inf 1.9999998 (2)
Primal infeasible - objective value 229.35305
PrimalInfeasible objective 229.3530527 - 1 iterations time 0.002, Presolve 0.00
In [22]:
LAP(D[:,1:2])   # Feasible to have extra car not assigned

Presolve 5 (0) rows, 6 (0) columns and 12 (0) elements
0  Obj 0 Primal inf 1.9999998 (2)
3  Obj 95.10604
Optimal - objective value 95.10604
Optimal objective 95.10604007 - 3 iterations time 0.002
Out[22]:
(95.10604007130576, [0.0 0.0; 0.0 1.0; 1.0 0.0])

Ex: Ten Random Cars and Passengers¶

In [23]:
Random.seed!(2342)

n = 10
XYc = randX([0 0; 100 100], n)                     # Car cooridinates
XYp = randX([0 0; 100 100], n)                     # Passenger cooridinates
D = D2(XYc, XYp)
X = greedyassignment(D)
TC = sum(D.*X)
Out[23]:
224.3157219777679
In [24]:
X = greedyassignment(D)
TC = sum(D.*X)
Out[24]:
224.3157219777679
In [25]:
using CairoMakie

dcf() = display(current_figure())

scatter(XYc[:,1],XYc[:,2], label="Cars")
scatter!(XYp[:,1],XYp[:,2], label="Passengers")

idx = findall(X .> 0)
for i in 1:length(idx)
    x = [XYc[idx[i][1], 1], XYp[idx[i][2], 1]]
    y = [XYc[idx[i][1], 2], XYp[idx[i][2], 2]]
    lines!(x, y)
end

axislegend(position=:rb)

dcf()
No description has been provided for this image
Out[25]:
CairoMakie.Screen{IMAGE}
In [26]:
TCᵒ, Xᵒ = LAP(D)

Presolve 20 (0) rows, 100 (0) columns and 200 (0) elements
0  Obj 0 Primal inf 9.999999 (10)
20  Obj 215.37037
Optimal - objective value 215.37037
Optimal objective 215.3703654 - 20 iterations time 0.002
Out[26]:
(215.37036543420638, [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.0])
In [27]:
scatter(XYc[:,1],XYc[:,2], label="Cars")
scatter!(XYp[:,1],XYp[:,2], label="Passengers")

idx = findall(Xᵒ .> 0)
for i in 1:length(idx)
    x = [XYc[idx[i][1], 1], XYp[idx[i][2], 1]]
    y = [XYc[idx[i][1], 2], XYp[idx[i][2], 2]]
    lines!(x, y)
end

axislegend(position=:rb)

dcf();
No description has been provided for this image

2. Transportation Problem¶

The transportation problem is a generalization of the linear assignment problem (LAP). As in the LAP, nodes are separated into supply and demand nodes with arcs representing the associated cost per unit of flow, connecting each supply node to all of the demand nodes. Unlike the LAP, in which each supply node can provide only one unit and each demand node can require only one unit, the transportation problem supply nodes can provide multiple units, and the demand nodes can require multiple units.

The transportation problem can be modeled as a network flow problem and then solving as an LP:

Transportation Problem:

$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{i \in S} \sum_{j \in D} c_{i,j} x_{i,j} \\ \quad \mbox{subject to} \quad \sum_{j \in D} x_{i,j} &\le& s_i\mbox{,} &\quad i \in S\\ \sum_{i \in S} x_{i,j} &=& d_j\mbox{,} &\quad j \in D\\ x_i &\ge& 0 \mbox{,} \end{eqnarray*} $

Ex: Transportation Problem with Excess Supply¶

image.png

In [28]:
C = [8 6 10 9; 9 12 13 7; 14 9 16 5]
d = [45, 20, 30, 30]
s = [55, 50, 40]
sum(s) - sum(d)
Out[28]:
20
In [29]:
using JuMP, Cbc

function trans(c,s,d)
    S, D = 1:length(s), 1:length(d)
    model = Model(Cbc.Optimizer)
    @variable(model, 0 <= x[S,D])
    @objective(model, Min, sum(c[i,j]*x[i,j] for i ∈ S, j ∈ D))
    @constraint(model, [i ∈ S], sum(x[i,j] for j ∈ D) <= s[i])
    @constraint(model, [j ∈ D], sum(x[i,j] for i ∈ S) == d[j])
    optimize!(model)
    if termination_status(model) == JuMP.OPTIMAL
        return objective_value(model), Array(value.(x))
    end
    println(termination_status(model))
end

TCᵒ, Xᵒ = trans(C,s,d)

Presolve 7 (0) rows, 12 (0) columns and 24 (0) elements
0  Obj 0 Primal inf 125 (4)
5  Obj 970
Optimal - objective value 970
Optimal objective 970 - 5 iterations time 0.002
Out[29]:
(970.0, [5.0 20.0 30.0 0.0; 40.0 0.0 0.0 0.0; 0.0 0.0 0.0 30.0])
In [30]:
Matrix([Xᵒ sum(Xᵒ, dims=2) s; sum(Xᵒ, dims=1) 0 0; d' 0 0 ])
Out[30]:
5×6 Matrix{Float64}:
  5.0  20.0  30.0   0.0  55.0  55.0
 40.0   0.0   0.0   0.0  40.0  50.0
  0.0   0.0   0.0  30.0  30.0  40.0
 45.0  20.0  30.0  30.0   0.0   0.0
 45.0  20.0  30.0  30.0   0.0   0.0

What if the total supply is less than the total demand?

In [31]:
s = [25, 50, 40]
@show sum(s) - sum(d)
trans(C,s,d)

sum(s) - sum(d) = -10
INFEASIBLE_OR_UNBOUNDED
Presolve 7 (0) rows, 12 (0) columns and 24 (0) elements
0  Obj 0 Primal inf 125 (4)
7  Obj 1050 Primal inf 9.9999999 (1)
Primal infeasible - objective value 1050
PrimalInfeasible objective 1050 - 7 iterations time 0.002

Ex: Vaccine Supply¶

Tomorrow, 130, 82, 215, and 175 doses of vaccine will be needed at clinics located at coordinates (3,2), (2,3), (4,2), and (3,4), respectively, and 180, 300, and 250 doses are available for shipment overnight from distribution centers located at coordinates (2,2), (1,1), and (2,4), respectively. Assuming the cost of supplying each clinic from a distribution center is proportional to the Euclidean distance between them, determine how much each distribution center should send to each clinic.

In [32]:
s = [180, 300, 250]
d = [130, 82, 215, 175]
XYdc = [2 2; 1 1; 2 4]
XYclinic = [3 2; 2 3; 4 2; 3 4]
D = D2(XYdc, XYclinic)
Out[32]:
3×4 Matrix{Float64}:
 1.0      1.0      2.0      2.23607
 2.23607  2.23607  3.16228  3.60555
 2.23607  1.0      2.82843  1.0
In [33]:
TCᵒ, Xᵒ = trans(D,s,d)

Presolve 7 (0) rows, 12 (0) columns and 24 (0) elements
0  Obj 0 Primal inf 602 (4)
7  Obj 1016.9118
Optimal - objective value 1016.9118
Optimal objective 1016.911758 - 7 iterations time 0.002
Out[33]:
(1016.9117575489613, [130.0 7.0 43.0 0.0; 0.0 0.0 172.0 0.0; 0.0 75.0 0.0 175.0])
In [34]:
Matrix([Xᵒ sum(Xᵒ, dims=2) s; sum(Xᵒ, dims=1) 0 0; d' 0 0 ])
Out[34]:
5×6 Matrix{Float64}:
 130.0   7.0   43.0    0.0  180.0  180.0
   0.0   0.0  172.0    0.0  172.0  300.0
   0.0  75.0    0.0  175.0  250.0  250.0
 130.0  82.0  215.0  175.0    0.0    0.0
 130.0  82.0  215.0  175.0    0.0    0.0

Ex: ALA Allocation¶

When NFs have capacity constraints, the allocation for the ALA procedure can be determined by solving the transportation problem. In Loc 4, the population of NC 100K cities was allocated to Raleigh and Charlotte based on closeness. All of the population of each city was allocated to either Raleigh or Charlotte. If the total population allocated was constrained, the population of some of the cities would have to be split between Raleigh and Charlotte. In the following, the capacity of Raleigh and Charlotte is 1.5 million people.

In [35]:
using Logjam.DataTools, DataFrames
using SparseArrays

function dgc(xy₁, xy₂; unit=:mi)
    length(xy₁) == length(xy₂) == 2 || error("Inputs must have length 2.")
    unit in [:mi, :km] || error("Unit must be :mi or :km")

    Δx, Δy = xy₂[1] - xy₁[1], xy₂[2] - xy₁[2]
    a = sind(Δy / 2)^2 + cosd(xy₁[2]) * cosd(xy₂[2]) * sind(Δx / 2)^2
    2 * asin(min(sqrt(a), 1.0)) * (unit == :mi ? 3958.75 : 6371.00)
end

Dgc(X₁, X₂) = [dgc(i, j) for i in eachrow(X₁), j in eachrow(X₂)]

function name2lonlat(name, df)
    idx = findfirst(r -> r.NAME == name, eachrow(df))   # Find first match
    if idx === nothing
        error("'$name' not found in $df")
    end
    return collect(df[idx, [:LON, :LAT]])
end

df = filter(r -> (r.STFIP == st2fips(:NC)) && (r.POP > 100_000), usplace())
select!(df, :NAME, :LAT, :LON, :POP)
P = hcat(df.LON, df.LAT)
xyC = name2lonlat("Charlotte", df)
xyR = name2lonlat("Raleigh", df)
X = vcat(xyC', xyR')
w = df.POP
Out[35]:
10-element Vector{Int64}:
 174721
 874579
 105240
 283506
 208501
 299035
 114059
 467665
 115451
 249545

Capacity Unconstrained¶

In [36]:
n, m = 2, 10

alloc(X) = sparse([argmin(c) for c in eachcol(Dgc(X, P))], 1:m, w, n, m)

W = alloc(X)'
Out[36]:
10×2 LinearAlgebra.Adjoint{Int64, SparseMatrixCSC{Int64, Int64}} with 10 stored entries:
      ⋅  174721
 874579       ⋅
 105240       ⋅
      ⋅  283506
      ⋅  208501
      ⋅  299035
 114059       ⋅
      ⋅  467665
      ⋅  115451
 249545       ⋅
In [37]:
sum(W, dims=1)
Out[37]:
1×2 SparseMatrixCSC{Int64, Int64} with 2 stored entries:
 1343423  1548879

Capacity Constrained¶

In [38]:
s = [1_500_000, 1_500_000]
d = w

alloc(X) = trans(Dgc(X, P), s, d)[2]

W = alloc(X)'

Presolve 2 (-10) rows, 10 (-10) columns and 20 (-20) elements
0  Obj 2.1079484e+08 Primal inf 1392302 (1) Dual inf 257.33497 (4)
1  Obj 80752829
Optimal - objective value 80752829
After Postsolve, objective 80752829, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 80752828.73 - 1 iterations time 0.002, Presolve 0.00
Out[38]:
10×2 adjoint(::Matrix{Float64}) with eltype Float64:
      0.0  174721.0
 874579.0       0.0
 105240.0       0.0
      0.0  283506.0
      0.0  208501.0
  48879.0  250156.0
 114059.0       0.0
      0.0  467665.0
      0.0  115451.0
 249545.0       0.0