Difference between solutions obtained from ODE model and ModelingToolkit model

Specific Domains Modelling & Simulations
1

Hi all, I am starting using OrdinaryDiffEq and ModelingToolkit and try to learn how to use these packages properly. To this end, I took the double pendulum as an example.

For the OrdinaryDiffEq model, I have implemented the following code :

using OrdinaryDiffEq, GLMakie

# First pendulum parameters
const m₁, l₁, L₁, I₁ = 0.1375, 0.3, 0.1536, 0.0044

# Second pendulum parameters
const m₂, l₂, L₂, I₂ = 0.1375, 0.3, 0.1536, 0.0044

# Gravity constant
const g = 9.81

# Equations of motion 
function double_pendulum!(dq, q, p, t)
    M24 = m₂*l₁*L₂*cos(q[3])
    M22 = I₁ + m₁*l₁^2 + M24
    M42 = I₂ + M24
    M44 = I₂

    M = [1.  0. 0. 0.;
        0. M22 0. M24;
        0. 0. 1. 0.;
        0. M42 0. M44]

    dq .= M\[
        q[2],
        m₂*l₁*L₂*(q[2] + q[4])^2*sin(q[3]) - (m₁*L₁ + m₂*l₁)*g*sin(q[1]),
        q[4],
        -m₂*l₁*L₂*q[2]^2*sin(q[3]) - m₂*L₂*g*sin(q[1] + q[3])
    ]
end

# System solution
q0 = [deg2rad(30.), 0., deg2rad(120.), 0.]
tspan = (0., 10.)

prob = ODEProblem(pendule_double!, q0, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)
teval = LinRange(tspan[1], tspan[2], 2000) 

# Results display
fig1 = Figure()
ax1 = Axis(fig1[1, 1], xlabel = "Time(s)", ylabel = "Angle (°)")
lines!(ax1, teval, rad2deg.(sol(teval)[1, :]), color = :blue, label = "alpha", linewidth = 2)
lines!(ax1, teval, rad2deg.(sol(teval)[3, :]), color = :red, label = "beta", linestyle = :dash, linewidth = 2)
xlims!(ax1, teval[1], teval[end])
axislegend(position = :lb)
display(GLMakie.Screen(), fig1)

The results obtained are the same as those I got with Simscape or by implementing it in Matlab or Python.

In another attempt, I tried to implement this problem using ModelingToolkit.The only difference with the previous code are given below :

using OrdinaryDiffEq, ModelingToolkit, GLMakie
using ModelingToolkit: t_nounits as t, D_nounits as D

# First pendulum parameters
const m₁, l₁, L₁, I₁ = 0.1375, 0.3, 0.1536, 0.0044

# Second pendulum parameters
const m₂, l₂, L₂, I₂ = 0.1375, 0.3, 0.1536, 0.0044

# Gravity constant
const g = 9.81

# Equations of motion 
@variables q₁(t), q₂(t), q₃(t), q₄(t)
eqs = [
    D(q₁) ~ q₂,
    (I₁ + m₂*l₁^2 + m₂*l₁*L₂*cos(q₃))*D(q₂) + m₂*l₁*L₂*cos(q₃)*D(q₄) ~ m₂*l₁*L₂*(q₂ + q₄)^2*sin(q₃) - (m₁*L₁ + m₂*l₁)*g*sin(q₁),
    D(q₃) ~ q₄,
    (I₂ + m₂*l₁*L₂*cos(q₃))*D(q₂) + I₂*D(q₄) ~ -m₂*l₁*L₂*q₂^2*sin(q₃) - m₂*L₂*g*sin(q₁ + q₃)
]

@named sys = ODESystem(eqs, t)
sys = structural_simplify(sys)

# System solution
u0 = [
    q₁ => deg2rad(30.),
    q₂ => 0.,
    q₃ => deg2rad(120.),
    q₄ => 0.
]
tspan = (0., 10.)

prob = ODEProblem(sys, u0, tspan, jac = true)
sol = solve(prob, Rodas5()) # For dealing with mass-matrix

Despite the solver provides a solution, the latter is far from that obtained with OrdinaryDiffEq.

What am I missing here ? Indeed, I think I have implemented something the wrong way.

Thank you for your help.

2

What exactly is the issue? Because I plotted them side by side and they look very similar. Can you describe what problem you’re seeing?

plots
plots600×400 56.5 KB

1 Like
3

After reading your answer, I figured out that there is actually no issue in the results. I didn’t pay attention to the way the variables are ordered after sys = structural_simplify(sys). I naively believed, that sol[1, :] = q1, sol[2, :] = q2, which is not the case, because in my cases I have :

julia> unknowns(sys)
5-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 q₁(t)
 q₃(t)
 q₂(t)
 q₄(t)
 q₂ˍt(t)

So in my case, to properly display the solution, I have to write :

fig1 = Figure()
ax1 = Axis(fig1[1, 1],
           xlabel = "Time (s)",
           ylabel = "Angle (°)")
lines!(ax1, teval, rad2deg.(sol(teval)[1, :]), color = :blue, label = "alpha", linewidth = 2)
lines!(ax1, teval, rad2deg.(sol(teval)[2, :]), color = :red, label = "beta", linestyle = :dash, linewidth = 2)
xlims!(ax1, teval[1], teval[end])
axislegend(position = :lb)
display(GLMakie.Screen(), fig1)

Or even better code (read the docs) :

fig1 = Figure()
ax1 = Axis(fig1[1, 1],
           xlabel = "Time (s)",
           ylabel = "Angle (°)")
lines!(ax1, teval, rad2deg.(sol(teval)[q₁]), color = :blue, label = "alpha", linewidth = 2)
lines!(ax1, teval, rad2deg.(sol(teval)[q₃]), color = :red, label = "beta", linestyle = :dash, linewidth = 2)
xlims!(ax1, teval[1], teval[end])
axislegend(position = :lb)
display(GLMakie.Screen(), fig1)

Thank you @ChrisRackauckas !

4

Yes use the symbolic indexing interface to be robust. Here the structural simplification is adding a dummy derivative to enforce the condtion, so it should be a more correct formulation for long term simulations where the other form can have numerical drifting