# Estimate the parameters of the Lorenz system from the dataset

Note: If data is generated with a fixed time step method and then is tested against with the same time step, there is a biased introduced since it's no longer about hitting the true solution, rather it's just about retreiving the same values that the ODE was first generated by! Thus this version uses adaptive timestepping for all portions so that way tests are against the true solution.

using ParameterizedFunctions, OrdinaryDiffEq, DiffEqParamEstim
using BlackBoxOptim, NLopt, Plots, QuadDIRECT
gr(fmt=:png)

Plots.GRBackend()

Xiang2015Bounds = Tuple{Float64, Float64}[(9, 11), (20, 30), (2, 3)] # for local optimizations
xlow_bounds = [9.0,20.0,2.0]
xhigh_bounds = [11.0,30.0,3.0]
LooserBounds = Tuple{Float64, Float64}[(0, 22), (0, 60), (0, 6)] # for global optimization
GloIniPar = [0.0, 0.5, 0.1] # for global optimizations
LocIniPar = [9.0, 20.0, 2.0] # for local optimization

g1 = @ode_def LorenzExample begin
dx = σ*(y-x)
dy = x*(ρ-z) - y
dz = x*y - β*z
end σ ρ β
p = [10.0,28.0,2.66] # Parameters used to construct the dataset
r0 = [1.0; 0.0; 0.0]                #[-11.8,-5.1,37.5] PODES Initial values of the system in space # [0.1, 0.0, 0.0]
tspan = (0.0, 30.0)                 # PODES sample of 3000 observations over the (0,30) timespan
prob = ODEProblem(g1, r0, tspan,p)
tspan2 = (0.0, 3.0)                 # Xiang test sample of 300 observations with a timestep of 0.01
prob_short = ODEProblem(g1, r0, tspan2,p)

dt = 30.0/3000
tf = 30.0
tinterval = 0:dt:tf
t  = collect(tinterval)

h = 0.01
M = 300
tstart = 0.0
tstop = tstart + M * h
tinterval_short = 0:h:tstop
t_short = collect(tinterval_short)

# Generate Data
data_sol_short = solve(prob_short,Vern9(),saveat=t_short,reltol=1e-9,abstol=1e-9)
data_short = convert(Array, data_sol_short) # This operation produces column major dataset obs as columns, equations as rows
data_sol = solve(prob,Vern9(),saveat=t,reltol=1e-9,abstol=1e-9)
data = convert(Array, data_sol)


Plot the data

plot(data_sol_short,vars=(1,2,3)) # the short solution
plot(data_sol,vars=(1,2,3)) # the longer solution
interpolation_sol = solve(prob,Vern7(),saveat=t,reltol=1e-12,abstol=1e-12)
plot(interpolation_sol,vars=(1,2,3))

xyzt = plot(data_sol_short, plotdensity=10000,lw=1.5)
xy = plot(data_sol_short, plotdensity=10000, vars=(1,2))
xz = plot(data_sol_short, plotdensity=10000, vars=(1,3))
yz = plot(data_sol_short, plotdensity=10000, vars=(2,3))
xyz = plot(data_sol_short, plotdensity=10000, vars=(1,2,3))
plot(plot(xyzt,xyz),plot(xy, xz, yz, layout=(1,3),w=1), layout=(2,1), size=(800,600))

xyzt = plot(data_sol, plotdensity=10000,lw=1.5)
xy = plot(data_sol, plotdensity=10000, vars=(1,2))
xz = plot(data_sol, plotdensity=10000, vars=(1,3))
yz = plot(data_sol, plotdensity=10000, vars=(2,3))
xyz = plot(data_sol, plotdensity=10000, vars=(1,2,3))
plot(plot(xyzt,xyz),plot(xy, xz, yz, layout=(1,3),w=1), layout=(2,1), size=(800,600))


## Find a local solution for the three parameters from a short data set

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short)
res1 = bboptimize(obj_short;SearchRange = LooserBounds, MaxSteps = 7e3)

Starting optimization with optimizer BlackBoxOptim.DiffEvoOpt{BlackBoxOptim
daptiveDiffEvoRandBin{3},BlackBoxOptim.RandomBound{BlackBoxOptim.RangePerDi
mSearchSpace}}
0.00 secs, 0 evals, 0 steps
0.50 secs, 1558 evals, 1457 steps, improv/step: 0.312 (last = 0.3116), fitn
ess=22.641385792
1.00 secs, 3269 evals, 3169 steps, improv/step: 0.289 (last = 0.2693), fitn
ess=0.011380680
1.50 secs, 5114 evals, 5015 steps, improv/step: 0.290 (last = 0.2931), fitn
ess=0.000000695
2.00 secs, 6965 evals, 6867 steps, improv/step: 0.286 (last = 0.2754), fitn
ess=0.000000000

