state_priority not working (object-oriented planar pendulum models)

[AndreasHofmann217]

Describe the bug 🐞

Coming from Modelica, we are currently writing some PlanarMechanics library for MTK.
We want to enforce states, since there is no dynamic state selection that manages systems becoming singular during simulation with an arbitrary valid set of selected states.
However, using the state_priority variable metadata does not have any effect applied to our models for the revolute joints in the system.
Using state_priority for the angle of Body at the end of the pendulum does work, but limits applicability for general purpose.

Expected behavior

Providing state_priority should hint the compiler to select states. For our pendulum examples this would be angles and angular velocities of the joints.

Minimal Reproducible Example 👇

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

# BASE COMPONENTS
@connector Frame begin
    x(t), [description = "x position resolved in world frame"]
    y(t), [description = "y position resolved in world frame"]
    phi(t), [description = "rotation angle resolved in world frame"]
    fx(t), [connect = Flow, description = "x component of force resolved in world frame"]
    fy(t), [connect = Flow, description = "y component of force resolved in world frame"]
    tau(t), [connect = Flow, description = "torque resolved in world frame"]
end

@mtkmodel Fixed begin
    @parameters begin
        r[1:2] = [0.0,0.0], [description = "Position resolved in world frame"]
        phi = 0.0, [description = "Fixed Angle resolved in world frame"]
    end
    @components begin
        frame_a = Frame()
    end
    @equations begin
        # position and rotation
        frame_a.x ~ r[1]
        frame_a.y ~ r[2]
        frame_a.phi ~ phi
    end
end

@mtkmodel FixedTranslation begin
    @parameters begin
        r[1:2] = [1.0,0.0], [description = "Fixed x,-y length of the rod resolved w.r.t. frame_a"]
    end
    @variables begin
        r0x(t), [description = "x component of the vector from frame_a to frame_b resolved in world frame"]
        r0y(t), [description = "y component of the vector from frame_a to frame_b resolved in world frame"]
    end
    @components begin
        frame_a = Frame()
        frame_b = Frame()
    end
    @equations begin
        # distance vector resolved in world frame
        r0x ~ cos(frame_a.phi) * r[1] - sin(frame_a.phi) * r[2]
        r0y ~ sin(frame_a.phi) * r[1] + cos(frame_a.phi) * r[2]
        # position and rotation
        frame_a.x + r0x ~ frame_b.x
        frame_a.y + r0y ~ frame_b.y
        frame_a.phi ~ frame_b.phi
        # forces
        frame_a.fx + frame_b.fx ~ 0.0
        frame_a.fy + frame_b.fy ~ 0.0
        frame_a.tau + frame_b.tau + r0x * frame_b.fy - r0y * frame_b.fx ~ 0.0
    end
end

@mtkmodel Body begin
    @parameters begin
        m = 1.0, [description = "Mass of Body"]
        J = 1.0, [description = "Inertia of Body"]
        g[1:2] = [0.0,-9.81], [description = "gravitational acceleration"]
    end
    @variables begin
        r1(t), [state_priority = -100.0, description = "x position resolved in world frame"]
        r2(t), [state_priority = -100.0, description = "y position resolved in world frame"]
        v1(t), [state_priority = -100.0, description = "x velocity resolved in world frame"]
        v2(t), [state_priority = -100.0, description = "y velocity resolved in world frame"]
        a1(t), [description = "x acceleration resolved in world frame"]
        a2(t), [description = "y acceleration resolved in world frame"]
        phi(t), [state_priority = 0.0, description = "rotation angle resolved in world frame"]
        w(t), [state_priority = 0.0, description = "angular velocity resolved in world frame"]
        z(t), [description = "angular acceleration resolved in world frame"]
    end
    @components begin
        frame_a = Frame()
    end
    @equations begin
        # derivatives
        D(r1) ~ v1
        D(r2) ~ v2
        D(v1) ~ a1
        D(v2) ~ a2
        D(phi) ~ w
        D(w) ~ z
        # position and rotation
        r1 ~ frame_a.x
        r2 ~ frame_a.y
        phi ~ frame_a.phi
        # Newton
        m * a1 ~ frame_a.fx + m * g[1]
        m * a2 ~ frame_a.fy + m * g[2]
        frame_a.tau ~ J * z
    end
