Net Mod 1: Linear Assignment and Transportation Problems

Packages Used: No new packages are used in this notebook.

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 1: 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.

C = [4 8 5; 2 1 9; 6 3 7]     # Car-to-passenger arc cost (distance) matrix

3×3 Matrix{Int64}:
 4  8  5
 2  1  9
 6  3  7

idx = argmin(C)               # Find closest passenger

CartesianIndex(2, 2)

C[idx]

1

C[idx[1],:] .= maximum(C)+1   # Assign large value so that assignment not made again
C[:,idx[2]] .= maximum(C)+1
C

3×3 Matrix{Int64}:
  4  11   5
 10  11  10
  6  11   7

idx = argmin(C)               # Find next closest passenger

CartesianIndex(1, 1)

C[idx]

4

C[idx[1],:] .= maximum(C)+1
C[:,idx[2]] .= maximum(C)+1
C

3×3 Matrix{Int64}:
 13  12  12
 13  11  10
 13  11   7

idx = argmin(C)

CartesianIndex(3, 3)

C[idx]

7

TC = 1 + 4 + 7                # Final cost of assignment

12

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

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))
    for i = 1:size(C,2)
        idx = argmin(C)
        X[idx] = 1
        C[idx[1],:] .= maximum(C)+1
        C[:,idx[2]] .= maximum(C)+1
    end
    return X
end
X = greedyassignment(C)

3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

TC = sum(sum(C.*X))

12.0

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

greedyassignment(C[1:2,:])   # Not feasible to have one car serve multiple passengers

2×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  1.0

greedyassignment(C[:,1:2])   # Feasible to have extra car not assigned

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 Cars} \sum_{j \in Pass} c_{ij} x_{i,j} \\ \quad \mbox{subject to} \quad \sum_{j \in Pass} x_{ij} &\le& 1\mbox{,} &\quad i \in Cars\mbox{,}\quad\mbox{: each car is assigned up to one passenger}\\ \sum_{i \in Cars} x_{ij} &=& 1\mbox{,} &\quad j \in Pass\mbox{,}\quad\mbox{: each passenger is assigned a car}\\ x_i &\ge& 0 \mbox{,} \end{eqnarray*} $

using JuMP, GLPK
Cars = 1:size(C,1)
Pass = 1:size(C,2)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= X[1:length(Cars),1:length(Pass)] )
@objective(model, Min, sum(C[i,j]*X[i,j] for i ∈ Cars, j ∈ Pass) )
for i ∈ Cars
    @constraint(model, sum(X[i,j] for j ∈ Pass) <= 1 )
end
for j ∈ Pass
    @constraint(model, sum(X[i,j] for i ∈ Cars) == 1 )
end
print(model)
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
Xᵒ = value.(X)
TCᵒ, 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.0\\ & X_{1,2} + X_{2,2} + X_{3,2} = 1.0\\ & X_{1,3} + X_{2,3} + X_{3,3} = 1.0\\ & X_{1,1} + X_{1,2} + X_{1,3} \leq 1.0\\ & X_{2,1} + X_{2,2} + X_{2,3} \leq 1.0\\ & X_{3,1} + X_{3,2} + X_{3,3} \leq 1.0\\ & X_{1,1} \geq 0.0\\ & X_{2,1} \geq 0.0\\ & X_{3,1} \geq 0.0\\ & X_{1,2} \geq 0.0\\ & X_{2,2} \geq 0.0\\ & X_{3,2} \geq 0.0\\ & X_{1,3} \geq 0.0\\ & X_{2,3} \geq 0.0\\ & X_{3,3} \geq 0.0\\ \end{aligned} $$

GLPK Simplex Optimizer 5.0
6 rows, 9 columns, 18 non-zeros
      0: obj =  0.000000000e+000 inf =  3.000e+000 (3)
      5: obj =  1.200000000e+001 inf =  0.000e+000 (0)
*     8: obj =  1.000000000e+001 inf =  0.000e+000 (0)
OPTIMAL LP SOLUTION FOUND

(10.0, [0.0 0.0 1.0; 1.0 0.0 0.0; 0.0 1.0 0.0])

Xᵒ

3×3 Matrix{Float64}:
 0.0  0.0  1.0
 1.0  0.0  0.0
 0.0  1.0  0.0

Ex 2: 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.

Side Issue: Combining Vector Outputs into a Matrix

myvecfun(i) = [i:i+3;]'                                      # Returns a row vector
[myvecfun(i) for i = 1:4:12]

3-element Vector{LinearAlgebra.Adjoint{Int64, Vector{Int64}}}:
 [1 2 3 4]
 [5 6 7 8]
 [9 10 11 12]

Using the splat operator ... to pass each element as a separate input to the vertical concatenation function vcat converts the vectors into a matrix:

vcat([myvecfun(i) for i = 1:4:12]...)                        # Returns a matrix

3×4 Matrix{Int64}:
 1   2   3   4
 5   6   7   8
 9  10  11  12

Distance Matrix

