|
| 1 | +""" |
| 2 | +$TYPEDSIGNATURES |
| 3 | +
|
| 4 | +Returns the optimal solution of the Stochastic VSP instance, by solving the associated compact MIP. |
| 5 | +Quadratic constraints are linearized using Mc Cormick linearization. |
| 6 | +Note: If you have Gurobi, use `grb_model` as `model_builder` instead of `highs_model`. |
| 7 | +""" |
| 8 | +function compact_linearized_mip(instance::Instance; model_builder=scip_model, silent=true) |
| 9 | + (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance |
| 10 | + nb_nodes = nv(graph) |
| 11 | + job_indices = 2:(nb_nodes - 1) |
| 12 | + nodes = 1:nb_nodes |
| 13 | + |
| 14 | + # Pre-processing |
| 15 | + ε = intrinsic_delays |
| 16 | + Rmax = maximum(sum(ε; dims=1)) |
| 17 | + nb_scenarios = size(ε, 2) |
| 18 | + Ω = 1:nb_scenarios |
| 19 | + |
| 20 | + # Model definition |
| 21 | + model = model_builder() |
| 22 | + silent && set_silent(model) |
| 23 | + |
| 24 | + # Variables and objective function |
| 25 | + @variable(model, y[u in nodes, v in nodes; has_edge(graph, u, v)], Bin) |
| 26 | + @variable(model, R[v in nodes, ω in Ω] >= 0) # propagated delay of job v |
| 27 | + @variable(model, yR[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)] >= 0) # yR[u, v] = y[u, v] * R[u, ω] |
| 28 | + |
| 29 | + @objective( |
| 30 | + model, |
| 31 | + Min, |
| 32 | + delay_cost * sum(sum(R[v, ω] for v in job_indices) for ω in Ω) / nb_scenarios # average total delay |
| 33 | + + |
| 34 | + vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles |
| 35 | + ) |
| 36 | + |
| 37 | + # Flow contraints |
| 38 | + @constraint( |
| 39 | + model, |
| 40 | + flow[i in job_indices], |
| 41 | + sum(y[j, i] for j in inneighbors(graph, i)) == |
| 42 | + sum(y[i, j] for j in outneighbors(graph, i)) |
| 43 | + ) |
| 44 | + @constraint( |
| 45 | + model, |
| 46 | + unit_demand[i in job_indices], |
| 47 | + sum(y[j, i] for j in inneighbors(graph, i)) == 1 |
| 48 | + ) |
| 49 | + |
| 50 | + # Delay propagation constraints |
| 51 | + @constraint(model, [ω in Ω], R[1, ω] == ε[1, ω]) |
| 52 | + @constraint(model, R_delay_1[v in job_indices, ω in Ω], R[v, ω] >= ε[v, ω]) |
| 53 | + @constraint( |
| 54 | + model, |
| 55 | + R_delay_2[v in job_indices, ω in Ω], |
| 56 | + R[v, ω] >= |
| 57 | + ε[v, ω] + sum( |
| 58 | + yR[u, v, ω] - y[u, v] * slacks[u, v][ω] for u in nodes if has_edge(graph, u, v) |
| 59 | + ) |
| 60 | + ) |
| 61 | + |
| 62 | + # Mc Cormick linearization constraints |
| 63 | + @constraint( |
| 64 | + model, |
| 65 | + R_McCormick_1[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)], |
| 66 | + yR[u, v, ω] >= R[u, ω] + Rmax * (y[u, v] - 1) |
| 67 | + ) |
| 68 | + @constraint( |
| 69 | + model, |
| 70 | + R_McCormick_2[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)], |
| 71 | + yR[u, v, ω] <= Rmax * y[u, v] |
| 72 | + ) |
| 73 | + |
| 74 | + # Solve model |
| 75 | + optimize!(model) |
| 76 | + solution = value.(y) |
| 77 | + |
| 78 | + sol = solution_from_JuMP_array(solution, graph) |
| 79 | + return objective_value(model), sol |
| 80 | +end |
| 81 | + |
| 82 | +""" |
| 83 | +$TYPEDSIGNATURES |
| 84 | +
|
| 85 | +Returns the optimal solution of the Stochastic VSP instance, by solving the associated compact quadratic MIP. |
| 86 | +Note: If you have Gurobi, use `grb_model` as `model_builder` instead of `highs_model`. |
| 87 | +
|
| 88 | +!!! warning |
| 89 | + You need to use a solver that supports quadratic constraints to use this method. |
| 90 | +""" |
| 91 | +function compact_mip(instance::Instance; model_builder=scip_model, silent=true) |
| 92 | + (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance |
| 93 | + nb_nodes = nv(graph) |
| 94 | + job_indices = 2:(nb_nodes - 1) |
| 95 | + nodes = 1:nb_nodes |
| 96 | + |
| 97 | + # Pre-processing |
| 98 | + ε = intrinsic_delays |
| 99 | + nb_scenarios = size(ε, 2) |
| 100 | + Ω = 1:nb_scenarios |
| 101 | + |
| 102 | + # Model definition |
| 103 | + model = model_builder() |
| 104 | + silent && set_silent(model) |
| 105 | + |
| 106 | + # Variables and objective function |
| 107 | + @variable(model, y[u in nodes, v in nodes; has_edge(graph, u, v)], Bin) |
| 108 | + @variable(model, R[v in nodes, ω in Ω] >= 0) # propagated delay of job v |
| 109 | + @variable(model, yR[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)] >= 0) # yR[u, v] = y[u, v] * R[u, ω] |
| 110 | + |
| 111 | + @objective( |
| 112 | + model, |
| 113 | + Min, |
| 114 | + delay_cost * sum(sum(R[v, ω] for v in job_indices) for ω in Ω) / nb_scenarios # average total delay |
| 115 | + + |
| 116 | + vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles |
| 117 | + ) |
| 118 | + |
| 119 | + # Flow contraints |
| 120 | + @constraint( |
| 121 | + model, |
| 122 | + flow[i in job_indices], |
| 123 | + sum(y[j, i] for j in inneighbors(graph, i)) == |
| 124 | + sum(y[i, j] for j in outneighbors(graph, i)) |
| 125 | + ) |
| 126 | + @constraint( |
| 127 | + model, |
| 128 | + unit_demand[i in job_indices], |
| 129 | + sum(y[j, i] for j in inneighbors(graph, i)) == 1 |
| 130 | + ) |
| 131 | + |
| 132 | + # Delay propagation constraints |
| 133 | + @constraint(model, [ω in Ω], R[1, ω] == ε[1, ω]) |
| 134 | + @constraint(model, R_delay_1[v in job_indices, ω in Ω], R[v, ω] >= ε[v, ω]) |
| 135 | + @constraint( |
| 136 | + model, |
| 137 | + R_delay_2[v in job_indices, ω in Ω], |
| 138 | + R[v, ω] >= |
| 139 | + ε[v, ω] + |
| 140 | + sum(y[u, v] * (R[u, ω] - slacks[u, v][ω]) for u in nodes if has_edge(graph, u, v)) |
| 141 | + ) |
| 142 | + |
| 143 | + # Solve model |
| 144 | + optimize!(model) |
| 145 | + solution = value.(y) |
| 146 | + |
| 147 | + sol = solution_from_JuMP_array(solution, graph) |
| 148 | + return objective_value(model), sol |
| 149 | +end |
0 commit comments