search

SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.39k stars 199 forks source link

Discrete-Time (Sub-)Systems #894

Closed asprionj closed 4 months ago

asprionj commented 3 years ago

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 commented 3 years ago

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: https://github.com/JuliaSymbolics/Symbolics.jl/issues/32 Discrete system: https://github.com/SciML/ModelingToolkit.jl/issues/571 Events: https://github.com/SciML/ModelingToolkit.jl/issues/38 Promotion: https://github.com/SciML/ModelingToolkit.jl/issues/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:

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 commented 3 years ago

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

jonniedie commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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 commented 2 years ago

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 commented 2 years ago

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.

ChrisRackauckas commented 2 years ago

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

The way to implement that well is... to put x_discr into p. In fact, it would work so well that it would even work with the tie-ins for adjoints of callbacks. https://github.com/SciML/DiffEqSensitivity.jl/blob/master/src/callback_tracking.jl#L220-L236 so that will make it reverse the discrete system when you reverse the solve. It's a rather beautiful solution if you look at all of the details.

JTaets commented 2 years ago

I like it then!

I'm a bit scared that using DifferentialEquations.jl and packages based on it will be a little bit more complicated to use due to the extra part in the parameters, which is hard to explain to first time users.

But I guess the inner workings can be hidden by deprecating prob.p and prob.u0 in favor of syntax like

And also remake(prob,p=...) acts on the first part of the parameters.

I think in such case, you have the best of both worlds. Easy implementation and still easily understandable outer workings I'm sure you guys will figure something out :)

ChrisRackauckas commented 2 years ago

But I guess the inner workings can be hidden by deprecating prob.p and prob.u0 in favor of syntax like

Yeah that's the plan. Actually, we already have sol[u] etc., so you shouldn't have to go poking into sol.u and remember array orderings anyways. We should just return the values as you would expect from a symbolic interface, and that gives us freedom to optimize the interns how we please. DifferentialEquations.jl is numerical so you control everything. ModelingToolkit is a modeling framework so it should take control.

ChrisRackauckas commented 4 months ago

See @baggepinnen 's various issues on the newer plans and versions of discrete continuous systems.