XY1 = [1 1; 4 4; 8 8]                                        # Car cooridinates
XY2 = [2 6; 5 2; 4 9]                                        # Passenger cooridinates
dist2(x,P) = sqrt.(sum((x' .- P).^2, dims = 2))              # Returns a row vector
D = vcat([dist2(XY1[i,:],XY2) for i in 1:size(XY1,1)]'...)   # Returns a matrix

3×3 Matrix{Float64}:
 5.09902  4.12311  8.544
 2.82843  2.23607  5.0
 6.32456  6.7082   4.12311

greedyassignment(D)

3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

sum(sum(D.*greedyassignment(D)))

11.458193116710234

using JuMP, GLPK
function LAP(C)
    N₁ = 1:size(C,1)
    N₂ = 1:size(C,2)
    model = Model(GLPK.Optimizer)
    @variable(model, 0 <= X[1:length(N₁), 1:length(N₂)] )
    @objective(model, Min, sum(C[i,j]*X[i,j] for i ∈ N₁, j ∈ N₂) )
    for i ∈ N₁
        @constraint(model, sum(X[i,j] for j ∈ N₂) <= 1 )
    end
    for j ∈ N₂
        @constraint(model, sum(X[i,j] for i ∈ N₁) == 1 )
    end
    set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
    optimize!(model)
    return value.(X)
end
Xᵒ = LAP(D)
TCᵒ = sum(sum(D.*Xᵒ))
TCᵒ, Xᵒ

GLPK Simplex Optimizer 5.0
6 rows, 9 columns, 18 non-zeros
      0: obj =  0.000000000e+000 inf =  3.000e+000 (3)
      5: obj =  1.145819312e+001 inf =  0.000e+000 (0)
*     8: obj =  1.107463838e+001 inf =  0.000e+000 (0)
OPTIMAL LP SOLUTION FOUND

(11.074638375981511, [0.0 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 1.0])

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

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

GLPK Simplex Optimizer 5.0
5 rows, 6 columns, 12 non-zeros
      0: obj =  0.000000000e+000 inf =  3.000e+000 (3)
      4: obj =  7.335087491e+000 inf =  1.000e+000 (1)
LP HAS NO PRIMAL FEASIBLE SOLUTION

2×3 Matrix{Float64}:
 0.0   2.0  0.0
 1.0  -1.0  1.0

LAP(D[:,1:2])   # Feasible to have extra car not assigned

GLPK Simplex Optimizer 5.0
5 rows, 6 columns, 12 non-zeros
      0: obj =  0.000000000e+000 inf =  2.000e+000 (2)
      3: obj =  7.335087491e+000 inf =  0.000e+000 (0)
*     6: obj =  6.951532750e+000 inf =  0.000e+000 (0)
OPTIMAL LP SOLUTION FOUND

3×2 Matrix{Float64}:
 0.0  1.0
 1.0  0.0
 0.0  0.0

Ex 3: Ten Random Cars and Passengers

using Random
Random.seed!(1215)
n = 10
XY1 = 10*rand(n,2)
XY2 = 10*rand(n,2)
D = vcat([dist2(XY1[i,:],XY2) for i in 1:size(XY1,1)]'...)

10×10 Matrix{Float64}:
 9.58611  4.51146   8.18951  6.09541   …  8.33746  4.08901   8.08572
 3.88209  1.88679   5.384    0.246301     4.24533  5.9087    6.42391
 3.36016  1.80746   7.13067  1.50139      5.96115  7.09449   8.13832
 6.51161  1.80711   7.86873  3.57492      7.31778  5.67638   8.37062
 8.11201  4.14016   3.8528   4.4492       4.29152  1.7083    3.80315
 9.5443   4.46973   7.75516  5.96929   …  7.96673  3.58464   7.61078
 1.33956  4.63668   6.84334  3.02714      5.21725  8.8013    8.21285
 7.047    4.3019   10.7824   5.64745      9.98565  8.75441  11.4185
 9.17575  7.48417   1.05827  6.681        2.65182  5.21653   1.02802
 8.48535  3.53646   6.30472  4.78033      6.46701  2.64633   6.27466

X = greedyassignment(D)

10×10 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  1.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  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  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  1.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  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  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0

TC = sum(sum(D.*X))

27.81613198823573

using Plots
idx = findall(X .> 0)
x = [[XY1[idx[i][1],1]; XY2[idx[i][2],1]] for i = 1:size(XY1,1)]
y = [[XY1[idx[i][1],2]; XY2[idx[i][2],2]] for i = 1:size(XY1,1)]
plot(x, y, label=false)
scatter!([XY1[:,1]],[XY1[:,2]],label="Cars")
scatter!([XY2[:,1]],[XY2[:,2]],label="Passengers", legend=:bottomright)

Xᵒ = LAP(D)

GLPK Simplex Optimizer 5.0
20 rows, 100 columns, 200 non-zeros
      0: obj =  0.000000000e+000 inf =  1.000e+001 (10)
     19: obj =  5.648289937e+001 inf =  0.000e+000 (0)
*    40: obj =  2.530998090e+001 inf =  0.000e+000 (0)
OPTIMAL LP SOLUTION FOUND

10×10 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  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  1.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  0.0  0.0  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  0.0  0.0  0.0  1.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  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  1.0  0.0  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

TCᵒ = sum(sum(D.*Xᵒ))

25.309980899746733

idx = findall(Xᵒ .> 0)
x = [[XY1[idx[i][1],1]; XY2[idx[i][2],1]] for i = 1:size(XY1,1)]
y = [[XY1[idx[i][1],2]; XY2[idx[i][2],2]] for i = 1:size(XY1,1)]
plot(x, y, label=false)
scatter!([XY1[:,1]],[XY1[:,2]],label="Cars")
scatter!([XY2[:,1]],[XY2[:,2]],label="Passengers", legend=:bottomright)

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 4: Transportation Problem with Excess Supply

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)

20

using JuMP, GLPK
function trans0(C,s,d)
    S = 1:length(s)
    D = 1:length(d)
    model = Model(GLPK.Optimizer)
    @variable(model, 0 <= X[1:length(S), 1:length(D)] )
    @objective(model, Min, sum(C[i,j]*X[i,j] for i ∈ S, j ∈ D) )
    for i ∈ S
        @constraint(model, sum(X[i,j] for j ∈ D) <= s[i] )
    end
    for j ∈ D
        @constraint(model, sum(X[i,j] for i ∈ S) == d[j] )
    end
    set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
    optimize!(model)
    return value.(X)
end
Xᵒ = trans0(C,s,d)
@show TCᵒ = sum(sum(C.*Xᵒ))
Xᵒ

GLPK Simplex Optimizer 5.0
7 rows, 12 columns, 24 non-zeros
      0: obj =  0.000000000e+000 inf =  1.250e+002 (4)
      6: obj =  1.100000000e+003 inf =  0.000e+000 (0)
*     9: obj =  9.700000000e+002 inf =  0.000e+000 (0)
OPTIMAL LP SOLUTION FOUND
TCᵒ = sum(sum(C .* Xᵒ)) = 970.0

3×4 Matrix{Float64}:
  5.0  20.0  30.0   0.0
 40.0   0.0   0.0   0.0
  0.0   0.0   0.0  30.0

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

s = [25, 50, 40]
@show sum(s) - sum(d)
trans0(C,s,d)

sum(s) - sum(d) = -10
GLPK Simplex Optimizer 5.0
7 rows, 12 columns, 24 non-zeros
      0: obj =  0.000000000e+000 inf =  1.250e+002 (4)
      6: obj =  1.170000000e+003 inf =  1.000e+001 (1)
LP HAS NO PRIMAL FEASIBLE SOLUTION

3×4 Matrix{Float64}:
 -5.0   0.0  30.0   0.0
 50.0   0.0   0.0   0.0
  0.0  20.0   0.0  30.0

Can add a check inside trans to detect supply being less than demand so that an error is thrown and nothing returned instead of having to check the solver message to determine that an erroneous solution is being returned when using trans0:

using JuMP, GLPK
function trans(C,s,d)
    if sum(s) < sum(d)
        error("Infeasible: Supply less than demand")
    end
    S = 1:length(s)
    D = 1:length(d)
    model = Model(GLPK.Optimizer)
    @variable(model, 0 <= X[1:length(S), 1:length(D)] )
    @objective(model, Min, sum(C[i,j]*X[i,j] for i ∈ S, j ∈ D) )
    for i ∈ S
        @constraint(model, sum(X[i,j] for j ∈ D) <= s[i] )
    end
    for j ∈ D
        @constraint(model, sum(X[i,j] for i ∈ S) == d[j] )
    end
    set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
    optimize!(model)
    return value.(X)
end
Xᵒ = trans(C,s,d)

Infeasible: Supply less than demand

Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:33
 [2] trans(C::Matrix{Int64}, s::Vector{Int64}, d::Vector{Int64})
   @ Main .\In[35]:4
 [3] top-level scope
   @ In[35]:21
 [4] eval
   @ .\boot.jl:360 [inlined]
 [5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base .\loading.jl:1116

Ex 5: 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.

s = [180, 300, 250]
d = [130, 82, 215, 175]
XY1 = [2 2; 1 1; 2 4]
XY2 = [3 2; 2 3; 4 2; 3 4]
dist2(x,P) = sqrt.(sum((x' .- P).^2, dims = 2))
D = vcat([dist2(XY1[i,:],XY2) for i in 1:size(XY1,1)]'...)

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

Xᵒ = trans(D,s,d)
TCᵒ = sum(sum(D.*Xᵒ))

GLPK Simplex Optimizer 5.0
7 rows, 12 columns, 24 non-zeros
      0: obj =  0.000000000e+000 inf =  6.020e+002 (4)
      6: obj =  1.244538090e+003 inf =  0.000e+000 (0)
*     9: obj =  1.016911758e+003 inf =  0.000e+000 (0)
OPTIMAL LP SOLUTION FOUND

1016.9117575489613

Matrix([Xᵒ sum(Xᵒ,dims=2) s; sum(Xᵒ,dims=1) 0 0; d' 0 0 ])

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