# ModelingToolkit compiler does not simplify as expected

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978>\
**Category:** Modelling & Simulations\
**Tags:** question, modelingtoolkit\
**Created:** [January 10, 2026, 10:41am UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978 "2026-01-10T10:41:39Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![PhysicsBender](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/physicsbender/32/220077_2.png) [@PhysicsBender](https://discourse.julialang.org/u/PhysicsBender)\
**Post date:** [January 10, 2026, 10:41am UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/1 "2026-01-10T10:41:39Z")

</div>

I have noticed that the ModelingToolkit compiler does not simplify systems as expected, creating more unknowns and differential equations. At first I thought it was because it was finding linear systems that it was unable to algebraically resolve but I was able to reproduce problems even with a simple Mass Spring system using Standard Library components. Is this intended behaviour? Coming from Modelica this seems to me like a drawback when considering that Modelica tools like dymola can reduce even multibody systems to the minimum configuration variables needed and only have problems when kinematic loops are present which is expected.

**Minimal Example**

```julia-auto
using ModelingToolkit
using ModelingToolkitStandardLibrary
@parameters t

D = Differential(t)

@component function TwoMassTranslationalSystem(; name )
    @named fixed = ModelingToolkitStandardLibrary.Mechanical.Translational.Fixed()
    @named spring = ModelingToolkitStandardLibrary.Mechanical.Translational.Spring(k = 1.0)
    @named mass = ModelingToolkitStandardLibrary.Mechanical.Translational.Mass(m = 1.0)
    eqs = Equation[
        connect(fixed.flange, spring.flange_a),
        connect(spring.flange_b, mass.flange),
    ]
    return System(eqs, t; name = name, systems = [fixed, spring, mass])
end

 @named spring_damper = TwoMassTranslationalSystem()
sys = mtkcompile(spring_damper)
unknowns(sys)
full_equations(sys)

```

**Julia Output**

```julia-auto

3-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 mass₊v(t)
 mass₊s(t)
 spring₊delta_s(t)

```

```julia-auto
3-element Vector{Equation}:
 Differential(t, 1)(mass₊v(t)) ~ mass₊g + (spring₊k*spring₊delta_s(t)) / mass₊m
 Differential(t, 1)(mass₊s(t)) ~ mass₊v(t)
 Differential(t, 1)(spring₊delta_s(t)) ~ -mass₊v(t)

```

---

<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 10, 2026, 10:54am UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/2 "2026-01-10T10:54:15Z")

</div>

This is because the translational library is using velcotities in the connectors instead of positions, there is an alternative implementation available called `ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica` which should not suffer from these issues, but it lacks some components compared to `Translational`. The dyad version of this library should also be free from this issue.

---

<div class="post-metadata">

