Constructing callbacks on observed variables with symbolic indexing of an MTK model

Specific Domains Modelling & Simulations
, ,
1

I am trying to figure out how to use the Symbolic indexing with callback functions. There is already the continuous events support within MTK which is quite nice and has worked for this MWE. When the variable goes to observed, however, I can no longer reference the variable within the event and must make it irreducible (way1). This can be undesirable since it makes the ODE into a DAE. I can also build the callback via ContinuousCallback as shown in the below MWE, however, it skips over the termination event for the first two cases. It does terminate for the last case.

I’d appreciate any input on personal experiences with building continuous callbacks on observed variables, and in the case of much more nonlinear terminations if there may be a better way to approach the problem. Thanks

MWE:

using ModelingToolkit
import ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq

function build_system(; name=:sys, way=:way1)
    ps = @parameters begin
        I = -30
        R_p = 2.0e-6
        F = 96485
        a = 723600
        L = 8.8e-5
        D_s = 3.9e-14
        c_s_max = 30555
    end

    states = @variables begin
        c_s(t)
        c_s_avg(t), [guess=0.0]
        SOC(t), [guess=(way==:way1 || way==:way2 || way==:way3) ? 0.0 : nothing, irreducible=(way==:way1 || way==:way2 || way==:way3)]
    end

    eqns = [D(c_s_avg) ~ -3 * (I / (F * a * L)) / R_p,
            c_s - c_s_avg ~ -R_p / D_s * (I / (F * a * L)) / 5,
            SOC ~ c_s_avg / c_s_max]

    # way 1 = continuous event + irreducible -> works
    cont_event_eq = [SOC ~ 1.0] => ((integrator, u, t, ctx) -> terminate!(integrator), [SOC], [c_s_max], [], nothing)

    return structural_simplify(ODESystem(eqns, t, [states...], ps, name=name, continuous_events= way==:way1 ? cont_event_eq : nothing))
end

sys = build_system()

function solve_problem(sys; way=:way1)
    prob = ODEProblem(sys, ModelingToolkit.get_guesses(sys), (0.0, 5000.0), ModelingToolkit.get_defaults(sys))

    # way 2 integ[:SOC] - 1.0 does not terminate
    # way 3 integ.u[ModelingToolkit.variable_index(integ.f.sys, :SOC)] - 1.0 does not terminate
    # way 4 integ.sol(t, idxs=:SOC) does terminate
    cont_callback = ContinuousCallback(
        (u, t, integ) -> way==:way2 ? (integ[:SOC] - 1.0) :
                         way == :way3 ? integ.u[ModelingToolkit.variable_index(integ.f.sys, :SOC)] - 1.0 :
                         integ.sol(t, idxs=:SOC) - 1.0,
        (integ) -> terminate!(integ)
    )
    
    sol = solve(prob, QNDF(), callback = way != :way1 ? cont_callback : nothing)
end

function build_and_solve(; way=:way1)
    sys = build_system(; way = way)
    sol = solve_problem(sys; way = way)
    return sol
end

sol = build_and_solve(way = :way4)
2 
  [0c46a032] DifferentialEquations v7.16.0
  [98e50ef6] JuliaFormatter v1.0.62
⌃ [961ee093] ModelingToolkit v9.64.0
  [429524aa] Optim v1.11.0
  [1dea7af3] OrdinaryDiffEq v6.91.0
  [d96e819e] Parameters v0.12.3
  [14b8a8f1] PkgTemplates v0.7.53
  [91a5bcdd] Plots v1.40.9
3

In the case with a strong nonlinearity, the final value of the callback is dictated by the degree of nonlinearity of the step function (in this code, the variable k controls sharpness of transition). When x_2 reaches x_2_max, the nonlinearity is effectively “activated” and the affect within the callback should, and does, occur when the value = 25.0

For different k values, however,
k = 50, end nonlinearity = 24.9999966995029
k = 500, end nonlinearity = 24.999966995029023
k = 500, end nonlinearity = 24.999669950290237

using ModelingToolkit
import ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq

function build_system(; name=:sys)#, way=:way1)
    ps = @parameters begin
        x_1_rate = .05
        x_1_x_2_diff = 5
        x_2_max = 200
        k = 500
        α = 0.0
    end

    states = @variables begin
        x_1(t), [guess=0.0]
        x_2(t)
        nonlinearity(t)
    end

    function nonlinearity_function(x, xmax, k, alpha)
        ratio = x / xmax
        return 50 - 50 / (1 + exp(-k * (ratio - (1 - alpha))))
    end

    eqns = [D(x_1) ~ x_1_rate,
            x_2 - x_1 ~ x_1_x_2_diff,
            nonlinearity ~ nonlinearity_function(x_2, x_2_max, k, α)]

    return structural_simplify(ODESystem(eqns, t, [states...], ps, name=name))
end

function solve_problem(sys)
    prob = ODEProblem(sys, ModelingToolkit.get_guesses(sys), (0.0, 10000.0), ModelingToolkit.get_defaults(sys))

    cont_callback = ContinuousCallback(
        (u, t, integ) -> integ.sol(t, idxs=:nonlinearity) - 25.0,
        (integ) -> terminate!(integ),
        save_positions = (true, true),
        initializealg=CheckInit(),
        reltol=1e-9, abstol=1e-9,
        interp_points=0
    )
    
    sol = solve(prob, QNDF(), callback = cont_callback)
end

function build_and_solve()
    sys = build_system()
    sol = solve_problem(sys)
    return sol
end

sol = build_and_solve()

sol[:nonlinearity][end]