I am interested in manipulating symbolic differential equations, but can’t seem to identify a package that provides the functionality I am looking for. I am particularly interested in this functionality for defining a dynamical system in the form x’(t) = f(x(t)). For example, given a set of ODEs, how would I use Julia to symbolically manipulate them to get them into the form x’(t) = f(x(t))?
1 Like
This is your go-to package, without discussion.
There doesn’t seem to be an explicit function for solving a set of symbolic equations. I feel like I am overlooking something. Is the idea behind ModelingToolkit that a set of equations is used to construct an ODESystem, which can then be used to construct an ODEProblem, etc.? Thanks for the quick response!
Not all ODEs can be written in that form
Not all ODEs can be written in such a separable form, e.g. x'(t) = sin(x(t) + t). Can you give an example of the input and output you want?
4 Likes
Good point. I should have specified that I am only interested in autonomous systems. It is clear now that my initial post lacked some detail which is probably important. I need to derive the equations of motion for a system symbolically, e.g. via Lagrangian mechanics, and solve the resulting symbolic equations so that I can then use the tools in DifferentialEquations, BoundaryValueDiffEq, etc., to solve numerical problems. The systems I am working with have ODEs which are rather unpleasant to look at, so I am looking for a way to make computing the equations of motion autonomous so that the output is a system of the form x'(t) = f(x(t)). The input, I suppose, would be a symbolic expression (a Lagrangian).
ModelingToolkit is working towards this but isn’t all of the way there for Legrangians right now. Generally that needs state selection. I think we’ll have demos along those lines in a month or so (and it would probably be nice to have a system type for this).
2 Likes
Here’s something I cooked up the other day to do what I think you’re describing. It’s probably not super robust and definitely not the “correct way” to do this, but it worked for the problem I tried and it gives me the correct results so ¯\(ツ)/¯
using Symbolics, ModelingToolkit, DifferentialEquations
∂_∂(x) = y -> expand_derivatives(Differential(x)(y))
function LagrangeEOM(L, q̇, q, p, t;
Q = zeros(length(q)),
default_u0 = [q̇..., q...] .=> 0.0,
kwargs...)
Q_vals = Q
inds = 1:length(q)
@variables q̈[inds](t) Q[inds](t)
d_dt = ∂_∂(t)
# Equations of motion
lagrange_eqs = map(q̇, q, Q) do q̇ᵢ, qᵢ, Qᵢ
d_dt(∂_∂(q̇ᵢ)(L)) - ∂_∂(qᵢ)(L) ~ Qᵢ
end
try
subs = d_dt.([q̇; q]) .=> [q̈; q̇]
new_eqs = [substitute(eq.lhs, subs) ~ eq.rhs for eq in lagrange_eqs]
lagrange_eqs = d_dt.(q̇) .~ Symbolics.solve_for(new_eqs, q̈)
catch e
@warn "Could not reduce equations due to: $e"
end
deriv_eqs = d_dt.(q) .~ q̇
force_eqs = Q .~ Q_vals
eqs = [lagrange_eqs; deriv_eqs; force_eqs]
return structural_simplify(ODESystem(eqs, t, [q̇..., q..., Q...], p; default_u0, kwargs...))
end
Now using it for a model with a sloshing liquid on a cart (basically just a pendulum on a cart), you have:
# Set up symbolic variables and parameters
@variables t x(t) v(t) θ(t) ω(t) F(t)
@parameters Lₐ Lₚ mₖ mₗ mₚ g
D = Differential(t)
## Inputs
# Make Lagrangian
V = -mₚ*Lₚ*g*cos(θ)
T = (mₖ+mₗ)*v^2/2 + mₚ*(Lₚ^2*ω^2 + 2*Lₚ*ω*v*cos(θ) +v^2)/2
L = T - V
# Cart input force
F = 1000sin(t)
# Generalized forces
Q = [F, 0]
# Make equations of motion
slosh_cart = LagrangeEOM(L, [v, ω], [x, θ], [Lₐ, Lₚ, mₖ, mₗ, mₚ, g], t; Q)
# Initial Conditions
ic = [
θ => deg2rad(10)
]
# Parameters
p = Dict([
Lₐ => 10
Lₚ => 0.5
mₖ => 100
mₗ => 0
mₚ => 25
g => 9.80665
])
## Simulation
sol = solve(ODEProblem(slosh_cart, ic, (0.0, 10.0), [p...]))
Like I said, no clue how robust this actually is. And I haven’t really set it up to handle constraints or anything. But it at least works for this simple model.
6 Likes
It’s the right way to do it, and would be a nice PR IMO (I would like to get HamiltonianSystems, Lagrange, etc.). The only issue is that this form can require state selection handling in the nonlinear tearing of structural_simplify, which @YingboMa is implementing but won’t complete until maybe the end of next week. After that, then this will be robust. Before that, there will be some odd cases that just won’t simplify all of the way on this kind of code.
But can you open an issue so we can discuss these kinds of alternative constructors?
5 Likes
I will try this out. Thank you!
1 Like
Here’s a more explicit version using @jonniedie 's example.
using ModelingToolkit, OrdinaryDiffEq
@variables t x(t) v(t) θ(t) ω(t) F(t)
@parameters Lₐ Lₚ mₖ mₗ mₚ g
D = Differential(t)
function L((v, ω), (x, θ), (Lₐ, Lₚ, mₖ, mₗ, mₚ, g), t)
V = -mₚ*Lₚ*g*cos(θ)
T = (mₖ+mₗ)*v^2/2 + mₚ*(Lₚ^2*ω^2 + 2*Lₚ*ω*v*cos(θ) + v^2)/2
T - V
end
function lagrangian2system(
L, q̇, q, p, t;
Q = zeros(length(q)),
defaults = [q̇; q] .=> 0.0,
kwargs...
)
Q_vals = Q
inds = eachindex(q)
@variables v[inds] x[inds] Q[inds](t)
sub = Base.Fix2(substitute, Dict([v.=>q̇; x.=>q]))
Lf = L(v, x, p, t)
F = ModelingToolkit.gradient(Lf, x) + Q
Lᵥ = ModelingToolkit.gradient(Lf, v)
rhs = sub.(F - ModelingToolkit.jacobian(Lᵥ, x) * q̇ - ModelingToolkit.derivative.(Lᵥ, (t,)))
M = sub.(ModelingToolkit.jacobian(Lᵥ, v))
D = Differential(t)
eqs = [
D.(q̇) .~ M \ rhs
D.(q) .~ q̇
Q .~ Q_vals
]
sys = ODESystem(eqs, t, [q̇; q; Q], p; defaults=defaults, kwargs...)
return structural_simplify(sys)
end
# Cart input force
F = 1000sin(t)
# Generalized forces
Q = [F, 0]
# Make equations of motion
slosh_cart = lagrangian2system(L, [v, ω], [x, θ], [Lₐ, Lₚ, mₖ, mₗ, mₚ, g], t; Q)
# Initial Conditions
ic = [
θ => deg2rad(10)
]
# Parameters
p = Dict([
Lₐ => 10
Lₚ => 0.5
mₖ => 100
mₗ => 0
mₚ => 25
g => 9.80665
])
## Simulation
prob = ODEProblem(slosh_cart, ic, (0.0, 10.0), [p...])
sol = solve(prob, Tsit5())
6 Likes
Hey,
here’s a n-pendulum example using @YingboMa 's lagrangian2system function:
## Define symbolic variables
n = 3
@variables t θ[1:n](t) ω[1:n](t)
@parameters g l m
## Define Lagrangian
function L(ω, θ, (g, l, m), t)
# Moments of inertia
J = m * l^2 /12
# Define vertical positions
y = vcat([l * cos(θ[1]) / 2], map(i -> sum(map(j -> l*cos(θ[j]), 1:i-1)) + l * cos(θ[i]) / 2, 2:n))
# Define velocities
vx = vcat([l * cos(θ[1]) * ω[1]/ 2], map(i -> sum(map(j -> l*cos(θ[j])*ω[j], 1:i-1)) + l * cos(θ[i]) * ω[i] / 2, 2:n))
vy = vcat([-l * sin(θ[1]) * ω[1] / 2], map(i -> sum(map(j -> -l*sin(θ[j])*ω[j], 1:i-1)) - l * sin(θ[i]) * ω[i] / 2, 2:n))
# Build Lagrangian
V = -sum(m.*y.*g)
T = 1/2 * sum(m.*(vx.^2 + vy.^2) + J.*(ω.^2))
T - V
end
# Friction terms
Q = zero(ω)
# Define ODESystem
n_pendulum = lagrangian2system(L, ω, θ, [g, l, m], t; Q)
# Initial Conditions
ic = vcat(
ω .=> zeros(n),
θ .=> deg2rad.(randn(n)*10)
)
# Parameters
p = Dict([
g => 9.806
l => 1.0
m => 1.0
])
## Simulation
prob = ODEProblem(n_pendulum, ic, (0.0, 10.0), [p...])
sol = solve(prob, Tsit5())
Plots.plot(sol)
The above works well, but as soon as I try to pass vector valued parameters m[1:n] and l[1:n] I get a “syntax: invalid let syntax” error in the ODEFunction call. Here’s the code for this:
using Symbolics, ModelingToolkit, DifferentialEquations
using Latexify, Plots
## Lagrangian functionality
function lagrangian2system(
L, q̇, q, p, t;
Q = zeros(length(q)),
defaults = [q̇; q] .=> 0.0,
kwargs...
)
Q_vals = Q
inds = eachindex(q)
@variables v[inds] x[inds] Q[inds](t)
sub = Base.Fix2(substitute, Dict([v.=>q̇; x.=>q]))
Lf = L(v, x, p, t)
F = ModelingToolkit.gradient(Lf, x) + Q
Lᵥ = ModelingToolkit.gradient(Lf, v)
rhs = sub.(F - ModelingToolkit.jacobian(Lᵥ, x) * q̇ - ModelingToolkit.derivative.(Lᵥ, (t,)))
M = sub.(ModelingToolkit.jacobian(Lᵥ, v))
D = Differential(t)
eqs = [
D.(q̇) .~ M \ rhs
D.(q) .~ q̇
Q .~ Q_vals
]
sys = ODESystem(eqs, t, [q̇; q; Q], p; defaults=defaults, kwargs...)
return structural_simplify(sys)
end
## Define symbolic variables
n = 3
@variables t θ[1:n](t) ω[1:n](t)
@parameters g l[1:n] m[1:n]
## Define Lagrangian
function L(ω, θ, (g, l, m), t)
# Moments of inertia
J = m .* l.^2 /12
# Define vertical positions
y = vcat([l[1] * cos(θ[1]) / 2], map(i -> sum(map(j -> l[j]*cos(θ[j]), 1:i-1)) + l[i] * cos(θ[i]) / 2, 2:n))
# Define velocities
vx = vcat([l[1] * cos(θ[1]) * ω[1]/ 2], map(i -> sum(map(j -> l[j]*cos(θ[j])*ω[j], 1:i-1)) + l[i] * cos(θ[i]) * ω[i] / 2, 2:n))
vy = vcat([-l[1] * sin(θ[1]) * ω[1] / 2], map(i -> sum(map(j -> -l[j]*sin(θ[j])*ω[j], 1:i-1)) - l[i] * sin(θ[i]) * ω[i] / 2, 2:n))
# Build Lagrangian
V = -sum(m.*y.*g)
T = 1/2 * sum(m.*(vx.^2 + vy.^2) + J.*(ω.^2))
T - V
end
# Friction terms
Q = zero(ω)
# Define ODESystem
n_pendulum = lagrangian2system(L, ω, θ, [g, l, m], t; Q)
# Initial Conditions
ic = vcat(
ω .=> zeros(n),
θ .=> deg2rad.(randn(n)*10)
)
# Parameters
p = Dict([
g => 9.806
l => repeat([1.0], n)
m => repeat([1.0], n)
])
## Simulation
prob = ODEProblem(n_pendulum, ic, (0.0, 10), [p...])
sol = solve(prob, Tsit5())
Plots.plot(sol)
Any suggestions how to properly pass those parameters?
3 Likes