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.

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

Benders documentation #975

Merged
merged 4 commits into from
Jul 4, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
71 changes: 52 additions & 19 deletions docs/src/api/benders.md
Original file line number Diff line number Diff line change
Expand Up @@ -74,20 +74,19 @@ The subproblems have the following form:

```math
\begin{aligned}
\min \quad& fy + \mathbf{1}z' + \mathbf{1}z'' &&& \\
\text{s.t.} \quad& Dy + z' \geq d - B\bar{x} && (5) \quad& {\color{blue}(\pi)} \\
& Ey + z'' \geq e && (6) \quad& {\color{blue}(\rho)} \\
\min \quad& fy + {\color{gray} \mathbf{1}z' + \mathbf{1}z''} &&& \\
\text{s.t.} \quad& Dy {\color{gray} + z'} \geq d - B\bar{x} && (5) \quad& {\color{blue}(\pi)} \\
& Ey {\color{gray} + z''} \geq e && (6) \quad& {\color{blue}(\rho)} \\
& l_2 \leq y \leq u_2 && (7) \quad& {\color{blue}(\sigma)}
\end{aligned}
```

where $y$ are second-stage variables, $z'$ and $z''$ are artificial variables,
where $y$ are second-stage variables, $z'$ and $z''$ are artificial variables (in grey because they are deactivated by default),
constraints (5) are the reformulation of linking constraints using the first-stage solution $\bar{x}$,
constraints (6) are the second-stage constraints,
and constraints (7) are the bounds on the second-stage variables.
In blue, we define the dual variables associated to these constraints.


**References**:

