Решение задачи Розенброка более чем 10 способами

Страница в процессе перевода.

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

Note В этом примере используется множество различных решателей из Optimization.jl. Подпакет каждого решателя необходимо устанавливать отдельно. Например, подробные сведения об установке и использовании пакета OptimizationOptimJL.jl см. на странице Optim.jl.

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

Интерфейс Optimization.jl предполагает, что функции определены с вектором аргументов оптимизации и вектором параметров , то есть:

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

Исходные области определения и записываются в вектор .

Для инициализации оптимизатора требуется начальная оценка местоположения минимума.

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

# Определяем решаемую задачу
using Optimization, ForwardDiff, Zygote

rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p = [1.0, 100.0]

f = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff())
l1 = rosenbrock(x0, _p)
prob = OptimizationProblem(f, x0, _p)

OptimizationProblem. In-place: true
u0: 2-element Vector{Float64}:
 0.0
 0.0

Решатели Optim.jl

Первый шаг: оптимизаторы без вычисления производных

using OptimizationOptimJL
sol = solve(prob, SimulatedAnnealing())
prob = OptimizationProblem(f, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8])
sol = solve(prob, SAMIN())

l1 = rosenbrock(x0, _p)
prob = OptimizationProblem(rosenbrock, x0, _p)
sol = solve(prob, NelderMead())

retcode: Success
u: 2-element Vector{Float64}:
 0.9999634355313174
 0.9999315506115275

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

optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff())
prob = OptimizationProblem(optf, x0, _p)
sol = solve(prob, BFGS())

retcode: Success
u: 2-element Vector{Float64}:
 0.9999999999373603
 0.99999999986862

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

sol = solve(prob, Newton())

retcode: Success
u: 2-element Vector{Float64}:
 0.9999999999999994
 0.9999999999999989

Оптимизатор второго порядка без гессиана

sol = solve(prob, Optim.KrylovTrustRegion())

retcode: Success
u: 2-element Vector{Float64}:
 0.999999999999108
 0.9999999999981819

Оптимизаторы на основе производных с различными ограничениями

cons = (res, x, p) -> res .= [x[1]^2 + x[2]^2]
optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); cons = cons)

prob = OptimizationProblem(optf, x0, _p, lcons = [-Inf], ucons = [Inf])
sol = solve(prob, IPNewton()) # Обратите внимание, что -Inf < x[1]^2 + x[2]^2 < Inf всегда верно

prob = OptimizationProblem(optf, x0, _p, lcons = [-5.0], ucons = [10.0])
sol = solve(prob, IPNewton()) # Опять-таки -5.0 < x[1]^2 + x[2]^2 < 10.0

prob = OptimizationProblem(optf, x0, _p, lcons = [-Inf], ucons = [Inf],
    lb = [-500.0, -500.0], ub = [50.0, 50.0])
sol = solve(prob, IPNewton())

prob = OptimizationProblem(optf, x0, _p, lcons = [0.5], ucons = [0.5],
    lb = [-500.0, -500.0], ub = [50.0, 50.0])
sol = solve(prob, IPNewton())

# Теперь обратите внимание, что x[1]^2 + x[2]^2 ≈ 0.5:
res = zeros(1)
cons(res, sol.u, _p)
println(res)

┌ Warning: The selected optimization algorithm requires second order derivatives, but `SecondOrder` ADtype was not provided.
│         So a `SecondOrder` with AutoForwardDiff() for both inner and outer will be created, this can be suboptimal and not work in some cases so
│         an explicit `SecondOrder` ADtype is recommended.
└ @ OptimizationBase ~/work/package/Optimization.jl/lib/OptimizationBase/src/cache.jl:51
┌ Warning: The selected optimization algorithm requires second order derivatives, but `SecondOrder` ADtype was not provided.
│         So a `SecondOrder` with AutoForwardDiff() for both inner and outer will be created, this can be suboptimal and not work in some cases so
│         an explicit `SecondOrder` ADtype is recommended.
└ @ OptimizationBase ~/work/package/Optimization.jl/lib/OptimizationBase/src/cache.jl:51
┌ Warning: The selected optimization algorithm requires second order derivatives, but `SecondOrder` ADtype was not provided.
│         So a `SecondOrder` with AutoForwardDiff() for both inner and outer will be created, this can be suboptimal and not work in some cases so
│         an explicit `SecondOrder` ADtype is recommended.
└ @ OptimizationBase ~/work/package/Optimization.jl/lib/OptimizationBase/src/cache.jl:51
┌ Warning: The selected optimization algorithm requires second order derivatives, but `SecondOrder` ADtype was not provided.
│         So a `SecondOrder` with AutoForwardDiff() for both inner and outer will be created, this can be suboptimal and not work in some cases so
│         an explicit `SecondOrder` ADtype is recommended.
└ @ OptimizationBase ~/work/package/Optimization.jl/lib/OptimizationBase/src/cache.jl:51
[0.49999999999999994]

