From e2c23fe3ac01030d1e3fc9b69f15f7eddb1fe32f Mon Sep 17 00:00:00 2001 From: Natacha Javerzat Date: Wed, 10 May 2023 14:31:43 +0200 Subject: [PATCH 1/6] init advanced demo --- docs/src/start/advanced-demo.jl | 547 ++++++++++++++++++++++++++++++++ 1 file changed, 547 insertions(+) create mode 100644 docs/src/start/advanced-demo.jl diff --git a/docs/src/start/advanced-demo.jl b/docs/src/start/advanced-demo.jl new file mode 100644 index 000000000..1f97fceb3 --- /dev/null +++ b/docs/src/start/advanced-demo.jl @@ -0,0 +1,547 @@ +# # Location Routing + +# We demonstrate the main features of Coluna on a variant of the Location Routing problem. +# In our variant, we consider customers that must be delivered by facilities. For simplicity, we assume that any customer's request can be fulfilled by any facility. To deliver a customer, we must open at least one route that passes through the customer's location. +# A route is defined as a vector of locations that satisfies the following rules: + +# - any route must start from a facility location +# - every route has a maximum length, i.e. the number of visited locations cannot exceed a fixed constant `nb_positions` + +# - the routes are said to be open, i.e. finish at last visited customer. +# - each route has a cost defined as the sum of the costs of its arcs. + +# Our objective is to select the best routes to open in order to deliver all customers while minimizing the opening costs of the routes. This cost depends on the original costs of the routes and the opening costs of the facilities. + + +# In this tutorial, we work on a small instance with only 2 facilities and 7 customers. The maximum length of a route is fixed to 6. +nb_positions = 6 +facilities_fixed_costs = [1200, 1500] +facilities = [1, 2] +customers = [3, 4, 5, 6, 7, 8, 9] +arc_costs = +[ + 0.0 253.241 253.924 253.706 353.612 373.716 319.347 246.488 342.445 ; + 253.241 0.0 212.017 161.689 270.646 268.338 178.325 166.935 231.827 ; + 253.924 212.017 0.0 142.475 233.867 238.415 182.768 170.346 216.402 ; + 253.706 161.689 142.475 0.0 285.9 288.07 226.365 156.479 294.774 ; + 353.612 270.646 233.867 285.9 0.0 420.728 304.439 248.598 391.131 ; + 373.716 268.338 238.415 288.07 420.728 0.0 324.366 295.006 381.664 ; + 319.347 178.325 182.768 226.365 304.439 324.366 0.0 224.823 307.48 ; + 246.488 166.935 170.346 156.479 248.598 295.006 224.823 0.0 213.504 ; + 342.445 231.827 216.402 294.774 391.131 381.664 307.48 213.504 0.0 + ] + +locations = vcat(facilities, customers) +nb_customers = length(customers) +nb_facilities = length(facilities) +positions = 1:length(nb_positions) + + + +# Let's define the model that will solve our problem; first we need to load some packages and dependencies: + +using JuMP, HiGHS, GLPK +using BlockDecomposition, Coluna + + +# We want to set an upper bound `nb_routes_per_locations` on the number of routes starting from a facility. +# This limit is calculated as follows: + +# We compute the minimum number of routes needed to visit all customers: +nb_routes = Int(ceil(nb_customers / nb_positions)) +# We define the upper bound `nb_routes_per_locations`: +nb_routes_per_locations = min(Int(ceil(nb_routes / nb_facilities)) * 2, nb_routes) + +# ## Direct model + +# First, we solve the problem by a direct approach, using the HiGHS solver. We start by initializing the solver: + +direct = JuMP.direct_model(HiGHS.Optimizer()) + + +# and we declare 3 types of binary variables: + +@variable(model, y[j in facilities], Bin) +@variable(model, z[u in locations, v in locations], Bin) +@variable(model, x[i in customers, j in facilities, k in routes_per_locations, p in positions], Bin) + +# variables `y_j` indicate the opening status of a given facility ; if `y_j = 1` then facility `j` is open, otherwise `y_j = 0` +# variables `z_(u,v)` indicate which arcs are used ; if `z_(u,v) = 1` then there is an arc between location `u` and location `v`, otherwise `z_(u,v) = 0` +# variables `x_(i, j, k, p)` are used to express cover constraints and to ensure the consistency of the routes ; `x_(i, j, k, p) = 1` if customer `i` is delivered from facility `j` at the position `p` of route `k`, otherwise `x_(i, j, k, p) = 0`. + +# Now, we add constraints to our model: + +# "each customer is visited by at least one route" +@constraint(model, cov[i in customers], + sum(x[i, j, k, p] for j in facilities, k in routes_per_locations, p in positions) >= 1) +# "for any route from any facility, its length does not exceed the fixed maximum length `nb_positions`" +@constraint(model, cardinality[j in facilities, k in routes_per_locations], + sum(x[i, j, k, p] for i in customers, p in positions) <= nb_positions * y[j]) +# "only one customer can be delivered by a given route at a given position" +@constraint(model, assign_setup[p in positions, j in facilities, k in routes_per_locations], + sum(x[i, j, k, p] for i in customers) <= y[j]) +# "a customer can only be delivered at position `p > 1` of a given route if there is a customer delivered at position `p-1` of the same route" +@constraint(model, open_route[j in facilities, k in routes_per_locations, p in positions; p > 1], + sum(x[i, j, k, p] for i in customers) <= sum(x[i, j, k, p-1] for i in customers)) +# "there is an arc between two customers whose demand is satisfied by the same route at consecutive positions" +@constraint(model, route_arc[i in customers, l in customers, indi in 1:nb_customers, indl in 1:nb_customers, j in facilities, k in routes_per_locations, p in positions; p > 1 && i != l && indi != indl], + z[customers[indi], customers[indl]] >= x[l, j, k, p] + x[i, j, k, p-1] - 1) +# "there is an arc between the facility `j` and the first customer visited by the route `k` from facility `j`" +@constraint(model, start_arc[i in customers, indi in 1:nb_customers, j in facilities, k in routes_per_locations], + z[facilities[j], customers[indi]] >= x[i, j, k, 1]) + + +# We set the objective function: +@objective(model, Min, + sum(arc_costs[u, v] * z[u, v] for u in locations, v in locations) + + + sum(facility_fixed_costs[j] * y[j] for j in facilities)) + +# ##TODO: run direct model + +# ## Decomp model + +# We can exploit the structure of the problem by generating routes starting from each facility. +# The most immediate decomposition is to consider each route traveled by a vehicle as a sub-problem. However, at a given facility, vehicles are identical and therefore any vehicle can travel on any route. So we have several identical subproblems at each facility. + +# We declare the axis: + +axis = collect(facilities) +@axis(Base_axis, axis) +@show Base_axis + +# and we set up the solver: + +# ##TODO: clean +coluna = optimizer_with_attributes( + Coluna.Optimizer, + "params" => Coluna.Params( + solver=Coluna.Algorithm.TreeSearchAlgorithm( + maxnumnodes = 1, + conqueralg = Coluna.ColCutGenConquer( + primal_heuristics = [ + #Coluna.ParameterizedHeuristic( + # Diva.Diving(), + # 1.0, + # 1.0, + # 1, + # 1, + # "Diving" + #) + ] + ) + ) # default branch-cut-and-price + ), + "default_optimizer" => GLPK.Optimizer # GLPK for the master & the subproblems +) + +decomp = BlockModel(coluna) + +# Let's declare the variables. We distinct the master variables `y` + +@variable(decomp, y[j in facilities], Bin) + +# from the sub-problem variables `x` and `z`: + +@variable(decomp, x[i in customers, j in Base_axis] <= 1) +@variable(decomp, z[u in locations, v in locations], Bin) + +# The information carried by the `x` variables may seem redundant with that of the `z` variables. +# The `x` variables are in fact introduced only in order to separate a family of robust valid inequalities. + +# We now declare our problem's constraints: +# The cover constraints are expressed w.r.t. the `x` variables: +@constraint(decomp, cov[i in customers], + sum(x[i, j] for j in facilities) >= 1) + +# We add a constraint to express that if a facility is not opened, there can be no arc between this facility and a customer. +@constraint(decomp, open_facility[j in facilities], + sum(z[j, i] for i in customers) <= y[j] * nb_routes_per_locations) + +# We do not need to add constraints ensuring the consistency of the routes because we solve our subproblems by pricing. It will therefore be the responsibility of the pricing callback to create consistent routes. ##TODO link with additional constraints of Direct model + +# We set the objective function: + +@objective(decomp, Min, +sum(arc_costs[u, v] * z[u, v] for u in locations, v in locations) ++ +sum(facilities_fixed_costs[j] * y[j] for j in facilities)) + +# ## Pricing callback + +# Each sub-problem could be solved by a MIP, provided the right sub-problem constraints are added. Here, we propose a resolution by enumeration within a pricing callback. + +# The general idea of enumeration is very simple: we enumerate the possible routes from a facility and keep the one with the lowest reduced cost, i.e. the one that improves the current solution the most. +# Enumerating all possible routes is very expensive. We improve the pricing efficiency a bit by pre-processing, for a given subset of customers and a given facility, the best order to visit the customers of the subset. This order depends only on the original cost of the arcs, so we need a method to compute it: + +struct EnumeratedRoute + length::Int + path::Vector{Int} +end + +function route_original_cost(costs, enroute::EnumeratedRoute) + route_cost = 0.0 + path = enroute.path + path_length = enroute.length + for i in 1:(path_length-1) + route_cost += costs[path[i], path[i+1]] + end + return route_cost +end + +function best_visit_order(costs, cust_subset, facility_id) + ## generate all the possible visit orders + set_size = size(cust_subset)[1] + all_paths = collect(multiset_permutations(cust_subset, set_size)) + all_enroutes = Vector{EnumeratedRoute}() + for path in all_paths + ## add the first index i.e. the facility id + enpath = vcat([facility_id], path) + ## length of the route = 1 + number of visited customers + enroute = EnumeratedRoute(set_size + 1, enpath) + push!(all_enroutes, enroute) + end + ## compute each route original cost + enroutes_costs = map(r -> + (r, route_original_cost(costs, r)), all_enroutes ) + ## keep the best visit order + tmp = argmin([c for (_, c) in enroutes_costs]) + (best_order, _) = enroutes_costs[tmp] + return best_order +end + +# We are now able to compute the best route for all the possible customers subsets, given a facility id: + +using Combinatorics + +function best_route_forall_cust_subsets(costs, customers, facility_id, max_size) + best_routes = Vector{EnumeratedRoute}() + all_subsets = Vector{Vector{Int}}() + for subset_size in 1:max_size + subsets = collect(combinations(customers, subset_size)) + for s in subsets + push!(all_subsets, s) + end + end + for s in all_subsets + route_s = best_visit_order(costs, s, facility_id) + push!(best_routes, route_s) + end + return best_routes +end + +# We store all the information given by the pre-computation in a dictionary. To each facility id we match a vector of routes that are the best visiting orders for each possible subset of customers. + +routes_per_facility = Dict( + j => best_route_forall_cust_subsets(arc_costs, customers, j, nb_positions) for j in facilities + ) + + +# We must also declare methods to calculate the contribution to the reduced cost of the two types of subproblem variables, `x` and `z`: + +function x_contribution(enroute::EnumeratedRoute, j::Int, x_red_costs) + x = 0.0 + visited_customers = enroute.path[2:enroute.length] + for i in visited_customers + x += x_red_costs["x_$(i)_$(j)"] + end + return x +end + +function z_contribution(enroute::EnumeratedRoute, z_red_costs) + z = 0.0 + for i in 1:(enroute.length-1) + current_position = enroute.path[i] + next_position = enroute.path[i+1] + z += z_red_costs["z_$(current_position)_$(next_position)"] + end + return z +end + +# We are now able to write our pricing callback: + +function my_pricing_callback(cbdata) + ## get the id of the facility + j = BlockDecomposition.indice(BlockDecomposition.callback_spid(cbdata, decomp)) + + ## retrieve variables reduced costs + 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 + ) + + ## keep route with minimum reduced cost. + red_costs_j = map(r -> ( + r, + x_contribution(r, j, x_red_costs) + z_contribution(r, z_red_costs) # the reduced cost of a route is the sum of the contribution of the variables + ), routes_per_facility[j] + ) + min_index = argmin([x for (_,x) in red_costs_j]) + (best_route, min_reduced_cost) = red_costs_j[min_index] + + ## Create the solution (send only variables with non-zero values) + + ## retrieve the route's arcs + 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])) + end + best_route_customers = best_route.path[2:best_route.length] + z_vars = [z[u, v] for (u,v) in best_route_arcs] + 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)) + sol_cost = min_reduced_cost + + ## Submit the solution of the subproblem to Coluna + MOI.submit(decomp, BlockDecomposition.PricingSolution(cbdata), sol_cost, sol_vars, sol_vals) + + ## Submit the dual bound to the solution of the subproblem + ## This bound is used to compute the contribution of the subproblem to the lagrangian + ## bound in column generation. + MOI.submit(decomp, BlockDecomposition.PricingDualBound(cbdata), sol_cost) # optimal solution + +end + + +# Set decomposition +#@dantzig_wolfe_decomposition(decomp, dec, Base_axis) +#subproblems = BlockDecomposition.getsubproblems(dec) +#specify!.(subproblems, lower_multiplicity=0, upper_multiplicity=nb_routes_per_locations, solver = my_pricing_callback) #set solver = nothing to test without my_pricing_callback +#subproblemrepresentative.(z, Ref(subproblems)) +##Solve the model +#JuMP.optimize!(decomp) + +# ## TODO: display "raw" decomp model output and comment, transition to next section + +# ## Strengthen with robust cuts (valid inequalities) + +# We want to add some constraints of the form `x_i_j <= y_j ∀i ∈ customers, ∀j ∈ facilities` in order to try to "improve" the integrality of `y_j`. However, there are `nb_customers x nb_facilities` such constraints. Instead of adding all of them, we want to consider a constraint if and only if it is violated, i.e. if the solution found does not respect the constraint (this is a cut generation approach). + +# We declare a structure representing an inequality `x_i_j <= y_j`for a fixed facility `j` and a fixed customer `i`: +struct OpenFacilityInequality + facility_id::Int + customer_id::Int +end + +# We write our valid inequalities callback: +function valid_inequalities_callback(cbdata) + ## Get variables valuations, store them into dictionaries + x_vals = Dict( + "x_$(i)_$(j)" => BlockDecomposition.callback_value(cbdata, x[i, j]) for i in customers, j in facilities + ) + y_vals = Dict( + "y_$(j)" => BlockDecomposition.callback_value(cbdata, y[j]) for j in facilities + ) + + ## Separate the valid inequalities i.e. retrieve the inequalities that are violated by the current solution. + inequalities = Vector{OpenFacilityInequality}() + + for j in facilities + for i in customers + x_i_j = x_vals["x_$(i)_$(j)"] + y_j = y_vals["y_$(j)"] + if x_i_j > y_j + push!(inequalities, OpenFacilityInequality(j, i)) + end + end + end + + ## Add the valid inequalities to the model. + for ineq in inequalities + constr = JuMP.@build_constraint(x[ineq.customer_id, ineq.facility_id] <= y[ineq.facility_id]) + MOI.submit(decomp, MOI.UserCut(cbdata), constr) + end +end + +MOI.set(decomp, MOI.UserCutCallback(), valid_inequalities_callback); + + + +#MOI.set(decomp, MOI.UserCutCallback(), valid_inequalities_callback); +#@dantzig_wolfe_decomposition(decomp, dec, Base_axis) +#subproblems = BlockDecomposition.getsubproblems(dec) +#specify!.(subproblems, lower_multiplicity=0, upper_multiplicity=nb_routes_per_locations, solver = my_pricing_callback) #set solver = nothing to test without my_pricing_callback +#subproblemrepresentative.(z, Ref(subproblems)) +#JuMP.optimize!(decomp) + +# ## TODO: comment on the improvement of the dual bound, fix display of both raw decomp and valid ineq decomp model + + + +# ## Strengthen with non-robust cuts (rank-one cuts) + +# ##TODO: describe the goal by looking at lambdas values -> highlight that the aim is to drive them towards integrality + +# Here, we implement special types of cuts called "rank-one cuts" (R1C). These cuts are non-robust in the sense that they can not be expressed only with the original variables of the model. In particular, they have to be expressed with the master columns variables `λ_k`. +# R1Cs are obtained by applying the Chvátal-Gomory procedure once, hence their name, on cover constraints. We must therefore be able to differentiate the cover constraints from the other constraints of the model. To do this, we exploit an advantage of Coluna that allows us to attach custom data to the constraints and variables of our model: + +# We create a special custom data with the only information we need to characterize our cover constraints: the customer id that corresponds to this constraint. +struct CoverConstrData <: BlockDecomposition.AbstractCustomData + customer::Int +end +# We declare our custom data to Coluna +BlockDecomposition.customconstrs!(decomp, CoverConstrData); +# And we attach one custom data to each cover constraint +for i in customers + customdata!(cov[i], CoverConstrData(i)) +end + + +# The rank-one cuts we are going to add are of the form: +# `sum(c_k λ_k) <= 1.0` +# for a fixed subset `r1c_cov_constrs` of cover constraints of size 3, with `λ_k` the master columns variables and c_k` s.t. +# `c_k = ⌊ 1/2 x |r1c_locations ∩ r1c_cov_constrs| ⌋` +# with `r1c_locations` the current solution (route) that corresponds to `λ_k`. +# e.g. if we consider cover constraints cov[3], cov[6] and cov[8] in our cut, then the route 1-4-6-7 gives a zero coefficient while the route 1-4-6-3 gives a coefficient equal to one. + +# But a problem arises: how to get the current solution `r1c_locations` that corresponds to a given `λ_k` ? To handle that difficulty, we use once again the custom data trick: + +# Each `λ_k` is associated to a `R1cVarData` structure that carries the current solution. +struct R1cVarData <: BlockDecomposition.AbstractCustomData + visited_locations::Vector{Int} +end + +# The rank-one cuts are associated with `R1cCutData` structures indicating which cover constraints are taken into account in the cut. +struct R1cCutData <: BlockDecomposition.AbstractCustomData + cov_constrs::Vector{Int} +end + +# We declare our custom data to Coluna: +BlockDecomposition.customvars!(decomp, R1cVarData) +BlockDecomposition.customconstrs!(decomp, [CoverConstrData, R1cCutData]); + +# This method is called by Coluna to compute the coefficients of the `λ_k` in the cuts: +function Coluna.MathProg.computecoeff( + ::Coluna.MathProg.Variable, var_custom_data::R1cVarData, + ::Coluna.MathProg.Constraint, constr_custom_data::R1cCutData +) + return floor(1/2 * length(var_custom_data.visited_locations ∩ constr_custom_data.cov_constrs)) +end + +# TODO: fix necessity to write computecoeff for cover constr or explain trick +function Coluna.MathProg.computecoeff( + ::Coluna.MathProg.Variable, ::R1cVarData, + ::Coluna.MathProg.Constraint, ::CoverConstrData) + return 0 +end + +# We are now able to write our rank-one cut callback: +function r1c_callback(cbdata) + original_sol = cbdata.orig_sol + master = Coluna.MathProg.getmodel(original_sol) + ## retrieve the cover constraints + cov_constrs = Int[] + for constr in values(Coluna.MathProg.getconstrs(master)) + if typeof(constr.custom_data) <: CoverConstrData + push!(cov_constrs, constr.custom_data.customer) + end + end + + ## retrieve the master columns λ + lambdas = Vector{Any}() + for (var_id, val) in original_sol + if Coluna.MathProg.getduty(var_id) <= Coluna.MathProg.MasterCol + push!(lambdas, (val, Coluna.MathProg.getvar(cbdata.form, var_id))) + end + end + + ## separate the valid R1Cs (i.e. those violated by the current solution) + ## for a fixed subset of cover constraints of size 3, iterate on the master columns and check violation: + subsets = collect(combinations(cov_constrs, 3)) + for cov_constr_subset in subsets + violation = 0 + for lambda in lambdas + (val, var) = lambda + if !isnothing(var.custom_data) + coeff = floor(1/2 * length(var.custom_data.visited_locations ∩ cov_constr_subset)) + violation += coeff * val + end + end + if violation > 1 + ## create the constraint and add it to the model, use custom data to keep information about the cut (= the subset of considered cover constraints) + MOI.submit(decomp, + MOI.UserCut(cbdata), + JuMP.ScalarConstraint(JuMP.AffExpr(0.0), MOI.LessThan(1.0)), + R1cCutData(cov_constr_subset) + ) + end + end + +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 model variables. Thus, their cost is not taken into account in the reduced cost calculations. We must therefore add it "manually" in the callback. + +# The contribution of R1Cs to the reduced cost computation is managed by the following method: +function r1c_contrib(enroute::EnumeratedRoute, custduals) + cost=0 + if !isempty(custduals) + for (r1c_cov_constrs, dual) in custduals + coeff = floor(1/2 * length(enroute.path ∩ r1c_cov_constrs)) + cost += coeff*dual + end + end + return cost +end + +# We re-write our pricing callback, with the additional contribution that corresponds to R1Cs cost: +function my_pricing_callback(cbdata) + j = BlockDecomposition.indice(BlockDecomposition.callback_spid(cbdata, decomp)) + 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 + ) + ## Get the dual values of the custom cuts to calculate contributions of + ## non-robust cuts to the cost of the solution: + custduals = Tuple{Vector{Int}, Float64}[] + for (_, constr) in Coluna.MathProg.getconstrs(cbdata.form.parent_formulation) + if typeof(constr.custom_data) == R1cCutData + push!(custduals, ( + constr.custom_data.cov_constrs, + Coluna.MathProg.getcurincval(cbdata.form.parent_formulation, constr) + )) + end + end + + ## Keep route with minimum reduced cost, + ## add variables contribution and also the non-robust cuts contribution + red_costs_j = map(r -> ( + r, + x_contribution(r, j, x_red_costs) + + z_contribution(r, z_red_costs) - #TODO: comment on sign ? + r1c_contrib(r, custduals) + ), routes_per_facility[j] + ) + min_index = argmin([x for (_,x) in red_costs_j]) + (best_route, min_reduced_cost) = red_costs_j[min_index] + + 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])) + end + best_route_customers = best_route.path[2:best_route.length] + z_vars = [z[u, v] for (u,v) in best_route_arcs] + x_vars = [x[i, j] for i in best_route_customers] + sol_vars = vcat(z_vars, x_vars) + sol_vals = vcat([1.0 for _ in z_vars], [1.0 for _ in x_vars]) + sol_cost = min_reduced_cost + + ## Submit the solution of the subproblem to Coluna + ## TODO: comment on custom data here + MOI.submit(decomp, BlockDecomposition.PricingSolution(cbdata), sol_cost, sol_vars, sol_vals, R1cVarData(best_route.path)) + MOI.submit(decomp, BlockDecomposition.PricingDualBound(cbdata), sol_cost) + +end + + +MOI.set(decomp, MOI.UserCutCallback(), r1c_callback); +@dantzig_wolfe_decomposition(decomp, dec, Base_axis) +subproblems = BlockDecomposition.getsubproblems(dec) +specify!.(subproblems, lower_multiplicity=0, upper_multiplicity=nb_routes_per_locations, solver = my_pricing_callback) +subproblemrepresentative.(z, Ref(subproblems)) +JuMP.optimize!(decomp) \ No newline at end of file From 95cdcc5588e613e4faf30ee51c33f03827dfd786 Mon Sep 17 00:00:00 2001 From: Natacha Javerzat Date: Wed, 10 May 2023 15:12:34 +0200 Subject: [PATCH 2/6] fix direct model, aux method to create model --- docs/src/start/advanced-demo.jl | 193 +++++++++++++++++--------------- 1 file changed, 100 insertions(+), 93 deletions(-) diff --git a/docs/src/start/advanced-demo.jl b/docs/src/start/advanced-demo.jl index 1f97fceb3..ab669db86 100644 --- a/docs/src/start/advanced-demo.jl +++ b/docs/src/start/advanced-demo.jl @@ -51,6 +51,7 @@ using BlockDecomposition, Coluna nb_routes = Int(ceil(nb_customers / nb_positions)) # We define the upper bound `nb_routes_per_locations`: nb_routes_per_locations = min(Int(ceil(nb_routes / nb_facilities)) * 2, nb_routes) +routes_per_locations = 1:nb_routes_per_locations # ## Direct model @@ -61,9 +62,9 @@ direct = JuMP.direct_model(HiGHS.Optimizer()) # and we declare 3 types of binary variables: -@variable(model, y[j in facilities], Bin) -@variable(model, z[u in locations, v in locations], Bin) -@variable(model, x[i in customers, j in facilities, k in routes_per_locations, p in positions], Bin) +@variable(direct, y[j in facilities], Bin) +@variable(direct, z[u in locations, v in locations], Bin) +@variable(direct, x[i in customers, j in facilities, k in routes_per_locations, p in positions], Bin) # variables `y_j` indicate the opening status of a given facility ; if `y_j = 1` then facility `j` is open, otherwise `y_j = 0` # variables `z_(u,v)` indicate which arcs are used ; if `z_(u,v) = 1` then there is an arc between location `u` and location `v`, otherwise `z_(u,v) = 0` @@ -72,30 +73,30 @@ direct = JuMP.direct_model(HiGHS.Optimizer()) # Now, we add constraints to our model: # "each customer is visited by at least one route" -@constraint(model, cov[i in customers], +@constraint(direct, cov[i in customers], sum(x[i, j, k, p] for j in facilities, k in routes_per_locations, p in positions) >= 1) # "for any route from any facility, its length does not exceed the fixed maximum length `nb_positions`" -@constraint(model, cardinality[j in facilities, k in routes_per_locations], +@constraint(direct, cardinality[j in facilities, k in routes_per_locations], sum(x[i, j, k, p] for i in customers, p in positions) <= nb_positions * y[j]) # "only one customer can be delivered by a given route at a given position" -@constraint(model, assign_setup[p in positions, j in facilities, k in routes_per_locations], +@constraint(direct, assign_setup[p in positions, j in facilities, k in routes_per_locations], sum(x[i, j, k, p] for i in customers) <= y[j]) # "a customer can only be delivered at position `p > 1` of a given route if there is a customer delivered at position `p-1` of the same route" -@constraint(model, open_route[j in facilities, k in routes_per_locations, p in positions; p > 1], +@constraint(direct, open_route[j in facilities, k in routes_per_locations, p in positions; p > 1], sum(x[i, j, k, p] for i in customers) <= sum(x[i, j, k, p-1] for i in customers)) # "there is an arc between two customers whose demand is satisfied by the same route at consecutive positions" -@constraint(model, route_arc[i in customers, l in customers, indi in 1:nb_customers, indl in 1:nb_customers, j in facilities, k in routes_per_locations, p in positions; p > 1 && i != l && indi != indl], +@constraint(direct, route_arc[i in customers, l in customers, indi in 1:nb_customers, indl in 1:nb_customers, j in facilities, k in routes_per_locations, p in positions; p > 1 && i != l && indi != indl], z[customers[indi], customers[indl]] >= x[l, j, k, p] + x[i, j, k, p-1] - 1) # "there is an arc between the facility `j` and the first customer visited by the route `k` from facility `j`" -@constraint(model, start_arc[i in customers, indi in 1:nb_customers, j in facilities, k in routes_per_locations], +@constraint(direct, start_arc[i in customers, indi in 1:nb_customers, j in facilities, k in routes_per_locations], z[facilities[j], customers[indi]] >= x[i, j, k, 1]) # We set the objective function: -@objective(model, Min, +@objective(direct, Min, sum(arc_costs[u, v] * z[u, v] for u in locations, v in locations) + - sum(facility_fixed_costs[j] * y[j] for j in facilities)) + sum(facilities_fixed_costs[j] * y[j] for j in facilities)) # ##TODO: run direct model @@ -106,10 +107,6 @@ direct = JuMP.direct_model(HiGHS.Optimizer()) # We declare the axis: -axis = collect(facilities) -@axis(Base_axis, axis) -@show Base_axis - # and we set up the solver: # ##TODO: clean @@ -117,55 +114,74 @@ coluna = optimizer_with_attributes( Coluna.Optimizer, "params" => Coluna.Params( solver=Coluna.Algorithm.TreeSearchAlgorithm( - maxnumnodes = 1, + maxnumnodes = 0, conqueralg = Coluna.ColCutGenConquer( primal_heuristics = [ #Coluna.ParameterizedHeuristic( - # Diva.Diving(), - # 1.0, - # 1.0, - # 1, - # 1, - # "Diving" - #) - ] - ) - ) # default branch-cut-and-price - ), - "default_optimizer" => GLPK.Optimizer # GLPK for the master & the subproblems -) - -decomp = BlockModel(coluna) - -# Let's declare the variables. We distinct the master variables `y` - -@variable(decomp, y[j in facilities], Bin) - -# from the sub-problem variables `x` and `z`: - -@variable(decomp, x[i in customers, j in Base_axis] <= 1) -@variable(decomp, z[u in locations, v in locations], Bin) - -# The information carried by the `x` variables may seem redundant with that of the `z` variables. -# The `x` variables are in fact introduced only in order to separate a family of robust valid inequalities. - -# We now declare our problem's constraints: -# The cover constraints are expressed w.r.t. the `x` variables: -@constraint(decomp, cov[i in customers], - sum(x[i, j] for j in facilities) >= 1) + # Diva.Diving(), + # 1.0, + # 1.0, + # 1, + # 1, + # "Diving" + #) + ] + ) + ) # default branch-cut-and-price + ), + "default_optimizer" => GLPK.Optimizer # GLPK for the master & the subproblems + ) + + + +# The following method creates the model according to the decomposition described: +function create_model() + + axis = collect(facilities) + @axis(Base_axis, axis) + @show Base_axis + + model = BlockModel(coluna) + + ## Let's declare the variables. We distinct the master variables `y` + + @variable(model, y[j in facilities], Bin) + + ## from the sub-problem variables `x` and `z`: + + @variable(model, x[i in customers, j in Base_axis] <= 1) + @variable(model, z[u in locations, v in locations], Bin) + + ## The information carried by the `x` variables may seem redundant with that of the `z` variables. + ## The `x` variables are in fact introduced only in order to separate a family of robust valid inequalities. + + ## We now declare our problem's constraints: + ## The cover constraints are expressed w.r.t. the `x` variables: + @constraint(model, cov[i in customers], + sum(x[i, j] for j in facilities) >= 1) + + ## We add a constraint to express that if a facility is not opened, there can be no arc between this facility and a customer. + @constraint(model, open_facility[j in facilities], + sum(z[j, i] for i in customers) <= y[j] * nb_routes_per_locations) + + ## We do not need to add constraints ensuring the consistency of the routes because we solve our subproblems by pricing. It will therefore be the responsibility of the pricing callback to create consistent routes. ##TODO link with additional constraints of Direct model + + ## We set the objective function: + + @objective(model, Min, + sum(arc_costs[u, v] * z[u, v] for u in locations, v in locations) + + + sum(facilities_fixed_costs[j] * y[j] for j in facilities)) -# We add a constraint to express that if a facility is not opened, there can be no arc between this facility and a customer. -@constraint(decomp, open_facility[j in facilities], - sum(z[j, i] for i in customers) <= y[j] * nb_routes_per_locations) + @dantzig_wolfe_decomposition(model, dec, Base_axis) + subproblems = BlockDecomposition.getsubproblems(dec) + specify!.(subproblems, lower_multiplicity=0, upper_multiplicity=nb_routes_per_locations, solver=my_pricing_callback) + subproblemrepresentative.(z, Ref(subproblems)) -# We do not need to add constraints ensuring the consistency of the routes because we solve our subproblems by pricing. It will therefore be the responsibility of the pricing callback to create consistent routes. ##TODO link with additional constraints of Direct model + return model, x, y, z, cov +end -# We set the objective function: -@objective(decomp, Min, -sum(arc_costs[u, v] * z[u, v] for u in locations, v in locations) -+ -sum(facilities_fixed_costs[j] * y[j] for j in facilities)) # ## Pricing callback @@ -258,11 +274,13 @@ function z_contribution(enroute::EnumeratedRoute, z_red_costs) return z end + + # We are now able to write our pricing callback: function my_pricing_callback(cbdata) ## get the id of the facility - j = BlockDecomposition.indice(BlockDecomposition.callback_spid(cbdata, decomp)) + j = BlockDecomposition.indice(BlockDecomposition.callback_spid(cbdata, model)) ## retrieve variables reduced costs z_red_costs = Dict( @@ -297,23 +315,19 @@ function my_pricing_callback(cbdata) sol_cost = min_reduced_cost ## Submit the solution of the subproblem to Coluna - MOI.submit(decomp, BlockDecomposition.PricingSolution(cbdata), sol_cost, sol_vars, sol_vals) + MOI.submit(model, BlockDecomposition.PricingSolution(cbdata), sol_cost, sol_vars, sol_vals) ## Submit the dual bound to the solution of the subproblem ## This bound is used to compute the contribution of the subproblem to the lagrangian ## bound in column generation. - MOI.submit(decomp, BlockDecomposition.PricingDualBound(cbdata), sol_cost) # optimal solution + MOI.submit(model, BlockDecomposition.PricingDualBound(cbdata), sol_cost) # optimal solution end - -# Set decomposition -#@dantzig_wolfe_decomposition(decomp, dec, Base_axis) -#subproblems = BlockDecomposition.getsubproblems(dec) -#specify!.(subproblems, lower_multiplicity=0, upper_multiplicity=nb_routes_per_locations, solver = my_pricing_callback) #set solver = nothing to test without my_pricing_callback -#subproblemrepresentative.(z, Ref(subproblems)) -##Solve the model -#JuMP.optimize!(decomp) +# Create the model: +(model, x, y, z, _) = create_model() +# Solve: +JuMP.optimize!(model) # ## TODO: display "raw" decomp model output and comment, transition to next section @@ -353,20 +367,14 @@ function valid_inequalities_callback(cbdata) ## Add the valid inequalities to the model. for ineq in inequalities constr = JuMP.@build_constraint(x[ineq.customer_id, ineq.facility_id] <= y[ineq.facility_id]) - MOI.submit(decomp, MOI.UserCut(cbdata), constr) + MOI.submit(model, MOI.UserCut(cbdata), constr) end end -MOI.set(decomp, MOI.UserCutCallback(), valid_inequalities_callback); - - - -#MOI.set(decomp, MOI.UserCutCallback(), valid_inequalities_callback); -#@dantzig_wolfe_decomposition(decomp, dec, Base_axis) -#subproblems = BlockDecomposition.getsubproblems(dec) -#specify!.(subproblems, lower_multiplicity=0, upper_multiplicity=nb_routes_per_locations, solver = my_pricing_callback) #set solver = nothing to test without my_pricing_callback -#subproblemrepresentative.(z, Ref(subproblems)) -#JuMP.optimize!(decomp) +# We re-declare the model and solve it with the inequalities_callback: +(model, x, y, z, _) = create_model() +MOI.set(model, MOI.UserCutCallback(), valid_inequalities_callback); +JuMP.optimize!(model) # ## TODO: comment on the improvement of the dual bound, fix display of both raw decomp and valid ineq decomp model @@ -383,8 +391,11 @@ MOI.set(decomp, MOI.UserCutCallback(), valid_inequalities_callback); struct CoverConstrData <: BlockDecomposition.AbstractCustomData customer::Int end + +(model, x, y, z, cov) = create_model() + # We declare our custom data to Coluna -BlockDecomposition.customconstrs!(decomp, CoverConstrData); +BlockDecomposition.customconstrs!(model, CoverConstrData); # And we attach one custom data to each cover constraint for i in customers customdata!(cov[i], CoverConstrData(i)) @@ -411,8 +422,8 @@ struct R1cCutData <: BlockDecomposition.AbstractCustomData end # We declare our custom data to Coluna: -BlockDecomposition.customvars!(decomp, R1cVarData) -BlockDecomposition.customconstrs!(decomp, [CoverConstrData, R1cCutData]); +BlockDecomposition.customvars!(model, R1cVarData) +BlockDecomposition.customconstrs!(model, [CoverConstrData, R1cCutData]); # This method is called by Coluna to compute the coefficients of the `λ_k` in the cuts: function Coluna.MathProg.computecoeff( @@ -463,7 +474,7 @@ function r1c_callback(cbdata) end if violation > 1 ## create the constraint and add it to the model, use custom data to keep information about the cut (= the subset of considered cover constraints) - MOI.submit(decomp, + MOI.submit(model, MOI.UserCut(cbdata), JuMP.ScalarConstraint(JuMP.AffExpr(0.0), MOI.LessThan(1.0)), R1cCutData(cov_constr_subset) @@ -489,7 +500,7 @@ end # We re-write our pricing callback, with the additional contribution that corresponds to R1Cs cost: function my_pricing_callback(cbdata) - j = BlockDecomposition.indice(BlockDecomposition.callback_spid(cbdata, decomp)) + 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 ) @@ -528,20 +539,16 @@ function my_pricing_callback(cbdata) z_vars = [z[u, v] for (u,v) in best_route_arcs] x_vars = [x[i, j] for i in best_route_customers] sol_vars = vcat(z_vars, x_vars) - sol_vals = vcat([1.0 for _ in z_vars], [1.0 for _ in x_vars]) + sol_vals = ones(Float64, length(z_vars) + length(x_vars)) sol_cost = min_reduced_cost ## Submit the solution of the subproblem to Coluna ## TODO: comment on custom data here - MOI.submit(decomp, BlockDecomposition.PricingSolution(cbdata), sol_cost, sol_vars, sol_vals, R1cVarData(best_route.path)) - MOI.submit(decomp, BlockDecomposition.PricingDualBound(cbdata), sol_cost) + MOI.submit(model, BlockDecomposition.PricingSolution(cbdata), sol_cost, sol_vars, sol_vals, R1cVarData(best_route.path)) + MOI.submit(model, BlockDecomposition.PricingDualBound(cbdata), sol_cost) end -MOI.set(decomp, MOI.UserCutCallback(), r1c_callback); -@dantzig_wolfe_decomposition(decomp, dec, Base_axis) -subproblems = BlockDecomposition.getsubproblems(dec) -specify!.(subproblems, lower_multiplicity=0, upper_multiplicity=nb_routes_per_locations, solver = my_pricing_callback) -subproblemrepresentative.(z, Ref(subproblems)) -JuMP.optimize!(decomp) \ No newline at end of file +MOI.set(model, MOI.UserCutCallback(), r1c_callback); +JuMP.optimize!(model) \ No newline at end of file From 77ae7f2135463d9a94039e092cb7f9dfc8afeddb Mon Sep 17 00:00:00 2001 From: Natacha Javerzat Date: Wed, 10 May 2023 15:45:33 +0200 Subject: [PATCH 3/6] mineure --- docs/src/start/advanced-demo.jl | 36 ++++++++++++++++----------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/docs/src/start/advanced-demo.jl b/docs/src/start/advanced-demo.jl index 5a0548e27..412fa4403 100644 --- a/docs/src/start/advanced-demo.jl +++ b/docs/src/start/advanced-demo.jl @@ -37,7 +37,7 @@ arc_costs = locations = vcat(facilities, customers) nb_customers = length(customers) nb_facilities = length(facilities) -positions = 1:length(nb_positions) +positions = 1:nb_positions # Let's define the model that will solve our problem; first we need to load some packages and dependencies: @@ -67,28 +67,28 @@ model = JuMP.direct_model(HiGHS.Optimizer()) @variable(model, z[u in locations, v in locations], Bin) @variable(model, x[i in customers, j in facilities, k in routes_per_locations, p in positions], Bin) -# variables `y_j` indicate the opening status of a given facility ; if `y_j = 1` then facility `j` is open, otherwise `y_j = 0` -# variables `z_(u,v)` indicate which arcs are used ; if `z_(u,v) = 1` then there is an arc between location `u` and location `v`, otherwise `z_(u,v) = 0` -# variables `x_(i, j, k, p)` are used to express cover constraints and to ensure the consistency of the routes ; `x_(i, j, k, p) = 1` if customer `i` is delivered from facility `j` at the position `p` of route `k`, otherwise `x_(i, j, k, p) = 0`. +# - variables `y_j` indicate the opening status of a given facility ; if `y_j = 1` then facility `j` is open, otherwise `y_j = 0` +# - variables `z_(u,v)` indicate which arcs are used ; if `z_(u,v) = 1` then there is an arc between location `u` and location `v`, otherwise `z_(u,v) = 0` +# - variables `x_(i, j, k, p)` are used to express cover constraints and to ensure the consistency of the routes ; `x_(i, j, k, p) = 1` if customer `i` is delivered from facility `j` at the position `p` of route `k`, otherwise `x_(i, j, k, p) = 0`. # Now, we add constraints to our model: -# "each customer is visited by at least one route" +# - "each customer is visited by at least one route" @constraint(model, cov[i in customers], sum(x[i, j, k, p] for j in facilities, k in routes_per_locations, p in positions) >= 1) -# "for any route from any facility, its length does not exceed the fixed maximum length `nb_positions`" +# - "for any route from any facility, its length does not exceed the fixed maximum length `nb_positions`" @constraint(model, cardinality[j in facilities, k in routes_per_locations], sum(x[i, j, k, p] for i in customers, p in positions) <= nb_positions * y[j]) -# "only one customer can be delivered by a given route at a given position" +# - "only one customer can be delivered by a given route at a given position" @constraint(model, assign_setup[p in positions, j in facilities, k in routes_per_locations], sum(x[i, j, k, p] for i in customers) <= y[j]) -# "a customer can only be delivered at position `p > 1` of a given route if there is a customer delivered at position `p-1` of the same route" +# - "a customer can only be delivered at position `p > 1` of a given route if there is a customer delivered at position `p-1` of the same route" @constraint(model, open_route[j in facilities, k in routes_per_locations, p in positions; p > 1], sum(x[i, j, k, p] for i in customers) <= sum(x[i, j, k, p-1] for i in customers)) -# "there is an arc between two customers whose demand is satisfied by the same route at consecutive positions" +# - "there is an arc between two customers whose demand is satisfied by the same route at consecutive positions" @constraint(model, route_arc[i in customers, l in customers, indi in 1:nb_customers, indl in 1:nb_customers, j in facilities, k in routes_per_locations, p in positions; p > 1 && i != l && indi != indl], z[customers[indi], customers[indl]] >= x[l, j, k, p] + x[i, j, k, p-1] - 1) -# "there is an arc between the facility `j` and the first customer visited by the route `k` from facility `j`" +# - "there is an arc between the facility `j` and the first customer visited by the route `k` from facility `j`" @constraint(model, start_arc[i in customers, indi in 1:nb_customers, j in facilities, k in routes_per_locations], z[facilities[j], customers[indi]] >= x[i, j, k, 1]) @@ -123,13 +123,13 @@ coluna = optimizer_with_attributes( conqueralg = Coluna.ColCutGenConquer( primal_heuristics = [ #Coluna.ParameterizedHeuristic( - # Diva.Diving(), - # 1.0, - # 1.0, - # 1, - # 1, - # "Diving" - #) + ## Diva.Diving(), + ## 1.0, + ## 1.0, + ## 1, + ## 1, + ## "Diving" + ##) ] ) ) # default branch-cut-and-price @@ -409,7 +409,7 @@ end # The rank-one cuts we are going to add are of the form: # `sum(c_k λ_k) <= 1.0` -# for a fixed subset `r1c_cov_constrs` of cover constraints of size 3, with `λ_k` the master columns variables and c_k` s.t. +# for a fixed subset `r1c_cov_constrs` of cover constraints of size 3, with `λ_k` the master columns variables and `c_k` s.t. # `c_k = ⌊ 1/2 x |r1c_locations ∩ r1c_cov_constrs| ⌋` # with `r1c_locations` the current solution (route) that corresponds to `λ_k`. # e.g. if we consider cover constraints cov[3], cov[6] and cov[8] in our cut, then the route 1-4-6-7 gives a zero coefficient while the route 1-4-6-3 gives a coefficient equal to one. From 25653f7835d5a0df93f304db1c11078bfc7c2cca Mon Sep 17 00:00:00 2001 From: Natacha Javerzat Date: Wed, 10 May 2023 15:58:25 +0200 Subject: [PATCH 4/6] maj TODOs --- docs/src/start/advanced-demo.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/start/advanced-demo.jl b/docs/src/start/advanced-demo.jl index 412fa4403..d492ec2fe 100644 --- a/docs/src/start/advanced-demo.jl +++ b/docs/src/start/advanced-demo.jl @@ -169,10 +169,10 @@ function create_model() @constraint(model, open_facility[j in facilities], sum(z[j, i] for i in customers) <= y[j] * nb_routes_per_locations) - ## We do not need to add constraints ensuring the consistency of the routes because we solve our subproblems by pricing. It will therefore be the responsibility of the pricing callback to create consistent routes. ##TODO link with additional constraints of Direct model - + + ## Contrary to the direct model, we are not obliged here to add constraints to ensure the consistency of the routes because we solve our subproblems by pricing. It will therefore be the responsibility of the pricing callback to create consistent routes. + ## We set the objective function: - @objective(model, Min, sum(arc_costs[u, v] * z[u, v] for u in locations, v in locations) + @@ -381,13 +381,13 @@ end MOI.set(model, MOI.UserCutCallback(), valid_inequalities_callback); JuMP.optimize!(model) -# ## TODO: comment on the improvement of the dual bound, fix display of both raw decomp and valid ineq decomp model +# ## TODO: comment on the improvement of the dual bound # ## Strengthen with non-robust cuts (rank-one cuts) -# ##TODO: describe the goal by looking at lambdas values -> highlight that the aim is to drive them towards integrality +# ##TODO: describe the goal by looking at lambdas values (if possible?) -> highlight that the aim is to drive them towards integrality # Here, we implement special types of cuts called "rank-one cuts" (R1C). These cuts are non-robust in the sense that they can not be expressed only with the original variables of the model. In particular, they have to be expressed with the master columns variables `λ_k`. # R1Cs are obtained by applying the Chvátal-Gomory procedure once, hence their name, on cover constraints. We must therefore be able to differentiate the cover constraints from the other constraints of the model. To do this, we exploit an advantage of Coluna that allows us to attach custom data to the constraints and variables of our model: From 048e99e07b9ababbda3135bce364987d278d977b Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Wed, 10 May 2023 16:00:59 +0200 Subject: [PATCH 5/6] Update docs/src/start/advanced-demo.jl --- docs/src/start/advanced-demo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/start/advanced-demo.jl b/docs/src/start/advanced-demo.jl index d492ec2fe..231661b98 100644 --- a/docs/src/start/advanced-demo.jl +++ b/docs/src/start/advanced-demo.jl @@ -132,7 +132,7 @@ coluna = optimizer_with_attributes( ##) ] ) - ) # default branch-cut-and-price + ) ## default branch-cut-and-price ), "default_optimizer" => GLPK.Optimizer # GLPK for the master & the subproblems ) From 99bc70e7bdfd1a5fd8cc10ea2a7f52d46c7626ba Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Wed, 10 May 2023 16:01:05 +0200 Subject: [PATCH 6/6] Update docs/src/start/advanced-demo.jl --- docs/src/start/advanced-demo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/start/advanced-demo.jl b/docs/src/start/advanced-demo.jl index 231661b98..56b123698 100644 --- a/docs/src/start/advanced-demo.jl +++ b/docs/src/start/advanced-demo.jl @@ -122,7 +122,7 @@ coluna = optimizer_with_attributes( maxnumnodes = 0, conqueralg = Coluna.ColCutGenConquer( primal_heuristics = [ - #Coluna.ParameterizedHeuristic( + ## Coluna.ParameterizedHeuristic( ## Diva.Diving(), ## 1.0, ## 1.0,