end

@mtkmodel RevoluteJointPassive begin
    @variables begin
        phi(t), [state_priority=100.0, description="Relative angle between frame_a and frame_b"]
        w(t), [state_priority=100.0, description="Angular velocity between frame_a and frame_b"]
        z(t), [description="Angular acceleration between frame_a and frame_b"]
    end
    @components begin
        frame_a = Frame()
        frame_b = Frame()
    end
    @equations begin
        # derivatives
        D(phi) ~ w
        D(w) ~ z
        # position and rotation
        frame_a.phi + phi ~ frame_b.phi
        frame_a.x ~ frame_b.x
        frame_a.y ~ frame_b.y
        # forces
        frame_a.fx + frame_b.fx ~ 0.0
        frame_a.fy + frame_b.fy ~ 0.0
        frame_a.tau ~ 0.0
        frame_a.tau + frame_b.tau ~ 0.0
    end
end

# SINGLE PENDULUM
@mtkmodel Pendulum begin
    @components begin
        fixed = Fixed()
        pivot1 = RevoluteJointPassive()
        rod1 = FixedTranslation()
        body = Body()   
    end
    @equations begin
        connect(fixed.frame_a, pivot1.frame_a)
        connect(pivot1.frame_b, rod1.frame_a)
        connect(rod1.frame_b, body.frame_a)
    end
end

@mtkcompile pendulum = Pendulum()
unknowns(pendulum)
equations(pendulum)


# DOUBLE PENDULUM MODEL
@mtkmodel DoublePendulum begin
    @components begin
        fixed = Fixed()
        pivot1 = RevoluteJointPassive()
        pivot2 = RevoluteJointPassive()
        rod1 = FixedTranslation()
        rod2 = FixedTranslation()
        body = Body()   
    end
    @equations begin
        connect(fixed.frame_a, pivot1.frame_a)
        connect(pivot1.frame_b, rod1.frame_a)
        connect(rod1.frame_b, pivot2.frame_a)
        connect(pivot2.frame_b, rod2.frame_a)
        connect(rod2.frame_b, body.frame_a)
    end
end

@mtkcompile doublependulum = DoublePendulum()
unknowns(doublependulum)
equations(doublependulum)

# DOUBLE PENDULUM MODEL
# does work if phi in body has a state_priority
@mtkmodel DoublePendulumTwoBody begin
    @components begin
        fixed = Fixed()
        pivot1 = RevoluteJointPassive()
        pivot2 = RevoluteJointPassive()
        rod1 = FixedTranslation()
        rod2 = FixedTranslation()
        body1 = Body()
        body2 = Body()    
    end
    @equations begin
        connect(fixed.frame_a, pivot1.frame_a)
        connect(pivot1.frame_b, rod1.frame_a)
        connect(rod1.frame_b, body1.frame_a)
        connect(rod1.frame_b, pivot2.frame_a)
        connect(pivot2.frame_b, rod2.frame_a)
        connect(rod2.frame_b, body2.frame_a)
    end
end

@mtkcompile doublependulumtwobody = DoublePendulumTwoBody()
unknowns(doublependulumtwobody)
equations(doublependulumtwobody)

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

(StateSelection) pkg> status
Status `C:\JuliaIssues\StateSelection\Project.toml`
⌃ [961ee093] ModelingToolkit v10.21.0
Info Packages marked with ⌃ have new versions available and may be upgradable.

Julia Version 1.11.6
Commit 9615af0f26 (2025-07-09 12:58 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 16 × 11th Gen Intel(R) Core(TM) i7-11850H @ 2.50GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, tigerlake)
Threads: 8 default, 0 interactive, 4 GC (on 16 virtual cores)
Environment:
  JULIA_NUM_THREADS = 8
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1
  JULIA_PKG_SERVER = https://juliahub.com