function con_c(res, x, p)
    res .= [x[1]^2 + x[2]^2]
end

optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); cons = con_c)
prob = OptimizationProblem(optf, x0, _p, lcons = [-Inf], ucons = [0.25^2])
sol = solve(prob, IPNewton()) # -Inf < cons_circ(sol.u, _p) = 0.25^2

retcode: Success
u: 2-element Vector{Float64}:
 0.24327905408863862
 0.05757865786675858

Решатели Evolutionary.jl

using OptimizationEvolutionary
sol = solve(prob, CMAES(μ = 40, λ = 100), abstol = 1e-15) # -Inf < cons_circ(sol.u, _p) = 0.25^2

retcode: Success
u: 2-element Vector{Float64}:
 0.24335547167668792
 0.057254819928253604

IPOPT через OptimizationMOI

using OptimizationMOI, Ipopt

function con2_c(res, x, p)
    res .= [x[1]^2 + x[2]^2, x[2] * sin(x[1]) - x[1]]
end

optf = OptimizationFunction(rosenbrock, Optimization.AutoZygote(); cons = con2_c)
prob = OptimizationProblem(optf, x0, _p, lcons = [-Inf, -Inf], ucons = [100.0, 100.0])
sol = solve(prob, Ipopt.Optimizer())

retcode: Success
u: 2-element Vector{Float64}:
 0.9999999999080653
 0.9999999998157655

Теперь перейдем к OptimizationOptimisers с автоматическим дифференцированием в обратном режиме

import OptimizationOptimisers
optf = OptimizationFunction(rosenbrock, Optimization.AutoZygote())
prob = OptimizationProblem(optf, x0, _p)
sol = solve(prob, OptimizationOptimisers.Adam(0.05), maxiters = 1000, progress = false)

retcode: Default
u: 2-element Vector{Float64}:
 0.999957450694706
 0.9999163934196296

Пробуем эволюционные методы из CMAEvolutionStrategy.jl

using OptimizationCMAEvolutionStrategy
sol = solve(prob, CMAEvolutionStrategyOpt())

retcode: Failure
u: 2-element Vector{Float64}:
 1.000000137136633
 1.0000002638470338

Теперь пробуем несколько решателей NLopt.jl с символьным дифференцированием через ModelingToolkit.jl

using OptimizationNLopt, ModelingToolkit
optf = OptimizationFunction(rosenbrock, Optimization.AutoSymbolics())
prob = OptimizationProblem(optf, x0, _p)

sol = solve(prob, Opt(:LN_BOBYQA, 2))
sol = solve(prob, Opt(:LD_LBFGS, 2))

retcode: Success
u: 2-element Vector{Float64}:
 0.9999999999894374
 0.9999999999844783

Добавляем прямоугольные ограничения и решаем задачу несколькими методами из NLopt.jl

prob = OptimizationProblem(optf, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8])
sol = solve(prob, Opt(:LD_LBFGS, 2))
sol = solve(prob, Opt(:G_MLSL_LDS, 2), local_method = Opt(:LD_LBFGS, 2), maxiters = 10000) #глобальный оптимизатор со случайными начальными точками локальной оптимизации

retcode: MaxIters
u: 2-element Vector{Float64}:
 0.8
 0.6400000000000001

Решатели BlackBoxOptim.jl

using OptimizationBBO
prob = Optimization.OptimizationProblem(rosenbrock, [0.0, 0.3], _p, lb = [-1.0, 0.2],
    ub = [0.8, 0.43])
sol = solve(prob, BBO_adaptive_de_rand_1_bin()) # -1.0 ≤ x[1] ≤ 0.8, 0.2 ≤ x[2] ≤ 0.43

retcode: MaxIters
u: 2-element Vector{Float64}:
 0.6577248384964002
 0.43

И это лишь малая часть того, что может предложить Optimization.jl!