Finding Maxima and Minima of ODEs Solutions

Setup

In this tutorial, we will show how to use Optimization.jl to find the maxima and minima of solutions. Let’s take a look at the double pendulum:

#Constants and setup
using OrdinaryDiffEq
initial = [0.01, 0.01, 0.01, 0.01]
tspan = (0.0, 100.0)

#Define the problem
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

#Pass to solvers
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]

In time, the solution looks like:

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

while it has the well-known phase-space plot:

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

Local Optimization

Let’s find out what some of the local maxima and minima are. Optim.jl can be used to minimize functions, and the solution type has a continuous interpolation which can be used. Let’s look for the local optima for the 4th variable around t=20. Thus, our optimization function is:

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

f (generic function with 1 method)

first(t) is the same as t[1] which transforms the array of size 1 into a number. idxs=4 is the same as sol(first(t))[4] but does the calculation without a temporary array and thus is faster. To find a local minimum, we can solve the optimization problem where the loss function is 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

From this printout, we see that the minimum is at t=18.63 and the value is -2.79e-2. We can get these in code-form via:

println(opt.u)

[18.63212745189591]

To get the maximum, we just minimize the negative of the function:

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

Let’s add the maxima and minima to the plots:

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

Global Optimization

If we instead want to find global maxima and minima, we need to look somewhere else. There are many choices for this. A pure Julia option are the BlackBoxOptim solvers within Optimization.jl, but I will continue the story with the OptimizationNLopt methods. To do this, we simply swap out to one of the global optimizers in the list. Let’s try 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")