DifferentialEquations.jl discrete delay problem?

In DifferentialEquations.jl (abbrev.: DiffEq.jl), DDE (delay differential equations) is supported.

But what I want is close to delay discrete equation; that is, get an external signal h(t), then the state x(t) = h(t - \tau) with delay constant (\tau > 0).

What is the best practice in this case?

Note: I would like to keep using DiffEq.jl rather than utilising other packages.

One remedy would be to literally define a function x(t) = h(t - tau) by receiving the external signal function h and calling x within ODE function.

Is there any other way?

The equation \dot x(t) = h(t - \tau) is an ordinary differential equation.
(Only if you access the state x in the past, i.e. x(t - \tau), you have a delay differential equation.)

Is h a function with a closed-form expression? Or is it something where you only know h(t_i) (with t_i being the timesteps?)

Not \dot{x}(t) = h(t-\tau) but x(t) = h(t - \tau).

It is possible to access the values of h at any instants, but it would be desirable if the access is limited only at the time instants when solver calls. For now, suppose that we can access the values of h any desired time instants.

the following remedy would work properly, but for the extensibility of my custom functions, I hope that the delayed state x(t) can be dealt with as the solution of an ODE.

That’s not like a differential equation at all?

It looks like just copying values? What is the complete problem you want to solve?

1 Like

Thanks for your attention, @SteffenPL !

Yes, it’s actually now a DE :slight_smile:
To avoid confusion, I’ll explain my situation below:

I’m trying to add a fault detection and isolation (FDI) module in dynamical system. Basically, this receives a true effectiveness signal h(t) and produces delayed (and possibly noisy) estimation of the effectiveness (namely, \hat{h}(t) ).
There are many models of the FDI. In my mind, the implementation of two simple FDI models are the primary task.

  1. (Simple delay model) \hat{h} (t) = h(t-\tau)
  2. (Low-pass-filter-like model) \dot{\hat{h}} (t) = (h(t) - \hat{h}(t)) / \tau

And I wanna deal with these models in the same manner. For the second model (and other parts of dynamical systems, e.g., multicopter in my case), my code is highly based on DifferentialEquations.jl.
In this regard, I wanna implement the first model as a part of differential equation.

That’s why I’m looking forward to ways of realising FDI models in ODE manner (although it would not be expressed as an ODE).

1 Like

Have you tried this? Discrete Problems · DifferentialEquations.jl
You could define the ODE function as

u0 = [0]
p = ( h(t) = t^2, 0.2 )
tspan = (0, 100)

function f(du, u, p, t) 
  h, tau = p
  du[:] .= h(t - tau)

prob = DiscreteProblem(f, u0, p, tspan)

(I haven’t tried the code, it’s just an idea)

Is it ok for you if you define the two different models with different RHS functions?

Can I compose DiscreteProblem and ODEProblem?
If so, it may be possible.

Roughly speaking,

  1. For the first model
x0 = [0]
tspan = (0, 100)
h = t -> t
function f(dx, x, p, t)
    dx .= (h(t) - x) / tau
prob = ODEProblem(f, x0, tspan)
  1. For the second model
x0 = [0]
tspan = (0, 100)
h = t -> t
function f(dx, x, p, t)
    dx .= h(t - tau)
prob = DiscreteProblem(f, x0, tspan)

Then, it seems hard to be dealt with the same manner; when changing models, I may have to write codes twice with different DEProblem

The discrete delay equation as you call it is usually known as a difference equation (assuming you are interested in continuous time; the DiscreteProblem appears to be for discrete time of I’m reading it correctly). My best guess for implementing it with DiffEq is to use DelayDiffEq with a singular mass matrix (you need to make sure that you pick a solver that is compatible with singular mass matrices). This will allow you to mix differential terms and difference terms directly.

I’m pretty sure there is an example in the DiffEq docs of using singular mass matrices. (I’m on my phone otherwise I’d look it up.)


@dawbarton Your solution is exactly what I’ve wanted!

I saw Handling Mass Matrices and adopt the code for my case.

  • Example code
using DifferentialEquations
using LinearAlgebra

function test()
    function dynamics!(du,u,p,t)
        du .= u - ones(size(u)) * t
    n = 3
    x0 = zeros(n)
    M = zeros(n) |> Diagonal |> Matrix
    f = ODEFunction(dynamics!, mass_matrix=M)
    tspan = (0, 10.0)
    prob = ODEProblem(f, x0, tspan)
    sol = solve(prob)
  • Result

julia> test()
retcode: Success
Interpolation: specialized 3rd order "free" stiffness-aware interpolation
t: 9-element Vector{Float64}:
u: 9-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0]
 [1.0e-6, 1.0e-6, 1.0e-6]
 [1.1e-5, 1.1e-5, 1.1e-5]
 [0.00011099999999999999, 0.00011099999999999999, 0.00011099999999999999]
 [0.0011109999999999998, 0.0011109999999999998, 0.0011109999999999998]
 [0.011110999999999996, 0.011110999999999996, 0.011110999999999996]
 [0.11111099999999996, 0.11111099999999996, 0.11111099999999996]
 [1.1111109999999995, 1.1111109999999995, 1.1111109999999995]
 [10.0, 10.0, 10.0]


  • This is independent of the initial state (you can choose any x0). Even nothing can be assigned.

For automatic composition, it’s coming through ModelingToolkit, and you can follow this issue for ongoing information:

2/4 chunks have been knocked out from it. I assume this is for optimal control, and yeah :sweat_smile: I’ll make ya’ll happy in a bit.

1 Like

Thanks a lot!

Yes, it is quite often necessary to compose discrete- and continuous-time dynamical systems (especially for discrete-input controller, e.g., model predictive control (MPC)).

As a student of control engineering, SciML is a very rich ecosystem for various kinds of simulation! :slight_smile:
I would be so pleased if the composition functionality is once provided.

EDIT: for other control engineers, I highly recommended you to take a look at ComponentArrays.jl to construct nested and complex dynamical systems.
Also, I’ve tried to find a good way to construct dynamical system and control simulator; FlightSims.jl.

Yeah the whole MPC story is coming together. The extra pieces are all about feeding the solver extra information, new solvers specifically designed for better handling of discontinuities, etc… next year will be a fun one around optimal control. But until then, I’d recommend using DiffEq directly, and yes ComponentArrays.jl is an amazing abstraction for making the “by hand” approach very clean (I plan to add ComponentArrays.jl to the DifferentialEquations.jl first ODE tutorial to promote it even more).


BTW, I have a plan of implementing discrete-time control via DiscreteCallback by considering parameters of parameterised ODE as control input.
Will there be a distinction more than the difference between state variable of ODE (upcoming supports as I understand) and parameters of parameterised ODE (my plan)?