div.ProseMirror

Finding Maxima and Minima of DiffEq Solutions

Chris Rackauckas

Setup

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

#Constants and setup
using DifferentialEquations
initial = [0.01, 0.01, 0.01, 0.01]
tspan = (0.,100.)
#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)
0.5s
Julia
→
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true timespan: (0.0, 100.0) u0: [0.01, 0.01, 0.01, 0.01]
sol = solve(poincare, Tsit5())
2.0s
Julia
→

In time, the solution looks like:

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

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

plot(sol, vars=(3,4), leg=false)
0.4s
Julia
→

Local Optimization

Let's fine 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(t,idxs=4)
0.3s
Julia
→
#3

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 minima, we can simply call Optim on this function. Let's find a local minimum:

using Pkg
pkg"add Optim"
using Optim
opt = optimize(f,18.0,22.0)
34.7s
Julia
→
Results of Optimization Algorithm * Algorithm: Brent's Method * Search Interval: [18.000000, 22.000000] * Minimizer: 1.863213e+01 * Minimum: -2.793164e-02 * Iterations: 11 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true * Objective Function Calls: 12

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.minimizer)
println(opt.minimum)
0.6s
Julia
→

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

f = (t) -> -sol(first(t),idxs=4)
opt2 = optimize(f,0.0,22.0)
0.3s
Julia
→
Results of Optimization Algorithm * Algorithm: Brent's Method * Search Interval: [0.000000, 22.000000] * Minimizer: 1.399975e+01 * Minimum: -2.269411e-02 * Iterations: 13 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true * Objective Function Calls: 14

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

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

Brent's method will locally minimize over the full interval. If we instead want a local maxima nearest to a point, we can use BFGS(). In this case, we need to optimize a vector [t], and thus dereference it to a number using first(t).

f = (t) -> -sol(first(t),idxs=4)
opt = optimize(f,[20.0],BFGS())
3.4s
Julia
→
* Status: success * Candidate solution Final objective value: -2.588588e-02 * Found with Algorithm: BFGS * Convergence measures |x - x'| = 1.11e-04 ≰ 0.0e+00 |x - x'|/|x'| = 4.78e-06 ≰ 0.0e+00 |f(x) - f(x')| = 1.68e-10 ≰ 0.0e+00 |f(x) - f(x')|/|f(x')| = 6.49e-09 ≰ 0.0e+00 |g(x)| = 8.41e-12 ≤ 1.0e-08 * Work counters Seconds run: 0 (vs limit Inf) Iterations: 4 f(x) calls: 16 ∇f(x) calls: 16

Global Optimization

If we instead want to find global maxima and minima, we need to look somewhere else. For this there are many choices. A pure Julia option is BlackBoxOptim.jl, but I will use NLopt.jl. Following the NLopt.jl tutorial but replacing their function with out own:

using Pkg
pkg"add NLopt ForwardDiff"
1.9s
Julia
→
import NLopt, ForwardDiff
count = 0 # keep track of # function evaluations
function g(t::Vector, grad::Vector)
  if length(grad) > 0
    #use ForwardDiff for the gradients
    grad[1] = ForwardDiff.derivative((t)->sol(first(t),idxs=4),t)
  end
  sol(first(t),idxs=4)
end
opt = NLopt.Opt(:GN_ORIG_DIRECT_L, 1)
NLopt.lower_bounds!(opt, [0.0])
NLopt.upper_bounds!(opt, [40.0])
NLopt.xtol_rel!(opt,1e-8)
NLopt.min_objective!(opt, g)
(minf,minx,ret) = NLopt.optimize(opt,[20.0])
println(minf," ",minx," ",ret)
NLopt.max_objective!(opt, g)
(maxf,maxx,ret) = NLopt.optimize(opt,[20.0])
println(maxf," ",maxx," ",ret)
4.7s
Julia
→
plot(sol, vars=(0,4), plotdensity=10000)
scatter!([minx],[minf],label="Global Min")
scatter!([maxx],[maxf],label="Global Max")
1.6s
Julia
→
Runtimes (1)