**Author:** ![PhysicsBender](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/physicsbender/32/220077_2.png) [@PhysicsBender](https://discourse.julialang.org/u/PhysicsBender)\
**Post date:** [January 10, 2026, 11:35am UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/3 "2026-01-10T11:35:31Z")

</div>

> [@baggepinnen](#):
>
> TranslationalModelica

I see. I switched to using the ModelicaComponent library and I got two equations. However I tried extending the system by adding an additional mass just to see if it would correcty deduce the system into the same unknowns albeit with added mass and I got an extra algebraic equation that looks easily solvable. What is tripping up the symbolic engine in this case ? Is it avoiding dividing by zero cause the mass could be parametrized as zero?

```julia-auto
@component function TwoMassTranslationalSystem(; name )
    @named fixed = ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica.Fixed()
    @named spring = ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica.Spring(c = 1.0)
    @named mass1 = ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica.Mass(m = 1.0)
    @named mass2 = ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica.Mass(m = 1.0)

    eqs = Equation[
        connect(fixed.flange, spring.flange_a),
        connect(spring.flange_b, mass1.flange_a),
        connect(spring.flange_b, mass2.flange_a),
    ]
    return System(eqs, t; name = name, systems = [fixed, spring, mass1, mass2])
end

 @named spring_damper = TwoMassTranslationalSystem()
sys = mtkcompile(spring_damper)

```

```julia-auto
julia> unknowns(sys)
3-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 mass2₊flange_a₊s(t)
 mass2₊flange_a₊sˍt(t)
 mass1₊flange_a₊sˍtt(t)

```

```julia-auto
julia> full_equations(sys)
3-element Vector{Equation}:
 Differential(t, 1)(mass2₊flange_a₊s(t)) ~ mass2₊flange_a₊sˍt(t)
 Differential(t, 1)(mass2₊flange_a₊sˍt(t)) ~ mass1₊flange_a₊sˍtt(t)
 0 ~ (-(-fixed₊s0 - spring₊s_rel0 + mass2₊flange_a₊s(t))*spring₊c - mass1₊m*mass1₊flange_a₊sˍtt(t)) / mass2₊m - mass1₊flange_a₊sˍtt(t)

```

---

<div class="post-metadata">

**Author:** ![PhysicsBender](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/physicsbender/32/220077_2.png) [@PhysicsBender](https://discourse.julialang.org/u/PhysicsBender)\
**Post date:** [January 10, 2026, 11:44am UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/4 "2026-01-10T11:44:26Z")

</div>

I tried a custom implementation with masses hardcoded to 1 and it’s still creates an 3rd unknown

---

<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 10, 2026, 12:25pm UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/5 "2026-01-10T12:25:38Z")

</div>

You need some new, not yet default, keyword args to `mtkcompile`

```julia-auto
julia> sys = mtkcompile(spring_damper; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))
Model spring_damper:
Equations (2):
  2 standard: see equations(spring_damper)
Unknowns (2): see unknowns(spring_damper)
  mass2₊flange_a₊s(t)
  mass2₊flange_a₊sˍt(t)

```

---

<div class="post-metadata">

**Author:** ![PhysicsBender](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/physicsbender/32/220077_2.png) [@PhysicsBender](https://discourse.julialang.org/u/PhysicsBender)\
**Post date:** [January 10, 2026, 12:26pm UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/6 "2026-01-10T12:26:57Z")

</div>

Thanks a lot. Problems like these have been bugging me for days.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 10, 2026, 12:28pm UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/7 "2026-01-10T12:28:08Z")

</div>

> [@PhysicsBender](#):
>
> have been bugging me for days

If I don’t find a solution to an issue myself within 20 minutes, I ask for help here.

---

<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 10, 2026, 12:28pm UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/8 "2026-01-10T12:28:13Z")

</div>

The options I just showed you there are really brand new, and there are still some kinks to be worked out. It’s the development of Multibody.jl that is driving these changes, so once a new release of Multibody.jl is out it’ll all be tested a bit more.

---

<div class="post-metadata">

**Author:** ![PhysicsBender](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/physicsbender/32/220077_2.png) [@PhysicsBender](https://discourse.julialang.org/u/PhysicsBender)\
**Post date:** [January 11, 2026, 12:37pm UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/9 "2026-01-11T12:37:54Z")

</div>

Hey @baggepinnen since I am a bit impatient I start experimenting with some custom multibody components I created and the compiler options you suggested and I had some intresting observations.

I created a custom `Frame()` connector with torque and force vectors as flow variables and euler angles and inertially expressed position vector as non flow variables ( I understand the limitations of using euler angles as part of the connector and for a proper solution that scales creating an orientation object that encodes both rotation matrix and angular velocity is propably beneficial but for the purposes of my toy example I think it should work) . I also created a simple `Body()` compoent that just writes the kinematics + Euler/Newton equations for dynamic equilibrium.

Depending on how I formulate the angular velocity kinematics in the model I get different response from the symbolic engine and compiler. I will try to list them into four scenarios, two of which contain only one body which should reduce to 12 unknowns, and another two which contain two bodies rigidly connected but should again reduce to 12 unknowns with summed inertia and mass.

1 . For the `one_body` system in which I explicitly wrote the kinematics that transform angular velocity to euler angle rates (`use_w = true`) the compile engine correctly takes the “simplest” approach and “roots” angular velocity picking them as the state which leads to symbolically lean equations.This creates 12 differential and 0 algebraic equations. Happy Days!

2 . For the `one_body_alt` system in which I explicitly wrote the inverse kinematics that transform euler angle rates to angular velocity (`use_w = false`) the compile engine picks the euler angle rates as states and then in order to calculate angular velocity components x and y creates two extra algebraic equations. I hypothesized this happens because of the singularity created at pitch angle equal to ± pi/2 but I could be wrong cause in the next example it does things that could lead to divisions by zero under certain circumastances. This creates 12 differential and 2 algebraic equations.

3 . For the `two_bodies` system in which `use_w = true` the compile engine strangely now picks the euler angle rates as states and then in order to calculate angular velocity components x and y creates two extra algebraic equations for each body similar to case 2. This creates 12 differential and 4 algebraic equations.

4 . For the `two_bodies_alt` system in which `use_w = false` the compile engine again picks the euler angle rates as states and manages to symbolically resolve the angular velocities without introducing any algebraic equations. However the expressions created are enormous, which for sure not lead to efficient code. This creates 12 differential and 0 algebraic equations.

It seems that depending one how you write the equations the same physics can lead to different sets of equations which is not always the simplest. It would be intresting if there was a way to force the compiler to pick certain variables as root/states just like you can do in Modelica. I tried using the `state_priority` metadata but it seems to be completely ignored by `mtkcompile`.

Below is the code showcasing all four cases.

```julia-auto
using ModelingToolkit
using LinearAlgebra

@parameters t
D = Differential(t)

@connector function Frame(; name , r_0 = nothing)
    @variables begin 
        r_0(t)[1:3] = r_0
        θ(t)[1:3] 
        f(t)[1:3],[connect = Flow]
        tau(t)[1:3],[connect = Flow]
    end
    (r_0, f, tau, θ) = collect.(( r_0, f, tau, θ))
    return System(Equation[], t, [r_0; f; tau; θ],[]; name = name)
end

@component function Body(; name, use_w = false)
    @named frame_a = Frame(r_0 = [0,0,0])
    @parameters begin
        m = 1, [description = "mass of the body (kg)"]
        Ixx = 0.1, [description = "moment of inertia around the body x-axis (kg⋅m²)"]
        Iyy = 0.1, [description = "moment of inertia around the body y-axis (kg⋅m²)"]
        Izz = 0.1, [description = "moment of inertia around the body z-axis (kg⋅m²)"]
        g = 9.81, [description = "gravity vector, resolved in inertial frame"]
    end
    g_0 = [0, 0, g]
    I = [Ixx 0 0; 0 Iyy 0; 0 0 Izz]
    @variables begin
        v_0(t)[1:3] 
        θ_t(t)[1:3] 
        w_a(t)[1:3] 
    end 
    (θ, r_0, f_a, tau_a, v_0, θ_t, w_a) = collect.((frame_a.θ, frame_a.r_0, frame_a.f, frame_a.tau, v_0, θ_t, w_a))
    R_a = [
        cos(θ[2])*cos(θ[3]) -cos(θ[1])*sin(θ[3]) + sin(θ[1])*sin(θ[2])*cos(θ[3]) sin(θ[1])*sin(θ[3]) + cos(θ[1])*sin(θ[2])*cos(θ[3]);
        cos(θ[2])*sin(θ[3]) cos(θ[1])*cos(θ[3]) + sin(θ[1])*sin(θ[2])*sin(θ[3]) -sin(θ[1])*cos(θ[3]) + cos(θ[1])*sin(θ[2])*sin(θ[3]);
        -sin(θ[2]) sin(θ[1])*cos(θ[2]) cos(θ[1])*cos(θ[2])
    ] # assuming ZYX Euler angles from local to global
    θ_t_exp = [
        w_a[1] + sin(θ[1]) * tan(θ[2]) * w_a[2] + cos(θ[1]) * tan(θ[2]) * w_a[3],
                             cos(θ[1]) * w_a[2] - sin(θ[1]) * w_a[3],
                 sin(θ[1]) / cos(θ[2]) * w_a[2] + cos(θ[1]) / cos(θ[2]) * w_a[3]
    ] # assuming ZYX Euler angles   
    w_a_exp = [
            θ_t[1] - θ_t[3]*sin(θ[2]),
            θ_t[2]*cos(θ[1]) + θ_t[3]*cos(θ[2])*sin(θ[1]),
           -θ_t[2]*sin(θ[1]) + θ_t[3]*cos(θ[2])*cos(θ[1])
    ] # assuming ZYX Euler angles
    ang_eqs = if use_w
        θ_t .~ θ_t_exp
     else
        w_a .~ w_a_exp
    end
    eqs = [
        ang_eqs
        D.(θ) .~ θ_t
        D.(r_0) .~ v_0     
        m*(D.(v_0) + g_0) .~ R_a*f_a # expressed in world frame
        I*D.(w_a) + cross(w_a, I*w_a) .~ tau_a # expressed in frame_a (body frame)
    ]
    return System(eqs, t; name = name, systems=[frame_a])
end

@component function TwoBodies(; name, use_w = false)
    @named B1 = Body(use_w = use_w)
    @named B2 = Body(use_w = use_w)
    eqs = Equation[
        connect(B1.frame_a, B2.frame_a)
    ]
    return System(eqs, t; name = name, systems=[B1, B2])
end

@named one_body = Body(use_w = true) # uses angular velocity w_a as state so equations are symbolically simpler and correctly reduced to 12 unknowns
@named one_body_alt = Body(use_w = false) # uses Euler angle rates θ_t as state which leads to the creation of two extra algebraic equations 14 in total to calculate w_a from θ_t
@named two_bodies = TwoBodies(use_w = true) # produces similar results to one_body_alt picking θ_t as state variables for both bodies and creates algebraic equations to calculate w_a 
@named two_bodies_alt = TwoBodies(use_w = false) # reduces to 12 unknowns per body by picking θ_t as state variables directly and since expressions are explicit no extra algebraic equations are created BUT the equations are more symbolically complex

one_body = mtkcompile(one_body; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))
two_bodies = mtkcompile(two_bodies; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))
one_body_alt = mtkcompile(one_body_alt; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))
two_bodies_alt = mtkcompile(two_bodies_alt; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))

```

---

<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 12, 2026, 7:01am UTC](https://discourse.julialang.org/t/modelingtoolkit-compiler-does-not-simplify-as-expected/134978/10 "2026-01-12T07:01:56Z")

</div>

> [@PhysicsBender](#):
>
> Depending on how I formulate the angular velocity kinematics in the model I get different response from the symbolic engine and compiler.

That’s expected, the equations involved are nonlinear and the symbolic machinery cannot solve them by itself.

> [@PhysicsBender](#):
>
> It would be intresting if there was a way to force the compiler to pick certain variables as root/states just like you can do in Modelica.

The Multibody.jl library is offering this feature, if you are impatient and do not want to wait for that you can implement this feature yourself using structural parameters.