I am trying to find multiple parameters for the optimized system which is governed by a system of equations.
In my code, I need to be able to use 4 or 5 equations and determine the parameters for the optimized system. I have simplified the code to a system of 3 equations for simplicity, (but it is 5 equations), as below:
function f1(du, u, p, t)
x, y = u
α, β, δ, γ = p
du[1] = dx = α * x - β * x * y
du[2] = dy = -δ * y + γ *x * y
du[3] = dz = α*x
end
u0 = [1.0; 1.0; 1.0]
tspan = (0.0, 10.0)
p = [1.5, 1.0, 3.0, 1.0]
prob = ODEProblem(f1, u0, tspan, p)
t = collect(range(0, stop = 10, length = 200)) # Creates the time vector
data = Array(solve(prob, Tsit5(), saveat = t, abstol = 1e-12, reltol = 1e-12))
bound = Tuple{Float64, Float64}[(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10), (0, 10), (0, 10)]
obj = multiple_shooting_objective(prob, Tsit5(), L2Loss(t, data),Optimization.AutoForwardDiff(); discontinuity_weight = 1.0, abstol = 1e-12, reltol = 1e-12)
optprob = Optimization.OptimizationProblem(obj, zeros(18), lb = first.(bound), ub = last.(bound))
optsol = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 10000)
Optimised_parameters = optsol.u[(end - 1):end]
println("Optimised parameters: ")
println("p₁: ", Optimised_parameters[1])
println("p₂: ", Optimised_parameters[2])
println(" ")
However, I am getting the error message:
ERROR: InexactError: Int64(4.666666666666667)
for the "OptSol" line. Which indicates that somewhere it is receiving an integer but expecting a Float. I am not sure how to fix this so I can extend the system to include more equations. Any help would be much appreciated!
To help, here is a working example (there is one less du term. I am trying to allow multiple more du terms):
using DifferentialEquations, RecursiveArrayTools, Plots, DiffEqParamEstim
using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationBBO
function f1(du, u, p, t)
x, y = u
α, β, δ, γ = p
du[1] = dx = α * x - β * x * y
du[2] = dy = -δ * y + γ *x * y
end
u0 = [1.0; 1.0]
tspan = (0.0, 10.0)
p = [1.5, 1.0, 3.0, 1.0]
prob = ODEProblem(f1, u0, tspan, p)
t = collect(range(0, stop = 10, length = 200)) # Creates the time vector
data = Array(solve(prob, Tsit5(), saveat = t, abstol = 1e-12, reltol = 1e-12))
bound = Tuple{Float64, Float64}[(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10), (0, 10), (0, 10)]
obj = multiple_shooting_objective(prob, Tsit5(), L2Loss(t, data),Optimization.AutoForwardDiff(); discontinuity_weight = 1.0, abstol = 1e-12, reltol = 1e-12)
optprob = Optimization.OptimizationProblem(obj, zeros(18), lb = first.(bound), ub = last.(bound))
optsol = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 10000)
Optimised_parameters = optsol.u[(end - 1):end]
println("Optimised parameters: ")
println("p₁: ", Optimised_parameters[1])
println("p₂: ", Optimised_parameters[2])
println(" ")
Reproducible. Below is some information from a debug session that might help someone explain how to correct the problem.
On Julia 1.10.7 (for Debugger.jl):
using Pkg
Pkg.add("Optimization"); import Optimization: Optimization, ODEProblem, solve
Pkg.add("OptimizationBBO"); import OptimizationBBO: BBO_adaptive_de_rand_1_bin_radiuslimited
Pkg.add("OrdinaryDiffEq"); import OrdinaryDiffEq: Tsit5
Pkg.add("DiffEqParamEstim"); import DiffEqParamEstim: L2Loss, multiple_shooting_objective
function f1(du, u, p, t)
x, y = u
α, β, δ, γ = p
du[1] = dx = α * x - β * x * y
du[2] = dy = -δ * y + γ *x * y
du[3] = dz = α*x
end
u0 = [1.0; 1.0; 1.0]
tspan = (0.0, 10.0)
p = [1.5, 1.0, 3.0, 1.0]
prob = ODEProblem(f1, u0, tspan, p)
t = collect(range(0, stop = 10, length = 200)) # Creates the time vector
data = Array(solve(prob, Tsit5(), saveat = t, abstol = 1e-12, reltol = 1e-12))
bound = Tuple{Float64, Float64}[(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10),
(0, 10), (0, 10), (0, 10), (0, 10), (0, 10), (0, 10)]
obj = multiple_shooting_objective(prob, Tsit5(), L2Loss(t, data), Optimization.AutoForwardDiff(); discontinuity_weight = 1.0, abstol = 1e-12, reltol = 1e-12)
optprob = Optimization.OptimizationProblem(obj, zeros(18), lb = first.(bound), ub = last.(bound))
Pkg.add("Debugger"); using Debugger; break_on(:error)
@run optsol = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 10000)
Produces
Breaking for error:
ERROR: InexactError: Int64(4.666666666666667)
Stacktrace:
[1] Int64(x::Float64)
@ Base float.jl:912
[2] (::DiffEqParamEstim.var"#43#48"{Nothing, Float64, DiffEqParamEstim.var"#1#2", @Kwargs{abstol::Float64, reltol::Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, typeof(f1), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, L2Loss{Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing}, Nothing})(p::Vector{Float64}, ::SciMLBase.NullParameters)
...
In Int64(x) at float.jl:908
908 function (::Type{$Ti})(x::$Tf)
909 if ($(Tf(typemin(Ti))) <= x < $(Tf(typemax(Ti)))) && isinteger(x)
910 return unsafe_trunc($Ti,x)
911 else
>912 throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
913 end
914 end
915 end
916 end
About to run: throw(InexactError(:Int64, Int64, 4.666666666666667))
At the debugger command prompt, moving up to the next frame:
1|debug> up
In #43(p, #unused#) at ~/.julia/packages/DiffEqParamEstim/6z5Ou/src/multiple_shooting_objective.jl:30
29 kwargs...)
30 cost_function = function (p, _ = nothing)
31 t0, tf = prob.tspan
32 P, N = length(prob.p), length(prob.u0)
>33 K = Int((length(p) - P) / N)
34 τ = range(t0, tf, length = K + 1)
35 sol = []
36 loss_val = 0
37 for k in 1:K
About to run: (Int64)(4.666666666666667)
Stack frame:
2|debug> fr
[2] #43(p, #unused#) at ~/.julia/packages/DiffEqParamEstim/6z5Ou/src/multiple_shooting_objective.jl:30
| p::Vector{Float64} = <[3.843240561664537, 8.092526443524417, 0.6805655220003193, 4.4963371509026695, 7.785492340246389, 3.6...>
| ::Int64 = 2
| N::Int64 = 3
| P::Int64 = 4
...
Evaluate a watch expression:
2|debug> w add length(p)
1] length(p): 18
2|debug> q
As N=3, P=4, length(p)=18, thus (length(p)-P)/N = 4.666666666666667
This suggests that (length(p) - P) must be a multiple of N. (A previous version used a floor here.)
Since the code above adds a variable u[3] compared to the working version (and the documentation example), N increases. This suggests 18 must be increased.
Where does 18 come from? The only explanation I found is "Why is it long 18 (tuples)?" "It’s the parameter bound. The first two are the real parameters, the rest are due to the inflection points. Then you have to ignore those in the output." (From here, where the comment notes DiffEqFlux.jl has an alternative.)
After increasing 18 to 19 it runs.
optprob = Optimization.OptimizationProblem(obj, zeros(19), lb = zeros(19), ub = 10 .* ones(19))
optsol = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 10000)
# retcode: MaxIters
# u: 19-element Vector{Float64} ...