Как прочитать соответствующее значение в указанное время в DiscreteCallback?

Подобно этому вопросу, я пытаюсь решить это ODE с зависящим от времени входным параметром. Он состоит из серии дискретные обратные вызовы. В определенные моменты изменяется параметр (не состояние!). Время и значения хранятся в nx2 Array. Но я не могу заставить функцию affect найти соответствующее значение параметра в указанное время. В приведенных примерах значение, присвоенное u[1], обычно является постоянным. Рассмотрим этот MWE (с подходом, очень похожим на Matlab), который работает правильно без обратного вызова:

using DifferentialEquations
using Plots

function odm2prod(dx, x, params, t)
    k_1, f_1, V_liq, X_in, Y_in, q_in = params

    rho_1 = k_1*x[1]
    q_prod = 0.52*f_1*x[1]
    # Differential Equations
    dx[1] = q_in/V_liq*(X_in - x[1]) - rho_1
    dx[2] = q_in/V_liq*(Y_in - x[2])
end

x0      = [3.15, 1.5]
tspan   = (0.0, 7.0)
params  = [0.22, 43, 155, 249, 58, 0]
prob    = ODEProblem(odm2prod, x0, tspan, params)

input   = [1.0 60; 1.1 0; 2.0 60; 2.3 0; 4.0 430; 4.05 0]
dosetimes = input[:,1]
function affect!(integrator)
    ind_t = findall(integrator.t == dosetimes)
    integrator.p[6] = input[ind_t, 2]
end
cb = PresetTimeCallback(dosetimes, affect!)
sol = solve(prob, Tsit5(), callback=cb, saveat=1/12)

plot(sol, vars=[1, 2])

Это не работает. Ошибка возникает в строке 22, поскольку сравнение вектора со скаляром, похоже, не определено в Julia, или существует особый синтаксис, о котором я не знаю.

Я знаю, что можно использовать зависящие от времени параметры в Julia, но я полагаю, что это будет работать только для непрерывных функций, а не для дискретных изменений !? Я не просматривал справку для interpolate, но не знаю, как использовать ее в моем конкретном случае.

Может кто-нибудь сказать мне, как заставить это работать, пожалуйста? Вероятно, потребуется всего несколько строк кода. Кроме того, я не обязательно хочу, чтобы dosetimes был частью sol.t, если они не совпадают.


julia differential-equations differentialequations.jl
person rotton    schedule 25.01.2020    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Вы используете findall неправильно, говорится в документации.

findall(f::Function, A)

Возвращает вектор I индексов или ключей A, где f(A[I]) возвращает true.

Тогда вы должны принять во внимание, что результатом поиска по всем является список. Поскольку вы ожидаете, что у него будет только один элемент, используйте только первый

function affect!(integrator)
    ind_t = findall(t -> t==integrator.t, dosetimes)
    integrator.p[6] = input[ind_t[1], 2]
end

и вы получите сюжет

сюжет решения оды

person Lutz Lehmann    schedule 25.01.2020
comment
На самом деле, я хотел имитировать функцию Matlab ismember с findall Джулии. Но тогда я должен был поставить аргументы наоборот ... В любом случае, ваши решения работают. Альтернативой может быть использование Interpolations.jl, но необходимая опция для previous интерполяции - все еще ожидает рассмотрения. - person rotton; 21.02.2020
comment
Я был удивлен, увидев эту работу, даже когда первый PresetTimeCallback происходит в начальный момент времени, например tspan[1]. Я работаю над проблемой, когда PresetTimeCallback изменяет состояние, а не параметр. В этом случае обратный вызов не запускается в начальный момент времени. Есть идеи, почему это может иметь значение? - person fabern; 02.07.2020
comment
В любом случае, если это кому-то поможет: во втором случае, когда affect! изменяет состояние, например integrator.u[1] += input[ind_t[1], 2] обходной путь может состоять из двух шагов: a) использовать PresetTimeCallback(dosetimes, affect!, initialize = (c,u,t,integrator) -> affect!(integrator)) (как указано в github.com/SciML/DifferentialEquations.jl/issues/). и б) изменить affect!, чтобы проверить if (integrator.t in dosetimes), иначе findall может вызвать ошибку. - person fabern; 02.07.2020
comment
Привет @fabern, я давно сюда не заглядывал. Я кодировал свои первые шаги в Джулии, чтобы посмотреть, смогу ли я перейти на нее навсегда. Но эта функция прыжка является важным аспектом в моем моделировании, поэтому мне нужно для нее действительно надежное решение. Кажется, что Pumas будет интересной альтернативой. , поскольку их DosageRegimen очень похож на корм, который я использую. В конце концов, я попробую. - person rotton; 03.08.2020