Symbolics and ModelingToolkit

General Usage
,
1

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[1] = x, q[2] = θ
params = @parameters m₁ m₂ l g F

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

t = vars[1]
qs = vars[2]

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[2]*l

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

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[1] => 0, 
        qs[2] => 0,
        q̇s[1] => 0,
        q̇s[2] => 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?
thanks in ahead.

2

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

3

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

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[1] = x, q[2] = θ
params = @parameters m₁ m₂ l g F

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

t = vars[1]
qs = vars[2]

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[2]*l

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

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[1] => 0, # sol.u[3] = x
        qs[2] => 0, # sol.u[1] = θ
        q̇s[1] => 0, # sol.u[4] = v
        q̇s[2] => 0] # sol.u[2] = ω


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)
4

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))[2]) ~ var"q(t)[2]ˍt"
 Differential(t)(var"q(t)[2]ˍt") ~ var"q(t)[2]ˍtt"
 Differential(t)((q(t))[1]) ~ var"q(t)[1]ˍt"
 Differential(t)(var"q(t)[1]ˍt") ~ var"q(t)[1]ˍtt"
 0 ~ (l*m₂*var"q(t)[1]ˍtt"*cos((q(t))[2]) + g*l*m₂*sin((q(t))[2])) / (-m₂*(l^2)) - var"q(t)[2]ˍ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[2])) => 0 to do this.

X_0 =  [qs[1] => 0, 
        qs[2] => 0,
        q̇s[1] => 0,
        q̇s[2] => 0,
        dt(dt(q[2])) => 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)

image
image1106×800 176 KB

5

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[1] = x, q[2] = θ
params = @parameters m₁ m₂ l g F

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

t = vars[1]
qs = vars[2]

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[2]*l

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

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[1])), dt(dt(qs[2]))])

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[1] => 0,
        qs[2] =>  π/6, 
        q̇s[1] => 0, 
        q̇s[2] => 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))[1])) ~ (l*m₂*Differential(t)(Differential(t)((q(t))[2]))*cos((q(t))[2]) - F*cos(t) - l*m₂*(Differential(t)((q(t))[2])^2)*sin((q(t))[2])) / (-m₁ - m₂)
 Differential(t)(Differential(t)((q(t))[2])) ~ (l*m₂*Differential(t)(Differential(t)((q(t))[1]))*cos((q(t))[2]) + g*l*m₂*sin((q(t))[2])) / (-m₂*(l^2))

which find the second derivative of each state as function of the second derivative of the other state
to over come this I add additional step :

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

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))[1])) ~ (-F*cos(t) - l*m₂*(Differential(t)((q(t))[2])^2)*sin((q(t))[2])) / (-m₁ - m₂) + (l*m₂*((g*sin((q(t))[2])) / l + (((-F*cos(t) - l*m₂*(Differential(t)((q(t))[2])^2)*sin((q(t))[2])) / (-m₁ - m₂))*cos((q(t))[2])) / l)*cos((q(t))[2])) / (((-m₂*(cos((q(t))[2])^2)) / (-m₁ - m₂) - 1)*(-m₁ - m₂))
 Differential(t)(Differential(t)((q(t))[2])) ~ ((g*sin((q(t))[2])) / l + (((-F*cos(t) - l*m₂*(Differential(t)((q(t))[2])^2)*sin((q(t))[2])) / (-m₁ - m₂))*cos((q(t))[2])) / l) / ((-m₂*(cos((q(t))[2])^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.

plot_2
plot_21024×683 103 KB

thank you all for the help.

1 Like