# Symbolics and ModelingToolkit

**URL:** <https://discourse.julialang.org/t/symbolics-and-modelingtoolkit/93010>\
**Category:** General Usage\
**Tags:** modelingtoolkit, symbolics\
**Created:** [January 16, 2023, 5:53am UTC](https://discourse.julialang.org/t/symbolics-and-modelingtoolkit/93010 "2023-01-16T05:53:13Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Daniel\_Lewis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_lewis/32/31250_2.png) [@Daniel\_Lewis](https://discourse.julialang.org/u/Daniel_Lewis)\
**Post date:** [January 16, 2023, 5:53am UTC](https://discourse.julialang.org/t/symbolics-and-modelingtoolkit/93010/1 "2023-01-16T05:53:14Z")

</div>

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

```plaintext
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.

---

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [January 16, 2023, 10:40am UTC](https://discourse.julialang.org/t/symbolics-and-modelingtoolkit/93010/2 "2023-01-16T10:40:04Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Daniel\_Lewis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_lewis/32/31250_2.png) [@Daniel\_Lewis](https://discourse.julialang.org/u/Daniel_Lewis)\
**Post date:** [January 17, 2023, 6:08am UTC](https://discourse.julialang.org/t/symbolics-and-modelingtoolkit/93010/3 "2023-01-17T06:08:09Z")

</div>

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`

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 17, 2023, 6:31am UTC](https://discourse.julialang.org/t/symbolics-and-modelingtoolkit/93010/4 "2023-01-17T06:31:10Z")

</div>

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
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.

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/f/6/f6cb481bf65e29c9b6f748128d61438261ec9502.png)

---

<div class="post-metadata">

**Author:** ![Daniel\_Lewis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_lewis/32/31250_2.png) [@Daniel\_Lewis](https://discourse.julialang.org/u/Daniel_Lewis)\
**Post date:** [January 17, 2023, 6:59am UTC](https://discourse.julialang.org/t/symbolics-and-modelingtoolkit/93010/5 "2023-01-17T06:59:13Z")

</div>

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

```julia
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:

```julia
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 :

```julia
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.

```julia
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

```julia
HO = structural_simplify(HO)

```

I get 4 equations and the following solution.

 ![plot_2](https://global.discourse-cdn.com/julialang/original/3X/5/e/5ec296a69717076b7478e1b273a2a9502901daab.jpeg)

thank you all for the help.