如何修复TypeError:在setindex中！在DifferentialEquations.jl中 [英] How to fix TypeError: in setindex! in DifferentialEquations.jl

问题描述

最近，我开始使用Julia的(v1.0.3)DifferentialEquations.jl软件包.我尝试解决一个简单的ODE系统，该系统的结构与我的真实模型相同，但是要小得多. 根据我使用的求解器，该示例可以解决或引发错误.考虑这个MWE，它是CSTR中连续/平行反应的化学工程模型:

Recently, I got started with Julia's (v1.0.3) DifferentialEquations.jl package. I tried solving a simple ODE system, with the same structure as my real model, but much smaller. Depending on the solver which I use, the example either solves or throws an error. Consider this MWE, a Chemical Engineering model of a consecutive / parallel reaction in a CSTR:

using DifferentialEquations
using Plots

# Modeling a consecutive / parallel reaction in a CSTR
# A --> 2B --> C, C --> 2B, B --> D
# PETERSEN-Matrix
#   No.     A       B       C       D       Rate
#   1      -1       2                       k1*A
#   2              -2       1               k2*B*B
#   3               2      -1               k3*C
#   4              -1               1       k4*B

function fpr(dx, x, params, t)
    k_1, k_2, k_3, k_4, q_in, V_liq, A_in, B_in, C_in, D_in = params
    # Rate equations
    rate = Array{Float64}(undef, 4)
    rate[1] = k_1*x[1]
    rate[2] = k_2*x[2]*x[2]
    rate[3] = k_3*x[3]
    rate[4] = k_4*x[2]

    dx[1] = -rate[1] + q_in/V_liq*(A_in - x[1])
    dx[2] = 2*rate[1] - 2*rate[2] + 2*rate[3] - rate[4] + q_in/V_liq*(B_in - x[2])
    dx[3] = rate[2] - rate[3] + q_in/V_liq*(C_in - x[3])
    dx[4] = rate[4] + q_in/V_liq*(D_in - x[4])
end 

u0 = [1.5, 0.1, 0, 0]
params = [1.0, 1.5, 0.75, 0.15, 3, 15, 0.5, 0, 0, 0]
tspan = (0.0, 15.0)
prob = ODEProblem(fpr, u0, tspan, params)
sol = solve(prob)
plot(sol)

这很好用. 但是，如果选择其他求解器，例如Rosenbrock23()或Rodas4()，则ODE无法解决，并且出现以下错误:

This works perfectly. However, if a choose a different solver, say Rosenbrock23() or Rodas4(), the ODE is not solved and I get the following error:

ERROR: LoadError: TypeError: in setindex!, in typeassert, expected Float64,
got ForwardDiff.Dual{Nothing,Float64,4}

我不会在此处粘贴整个堆栈跟踪，因为它很长，但是您可以通过将sol = solve(prob)更改为sol = solve(prob, Rosenbrock23())来轻松地重现该堆栈.在我看来，当求解程序尝试导出Jacobian时会发生错误，但是我不知道为什么.以及为什么默认求解器可以工作，而其他求解器却不能工作?

I won't paste the whole stacktrace here, since it is very long, but you can easily reproduce this by changing sol = solve(prob) into sol = solve(prob, Rosenbrock23()). It seems to me that the error occurs when the solver tries to derive Jacobians, but I have no clue why. And why does the default solver work, but others don't?

请，有人可以告诉我为什么会发生此错误，以及如何解决该错误吗?

Please, could anyone tell me why this error occurs and how it can be fixed?

推荐答案

自动区分的工作原理是通过函数传递Dual类型，而不是通常使用浮点数的函数.因此出现问题是因为您将内部值rate固定为类型Vector{Float64}(请参见第三点

Automatic differentiation works by passing Dual types through your function, instead of the floats you would normally use it with. So the problem arises because you fix the internal value rate to be of type Vector{Float64} (see the third point here, and this advice). Fortunately, that's easy to fix (and even better looking, IMHO):

julia> function fpr(dx, x, params, t)
           k_1, k_2, k_3, k_4, q_in, V_liq, A_in, B_in, C_in, D_in = params
           # Rate equations
           # should actually be rate = [k_1*x[1], k_2*x[2]*x[2], k_3*x[3], k_4*x[2]], as per @LutzL's comment
           rate = [k_1*x[1], k_2*x[2], k_3*x[3], k_4*x[2]]

           dx[1] = -rate[1] + q_in/V_liq*(A_in - x[1])
           dx[2] = 2*rate[1] - 2*rate[2] + 2*rate[3] - rate[4] + q_in/V_liq*(B_in - x[2])
           dx[3] = rate[2] - rate[3] + q_in/V_liq*(C_in - x[3])
           dx[4] = rate[4] + q_in/V_liq*(D_in - x[4])
       end

对Rosenbrock23和Rodas4均适用.

或者，您可以使用Rosenbrock23(autodiff=false)(我认为它将使用有限差分代替)关闭AD，或者提供雅可比行列式.

Alternatively, you can turn off AD with Rosenbrock23(autodiff=false) (which, I think, will use finite differences instead), or supply a Jacobian.

这篇关于如何修复TypeError:在setindex中！在DifferentialEquations.jl中的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！