Нелинейная оптимизация с ограничениями

Интерфейс нелинейной оптимизации с ограничениями в Optim предполагает, что пользователь может написать задачу оптимизации следующим образом.

Для ограничений в виде равенства для или вы устанавливаете эти конкретные записи границ равными, . Аналогично, задание или означает, что ограничение не ограничено снизу или сверху, соответственно.

Оптимизация с ограничениями с `IPNewton`

Мы рассмотрим примеры использования интерфейса ограничений с алгоритмом оптимизации Ньютона с внутренней точкой IPNewton.

В этих примерах мы работаем со стандартной функцией Розенброка. Цель и ее производные определяются следующим образом:

fun(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end

function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;

Интерфейс оптимизации

Для решения задачи оптимизации с ограничениями мы вызываем метод optimize.

optimize(d::AbstractObjective, constraints::AbstractConstraints, initial_x::Tx, method::ConstrainedOptimizer, options::Options)

Мы можем создать экземпляры AbstractObjective и AbstractConstraints, используя типы TwiceDifferentiable и TwiceDifferentiableConstraints из пакета NLSolversBase.jl.

Минимизация блока

Мы хотим оптимизировать функцию Розенброка в блоке , начиная с точки . Ограничения блока заданы, например, с помощью TwiceDifferentiableConstraints(lx, ux).

x0 = [0.0, 0.0]
df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)

lx = [-0.5, -0.5]; ux = [0.5, 0.5]
dfc = TwiceDifferentiableConstraints(lx, ux)

res = optimize(df, dfc, x0, IPNewton())

 * Status: success

 * Candidate solution
    Final objective value:     2.500000e-01

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 4.39e-10 ≰ 0.0e+00
    |x - x'|/|x'|          = 8.79e-10 ≰ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 1.00e+00 ≰ 1.0e-08

 * Work counters
    Seconds run:   1  (vs limit Inf)
    Iterations:    43
    f(x) calls:    68
    ∇f(x) calls:   68

Как и в случае с остальной частью Optim, вы также можете использовать autodiff=:forward и просто передать fun.

