Julia 优化脚本期间出现“Int64”错误消息

问题描述 投票：0回答：1

我正在尝试为由方程组控制的优化系统找到多个参数。

在我的代码中，我需要能够使用 4 或 5 个方程并确定优化系统的参数。为了简单起见，我将代码简化为 3 个方程组（但实际上是 5 个方程），如下所示：


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

但是，我收到错误消息：

ERROR: InexactError: Int64(4.666666666666667)

“OptSol”系列。 这表明某个地方正在接收一个整数，但需要一个浮点数。我不知道如何解决这个问题，以便我可以扩展系统以包含更多方程。 任何帮助将不胜感激！

为了提供帮助，这里有一个工作示例（少了一个 du 术语。我试图允许更多的 du 术语）：

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(" ")
optimization julia ode
1个回答
0
投票

可重复。 以下是调试会话中的一些信息，可能有助于某人解释如何纠正问题。

在 Julia 1.10.7 上（对于 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)

产品

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

在调试器命令提示符处，将 

up
移动到下一帧：

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)

堆栈

fr
名称：

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

评估 

w
atch 表达式：

2|debug> w add length(p)
1] length(p): 18

2|debug> q

当N=3，P=4，长度(p)=18，因此(长度(p)-P)/N = 4.666666666666667

这表明 (length(p) - P) 必须是 N 的倍数。（以前的版本在这里使用了 

floor
。）

由于与工作版本（以及文档示例）相比，上面的代码添加了变量 u[3]，因此 N 增加。 这表明必须增加 18。

18从哪里来？ 我找到的唯一解释是“为什么长18（元组）？” “这是参数绑定。前两个是真正的参数，其余的是拐点。然后你必须忽略输出中的那些。” （从这里，评论指出 DiffEqFlux.jl 有一个替代方案。）

将18增加到19后运行。

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} ...
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.