```@docs
Expand Down Expand Up @@ -127,6 +126,11 @@ Coluna.Benders.new_output

## Benders cut generation iteration

This is a description of how the `Coluna.Benders.run_benders_iteration!` generic function behaves with the default implementation.

These are the main steps of a Benders cut generation iteration without stabilization.
Click on the step to go to the corresponding section.

```mermaid
flowchart TB;
id1(Optimize master)
Expand Down Expand Up @@ -158,12 +162,15 @@ flowchart TB;
click id4 href "#Subproblem-iterator" "Link to doc"
click id5 href "#Separation-subproblem-optimization" "Link to doc"
click id6 href "#Set-of-generated-cuts" "Link to doc"
click id9 href "#Unboundeness-check" "Link to doc"
click id9 href "#Unboundedness-check" "Link to doc"
click id11 href "#Current-primal-solution" "Link to doc"
click id7 href "#Cuts-insertion" "Link to doc"
click id8 href "#Iteration-output" "Link to doc"
```

In the default implementation, some sections may have different behaviors depending on the
result of previous steps.

### Master optimization

This operation consists in optimize the master problem in order to find a first-level
Expand Down Expand Up @@ -204,7 +211,6 @@ If the solver does not provide a dual infeasibility certificate, the implementat
has an "emergency" routine to provide a first-stage feasible solution by solving the master LP with cost of second stage variables set to zero.
We recommend using a solver that provides a dual infeasibility certificate and avoiding the "emergency" routine.


**References**:

```@docs
Expand All @@ -215,19 +221,34 @@ Go back to the [cut generation iteration diagram](#Benders-cut-generation-iterat

### Setup separation subproblems

The separation subproblems differs depending on wether the master is unbounded or not.
!!! info
The separation subproblems differs depending on whether the restricted master is unbounded or not:
- if the restricted master is optimal, the generic function calls `Coluna.Benders.update_sp_rhs!`
- if the restricted master is unbounded, the generic function calls `Coluna.Benders.setup_separation_for_unbounded_master_case!`

If the master is unbounded, the generic function calls `Coluna.Benders.setup_separation_for_unbounded_master_case!`; otherwise, it calls `Coluna.Benders.update_sp_rhs!`.
Default implementation of `Coluna.Benders.update_sp_rhs!` updates the right-hand side of the linking constraints (5).

**Reference**:
```@docs
Coluna.Benders.update_sp_rhs!
```

Default implementation of `Coluna.Benders.setup_separation_for_unbounded_master_case!`
gives raise to the formulation proposed in Lemma 2 of Bonami et al.
gives raise to the formulation proposed in Lemma 2 of Bonami et al:

Default implementation of `Coluna.Benders.update_sp_rhs!` updates the right-hand side of the linking constraints (5).
```math
\begin{aligned}
(SepB) \equiv \min \quad& fy + {\color{gray} \mathbf{1}z' + \mathbf{1}z''} &&& \\
\text{s.t.} \quad& Dy {\color{gray} + z'} \geq -B\bar{x} && (5a) \quad& {\color{blue}(\pi)} \\
& Ey {\color{gray} + z''} \geq 0 && (6a) \quad& {\color{blue}(\rho)} \\
& y \geq 0 && (7a) \quad& {\color{blue}(\sigma)}
\end{aligned}
```
where $y$ are second-stage variables, $z'$ and $z''$ are artificial variables (in grey because they are deactivated by default), and $\bar{x}$ is the an unbounded ray of the restricted master.

**References**:
**Reference**:

```@docs
Coluna.Benders.update_sp_rhs!
Coluna.Benders.setup_separation_for_unbounded_master_case!
```

Expand Down Expand Up @@ -288,6 +309,11 @@ Coluna.Algorithm.CutsSet
Coluna.Algorithm.SepSolSet
```

The default implementation of `push_in_set!` has the responsibility to check if the cut is
violated. Given $\bar{eta}_k$ solution to the restricted master and $\bar{y}$ solution to the separation problem, the cut is considered as violated when:
- the separation subproblem was infeasible
- or $\bar{\eta}_k \geq f\bar{y}$

**References**:

```@docs
Expand All @@ -298,9 +324,18 @@ Coluna.Benders.push_in_set!

Go back to the [cut generation iteration diagram](#Benders-cut-generation-iteration).

### Unboundness check
### Unboundedness check

Lorem ipsum.
!!! info
This check is performed only when the restricted master is unbounded.

To perform this check, we need a solution to each separation problem.

Let $(\bar{\eta}_k)_{k \in K}$ be the value of second stage variables in the dual infeasibility certificate of the restricted master.
Let $\bar{y}$ be an optimal solution to the separation problem **(SepB)**.

As indicated by Bonami et al., if $f\bar{y} \leq \sum_{k \in K} \bar{\eta}_k$, then the
original problem is unbounded (by definition of an unbounded ray of the original problem).

**References**:

Expand All @@ -310,11 +345,9 @@ Coluna.Benders.master_is_unbounded

### Cuts insertion

The default implementation inserts into the master all the cuts stored in the `` object.
The default implementation inserts into the master all the cuts stored in the `CutsSet` object.

**References**:

**References**:
**Reference**:

```@docs
Coluna.Benders.insert_cuts!
Expand Down
105 changes: 84 additions & 21 deletions src/Algorithm/benders/default.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,8 @@
"""
BendersContext(reformulation, algo_params) -> BendersContext

Default implementation of the Benders algorithm.
"""
struct BendersContext <: Benders.AbstractBendersContext
reform::Reformulation
optim_sense
Expand Down Expand Up @@ -47,6 +52,17 @@ function Benders.setup_reformulation!(reform::Reformulation, env)
return
end


"""
Output of the default implementation of the `Benders.optimize_master_problem!` method.

It contains:
- `ip_solver`: `true` if the master problem is solved with a MIP solver and involves integral variables, `false` otherwise.
- `result`: the result of the master problem optimization stored in an `OptimizationState` object.
- `infeasible`: `true` if the master at the current iteration is infeasible; `false` otherwise.
- `unbounded`: `true` if the master at the current iteration is unbounded; `false` otherwise.
- `certificate`: `true` if the master at the current iteration is unbounded and if the current result is a dual infeasibility certificate, `false` otherwise.
"""
struct BendersMasterResult{F}
ip_solver::Bool
result::OptimizationState{F}
Expand Down Expand Up @@ -202,13 +218,27 @@ function Benders.setup_separation_for_unbounded_master_case!(ctx::BendersContext
return
end

"""
Solution to the separation problem together with its corresponding benders cut.

It contains:
- `min_sense`: `true` if it's a minimization problem; `false` otherwise.
- `lhs`: the left-hand side of the cut.
- `rhs`: the right-hand side of the cut.
- `dual_sol`: an optimal dual solution to the separation problem.
"""
struct GeneratedCut{F}
min_sense::Bool
lhs::Dict{VarId, Float64}
rhs::Float64
dual_sol::DualSolution{F}
end

"""
Stores a collection of cuts.

It contains `cuts` a vector of `GeneratedCut` objects.
"""
struct CutsSet
cuts::Vector{GeneratedCut}
CutsSet() = new(GeneratedCut[])
Expand All @@ -218,28 +248,48 @@ Base.iterate(set::CutsSet, state) = iterate(set.cuts, state)

Benders.set_of_cuts(::BendersContext) = CutsSet()

"""
Primal solutions to the separation problems optimized at the current iteration.
This is used to build a primal solution.

It contains `sols` a vector of primal solutions.
"""
struct SepSolSet{F}
sols::Vector{MathProg.PrimalSolution{F}}
end
SepSolSet{F}() where {F} = SepSolSet{F}(MathProg.PrimalSolution{F}[])
Benders.set_of_sep_sols(::BendersContext) = SepSolSet{MathProg.Formulation{MathProg.BendersSp}}()

"""
Output of the default implementation of the `Benders.optimize_separation_problem!` and
`Benders.treat_infeasible_separation_problem_case!` methods.

It contains:
- `second_stage_estimation_in_master`: the value of the second stage cost variable in the solution to the master problem.
- `second_stage_cost`: the value of the second stage cost variable in the solution to the separation problem.
- `lp_primal_sol`: the primal solution to the separation problem.
- `infeasible`: `true` if the current separation problem is infeasible; `false` otherwise.
- `unbounded`: `true` if the current separation problem is unbounded; `false` otherwise.
- `cut`: the cut generated by the separation problem.
- `infeasible_treatment`: `true` if this objects is an output of the `Benders.treat_infeasible_separation_problem_case!` method; `false` otherwise.
- `unbounded_master`: `true` if the separation subproblem has the form of Lemma 2 to separate a cut to truncate an unbounded ray of the restricted master problem; `false` otherwise.
"""
struct BendersSeparationResult{F}
second_stage_estimation_in_master::Float64
second_stage_cost::Union{Nothing,Float64}
lp_primal_sol::Union{Nothing,MathProg.PrimalSolution{F}}
infeasible::Bool
unbounded::Bool
certificate::Union{Nothing,MathProg.DualSolution{F}}
cut::Union{Nothing,GeneratedCut{F}}
infeasible_separation::Bool
infeasible_treatment::Bool
unbounded_master::Bool
end

Benders.get_obj_val(res::BendersSeparationResult) = res.second_stage_cost
Benders.get_primal_sol(res::BendersSeparationResult) = res.lp_primal_sol
Benders.is_infeasible(res::BendersSeparationResult) = res.infeasible
Benders.is_unbounded(res::BendersSeparationResult) = res.unbounded
Benders.get_dual_sol(res::BendersSeparationResult) = res.cut.dual_sol

## original MIP:
## min cx + dy s.t.
Expand Down Expand Up @@ -310,11 +360,11 @@ function Benders.optimize_separation_problem!(ctx::BendersContext, sp::Formulati
opt_state = run!(ctx.separation_solve_alg, env, sp, input)

if getterminationstatus(opt_state) == UNBOUNDED
return BendersSeparationResult{Formulation{BendersSp}}(estimated_cost, nothing, get_best_lp_primal_sol(opt_state), false, true, nothing, nothing, false, unbounded_master)
return BendersSeparationResult{Formulation{BendersSp}}(estimated_cost, nothing, get_best_lp_primal_sol(opt_state), false, true, nothing, false, unbounded_master)
end

if getterminationstatus(opt_state) == INFEASIBLE ## we then enter treat_infeasible_separation_problem_case! (phase 1)
return BendersSeparationResult{Formulation{BendersSp}}(estimated_cost, nothing, get_best_lp_primal_sol(opt_state), true, false, nothing, nothing, false, unbounded_master)
return BendersSeparationResult{Formulation{BendersSp}}(estimated_cost, nothing, get_best_lp_primal_sol(opt_state), true, false, nothing, false, unbounded_master)
end
## create and add cuts to the result
dual_sol = get_best_lp_dual_sol(opt_state)
Expand All @@ -324,10 +374,14 @@ function Benders.optimize_separation_problem!(ctx::BendersContext, sp::Formulati
cut_rhs = _compute_cut_rhs_contrib(ctx, sp, dual_sol)

cut = GeneratedCut(min_sense, cut_lhs, cut_rhs, dual_sol)
return BendersSeparationResult(estimated_cost, cost, get_best_lp_primal_sol(opt_state), false, false, dual_sol, cut, false, unbounded_master)
return BendersSeparationResult(estimated_cost, cost, get_best_lp_primal_sol(opt_state), false, false, cut, false, unbounded_master)
end

function Benders.master_is_unbounded(ctx::BendersContext, second_stage_cost, unbounded_master_case)
if !unbounded_master_case
return false
end

estimated_cost = 0
for (spid, sp) in Benders.get_benders_subprobs(ctx)
second_stage_cost_var = sp.duty_data.second_stage_cost_var
Expand All @@ -336,7 +390,7 @@ function Benders.master_is_unbounded(ctx::BendersContext, second_stage_cost, unb

min_sense = Benders.is_minimization(ctx)
sc = min_sense ? 1.0 : - 1.0
return sc * second_stage_cost < sc * estimated_cost + 1e-5 && unbounded_master_case
return sc * second_stage_cost < sc * estimated_cost + 1e-5
end

## it is a phase 1: add artificial variables in order to find a feasible solution
Expand All @@ -362,9 +416,8 @@ function Benders.treat_infeasible_separation_problem_case!(ctx::BendersContext,
end
_deactivate_art_vars(sp)

if getterminationstatus(opt_state) == INFEASIBLE # should not happen
error("A")
return BendersSeparationResult(estimated_cost, nothing, get_best_lp_primal_sol(opt_state), false, true, nothing, nothing, true, unbounded_master_case)
if getterminationstatus(opt_state) == INFEASIBLE
error("A") # should not happen
end

dual_sol = get_best_lp_dual_sol(opt_state)
Expand All @@ -379,14 +432,7 @@ function Benders.treat_infeasible_separation_problem_case!(ctx::BendersContext,
cut_lhs = _compute_cut_lhs(ctx, sp, dual_sol, true)
cut_rhs = _compute_cut_rhs_contrib(ctx, sp, dual_sol)
cut = GeneratedCut(min_sense, cut_lhs, cut_rhs, dual_sol)
return BendersSeparationResult(estimated_cost, cost, get_best_lp_primal_sol(opt_state), false, false, dual_sol, cut, true, unbounded_master_case)
end

function Benders.get_dual_sol(res::BendersSeparationResult)
if res.infeasible
return res.certificate
end
return res.cut.dual_sol
return BendersSeparationResult(estimated_cost, cost, get_best_lp_primal_sol(opt_state), false, false, cut, true, unbounded_master_case)
end

function Benders.push_in_set!(ctx::BendersContext, set::CutsSet, sep_result::BendersSeparationResult)
Expand All @@ -398,9 +444,8 @@ function Benders.push_in_set!(ctx::BendersContext, set::CutsSet, sep_result::Ben
eq = abs(sep_result.second_stage_cost - sep_result.second_stage_estimation_in_master) < 1e-5
gt = sc * sep_result.second_stage_cost + 1e-5 > sc * sep_result.second_stage_estimation_in_master


# if cost of separation result > second cost variable in master result
if !eq && gt || sep_result.infeasible_separation
if !eq && gt || sep_result.infeasible_treatment
push!(set.cuts, sep_result.cut)
return true
end
Expand Down Expand Up @@ -451,8 +496,6 @@ function Benders.insert_cuts!(reform, ctx::BendersContext, cuts)
warning = CutAlreadyInsertedBendersWarning(
in_master, is_active, cut_id, master, spform
)
# @warn warning
# push!(cuts_to_insert, cut)
throw(warning) # TODO: parameter
end
else
Expand Down Expand Up @@ -527,6 +570,17 @@ function Benders.build_primal_solution(context::BendersContext, mast_primal_sol,
)
end

"""
Output of the default implementation of an iteration of the Benders algorithm.

It contains:
- `min_sense`: the original problem is a minimization problem
- `nb_new_cuts`: the number of new cuts added to the master problem
- `ip_primal_sol`: the primal solution to the original problem found during this iteration
- `infeasible`: the original problem is infeasible
- `time_limit_reached`: the time limit was reached
- `master`: the solution value to the master problem
"""
struct BendersIterationOutput <: Benders.AbstractBendersIterationOutput
min_sense::Bool
nb_new_cuts::Int
Expand Down Expand Up @@ -560,6 +614,15 @@ function Benders.new_iteration_output(
)
end

"""
Output of the default implementation of the Benders algorithm.

It contains:
- `infeasible`: the original problem is infeasible
- `time_limit_reached`: the time limit was reached
- `mlp`: the final bound obtained with the Benders cut algorithm
- `ip_primal_sol`: the best primal solution to the original problem found by the Benders cut algorithm
"""
struct BendersOutput <: Benders.AbstractBendersOutput
infeasible::Bool
time_limit_reached::Bool
Expand Down
6 changes: 6 additions & 0 deletions src/Algorithm/benders/printer.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
"""
BendersPrinterContext(reformulation, algo_params) -> BendersPrinterContext

Creates a context to run the default implementation of the Benders algorithm
together with a printer that prints information about the algorithm execution.
"""
mutable struct BendersPrinterContext
inner::BendersContext
sp_elapsed_time::Float64
Expand Down
2 changes: 1 addition & 1 deletion src/Algorithm/colgen/default.jl
Original file line number Diff line number Diff line change
Expand Up @@ -569,7 +569,7 @@ struct GeneratedColumn
end

"""
Columns generated at the current iterations that forms the "current primal solution".
Columns generated at the current iteration that forms the "current primal solution".
This is used to compute the Lagragian dual bound.

It contains:
Expand Down
Loading