Optimization stopped after 7001 steps and 2.0429041385650635 seconds
Termination reason: Max number of steps (7000) reached
Steps per second = 3426.984099663876
Function evals per second = 3474.955024069969
Improvements/step = 0.28614285714285714
Total function evaluations = 7099

Best candidate found: [10.0, 28.0, 2.66]

Fitness: 0.000000000

# Tolernace is still too high to get close enough

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9)
res1 = bboptimize(obj_short;SearchRange = LooserBounds, MaxSteps = 7e3)

Starting optimization with optimizer BlackBoxOptim.DiffEvoOpt{BlackBoxOptim
daptiveDiffEvoRandBin{3},BlackBoxOptim.RandomBound{BlackBoxOptim.RangePerDi
mSearchSpace}}
0.00 secs, 0 evals, 0 steps
0.50 secs, 1254 evals, 1095 steps, improv/step: 0.283 (last = 0.2831), fitn
ess=80.049861865
1.00 secs, 2503 evals, 2343 steps, improv/step: 0.277 (last = 0.2708), fitn
ess=0.255816863
1.50 secs, 3798 evals, 3638 steps, improv/step: 0.280 (last = 0.2849), fitn
ess=0.000855484
2.00 secs, 5083 evals, 4924 steps, improv/step: 0.286 (last = 0.3025), fitn
ess=0.000001922
2.50 secs, 6363 evals, 6205 steps, improv/step: 0.285 (last = 0.2849), fitn
ess=0.000000006

Optimization stopped after 7001 steps and 2.8151230812072754 seconds
Termination reason: Max number of steps (7000) reached
Steps per second = 2486.925011107364
Function evals per second = 2543.0504434391687
Improvements/step = 0.2847142857142857
Total function evaluations = 7159

Best candidate found: [10.0, 28.0, 2.66]

Fitness: 0.000000000

# With the tolerance lower, it achieves the correct solution in 3.5 seconds.

