New calculate_statespace function for AbstractODESystem

[cadojo]

As mentioned in #1059, we can now generate the A and B matrices of linear state space models for any AbstractODESystem because we have the new calculate_control_jacobian function.

Users can generate_jacobian and generate_control_jacobian for fast functions which return the system's A and B matrices. Still, it would be nice to have a convenience function which simply returns the linearizations for any AbstractODESystem. Users could call calculate_statespace(sys::AbstractODESystem, subs::Dict) to return the linearization evaluated at the points in the Dict subs.

That's what this PR adds. One additional test has been added, which passes.

The reason I did not use generate_jacobian and generate_control_jacobian in calculate_statespace is because I did not want calculate_statespace to have to compile functions with every call, and it doesn't seem right to cache RuntimeGeneratedFunctions in every concrete AbstractODESystem type.

Anyone have any opinions / thoughts on whether this should be included? I'd find it useful for PolynomialGTM.jl.

[baggepinnen]

Would it make sense to include functionality for the user to specify outputs also? So that one can obtain a system on the form

x' = Ax + Bu
y = Cx + Du

The controls already indicate what u is so that B can be obtained, perhaps measured outputs could be indicated similarly?

[cadojo]

Would it make sense to include functionality for the user to specify outputs also?

I'd be happy implementing this, though I'm a controls engineer, and I worry about us GN&C folks overtaking a modeling package which is designed for a different purpose (scientific machine learning).

Right now, I don't see any easy way to specify outputs in an AbstractODESystem, I think we'd need to change eqs in ODESystem (and similar systems) to state_eqs, and then we'd need to add output_eqs and outputs to ODESystem. Like I said, I'm happy to do so, but what are other peoples thoughts?

[ChrisRackauckas]

The reason I did not use generate_jacobian and generate_control_jacobian in calculate_statespace is because I did not want calculate_statespace to have to compile functions with every call, and it doesn't seem right to cache RuntimeGeneratedFunctions in every concrete AbstractODESystem type.

Makes sense to me.

Would it make sense to include functionality for the user to specify outputs also?

I think so. It should somehow carry over the observed equations from the original system, right?

[cadojo]

I think so. It should somehow carry over the observed equations from the original system, right?

I hadn't played around with observed before in MTK – I just did. Here's the familiar spring-mass-damper example.

# On MTK's master branch as of this post...

using ModelingToolkit

sys = let
    @parameters t f k d
    @variables x(t) ẋ(t) y(t)
    δ = Differential(t)
           
    eqs = [
        δ(x) ~ ẋ,
        δ(ẋ) ~ f - k*x - d*ẋ
    ]
           
    @named HarmonicOscillator = ODESystem(eqs, t, [x, ẋ], [f, k, d]; controls=[f], observed=[y ~ x + ẋ])
end |> structural_simplify

We could have also left observed blank in the ODESystem constructor, and simply added y ~ x + ẋ to eqs. The structural_simplify function would have moved y from equations to observed.

I agree it makes sense that outputs are somehow associated with the existing observed functionality. But then, should all observed variables be assumed outputs? I think not, but then if we linearize the system, we'd need to substitute the symbolic observed equations with numerical values, and those observed equations fall out of the linearization. Maybe that's okay? I don't feel I'm experienced enough to say, after looking at this for a couple minutes.

Also, if outputs are a subset of observed, then the usage feels a bit clunky. You'd specify every control, and every output twice. If we do add outputs, should I go ahead and separate parameters and controls too? I think that might be the right call anyway.

[ChrisRackauckas]

I agree it makes sense that outputs are somehow associated with the existing observed functionality. But then, should all observed variables be assumed outputs? I think not, but then if we linearize the system, we'd need to substitute the symbolic observed equations with numerical values, and those observed equations fall out of the linearization. Maybe that's okay? I don't feel I'm experienced enough to say, after looking at this for a couple minutes.

The difference though is that the observed equations, if left fully intact, will not be y = Cx + Du but something nonlinear. @baggepinnen do those need to be linear? Would it make sense to linearize them automatically?

