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 OrdinaryDiffEq
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]
     = u[2]
    β  = u[3]
     = u[4]
    udot .=
    [2(-(1+cos(β)))/(3-cos(2β)),
    -2sin(α) - sin(α+β),
    2(-(1+cos(β)) + (3+2cos(β)))/(3-cos(2β)),
    -sin(α+β) - 2sin(β)*(((-))/(3-cos(2β))) + 2sin(2β)*((^2 - 2(1+cos(β))* + (3+2cos(β))^2)/(3-cos(2β))^2)]
end

#Pass to solvers
poincare = ODEProblem(double_pendulum_hamiltonian, initial, tspan)
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())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 193-element Array{Float64,1}:
   0.0                
   0.08332584852065579
   0.24175280271811872
   0.4389536500504315 
   0.6797322542488147 
   0.964763376337819  
   1.3179449556841032 
   1.7031210236280163 
   2.0678477932001846 
   2.471782525434673  
   ⋮                  
  95.84571836675003   
  96.35777612654726   
  96.9291238553263    
  97.4467872981331    
  97.9624744296349    
  98.51182496995675   
  99.06081878698582   
  99.58283477685029   
 100.0                
u: 193-element Array{Array{Float64,1},1}:
 [0.01, 0.01, 0.01, 0.01]                          
 [0.00917069, 0.006669, 0.0124205, 0.00826641]     
 [0.00767328, 0.000374625, 0.0164426, 0.00463683]  
 [0.00612597, -0.00730546, 0.0199674, -0.000336506]
 [0.0049661, -0.0163086, 0.0214407, -0.00670509]   
 [0.00479557, -0.0262381, 0.0188243, -0.0139134]   
 [0.00605469, -0.0371246, 0.0100556, -0.0210382]   
 [0.00790078, -0.046676, -0.00267353, -0.025183]   
 [0.00827652, -0.0527843, -0.0127315, -0.0252581]  
 [0.00552358, -0.0552525, -0.0168439, -0.021899]   
 ⋮                                                 
 [-0.0148868, 0.0423324, 0.0136282, 0.0180291]     
 [-0.00819054, 0.0544225, 0.00944831, 0.0177401]   
 [0.00412448, 0.0567489, -0.00515392, 0.017597]    
 [0.0130796, 0.0480772, -0.0137706, 0.0182866]     
 [0.0153161, 0.0316313, -0.00895722, 0.0171185]    
 [0.0111156, 0.00992938, 0.0072972, 0.0103535]     
 [0.00571392, -0.0117872, 0.020508, -0.00231029]   
 [0.00421143, -0.0299109, 0.0187506, -0.0156505]   
 [0.00574124, -0.0416539, 0.00741327, -0.023349]

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

using Optim
opt = optimize(f,18.0,22.0)
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)
18.632126799595834
println(opt.minimum)
-0.027931635264246277

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)
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")

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())
Results of Optimization Algorithm
 * Algorithm: BFGS
 * Starting Point: [20.0]
 * Minimizer: [23.297607288716723]
 * Minimum: -2.588588e-02
 * Iterations: 4
 * Convergence: true
   * |x - x'| ≤ 0.0e+00: false 
     |x - x'| = 1.11e-04 
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: false
     |f(x) - f(x')| = -6.49e-09 |f(x)|
   * |g(x)| ≤ 1.0e-08: true 
     |g(x)| = 8.41e-12 
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 16
 * Gradient 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:

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)
-0.027931635264246215 [18.6321] XTOL_REACHED
NLopt.max_objective!(opt, g)
(maxf,maxx,ret) = NLopt.optimize(opt,[20.0])
println(maxf," ",maxx," ",ret)
0.027968571933041936 [6.5537] XTOL_REACHED
plot(sol, vars=(0,4), plotdensity=10000)
scatter!([minx],[minf],label="Global Min")
scatter!([maxx],[maxf],label="Global Max")

Appendix

This tutorial is part of the DiffEqTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

using DiffEqTutorials
DiffEqTutorials.weave_file("ode_extras","03-ode_minmax.jmd")

Computer Information:

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, ivybridge)

Package Information:

Status `~/.julia/environments/v1.1/Project.toml`
[7e558dbc-694d-5a72-987c-6f4ebed21442] ArbNumerics 0.5.4
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.4.2
[be33ccc6-a3ff-5ff2-a52e-74243cff1e17] CUDAnative 2.2.0
[3a865a2d-5b23-5a0f-bc46-62713ec82fae] CuArrays 1.0.2
[55939f99-70c6-5e9b-8bb0-5071ed7d61fd] DecFP 0.4.8
[abce61dc-4473-55a0-ba07-351d65e31d42] Decimals 0.4.0
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 1.1.0
[eb300fae-53e8-50a0-950c-e21f52c2b7e0] DiffEqBiological 3.8.2
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.9.0
[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.6.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.1.0
[6d1b261a-3be8-11e9-3f2f-0b112a9a8436] DiffEqTutorials 0.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.4.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.20.0
[497a8b3b-efae-58df-a0af-a86822472b78] DoubleFloats 0.9.1
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.3
[c91e804a-d5a3-530f-b6f0-dfbca275c004] Gadfly 1.0.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.18.1
[4138dd39-2aa7-5051-a626-17a0bb65d9c8] JLD 0.9.1
[23fbe1c1-3f47-55db-b15f-69d7ec21a316] Latexify 0.8.2
[eff96d63-e80a-5855-80a2-b1b0885c5ab7] Measurements 2.0.0
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 0.2.0
[76087f3c-5699-56af-9a33-bf431cd00edd] NLopt 0.5.1
[2774e3e8-f4cf-5e23-947b-6d7e65073b56] NLsolve 4.0.0
[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.18.1
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.8.1
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.1.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.25.1
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 0.20.0
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.11.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.11.0
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.6.1
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 0.15.0
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.0
[b77e0a4c-d291-57a0-90e8-8db25a27a240] InteractiveUtils
[37e2e46d-f89d-539d-b4ee-838fcccc9c8e] LinearAlgebra
[44cfe95a-1eb2-52ea-b672-e2afdf69b78f] Pkg