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")
```