# Finding Maxima and Minima of DiffEq Solutions

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


### 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


### 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

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
[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
[44cfe95a-1eb2-52ea-b672-e2afdf69b78f] Pkg