Нахождение максимумов и минимумов решений обыкновенных дифференциальных уравнение

Настройка

В этом руководстве мы покажем, как с помощью пакета Optimization.jl находить максимумы и минимумы решения. Давайте рассмотрим двойной маятник.

#Константы и настройка
using OrdinaryDiffEq
initial = [0.01, 0.01, 0.01, 0.01]
tspan = (0.0, 100.0)

#Определение задачи
function double_pendulum_hamiltonian(udot, u, p, t)
    α = u[1]
    lα = u[2]
    β = u[3]
    lβ = u[4]
    udot .= [2(lα - (1 + cos(β))lβ) / (3 - cos(2β)),
        -2sin(α) - sin(α + β),
        2(-(1 + cos(β))lα + (3 + 2cos(β))lβ) / (3 - cos(2β)),
        -sin(α + β) - 2sin(β) * (((lα - lβ)lβ) / (3 - cos(2β))) +
        2sin(2β) * ((lα^2 - 2(1 + cos(β))lα * lβ + (3 + 2cos(β))lβ^2) / (3 - cos(2β))^2)]
end

#Передача в решатели
poincare = ODEProblem(double_pendulum_hamiltonian, initial, tspan)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 4-element Vector{Float64}:
 0.01
 0.01
 0.01
 0.01

sol = solve(poincare, Tsit5())

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 193-element Vector{Float64}:
   0.0
   0.08332584852065579
   0.24175300587841853
   0.4389533535703127
   0.6797301355043014
   0.9647629844957191
   1.3179426427778789
   1.7031226589934565
   2.0678504564408504
   2.4717900471019023
   ⋮
  95.8457265154543
  96.35778609875159
  96.92912892765023
  97.44679036707718
  97.96248017978454
  98.5118292417693
  99.06082273312332
  99.58284082787185
 100.0
u: 193-element Vector{Vector{Float64}}:
 [0.01, 0.01, 0.01, 0.01]
 [0.00917068738040534, 0.0066690004553842845, 0.012420525490765832, 0.008266408515192912]
 [0.007673275265972518, 0.00037461737897660346, 0.016442590227730366, 0.004636827483318282]
 [0.0061259744192393014, -0.007305450189721184, 0.01996737108423187, -0.00033649798308967233]
 [0.0049661106627111465, -0.01630851653373806, 0.021440659476204688, -0.006705037098400459]
 [0.004795568366455252, -0.026238104231339168, 0.01882432482705364, -0.01391336508452692]
 [0.006054680216573254, -0.03712455410888037, 0.010055700343225163, -0.021038128751556112]
 [0.007900784624609724, -0.04667607081473084, -0.002673583699559344, -0.025183036571892782]
 [0.008276510353924267, -0.052784334378941124, -0.01273154768229481, -0.02525804011008269]
 [0.005523496102562292, -0.055252504127907325, -0.01684388181977143, -0.021898962485296523]
 ⋮
 [-0.014886778911438145, 0.042332620369772034, 0.013628258608953085, 0.01802907676255885]
 [-0.008190351596162247, 0.05442260545036836, 0.009448112771214356, 0.017740074304871352]
 [0.00412458896341776, 0.05674882926735902, -0.005154041839394905, 0.01759697731920902]
 [0.013079674415301898, 0.04807713954653274, -0.013770643022910575, 0.018286645951147314]
 [0.015316052216181856, 0.03163113020569826, -0.008957094529053966, 0.017118433263410053]
 [0.011115540528354754, 0.009929208333814447, 0.007297331934230443, 0.010353457826205819]
 [0.005713897643562407, -0.011787329261525195, 0.020508032523607597, -0.002310390783307738]
 [0.004211435541855515, -0.02991110357130916, 0.018750502558174825, -0.015650642104010095]
 [0.0057412395742532365, -0.04165385988427417, 0.007413270264835161, -0.023348978443051394]

Во времени решение выглядит следующим образом.

using Plots;
gr();
plot(sol, vars = [(0, 3), (0, 4)], leg = false, plotdensity = 10000)

Оно имеет известный график фазового пространства:

plot(sol, vars = (3, 4), leg = false)

Локальная оптимизация

Давайте найдем локальные максимумы и минимумы. С помощью Optim.jl можно минимизировать функции, причем тип решения имеет непрерывную интерполяцию, чем можно воспользоваться. Давайте найдем локальные оптимумы 4-й переменной в окрестности точки t=20. Функция оптимизации имеет следующий вид:

f(t, _) = sol(first(t), idxs = 4)

f (generic function with 1 method)

Выражение first(t) равносильно t[1] и преобразует массив размера 1 в число. Выражение idxs=4 равносильно sol(first(t))[4], но вычисление выполняется без временного массива, и поэтому производится быстрее. Чтобы найти локальный минимум, можно решить задачу оптимизации с функцией потерь f:

using Optimization, OptimizationNLopt, ForwardDiff
optf = OptimizationFunction(f, Optimization.AutoForwardDiff())
min_guess = 18.0
optprob = OptimizationProblem(optf, [min_guess], lb = [0.0], ub = [100.0])
opt = solve(optprob, NLopt.LD_LBFGS())

u: 1-element Vector{Float64}:
 18.63212745189591

Из этих выходных данных видно, что минимум находится в точке t=18.63 и имеет значение -2.79e-2. С помощью кода этот результат можно получить так:

println(opt.u)

[18.63212745189591]

Чтобы получить максимум, мы просто минимизируем функцию с обратным знаком:

fminus(t, _) = -sol(first(t), idxs = 4)

optf = OptimizationFunction(fminus, Optimization.AutoForwardDiff())
min_guess = 22.0
optprob2 = OptimizationProblem(optf, [min_guess], lb = [0.0], ub = [100.0])
opt2 = solve(optprob2, NLopt.LD_LBFGS())

u: 1-element Vector{Float64}:
 23.297606513035095

Давайте добавим максимумы и минимумы на графики:

plot(sol, vars = (0, 4), plotdensity = 10000)
scatter!([opt.u], [opt.minimum], label = "Local Min")
scatter!([opt2.u], [-opt2.minimum], label = "Local Max")

Глобальная оптимизация

Если вместо этого нужно найти глобальные максимумы и минимумы, потребуется иной подход. Вариантов много. В Julia для этого есть свое собственное средство — решатели BlackBoxOptim из пакета Optimization.jl. Однако мы вместо этого воспользуемся методами OptimizationNLopt. Для этого мы просто переключимся на один из глобальных оптимизаторов из списка. Давайте попробуем GN_ORIG_DIRECT_L:

gopt = solve(optprob, NLopt.GN_ORIG_DIRECT_L())
gopt2 = solve(optprob2, NLopt.GN_ORIG_DIRECT_L())

@show gopt.u, gopt2.u

([76.29050924282319], [47.4603033758301])

plot(sol, vars = (0, 4), plotdensity = 10000)
scatter!([gopt.u], [gopt.minimum], label = "Global Min")
scatter!([gopt2.u], [-gopt2.minimum], label = "Global Max")