Callbacks causing wrong solution mixed ODE and Linear system

[Lode]

Hello,

I am new in Julia and tried to model a system with a discrete stepping feedback controller.

The system was modelled using ‘ModelingToolkit’ and thereafter converter to an ‘ODEProblem’.

The discrete feedback controller was given as a function and insert in the system by a ‘DiscreteCalback’.

A part of the variables of the system is equal to the feedback of the controller another part is given by differential equations.

When testing it, I noticed that the solution of the variables (of the system) equal to the feedback of the controller are returned as constant (equal to the value of the last time of the simulations) over the simulations time, but the other variable (given by an ODE relations to this variables) are correctly solved taken into account this stepping effect.

As given in the (simplified) code example below .

Thank you for your reply

using ModelingToolkit
using OrdinaryDiffEq
using Plots
# define basic ODEProblem
@parameters t p q
@variables u1(t) u2(t)
D = Differential(t)
eq= [D(u1) ~ -0.5*u1*p - u2,
      u2 ~ q]
@named f_test = ODESystem(eq, t; name =:f_test)
smp_sys = structural_simplify(f_test)
parameters(smp_sys)
u0_sys = [u1 => 10.0,
          u2 => 10.0]
p_sys = [p => 0.0, q => 0.0]
prob_sys = ODEProblem(smp_sys, u0_sys, (0.0, 10.0), p_sys)

# define callbacks 
const tstop1 = [5.0]
const tstop2 = [8.0]

function condition(u, t, integrator)
    t in tstop1
end

function condition2(u, t, integrator)
    t in tstop2
end

pos = 2
function affect!(integrator)
    integrator.p[pos] = 5
end

function affect2!(integrator)
    integrator.p[2] = -5
end

save_positions = (true, true)

cb = DiscreteCallback(condition, affect!, save_positions = save_positions)

save_positions = (false, true)

cb2 = DiscreteCallback(condition2, affect2!, save_positions = save_positions)

cbs = CallbackSet(cb, cb2)
const tstop = [5.0; 8.0]
# solve combined system
sol = solve(prob_sys, Tsit5(), callback = cbs, tstops = tstop)
t_sol = sol[t]
u1_sol = sol[u1]
u2_sol = sol[u2]
# Ploting
t_u2_expected = [0.0, 4.99999, 5.00001, 7.9999, 8.0001, 10]
function u2_expected_calc(t_sol)
    sol = zeros(length(t_sol))
    for i in 1:length(t_sol)
        sol[i] = (t_sol[i] < tstop1[1]) ? 0 : (t_sol[i] < tstop2[1] ? 5 : -5)
    end
    return sol
end

u2_expected = u2_expected_calc(t_u2_expected)
plot(t_sol, [u1_sol, u2_sol], title = "Wrong u₂ solution", label = ["u₁ solution solver" "u₂ solutions solver"]);
plot!(t_u2_expected, [u2_expected], label = "u₂ 'real' solution")

[Sevi]

Hi @Lode and welcome to the forum!

Could you please clarify what system you are trying to simulate? From the code, it looks like you have these equations you are solving:

\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} u_1 &= - \frac{1}{2} u_1 p - u_2 \\ u_2 &= q \end{align}

Is this correct (note that there is just u_2 and not the derivative in the second equation)? I would expect the solution to look like the blue and green plot that you showed (the parameter p stays zero and u_2 is set to a constant value at some point)

[Lode]

Thank you for your reply @Sevi,

The system here doesn’t represent anything, it is more like a try-out case to see verry quick if the code blocks works as expected since the main code doesn’t work as expected. In practice, I am making a small toolbox to simulate batteries based on equivalent schemes (in a nutshell this like an electric circuit where the components have changing values, for more detail see here). This system will be mixed of (non-)linear and differential equations.

The controller will be directly control the charge currents in order to keep the battery in his design limits (maximum/minimum voltage, current, power etc.). It must be easily to change since there are numberless strategies. The charge current is a variable that also must be saved.

The write out equations as in previous post are indeed correct and p is fixed to zero.

But my main problem lies in the solver it-self, as in previous post already said you would expect that he return the green and blue curve. But instead for u₂, he is given back the orange curve instead of the green curve. The solver knows the changes of q since he adopt them in the u₁ solution, this is easily visible in the moments between 0 and 5 where u₁ is fixed due to well-chosen values of p and q. If the solver would use the orange curve of u₂ in all moments u₁ would not fixed, but instead continous linear decreasing like in the last moments.

[Sevi]

Apologies, I mixed up the plots. Yeah, it is strange that the “solution” looks like that.

It seems like the mechanism that turns u2 into an “observed” variable (because it’s not really a dynamic variable, it is not part of the variables of the ODE system in the end), doesn’t actually record the values, but assumes that the value is always the one at the final time

At least that would make sense if the observed variable is just equal to a parameter which is (?) assumed to be constant.

Here’s what I did to find out more about what’s going on, in case it helps you or someone else.

Details

The code was run with Julia 1.10.2 and

  [961ee093] ModelingToolkit v9.5.0
  [1dea7af3] OrdinaryDiffEq v6.74.0
  [91a5bcdd] Plots v1.40.2

• check what are the actual variables and observed quantities of the system (the variable u2 is correctly identified as one that is not actually dynamic)

julia> unknowns(smp_sys)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 u1(t)

julia> observed(smp_sys)
1-element Vector{Equation}:
 u2(t) ~ q

• use Chtulhu to “descend” into the function call that extracts the solution (I find this very useful to quickly check which functions are called, especially when there is a lot of dispatch going on – since the types here have a lot of parameters, you might want to turn off the type information by pressing t while in the descend mode)

julia> using Cthulhu
julia> @descend sol[u2]
# lots of output here, press `t` to turn off all the type info and see the source code

# after a few more descends, we arrive at

function _getindex(A::AbstractDiffEqArray, ::ScalarSymbolic, sym)
...
     elseif is_observed(A, sym)
         return observed(A, sym).(A.u, (parameter_values(A),), A.t)
...
end

• the function above looks like a broadcastet call of whatever observed(A, sym) returns over the values of u and the times. The parameter values in (parameter_values(A),) are just the values at the end of the integration though

My best guess is that dynamically changing parameters should also be labelled as dynamic in some way. Otherwise, the values of the “observed” quantities is assumed to only depend on the values of the variables and the (constant) parameters.

Perhaps it is possible to use a time-dependent function instead of the parameter on the right-hand side of your ODE?