Jacobian order

[baggepinnen]

There is currently a potential mismatch between the order in which states and equations are stored, in the following case, the dynamics for the spring deflection and damper velocity are equations 3 and 4, but the corresponding states are numbered 3 and 6. There is no error here, but when the jacobian is calculated, it will not be on the form

x' = Ax

but rather

xp = states(model)[permutation]
xp' = A*xp

where permutation is determined by the equation order.

julia> equations(model)
6-element Vector{Equation}:
 Differential(t)(mass1₊vel(t)) ~ (u(t) - sd₊damper₊c*sd₊damper₊vel(t) - sd₊spring₊k*sd₊spring₊x(t)) / mass1₊m
 Differential(t)(mass2₊vel(t)) ~ (sd₊damper₊c*sd₊damper₊vel(t) + sd₊spring₊k*sd₊spring₊x(t)) / mass2₊m
 sd₊spring₊x(t) ~ mass1₊pos(t) - mass2₊pos(t)
 sd₊damper₊vel(t) ~ mass1₊vel(t) - mass2₊vel(t)
 Differential(t)(mass1₊pos(t)) ~ mass1₊vel(t)
 Differential(t)(mass2₊pos(t)) ~ mass2₊vel(t)

julia> states(model)
6-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 mass1₊vel(t)
 mass2₊vel(t)
 sd₊spring₊x(t)
 mass1₊pos(t)
 mass2₊pos(t)
 sd₊damper₊vel(t)

julia> calculate_jacobian(model)
6×6 Matrix{Num}:
 0   0  (-sd₊spring₊k) / mass1₊m  0   0  (-sd₊damper₊c) / mass1₊m
 0   0     sd₊spring₊k / mass2₊m  0   0    sd₊damper₊c / mass2₊m
 0   0                         0  1  -1          0
 1  -1                         0  0   0          0
 1   0                         0  0   0          0
 0   1                         0  0   0          0

The model above is not minimal, in the sense that some states can be computed from others. The problem goes away after strucutral_simplify

julia> equations(model)
4-element Vector{Equation}:
 Differential(t)(mass1₊vel(t)) ~ (u(t) - sd₊spring₊k*(mass1₊pos(t) - mass2₊pos(t)) - sd₊damper₊c*(mass1₊vel(t) - mass2₊vel(t))) / mass1₊m
 Differential(t)(mass2₊vel(t)) ~ (sd₊spring₊k*(mass1₊pos(t) - mass2₊pos(t)) + sd₊damper₊c*(mass1₊vel(t) - mass2₊vel(t))) / mass2₊m
 Differential(t)(mass1₊pos(t)) ~ mass1₊vel(t)
 Differential(t)(mass2₊pos(t)) ~ mass2₊vel(t)

julia> calculate_jacobian(model)
4×4 Matrix{Num}:
 (-sd₊damper₊c) / mass1₊m    sd₊damper₊c / mass1₊m   (-sd₊spring₊k) / mass1₊m     sd₊spring₊k / mass1₊m
   sd₊damper₊c / mass2₊m   (-sd₊damper₊c) / mass2₊m     sd₊spring₊k / mass2₊m  (-sd₊spring₊k) / mass2₊m
         1                         0                                        0                         0
         0                         1                                        0                         0

Is the ordering after simplification a lucky coincidence, or can the ordering of states and equations always be made the same?

[Keno]

Your request is to have the D(x_i) ~ rhs to have the same index in equations as x_i does in states? For a general system there's no unique association between the equations and the x_i states, so it's probably not possible in general. Why do you need this?

[baggepinnen]

Your request is to have the D(x_i) ~ rhs to have the same index in equations as x_i does in states?

Yes! It's not strictly needed for anything, other than to match the textbook definition of a linear system of ODE's. When the user calls the (yet to be defined) linearize function to get a linear statespace system back

x' = Ax + Bu
y  = Cx + Du

they expect x' and x to have the relation x' = dx/dt.

I realize this can't hold in general, in particular if there are algebraic equations left that can not be simplified. In such a case you could potentially linearize the system to

Ex' = Ax + Bu
 y  = Cx + Du

in which case it's common to place the algebraic equations last, i.e., E = blkdiag(I(n_differential), 0*I(n_algebraic))

[crlaugh]

This practice is followed in Dymola and all other Modelica tools with which I am familiar - the order of D(x) is always the same as that of x which multiplies the Jacobian A, in my experience. Because the standard practice in these tools is to reduce to index 0, these tools always end up with a linear system of ODEs. If we only reduce the index to 1, using a canonical form like

Ex' = Ax + Bu 
y  = Cx + Du

would be helpful.

To add to @baggepinnen's comments, I can think of two reasons to do this:

1. Convention among control engineers/mathematical modelers, and
2. Looking at the structure of A with knowledge of the states can be very helpful in model diagnostics. Since the structure of A tells the model designer which terms are coupled to the derivatives of which states, the model designer only has to remember one order of states in the vector when looking at the details of the model construction. In general, the BLT structure of A brings out the fact that most couplings are fairly local for many systems (with some notable exceptions), so one can quickly ascertain the couplings between terms via the size and makeup of the blocks.

I have used this technique a number of times when looking at large-scale models because missing (or extra) terms in the model are easy to see this way.

[ChrisRackauckas]

It matches states(sys) which is the canonical ordering in the model. Specify the ordering there if you want it different: the ODESystem constructor can take in the array of states in the order you wish for this purpose. The equations cannot give a canonical ordering unless it's an ODE system, since then it's not necessarily defined which equations define which variables in the algebraic parts.