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.