From e30cb36adbf97f4f81c500dab52f2c4603e235ac Mon Sep 17 00:00:00 2001 From: Natacha Javerzat Date: Fri, 12 May 2023 16:03:43 +0200 Subject: [PATCH] re-organize demo, fix typos --- docs/src/start/advanced-demo.jl | 218 ++++++++++++++++---------------- 1 file changed, 112 insertions(+), 106 deletions(-) diff --git a/docs/src/start/advanced-demo.jl b/docs/src/start/advanced-demo.jl index ddf004af3..6816f92ca 100644 --- a/docs/src/start/advanced-demo.jl +++ b/docs/src/start/advanced-demo.jl @@ -174,7 +174,7 @@ function create_model(optimizer, pricing_algorithms) subproblemrepresentative.(z, Ref(subproblems)) return model, x, y, z, cov -end +end; # Note that contrary to the direct model, we don't have to add constraints to ensure the # consistency of the routes because we solve our subproblems using a pricing callback. @@ -218,7 +218,7 @@ function route_original_cost(costs, route::Route) route_cost += costs[path[i], path[i+1]] end return route_cost -end +end; # Since the cost of the route depends only on the arcs and there is no resource consumption # constraints on the route or master constraints that may have an effect of the customer visit order, @@ -244,7 +244,7 @@ function best_visit_order(costs, cust_subset, facility_id) tmp = argmin([c for (_, c) in routes_costs]) (best_order, _) = routes_costs[tmp] return best_order -end +end; # We are now able to compute the best route for all the possible customers subsets, # given a facility id: @@ -265,7 +265,7 @@ function best_route_forall_cust_subsets(costs, customers, facility_id, max_size) push!(best_routes, route_s) end return best_routes -end +end; # We store all the information given by this pre-processing phase in a dictionary. # To each facility id we match a vector of routes that are the best visiting orders @@ -287,7 +287,7 @@ function x_contribution(route::Route, j::Int, x_red_costs) x += x_red_costs["x_$(i)_$(j)"] end return x -end +end; function z_contribution(route::Route, z_red_costs) z = 0.0 @@ -297,7 +297,7 @@ function z_contribution(route::Route, z_red_costs) z += z_red_costs["z_$(current_position)_$(next_position)"] end return z -end +end; # We are now able to write our pricing callback: @@ -345,7 +345,7 @@ function pricing_callback(cbdata) ## bound in column generation. MOI.submit(model, BlockDecomposition.PricingDualBound(cbdata), sol_cost) ## optimal solution -end +end; # Create the model: model, x, y, z, _ = create_model(coluna, pricing_callback); @@ -354,8 +354,6 @@ model, x, y, z, _ = create_model(coluna, pricing_callback); JuMP.optimize!(model) -# TODO: display "raw" decomp model output and comment, transition to next section - # ### Strengthen with robust cuts (valid inequalities) # We introduce of first type of classic valid inequalities that tries to improve the @@ -403,14 +401,14 @@ function valid_inequalities_callback(cbdata) constr = JuMP.@build_constraint(x[ineq.customer_id, ineq.facility_id] <= y[ineq.facility_id]) MOI.submit(model, MOI.UserCut(cbdata), constr) end -end +end; # We re-declare the model and optimize it with these valid inequalites: model, x, y, z, _ = create_model(coluna, pricing_callback); MOI.set(model, MOI.UserCutCallback(), valid_inequalities_callback); JuMP.optimize!(model) -# TODO: comment on the improvement of the dual bound + # ### Strengthen with non-robust cuts (rank-one cuts) @@ -496,7 +494,7 @@ BlockDecomposition.customconstrs!(model, [CoverConstrData, R1cCutData]); # contribution (i.e. linear expression of subproblem variable valuations) and the non-robust # contribution retuned by the `computecoeff` method. -# So we define this first method that Coluna call sget the coefficients of the columns +# So we define this first method that Coluna calls to get the coefficients of the columns # in the R1C: function Coluna.MathProg.computecoeff( ::Coluna.MathProg.Variable, var_custom_data::R1cVarData, @@ -555,11 +553,11 @@ function r1c_callback(cbdata) ) end end -end +end; # You should find disturbing the way we add the non-robust valid inequalities because we # literally add the constraint `0 <= 1` to the model. -# This is because you don't have access to the columns from the seperation callbacl and how +# This is because you don't have access to the columns from the separation k and how # Coluna computes the coefficient of the columns in the constraints. # You can only express the robust-part of the inequality. The non-robust part is calculated # in the `compute_coeff` method. @@ -567,7 +565,7 @@ end # The last thing we need to do to complete the implementation of R1Cs is to update our # pricing callback. Unlike valid inequalities, R1Cs are not expressed directly with the # subproblem variables. -# Thus, their contribution to the redcued cost of a column is not captured by the reduced cost +# Thus, their contribution to the reduced cost of a column is not captured by the reduced cost # of subproblem variables. # We must therefore take this contribution into account "manually". @@ -581,7 +579,7 @@ function r1c_contrib(route::Route, custduals) end end return cost -end +end; # We re-write our pricing callback to: # - retrieve the dual cost of the R1Cs @@ -677,20 +675,10 @@ function add_nearest_neighbor(route::Route, customers, costs) push!(route.path, nearest) route.length += 1 end -end - -# Then we define our inexact pricing callback: -function approx_pricing(cbdata) - - j = BlockDecomposition.indice(BlockDecomposition.callback_spid(cbdata, model)) - z_red_costs = Dict( - "z_$(u)_$(v)" => BlockDecomposition.callback_reduced_cost(cbdata, z[u, v]) for u in locations, v in locations - ) - x_red_costs = Dict( - "x_$(i)_$(j)" => BlockDecomposition.callback_reduced_cost(cbdata, x[i, j]) for i in customers - ) - +end; +# We then define our heuristic for the enumeration of the routes, the method returns the best route found by he heuristic together with its cost: +function enumeration_heuristic(x_red_costs, z_red_costs, j) ## initialize our "greedy best route" best_route = Route(1, [j]) ## initialize the route's cost to zero @@ -708,7 +696,24 @@ function approx_pricing(cbdata) break end end + return (best_route, current_redcost) +end + +# We can now define our inexact pricing callback: +function approx_pricing(cbdata) + j = BlockDecomposition.indice(BlockDecomposition.callback_spid(cbdata, model)) + z_red_costs = Dict( + "z_$(u)_$(v)" => BlockDecomposition.callback_reduced_cost(cbdata, z[u, v]) for u in locations, v in locations + ) + x_red_costs = Dict( + "x_$(i)_$(j)" => BlockDecomposition.callback_reduced_cost(cbdata, x[i, j]) for i in customers + ) + + ## call the heuristic to elect the "greedy best route": + best_route, sol_cost = enumeration_heuristic(x_red_costs, z_red_costs, j) + + ## build the solution: best_route_arcs = Vector{Tuple{Int, Int}}() for i in 1:(best_route.length - 1) push!(best_route_arcs, (best_route.path[i], best_route.path[i+1])) @@ -719,8 +724,6 @@ function approx_pricing(cbdata) x_vars = [x[i, j] for i in best_route_customers] sol_vars = vcat(z_vars, x_vars) sol_vals = ones(Float64, length(z_vars) + length(x_vars)) - ## take the eventual rank-one cuts contribution into account - sol_cost = current_redcost MOI.submit(model, BlockDecomposition.PricingSolution(cbdata), sol_cost, sol_vars, sol_vals) ## as the procedure is inexact, no dual bound can be computed, we set it to -Inf @@ -748,84 +751,9 @@ for i in customers end # Optimize: - JuMP.optimize!(model) -# ## Benders decomposition - -# In this last section, we show you how you solve the linear relaxation of the problem using Benders. - -# At the first stage, we chose a subset of facilities to open. -# At the second stage, we solve the linear rekaxation of the routing problem for each facility that is open. -# Since we play with a toy instance and we can enumerate all the routes, the second-level will -# return select the routes to use. - -# We introduce the variables. -# Let `y[j]` equals 1 if the facility `j` is open and 0 otherwise. -# Let `λ[j,k]` equals 1 if route `k` starting from facility `j` is selected and 0 otherwise. - -# Since there is only one subproblem at the second stage, we introduce a fake axis that contains -# only one element. -# In the case of Benders, when the second-stage is decomposable, you can imagine the axis as the set of scenarios for instance. - -fake = 1 -@axis(axis, collect(fake:fake)) - -coluna = JuMP.optimizer_with_attributes( - Coluna.Optimizer, - "params" => Coluna.Params( - - solver = Coluna.Algorithm.BendersCutGeneration() - ), - "default_optimizer" => GLPK.Optimizer -) - -model = BlockModel(coluna); - - -# We introduce auxiliary structures to make the code readable. - -## routes covering customer i from facility j. -covering_routes = Dict( - (j, i) => findall(r -> (i in r.path), routes_per_facility[j]) for i in customers, j in facilities -) -## routes costs from facility j. -routes_costs = Dict( - j => [route_original_cost(arc_costs, r) for r in routes_per_facility[j]] for j in facilities -) - -# We declare the variables. -@variable(model, 0 <= y[j in facilities] <= 1) ## 1st stage -@variable(model, 0 <= λ[f in axis, j in facilities, k in 1:length(routes_per_facility[j])] <= 1) ## 2nd stage - -# We declare the constraints. - -## Linking constraints -@constraint(model, open[fake in axis, j in facilities, k in 1:length(routes_per_facility[j])], - y[j] >= λ[fake, j, k]) - -## Second-stage constraints -@constraint(model, cover[fake in axis, i in customers], - sum(λ[fake, j, k] for j in facilities, k in covering_routes[(j,i)]) >= 1) - -## Second-stage constraints -@constraint(model, limit_nb_routes[fake in axis, j in facilities], - sum(λ[fake, j, q] for q in 1:length(routes_per_facility[j])) <= nb_routes_per_facility -) - -## First-stage constraint -@constraint(model, min_opening, - sum(y[j] for j in facilities) >= 1) - -@objective(model, Min, - sum(facilities_fixed_costs[j] * y[j] for j in facilities) + - sum(routes_costs[j][k] * λ[fake, j, k] for j in facilities, k in 1:length(routes_per_facility[j]))) - -# We perform the decomposition over the axis and we optimize the problem. -@benders_decomposition(model, dec, axis) -JuMP.optimize!(model) - # ## Example of comparison of the dual bounds obtained on a larger instance. @@ -922,3 +850,81 @@ attach_data(model, cov) MOI.set(model, MOI.UserCutCallback(), cuts_callback); # dual bound found after optimization = 1600.63 + + +# ## Benders decomposition + +# In this last section, we show you how you solve the linear relaxation of the problem using Benders. + +# At the first stage, we chose a subset of facilities to open. +# At the second stage, we solve the linear relaxation of the routing problem for each facility that is open. +# Since we play with a toy instance and we can enumerate all the routes, the second-level will +# return select the routes to use. + +# We introduce the variables. +# Let `y[j]` equals 1 if the facility `j` is open and 0 otherwise. +# Let `λ[j,k]` equals 1 if route `k` starting from facility `j` is selected and 0 otherwise. + +# Since there is only one subproblem at the second stage, we introduce a fake axis that contains +# only one element. +# In the case of Benders, when the second-stage is decomposable, you can imagine the axis as the set of scenarios for instance. + +fake = 1 +@axis(axis, collect(fake:fake)) + +coluna = JuMP.optimizer_with_attributes( + Coluna.Optimizer, + "params" => Coluna.Params( + + solver = Coluna.Algorithm.BendersCutGeneration() + ), + "default_optimizer" => GLPK.Optimizer +) + +model = BlockModel(coluna); + + +# We introduce auxiliary structures to make the code readable. + +## routes covering customer i from facility j. +covering_routes = Dict( + (j, i) => findall(r -> (i in r.path), routes_per_facility[j]) for i in customers, j in facilities +) +## routes costs from facility j. +routes_costs = Dict( + j => [route_original_cost(arc_costs, r) for r in routes_per_facility[j]] for j in facilities +) + +# We declare the variables. +@variable(model, 0 <= y[j in facilities] <= 1) ## 1st stage +@variable(model, 0 <= λ[f in axis, j in facilities, k in 1:length(routes_per_facility[j])] <= 1) ## 2nd stage + +# We declare the constraints. + +## Linking constraints +@constraint(model, open[fake in axis, j in facilities, k in 1:length(routes_per_facility[j])], + y[j] >= λ[fake, j, k]) + +## Second-stage constraints +@constraint(model, cover[fake in axis, i in customers], + sum(λ[fake, j, k] for j in facilities, k in covering_routes[(j,i)]) >= 1) + +## Second-stage constraints +@constraint(model, limit_nb_routes[fake in axis, j in facilities], + sum(λ[fake, j, q] for q in 1:length(routes_per_facility[j])) <= nb_routes_per_facility +) + +## First-stage constraint +@constraint(model, min_opening, + sum(y[j] for j in facilities) >= 1) + +@objective(model, Min, + sum(facilities_fixed_costs[j] * y[j] for j in facilities) + + sum(routes_costs[j][k] * λ[fake, j, k] for j in facilities, k in 1:length(routes_per_facility[j]))) + +# We perform the decomposition over the axis and we optimize the problem. +@benders_decomposition(model, dec, axis) +JuMP.optimize!(model) + + +