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

Question

Аналогично этот вопрос , я пытаюсь решить этот 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 , но я полагаю, что это будет работать только для непрерывных функций, а не для дискретных изменений !? Я не посмотрел справку по interpolate, но не уверен, как использовать ее для моего конкретного c случая.

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

Answer 1

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

findall(f::Function, A)

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

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

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

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