MTK equation explorations: mechanical example

mechanical-example.md

I'm trying to implement this Modelica example from Michael Tiller's Modelica by Design. I don't understand the rules for writing DAE equations in MTK.

Here's a version in semi-explicit form with the derivative terms alone on the left-hand side:

using ModelingToolkit, OrdinaryDiffEq
const D = Differential(t)
@variables t

@variables phi1(t) phi2(t)=1.0 omega1(t) omega2(t)
J1 = 0.4; J2 = 1.0; c1 = 11.0; c2 = 5.0; d1 = 0.2; d2 = 1.0
eq = [
    ## Equations for inertia 1
    D(phi1) ~ omega1
    D(omega1) ~ (c1*(phi2 - phi1) + d1*(omega2 - omega1)) / J1
    ## Equations for inertia 2
    D(phi2) ~ omega2
    D(omega2) ~ (c1*(phi1 - phi2) + d1*(omega1 - omega2) - c2*phi2 - d2*omega2) / J2
]
@named sys = ODESystem(eq)
sys = structural_simplify(sys)
prob = ODAEProblem(sys, [0.0, 0.0, 1.0, 0.0], (0.0, 5.0))
sol = solve(prob, Tsit5())
#plot(sol)

If the J variables are moved to the LHS, the model still compiles and solves, but the results are wrong.

@variables phi1(t) phi2(t)=1.0 omega1(t) omega2(t)
J1 = 0.4; J2 = 1.0; c1 = 11.0; c2 = 5.0; d1 = 0.2; d2 = 1.0
eq = [
    ## Equations for inertia 1
    D(phi1) ~ omega1
    J1 * D(omega1) ~ c1*(phi2 - phi1) + d1*(omega2 - omega1)
    ## Equations for inertia 2
    D(phi2) ~ omega2
    J2 * D(omega2) ~ c1*(phi1 - phi2) + d1*omega1 - d1*omega2 - c2*phi2 - d2*omega2
]
@named sys = ODESystem(eq)
sys = structural_simplify(sys)
prob = ODAEProblem(sys, [0.0, 1.0, 0.0, 0.0], (0.0, 5.0))
sol = solve(prob, Tsit5())

Here is the result. Note that the first and last states (phi1 and omega in this case) are zero in the solution.

julia> sol = solve(prob, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 33-element Vector{Float64}:
 0.0
 6.243756243756244e-5
 0.0006868131868131868
 0.00693056943056943
 0.0348374538259648
 0.0893616793385086
 0.1653307171340363
 0.26240965596179794
 ⋮
 3.5528948115882653
 3.7836471920118484
 3.987921039780141
 4.210449676082234
 4.428818087102506
 4.649549476064513
 4.869944997278494
 5.0
u: 33-element Vector{Vector{Float64}}:
 [0.0, 1.0, 0.0, 0.0]
 [0.0, 0.9999999688131854, -0.0009989635644382343, 0.0]
 [0.0, 0.9999962273400609, -0.01098446997686435, 0.0]
 [0.0, 0.99961682523775, -0.11041512769799762, 0.0]
 [0.0, 0.990440150823186, -0.5441439014764878, 0.0]
 [0.0, 0.9389885280051328, -1.3271124111832413, 0.0]
 [0.0, 0.8023591682935984, -2.228408563735211, 0.0]
 [0.0, 0.5457558367434211, -2.976889984759295, 0.0]
 ⋮
 [0.0, 0.028072734697371024, -0.4779846928236069, 0.0]
 [0.0, -0.06537025888426529, -0.2828116395253121, 0.0]
 [0.0, -0.0920292458364988, 0.023670804272775167, 0.0]
 [0.0, -0.05667651977173833, 0.261888643671073, 0.0]
 [0.0, 0.006128488089785196, 0.2756482646292357, 0.0]
 [0.0, 0.050855520312849734, 0.1104192952305579, 0.0]
 [0.0, 0.0525687615449217, -0.08695175519248523, 0.0]
 [0.0, 0.03595629036784663, -0.1607920150668361, 0.0]

It seems like a model with poorly formed equations should show a warning (or even fail). I couldn't find anything in the documentation on allowed forms for DAE equations.