1

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(" ")
4
  • 2
    Hi @AndyDufresne, providing a working example would make the helping you easier. Commented Dec 3, 2024 at 12:36
  • connected to discourse.julialang.org/t/… Commented Dec 3, 2024 at 12:36
  • 1
    Hi @Tanj, thanks. I provided an example. Commented Dec 4, 2024 at 1:02
  • Could you also provide the package versions? julia> using Pkg; Pkg.status() Commented Dec 6, 2024 at 16:16

1 Answer 1

2

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} ...
Sign up to request clarification or add additional context in comments.

1 Comment

Thanks so much! Such a detailed, helpful answer.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.