Discrete-Time (Sub-)Systems

[asprionj]

There is a discussion about how to implement hybrid system, i.e. a combination of 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. 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.

[ChrisRackauckas]

It makes sense, and... close to duplicate but more so it's a combination of issues with a small cherry on top. And indeed, that Discourse issue describes why I don't want to waste any more time on DEDataArray: it's a dead end with some fundamental issues because of how the discrete data is "supposed" to propagate in broadcast, but that ends up making the propagation of controls information non-deterministic. We can make Tsit5 work with it and keep it tested, but it's not a good general solution. The general solution has to be MTK, with its observation equations, and all of that. So I want to focus the discussion on getting a real solution, and real workarounds until that's done.

(Don't feel bad about the duplicate issues though: it probably is good for populating the search term with different titles.)

The issues that you're touching on here are:

Difference operator: JuliaSymbolics/Symbolics.jl#32
Discrete system: #571
Events: #38
Promotion: #817

So there's a few steps to get there, at least to do it correctly. We want to make sure the whole symbolic system is directly supportive of it. You can hack it in today, but that won't mean you'll get great symbolics on it! So the steps are:

• Define a differencing operator so that we can symbolically represent a discrete-time equation. Define some symbolics for expanding and solving them.
• Define DiscreteSystem so that we have a symbolic analog to DiscreteProblem. Define all of the conversions. This is fairly routine and a great first issue BTW: it's just a bit more than a copy-paste of some ODESystem stuff.
• Define a way to symbolically specify and build events. @YingboMa is tackling that one, but it'll take a few months.
• Define a promotion system on subsystems, so you can put a NonlinearSystem inside of an ODESystem
• Inside of ODESystem, promote differencing equations to DiscreteCallbacks

But yes, let me indeed state that we want to solve the grand problems mentioned in https://www.mathworks.com/videos/accelerating-the-pace-and-scope-of-control-system-design-1493659266810.html in a way that's fully generalized to acausality, stochasticity, delays, PDEs, and more. So making it nice to do control systems is a big issue.

The workaround you should do today is just add an extra state with D(x) ~ 0, and then directly add a DiscreteCallback on the generated code.

Now, there is a very interesting question here of how to represent the control data in a way that can be accessed for post plotting. Here's a few things:

1. You could extend the system with a D(x) ~ 0, but you can't simplify that out if you want to change it with a callback.
2. You could add it to the parameters and modify that parameter. Then you don't have a way to get the time series.
3. You could have different parts of the parameters. "Parameters" and "caches". Controls could create a ZeroSpline (from https://github.com/PumasAI/DataInterpolations.jl) which it appends to in the callback while it does the change.

@YingboMa this is a big one to think about.

[ChrisRackauckas]

From a discussion with Martin and Hilding we settled on (3).

[jonniedie]

I'm not sure if there's anything set up for this already, but is there a plan in place for handling discrete samples from continuous variables as well?. What I think would be really neat syntax for this is something like:

@variables t
D = Differential(t)
k = SampledTime(t, 0.01)

@variables x(t), u(k)

eqs = [
    D(D(x)) = x + u(k)
    u(k+1) = -3x(k) - 2x(k-1) + 5u(k) - 3u(k-1)
]

This also helps with the fact that different discrete time variables might have different sample rates. You could have different SampledTimes for any discrete control variable in your system each one of them would be able to pull previous sampled values of the continuous variables at the time points they need to.

[ChrisRackauckas]

See step Define a differencing operator so that we can symbolically represent a discrete-time equation. Define some symbolics for expanding and solving them.

[asprionj]

From a discussion with Martin and Hilding we settled on (3).

Will you start this soon, or already did? Can I help somehow?
Guess this feature interrelates with quite some internals of MTK (and DiffEq?), so probably not a "good first issue"...?

[ChrisRackauckas]

This is a deep one. I'm going to probably be doing this one myself, but will post what the internals will work like hopefully soon.

[pepijndevos]

In https://github.com/JuliaComputing/CxxrtlCosim we're trying to use a digital simulator for cosimulation with MTK.
This works great with DiscreteSystem, but what would be even more amazing is if promotion from DiscreteSystem to ODESystem works, which would allow full mixed signal simulation.

Inside of ODESystem, promote differencing equations to DiscreteCallbacks

[JTaets]

Since working some more with ModelingToolkit, I would like to express my opinion on this.

I think that D(x) ~ 0 is still best, having discrete variables in prob.p seems wrong.
Theoretically the state means the information required to fully specify the system and that contains the discrete states.
Putting it inside the parameters makes less sense to me, because parameters should be in my opinion constants or functions of time, not functions of other states.

I think the problem of solving ODEs being ~O(n^3) can be lowered for discrete states x_discr in some way.
Maybe only simulate the continuous 'substate' x_cont in the ODE-solvers like dx_cont = f(vcat(x_cont,x_discr),p,t) (but its inplace version) or even dx_cont = f(x_cont,x_discr,p,t) . The only place where the discrete state comes into play is in the callbacks and at the input of f

Don't know how feasible this is and I could be horribly wrong.