Force MTK to generate explicit ODEs?

Specific Domains Modelling & Simulations
, , , ,
1

Hi!

I am using ModelingToolkit.jl to model a physical system. The compiler generates an index-1 DAE, which works great for simulation.

However, I need to use a standard Kalman Filter, which requires a minimal state-space representation (\dot{x} = Ax + Bu) with a full-rank observability matrix.

The Problem:

When I use ModelingToolkit.linearize(sys, inputs, outputs), the resulting state-space matrices include the dummy (algebraic) states.

  • Total states: 5 (3 true states + 2 dummy/algebraic states).
  • Observability Matrix Rank: 3.
  • Result: The system is unobservable and the Kalman covariance update fails.

My Question:

Is there a way to force MTK to perform full substitution/elimination of these algebraic variables to generate a purely explicit ODE system before or during linearization? I essentially need to convert the Mass Matrix DAE into a standard ODE so I can get a minimal A matrix.

1 Like
2

You can use model-order index reduction: Run only the algebraic part of the solver - #15 by baggepinnen

1 Like
3

See Home · ControlSystemsMTK Documentation for how to get a normal statespace system

1 Like
4

See also GitHub - baggepinnen/LowLevelParticleFiltersMTK.jl: An interface for state estimation using LowLevelParticleFilters on ModelingToolkit models which might be of interest

1 Like
5

Thanks for the link!
I got it to work, but, interestingly, not with named_ss(), though i suspect the issue is linearize() which is called by named_ss.

For reference, the DAE produced by mtk consists of 3 differential variables (position s, current i, velocity v) and 2 algebraic variables (force f, acceleration a).

The equations(sys) are:

 Differential(t)(axisâ‚Šprismaticâ‚Šs(t)) ~ axisâ‚Šprismaticâ‚Šv(t)
 Differential(t)(axisâ‚Šmotorâ‚Šinductorâ‚Ši(t)) ~ axisâ‚Šmotorâ‚Šinductorâ‚Šv(t) / axisâ‚Šmotorâ‚Šinductorâ‚ŠL
 Differential(t)(axisâ‚Šprismaticâ‚Šv(t)) ~ axisâ‚Šfrictionâ‚Šflange_bâ‚ŠsËŤtt(t) # <- This is an algebraic state
 0 ~ axisâ‚Šprismaticâ‚Šf(t) + axisâ‚Štransâ‚Šflange_transâ‚Šf(t)
 0 ~ axisâ‚Šbody_movâ‚Šframe_aâ‚ŠxËŤtt(t) - axisâ‚Šbody_movâ‚Šr1ËŤtt(t)

Using named_ss produced the same ss model as linearize with all 5 states.

Here is the A matrix produced by linearize:

5Ă—5 Matrix{Float64}:
  0.0    0.0         1.0      0.0   0.0
 -0.0  -20.0      -251.327   -0.0  -0.0
  0.0    0.0         0.0      0.0   1.0
  0.0  -24.5417   -308.4      0.0   0.0
  0.0   -4.90835   -61.6801   0.0   0.0

Notice that the rows 2, 4, 5 (i, f, a) are all linearly dependent.
That feels like a problem?

I tried bypassing linearize by manually building the Descriptor System, using calculate_massmatrix and ForwardDiff.jacobian on the compiled ODEProblem.

# 1. Extract Mass Matrix (E) from the working ODEProblem
E = prob.f.mass_matrix

# 2. Compute A numerically using ForwardDiff on the ODE function
f_func = (u) -> (out=similar(u); prob.f(out, u, prob.p, t); out)
A_numeric = ForwardDiff.jacobian(f_func, prob.u0)

# 3. Convert manually
sys_dss = dss(A_numeric, E, B, C, D)
sys_final = dss2ss(sys_dss, fast=false)

The A matrix now looks like this:

5Ă—5 Matrix{Float64}:
  0.0    0.0         1.0     0.0   0.0
 -0.0  -20.0      -251.327  -0.0  -0.0
  0.0    0.0         0.0     0.0   1.0
  0.0    1.25664     0.0    -1.0  -0.120409
  0.0    0.0         0.0    -0.2   1.0

sys_final now is the correct ss model with 3 states.

Am i missing something or does linearize actually fail to linearize my model correctly?

6

Would you mind posting the complete model and your call to linearize, I’ll look into it in the morning?

7

Sure, i can do that.
Might take me some time though, i need to extract the components first so you wont have to deal with a hundred files.

8

You can also try calling minreal on the system returned by named_ss to see if you have gotten a model with pole-zero cancelations

9

Sorry for the delay.
I created a test model and a test script, but im not allowed to upload files or post links yet.
Is there another way to post code?
But maybe it won’t be necessary, because the issue is probably that I don’t understand named_ss properly.

I made a mistake when testing linearize a few days back.
I beliefe linearize works just fine (simulation of the linearized model is correct), but named_ss still doesn’t remove the dummy states.

My question now is, what is named_ss trying to do?
I demonstrated in a previous reply that dss2ss is capable of producing a reduced state space representation of my model after manually building the descriptor system.
As i understand it, named_ss also uses dss2ss, but retrieves the linearized model using linearize, which (always?) returns an explicit state space model containing the dummy states, instead of a descriptor system.
So does named_ss try to convert the ss model into descriptor form before calling dss2ss?

To answer your last question: minreal did work and produced an observable system. :+1:

10

Yeah I suspected that would work. linearize converts an index 1 DAE into an ODE, but neither linearize nor named_ss do any reductions, dss2ss does. The only thing named_ss does over linearize is that it can handle situations in which MTK fails to recognize a proper system, producing equations that appear to require the derivative of inputs, named_ss performs a transformation to proper form in this case. In addition to this, named_ss wraps the result of linearization in a statespace object.

minreal computes a minimal realization, that is, it removes modes that are unobservable or uncontrollable, cancelling poles and zeros. dss2ss also did this for you. Calling named_ss |> minreal, or named_ss |> sminreal is thus recommended when you have a DAE. Note, the result is not incorrect with the additional variables, it’s just a larger representation than you need to represent the transfer function from the input to the output.

It’s no longer needed now that we know that the issue is :slight_smile:

1 Like