Computing the duals of a mixed-integer program

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

This tutorial explains how to compute dual variables (shadow prices) for a mixed-integer linear program by fixing the discrete variables to their optimal values and resolving as a linear program. Dual information is not available directly from MIP solvers.

Learning intentions:

• Understand why MIP solvers cannot return dual values, and why fixing integer variables to their optimal values and re-solving as an LP is a valid workaround
• Manually relax integer variables using fix and unset_binary, re-solve the resulting LP, and recover the dual values
• Use fix_discrete_variables to automate the relaxation, and call the returned undo function to restore the original MIP formulation

Required packages

This tutorial uses the following packages:

using JuMP
import HiGHS

The model

Our example model is the unit commitment example from Unit commitment. The details are unimportant, other than to note that there are two types of continuous variables, g and w, representing the quantity of generation from thermal and wind plants, and a discrete variable dispatch, which is 1 if plant i is operating, and 0 if not.

We are interested in the "dual" of the power_balance constraint, because it represents the marginal price of electricity that consumers should pay for their consumption.

generators = [
    (min = 0.0, max = 1000.0, fixed_cost = 1000.0, variable_cost = 50.0),
    (min = 300.0, max = 1000.0, fixed_cost = 0.0, variable_cost = 100.0),
]
N = length(generators)
model = Model(HiGHS.Optimizer)
set_silent(model)
@variables(model, begin
    generators[i].min <= g[i=1:N] <= generators[i].max
    0 <= w <= 200
    dispatch[i=1:N], Bin
end)
@constraints(model, begin
    power_balance, sum(g[i] for i in 1:N) + w == 1500
    [i = 1:N], g[i] <= generators[i].max * dispatch[i]
    [i = 1:N], g[i] >= generators[i].min * dispatch[i]
end)
@objective(
    model,
    Min,
    sum(
        generators[i].fixed_cost * dispatch[i] +
        generators[i].variable_cost * g[i] for i in 1:N
    )
)
print(model)

Min 1000 dispatch[1] + 50 g[1] + 100 g[2]
Subject to
 power_balance : g[1] + g[2] + w = 1500
 g[1] ≥ 0
 g[2] - 300 dispatch[2] ≥ 0
 g[1] - 1000 dispatch[1] ≤ 0
 g[2] - 1000 dispatch[2] ≤ 0
 g[1] ≥ 0
 g[2] ≥ 300
 w ≥ 0
 g[1] ≤ 1000
 g[2] ≤ 1000
 w ≤ 200
 dispatch[1] binary
 dispatch[2] binary

Manually fix the variables

If we optimize this model, we obtain a dual_status of NO_SOLUTION:

optimize!(model)
assert_is_solved_and_feasible(model)
dual_status(model)

NO_SOLUTION::ResultStatusCode = 0

This is because HiGHS cannot compute the duals of a mixed-integer program. We can work around this problem by fixing the integer variables to their optimal solution, relaxing integrality, and re-solving as a linear program.

discrete_values = round.(Int, value.(dispatch))
fix.(dispatch, discrete_values; force = true)
unset_binary.(dispatch)
print(model)

Min 1000 dispatch[1] + 50 g[1] + 100 g[2]
Subject to
 power_balance : g[1] + g[2] + w = 1500
 g[1] ≥ 0
 g[2] - 300 dispatch[2] ≥ 0
 g[1] - 1000 dispatch[1] ≤ 0
 g[2] - 1000 dispatch[2] ≤ 0
 dispatch[1] = 1
 dispatch[2] = 1
 g[1] ≥ 0
 g[2] ≥ 300
 w ≥ 0
 g[1] ≤ 1000
 g[2] ≤ 1000
 w ≤ 200

Now if we re-solve the problem, we obtain a FEASIBLE_POINT for the dual:

optimize!(model)
assert_is_solved_and_feasible(model)
dual_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

and a marginal price of electricity of $100/MWh:

dual(power_balance)

100.0

To reset the problem back to a mixed-integer linear program, we need to unfix and call set_binary:

unfix.(dispatch)
set_binary.(dispatch)
print(model)

Min 1000 dispatch[1] + 50 g[1] + 100 g[2]
Subject to
 power_balance : g[1] + g[2] + w = 1500
 g[1] ≥ 0
 g[2] - 300 dispatch[2] ≥ 0
 g[1] - 1000 dispatch[1] ≤ 0
 g[2] - 1000 dispatch[2] ≤ 0
 g[1] ≥ 0
 g[2] ≥ 300
 w ≥ 0
 g[1] ≤ 1000
 g[2] ≤ 1000
 w ≤ 200
 dispatch[1] binary
 dispatch[2] binary

Use fix_discrete_variables

Manually choosing the variables to relax and fix works for our small example, but it becomes more difficult in problems with a larger number of binary and integer variables. To automate the process we just did manually, JuMP provides the fix_discrete_variables function:

optimize!(model)
assert_is_solved_and_feasible(model)
dual_status(model)

NO_SOLUTION::ResultStatusCode = 0

undo = fix_discrete_variables(model);

Here undo is a function that, when called with no arguments, returns the model to the original mixed-integer formulation.

print(model)

Min 1000 dispatch[1] + 50 g[1] + 100 g[2]
Subject to
 power_balance : g[1] + g[2] + w = 1500
 g[1] ≥ 0
 g[2] - 300 dispatch[2] ≥ 0
 g[1] - 1000 dispatch[1] ≤ 0
 g[2] - 1000 dispatch[2] ≤ 0
 dispatch[1] = 1
 dispatch[2] = 1
 g[1] ≥ 0
 g[2] ≥ 300
 w ≥ 0
 g[1] ≤ 1000
 g[2] ≤ 1000
 w ≤ 200

optimize!(model)
assert_is_solved_and_feasible(model)
dual_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

dual(power_balance)

100.0

Finally, call undo to revert the reformulation

undo()
print(model)

Min 1000 dispatch[1] + 50 g[1] + 100 g[2]
Subject to
 power_balance : g[1] + g[2] + w = 1500
 g[1] ≥ 0
 g[2] - 300 dispatch[2] ≥ 0
 g[1] - 1000 dispatch[1] ≤ 0
 g[2] - 1000 dispatch[2] ≤ 0
 g[1] ≥ 0
 g[2] ≥ 300
 w ≥ 0
 g[1] ≤ 1000
 g[2] ≤ 1000
 w ≤ 200
 dispatch[1] binary
 dispatch[2] binary