Duality handlers

This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.

The purpose of this tutorial is to demonstrate trivial examples that expose the strengths and weaknesses of each duality handler.

For more information on SDDP.jl's duality handlers, see Integrality.

This tutorial uses the following packages:

using SDDP
import HiGHS

Note that these trivial examples exposed a bug in HiGHS, which we work-around by turning off presolve. In real examples, you should probably leave the default options.

Optimizer = optimizer_with_attributes(HiGHS.Optimizer, "presolve" => "off")

MathOptInterface.OptimizerWithAttributes(HiGHS.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("presolve") => "off"])

Training function

First, we need a function to simplify our testing:

function train_and_evaluate_bounds(
    model_fn::Function,
    duality_handler::SDDP.AbstractDualityHandler;
    print_level::Int = 0,
    kwargs...,
)
    model = model_fn()
    SDDP.train(model; print_level, duality_handler, kwargs...)
    simulations = SDDP.simulate(model, 1)
    lower_bound = SDDP.calculate_bound(model)
    println("lower_bound: $(lower_bound)")
    upper_bound = sum(data[:stage_objective] for data in only(simulations))
    println("upper_bound: $(upper_bound)")
    return
end

train_and_evaluate_bounds (generic function with 1 method)

This function builds a new model, trains it using the provided duality_handler, and then prints the lower and upper bounds. We'll assume that the models we pass in are deterministic so we need to conduct only a single simulation to evaluate the upper bound.

ContinuousConicDuality

The default duality handler in SDDP.jl is ContinuousConicDuality. To compute a cut, it solves the continuous relaxation of the MIP.

In the same way that solution to a relaxed linear program may be far from the optimal MIP solution, the biggest downside to ContinuousConicDuality is that many models have large gaps between the lower and upper bound:

function model_1()
    return SDDP.LinearPolicyGraph(;
        stages = 2,
        lower_bound = 0.0,
        optimizer = Optimizer,
    ) do sp, t
        @variable(sp, x, Bin, SDDP.State, initial_value = 1.0)
        @variable(sp, y, Bin)
        @constraint(sp, x.out == x.in)
        if t == 1
            @stageobjective(sp, x.out)
        else
            @stageobjective(sp, y)
            @constraint(sp, y >= x.in - 0.5)
        end
    end
end

train_and_evaluate_bounds(model_1, SDDP.ContinuousConicDuality())

lower_bound: 1.5
upper_bound: 2.0

StrengthenedConicDuality

One technique to improve upon ContinuousConicDuality is StrengthenedConicDuality. Without going into the technical details, both use the continuous relaxation to compute a valid subgradient for the cut. StrengthenedConicDuality then tries to improve the cut by solving an additional integer program.

In this example, StrengthenedConicDuality can improve upon the lower bound and prove that the upper bound of 2.0 is optimal:

train_and_evaluate_bounds(model_1, SDDP.StrengthenedConicDuality())

lower_bound: 2.0
upper_bound: 2.0

Sometimes, however, StrengthenedConicDuality cannot improve upon ContinuousConicDuality:

function model_2()
    return SDDP.LinearPolicyGraph(;
        stages = 2,
        lower_bound = 0.0,
        optimizer = Optimizer,
    ) do sp, t
        @variable(sp, x, SDDP.State, initial_value = 0.1)
        @variable(sp, y, Int)
        @constraint(sp, x.out == x.in)
        if t == 1
            @stageobjective(sp, x.out)
        else
            @stageobjective(sp, y)
            @constraint(sp, y >= x.in + 0.1)
            @constraint(sp, y >= -x.in + 0.1)
        end
    end
end

train_and_evaluate_bounds(model_2, SDDP.ContinuousConicDuality())

lower_bound: 0.30000000000000004
upper_bound: 1.1

train_and_evaluate_bounds(model_2, SDDP.StrengthenedConicDuality())

lower_bound: 0.3
upper_bound: 1.1

Even though it is sometimes tighter than ContinuousConicDuality and it can never be worse, StrengthenedConicDuality is not the default duality handler because it is more expensive to compute; it solves a mixed-integer program whereas ContinuousConicDuality solves a continuous relaxation.

LagrangianDuality

A technique to improve upon StrengthenedConicDuality is LagrangianDuality. Without going into the technical details, both use the continuous relaxation to compute a valid subgradient for the cut, but, where StrengthenedConicDuality tries to improve the cut by solving a single additional integer program, LagrangianDuality may solve many integer programs.

LagrangianDuality finds the optimal policy for model_1:

train_and_evaluate_bounds(model_1, SDDP.LagrangianDuality())

lower_bound: 2.0
upper_bound: 2.0

and also for model_2:

train_and_evaluate_bounds(model_2, SDDP.LagrangianDuality())

lower_bound: 1.0999999999999999
upper_bound: 1.1

Sometimes, however, LagrangianDuality does not close the gap:

function model_3()
    return SDDP.LinearPolicyGraph(;
        stages = 2,
        lower_bound = -1.0,
        optimizer = Optimizer,
    ) do sp, t
        @variable(sp, -1 <= x <= 0.5, SDDP.State, initial_value = 0.0)
        @variable(sp, y)
        @stageobjective(sp, y)
        if t == 1
            @constraint(sp, y >= x.out)
            @constraint(sp, y >= -x.out)
        else
            @variable(sp, z, Bin)
            @constraint(sp, y >= 1 - x.in - 3 * z)
            @constraint(sp, y >= 1 + x.in - 3 * (1 - z))
        end
    end
end

train_and_evaluate_bounds(model_3, SDDP.LagrangianDuality())

