# 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 the duals of a mixed-integer linear program by fixing the discrete variables to their optimal solution and resolving as a linear program.

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 = 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.

Tip

After calling fix_discrete_variables, you can set a new solver with set_optimizer if your mixed-integer solver does not support computing a dual solution.

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