The set covering problem refers to minimizing the number of sets whose union includes all of the objects to be covered. Each set can include (or cover) a different subset of objects, and each object can be covered by multiple sets. There are many applications of set covering; for example, each set might represent the customers that are within a one hour drive of a potential site for a service center, and a solution to the set covering problem would correspond to the minimum number of sites required so that each customer can be serviced within one hour.
General formulation of the problem:
$ \begin{eqnarray*} \quad M &=& \bigl\{ 1, \dots, m \bigr\}, & \quad \mbox{objects to be covered} \\ N &=& \bigl\{ 1, \dots, n \bigr\}, & \quad \mbox{indices of } n \mbox{ subsets of } M \\ M_i &\subseteq& M, i \in N, & \quad n \mbox{ subsets of } M \\ I &\subseteq& N, & \quad \mbox{covering of } M \\ I^* &=& \mathrm{arg}\underset{I}{\operatorname{min}} \bigl\{ \lvert I \rvert : \bigcup_{i \in I} M_i = M \bigr\}, & \quad \mbox{min cost covering of } M. \end{eqnarray*}$
BIP formulation of the problem:
$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{i \in N} x_i \\ \quad \mbox{subject to} \quad \sum_{i \in N} a_{ji} x_i &\ge& 1, &\quad j \in M\\ x_i &\in& \bigl\{ 0, 1 \bigr\}, &\quad i \in N \end{eqnarray*} $
where
$ \begin{eqnarray*} \quad x_i &=& \begin{cases} 1, & \mbox{if } M_i \mbox{ is in cover} \\ 0, & \mbox{otherwise} \end{cases} \\ a_{ji} &=& \begin{cases} 1, & \mbox{if } j \in M_i \\ 0, & \mbox{otherwise.} \end{cases} \\ \end{eqnarray*} $
$ \begin{eqnarray*} \quad M &=& \bigl\{ 1, \dots, 6 \bigr\} \\ N &=& \bigl\{ 1, \dots, 5 \bigr\} \\ M_1 &=& \bigl\{ 1, 2 \bigr\}, M_2 = \bigl\{ 1, 4, 5 \bigr\}, M_3 = \bigl\{ 3, 5 \bigr\}, M_4 = \bigl\{ 2, 3, 6 \bigr\}, M_5 = \bigl\{ 6 \bigr\} \\ I^* &=& \mathrm{arg}\underset{I}{\operatorname{min}} \bigl\{\sum_{i \in I} x_i : \bigcup_{i \in I} M_i = M \bigr\} = \bigl\{ 2, 4 \bigr\} \\ \sum_{i \in I^*} x_i &=& 2 \end{eqnarray*}$
M = 1:6
N = 1:5
Mi = [[1,2],[1,4,5],[3,5],[2,3,6],[6]]
A = zeros(length(M),length(N)) # A = objects x subsets
for i = 1:size(A,2)
A[Mi[i],i] .= 1
end
A
6×5 Matrix{Float64}: 1.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0
using JuMP, Cbc
model = Model(Cbc.Optimizer)
@variable(model, x[1:length(N)], Bin )
@objective(model, Min, sum(x[i] for i ∈ N) )
@constraint(model, [j ∈ M], sum(A[j,i]*x[i] for i ∈ N) >= 1 )
set_optimizer_attribute(model, "logLevel", 1)
optimize!(model)
TCᵒ, xᵒ = objective_value(model), value.(x)
Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Jan 1 1970 command line - Cbc_C_Interface -logLevel 1 -solve -quit (default strategy 1) Continuous objective value is 2 - 0.00 seconds Cgl0004I processed model has 3 rows, 4 columns (4 integer (4 of which binary)) and 6 elements Cutoff increment increased from 1e-05 to 0.9999 Cbc0038I Initial state - 0 integers unsatisfied sum - 0 Cbc0038I Solution found of 2 Cbc0038I Before mini branch and bound, 4 integers at bound fixed and 0 continuous Cbc0038I Mini branch and bound did not improve solution (0.00 seconds) Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of 2 - took 0.00 seconds Cbc0012I Integer solution of 2 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds) Cbc0001I Search completed - best objective 2, took 0 iterations and 0 nodes (0.00 seconds) Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost Cuts at root node changed objective from 2 to 2 Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 2.00000000 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.00 Time (Wallclock seconds): 0.00 Total time (CPU seconds): 0.00 (Wallclock seconds): 0.00
(2.0, [0.0, 1.0, 0.0, 1.0, 0.0])
solution_summary(model).termination_status
OPTIMAL::TerminationStatusCode = 1
idx = findall(value.(x) .> 0)
2-element Vector{Int64}: 2 4
using JuMP, Cbc
function setcover(A)
M = 1:size(A,1)
N = 1:size(A,2)
model = Model(Cbc.Optimizer)
@variable(model, x[1:length(N)], Bin )
@objective(model, Min, sum(x[i] for i ∈ N) )
@constraint(model, [j ∈ M], sum(A[j,i]*x[i] for i ∈ N) >= 1 )
set_optimizer_attribute(model, "logLevel", 0)
set_optimizer_attribute(model, "seconds", 60.0)
optimize!(model)
println(solution_summary(model).termination_status)
return findall(value.(x) .> 0)
end
setcover(A)
OPTIMAL
2-element Vector{Int64}: 2 4
using LightGraphs, SimpleWeightedGraphs
i = [1, 1, 2, 2, 3, 3, 4, 4, 5]
j = [2, 3, 3, 4, 4, 5, 5, 6, 6]
w = [4, 2, 1, 5, 8, 10, 2, 6, 3]
G = SimpleWeightedGraph(i, j, w)
dijkstra_shortest_paths(G,1).dists
6-element Vector{Int64}: 0 3 2 8 10 13
D = vcat([dijkstra_shortest_paths(G,i).dists for i ∈ vertices(G)]'...)
6×6 Matrix{Int64}: 0 3 2 8 10 13 3 0 1 5 7 10 2 1 0 6 8 11 8 5 6 0 2 5 10 7 8 2 0 3 13 10 11 5 3 0
A = D .< 5
6×6 BitMatrix: 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 1 1 0 0 0 0 1 1 1 0 0 0 0 1 1
setcover(A)
OPTIMAL
2-element Vector{Int64}: 3 5
using DataFrames, CSV
df = DataFrame(CSV.File("Net_Mod-2-Data-arc.csv"))
i, j, d = Int.(df.i), Int.(df.j), Float64.(df.d)
G = SimpleWeightedGraph(i, j, d)
{847, 1210} undirected simple Int64 graph with Float64 weights
D = vcat([dijkstra_shortest_paths(G,i).dists for i ∈ vertices(G)]'...)
847×847 Matrix{Float64}: 0.0 590.4 1594.0 354.1 2863.2 … 1746.1 1917.5 1744.2 1744.7 590.4 0.0 2136.3 704.3 3425.3 2288.4 2459.8 2286.5 2287.0 1594.0 2136.3 0.0 1785.4 2166.6 1067.1 911.4 1065.2 1065.7 354.1 704.3 1785.4 0.0 3074.4 1937.5 2108.9 1935.6 1936.1 2863.2 3425.3 2166.6 3074.4 0.0 1191.2 1483.8 1194.3 1194.8 410.8 572.7 1951.9 519.9 3240.9 … 2104.0 2275.4 2102.1 2102.6 388.4 538.0 1929.5 497.5 3218.5 2081.6 2253.0 2079.7 2080.2 2536.3 3078.6 1251.3 2727.7 1449.1 945.1 652.5 948.2 948.7 964.2 1506.5 1345.4 1155.6 2633.2 1496.3 1668.9 1494.4 1494.9 681.5 1223.8 1136.9 872.9 2424.7 1287.8 1460.4 1285.9 1286.4 2972.9 3515.2 1687.9 3164.3 1380.4 … 1343.1 1089.1 1346.2 1346.7 417.8 826.1 1358.0 475.2 2647.0 1510.1 1681.5 1508.2 1508.7 351.2 777.7 1574.0 426.8 2863.0 1726.1 1897.5 1724.2 1724.7 ⋮ ⋱ ⋮ 892.0 1434.3 702.0 1083.4 2069.5 … 932.6 1025.5 930.7 931.2 939.1 1481.4 654.9 1130.5 2116.6 979.7 978.4 977.8 978.3 968.5 1510.8 625.5 1159.9 2142.9 1006.0 949.0 1004.1 1004.6 1031.9 1574.2 562.1 1223.3 2079.5 942.6 885.6 940.7 941.2 1071.2 1613.5 728.2 1262.6 2245.6 1108.7 1051.7 1106.8 1107.3 1343.2 1885.5 533.9 1534.6 2163.4 … 1026.5 969.5 1024.6 1025.1 1743.6 2285.9 1064.6 1935.0 1193.7 2.5 295.1 0.6 1.1 1744.7 2287.0 1065.7 1936.1 1194.8 3.6 296.2 0.5 1.0 1746.1 2288.4 1067.1 1937.5 1191.2 0.0 292.6 3.1 3.6 1917.5 2459.8 911.4 2108.9 1483.8 292.6 0.0 295.7 296.2 1744.2 2286.5 1065.2 1935.6 1194.3 … 3.1 295.7 0.0 0.5 1744.7 2287.0 1065.7 1936.1 1194.8 3.6 296.2 0.5 0.0
A = D .< 450
idx = setcover(A)
OPTIMAL
14-element Vector{Int64}: 32 79 109 143 367 390 452 486 693 715 742 771 804 820
using Plots
df = DataFrame(CSV.File("Net_Mod-2-Data-node.csv"))
x, y = Float64.(df.x), Float64.(df.y)
plot([x[i]'; x[j]'], [y[i]'; y[j]'], label=false)
scatter!([x[idx]], [y[idx]], label="Terminal")
BIP formulation of the set packing problem:
$ \begin{eqnarray*} \quad \mbox{Maximize} \quad \sum_{i \in N} x_i \\ \quad \mbox{subject to} \quad \sum_{i \in N} a_{ji} x_i &\le& 1, &\quad j \in M\\ x_i &\in& \bigl\{ 0, 1 \bigr\}. &\quad i \in N \end{eqnarray*} $
M = 1:6
N = 1:5
Mi = [[1,2],[1,4,5],[3,5],[2,3,6],[6]]
A = zeros(length(M),length(N)) # A = objects x subsets
for i = 1:size(A,2)
A[Mi[i],i] .= 1
end
using JuMP, Cbc
model = Model(Cbc.Optimizer)
@variable(model, x[1:length(N)], Bin )
@objective(model, Max, sum(x[i] for i ∈ N) )
@constraint(model, [j ∈ M], sum(A[j,i]*x[i] for i ∈ N) <= 1 )
set_optimizer_attribute(model, "logLevel", 1)
optimize!(model)
TCᵒ, xᵒ = objective_value(model), value.(x)
Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Jan 1 1970 command line - Cbc_C_Interface -logLevel 1 -solve -quit (default strategy 1) Continuous objective value is 3 - 0.00 seconds Cgl0004I processed model has 0 rows, 0 columns (0 integer (0 of which binary)) and 0 elements Cbc3007W No integer variables - nothing to do Cuts at root node changed objective from -3 to -1.79769e+308 Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 3.00000000 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.00 Time (Wallclock seconds): 0.00 Total time (CPU seconds): 0.00 (Wallclock seconds): 0.00
(3.0, [1.0, 0.0, 1.0, 0.0, 1.0])
The bin packing problem involves determining the minimum number of equal-capacity bins required to pack different size objects so that all of the objects assigned to a bin do not exceed its capacity. The 1-D bin packing problem refers to the bins having a single scalar capacity and each object having a single scalar size. The 2-D bin packing problem can, for example, refer to each bin having both weight and cubic volume restrictions on its capacity and each object having a weight and cubic volume.
General formulation of the 1-D bin packing problem (where the capacity restruction is assumed to be cubic volume):
$ \begin{eqnarray*} \quad M &=& \bigl\{ 1, \dots, m \bigr\}, \quad \mbox{objects to be packed into bins} \\ v_j &=& \mbox{volume of object } n \\ V &=& \mbox{volume of each bin } B_i \bigl( \max v_j \le V \bigr) \\ B^* &=& \mathrm{arg}\underset{B}{\operatorname{min}} \bigl\{ \lvert B \rvert : \sum_{j \in B_i} v_j \le V, \bigcup_{B_i \in B} = M \bigr\}, \quad \mbox{min cost bin packing of } M. \end{eqnarray*}$
BIP formulation of the problem:
$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{i \in N} y_i \\ \quad \mbox{subject to} \quad V y_i &\ge& \sum_{j \in M} v_j x_{ij}, &\quad i \in N &\quad (a) \\ \sum_{i \in N} x_{ij} &=& 1, &\quad j \in M &\quad (b) \\ y_i &\in& \bigl\{ 0, 1 \bigr\}, &\quad i \in N \\ x_{ij} &\in& \bigl\{ 0, 1 \bigr\}, &\quad i \in N, j \in M \end{eqnarray*} $
where
$ \begin{eqnarray*} \quad y_i &=& \begin{cases} 1, & \mbox{if bin } B_i \mbox{ is used in packing} \\ 0, & \mbox{otherwise} \end{cases} \\ x_{ij} &=& \begin{cases} 1, & \mbox{if object } j \mbox{ is packed into bin } B_i \\ 0, & \mbox{otherwise.} \end{cases} \\ \end{eqnarray*} $
In the formulation, there are a set N of potential bins. Constraints (a) ensure that, for each bin, the total volume of the assigned objects do not exceed its capacity, and constraints (b) ensure that, for each object, it is assigned to one bin.
Pack 20 objects ranging randomly in size from 1 to 5 into min number of bins of capacity 10:
using Random
m = 20
M = 1:m
V = 10
Random.seed!(1244)
v = rand(1:5,m)
20-element Vector{Int64}: 3 4 2 3 3 4 1 5 2 1 3 5 4 4 5 3 4 4 3 3
LB = ceil(sum(v)/V)
7.0
using JuMP, Cbc
model = Model(Cbc.Optimizer)
@variable(model, y[1:m], Bin )
@variable(model, X[1:m,1:m], Bin )
@objective(model, Min, sum(y[i] for i ∈ M) )
@constraint(model, [i ∈ M], V*y[i] >= sum(v[j]*X[i,j] for j ∈ M) )
@constraint(model, [j ∈ M], sum(X[i,j] for i ∈ M) == 1 )
set_optimizer_attribute(model, "logLevel", 1)
optimize!(model)
TCᵒ, yᵒ, Xᵒ = objective_value(model), value.(y), value.(X)
Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Jan 1 1970 command line - Cbc_C_Interface -logLevel 1 -solve -quit (default strategy 1) Continuous objective value is 6.6 - 0.00 seconds Cgl0004I processed model has 40 rows, 420 columns (420 integer (420 of which binary)) and 820 elements Cutoff increment increased from 1e-05 to 0.9999 Cbc0038I Initial state - 11 integers unsatisfied sum - 3.53333 Cbc0038I Pass 1: suminf. 1.90000 (5) obj. 6.9 iterations 18 Cbc0038I Pass 2: suminf. 1.30000 (5) obj. 8 iterations 20 Cbc0038I Pass 3: suminf. 1.00000 (3) obj. 8 iterations 3 Cbc0038I Pass 4: suminf. 1.00000 (3) obj. 8 iterations 0 Cbc0038I Pass 5: suminf. 0.00000 (0) obj. 8 iterations 12 Cbc0038I Solution found of 8 Cbc0038I Before mini branch and bound, 405 integers at bound fixed and 0 continuous Cbc0038I Full problem 40 rows 420 columns, reduced to 0 rows 0 columns Cbc0038I Mini branch and bound did not improve solution (0.01 seconds) Cbc0038I Round again with cutoff of 6.96009 Cbc0038I Pass 6: suminf. 1.90000 (5) obj. 6.9 iterations 0 Cbc0038I Pass 7: suminf. 1.90000 (5) obj. 6.9 iterations 4 Cbc0038I Pass 8: suminf. 1.70000 (5) obj. 6.9 iterations 6 Cbc0038I Pass 9: suminf. 2.23991 (9) obj. 6.96009 iterations 22 Cbc0038I Pass 10: suminf. 1.03991 (5) obj. 6.96009 iterations 13 Cbc0038I Pass 11: suminf. 1.10000 (5) obj. 6.7 iterations 6 Cbc0038I Pass 12: suminf. 0.83991 (5) obj. 6.96009 iterations 6 Cbc0038I Pass 13: suminf. 1.03991 (5) obj. 6.96009 iterations 6 Cbc0038I Pass 14: suminf. 2.93991 (10) obj. 6.96009 iterations 52 Cbc0038I Pass 15: suminf. 2.93991 (10) obj. 6.96009 iterations 39 Cbc0038I Pass 16: suminf. 3.03991 (10) obj. 6.96009 iterations 25 Cbc0038I Pass 17: suminf. 2.03991 (9) obj. 6.96009 iterations 9 Cbc0038I Pass 18: suminf. 2.43991 (9) obj. 6.96009 iterations 14 Cbc0038I Pass 19: suminf. 2.03991 (9) obj. 6.96009 iterations 11 Cbc0038I Pass 20: suminf. 2.33991 (8) obj. 6.96009 iterations 65 Cbc0038I Pass 21: suminf. 1.33991 (7) obj. 6.96009 iterations 34 Cbc0038I Pass 22: suminf. 2.16009 (7) obj. 6.96009 iterations 27 Cbc0038I Pass 23: suminf. 1.63991 (7) obj. 6.96009 iterations 19 Cbc0038I Pass 24: suminf. 1.96009 (7) obj. 6.96009 iterations 22 Cbc0038I Pass 25: suminf. 1.63991 (7) obj. 6.96009 iterations 14 Cbc0038I Pass 26: suminf. 1.23991 (6) obj. 6.96009 iterations 15 Cbc0038I Pass 27: suminf. 1.23991 (7) obj. 6.96009 iterations 8 Cbc0038I Pass 28: suminf. 2.43991 (8) obj. 6.96009 iterations 23 Cbc0038I Pass 29: suminf. 2.16009 (7) obj. 6.96009 iterations 10 Cbc0038I Pass 30: suminf. 1.63991 (7) obj. 6.96009 iterations 14 Cbc0038I Pass 31: suminf. 1.23991 (6) obj. 6.96009 iterations 18 Cbc0038I Pass 32: suminf. 1.23991 (7) obj. 6.96009 iterations 9 Cbc0038I Pass 33: suminf. 2.43991 (8) obj. 6.96009 iterations 21 Cbc0038I Pass 34: suminf. 2.16009 (7) obj. 6.96009 iterations 10 Cbc0038I Pass 35: suminf. 1.63991 (7) obj. 6.96009 iterations 12 Cbc0038I No solution found this major pass Cbc0038I Before mini branch and bound, 345 integers at bound fixed and 0 continuous Cbc0038I Full problem 40 rows 420 columns, reduced to 26 rows 66 columns Cbc0038I Mini branch and bound improved solution from 8 to 7 (0.07 seconds) Cbc0038I After 0.07 seconds - Feasibility pump exiting with objective of 7 - took 0.06 seconds Cbc0012I Integer solution of 7 found by feasibility pump after 0 iterations and 0 nodes (0.07 seconds) Cbc0001I Search completed - best objective 7, took 0 iterations and 0 nodes (0.07 seconds) Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost Cuts at root node changed objective from 6.6 to 6.6 Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 7.00000000 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.07 Time (Wallclock seconds): 0.07 Total time (CPU seconds): 0.07 (Wallclock seconds): 0.07
(7.0, [1.0, 1.0, 1.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, 1.0, 1.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])
Xᵒ.!=0
20×20 BitMatrix: 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 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 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 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 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 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 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 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 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 1 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 1 0 1 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
B = [findall(Xᵒ[i,:] .!= 0) for i ∈ findall(yᵒ.>0)]
7-element Vector{Vector{Int64}}: [1, 9, 10, 17] [4, 11, 18] [7, 13, 15] [6, 16, 20] [8, 12] [2, 19] [3, 5, 14]
[v[B[i]] for i=1:length(B)]
7-element Vector{Vector{Int64}}: [3, 2, 1, 4] [3, 3, 4] [1, 4, 5] [4, 3, 3] [5, 5] [4, 3] [2, 3, 4]
[sum(v[B[i]]) for i=1:length(B)]
7-element Vector{Int64}: 10 10 10 10 10 7 9
Pack 200 objects ranging randomly in size from 1 to 10 into min number of bins of capacity 20:
using Random
m = 200
M = 1:m
V = 20
Random.seed!(57839)
v = rand(1:10,m)
LB = ceil(sum(v)/V)
55.0
using JuMP, Cbc
model = Model(Cbc.Optimizer)
@variable(model, y[1:m], Bin )
@variable(model, X[1:m,1:m], Bin )
@objective(model, Min, sum(y[i] for i ∈ M) )
@constraint(model, [i ∈ M], V*y[i] >= sum(v[j]*X[i,j] for j ∈ M) )
@constraint(model, [j ∈ M], sum(X[i,j] for i ∈ M) == 1 )
set_optimizer_attribute(model, "logLevel", 1)
set_optimizer_attribute(model, "seconds", 30.0)
optimize!(model)
Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Jan 1 1970 command line - Cbc_C_Interface -seconds 30.0 -logLevel 1 -solve -quit (default strategy 1) seconds was changed from 1e+100 to 30 Continuous objective value is 54.85 - 0.43 seconds Cgl0004I processed model has 400 rows, 40200 columns (40200 integer (40200 of which binary)) and 80200 elements Cutoff increment increased from 1e-05 to 0.9999 Cbc0038I Initial state - 101 integers unsatisfied sum - 30.2254 Cbc0038I Pass 1: suminf. 19.26825 (77) obj. 59.55 iterations 504 Cbc0038I Pass 2: suminf. 19.06825 (73) obj. 59.55 iterations 43 Cbc0038I Pass 3: suminf. 12.11825 (56) obj. 66.35 iterations 144 Cbc0038I Pass 4: suminf. 10.31825 (51) obj. 66.35 iterations 38 Cbc0038I Pass 5: suminf. 6.90794 (38) obj. 73 iterations 79 Cbc0038I Pass 6: suminf. 6.10794 (36) obj. 73 iterations 9 Cbc0038I Pass 7: suminf. 4.24127 (28) obj. 73 iterations 44 Cbc0038I Pass 8: suminf. 4.01905 (26) obj. 73 iterations 12 Cbc0038I Pass 9: suminf. 3.78333 (20) obj. 75.45 iterations 65 Cbc0038I Pass 10: suminf. 3.58333 (21) obj. 75.45 iterations 18 Cbc0038I Pass 11: suminf. 3.13889 (19) obj. 75.45 iterations 35 Cbc0038I Pass 12: suminf. 2.93889 (20) obj. 75.45 iterations 19 Cbc0038I Pass 13: suminf. 3.33333 (13) obj. 82 iterations 64 Cbc0038I Pass 14: suminf. 2.93333 (11) obj. 82 iterations 4 Cbc0038I Pass 15: suminf. 2.06667 (8) obj. 82 iterations 22 Cbc0038I Pass 16: suminf. 2.06667 (9) obj. 82 iterations 3 Cbc0038I Pass 17: suminf. 2.48889 (9) obj. 82 iterations 26 Cbc0038I Pass 18: suminf. 2.48889 (9) obj. 82 iterations 2 Cbc0038I Pass 19: suminf. 12.94048 (43) obj. 82 iterations 460 Cbc0038I Pass 20: suminf. 7.38810 (28) obj. 82 iterations 615 Cbc0038I Pass 21: suminf. 6.78810 (26) obj. 82 iterations 10 Cbc0038I Pass 22: suminf. 3.02937 (16) obj. 83.15 iterations 48 Cbc0038I Pass 23: suminf. 2.03571 (12) obj. 83.15 iterations 26 Cbc0038I Pass 24: suminf. 2.80714 (10) obj. 83.15 iterations 37 Cbc0038I Pass 25: suminf. 2.52143 (10) obj. 83.15 iterations 7 Cbc0038I Pass 26: suminf. 2.45714 (6) obj. 85 iterations 37 Cbc0038I Pass 27: suminf. 1.97143 (6) obj. 85 iterations 17 Cbc0038I Pass 28: suminf. 12.84048 (40) obj. 82.6 iterations 510 Cbc0038I Pass 29: suminf. 7.30556 (26) obj. 82.6 iterations 691 Cbc0038I Pass 30: suminf. 6.30556 (22) obj. 82.6 iterations 9 Cbc0038I No solution found this major pass Cbc0038I Before mini branch and bound, 39664 integers at bound fixed and 0 continuous Cbc0038I Full problem 400 rows 40200 columns, reduced to 160 rows 430 columns Cbc0038I Mini branch and bound improved solution from 1.79769e+308 to 64 (2.09 seconds) Cbc0038I Round again with cutoff of 62.1851 Cbc0038I Pass 30: suminf. 19.26825 (77) obj. 59.55 iterations 0 Cbc0038I Pass 31: suminf. 17.76825 (70) obj. 61.05 iterations 51 Cbc0038I Pass 32: suminf. 17.11650 (73) obj. 62.1851 iterations 186 Cbc0038I Pass 33: suminf. 14.81650 (65) obj. 62.1851 iterations 66 Cbc0038I Pass 34: suminf. 13.84507 (61) obj. 62.1851 iterations 51 Cbc0038I Pass 35: suminf. 13.34507 (62) obj. 62.1851 iterations 18 Cbc0038I Pass 36: suminf. 13.34507 (62) obj. 62.1851 iterations 3 Cbc0038I Pass 37: suminf. 22.42681 (85) obj. 62.1851 iterations 523 Cbc0038I Pass 38: suminf. 14.35618 (70) obj. 62.1851 iterations 666 Cbc0038I Pass 39: suminf. 13.17840 (68) obj. 62.1851 iterations 24 Cbc0038I Pass 40: suminf. 12.55618 (65) obj. 62.1851 iterations 59 Cbc0038I Pass 41: suminf. 14.52840 (70) obj. 62.1851 iterations 126 Cbc0038I Pass 42: suminf. 14.32840 (70) obj. 62.1851 iterations 13 Cbc0038I Pass 43: suminf. 13.61729 (66) obj. 62.1851 iterations 156 Cbc0038I Pass 44: suminf. 13.21729 (67) obj. 62.1851 iterations 64 Cbc0038I Pass 45: suminf. 13.63396 (67) obj. 62.1851 iterations 54 Cbc0038I Pass 46: suminf. 13.43951 (66) obj. 62.1851 iterations 7 Cbc0038I Pass 47: suminf. 12.43951 (63) obj. 62.1851 iterations 18 Cbc0038I Pass 48: suminf. 12.31729 (61) obj. 62.1851 iterations 105 Cbc0038I Pass 49: suminf. 12.01729 (61) obj. 62.1851 iterations 11 Cbc0038I Pass 50: suminf. 11.81729 (58) obj. 62.1851 iterations 52 Cbc0038I Pass 51: suminf. 11.51729 (59) obj. 62.1851 iterations 16 Cbc0038I Pass 52: suminf. 11.15062 (59) obj. 62.1851 iterations 31 Cbc0038I Pass 53: suminf. 10.85062 (57) obj. 62.1851 iterations 19 Cbc0038I Pass 54: suminf. 11.15062 (59) obj. 62.1851 iterations 46 Cbc0038I Pass 55: suminf. 11.15062 (59) obj. 62.1851 iterations 1 Cbc0038I Pass 56: suminf. 11.13158 (55) obj. 62.1851 iterations 147 Cbc0038I Pass 57: suminf. 10.83158 (56) obj. 62.1851 iterations 25 Cbc0038I Pass 58: suminf. 10.53158 (55) obj. 62.1851 iterations 79 Cbc0038I Pass 59: suminf. 10.53158 (55) obj. 62.1851 iterations 15 Cbc0038I No solution found this major pass Cbc0038I Before mini branch and bound, 39770 integers at bound fixed and 0 continuous Cbc0038I Full problem 400 rows 40200 columns, reduced to 123 rows 287 columns Cbc0038I Mini branch and bound did not improve solution (2.73 seconds) Cbc0038I After 2.73 seconds - Feasibility pump exiting with objective of 64 - took 1.52 seconds Cbc0012I Integer solution of 64 found by feasibility pump after 0 iterations and 0 nodes (2.74 seconds) Cbc0038I Full problem 400 rows 40200 columns, reduced to 51 rows 86 columns Cbc0031I 66 added rows had average density of 604.75758 Cbc0013I At root node, 66 cuts changed objective from 54.85 to 54.85 in 10 passes Cbc0014I Cut generator 0 (Probing) - 21 row cuts average 2.0 elements, 0 column cuts (0 active) in 1.756 seconds - new frequency is -100 Cbc0014I Cut generator 1 (Gomory) - 353 row cuts average 2043.2 elements, 0 column cuts (0 active) in 1.347 seconds - new frequency is -100 Cbc0014I Cut generator 2 (Knapsack) - 3 row cuts average 2.3 elements, 0 column cuts (0 active) in 0.010 seconds - new frequency is -100 Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.003 seconds - new frequency is -100 Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 178 row cuts average 505.1 elements, 0 column cuts (0 active) in 0.130 seconds - new frequency is -100 Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.035 seconds - new frequency is -100 Cbc0014I Cut generator 6 (TwoMirCuts) - 124 row cuts average 172.8 elements, 0 column cuts (0 active) in 0.206 seconds - new frequency is -100 Cbc0014I Cut generator 7 (ZeroHalf) - 121 row cuts average 216.9 elements, 0 column cuts (0 active) in 0.246 seconds - new frequency is -100 Cbc0010I After 0 nodes, 1 on tree, 64 best solution, best possible 54.85 (9.88 seconds) Cbc0038I Full problem 400 rows 40200 columns, reduced to 98 rows 249 columns Cbc0038I Full problem 400 rows 40200 columns, reduced to 101 rows 257 columns Cbc0010I After 100 nodes, 61 on tree, 64 best solution, best possible 54.85 (21.43 seconds) Cbc0038I Full problem 400 rows 40200 columns, reduced to 95 rows 245 columns Cbc0010I After 200 nodes, 112 on tree, 64 best solution, best possible 54.85 (23.84 seconds) Cbc0038I Full problem 400 rows 40200 columns, reduced to 95 rows 247 columns Cbc0010I After 300 nodes, 166 on tree, 64 best solution, best possible 54.85 (26.06 seconds) Cbc0038I Full problem 400 rows 40200 columns, reduced to 94 rows 249 columns Cbc0010I After 400 nodes, 210 on tree, 64 best solution, best possible 54.85 (28.20 seconds) Cbc0020I Exiting on maximum time Cbc0005I Partial search - best objective 64 (best possible 54.85), took 8730 iterations and 457 nodes (29.44 seconds) Cbc0032I Strong branching done 1602 times (17662 iterations), fathomed 0 nodes and fixed 2 variables Cbc0035I Maximum depth 128, 0 variables fixed on reduced cost Cuts at root node changed objective from 54.85 to 54.85 Probing was tried 10 times and created 21 cuts of which 0 were active after adding rounds of cuts (1.756 seconds) Gomory was tried 10 times and created 353 cuts of which 0 were active after adding rounds of cuts (1.347 seconds) Knapsack was tried 10 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.010 seconds) Clique was tried 10 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.003 seconds) MixedIntegerRounding2 was tried 10 times and created 178 cuts of which 0 were active after adding rounds of cuts (0.130 seconds) FlowCover was tried 10 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.035 seconds) TwoMirCuts was tried 10 times and created 124 cuts of which 0 were active after adding rounds of cuts (0.206 seconds) ZeroHalf was tried 10 times and created 121 cuts of which 0 were active after adding rounds of cuts (0.246 seconds) ImplicationCuts was tried 8 times and created 15 cuts of which 0 were active after adding rounds of cuts (0.036 seconds) Result - Stopped on time limit Objective value: 64.00000000 Lower bound: 54.850 Gap: 0.17 Enumerated nodes: 457 Total iterations: 8730 Time (CPU seconds): 29.56 Time (Wallclock seconds): 29.56 Total time (CPU seconds): 29.56 (Wallclock seconds): 29.56
has_values(model)
true
TCᵒ, yᵒ, Xᵒ = objective_value(model), value.(y), value.(X)
(64.0, [1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.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])
B = [findall(Xᵒ[i,:] .!= 0) for i ∈ findall(yᵒ.>0)]
64-element Vector{Vector{Int64}}: [1, 2, 7, 27] [3, 85, 120, 185] [5, 11] [14, 79, 196] [13, 62, 167] [23, 30, 132] [28, 36, 137] [22, 39, 41] [38, 42, 76] [21, 74, 183] [16, 24, 50] [47, 49, 86] [43, 107] ⋮ [112, 171, 180] [4, 106, 188] [146, 168, 198] [128, 159, 165, 194] [113, 125, 143, 184, 186, 197] [32, 105, 187] [25, 151, 178] [20] [145, 163, 190] [17, 124, 182] [15, 148, 191] [34, 59, 139]
objective_value(model)
64.0
objective_bound(model)
54.85
absgap = objective_value(model)/objective_bound(model)
1.1668185961713764
relgap = (objective_value(model) - objective_bound(model))/objective_value(model)
0.14296874999999998
MOI.get(model, MOI.RelativeGap())
0.14296874999999998
using JuMP, Gurobi
model = Model(Gurobi.Optimizer)
@variable(model, y[1:m], Bin )
@variable(model, X[1:m,1:m], Bin )
@objective(model, Min, sum(y[i] for i ∈ M) )
@constraint(model, [i ∈ M], V*y[i] >= sum(v[j]*X[i,j] for j ∈ M) )
@constraint(model, [j ∈ M], sum(X[i,j] for i ∈ M) == 1 )
set_optimizer_attribute(model, "OutputFlag", 1)
set_optimizer_attribute(model, "TimeLimit", 30)
optimize!(model)
Academic license - for non-commercial use only - expires 2022-07-31 Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64) Thread count: 4 physical cores, 8 logical processors, using up to 8 threads Optimize a model with 400 rows, 40200 columns and 80200 nonzeros Model fingerprint: 0x0ede761d Variable types: 0 continuous, 40200 integer (40200 binary) Coefficient statistics: Matrix range [1e+00, 2e+01] Objective range [1e+00, 1e+00] Bounds range [0e+00, 0e+00] RHS range [1e+00, 1e+00] Found heuristic solution: objective 129.0000000 Presolve time: 0.19s Presolved: 400 rows, 40200 columns, 80200 nonzeros Variable types: 0 continuous, 40200 integer (40200 binary) Root relaxation: objective 5.485000e+01, 1669 iterations, 0.15 seconds Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time 0 0 54.85000 0 1 129.00000 54.85000 57.5% - 0s H 0 0 55.0000000 54.85000 0.27% - 0s 0 0 54.85000 0 1 55.00000 54.85000 0.27% - 0s Explored 1 nodes (13360 simplex iterations) in 0.62 seconds Thread count was 8 (of 8 available processors) Solution count 2: 55 129 Optimal solution found (tolerance 1.00e-04) Best objective 5.500000000000e+01, best bound 5.500000000000e+01, gap 0.0000% User-callback calls 97, time in user-callback 0.00 sec
LB = ceil(sum(v)/V)
55.0
The graph coloring problem is to assign colors to the vertices of a graph so that any two vertices connected by an edge are assigned different colors. The edges of a graph and faces of a planar graph could be colored instead of vertices but it is always possible to convert these coloring problems to a vertex covering problem. The graph covering problem has many applications associated with assigning resources to activities where each activity is represented by a vertex, and each edge represents pairs of activities that cannot be assigned to the same resource due to some conflict. Each resource is represented by a different color, and it is usually desired to find the minimum number of colors (resources) required for all of the activities.
General formulation of the problem:
$ \begin{eqnarray*} \quad M &=& \bigl\{ 1, \dots, m \bigr\}, \quad \mbox{set of vertices in graph } G \\ W &=& \mbox{adjacency matrix of } G \\ K^* &=& \mathrm{arg}\underset{K}{\operatorname{min}} \bigl\{ \lvert K \rvert : (i,j) \in W, i \in K_r, j \in K_s, \bigcup_{K_r,K_s \in K} = M \bigr\}, \quad \mbox{min coloring of } G. \end{eqnarray*}$
BIP formulation of the problem:
$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{k \in N} y_k \\ \quad \mbox{subject to} \quad \sum_{k \in N} x_{ik} &=& 1, &\quad i \in M &\quad (a) \\ x_{ik} + x_{jk} &\le& 1, &\quad (i,j) \in W; k \in N &\quad (b) \\ x_{ik} &\le& y_k, &\quad i \in M, k \in N &\quad (c) \\ y_i &\in& \bigl\{ 0, 1 \bigr\}, &\quad i \in M \\ x_{ij} &\in& \bigl\{ 0, 1 \bigr\}, &\quad i \in M; j \in M \end{eqnarray*} $
where
$ \begin{eqnarray*} \quad y_k &=& \begin{cases} 1, & \mbox{if color } k \mbox{ is used} \\ 0, & \mbox{otherwise} \end{cases} \\ x_{ik} &=& \begin{cases} 1, & \mbox{if vertex } i \mbox{ is colored } k \\ 0, & \mbox{otherwise.} \end{cases} \\ \end{eqnarray*} $
In the formulation, there are a set N of potential colors. Constraints (a) ensure that, for each vertex, it is assigned to one color. Constraints (b) ensure that, for each edge (i,j) of W, if vertex i is assigned to color k then vertex j is not assigned to k (i.e., If X then not Y). Constraints (c) ensure that a vertex can be colored k only if color k is being used.
W = [0 1 1 1 1; 0 0 1 1 0; 0 0 0 1 0; 0 0 0 0 1; 0 0 0 0 0]
W = W + W'
5×5 Matrix{Int64}: 0 1 1 1 1 1 0 1 1 0 1 1 0 1 0 1 1 1 0 1 1 0 0 1 0
using LightGraphs, SimpleWeightedGraphs
G = SimpleWeightedGraph(W)
b = src.(edges(G))
e = dst.(edges(G))
[b e]
8×2 Matrix{Int64}: 1 2 1 3 2 3 1 4 2 4 3 4 1 5 4 5
zip(b,e)
zip([1, 1, 2, 1, 2, 3, 1, 4], [2, 3, 3, 4, 4, 4, 5, 5])
[(i,j) for (i,j) in zip(b,e)]
8-element Vector{Tuple{Int64, Int64}}: (1, 2) (1, 3) (2, 3) (1, 4) (2, 4) (3, 4) (1, 5) (4, 5)
[i+j for (i,j) in zip(b,e)]
8-element Vector{Int64}: 3 4 5 5 6 7 6 9
using JuMP, Cbc
model = Model(Cbc.Optimizer)
M = 1:size(W,1)
@variable(model, y[1:length(M)], Bin )
@variable(model, X[1:length(M),1:length(M)], Bin )
@objective(model, Min, sum(y[i] for i ∈ M ))
@constraint(model, [i ∈ M], sum(X[i,k] for k ∈ M) == 1 )
@constraint(model, [(i,j) ∈ zip(b,e), k ∈ M], X[i,k] + X[j,k] <= 1 )
@constraint(model, [i ∈ M, k ∈ M], X[i,k] <= y[k] )
set_optimizer_attribute(model, "logLevel", 1)
set_optimizer_attribute(model, "seconds", 30.0)
optimize!(model)
Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Jan 1 1970 command line - Cbc_C_Interface -seconds 30.0 -logLevel 1 -solve -quit (default strategy 1) seconds was changed from 1e+100 to 30 Continuous objective value is 1 - 0.00 seconds Cgl0003I 0 fixed, 0 tightened bounds, 40 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 41 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 33 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 17 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 10 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 5 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 5 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 3 strengthened rows, 0 substitutions Cgl0004I processed model has 15 rows, 30 columns (30 integer (30 of which binary)) and 70 elements Cutoff increment increased from 1e-05 to 0.9999 Cbc0038I Initial state - 0 integers unsatisfied sum - 0 Cbc0038I Solution found of 4 Cbc0038I Before mini branch and bound, 30 integers at bound fixed and 0 continuous Cbc0038I Mini branch and bound did not improve solution (0.01 seconds) Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of 4 - took 0.00 seconds Cbc0012I Integer solution of 4 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds) Cbc0001I Search completed - best objective 4, took 0 iterations and 0 nodes (0.01 seconds) Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost Cuts at root node changed objective from 4 to 4 Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 4.00000000 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.01 Time (Wallclock seconds): 0.01 Total time (CPU seconds): 0.01 (Wallclock seconds): 0.01
TCᵒ, yᵒ, Xᵒ = objective_value(model), value.(y), value.(X)
@show TCᵒ
@show yᵒ
Xᵒ
TCᵒ = 4.0 yᵒ = [1.0, 0.0, 1.0, 1.0, 1.0]
5×5 Matrix{Float64}: 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.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 0.0 1.0 0.0 0.0
Kᵒ = [findall(Xᵒ[:,i] .!= 0) for i ∈ findall(yᵒ.>0)]
4-element Vector{Vector{Int64}}: [3] [2, 5] [1] [4]
using DataFrames, CSV
df = DataFrame(CSV.File("MIP-3-Data-state_adj.csv", header=true))
first(df,5)
5 rows × 54 columns (omitted printing of 45 columns)
STATEFP | STNAME | STABBR | POPULATION | LATITUDE | LONGITUDE | AL | AZ | AR | |
---|---|---|---|---|---|---|---|---|---|
Int64 | String | String | Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | |
1 | 1 | Alabama | AL | 4779736 | 33.0081 | -86.7568 | 0 | 0 | 0 |
2 | 4 | Arizona | AZ | 6392017 | 33.3683 | -111.864 | 0 | 0 | 0 |
3 | 5 | Arkansas | AR | 2915918 | 35.1426 | -92.6552 | 0 | 0 | 0 |
4 | 6 | California | CA | 37253956 | 35.4636 | -119.325 | 0 | 1 | 0 |
5 | 8 | Colorado | CO | 5029196 | 39.5134 | -105.208 | 0 | 0 | 0 |
W = Matrix(DataFrame(CSV.File("MIP-3-Data-state_adj.csv", header=true)))
W = Int.(W[:,7:end])
48×48 Matrix{Int64}: 0 0 0 0 0 0 0 1 1 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 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 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 0 0 0 0 0 0 1 0 0 0 0 0 1 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 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 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 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 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 … 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 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 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 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 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 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 1 0 0 0 0 1 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 0 0 0 0 0 0 1 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 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 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0
using LightGraphs, SimpleWeightedGraphs
G = SimpleWeightedGraph(W)
b = src.(edges(G))
e = dst.(edges(G));
using Plots
lon, lat, str = df.LONGITUDE, df.LATITUDE, df.STABBR
plot([lon[b]'; lon[e]'], [lat[b]'; lat[e]'], label=false)
scatter!([lon], [lat], label=false)
annotate!([(lon[i],lat[i],(" "*str[i],7,:left)) for i in 1:length(str)])
using JuMP, Cbc
model = Model(Cbc.Optimizer)
M = 1:size(W,1)
@variable(model, y[1:length(M)], Bin )
@variable(model, X[1:length(M),1:length(M)], Bin )
@objective(model, Min, sum(y[i] for i ∈ M ))
@constraint(model, [i ∈ M], sum(X[i,k] for k ∈ M) == 1 )
@constraint(model, [(i,j) ∈ zip(b,e), k ∈ M], X[i,k] + X[j,k] <= 1 )
@constraint(model, [i ∈ M, k ∈ M], X[i,k] <= y[k] )
set_optimizer_attribute(model, "logLevel", 1)
set_optimizer_attribute(model, "seconds", 30.0)
optimize!(model)
Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Jan 1 1970 command line - Cbc_C_Interface -seconds 30.0 -logLevel 1 -solve -quit (default strategy 1) seconds was changed from 1e+100 to 30 Continuous objective value is 1 - 0.07 seconds Cgl0003I 0 fixed, 0 tightened bounds, 5099 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 4688 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 2589 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 1410 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 342 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 54 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 39 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 24 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 5 strengthened rows, 0 substitutions Cgl0004I processed model has 2887 rows, 2352 columns (2352 integer (2352 of which binary)) and 13617 elements Cutoff increment increased from 1e-05 to 0.9999 Cbc0045I 48 integer variables out of 2352 objects (2352 integer) have cost of 1 - high priority Cbc0045I branch on satisfied N create fake objective Y random cost Y Cbc0038I Initial state - 248 integers unsatisfied sum - 38.7673 Cbc0038I Pass 1: suminf. 22.00000 (44) obj. 4 iterations 1914 Cbc0038I Pass 2: suminf. 0.00000 (0) obj. 6 iterations 596 Cbc0038I Solution found of 6 Cbc0038I Before mini branch and bound, 2063 integers at bound fixed and 0 continuous Cbc0038I Full problem 2887 rows 2352 columns, reduced to 342 rows 287 columns Cbc0038I Mini branch and bound improved solution from 6 to 4 (3.77 seconds) Cbc0038I After 3.77 seconds - Feasibility pump exiting with objective of 4 - took 0.35 seconds Cbc0012I Integer solution of 4 found by feasibility pump after 0 iterations and 0 nodes (3.77 seconds) Cbc0001I Search completed - best objective 4, took 0 iterations and 0 nodes (3.77 seconds) Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost Cuts at root node changed objective from 4 to 4 Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 4.00000000 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 3.86 Time (Wallclock seconds): 3.86 Total time (CPU seconds): 3.86 (Wallclock seconds): 3.86
yᵒ, Xᵒ = value.(y), value.(X)
Kᵒ = [findall(Xᵒ[:,i] .!= 0) for i ∈ findall(yᵒ.>0)]
4-element Vector{Vector{Int64}}: [9, 15, 19, 22, 26, 29, 36, 48] [1, 2, 3, 13, 14, 20, 27, 31, 32, 35, 37, 46] [4, 5, 7, 8, 10, 12, 17, 23, 30, 38, 39, 41, 44, 47] [6, 11, 16, 18, 21, 24, 25, 28, 33, 34, 40, 42, 43, 45]
using Plots
plot([lon[b]'; lon[e]'], [lat[b]'; lat[e]'], label=false)
for i = 1:length(Kᵒ)
scatter!([lon[Kᵒ[i]]], [lat[Kᵒ[i]]], label="Color "*string(i) )
end
annotate!([(lon[i],lat[i],(" "*str[i],7,:left)) for i in 1:length(str)], legend=:bottomleft )
In 2019, a total of 30 PhD students took either the OR or ISE Qualifying Exam. Each student selected four areas to be tested in from eleven available areas. The portion of the exam for each area is offered on a different day, and multiple areas can be scheduled for the same day as long as no student is taking both areas. The objective is to determine the minimum number of days required for the exam.
L = Matrix(DataFrame(CSV.File("MIP-3-Data-QE_2019,csv", header=false)))
30×4 Matrix{Int64}: 1 3 4 5 1 2 4 8 1 5 7 8 5 6 7 8 4 6 7 8 5 6 7 8 1 5 6 8 1 2 4 6 3 4 5 6 1 3 5 6 4 6 7 8 4 6 7 8 4 6 7 8 ⋮ 1 3 4 6 1 2 3 4 6 7 8 9 7 8 9 10 7 8 9 11 5 7 8 9 6 7 8 9 7 8 9 10 7 8 9 10 6 7 8 9 7 8 9 10 10 9 8 7
m = maximum(L)
11
W = Int.(zeros(m,m))
for k in eachrow(L)
for i = 1:length(k)-1, j = i+1:length(k)
W[k[i],k[j]] = 1
end
end
W
11×11 Matrix{Int64}: 0 1 1 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
W = W+W'
11×11 Matrix{Int64}: 0 1 1 1 1 1 1 1 0 0 0 1 0 1 1 1 1 0 1 0 0 0 1 1 0 1 1 1 0 0 0 0 0 1 1 1 0 1 1 1 1 0 0 0 1 1 1 1 0 1 1 1 1 0 0 1 1 1 1 1 0 1 1 1 0 0 1 0 0 1 1 1 0 2 2 2 1 1 1 0 1 1 1 2 0 2 2 1 0 0 0 0 1 1 2 2 0 2 1 0 0 0 0 0 0 2 2 2 0 0 0 0 0 0 0 0 1 1 1 0 0
W = Int.(zeros(m,m))
for k in eachrow(L)
k = sort(k)
for i = 1:length(k)-1, j = i+1:length(k)
W[k[i],k[j]] = 1
end
end
W = W+W'
11×11 Matrix{Int64}: 0 1 1 1 1 1 1 1 0 0 0 1 0 1 1 1 1 0 1 0 0 0 1 1 0 1 1 1 0 0 0 0 0 1 1 1 0 1 1 1 1 0 0 0 1 1 1 1 0 1 1 1 1 0 0 1 1 1 1 1 0 1 1 1 0 0 1 0 0 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 0 0 0 0 1 1 1 1 0 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 1 1 1 0 0
using LightGraphs, SimpleWeightedGraphs
G = SimpleWeightedGraph(W)
b = src.(edges(G))
e = dst.(edges(G));
using JuMP, Cbc
model = Model(Cbc.Optimizer)
M = 1:size(W,1)
@variable(model, y[1:length(M)], Bin )
@variable(model, X[1:length(M),1:length(M)], Bin )
@objective(model, Min, sum(y[i] for i ∈ M ))
@constraint(model, [i ∈ M], sum(X[i,k] for k ∈ M) == 1 )
@constraint(model, [(i,j) ∈ zip(b,e), k ∈ M], X[i,k] + X[j,k] <= 1 )
@constraint(model, [i ∈ M, k ∈ M], X[i,k] <= y[k] )
set_optimizer_attribute(model, "logLevel", 1)
set_optimizer_attribute(model, "seconds", 30.0)
optimize!(model)
Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Jan 1 1970 command line - Cbc_C_Interface -seconds 30.0 -logLevel 1 -solve -quit (default strategy 1) seconds was changed from 1e+100 to 30 Continuous objective value is 1 - 0.00 seconds Cgl0003I 0 fixed, 0 tightened bounds, 378 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 419 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 392 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 290 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 179 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 101 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 60 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 37 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 34 strengthened rows, 0 substitutions Cgl0004I processed model has 160 rows, 132 columns (132 integer (132 of which binary)) and 1058 elements Cutoff increment increased from 1e-05 to 0.9999 Cbc0045I 11 integer variables out of 132 objects (132 integer) have cost of 1 - high priority Cbc0045I branch on satisfied N create fake objective Y random cost Y Cbc0038I Initial state - 0 integers unsatisfied sum - 0 Cbc0038I Solution found of 6 Cbc0038I Before mini branch and bound, 132 integers at bound fixed and 0 continuous Cbc0038I Mini branch and bound did not improve solution (0.08 seconds) Cbc0038I After 0.08 seconds - Feasibility pump exiting with objective of 6 - took 0.00 seconds Cbc0012I Integer solution of 6 found by feasibility pump after 0 iterations and 0 nodes (0.08 seconds) Cbc0001I Search completed - best objective 6, took 0 iterations and 0 nodes (0.08 seconds) Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost Cuts at root node changed objective from 6 to 6 Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 6.00000000 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.08 Time (Wallclock seconds): 0.08 Total time (CPU seconds): 0.08 (Wallclock seconds): 0.08
yᵒ, Xᵒ = value.(y), value.(X)
Kᵒ = [findall(Xᵒ[:,i] .!= 0) for i ∈ findall(yᵒ.>0)]
6-element Vector{Vector{Int64}}: [4] [3, 8] [2, 7] [1, 9] [5, 11] [6, 10]
Issue: What if each student took a different number of exams so that each row of L
did not have the same number of elements?
Instead of reading a CSV, can create a string of data for three students taking 4, 2, and 3 exams, respectively:
data = """
1,3,2,4
4,5
2,3,5
"""
"1,3,2,4\n4,5\n2,3,5\n"
Matrix(DataFrame(CSV.File(IOBuffer(data); header=false)))
┌ Warning: thread = 1 warning: only found 2 / 4 columns around data row: 2. Filling remaining columns with `missing` └ @ CSV C:\Users\kay\.julia\packages\CSV\Zl2ww\src\file.jl:612 ┌ Warning: thread = 1 warning: only found 3 / 4 columns around data row: 3. Filling remaining columns with `missing` └ @ CSV C:\Users\kay\.julia\packages\CSV\Zl2ww\src\file.jl:612
3×4 Matrix{Union{Missing, Int64}}: 1 3 2 4 4 5 missing missing 2 3 5 missing
coalesce.(Matrix(DataFrame(CSV.File(IOBuffer(data); header=false, silencewarnings=true))), 0)
3×4 Matrix{Int64}: 1 3 2 4 4 5 0 0 2 3 5 0
Convert matrix into a nested array:
L = [r[r.!=0] for r in eachrow(coalesce.(
Matrix(DataFrame(CSV.File(IOBuffer(data); header=false, silencewarnings=true))), 0))]
3-element Vector{Vector{Int64}}: [1, 3, 2, 4] [4, 5] [2, 3, 5]
m = maximum(vcat(L...)) # Using ... and vcat to convert to non-nested array
5
W = Int.(zeros(m,m))
for k in L # Using L instead of eachrow(L) since nested array not matrix
k = sort(k)
for i = 1:length(k)-1, j = i+1:length(k)
W[k[i],k[j]] = 1
end
end
W = W+W'
5×5 Matrix{Int64}: 0 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 0 1 0 1 1 1 0