ModelingToolkit.jl
ModelingToolkit.jl copied to clipboard
SciML
Rules seems to remove some metadata when rewriting equations
Hi, John here again sorry for a long post but I include everything for replication and (For context) (Title stolen from one of Yingbos comments) I have the following Modelica model:
model Pendulum
parameter Real x0 = 10;
parameter Real y0 = 10;
parameter Real g = 9.81;
parameter Real L = sqrt(x0^2 + y0^2);
Real x(start = x0);
Real y(start = y0);
Real vx;
Real vy;
Real phi(start = 1., fixed = true);
Real phid;
equation
x = L * sin(phi);
y = -L * cos(phi);
der(phi) = phid;
der(x) = vx;
der(y) = vy;
der(phid) = -g / L * sin(phi);
end Pendulum;
From this I generate this model:
using ModelingToolkit
using DifferentialEquations
begin
begin
saved_values_Pendulum = SavedValues(Float64, Tuple{Float64,Array})
function PendulumCallbackSet(aux)
local p = aux[1]
local reals = aux[2]
nothing
return CallbackSet()
end
end
function PendulumModel(tspan = (0.0, 1.0))
ModelingToolkit.@variables t #= C:\Users\johti17\Projects\Programming\JuliaPackages\OM.jl\OMBackend.jl\src\CodeGeneration\MTK_CodeGeneration.jl:235 =#
der = ModelingToolkit.Differential(t)
parameters = ModelingToolkit.@parameters((x0, y0, g, L)) #= C:\Users\johti17\Projects\Programming\JuliaPackages\OM.jl\OMBackend.jl\src\CodeGeneration\MTK_CodeGeneration.jl:237 =#
begin
function generateStateVariables()
(:t, :(x(t)), :(y(t)), :(phi(t)), :(phid(t)))
end
function generateAlgebraicVariables()
(:t, :(vx(t)), :(vy(t)))
end
variableConstructors =
Function[generateStateVariables, generateAlgebraicVariables]
end
componentVars = []
for constructor in variableConstructors
res = eval(
ModelingToolkit.Symbolics._parse_vars(
"CustomCall",
Real,
constructor(),
),
)
push!(componentVars, res[2:end])
end
vars = collect(Iterators.flatten(componentVars))
pars = Dict(
x0 => float(10.0),
y0 => float(10.0),
g => float(9.81),
L => float(sqrt(x0^2.0 + y0^2.0)),
)
startEquationComponents = []
begin
startEquationConstructors = Function[]
begin
function generateStartEquations0()
[
vx => 0.0,
vy => 0.0,
x => pars[x0],
y => pars[y0],
phi => 1.0,
phid => 0.0,
]
end
push!(startEquationConstructors, generateStartEquations0)
end
end
for constructor in startEquationConstructors
push!(startEquationComponents, constructor())
end
initialValues = collect(Iterators.flatten(startEquationComponents))
equationComponents = []
begin
equationConstructors = Function[]
begin
x0 = float(10.0)
y0 = float(10.0)
g = float(9.81)
L = float(sqrt(x0^2.0 + y0^2.0))
#=Note I rewrite the equations by manipulating the Julia ast here (my old solution)=#
function generateEquations0()
[
0 ~ x - L * sin(phi),
0 ~ y - -(L * cos(phi)),
der(phi) ~ phid,
der(x) ~ vx,
der(y) ~ vy,
der(phid) ~ - ((g / L) * sin(phi)),
]
end
push!(equationConstructors, generateEquations0)
end
end
for constructor in equationConstructors
push!(equationComponents, constructor())
end
eqs = collect(Iterators.flatten(equationComponents))
nonLinearSystem = ModelingToolkit.ODESystem(
eqs,
t,
vars,
parameters,
name = :($(Symbol("Pendulum"))),
)
firstOrderSystem = ModelingToolkit.ode_order_lowering(nonLinearSystem)
reducedSystem = structural_simplify(firstOrderSystem)
local event_p = [10.0, 10.0, 9.81, 0]
local discreteVars = collect(values(Dict([])))
local event_vars = vcat(collect(values(Dict([initialValues...]))), discreteVars)
local aux = Array{Array{Float64}}(undef, 2)
aux[1] = event_p
aux[2] = event_vars
problem = ModelingToolkit.ODEProblem(
reducedSystem,
initialValues,
tspan,
pars,
callback = PendulumCallbackSet(aux),
)
return (problem, initialValues, reducedSystem, tspan, pars, vars)
end
end
(PendulumModel_problem, ivs, rs, tspan, pars, vars) = PendulumModel()
function PendulumSimulate(tspan)
solve(PendulumModel_problem, tspan = tspan)
end
function PendulumSimulate(tspan = (0.0, 1.0); solver = Rodas5())
solve(PendulumModel_problem, tspan = tspan, solver)
end
julia> full_equations(rs)
4-element Vector{Symbolics.Equation}:
Differential(t)(phi(t)) ~ phid(t)
Differential(t)(x(t)) ~ 14.142135623730951phid(t)*cos(phi(t))
Differential(t)(y(t)) ~ 14.142135623730951phid(t)*sin(phi(t))
Differential(t)(phid(t)) ~ -0.6936717523440031sin(phi(t))
Everything seems to work fine:) However, the way I used to do the rewriting of the equation was slow, so I "improved it" by stealing some code from #1164
However,
"""
Temporary rewrite function. Not very pretty...
Original code by Chris R. Modifed it to fix terms of type X * D(Y).
The solution to solve it is not pretty and is probably quite flaky.
"""
function move_diffs(eq::Equation; rewrite)
# Do not modify `D(x) ~ ...`, already correct
# Ignore `δ(x) ~ ...` for now
res = if !(eq.lhs isa Term && operation(eq.lhs) isa Differential) &&
!(eq.lhs isa Term && operation(eq.lhs) isa Difference)
_eq = eq.rhs-eq.lhs
rhs = rewrite(_eq)
if rhs === nothing
eq
end
lhs = _eq - rhs
if !(lhs isa Number) && (operation(lhs) isa Differential)
lhs ~ -rhs
elseif !(lhs isa Number) && (operation(lhs) == *)
#= This code is probably quite flaky, however, it should not be needed after similar things are introduced in MTK. =#
newRhs = substitute(rhs, rhs => rhs / arguments(lhs)[1])
newLhs = substitute(lhs, lhs => arguments(lhs)[2])
#lhs = lhs_
~(newLhs, -newRhs)
else
-lhs ~ rhs
end
else
eq
end
return res
end
This function is called with:
function makeODESystem(deqs, iv, vars, pars; name)
# Equation is not in the canonical form
# Use a rule to find all g(u)*D(u) terms
# and move those to the lhs
D = Differential(iv)
r1 = SymbolicUtils.@rule ~~a * D(~~b) * ~~c => 0
r2 = SymbolicUtils.@rule D(~~b) => 0
# debugRewrite(deqs, iv, vars, parameters)
remove_diffs = SymbolicUtils.Postwalk(SymbolicUtils.Chain([r1,r2]))
deqs = [move_diffs(eq, rewrite = remove_diffs) for eq in deqs]
# debugRewrite(deqs, iv, vars, parameters)
res = ModelingToolkit.ODESystem(deqs, iv, vars, pars; name = name)
return res
end
This works fine for most of my tests, however, when I instead rewrite my system using these rules things go wrong..
This is before my new transformation:
@variables(x(t),y(t),phi(t),phid(t),vx(t),vy(t),)
eqs =[
0 ~ x(t) - 14.142135623730951sin(phi(t)),
0 ~ 14.142135623730951cos(phi(t)) + y(t),
0 ~ Differential(t)(phi(t)) - phid(t),
0 ~ Differential(t)(x(t)) - vx(t),
0 ~ Differential(t)(y(t)) - vy(t),
0 ~ 0.6936717523440031sin(phi(t)) + Differential(t)(phid(t)),
]
@parameters(x0,y0,g,L,)
and this is after:
@variables(x(t),y(t),phi(t),phid(t),vx(t),vy(t),)
eqs =[
0 ~ x(t) - 14.142135623730951sin(phi(t)),
0 ~ 14.142135623730951cos(phi(t)) + y(t),
Differential(t)(phi(t)) ~ phid(t),
Differential(t)(x(t)) ~ vx(t),
Differential(t)(y(t)) ~ vy(t),
Differential(t)(phid(t)) ~ -0.6936717523440031sin(phi(t)),
]
@parameters(x0,y0,g,L,)
However, this particular model needs to have its index reduced, and my compiler calls your index reduction routine if it is needed, using this new technique I get:
[ Info: Interactive evaluation failed: ModelingToolkit.InvalidSystemException("The system is missing an equation for vx(t).") with mode: MTK_MODE
[ Info: ModelingToolkit.InvalidSystemException("The system is missing an equation for vx(t).")
ERROR: InvalidSystemException: The system is missing an equation for vx(t).
Stacktrace:
[1] simulateModel(modelName::String; MODE::OMBackend.BackendMode, tspan::Tuple{Float64, Float64}, solver::Expr)
@ OMBackend C:\Users\johti17\Projects\Programming\JuliaPackages\OM.jl\OMBackend.jl\src\backendAPI.jl:343
[2] #runModelFM#2
@ C:\Users\johti17\Projects\Programming\JuliaPackages\OM.jl\src\OM.jl:97 [inlined]
[3] macro expansion
@ .\timing.jl:220 [inlined]
[4] runModelMTK(model::String, file::String; timeSpan::Tuple{Float64, Float64})
@ Main C:\Users\johti17\Projects\Programming\JuliaPackages\OM.jl\test\runtests.jl:120
[5] top-level scope
@ .\timing.jl:220 [inlined]
[6] top-level scope
@ .\REPL[78]:0
caused by: InvalidSystemException: The system is missing an equation for vx(t).
Interestingly enough, this error only occurs when conducting index reduction, either via dae_index_lowering or via structural_simplify for models that do not require index reduction the models are simulating just fine
The complete model using the new rewrite scheme:
using ModelingToolkit
using DifferentialEquations
begin
begin
saved_values_Pendulum = SavedValues(Float64, Tuple{Float64,Array})
function PendulumCallbackSet(aux)
local p = aux[1]
local reals = aux[2]
nothing
return CallbackSet()
end
end
function PendulumModel(tspan = (0.0, 1.0))
ModelingToolkit.@variables t
der = ModelingToolkit.Differential(t)
parameters = ModelingToolkit.@parameters((x0, y0, g, L))
begin
function generateStateVariables()
(:t, :(x(t)), :(y(t)), :(phi(t)), :(phid(t)))
end
function generateAlgebraicVariables()
(:t, :(vx(t)), :(vy(t)))
end
variableConstructors =
Function[generateStateVariables, generateAlgebraicVariables]
end
componentVars = []
for constructor in variableConstructors
res = eval(
ModelingToolkit.Symbolics._parse_vars(
"CustomCall",
Real,
constructor(),
),
)
push!(componentVars, res[2:end])
end
vars = collect(Iterators.flatten(componentVars))
pars = Dict(
x0 => float(10.0),
y0 => float(10.0),
g => float(9.81),
L => float(sqrt(x0^2.0 + y0^2.0)),
)
startEquationComponents = []
begin
startEquationConstructors = Function[]
begin
function generateStartEquations0()
[
vx => 0.0,
vy => 0.0,
x => pars[x0],
y => pars[y0],
phi => 1.0,
phid => 0.0,
]
end
push!(startEquationConstructors, generateStartEquations0)
end
end
for constructor in startEquationConstructors
push!(startEquationComponents, constructor())
end
initialValues = collect(Iterators.flatten(startEquationComponents))
equationComponents = []
begin
equationConstructors = Function[]
begin
x0 = float(10.0)
y0 = float(10.0)
g = float(9.81)
L = float(sqrt(x0^2.0 + y0^2.0))
function generateEquations0()
[
0 ~ x - L * sin(phi),
0 ~ y - -(L * cos(phi)),
0 ~ der(phi) - phid,
0 ~ der(x) - vx,
0 ~ der(y) - vy,
0 ~ der(phid) - -((g / L) * sin(phi)),
]
end
push!(equationConstructors, generateEquations0)
end
end
for constructor in equationConstructors
push!(equationComponents, constructor())
end
eqs = collect(Iterators.flatten(equationComponents))
nonLinearSystem = OMBackend.CodeGeneration.makeODESystem(
eqs,
t,
vars,
parameters,
name = :($(Symbol("Pendulum"))),
)
firstOrderSystem = nonLinearSystem
reducedSystem =
ModelingToolkit.structural_simplify(firstOrderSystem; simplify = true)
local event_p = [10.0, 10.0, 9.81, 0]
local discreteVars = collect(values(Dict([])))
local event_vars = vcat(collect(values(Dict([initialValues...]))), discreteVars)
local aux = Array{Array{Float64}}(undef, 2)
aux[1] = event_p
aux[2] = event_vars
problem = ModelingToolkit.ODEProblem(
reducedSystem,
initialValues,
tspan,
pars,
callback = PendulumCallbackSet(aux),
)
return (problem, initialValues, reducedSystem, tspan, pars, vars)
end
end
(PendulumModel_problem, ivs, PendulumModel_ReducedSystem, tspan, pars, vars) =
PendulumModel()
function PendulumSimulate(tspan)
solve(PendulumModel_problem, tspan = tspan)
end
function PendulumSimulate(tspan = (0.0, 1.0); solver = Rodas5())
solve(PendulumModel_problem, tspan = tspan, solver)
end
Note the call to:
nonLinearSystem = OMBackend.CodeGeneration.makeODESystem(
eqs,
t,
vars,
parameters,
name = :($(Symbol("Pendulum"))))
This could be replaced with the rewrite method above, while this rewriting does seem to work for most cases structural_simplify fails and I can't understand why full equations returns the same results for both systems
Cheers,
John
Some more context, For now, I have managed to resolve the issue, but not in a "pretty" way:
function makeODESystem(deqs, iv, vars, pars, idxReduction; name)
# Equation is not in the canonical form
# Use a rule to find all g(u)*D(u) terms
# and move those to the lhs
D = Differential(iv)
r1 = SymbolicUtils.@rule ~~a * D(~~b) * ~~c => 0
r2 = SymbolicUtils.@rule D(~~b) => 0
debugRewrite(deqs, iv, vars, pars)
remove_diffs = SymbolicUtils.Postwalk(SymbolicUtils.Chain([r1,r2]))
deqs = [move_diffs(eq, rewrite = remove_diffs) for eq in deqs]
#= Special computationally heavy routine for systems that need index reduction.. =#
if idxReduction
res = debugRewrite(deqs, iv, vars, pars)
println(res)
expr = Meta.parse(res)
eval(expr)
res = ModelingToolkit.ODESystem(eqs2, iv, vars2, pars; name = name)
res2 = ModelingToolkit.dae_index_lowering(res)
return res2
end
return ModelingToolkit.ODESystem(deqs, iv, vars, pars; name = name)
end
Basically, since this only seems to fail for the system that needs index reduction I convert the system to a string and reparse and reevaluate it, this made all my tests pass, still, it seems that MTK system lose some kind of information when you apply rewrite rules on them
