# Computing the duals of a mixed-integer program

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.0
g[1] ≥ 0.0
g[2] - 300 dispatch[2] ≥ 0.0
g[1] - 1000 dispatch[1] ≤ 0.0
g[2] - 1000 dispatch[2] ≤ 0.0
g[1] ≥ 0.0
g[2] ≥ 300.0
w ≥ 0.0
g[1] ≤ 1000.0
g[2] ≤ 1000.0
w ≤ 200.0
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)
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.0
g[1] ≥ 0.0
g[2] - 300 dispatch[2] ≥ 0.0
g[1] - 1000 dispatch[1] ≤ 0.0
g[2] - 1000 dispatch[2] ≤ 0.0
dispatch[1] = 1.0
dispatch[2] = 1.0
g[1] ≥ 0.0
g[2] ≥ 300.0
w ≥ 0.0
g[1] ≤ 1000.0
g[2] ≤ 1000.0
w ≤ 200.0
```

Now if we re-solve the problem, we obtain a `FEASIBLE_POINT`

for the dual:

```
optimize!(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.0
g[1] ≥ 0.0
g[2] - 300 dispatch[2] ≥ 0.0
g[1] - 1000 dispatch[1] ≤ 0.0
g[2] - 1000 dispatch[2] ≤ 0.0
g[1] ≥ 0.0
g[2] ≥ 300.0
w ≥ 0.0
g[1] ≤ 1000.0
g[2] ≤ 1000.0
w ≤ 200.0
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)
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.

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.0
g[1] ≥ 0.0
g[2] - 300 dispatch[2] ≥ 0.0
g[1] - 1000 dispatch[1] ≤ 0.0
g[2] - 1000 dispatch[2] ≤ 0.0
dispatch[1] = 1.0
dispatch[2] = 1.0
g[1] ≥ 0.0
g[2] ≥ 300.0
w ≥ 0.0
g[1] ≤ 1000.0
g[2] ≤ 1000.0
w ≤ 200.0
```

```
optimize!(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.0
g[1] ≥ 0.0
g[2] - 300 dispatch[2] ≥ 0.0
g[1] - 1000 dispatch[1] ≤ 0.0
g[2] - 1000 dispatch[2] ≤ 0.0
g[1] ≥ 0.0
g[2] ≥ 300.0
w ≥ 0.0
g[1] ≤ 1000.0
g[2] ≤ 1000.0
w ≤ 200.0
dispatch[1] binary
dispatch[2] binary
```

This tutorial was generated using Literate.jl. View the source `.jl`

file on GitHub.