Как исправить ошибку TypeError: в setindex!в DifferentialEquations.jl - PullRequest
Новая
       118

Как исправить ошибку TypeError: в setindex!в DifferentialEquations.jl

0 голосов
/ 29 января 2019

Недавно я начал с пакета Julia (v1.0.3) DifferentialEquations.jl.Я попытался решить простую систему ODE, с той же структурой, что и моя реальная модель, но гораздо меньшеВ зависимости от решателя, который я использую, пример либо решает, либо выдает ошибку.Рассмотрим MWE, модель химической инженерии последовательной / параллельной реакции в 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 не будет решен, и я получу следующую ошибку: 

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

Я не буду вставлять сюда всю трассировку стека, так какэто очень долго, но вы можете легко воспроизвести это, изменив sol = solve(prob) на sol = solve(prob, Rosenbrock23()).Мне кажется, что ошибка возникает, когда решатель пытается вывести якобиан, но я понятия не имею, почему.И почему решатель по умолчанию работает, а другие нет?

Может, кто-нибудь подскажет, почему возникает эта ошибка и как ее можно исправить?

1 Ответ

0 голосов
/ 29 января 2019

Автоматическое дифференцирование работает, пропуская Dual типы через вашу функцию, а не с плавающей точкой, с которой вы обычно ее используете.Таким образом, проблема возникает из-за того, что вы фиксируете внутреннее значение rate для типа Vector{Float64} (см. Третью точку здесь и этот совет ).К счастью, это легко исправить (и даже лучше выглядит, ИМХО): 

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.

В качестве альтернативы, вы можете отключить AD с помощьюRosenbrock23(autodiff=false) (который, я думаю, вместо этого будет использовать конечные различия) или снабдит якобиана.

...