Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

re-organize demo, fix typos #900

Merged
merged 1 commit into from
May 12, 2023
Merged

re-organize demo, fix typos #900

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
218 changes: 112 additions & 106 deletions docs/src/start/advanced-demo.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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,
Expand All @@ -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:
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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:

Expand Down Expand Up @@ -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);
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -555,19 +553,19 @@ 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.

# 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".

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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]))
Expand All @@ -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
Expand Down Expand Up @@ -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.

Expand Down Expand Up @@ -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)