Ex 1: Unter's Passenger-to-Car Assignment (Ex 1 in Net Mod 1)¶
Packages Used: Functions from the following packages are used in this notebook for the first time:
Combinatorial optimization refers to procedures used to find solutions from a collection of finite objects that characterize possible solution instances. These procedures fall into three broad categories:
Brute-Force Optimization
Heuristics and optimal procedures efficiently look at a small fraction of the possible solution instances to find a locally or globally optimal solution, but it is possible in many cases if the number of solutions is not too big to brute-force generate all possible solutions to find the global optimum. The ability to generate solution instances is useful for many reasons:
The Combinatorics package has several functions that generate combinatorial objects that can be used to represent instances of many different types of problems:
using Combinatorics
@show names(Combinatorics);
names(Combinatorics) = [:Combinatorics, :CoolLexCombinations, :Partition, :SkewDiagram, :bellnum, :catalannum, :character, :combinations, :derangement, :doublefactorial, :factorial, :fibonaccinum, :gamma, :hyperfactorial, :integer_partitions, :isrimhook, :jacobisymbol, :lassallenum, :legendresymbol, :leglength, :levicivita, :lobbnum, :lucasnum, :multiexponents, :multifactorial, :multinomial, :multiset_combinations, :multiset_permutations, :narayana, :ncpartitions, :nthperm, :nthperm!, :parity, :partialderangement, :partitions, :partitionsequence, :permutations, :powerset, :prevprod, :primorial, :stirlings1, :stirlings2, :subfactorial, :with_replacement_combinations]
A permutation is listing of the indices from 1 to n where the order does matter. The number of permutations is n factorial:
$ \quad n! = n \cdot (n - 1) \cdot (n - 2) \dotsb 2 \cdot 1 \implies 3! = 6 $
n = 3
factorial(n)
6
Need to use collect to generate all permutations because permutations only iterates one permutation at a time:
using Combinatorics
collect(permutations(1:n))
6-element Vector{Vector{Int64}}:
[1, 2, 3]
[1, 3, 2]
[2, 1, 3]
[2, 3, 1]
[3, 1, 2]
[3, 2, 1]
The function randperm returns a random ordering the indices 1 to n:
using Random
Random.seed!(234)
randperm(n)
3-element Vector{Int64}:
2
1
3
randperm(n)
3-element Vector{Int64}:
1
3
2
Often used to randomly reorder the elements in an n-element array:
v = [1,13,5,27,9]
v[randperm(length(v))]
5-element Vector{Int64}:
5
13
1
27
9
The function shuffle(v) is equavilent to v[randperm(length(v))], which is both more concise and does not require that the array v first be created:
shuffle(v)
5-element Vector{Int64}:
5
1
9
27
13
The number of permutations, $n!$, for $n = 12$ and $n = 50$:
factorial(12)
479001600
log10(factorial(12)) # Number of digits
8.680336964082644
factorial(big(50))
30414093201713378043612608166064768844377641568960512000000000000
log10(factorial(big(50))) # Number of digits
64.48307487247203535631480253445628010182653214430050816125726166556077150092167
n = 3
C = [4 8 5; 2 1 9; 6 3 7]
Xᵒ = [0 0 1; 1 0 0; 0 1 0] # Optimal solution
TCᵒ = sum(sum(C.*Xᵒ))
10
Xᵒ # Xᵒ termed a permutation matrix
3×3 Matrix{Int64}:
0 0 1
1 0 0
0 1 0
α = [3,1,2] # α termed a permutation vector
[(i,j) for (i,j) in enumerate(α)]
3-element Vector{Tuple{Int64, Int64}}:
(1, 3)
(2, 1)
(3, 2)
sum([C[i,j] for (i,j) in enumerate(α)])
10
TCᵒ, αᵒ = Inf, []
for α in permutations(1:n)
TC = sum([C[i,j] for (i,j) in enumerate(α)])
if TC < TCᵒ
TCᵒ, αᵒ = TC, α
end
println(α,": ",TC)
end
[1, 2, 3]: 12 [1, 3, 2]: 16 [2, 1, 3]: 17 [2, 3, 1]: 23 [3, 1, 2]: 10 [3, 2, 1]: 12
A combination is a selection of items from a collection where the order does not matter. The number of combinations of n items taken k at a time without repetition (also referred to as "n choose k" and the binomial coefficient) is:
$ \quad \dbinom{n}{k} = \dfrac{n!}{(n - k)!\,k!} \implies \dbinom{4}{3} = 4 $
n = 4
binomial(n,3)
4
collect(combinations(1:n,3))
4-element Vector{Vector{Int64}}:
[1, 2, 3]
[1, 2, 4]
[1, 3, 4]
[2, 3, 4]
Given an array of length $n$, the number of pairwise interchanges of its elements is
$ \quad \dbinom{n}{2} = \dfrac{n(n - 1)}{2} \implies \dbinom{4}{2} = 6 $
When the array represents a permutation, it is often computationally infeasible to generate all of the permutations; instead, just looking at the pairwise interchanges of the array allows it to be altered in an efficient manner.
collect(combinations(1:n,2))
6-element Vector{Vector{Int64}}:
[1, 2]
[1, 3]
[1, 4]
[2, 3]
[2, 4]
[3, 4]
for i = 1:n-1, j = i+1:n
println([i,j])
end
[1, 2] [1, 3] [1, 4] [2, 3] [2, 4] [3, 4]
collect(permutations(1:n))
24-element Vector{Vector{Int64}}:
[1, 2, 3, 4]
[1, 2, 4, 3]
[1, 3, 2, 4]
[1, 3, 4, 2]
[1, 4, 2, 3]
[1, 4, 3, 2]
[2, 1, 3, 4]
[2, 1, 4, 3]
[2, 3, 1, 4]
[2, 3, 4, 1]
[2, 4, 1, 3]
[2, 4, 3, 1]
[3, 1, 2, 4]
[3, 1, 4, 2]
[3, 2, 1, 4]
[3, 2, 4, 1]
[3, 4, 1, 2]
[3, 4, 2, 1]
[4, 1, 2, 3]
[4, 1, 3, 2]
[4, 2, 1, 3]
[4, 2, 3, 1]
[4, 3, 1, 2]
[4, 3, 2, 1]
using Random
Random.seed!(5223)
αᵒ = randperm(n)
println(" : ",αᵒ)
for (i,j) in combinations(1:n,2)
α = copy(αᵒ) # Need to make copy so αᵒ does not get changed
α[[j,i]] .= α[[i,j]]
println(i,",",j,": ",α)
end
: [3, 4, 2, 1] 1,2: [4, 3, 2, 1] 1,3: [2, 4, 3, 1] 1,4: [1, 4, 2, 3] 2,3: [3, 2, 4, 1] 2,4: [3, 1, 2, 4] 3,4: [3, 4, 1, 2]
factorial(12),binomial(12,2)
(479001600, 66)
factorial(big(50)),binomial(50,2)
(30414093201713378043612608166064768844377641568960512000000000000, 1225)
using LightGraphs, SimpleWeightedGraphs
i = [1, 1, 2, 2, 3, 3, 4, 4, 5] # Index of starting node of each arc
j = [2, 3, 3, 4, 4, 5, 5, 6, 6] # Index of ending node of each arc
w = [4.,2.,1.,5.,8.,10.,2.,6.,3.] # Weight (e.g., cost, distance, time) of each arc
G = SimpleWeightedGraph(i, j, w) # Create undirected network
ds = dijkstra_shortest_paths(G, 1) # Shortest paths from node 1 to other all nodes
@show ds.dists[6] # Distance from node 1 to 6
enumerate_paths(ds, 6) # Path from node 1 to node 6
ds.dists[6] = 13.0
6-element Vector{Int64}:
1
3
2
4
5
6
The number of possible paths in an undirected network is:
$ \begin{eqnarray} \quad \sum_{k=0}^{n-2} \dbinom{n-2}{k} k! &\implies& \sum_{k=0}^{6-2} \dbinom{6-2}{k} k! = 65, \mbox{ for } n = 6 \end{eqnarray}$
n = 6
factorial(n), sum([binomial(n-2,k)*factorial(k) for k=0:(n-2)])
(720, 65)
cnt = 0
for k = 0:n-2
for i in combinations(1:n-2,k) # Don't include nodes 1 and 6
for j in permutations(i)
cnt += 1
p = hcat(1, j.+1..., 6) # Add 1 so j starts at 2
println(cnt,": ",p)
end
end
end
1: [1 6] 2: [1 2 6] 3: [1 3 6] 4: [1 4 6] 5: [1 5 6] 6: [1 2 3 6] 7: [1 3 2 6] 8: [1 2 4 6] 9: [1 4 2 6] 10: [1 2 5 6] 11: [1 5 2 6] 12: [1 3 4 6] 13: [1 4 3 6] 14: [1 3 5 6] 15: [1 5 3 6] 16: [1 4 5 6] 17: [1 5 4 6] 18: [1 2 3 4 6] 19: [1 2 4 3 6] 20: [1 3 2 4 6] 21: [1 3 4 2 6] 22: [1 4 2 3 6] 23: [1 4 3 2 6] 24: [1 2 3 5 6] 25: [1 2 5 3 6] 26: [1 3 2 5 6] 27: [1 3 5 2 6] 28: [1 5 2 3 6] 29: [1 5 3 2 6] 30: [1 2 4 5 6] 31: [1 2 5 4 6] 32: [1 4 2 5 6] 33: [1 4 5 2 6] 34: [1 5 2 4 6] 35: [1 5 4 2 6] 36: [1 3 4 5 6] 37: [1 3 5 4 6] 38: [1 4 3 5 6] 39: [1 4 5 3 6] 40: [1 5 3 4 6] 41: [1 5 4 3 6] 42: [1 2 3 4 5 6] 43: [1 2 3 5 4 6] 44: [1 2 4 3 5 6] 45: [1 2 4 5 3 6] 46: [1 2 5 3 4 6] 47: [1 2 5 4 3 6] 48: [1 3 2 4 5 6] 49: [1 3 2 5 4 6] 50: [1 3 4 2 5 6] 51: [1 3 4 5 2 6] 52: [1 3 5 2 4 6] 53: [1 3 5 4 2 6] 54: [1 4 2 3 5 6] 55: [1 4 2 5 3 6] 56: [1 4 3 2 5 6] 57: [1 4 3 5 2 6] 58: [1 4 5 2 3 6] 59: [1 4 5 3 2 6] 60: [1 5 2 3 4 6] 61: [1 5 2 4 3 6] 62: [1 5 3 2 4 6] 63: [1 5 3 4 2 6] 64: [1 5 4 2 3 6] 65: [1 5 4 3 2 6]
W = adjacency_matrix(G)
Matrix(W)
6×6 Matrix{Float64}:
0.0 4.0 2.0 0.0 0.0 0.0
4.0 0.0 1.0 5.0 0.0 0.0
2.0 1.0 0.0 8.0 10.0 0.0
0.0 5.0 8.0 0.0 2.0 6.0
0.0 0.0 10.0 2.0 0.0 3.0
0.0 0.0 0.0 6.0 3.0 0.0
W = Matrix(replace!(W, 0. => Inf))
6×6 Matrix{Float64}:
Inf 4.0 2.0 Inf Inf Inf
4.0 Inf 1.0 5.0 Inf Inf
2.0 1.0 Inf 8.0 10.0 Inf
Inf 5.0 8.0 Inf 2.0 6.0
Inf Inf 10.0 2.0 Inf 3.0
Inf Inf Inf 6.0 3.0 Inf
p = enumerate_paths(ds, 6)
6-element Vector{Int64}:
1
3
2
4
5
6
[(i,j) for (i,j) in zip(p[1:end-1],p[2:end])]
5-element Vector{Tuple{Int64, Int64}}:
(1, 3)
(3, 2)
(2, 4)
(4, 5)
(5, 6)
sum([W[i,j] for (i,j) in zip(p[1:end-1],p[2:end])])
13.0
TCᵒ, pᵒ, cnt = Inf, 0, 0
for k = 0:n-2
for i in combinations(1:n-2,k)
for j in permutations(i)
p = hcat(1, j.+1..., 6)
TC = sum([W[i,j] for (i,j) in zip(p[1:end-1],p[2:end])])
if TC < TCᵒ
TCᵒ, pᵒ = TC, p
end
if ~isinf(TC)
cnt += 1
println(cnt,": ",p," = ",TC)
end
end
end
end
pᵒ, TCᵒ
1: [1 2 4 6] = 15.0 2: [1 3 4 6] = 16.0 3: [1 3 5 6] = 15.0 4: [1 2 3 4 6] = 19.0 5: [1 3 2 4 6] = 14.0 6: [1 2 3 5 6] = 18.0 7: [1 2 4 5 6] = 14.0 8: [1 3 4 5 6] = 15.0 9: [1 3 5 4 6] = 20.0 10: [1 2 3 4 5 6] = 18.0 11: [1 2 3 5 4 6] = 23.0 12: [1 2 4 3 5 6] = 30.0 13: [1 3 2 4 5 6] = 13.0
([1 3 … 5 6], 13.0)
@doc powerset
powerset(a, min=0, max=length(a))Generate all subsets of an indexable object a including the empty set, with cardinality bounded by min and max. Because the number of subsets can be very large, this function returns an iterator object. Use collect(powerset(a, min, max)) to get an array of all subsets.
n = 3
2^n
8
collect(powerset(1:n))
8-element Vector{Vector{Int64}}:
[]
[1]
[2]
[3]
[1, 2]
[1, 3]
[2, 3]
[1, 2, 3]
collect(powerset(1:n,1,n)) # 2^n - 1 non-empty subsets
7-element Vector{Vector{Int64}}:
[1]
[2]
[3]
[1, 2]
[1, 3]
[2, 3]
[1, 2, 3]
collect(powerset(1:n,2,2))
3-element Vector{Vector{Int64}}:
[1, 2]
[1, 3]
[2, 3]
[sum(i) for i in powerset(1:n,2,2)] # Can iterate over subsets
3-element Vector{Int64}:
3
4
5
Instead of solving Ex 1 in MIP 2 as a MILP, all non-empty subsets of the investment opportunities could be considered the fine the subset that maximizes profit while statisfying all of the constraints:
revenue = [9, 12, 7, 5, 2]
cost = [5, 7, 4, 3, 1]
cmax = 14
m = length(revenue)
collect(powerset(1:m,1,m))
31-element Vector{Vector{Int64}}:
[1]
[2]
[3]
[4]
[5]
[1, 2]
[1, 3]
[1, 4]
[1, 5]
[2, 3]
[2, 4]
[2, 5]
[3, 4]
⋮
[1, 3, 5]
[1, 4, 5]
[2, 3, 4]
[2, 3, 5]
[2, 4, 5]
[3, 4, 5]
[1, 2, 3, 4]
[1, 2, 3, 5]
[1, 2, 4, 5]
[1, 3, 4, 5]
[2, 3, 4, 5]
[1, 2, 3, 4, 5]
for M in powerset(1:m,1,m)
x = zeros(m)
x[M] .= 1
println(x)
end
[1.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, 1.0, 0.0] [0.0, 0.0, 0.0, 0.0, 1.0] [1.0, 1.0, 0.0, 0.0, 0.0] [1.0, 0.0, 1.0, 0.0, 0.0] [1.0, 0.0, 0.0, 1.0, 0.0] [1.0, 0.0, 0.0, 0.0, 1.0] [0.0, 1.0, 1.0, 0.0, 0.0] [0.0, 1.0, 0.0, 1.0, 0.0] [0.0, 1.0, 0.0, 0.0, 1.0] [0.0, 0.0, 1.0, 1.0, 0.0] [0.0, 0.0, 1.0, 0.0, 1.0] [0.0, 0.0, 0.0, 1.0, 1.0] [1.0, 1.0, 1.0, 0.0, 0.0] [1.0, 1.0, 0.0, 1.0, 0.0] [1.0, 1.0, 0.0, 0.0, 1.0] [1.0, 0.0, 1.0, 1.0, 0.0] [1.0, 0.0, 1.0, 0.0, 1.0] [1.0, 0.0, 0.0, 1.0, 1.0] [0.0, 1.0, 1.0, 1.0, 0.0] [0.0, 1.0, 1.0, 0.0, 1.0] [0.0, 1.0, 0.0, 1.0, 1.0] [0.0, 0.0, 1.0, 1.0, 1.0] [1.0, 1.0, 1.0, 1.0, 0.0] [1.0, 1.0, 1.0, 0.0, 1.0] [1.0, 1.0, 0.0, 1.0, 1.0] [1.0, 0.0, 1.0, 1.0, 1.0] [0.0, 1.0, 1.0, 1.0, 1.0] [1.0, 1.0, 1.0, 1.0, 1.0]
TPᵒ, xᵒ, cnt = -Inf, [], 0
for M in powerset(1:m,1,m)
x = zeros(m)
x[M] .= 1
TP = sum(revenue[M] .- cost[M])
if TP > TPᵒ &&
sum([cost[i]*x[i] for i ∈ M]) <= cmax # Budget constraint
TPᵒ, xᵒ = TP, x
end
cnt += 1
println(cnt,": ",TPᵒ,", ",TP,", ",sum(cost[M]),", ", x)
end
TPᵒ, xᵒ
1: 4, 4, 5, [1.0, 0.0, 0.0, 0.0, 0.0] 2: 5, 5, 7, [0.0, 1.0, 0.0, 0.0, 0.0] 3: 5, 3, 4, [0.0, 0.0, 1.0, 0.0, 0.0] 4: 5, 2, 3, [0.0, 0.0, 0.0, 1.0, 0.0] 5: 5, 1, 1, [0.0, 0.0, 0.0, 0.0, 1.0] 6: 9, 9, 12, [1.0, 1.0, 0.0, 0.0, 0.0] 7: 9, 7, 9, [1.0, 0.0, 1.0, 0.0, 0.0] 8: 9, 6, 8, [1.0, 0.0, 0.0, 1.0, 0.0] 9: 9, 5, 6, [1.0, 0.0, 0.0, 0.0, 1.0] 10: 9, 8, 11, [0.0, 1.0, 1.0, 0.0, 0.0] 11: 9, 7, 10, [0.0, 1.0, 0.0, 1.0, 0.0] 12: 9, 6, 8, [0.0, 1.0, 0.0, 0.0, 1.0] 13: 9, 5, 7, [0.0, 0.0, 1.0, 1.0, 0.0] 14: 9, 4, 5, [0.0, 0.0, 1.0, 0.0, 1.0] 15: 9, 3, 4, [0.0, 0.0, 0.0, 1.0, 1.0] 16: 9, 12, 16, [1.0, 1.0, 1.0, 0.0, 0.0] 17: 9, 11, 15, [1.0, 1.0, 0.0, 1.0, 0.0] 18: 10, 10, 13, [1.0, 1.0, 0.0, 0.0, 1.0] 19: 10, 9, 12, [1.0, 0.0, 1.0, 1.0, 0.0] 20: 10, 8, 10, [1.0, 0.0, 1.0, 0.0, 1.0] 21: 10, 7, 9, [1.0, 0.0, 0.0, 1.0, 1.0] 22: 10, 10, 14, [0.0, 1.0, 1.0, 1.0, 0.0] 23: 10, 9, 12, [0.0, 1.0, 1.0, 0.0, 1.0] 24: 10, 8, 11, [0.0, 1.0, 0.0, 1.0, 1.0] 25: 10, 6, 8, [0.0, 0.0, 1.0, 1.0, 1.0] 26: 10, 14, 19, [1.0, 1.0, 1.0, 1.0, 0.0] 27: 10, 13, 17, [1.0, 1.0, 1.0, 0.0, 1.0] 28: 10, 12, 16, [1.0, 1.0, 0.0, 1.0, 1.0] 29: 10, 10, 13, [1.0, 0.0, 1.0, 1.0, 1.0] 30: 10, 11, 15, [0.0, 1.0, 1.0, 1.0, 1.0] 31: 10, 15, 20, [1.0, 1.0, 1.0, 1.0, 1.0]
(10, [1.0, 1.0, 0.0, 0.0, 1.0])
TPᵒ, xᵒ, cnt = -Inf, [], 0
for M in powerset(1:m,1,m)
x = zeros(m)
x[M] .= 1
TP = sum(revenue[M] .- cost[M])
if TP > TPᵒ &&
sum(cost[i]*x[i] for i ∈ M) <= cmax && # Budget constraint
sum(x[i] for i ∈ M) <= 2 && # 1. Can only make two investments
x[2] <= x[5] && # 2. If 2 is made, 5 must also be made
x[1] + x[3] == 1 && # 3. Both 1 and 3 cannot be made
x[5] >= x[1] && x[5] >= x[4] && # 4. If 1 or 4 is made then 5 must be made
x[2] + x[3] + x[4] >= 1 # 5. At least one of 2,3, or 4 must be made
TPᵒ, xᵒ = TP, x
end
cnt += 1
println(cnt,": ",TPᵒ,", ",TP,", ",sum(cost[M]),", ", x)
end
TPᵒ, xᵒ
1: -Inf, 4, 5, [1.0, 0.0, 0.0, 0.0, 0.0] 2: -Inf, 5, 7, [0.0, 1.0, 0.0, 0.0, 0.0] 3: 3, 3, 4, [0.0, 0.0, 1.0, 0.0, 0.0] 4: 3, 2, 3, [0.0, 0.0, 0.0, 1.0, 0.0] 5: 3, 1, 1, [0.0, 0.0, 0.0, 0.0, 1.0] 6: 3, 9, 12, [1.0, 1.0, 0.0, 0.0, 0.0] 7: 3, 7, 9, [1.0, 0.0, 1.0, 0.0, 0.0] 8: 3, 6, 8, [1.0, 0.0, 0.0, 1.0, 0.0] 9: 3, 5, 6, [1.0, 0.0, 0.0, 0.0, 1.0] 10: 3, 8, 11, [0.0, 1.0, 1.0, 0.0, 0.0] 11: 3, 7, 10, [0.0, 1.0, 0.0, 1.0, 0.0] 12: 3, 6, 8, [0.0, 1.0, 0.0, 0.0, 1.0] 13: 3, 5, 7, [0.0, 0.0, 1.0, 1.0, 0.0] 14: 4, 4, 5, [0.0, 0.0, 1.0, 0.0, 1.0] 15: 4, 3, 4, [0.0, 0.0, 0.0, 1.0, 1.0] 16: 4, 12, 16, [1.0, 1.0, 1.0, 0.0, 0.0] 17: 4, 11, 15, [1.0, 1.0, 0.0, 1.0, 0.0] 18: 4, 10, 13, [1.0, 1.0, 0.0, 0.0, 1.0] 19: 4, 9, 12, [1.0, 0.0, 1.0, 1.0, 0.0] 20: 4, 8, 10, [1.0, 0.0, 1.0, 0.0, 1.0] 21: 4, 7, 9, [1.0, 0.0, 0.0, 1.0, 1.0] 22: 4, 10, 14, [0.0, 1.0, 1.0, 1.0, 0.0] 23: 4, 9, 12, [0.0, 1.0, 1.0, 0.0, 1.0] 24: 4, 8, 11, [0.0, 1.0, 0.0, 1.0, 1.0] 25: 4, 6, 8, [0.0, 0.0, 1.0, 1.0, 1.0] 26: 4, 14, 19, [1.0, 1.0, 1.0, 1.0, 0.0] 27: 4, 13, 17, [1.0, 1.0, 1.0, 0.0, 1.0] 28: 4, 12, 16, [1.0, 1.0, 0.0, 1.0, 1.0] 29: 4, 10, 13, [1.0, 0.0, 1.0, 1.0, 1.0] 30: 4, 11, 15, [0.0, 1.0, 1.0, 1.0, 1.0] 31: 4, 15, 20, [1.0, 1.0, 1.0, 1.0, 1.0]
(4, [0.0, 0.0, 1.0, 0.0, 1.0])