# Symbolics and ModelingToolkit

hi,
i’m trying to simulate a cart pole using the code example taken from this topic.
my code looks like this:

``````using Symbolics, Plots, DifferentialEquations, ModelingToolkit

n = 2 # number of general coordinates

vars = @variables t (q(t))[1:n] # q = x, q = θ
params = @parameters m₁ m₂ l g F

Q = [F*cos(t), 0]# general forces

t = vars
qs = vars

dt = Differential(t)
dqs = [Differential(q) for q in qs]
dq̇s = [Differential(dt(q)) for q in qs]

q̇s = [dt(q) for q in qs]

v = q̇s*l

T = 0.5*(m₁+m₂)*q̇s^2 + 0.5*m₂*v^2 + m₂*l*q̇s*q̇s*cos(qs)
U = -m₂*g*l*cos(qs)

L = T - U

Eᵢ = Vector{Num}()
Eᵢᵢ = Vector{Num}()

for (ind, (q, dq̇, dq)) in enumerate(zip(qs, dq̇s, dqs))
push!(Eᵢ, expand_derivatives(dt(expand_derivatives(dq̇(L)))) - expand_derivatives(dq(L)))
push!(Eᵢᵢ, Symbolics.solve_for(Eᵢ[ind] ~ Q[ind], dt(dt(q))))
end

system = [dt(dt(q)) ~ Eᵢᵢ[ind] for (ind, q) in enumerate(qs)]

@named HO = ODESystem(system, t, qs, params)

HO = structural_simplify(HO) # HERE IS THE PROBLEM : Model HO with 5 equations

X_0 =  [qs => 0,
qs => 0,
q̇s => 0,
q̇s => 0]

M = [   m₁ => 30, # kg
m₂ => 10, # kg
l => 30, # m
g => 9.81, # m/s^2
F => 0] # N

tspan = (0.0, 100.0)

prob = ODEProblem(HO, X_0, tspan, M, abstol = 1e-8, reltol = 1e-8)
sol = solve(prob)
``````

my problem is this: when I call the `structural_simplify(HO)` I get 5 equations which souled be 4.
what am I missing and what did I do wrong?

Just a quick look. With respect to the earlier solved problem, you should have some “expand_derivatives” terms in T. Can you check that?

hi,
thanks for the suggestion but unfortunately, after trying to call `expand_derivatives` in any possible place and combination I still get 5 equations .

my updated code with all the newly placed `expand_derivatives`

``````using Symbolics, Plots, DifferentialEquations, ModelingToolkit

n = 2 # number of general coordinates

vars = @variables t (q(t))[1:n] # q = x, q = θ
params = @parameters m₁ m₂ l g F

Q = [F*cos(t), 0]# general forces

t = vars
qs = vars

dt = Differential(t)
dqs = [Differential(q) for q in qs]
dq̇s = [Differential(expand_derivatives(dt(q))) for q in qs]

q̇s = [expand_derivatives(dt(q)) for q in qs]

v = q̇s*l

T = expand_derivatives(0.5*(m₁+m₂)*q̇s^2 + 0.5*m₂*v^2 + m₂*l*q̇s*q̇s*cos(qs))
U = -m₂*g*l*cos(qs)

L = T - U

Eᵢ = Vector{Num}()
Eᵢᵢ = Vector{Num}()

for (ind, (q, dq̇, dq)) in enumerate(zip(qs, dq̇s, dqs))
push!(Eᵢ, expand_derivatives(dt(expand_derivatives(dq̇(L)))) - expand_derivatives(dq(L)))
push!(Eᵢᵢ, Symbolics.solve_for(Eᵢ[ind] ~ Q[ind], expand_derivatives(dt(expand_derivatives(dt(q))))))
end

system = [expand_derivatives(dt(expand_derivatives(dt(q)))) ~ Eᵢᵢ[ind] for (ind, q) in enumerate(qs)]

@named HO = ODESystem(system, t, qs, params)

HO = structural_simplify(HO)

X_0 =  [qs => 0, # sol.u = x
qs => 0, # sol.u = θ
q̇s => 0, # sol.u = v
q̇s => 0] # sol.u = ω

M = [   m₁ => 30, # kg
m₂ => 10, # kg
l => 30, # m
g => 9.81, # m/s^2
F => 0] # N

tspan = (0.0, 100.0)

prob = ODEProblem(HO, X_0, tspan, M, abstol = 1e-8, reltol = 1e-8)
sol = solve(prob)
``````

You get a system with 4 differential equations and one algebraic equation, so you actually have the number of differential states that you expect.

``````julia> equations(HO)
5-element Vector{Equation}:
Differential(t)((q(t))) ~ var"q(t)ˍt"
Differential(t)(var"q(t)ˍt") ~ var"q(t)ˍtt"
Differential(t)((q(t))) ~ var"q(t)ˍt"
Differential(t)(var"q(t)ˍt") ~ var"q(t)ˍtt"
0 ~ (l*m₂*var"q(t)ˍtt"*cos((q(t))) + g*l*m₂*sin((q(t)))) / (-m₂*(l^2)) - var"q(t)ˍtt"
``````

