Документация Engee

Вычисление двойственных решений частично целочисленной программы

В этом руководстве объясняется, как вычислить двойственные решения частично целочисленной линейной программы путем фиксации дискретных переменных в точке оптимального решения и решения линейной программы.

В этом руководстве используются следующие пакеты.

using JuMP
import HiGHS

Модель

В качестве модели взят пример составления графика нагрузки энергоблоков из раздела Unit commitment. Подробности не важны, за исключением того, что есть два типа непрерывных переменных, g и w, представляющих количество вырабатываемой тепловыми и ветровыми электростанциями энергии, и дискретная переменная dispatch, которая равна 1, если электростанция i работает, и 0 в противном случае.

Нас интересует «двойственное» решение ограничения power_balance, так как оно представляет предельную цену потребляемой электроэнергии.

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

Фиксация переменных вручную

Если оптимизировать эту модель, мы получим значение NO_SOLUTION метода dual_status:

optimize!(model)
@assert is_solved_and_feasible(model)
dual_status(model)
NO_SOLUTION::ResultStatusCode = 0

Причина в том, что HiGHS не может вычислять двойственные решения частично целочисленной программы. Эту проблему можно обойти, зафиксировав целочисленные переменные в точке оптимального решения, ослабив целочисленность и повторно решив задачу как линейную программу.

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

Если теперь решить задачу заново, мы получим FEASIBLE_POINT для двойственного решения:

optimize!(model)
@assert is_solved_and_feasible(model)
dual_status(model)
FEASIBLE_POINT::ResultStatusCode = 1

и предельную цену электроэнергии 100 долл./МВт·ч:

dual(power_balance)
100.0

Чтобы вернуться к частично целочисленной линейной программе, нужно вызвать unfix и 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

Использование fix_discrete_variables

Выбор переменных вручную для ослабления и фиксации подходит для нашего небольшого примера, но затруднительнее в задачах с большим количеством двоичных и целочисленных переменных. Для автоматизации процесса, который мы только что выполнили вручную, JuMP предоставляет функцию fix_discrete_variables:

optimize!(model)
@assert is_solved_and_feasible(model)
dual_status(model)
NO_SOLUTION::ResultStatusCode = 0
undo = fix_discrete_variables(model);

Здесь undo — это функция, которая при вызове без аргументов возвращает модель к исходной частично целочисленной формулировке.

После вызова fix_discrete_variables можно задать новый решатель с помощью set_optimizer, если ваш частично целочисленный решатель не поддерживает вычисление двойственного решения.

 
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

Наконец, вызовем undo, чтобы отменить переформулировку:

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