Вложенные задачи оптимизации

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

Для моделирования задачи определяется пользовательский оператор, который выполняет разложение нижней задачи внутри верхней. Наконец, мы покажем, как повысить производительность за счет кэша, который позволяет избежать решения нижней задачи.

Более простой пример написания пользовательского оператора см. в руководстве User-defined Hessians.

Для моделирования и решения двухуровневых задач оптимизации также можно использовать расширение JuMP BilevelJuMP.jl.

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

using JuMP
import Ipopt

В оставшейся части руководства нашей целью будет решение двухуровневой задачи оптимизации:

Эта двухуровневая задача состоит из двух вложенных задач оптимизации: верхнего уровня, включающего переменные , и нижнего уровня, включающего переменные . С точки зрения задачи нижнего уровня значения являются фиксированными параметрами, поэтому модель оптимизирует с учетом этих фиксированных параметров. В то же время задача верхнего уровня оптимизирует и с учетом отклика .

Разложение

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

function solve_lower_level(x...)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, y[1:2])
    @objective(
        model,
        Max,
        x[1]^2 * y[1] + x[2]^2 * y[2] - x[1] * y[1]^4 - 2 * x[2] * y[2]^4,
    )
    @constraint(model, (y[1] - 10)^2 + (y[2] - 10)^2 <= 25)
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return objective_value(model), value.(y)
end

solve_lower_level (generic function with 1 method)

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

function V(x...)
    f, _ = solve_lower_level(x...)
    return f
end

V (generic function with 1 method)

Затем можно подставить в полную задачу, чтобы получилось следующее.

Это похоже на задачу нелинейной оптимизации с пользовательским оператором ! Однако поскольку решает задачу оптимизации на внутреннем уровне, мы не можем использовать автоматическое дифференцирование для вычисления первой и второй производных. Вместо этого можно прибегнуть к имеющейся в JuMP возможности передавать функции обратного вызова для градиента и гессиана.

Сначала нужно определить градиент относительно . В общем случае вычислить его может быть трудно, но поскольку встречается только в целевой функции, мы можем просто дифференцировать целевую функцию по , получив следующее.

function ∇V(g::AbstractVector, x...)
    _, y = solve_lower_level(x...)
    g[1] = 2 * x[1] * y[1] - y[1]^4
    g[2] = 2 * x[2] * y[2] - 2 * y[2]^4
    return
end

∇V (generic function with 1 method)

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

function ∇²V(H::AbstractMatrix, x...)
    _, y = solve_lower_level(x...)
    H[1, 1] = 2 * y[1]
    H[2, 2] = 2 * y[2]
    return
end

∇²V (generic function with 1 method)

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

Теперь есть все необходимое для определения двухуровневой задачи оптимизации:

model = Model(Ipopt.Optimizer)
@variable(model, x[1:2] >= 0)
@operator(model, op_V, 2, V, ∇V, ∇²V)
@objective(model, Min, x[1]^2 + x[2]^2 + op_V(x[1], x[2]))
optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model)

* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : -4.18983e+05
  Dual objective value : 0.00000e+00

* Work counters
  Solve time (sec)   : 1.02851e+00
  Barrier iterations : 34

Оптимальное целевое значение:

objective_value(model)

-418983.4868064087

Оптимальные переменные решения верхнего уровня :

value.(x)

2-element Vector{Float64}:
 154.97862339279123
 180.00961428934832

Чтобы вычислить оптимальные переменные решения нижнего уровня, нужно вызвать solve_lower_level с оптимальными переменными решения верхнего уровня:

_, y = solve_lower_level(value.(x)...)
y

2-element Vector{Float64}:
 7.0725939609898285
 5.946569892949626

Повышение производительности

Такой подход работает, но у него есть проблема с производительностью: каждый раз, когда нужно вычислить значение, градиент или гессиан , приходится заново решать задачу оптимизации нижнего уровня. На это расходуются ресурсы, так как градиент и гессиан часто будут вызываться в одной точке, а значит, решение задачи дважды с одними и теми же входными данными приводит к ненужным повторениям.

Эту проблему можно обойти, используя кэш:

mutable struct Cache
    x::Any
    f::Float64
    y::Vector{Float64}
end

с функцией обновления кэша при необходимости:

function _update_if_needed(cache::Cache, x...)
    if cache.x !== x
        cache.f, cache.y = solve_lower_level(x...)
        cache.x = x
    end
    return
end

_update_if_needed (generic function with 1 method)

Затем мы определяем кэшированные версии наших трех функций, которые вызывают _updated_if_needed и возвращают значения из кэша.

function cached_f(cache::Cache, x...)
    _update_if_needed(cache, x...)
    return cache.f
end

function cached_∇f(cache::Cache, g::AbstractVector, x...)
    _update_if_needed(cache, x...)
    g[1] = 2 * x[1] * cache.y[1] - cache.y[1]^4
    g[2] = 2 * x[2] * cache.y[2] - 2 * cache.y[2]^4
    return
end

function cached_∇²f(cache::Cache, H::AbstractMatrix, x...)
    _update_if_needed(cache, x...)
    H[1, 1] = 2 * cache.y[1]
    H[2, 2] = 2 * cache.y[2]
    return
end

cached_∇²f (generic function with 1 method)

Теперь все готово для настройки и решения задачи оптимизации верхнего уровня:

model = Model(Ipopt.Optimizer)
@variable(model, x[1:2] >= 0)
cache = Cache(Float64[], NaN, Float64[])
@operator(
    model,
    op_cached_f,
    2,
    (x...) -> cached_f(cache, x...),
    (g, x...) -> cached_∇f(cache, g, x...),
    (H, x...) -> cached_∇²f(cache, H, x...),
)
@objective(model, Min, x[1]^2 + x[2]^2 + op_cached_f(x[1], x[2]))
optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model)

* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : -4.18983e+05
  Dual objective value : 0.00000e+00

* Work counters
  Solve time (sec)   : 2.68806e-01
  Barrier iterations : 34

Можно проверить, получается ли то же самое целевое значение:

objective_value(model)

-418983.4868064087

и переменная решения верхнего уровня :

value.(x)

2-element Vector{Float64}:
 154.97862339279123
 180.00961428934832