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