Чтобы задать только нижние границы, используйте ux = fill(Inf, 2).

ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

 * Status: success (objective increased between iterations)

 * Candidate solution
    Final objective value:     7.987239e-20

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 3.54e-10 ≰ 0.0e+00
    |x - x'|/|x'|          = 3.54e-10 ≰ 0.0e+00
    |f(x) - f(x')|         = 2.40e-19 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 3.00e+00 ≰ 0.0e+00
    |g(x)|                 = 8.83e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    35
    f(x) calls:    63
    ∇f(x) calls:   63

Определение задач без ограничений

Задачу без ограничений можно определить, передав либо границы Inf, либо пустые массивы. Обратите внимание, что мы должны передать правильную информацию о типе в пустые lx и ux

lx = fill(-Inf, 2); ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

lx = Float64[]; ux = Float64[]
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

 * Status: success

 * Candidate solution
    Final objective value:     5.998937e-19

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 1.50e-09 ≰ 0.0e+00
    |x - x'|/|x'|          = 1.50e-09 ≰ 0.0e+00
    |f(x) - f(x')|         = 1.80e-18 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 3.00e+00 ≰ 0.0e+00
    |g(x)|                 = 7.92e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    34
    f(x) calls:    63
    ∇f(x) calls:   63

Универсальные нелинейные ограничения

Теперь рассмотрим задачу Розенброка с ограничением в

Мы передаем информацию об ограничениях в optimize, определяя вектор-функцию c(x) и ее якобиан J(x).

Информация о гессиане обрабатывается по-другому, и при этом принимается во внимание лагранжиан соответствующей преобразованной задачи оптимизации со слабыми переменными. Процесс похож на работу библиотеки CUTEst. Пусть представляет собой гессиан -й компоненты универсальных ограничений, а — соответствующую двойственную переменную в лагранжиане. Затем нам нужно, чтобы объект constraint добавил значения к гессиану цели, взвешенному с помощью .

Далее форма Julia для заданной функции и информация о производной добавляются следующим образом.

con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)
function con_jacobian!(J, x)
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    J
end
function con_h!(h, x, λ)
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
end;

Обратите внимание, что con_h! добавляет λ-взвешенное значение гессиана каждого элемента c(x) к гессиану fun.

Затем мы можем оптимизировать функцию Розенброка внутри шара радиусом .

lx = Float64[]; ux = Float64[]
lc = [-Inf]; uc = [0.5^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

 * Status: success

 * Candidate solution
    Final objective value:     2.966216e-01

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 7.71e-01 ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    32
    f(x) calls:    111
    ∇f(x) calls:   111

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

lc = [0.1^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

 * Status: success

 * Candidate solution
    Final objective value:     2.966216e-01

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 7.71e-01 ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    36
    f(x) calls:    162
    ∇f(x) calls:   162

Обратите внимание, что алгоритм предупреждает, что начальное предположение не является внутренней точкой. Часто IPNewton может справиться с этим, однако, если начальная догадка такова, что c(x) = u_c, то алгоритм в настоящий момент завершается сбоем. Возможно, в будущем мы это исправим.

Несколько ограничений

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

function con2_c!(c, x)
    c[1] = x[1]^2 + x[2]^2     ## Первое ограничение
    c[2] = x[2]*sin(x[1])-x[1] ## Второе ограничение
    c
end
function con2_jacobian!(J, x)
    # Первое ограничение
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    # Второе ограничение
    J[2,1] = x[2]*cos(x[1])-1.0
    J[2,2] = sin(x[1])
    J
end
function con2_h!(h, x, λ)
    # Первое ограничение
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
    # Второе ограничение
    h[1,1] += λ[2]*x[2]*-sin(x[1])
    h[1,2] += λ[2]*cos(x[1])
    # Симметризировать h
    h[2,1]  = h[1,2]
    h
end;

Мы создаем объекты ограничений и вызываем IPNewton с начальным предположением .

x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
dfc = TwiceDifferentiableConstraints(con2_c!, con2_jacobian!, con2_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

 * Status: success

 * Candidate solution
    Final objective value:     1.000000e+00

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 6.90e-10 ≰ 0.0e+00
    |x - x'|/|x'|          = 3.55e+08 ≰ 0.0e+00
    |f(x) - f(x')|         = 1.38e-09 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.38e-09 ≰ 0.0e+00
    |g(x)|                 = 2.00e+00 ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    29
    f(x) calls:    219
    ∇f(x) calls:   219

Простая программа

Ниже приводится версия программы без комментариев. Файл также доступен здесь: ipnewton_basics.jl

using Optim, NLSolversBase #hide
import NLSolversBase: clear! #hide

fun(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end

function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;

x0 = [0.0, 0.0]
df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)

lx = [-0.5, -0.5]; ux = [0.5, 0.5]
dfc = TwiceDifferentiableConstraints(lx, ux)

res = optimize(df, dfc, x0, IPNewton())

ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

lx = fill(-Inf, 2); ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

lx = Float64[]; ux = Float64[]
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)
function con_jacobian!(J, x)
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    J
end
function con_h!(h, x, λ)
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
end;

lx = Float64[]; ux = Float64[]
lc = [-Inf]; uc = [0.5^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

lc = [0.1^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

function con2_c!(c, x)
    c[1] = x[1]^2 + x[2]^2     ## Первое ограничение
    c[2] = x[2]*sin(x[1])-x[1] ## Второе ограничение
    c
end
function con2_jacobian!(J, x)
    # Первое ограничение
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    # Второе ограничение
    J[2,1] = x[2]*cos(x[1])-1.0
    J[2,2] = sin(x[1])
    J
end
function con2_h!(h, x, λ)
    # Первое ограничение
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
    # Второе ограничение
    h[1,1] += λ[2]*x[2]*-sin(x[1])
    h[1,2] += λ[2]*cos(x[1])
    # Симметризировать h
    h[2,1]  = h[1,2]
    h
end;

x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
dfc = TwiceDifferentiableConstraints(con2_c!, con2_jacobian!, con2_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

# Этот файл был создан с помощью Literate.jl, https://github.com/fredrikekre/Literate.jl.

Эта страница была создана с помощью Literate.jl.