Смешанные задачи комплементарности

Это руководство представляет собой набор смешанных задач комплементарности.

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

using JuMP
import PATHSolver

Вводная информация

Смешанная задача комплементарности имеет следующую форму:

где ограничение обеспечивает следующие отношения:

Если , то .
Если , то .
Если , то .

Возможно, вам приходилось видеть такую запись задачи комплементарности: . Это частный случай смешанной задачи комплементарности, в котором и .

Важно отметить, что в смешанной задаче комплементарности нет целевой функции и ограничений других типов.

Линейная комплементарность

Смешанная задача комплементарности формулируется с помощью символа перпендикулярности ⟂ (в REPL введите \perp<tab>).

M = [0 0 -1 -1; 0 0 1 -2; 1 -1 2 -2; 1 2 -2 4]
q = [2, 2, -2, -6]
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, 0 <= x[1:4] <= 10, start = 0)
@constraint(model, M * x + q ⟂ x)
optimize!(model)
@assert is_solved_and_feasible(model)
value.(x)

4-element Vector{Float64}:
 2.8
 0.0
 0.8
 1.2

Другие способы записи линейных задач комплементарности

Вам не обязательно использовать один вектор переменных, а ограничения комплементарности могут задаваться в любом порядке. Кроме того, вы можете использовать символ перпендикулярности, синтаксис complements(F, x) или множество MOI.Complements.

model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, 0 <= w <= 10, start = 0)
@variable(model, 0 <= x <= 10, start = 0)
@variable(model, 0 <= y <= 10, start = 0)
@variable(model, 0 <= z <= 10, start = 0)
@constraint(model, complements(y - 2z + 2, x))
@constraint(model, [-y - z + 2, w] in MOI.Complements(2))
@constraint(model, w + 2x - 2y + 4z - 6 ⟂ z)
@constraint(model, w - x + 2y - 2z - 2 ⟂ y)
optimize!(model)
@assert is_solved_and_feasible(model)
value.([w, x, y, z])

4-element Vector{Float64}:
 2.8
 0.0
 0.8
 1.2

Перевозки

Этот пример представляет собой переформулировку задачи о перевозках из главы 3.3 книги Dantzig, G.B. (1963). Linear Programming and Extensions. Princeton University Press, Princeton, New Jersey. Он основан на модели GAMS gamslib_transmcp.

capacity = Dict("seattle" => 350, "san-diego" => 600)
demand = Dict("new-york" => 325, "chicago" => 300, "topeka" => 275)
cost = Dict(
    ("seattle" => "new-york") => 90 * 2.5 / 1_000,
    ("seattle" => "chicago") => 90 * 1.7 / 1_000,
    ("seattle" => "topeka") => 90 * 1.8 / 1_000,
    ("san-diego" => "new-york") => 90 * 2.5 / 1_000,
    ("san-diego" => "chicago") => 90 * 1.8 / 1_000,
    ("san-diego" => "topeka") => 90 * 1.4 / 1_000,
)
plants, markets = keys(capacity), keys(demand)
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, w[i in plants] >= 0)
@variable(model, p[j in markets] >= 0)
@variable(model, x[i in plants, j in markets] >= 0)
@constraints(
    model,
    begin
        [i in plants, j in markets], w[i] + cost[i=>j] - p[j] ⟂ x[i, j]
        [i in plants], capacity[i] - sum(x[i, :]) ⟂ w[i]
        [j in markets], sum(x[:, j]) - demand[j] ⟂ p[j]
    end
)
optimize!(model)
@assert is_solved_and_feasible(model)
value.(p)

1-dimensional DenseAxisArray{Float64,1,...} with index sets:
    Dimension 1, ["new-york", "chicago", "topeka"]
And data, a 3-element Vector{Float64}:
 0.22500000000033224
 0.15299999994933886
 0.126

Ожидаемая полезность страхования

Этот пример взят из лекции курса AAE706, который читал Томас Ф. Резерфорд (Thomas F. Rutherford) в Висконсинском университете в Мадисоне. В нем моделируется ожидаемое страховое покрытие K, которое рациональному субъекту следует выбрать для страхования риска, возникающего с вероятностью pi и приводящего к убыткам в размере L.