[ChrisRackauckas]

I'd be happy implementing this, though I'm a controls engineer, and I worry about us GN&C folks overtaking a modeling package which is designed for a different purpose (scientific machine learning).

Don't worry about it. In fact, the whole purpose of ModelingToolkit is to be a general enough symbolic specification of these systems that all domains can feed into it. I think controls functionality is very necessary for a lot of domains, even if in some cases it's called "input functions" or "dosing". So observed variables, the ability to specify controls, transform this into optimization problems and other things: this is all very necessary to make first class in the ODESystem.

[cadojo]

The difference though is that the observed equations, if left fully intact, will not be y = Cx + Du but something nonlinear.

Might be misunderstanding the question here, but I figured the linearization would work the same way the Ax + Bu linearization works. The A and B matrices are found by taking the partial derivatives of the state equations with respect to states x, and controls u.

To find C and D, we'd linearize the observed equations with respect to the outputs y, and the controls u.

So it's alright for the observed equations to be nonlinear in an ODESystem the same way the state equations are allowed to be nonlinear. Right?

Edit. Per my earlier comment...

I agree it makes sense that outputs are somehow associated with the existing observed functionality. But then, should all observed variables be assumed outputs? I think not, but then if we linearize the system, we'd need to substitute the symbolic observed equations with numerical values, and those observed equations fall out of the linearization. Maybe that's okay? I don't feel I'm experienced enough to say, after looking at this for a couple minutes.

After writing out the logic myself in this comment, I think it is fine to specify outputs, and have them live alongside other observed quantities.

[ChrisRackauckas]

yeah, I guess linearize the observed equations and that will give that form.

[cadojo]

It would be fine to have non-output quantities in the left-hand side of the observed equations in the same way it's fine to have non-control quantities in the parameters, right?

Edit. I guess we need to decide if the number of observed quantities is equivalent to the number of outputs, or if the output equations are a subset of the observed equations.

[baggepinnen]

As I understand observed it's a list of expressions that can be derived from other states without integration? This can be different from a set of measured outputs in general. A measurement equation on a general form is any function of both states and inputs y = g(x,u), (possibly also time g(x,u,t) ?) with a linearized form y = Cx + Du. For calculate_statespace it makes sense to return the form y = Cx + Du, but before such a linearization is requested I think it may be valuable to keep any original expression for y provided by the user.

---

As a side, perhaps calculate_statespace should really be called linearize_dynamics or something like that, since a statespace system can just as well be nonlinear

x' = f(x,u,t)
y  = g(x,u,t)

---

Unrelated to the current discussion, but touching on

the whole purpose of ModelingToolkit is to be a general enough symbolic specification of these systems that all domains can feed into it.

It's common in controls to also have a noise/disturbance specification. Such a specification can take the form

x' = Ax + Bu + Ew where w ~ N(0, I)
y = Cx + Du + Fe where e ~ N(0, I)

This form is used to specify how disturbances enter the system. Commonly, measurement noise, periodic noise (think 50/60Hz measurement disturbances), or low-frequency model errors like Coulomb friction acting in the same place as the input u.
Disturbances are often treated as a second set of inputs, but they are often not known. The most typical references on disturbance models would be

• Robust and Optimal Control, Zhou , Doyle, Glover
• Introduction to Stochastic Control Theory, Karl J Astrom

Maybe this kind of disturbance specification can be left for later, but it would be good to keep in the back of the mind while designing interfaces etc. In typical control-systems design ala matlab or ControlSystems.jl, noise models are often added to a system model, some simple examples are here (mostly undocumented unfortunately)
https://github.com/JuliaControl/RobustAndOptimalControl.jl/blob/master/src/model_augmentation.jl
after which an observer (e.g., a kalman filter) is constructed to estimate the states of both the system and the noise model.

[ChrisRackauckas]

It's common in controls to also have a noise/disturbance specification. Such a specification can take the form

That's an SDESystem, which is an AbstractODESystem which all of this should apply to.

[YingboMa]

Closing in favor of the #1111 approach