lower_bound: 0.3333333333333333
upper_bound: 1.0

but it may still be better than StrengthenedConicDuality:

train_and_evaluate_bounds(model_3, SDDP.StrengthenedConicDuality())

lower_bound: -0.5
upper_bound: 1.0

The algorithm behind LagrangianDuality is significantly more complicated than StrengthenedConicDuality. For some models it can be helpful, for others, the increased computational cost is not worth the improvement in the tightness of the value function.

Different policies

So far, the different duality handlers have led to different lower bounds, but identical upper bounds. This is an artifact of our trivial examples. Using a more sophisticated duality handler can improve the lower bound and lead to a cheaper policy (the upper bound). Here's an example:

function model_4()
    return SDDP.LinearPolicyGraph(;
        stages = 2,
        lower_bound = -1.0,
        optimizer = Optimizer,
    ) do sp, t
        @variable(sp, -1 <= x <= 0.5, SDDP.State, initial_value = 0.0)
        @variable(sp, y)
        @variable(sp, z, Bin)
        if t == 1
            @stageobjective(sp, -0.1 * x.out)
        else
            @stageobjective(sp, y)
            @constraint(sp, y >= 1 - x.in - 3 * z)
            @constraint(sp, y >= 1 + x.in - 3 * (1 - z))
        end
    end
end

model_4 (generic function with 1 method)

ContinuousConicDuality finds in a policy that costs 0.45:

train_and_evaluate_bounds(model_4, SDDP.ContinuousConicDuality())

lower_bound: -0.55
upper_bound: 0.45

whereas LagrangianDuality finds a policy that costs 0.1:

train_and_evaluate_bounds(model_4, SDDP.LagrangianDuality())

lower_bound: 0.1
upper_bound: 0.1

This relationship is not guaranteed to hold. In some models ContinuousConicDuality may find a cheaper policy than LagrangianDuality, even though the latter finds a tighter lower bound. In general, you should experiment with different duality handlers to see what works best for your problem.

BanditDuality

The trade-off between the computational cost and the tightness of a formulation can be tricky to manage. SDDP.jl includes BanditDuality, which is an algorithm that does not appear in the published academic literature. The BanditDuality duality handler treats the problem of choosing a duality handler for each iteration of the SDDP algorithm as a multi-armed bandit problem, where the reward is the change in the lower bound after each iteration per second of computation time. The multi-armed bandit problem allows us to trade off many fast but weak iterations of ContinuousConicDuality against a small number of relatively strong iterations of LagrangianDuality.

duality_handler = SDDP.BanditDuality(
    SDDP.ContinuousConicDuality(),
    SDDP.StrengthenedConicDuality(),
    SDDP.LagrangianDuality(),
)
train_and_evaluate_bounds(model_1, duality_handler)

lower_bound: 2.0
upper_bound: 2.0

train_and_evaluate_bounds(model_2, duality_handler)

lower_bound: 1.0999999999999999
upper_bound: 1.1

train_and_evaluate_bounds(model_3, duality_handler)

lower_bound: 0.3333333333333333
upper_bound: 1.0

The BanditDuality is often a very good choice to use in practice. It is not the default because a tighter lower bound does not always lead to a better policy, so we opt for the simplest and fastest duality handler as the default.

FixedDiscreteDuality

ContinuousConicDuality, StrengthenedConicDuality, and LagrangianDuality) all have the useful feature that their cuts produce valid lower bounds on the optimal policy cost.

FixedDiscreteDuality is a duality handler that does NOT guarantee a valid lower bound.

It works by first solving the mixed-integer problem, fixing the discrete variables to their optimal value, and then solving the continuous relaxation. This has the same computational cost as StrengthenedConicDuality.

The resulting cuts are not guaranteed to be globally valid, and they may cut off part of the true value function. This means that the "lower" bound (computed using calculate_bound) may be greater than the optimal policy cost.

Sometimes FixedDiscreteDuality works and finds the same solutions as StrengthenedConicDuality:

train_and_evaluate_bounds(model_1, SDDP.FixedDiscreteDuality())

lower_bound: 2.0
upper_bound: 2.0

train_and_evaluate_bounds(model_1, SDDP.StrengthenedConicDuality())

lower_bound: 2.0
upper_bound: 2.0

For other models, FixedDiscreteDuality can find solutions that are tighter than StrengthenedConicDuality:

train_and_evaluate_bounds(model_2, SDDP.FixedDiscreteDuality())

lower_bound: 1.1
upper_bound: 1.1

train_and_evaluate_bounds(model_2, SDDP.StrengthenedConicDuality())

lower_bound: 0.3
upper_bound: 1.1

In other cases, FixedDiscreteDuality finds a solution that looks optimal, but because the lower bound lies, we have no way of formally verifying if the upper bound is optimal:

train_and_evaluate_bounds(model_3, SDDP.FixedDiscreteDuality())

lower_bound: 1.0
upper_bound: 1.0

As another example, here's a model that also appears to be optimal:

train_and_evaluate_bounds(model_4, SDDP.FixedDiscreteDuality())

lower_bound: 0.45
upper_bound: 0.45

But LagrangianDuality can find (and prove) that there is a better solution:

train_and_evaluate_bounds(model_4, SDDP.LagrangianDuality())

lower_bound: 0.1
upper_bound: 0.1

We implemented FixedDiscreteDuality because it may be useful. In some cases, it can find solutions that are tighter than StrengthenedConicDuality for the same computational cost. However, because the cuts are not guaranteed to be globally valid, you should exercise great caution when interpreting the solution. It may be a useful heuristic, or it may not. Use FixedDiscreteDuality at your own risk.