To simulate this system, you unfortunately have to provide an initial condition also for one dummy derivative, the code below uses `dt(dt(q)) => 0` to do this.

``````X_0 =  [qs => 0,
qs => 0,
q̇s => 0,
q̇s => 0,
dt(dt(q)) => 0]

M = [   m₁ => 30, # kg
m₂ => 10, # kg
l => 30, # m
g => 9.81, # m/s^2
F => 1] # N

tspan = (0.0, 100.0)

prob = ODEProblem(HO, X_0, tspan, M, abstol = 1e-8, reltol = 1e-8)
sol = solve(prob, Rodas5P())
plot(sol)
``````

thank you for the insight to understand this.
after doing some more exploration I found the following code to over come the additional equation.

``````using Symbolics, Plots, DifferentialEquations, ModelingToolkit

n = 2 # number of general coordinates

vars = @variables t (q(t))[1:n] # q = x, q = θ
params = @parameters m₁ m₂ l g F

Q = [F*cos(t), 0]# general forces

t = vars
qs = vars

dt = Differential(t)
dqs = [Differential(q) for q in qs]
dq̇s = [Differential(dt(q)) for q in qs]

q̇s = [dt(q) for q in qs]

v = q̇s*l

T = 0.5*(m₁+m₂)*q̇s^2 + 0.5*m₂*v^2 + m₂*l*q̇s*q̇s*cos(qs)
U = -m₂*g*l*cos(qs)

L = T - U

Eᵢ = Vector{Num}()
Eᵢᵢ = Vector{Num}()

for (ind, (q, dq̇, dq)) in enumerate(zip(qs, dq̇s, dqs))
push!(Eᵢ, expand_derivatives(dt(expand_derivatives(dq̇(L)))) - expand_derivatives(dq(L)))
push!(Eᵢᵢ, Symbolics.solve_for(Eᵢ[ind] ~ Q[ind], dt(dt(q))))
end

system = [dt(dt(q)) ~ E for (E, q) in zip(Eᵢᵢ,qs)]

tmp = Symbolics.solve_for(system, [dt(dt(qs)), dt(dt(qs))])

system = [dt(dt(q)) ~ e for (e, q) in zip(tmp,qs)]

@named HO = ODESystem(system, t, qs, params)

HO = structural_simplify(HO)

X_0 =  [qs => 0,
qs =>  π/6,
q̇s => 0,
q̇s => 0]

M = [   m₁ => 30, # kg
m₂ => 10, # kg
l => 30, # m
g => 9.81, # m/s^2
F => 0] # N

tspan = (0.0, 100.0)

prob = ODEProblem(HO, X_0, tspan, M, abstol = 1e-8, reltol = 1e-8)
sol = solve(prob)
``````

As I understand it `system = [dt(dt(q)) ~ E for (E, q) in zip(Eᵢᵢ,qs)]` produce the following equations:

``````2-element Vector{Equation}:
Differential(t)(Differential(t)((q(t)))) ~ (l*m₂*Differential(t)(Differential(t)((q(t))))*cos((q(t))) - F*cos(t) - l*m₂*(Differential(t)((q(t)))^2)*sin((q(t)))) / (-m₁ - m₂)
Differential(t)(Differential(t)((q(t)))) ~ (l*m₂*Differential(t)(Differential(t)((q(t))))*cos((q(t))) + g*l*m₂*sin((q(t)))) / (-m₂*(l^2))
``````

which find the second derivative of each state as function of the second derivative of the other state

``````tmp = Symbolics.solve_for(system, [dt(dt(qs)), dt(dt(qs))])

system = [dt(dt(q)) ~ e for (e, q) in zip(tmp,qs)]
``````

This step created 2 equations in which the second derivative of each state is a function of lower order derivates of all the other states, which is the “normal” way to write those set of equations.

``````2-element Vector{Equation}:
Differential(t)(Differential(t)((q(t)))) ~ (-F*cos(t) - l*m₂*(Differential(t)((q(t)))^2)*sin((q(t)))) / (-m₁ - m₂) + (l*m₂*((g*sin((q(t)))) / l + (((-F*cos(t) - l*m₂*(Differential(t)((q(t)))^2)*sin((q(t)))) / (-m₁ - m₂))*cos((q(t)))) / l)*cos((q(t)))) / (((-m₂*(cos((q(t)))^2)) / (-m₁ - m₂) - 1)*(-m₁ - m₂))
Differential(t)(Differential(t)((q(t)))) ~ ((g*sin((q(t)))) / l + (((-F*cos(t) - l*m₂*(Differential(t)((q(t)))^2)*sin((q(t)))) / (-m₁ - m₂))*cos((q(t)))) / l) / ((-m₂*(cos((q(t)))^2)) / (-m₁ - m₂) - 1)
``````

(A lot of text but the important thing is that there aren’t any second derivative at the right hand side of each equation).

Now when I call

``````HO = structural_simplify(HO)
``````

I get 4 equations and the following solution.

thank you all for the help.

1 Like