pi = 0.01  # Вероятность плохого исхода
L = 0.5    # Убытки при плохом исходе
γ = 0.02   # Премия за покрытие
σ = 0.5    # Эластичность
ρ = -1     # Экспонента риска
U(C) = C^ρ / ρ
MU(C) = C^(ρ - 1)
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, EU, start = 1)   # Ожидаемая полезность
@variable(model, EV, start = 1)   # Эквивалентное изменение дохода
@variable(model, C_G, start = 1)  # Потребление в хороший день
@variable(model, C_B, start = 1)  # Потребление в плохой день
@variable(model, K, start = 1)    # Покрытие
@constraints(
    model,
    begin
        (1 - pi) * U(C_G) + pi * U(C_B) - EU ⟂ EU
        100 * (((1 - pi) * C_G^ρ + pi * C_B^ρ)^(1 / ρ) - 1) - EV ⟂ EV
        1 - γ * K - C_G ⟂ C_G
        1 - L + (1 - γ) * K - C_B ⟂ C_B
        γ * ((1 - pi) * MU(C_G) + pi * MU(C_B)) - pi * MU(C_B) ⟂ K
    end
)
optimize!(model)
@assert is_solved_and_feasible(model)
value(K)

0.20474003534537774

Потребление электроэнергии

Этот пример представляет собой формулировку в виде смешанной задачи комплементарности примера 3.3.1 из работы DAertrycke2017.

В примере моделируется нейтральное к риску конкурентное равновесие между производителем и потребителем электроэнергии.

В примере предполагается, что производитель планирует инвестировать в новую электростанцию мощностью (МВт). Годовые капитальные затраты этого завода составляют (€/МВт), а эксплуатационные расходы — (€/МВт·ч). В году 8760 часов.

После капитальных вложений есть пять возможных сценариев потребления , которые произойдут с вероятностью . В каждом сценарии производитель производит МВт электроэнергии.

В модели имеется один потребитель с квадратичной функцией полезности .

Теперь мы построим и решим смешанную задачу комплементарности, дав ряд кратких комментариев. Экономическое обоснование модели потребовало бы более объемного руководства. Подробные сведения см. в работе DAertrycke2017.

I = 90_000                     # Годовые капитальные затраты
C = 60                         # Эксплуатационные расходы на МВт·ч
τ = 8_760                      # Количество часов в году
θ = [0.2, 0.2, 0.2, 0.2, 0.2]  # Вероятности сценариев
A = [300, 350, 400, 450, 500]  # Коэффициенты функции полезности
B = 1                          # Коэффициенты функции полезности
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, x >= 0, start = 1)           # Установленная мощность
@variable(model, Q[ω = 1:5] >= 0, start = 1)  # Потребление
@variable(model, Y[ω = 1:5] >= 0, start = 1)  # Производство
@variable(model, P[ω = 1:5], start = 1)       # Цена электроэнергии
@variable(model, μ[ω = 1:5] >= 0, start = 1)  # Маржа дефицита капитала
# Стоимость паевых инвестиций равна годовому дефициту прибыли, либо инвестиции равны 0
@constraint(model, I - τ * θ' * μ ⟂ x)
# Разница между ценой и маржой дефицита равна эксплуатационным расходам
@constraint(model, [ω = 1:5], C - (P[ω] - μ[ω]) ⟂ Y[ω])
# Цена равна предельной полезности для потребителя
@constraint(model, [ω = 1:5], P[ω] - (A[ω] - B * Q[ω]) ⟂ Q[ω])
# Производство равно потреблению
@constraint(model, [ω = 1:5], Y[ω] - Q[ω] ⟂ P[ω])
# Производство не превышает мощности
@constraint(model, [ω = 1:5], x - Y[ω] ⟂ μ[ω])
optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model)

* Solver : Path 5.0.03

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "The problem was solved"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : 0.00000e+00

* Work counters
  Solve time (sec)   : 2.07000e-04

Равновесным решением является мощность 389 МВт:

value(x)

389.31506849315065

Производство в каждом сценарии:

value.(Q)

5-element Vector{Float64}:
 240.0000000000001
 289.9999999999999
 340.0
 389.31506849315065
 389.31506849315065

Цена в каждом сценарии:

value.(P)

5-element Vector{Float64}:
  59.999999999999886
  60.0
  59.99999999999994
  60.68493150684928
 110.68493150684935