obj_short = build_loss_objective(prob_short,Vern9(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9,abstol=1e-9)
res1 = bboptimize(obj_short;SearchRange = LooserBounds, MaxSteps = 7e3)

Starting optimization with optimizer BlackBoxOptim.DiffEvoOpt{BlackBoxOptim
daptiveDiffEvoRandBin{3},BlackBoxOptim.RandomBound{BlackBoxOptim.RangePerDi
mSearchSpace}}
0.00 secs, 0 evals, 0 steps
0.50 secs, 1569 evals, 1427 steps, improv/step: 0.315 (last = 0.3146), fitn
ess=5.877435497
1.00 secs, 3169 evals, 3029 steps, improv/step: 0.294 (last = 0.2765), fitn
ess=0.000432088
1.50 secs, 4736 evals, 4596 steps, improv/step: 0.291 (last = 0.2846), fitn
ess=0.000002200
2.00 secs, 6289 evals, 6149 steps, improv/step: 0.294 (last = 0.3007), fitn
ess=0.000000002

Optimization stopped after 7001 steps and 2.2673919200897217 seconds
Termination reason: Max number of steps (7000) reached
Steps per second = 3087.6885191171395
Function evals per second = 3148.9924334375623
Improvements/step = 0.2915714285714286
Total function evaluations = 7140

Best candidate found: [10.0, 28.0, 2.66]

Fitness: 0.000000000

# With the more accurate solver Vern9 in the solution of the ODE, the convergence is less efficient!

# Fastest BlackBoxOptim: 3.5 seconds


# Using NLopt

First, the global optimization algorithms

obj_short = build_loss_objective(prob_short,Vern9(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9,abstol=1e-9)

opt = Opt(:GN_ORIG_DIRECT_L, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Accurate 3.2 seconds

1.258503 seconds (5.24 M allocations: 498.290 MiB, 9.85% gc time)
(7.389526220984814e-18, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:GN_CRS2_LM, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Accurate 3.0 seconds

1.119735 seconds (4.73 M allocations: 450.375 MiB, 9.81% gc time)
(2.7919617960803954e-18, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:GN_ISRES, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Accurate to single precision 8.2 seconds

3.123742 seconds (13.41 M allocations: 1.247 GiB, 9.04% gc time)
(0.0020472369723499527, [9.99773, 28.0006, 2.66001], :MAXEVAL_REACHED)

opt = Opt(:GN_ESCH, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Approximatively accurate, good starting values for local optimization

3.150222 seconds (13.41 M allocations: 1.247 GiB, 9.21% gc time)
(413.36121809892586, [10.5187, 27.6462, 2.80337], :MAXEVAL_REACHED)


Next, the local optimization algorithms that could be used after the global algorithms as a check on the solution and its precision. All the local optimizers are started from LocIniPar and with the narrow bounds of the Xiang2015Paper.

opt = Opt(:LN_BOBYQA, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # 0.1 seconds

0.032699 seconds (148.87 k allocations: 14.175 MiB)
(2.7636309213683456e-18, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:LN_NELDERMEAD, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.29 sec

0.110738 seconds (439.58 k allocations: 41.814 MiB, 8.86% gc time)
(2.767990920967993e-18, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:LD_SLSQP, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.21 sec

0.123957 seconds (448.49 k allocations: 34.165 MiB, 11.93% gc time)
(1.1116793161130144e-15, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:LN_COBYLA, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 1.84 sec

0.726348 seconds (3.10 M allocations: 294.844 MiB, 8.94% gc time)
(2.8465209128971724e-18, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:LN_NEWUOA_BOUND, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.18 sec ROUNDOFF LIMITED

0.096803 seconds (254.81 k allocations: 24.262 MiB, 6.18% gc time)
(2.514234891513546e-9, [10.0, 28.0, 2.66], :SUCCESS)

opt = Opt(:LN_PRAXIS, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.18 sec

0.072654 seconds (278.94 k allocations: 26.561 MiB, 14.07% gc time)
(2.763441276332588e-18, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:LN_SBPLX, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.65 sec

0.237722 seconds (1.00 M allocations: 95.515 MiB, 8.21% gc time)
(2.779471170077341e-18, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:LD_MMA, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.7 sec

0.262558 seconds (1.11 M allocations: 105.400 MiB, 9.71% gc time)
(2.506182282152645e-16, [10.0, 28.0, 2.66], :XTOL_REACHED)

opt = Opt(:LD_LBFGS, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.12 sec

0.040653 seconds (168.59 k allocations: 16.078 MiB, 13.50% gc time)
(1.116061455658539e-15, [10.0, 28.0, 2.66], :SUCCESS)

opt = Opt(:LD_TNEWTON_PRECOND_RESTART, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.15 sec

0.053640 seconds (215.41 k allocations: 20.545 MiB, 9.43% gc time)
(1.1161902238149476e-15, [10.0, 28.0, 2.66], :SUCCESS)


## Now let's solve the longer version for a global solution

Notice from the plotting above that this ODE problem is chaotic and tends to diverge over time. In the longer version of parameter estimation, the dataset is increased to 3000 observations per variable with the same integration time step of 0.01. Vern9 solver with reltol=1e-9 and abstol=1e-9 has been established to be accurate on the time interval [0,50]

# BB with Vern9 converges very slowly. The final values are within the NarrowBounds.
obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),tstops=t,reltol=1e-9,abstol=1e-9)

res1 = bboptimize(obj;SearchRange = LooserBounds, MaxSteps = 4e3) # Default adaptive_de_rand_1_bin_radiuslimited 33 sec [10.2183, 24.6711, 2.28969]

Starting optimization with optimizer BlackBoxOptim.DiffEvoOpt{BlackBoxOptim
daptiveDiffEvoRandBin{3},BlackBoxOptim.RandomBound{BlackBoxOptim.RangePerDi
mSearchSpace}}
0.00 secs, 0 evals, 0 steps
0.50 secs, 149 evals, 92 steps, improv/step: 0.554 (last = 0.5543), fitness
=586639.804024005
1.01 secs, 305 evals, 211 steps, improv/step: 0.374 (last = 0.2353), fitnes
s=526920.067405802
1.51 secs, 465 evals, 350 steps, improv/step: 0.346 (last = 0.3022), fitnes
s=526920.067405802
2.01 secs, 625 evals, 504 steps, improv/step: 0.349 (last = 0.3571), fitnes
s=526920.067405802
2.51 secs, 785 evals, 661 steps, improv/step: 0.324 (last = 0.2420), fitnes
s=526920.067405802
3.01 secs, 942 evals, 818 steps, improv/step: 0.312 (last = 0.2611), fitnes
s=507855.414342085
3.52 secs, 1100 evals, 976 steps, improv/step: 0.288 (last = 0.1646), fitne
ss=507855.414342085
4.02 secs, 1254 evals, 1130 steps, improv/step: 0.273 (last = 0.1753), fitn
ess=507855.414342085
4.52 secs, 1410 evals, 1286 steps, improv/step: 0.259 (last = 0.1603), fitn
ess=507855.414342085
5.03 secs, 1571 evals, 1447 steps, improv/step: 0.254 (last = 0.2174), fitn
ess=507855.414342085
5.53 secs, 1730 evals, 1606 steps, improv/step: 0.251 (last = 0.2201), fitn
ess=507855.414342085
6.03 secs, 1889 evals, 1765 steps, improv/step: 0.241 (last = 0.1447), fitn
ess=507855.414342085
6.53 secs, 2048 evals, 1924 steps, improv/step: 0.234 (last = 0.1572), fitn
ess=507855.414342085
7.03 secs, 2207 evals, 2083 steps, improv/step: 0.227 (last = 0.1321), fitn
ess=470702.977254070
7.53 secs, 2365 evals, 2241 steps, improv/step: 0.221 (last = 0.1456), fitn
ess=470702.977254070
8.04 secs, 2523 evals, 2399 steps, improv/step: 0.217 (last = 0.1646), fitn
ess=470702.977254070
8.54 secs, 2678 evals, 2554 steps, improv/step: 0.212 (last = 0.1355), fitn
ess=470702.977254070
9.04 secs, 2837 evals, 2713 steps, improv/step: 0.204 (last = 0.0755), fitn
ess=470702.977254070
9.54 secs, 2993 evals, 2869 steps, improv/step: 0.196 (last = 0.0577), fitn
ess=459341.294198028
10.04 secs, 3142 evals, 3018 steps, improv/step: 0.189 (last = 0.0537), fit
ness=459341.294198028
10.55 secs, 3299 evals, 3175 steps, improv/step: 0.183 (last = 0.0637), fit
ness=459341.294198028
11.05 secs, 3458 evals, 3334 steps, improv/step: 0.177 (last = 0.0629), fit
ness=459341.294198028
11.55 secs, 3619 evals, 3495 steps, improv/step: 0.171 (last = 0.0311), fit
ness=459341.294198028
12.05 secs, 3779 evals, 3655 steps, improv/step: 0.164 (last = 0.0250), fit
ness=459341.294198028
12.55 secs, 3940 evals, 3816 steps, improv/step: 0.159 (last = 0.0373), fit
ness=459341.294198028
13.05 secs, 4100 evals, 3976 steps, improv/step: 0.155 (last = 0.0625), fit
ness=459341.294198028

Optimization stopped after 4001 steps and 13.135531902313232 seconds
Termination reason: Max number of steps (4000) reached
Steps per second = 304.5936799327787
Function evals per second = 314.0337239996781
Improvements/step = 0.15425
Total function evaluations = 4125

Best candidate found: [10.9959, 24.9153, 2.42511]

Fitness: 459341.294198028

#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :adaptive_de_rand_1_bin, MaxSteps = 4e3) # Method 32 sec [13.2222, 25.8589, 2.56176]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :dxnes, MaxSteps = 2e3) # Method dxnes 119 sec  [16.8648, 24.393, 2.29119]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :xnes, MaxSteps = 2e3) # Method xnes 304 sec  [19.1647, 24.9479, 2.39467]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :de_rand_1_bin_radiuslimited, MaxSteps = 2e3) # Method 44 sec  [13.805, 24.6054, 2.37274]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :generating_set_search, MaxSteps = 2e3) # Method 195 sec [19.1847, 24.9492, 2.39412]

# using Evolutionary
# N = 3
# @time result, fitness, cnt = cmaes(obj, N; μ = 3, λ = 12, iterations = 1000) # cmaes( rastrigin, N; μ = 15, λ = P, tol = 1e-8)

opt = Opt(:GN_ORIG_DIRECT_L, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Fail to converge

6.717586 seconds (25.97 M allocations: 2.235 GiB, 7.93% gc time)
(470298.735717923, [7.04666, 23.6661, 1.8066], :XTOL_REACHED)

opt = Opt(:GN_CRS2_LM, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 20000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Hit and miss. converge approximately accurate values for local opt.91 seconds

63.291080 seconds (243.03 M allocations: 20.916 GiB, 8.02% gc time)
(526453.8669855164, [22.0, 24.9404, 2.33334], :MAXEVAL_REACHED)

opt = Opt(:GN_ISRES, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 50000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Approximately accurate within local bounds

155.589962 seconds (607.55 M allocations: 52.287 GiB, 7.96% gc time)
(504539.26387804036, [15.466, 25.0421, 2.35805], :MAXEVAL_REACHED)

opt = Opt(:GN_ESCH, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 20000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Approximately accurate

62.609640 seconds (243.02 M allocations: 20.915 GiB, 7.90% gc time)
(496453.9988525161, [11.3667, 23.5838, 2.21621], :MAXEVAL_REACHED)


This parameter estimation on the longer sample proves to be extremely challenging for the global optimizers. BlackBoxOptim is best in optimizing the objective function. All of the global algorithms produces final parameter estimates that could be used as starting values for further refinement with the local optimization algorithms.

opt = Opt(:LN_BOBYQA, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Claims SUCCESS but does not iterate to the true values.

0.352244 seconds (1.36 M allocations: 119.935 MiB, 7.57% gc time)
(588113.2784194095, [9.86259, 20.5811, 2.0], :SUCCESS)

opt = Opt(:LN_NELDERMEAD, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-9)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Inaccurate final values

31.582817 seconds (121.51 M allocations: 10.457 GiB, 7.84% gc time)
(404754.5095397624, [9.67892, 23.5168, 2.16107], :MAXEVAL_REACHED)

opt = Opt(:LD_SLSQP, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Inaccurate final values

0.076167 seconds (267.27 k allocations: 23.557 MiB, 6.90% gc time)
(575377.5788505444, [9.6232, 23.116, 2.3116], :FAILURE)


No local optimizer can improve the global solution to the true values.

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short)
lower = [0.0,0.0,0.0]
upper = [50.0,50.0,50.0]
splits = ([1.0,5.0,15.0],[0,10,20],[0,10,20])
@time root, x0 = analyze(obj_short,splits,lower,upper)

10.791007 seconds (30.06 M allocations: 1.593 GiB, 9.30% gc time)
(BoxRoot@[NaN, NaN, NaN], [5.0, 10.0, 10.0])

minimum(root)

Box1.0816479008038839e-6@[9.99994, 28.0, 2.66]

obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),tstops=t,reltol=1e-9,abstol=1e-9)
lower = [0.0,0.0,0.0]
upper = [50.0,50.0,50.0]
splits = ([0,5.0,15.0],[0,15,30],[0,2,5])
@time root, x0 = analyze(obj,splits,lower,upper)

3.087581 seconds (9.39 M allocations: 735.529 MiB, 7.90% gc time)
(BoxRoot@[NaN, NaN, NaN], [5.0, 15.0, 2.0])

minimum(root)

Box528522.0382483905@[26.702, 24.7704, 2.24937]


# Conclusion:

1. As expected the Lorenz system is extremely sensitive to initial space values. Starting the integration from r0 = [0.1,0.0,0.0] produces convergence with the short sample of 300 observations. This can be achieved by all the global optimizers as well as most of the local optimizers. Instead starting from r0= [-11.8,-5.1,37.5], as in PODES, with the shorter sample shrinks the number of successful algorithms to 3: BBO, :GN_CRS2_LMand :LD_SLSQP. For the longer sample, all the algorithms fail.

2. When trying to hit the real data, having a low enough tolerance on the numerical solution is key. If the numerical solution is too rough, then we can never actually hone in on the true parameters since even with the true parameters we will erroneously induce numerical error. Maybe this could be adaptive?

3. Excessively low tolerance in the numerical solution is inefficient and delays the convergence of the estimation.

4. The estimation method and the global versus local optimization make a huge difference in the timings. Here, BBO always find the correct solution for a global optimization setup. For local optimization, most methods in NLopt, like :LN_BOBYQA, solve the problem in <0.05 seconds. This is an algorithm that can scale a local optimization but we are aiming to scale a global optimization.

5. QuadDIRECT performs very well on the shorter problem but doesn't give very great results for the longer in the Lorenz case, more can be read about the algorithm here.

6. Fitting shorter timespans is easier... maybe this can lead to determining a minimal sample size for the optimizers and the estimator to succeed.