# DiffEqs: Hybrid (continuous+discrete) system / periodic callback

**URL:** <https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791>\
**Category:** Modelling & Simulations\
**Tags:** diffeq\
**Created:** [May 2, 2019, 7:35pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791 "2019-05-02T19:35:31Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 14, 2020, 11:29pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/21 "2020-04-14T23:29:03Z")

</div>

Yeah no worries, I thought I’d lay it all out there because you seem fairly interested in DEDataArrays, and so it’s probably good for me to write down somewhere why I think they may not be a long term solution. But I don’t have a better long term solution, and maybe a simple change is all that’s needed to fix 99.99% of cases, so I’ll probably play around with those ideas I proposed tonight. I really wanted this text somewhere because after I do this PR, I assume that I’ll need to have a post to point to in order to explain to future contributors why it’s so crucial that interpolation uses `similar` on the right and not the left end point, even though it overwrites all of the values. Writing generic code is fun 😉

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 15, 2020, 6:10am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/22 "2020-04-15T06:10:31Z")

</div>

You’re too strict on your typing. Always be careful with types if you want AD to work. Here, AD only occurs on continuous variables, so the types of continuous and discrete variables should be separated.

```julia
mutable struct CtrlSimTypeScalar{T,T2} <: DEDataVector{T}
    x::Vector{T}
    ctrl_x::T2 # controller state
    ctrl_u::T2 # controller output
end

```

with this, your test works now, and it also works on out-of-place and implicit methods:

```julia
  function f(x,p,t)
      x.ctrl_u .- x
  end
  function f(dx,x,p,t)
      dx[1] = x.ctrl_u - x[1]
  end

  ctrl_f(e, x) = 0.85(e + 1.2x)

  x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
  prob = ODEProblem{false}(f, x0, (0.0, 8.0))

  function ctrl_fun(int)
      e = 1 - int.u[1] # control error

      # pre-calculate values to avoid overhead in iteration over cache
      x_new = int.u.ctrl_x + e
      u_new = ctrl_f(e, x_new)

      # iterate over cache...
      if DiffEqBase.isinplace(int.sol.prob)
        for c in full_cache(int)
            c.ctrl_x = x_new
            c.ctrl_u = u_new
        end
      end
  end

  integrator = init(prob, Tsit5())

  # take 8 steps with a sampling time of 1s
  Ts = 1.0
  for i in 1:8
      ctrl_fun(integrator)
      DiffEqBase.step!(integrator,Ts,true)
  end

  sol = integrator.sol
  @test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]

  prob = ODEProblem{true}(f, x0, (0.0, 8.0))

  integrator = init(prob, Rosenbrock23())

  # take 8 steps with a sampling time of 1s
  Ts = 1.0
  for i in 1:8
      ctrl_fun(integrator)
      DiffEqBase.step!(integrator,Ts,true)
  end

  sol = integrator.sol
  @test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]

```

with PRs:

[https://github.com/SciML/DiffEqBase.jl/pull/481](https://github.com/SciML/DiffEqBase.jl/pull/481)  
[https://github.com/SciML/OrdinaryDiffEq.jl/pull/1110](https://github.com/SciML/OrdinaryDiffEq.jl/pull/1110)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 15, 2020, 3:35pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/23 "2020-04-15T15:35:14Z")

</div>

@asprionj your issues are now fixed on the latest release.

---

<div class="post-metadata">

**Author:** ![asprionj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asprionj/32/6856_2.png) [@asprionj](https://discourse.julialang.org/u/asprionj)\
**Post date:** [April 15, 2020, 9:06pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/24 "2020-04-15T21:06:29Z")

</div>

Awesome(ly quick), thanks! 👍

---

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [April 30, 2020, 5:23am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/25 "2020-04-30T05:23:32Z")

</div>

The out-of-place example doesn’t actually do anything - everything stays at 0.0 and `ctrl_fun` has no effect on the integrator.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 30, 2020, 9:42am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/26 "2020-04-30T09:42:17Z")

</div>

which example?

---

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [April 30, 2020, 4:18pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/27 "2020-04-30T16:18:11Z")

</div>

The out-of-place test in the post to which I replied:

```julia
mutable struct CtrlSimTypeScalar{T,T2} <: DEDataVector{T}
    x::Vector{T}
    ctrl_x::T2 # controller state
    ctrl_u::T2 # controller output
end

function f(x,p,t)
      x.ctrl_u .- x
end

ctrl_f(e, x) = 0.85(e + 1.2x)

x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
prob = ODEProblem{false}(f, x0, (0.0, 8.0))

function ctrl_fun(int)
    e = 1 - int.u[1] # control error

    # pre-calculate values to avoid overhead in iteration over cache
    x_new = int.u.ctrl_x + e
    u_new = ctrl_f(e, x_new)

    # iterate over cache...
    if DiffEqBase.isinplace(int.sol.prob)
        for c in full_cache(int)
            c.ctrl_x = x_new
            c.ctrl_u = u_new
        end
    end
end

integrator = init(prob, Tsit5())

# take 8 steps with a sampling time of 1s
Ts = 1.0
for i in 1:8
    ctrl_fun(integrator)
    DiffEqBase.step!(integrator,Ts,true)
end

sol = integrator.sol

# nothing happened
getfield.(sol.u,:x) == [[0.0] for i in 1:length(sol)]

```

I know it’s a test for the interpolation of the discrete variables, but what would make this a fully functional example?

Adding the `else` statement to ctrl\_fun almost makes it work

```julia
function ctrl_fun(int)
  e = 1 - int.u[1] # control error

  # pre-calculate values to avoid overhead in iteration over cache
  x_new = int.u.ctrl_x + e
  u_new = ctrl_f(e, x_new)

  # iterate over cache...
  if DiffEqBase.isinplace(int.sol.prob)
    for c in full_cache(int)
        c.ctrl_x = x_new
        c.ctrl_u = u_new
    end
  else
    int.u.ctrl_x = x_new
    int.u.ctrl_u = u_new
    u_modified!(int,true)
  end
end

```

but the effect is delayed until 1 second.  
 ![outofplace](https://global.discourse-cdn.com/julialang/original/3X/d/f/df12bfdf60b1b7d1a8129e1ecd7ceef1fbe66327.png)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 1, 2020, 11:16am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/28 "2020-05-01T11:16:04Z")

</div>

```julia
mutable struct CtrlSimTypeScalar{T,T2} <: DEDataVector{T}
    x::Vector{T}
    ctrl_x::T2 # controller state
    ctrl_u::T2 # controller output
end

function f(x,p,t)
    x.ctrl_u .- x
end
function f(dx,x,p,t)
    dx[1] = x.ctrl_u - x[1]
end

ctrl_f(e, x) = 0.85(e + 1.2x)

x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
prob = ODEProblem{false}(f, x0, (0.0, 8.0))

function ctrl_fun(int)
    e = 1 - int.u[1] # control error

    # pre-calculate values to avoid overhead in iteration over cache
    x_new = int.u.ctrl_x + e
    u_new = ctrl_f(e, x_new)

    # iterate over cache...
    if DiffEqBase.isinplace(int.sol.prob)
      for c in full_cache(int)
          c.ctrl_x = x_new
          c.ctrl_u = u_new
      end
    end
end

integrator = init(prob, Tsit5())

# take 8 steps with a sampling time of 1s
Ts = 1.0
for i in 1:8
    ctrl_fun(integrator)
    DiffEqBase.step!(integrator,Ts,true)
end

sol = integrator.sol
@test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]

```

passes

---

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [May 1, 2020, 3:39pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/29 "2020-05-01T15:39:11Z")

</div>

The test of the interpolation of `ctrl_u` passes, but the integration itself is incorrect, i.e.

```julia
@test all(u[1] == 0 for u in sol.u)

```

also passes, which is obviously not what we want. The function `ctrl_fun` in its current form does not currently update the integrator when the problem is out-of-place—what is the correct way to do that?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 3, 2020, 6:03am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/30 "2020-05-03T06:03:10Z")

</div>

```julia
using OrdinaryDiffEq, Test

mutable struct CtrlSimTypeScalar{T,T2} <: DEDataVector{T}
    x::Vector{T}
    ctrl_x::T2 # controller state
    ctrl_u::T2 # controller output
end

function f(x,p,t)
    x.ctrl_u .- x
end
function f(dx,x,p,t)
    dx[1] = x.ctrl_u - x[1]
end

ctrl_f(e, x) = 0.85(e + 1.2x)

e = 1 - 0.0
x_new = 0.0 + e
x0 = CtrlSimTypeScalar([0.0], x_new, ctrl_f(e, x_new))
prob = ODEProblem{false}(f, x0, (0.0, 8.0))

function ctrl_fun(int)
    e = 1 - int.u[1] # control error

    # pre-calculate values to avoid overhead in iteration over cache
    x_new = int.u.ctrl_x + e
    u_new = ctrl_f(e, x_new)

    # iterate over cache...
    if DiffEqBase.isinplace(int.sol.prob)
      for c in full_cache(int)
          c.ctrl_x = x_new
          c.ctrl_u = u_new
      end
    else
        int.u = CtrlSimTypeScalar(int.u.x,x_new,u_new)
    end
    u_modified!(int,true)
end

integrator = init(prob, Tsit5())

# take 8 steps with a sampling time of 1s
Ts = 1.0
for i in 1:8
    ctrl_fun(integrator)
    DiffEqBase.step!(integrator,Ts,true)
end

sol = integrator.sol
@test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]
@test any(u[1] != 0 for u in sol.u)

using Plots
# extract solution, plot...
sol = integrator.sol
plot(sol, label="system")

# actual full steps, with steppre line-style ("correct")
plot!(sol.t, [u.ctrl_u for u in sol.u], linetype=:steppre, color=:black, marker=true, label="u[.].ctrl_u")

# interpolated control signal ("wrong")
t_ctrl = 0:0.01:sol.t[end]
plot!(t_ctrl, [sol(t).ctrl_u for t in t_ctrl], color=:red, label="u(t).ctrl_u")

```

is what you want, and at the same time, further highlights why you don’t want to do this.

---

<div class="post-metadata">

**Author:** ![jt1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jt1/32/38331_2.png) [@jt1](https://discourse.julialang.org/u/jt1)\
**Post date:** [March 5, 2021, 8:29pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/31 "2021-03-05T20:29:36Z")

</div>

Hey @ChrisRackauckas  
Are the methods talked about here, the most efficient ones?

I benchmarked 2 code snippets from this thread in the code underneath and they seem to have many more allocations than expected (more than 2000 for only 34 timesteps). This while the integration and callback functions have no allocations.

```julia
using DifferentialEquations
using BenchmarkTools

function normal_integration()
    f(x, p, t) = p .- x # the system ("unit first-order lag")
    ctrl_f(e, x) = 0.85(e + 1.2x) # the controller (proportional-integral)
    ctrl_x = [0.0] # array to record (discrete) state of controller

    function ctrl_cb(int) # the callback function
        e = 1 - int.u[1] # control error
        push!(ctrl_x, ctrl_x[end] + e) # integrate the control error
        int.p = ctrl_f(e, ctrl_x[end]) # update controller output
    end

    pcb = PeriodicCallback(ctrl_cb, 1.0)
    prob = ODEProblem(f, [0.0], (0.0, 8.0), 0.0, callback=pcb)
    sol = solve(prob)
end

## With DEDataArray

mutable struct CtrlSimTypeScalar{T} <: DEDataVector{T}
    x::Array{T,1}
    ctrl_x::T # controller state
    ctrl_u::T # controller output
end

function dedata_integration()
    function f(dx,x,p,t)
        dx[1] = x.ctrl_u - x[1]
    end

    ctrl_f(e, x) = 0.85(e + 1.2x)

    x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
    prob = ODEProblem(f, x0, (0.0, 8.0))

    function ctrl_fun(int)
        e = 1 - int.u[1] # control error

        # pre-calculate values to avoid overhead in iteration over cache
        x_new = int.u.ctrl_x + e
        u_new = ctrl_f(e, x_new)

        # iterate over cache...
        for c in full_cache(int)
            c.ctrl_x = x_new
            c.ctrl_u = u_new
        end
    end

    integrator = init(prob, Tsit5()) # only Tsit5 seems to work

    # take 8 steps with a sampling time of 1s
    Ts = 1.0
    for i in 1:8
        ctrl_fun(integrator)
        step!(integrator,Ts,true)
    end
    integrator.sol
end

@btime normal_integration() #108.900 μs (2026 allocations: 185.73 KiB)
@btime dedata_integration() #160.600 μs (2140 allocations: 186.73 KiB)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 8, 2021, 6:03am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/32 "2021-03-08T06:03:18Z")

</div>

It’s not going to be optimal because it’s using a mutable struct in a way that will allocate. ModelingToolkit.jl with removing of observables will be better.

---

<div class="post-metadata">

**Author:** ![asprionj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asprionj/32/6856_2.png) [@asprionj](https://discourse.julialang.org/u/asprionj)\
**Post date:** [March 10, 2021, 4:08pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/33 "2021-03-10T16:08:04Z")

</div>

The logging of the signals seems not to be the time-consuming part. Modyfing `normal_integration()` to just use a scalar to store the controller state:

```julia
no_log_integration()
...
    ctrl_x = 0.0 # state of controller (error integrator)

    function ctrl_cb(int) # the callback function
        e = 1 - int.u[1] # control error
        ctrl_x += e
        int.p = ctrl_f(e, ctrl_x) # update controller output
    end
...
end

```

On my machine gives:

```julia
@btime normal_integration() #137.600 μs (2044 allocations: 185.88 KiB)
@btime dedata_integration() #187.400 μs (2305 allocations: 205.55 KiB)
@btime no_log_integration() #140.499 μs (2077 allocations: 186.22 KiB)

```

Running the `solve` with `save_on = false` does not drastically reduce the runtime and allocations:

```julia
132.301 μs (1858 allocations: 164.41 KiB)

```

But overall, `DEDataArray` indeed seems to add quite some overhead.  
How would the example here be implemented using ModelingToolkit.jl? Would be interesting to compare that as well! (The implementation itself, and the performance.)

**Edit** : @ChrisRackauckas, AFAIK there is no integration of such discrete events / systems in ModelingToolkit, correct? So MTK would just be used to model the “physical” / continuous part of the system, then use one of the above approaches to “buckle” the discrete-time controller onto that model?

---

<div class="post-metadata">

**Author:** ![jt1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jt1/32/38331_2.png) [@jt1](https://discourse.julialang.org/u/jt1)\
**Post date:** [March 15, 2021, 8:54am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/34 "2021-03-15T08:54:53Z")

</div>

I think maybe a more natural way of mixing discrete and continuous systems than the above methods is to add the discrete state as an extra state equation to the physical system like: `dctrl_x = 0`  
and then update the `ctrl_x` and `int.p` using a `PeriodicCallback`. Then atleast the `ctrl_x` is also stored in `sol.u` and you can use `EnsembeModels` etc. with this system more easily.

This isn’t very modular though as changing the amount of states in the controller requires you to change the equations of the physical system.  
But i guess `ModelingToolkit` might let you describe difference equations for the discrete systems at some point and do something like this automatically.

---

<div class="post-metadata">

**Author:** ![asprionj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asprionj/32/6856_2.png) [@asprionj](https://discourse.julialang.org/u/asprionj)\
**Post date:** [March 15, 2021, 10:21am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/35 "2021-03-15T10:21:58Z")

</div>

Adding additional states: wouldn’t this just result in the same overhead and rather cumbersome “update” as with `DEDataArray`s? I.e. having to loop over some cache.

@zekeriya.sari: How are discrete (sub-)systems integrated into the overall simulation in [Causal.jl](https://github.com/zekeriyasari/Causal.jl)? You also rely on DiffEqs.jl for the solution, and allow multiple discrete subsystems with different sampling times, right?

ModelingToolkit.jl: Yes, I’d also really like to know whether there is…

- Already a way of doing this,
- Something planned, or
- An Idea how this could / should be realized, as a starting point?

---

<div class="post-metadata">

**Author:** ![jt1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jt1/32/38331_2.png) [@jt1](https://discourse.julialang.org/u/jt1)\
**Post date:** [March 15, 2021, 12:57pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/36 "2021-03-15T12:57:42Z")

</div>

Yeah you are right, it is indeed slower (see code). But that might also partly be because in this case, it changes a scalar ODE to a non-scalar ODE and requires indexing for the system state where it didn’t before.

And thank you @asprionj for making me aware of [Causal.jl](https://github.com/zekeriyasari/Causal.jl). I didn’t know it existed.

```julia
using DifferentialEquations
using BenchmarkTools

function normal_integration()
    f(x, p, t) = p .- x # the system ("unit first-order lag")
    ctrl_f(e, x) = 0.85(e + 1.2x) # the controller (proportional-integral)
    ctrl_x = [0.0] # array to record (discrete) state of controller

    function ctrl_cb(int) # the callback function
        e = 1 - int.u[1] # control error
        push!(ctrl_x, ctrl_x[end] + e) # integrate the control error
        int.p = ctrl_f(e, ctrl_x[end]) # update controller output
    end

    pcb = PeriodicCallback(ctrl_cb, 1.0)
    prob = ODEProblem(f, [0.0], (0.0, 8.0), 0.0, callback=pcb)
    sol = solve(prob)
end

## With DEDataArray

mutable struct CtrlSimTypeScalar{T} <: DEDataVector{T}
    x::Array{T,1}
    ctrl_x::T # controller state
    ctrl_u::T # controller output
end

function dedata_integration()
    function f(dx,x,p,t)
        dx[1] = x.ctrl_u - x[1]
    end

    ctrl_f(e, x) = 0.85(e + 1.2x)

    x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
    prob = ODEProblem(f, x0, (0.0, 8.0))

    function ctrl_fun(int)
        e = 1 - int.u[1] # control error

        # pre-calculate values to avoid overhead in iteration over cache
        x_new = int.u.ctrl_x + e
        u_new = ctrl_f(e, x_new)

        # iterate over cache...
        for c in full_cache(int)
            c.ctrl_x = x_new
            c.ctrl_u = u_new
        end
    end

    integrator = init(prob, Tsit5()) # only Tsit5 seems to work

    # take 8 steps with a sampling time of 1s
    Ts = 1.0
    for i in 1:8
        ctrl_fun(integrator)
        DifferentialEquations.step!(integrator,Ts,true)
    end
    integrator.sol
end

## Augmenting the state

function augmented_integration() 
    function f(dx, x, p, t)# the system ("unit first-order lag")
        dx[1] = p .- x[1]
        dx[2] = 0.0 #This isn't really needed i guess if x0[2] = 0
        dx
    end

    ctrl_f(e, x) = 0.85(e + 1.2x) # the controller (proportional-integral)

    function ctrl_cb(int) # the callback function
        e = 1 - int.u[1] # control error
        next_ctrl_x = int.u[2] + e # integrate the control error
        int.p = ctrl_f(e, next_ctrl_x) # update controller output
        for c in full_cache(int)
            c[2] = next_ctrl_x
        end
    end

    pcb = PeriodicCallback(ctrl_cb, 1.0)
    prob = ODEProblem(f, [0.0,0.0], (0.0, 8.0), 0.0, callback=pcb)
    sol = solve(prob)
end

@btime normal_integration() #108.900 μs (2026 allocations: 185.73 KiB)
@btime dedata_integration() #160.600 μs (2140 allocations: 186.73 KiB)
@btime augmented_integration() #170.900 μs (1343 allocations: 120.23 KiB)

```

---

<div class="post-metadata">

**Author:** ![zekeriya.sari](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zekeriya.sari/32/13695_2.png) [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)\
**Post date:** [March 16, 2021, 7:39am UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/37 "2021-03-16T07:39:25Z")

</div>

> [@asprionj](#):
>
> @zekeriya.sari: How are discrete (sub-)systems integrated into the overall simulation in [Causal.jl](https://github.com/zekeriyasari/Causal.jl)? You also rely on DiffEqs.jl for the solution, and allow multiple discrete subsystems with different sampling times, right?

Yes @asprionj , Causal.jl depends on DifferentialEquations.jl for diffeq solving. But, the thing is that Causal.jl is not ready for hybrid systems yet. In its current version, it expects the components of the models be complete continuous-time or discrete-time. The main obstacle in the simulation of hybrid systems in Causal.jl is that in Causal.jl the components of the model share a single clock while evolving. In case of fully discrete-time or fully continuous-time systems, this is what we want: all the components evolve with respect to a unique time reference, i.e. the model clock. However, in hybrid systems, the discrete events (modeled by discrete-time systems) need a separate time reference, which implies a single clock can not be used and the discrete-time and continuous time components need their separate clocks. Simulation of hybrid systems is definitely a feature that must be included in Causal.jl and I’m working on it.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 16, 2021, 3:08pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/38 "2021-03-16T15:08:32Z")

</div>

> [@zekeriya.sari](#):
>
> However, in hybrid systems, the discrete events (modeled by discrete-time systems) need a separate time reference, which implies a single clock can not be used and the discrete-time and continuous time components need their separate clocks. Simulation of hybrid systems is definitely a feature that must be included in Causal.jl and I’m working on it.

That can be handled via the callback systems, syncing the clock by having the discrete systems add tstops to the others.

> [@asprionj](#):
>
> - Already a way of doing this,
> - Something planned, or
> - An Idea how this could / should be realized, as a starting point?

I did a full write up here:

> <https://github.com/SciML/ModelingToolkit.jl/issues/894#issuecomment-799599521>
>
> There is a discussion about how to implement hybrid system, i.e. a combination o…f a continuous system (in my domain of control-systems engineering often called the "plant" or the "physical system") with a discrete-time system ("controller"), \[on Discourse\](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791). The question came up whether MTK does already, will in the future, or should at all support discrete-time systems.
> 
> Please don't cosider this a "request to implement it", but a discussion on whether it makes sense, how it could be done, what effort it would take, etc. I know too little about the inner workings of MTK (and DiffEqs.jl) to judge any of this.

Basically, there are ways to work around it right now. We will be fully automating it, but it will take a little bit. And I put a suggestion in there for how to hand code it in a way that’s the most performant

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [August 16, 2021, 2:25pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/39 "2021-08-16T14:25:11Z")

</div>

Sorry to resurrect an old topic. I just want to register my findings.

I use `DEDataVector` a lot. In my case, solutions like extending the state with dummy variables would not be feasible, and that’s why:

The main project here using this framework is a satellite simulator. In this case, we use `DiscreteCallback` to simulate the full attitude control software. We have more than 2,000 parameters, including variables that are boolean, integers, etc. The idea is that we modeled all the embedded software functionality as precisely as we can. Thus, we have a workspace, which is `<: DEDataVector`, containing the continuous variables and a lot of other structures like `sensors`, `actuators`, `controller`, etc. that models the embedded software.

Since `DEDataVector` was marked as deprecated, I was searching for a possible replacement. Right now, ModelingToolkit seems to have some issues that will not allow me to migrate everything to it [1-3]. `DEDataVector` is working fine as long as you follow the instructions in the documentation and perform a full cache update. I tried with many solvers.

Thinking into the future, if `DEDataArrays` is really removed from the ecosystem (I might have problems internally way before this happens since this package is marked as deprecated), then I found one solution that seems to work for me. In the parameters, I pass two variables `(w::Workspace, vw::Vector{Workspace})`. The first one is the latest set of variables updated by the onboard software, and the latter is a vector containing the history of the workspace. Hence, instead of updating the cache at the end of the callback, I push a copy of `w` into `vw`, leading to a value for every `t` in my `tstops`. I did not see any modification in the solution, and it is **way** faster, probably due to reduced allocation. The only downside is that the plotting of discrete variables is not as nice as we get using `DEDataArrays`, which really shows that discrete behavior.

[1] [https://github.com/SciML/ModelingToolkit.jl/issues/1194](https://github.com/SciML/ModelingToolkit.jl/issues/1194)  
[2] [https://github.com/SciML/ModelingToolkit.jl/issues/1196](https://github.com/SciML/ModelingToolkit.jl/issues/1196)  
[3] [https://github.com/SciML/ModelingToolkit.jl/issues/1203](https://github.com/SciML/ModelingToolkit.jl/issues/1203)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 17, 2021, 9:26pm UTC](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/40 "2021-08-17T21:26:03Z")

</div>

The issue with DEDataArrays is more about correctness. They have embedded data inside of there, so what should happen with `x .= y .+ z`? This is well-defined on the iterated values, not well-defined on the extra data. So what it does is it just keeps the data of `x`. There are many integrators for which this setup works for (basically all of the non-stiff ones), which is why it “works”. But it’s not a good general solution because it’s essentially letting the correctness of the program be decided by undefined behavior. For example what should `A*x` for `A` a matrix and `x` a DEDataArray be? Maybe it should propagate the same metadata? But if it’s `A` a DEDataMatrix and `x` a `DEDataVector`, then which metadata do you choose? There are so many of these corner cases that this design is just simply flawed. By “deprecated” I mean, I’ll continue to keep the repo working and the current tests passing, but there won’t be any more development of them. ModelingToolkit has better defined semantics which will give more optimizations and more correctness guarantees, that’s just overall better, and thus that’s the direction this should all be going. That said, anyone is free to continue using DEDataArrays via the DEDataArrays.jl package, but I don’t plan on tracking down hacks for every combination and reordering ODE solver evaluations to make the undefined behavior work in every novel case (in fact, I think that it there’s some cases where there is no ordering which always gives the results someone wants, so it’s just not doable in general).

[Previous page](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791.md?page=1)

[Next page](https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791.md?page=3)