Initialization system fails to structural simplify

[bradcarman]

Describe the bug 🐞

Initialization system fails for the given example. The same system is able to solve in Modelica.

Expected behavior

Initialization system should be able to run structural_simplify

Minimal Reproducible Example 👇

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

function System(;name)

    pars = @parameters begin
        m = 1.0
        L = 0.5
        g = 9.81
    end

    vars = @variables begin
        x1(t) = L
        y1(t) = 0
        dx1(t), [guess = 0]
        dy1(t), [guess = 0]
        ddx1(t), [guess = 0]
        ddy1(t), [guess = -g]
        F1(t), [guess = 0]
        θ1(t), [state_priority=10, guess = 0]
        dθ1(t), [guess = 0]
        ddθ1(t), [guess = 0]

        x2(t) = 2L
        y2(t) = 0
        dx2(t), [guess = 0]
        dy2(t), [guess = 0]
        ddx2(t), [guess = 0]
        ddy2(t), [guess = -g]
        F2(t), [guess = 0]
        θ2(t), [state_priority=10, guess = 0]
        dθ2(t), [guess = 0]
        ddθ2(t), [guess = 0]

        x3(t) = 3L
        y3(t) = 0
        dx3(t), [guess = 0]
        dy3(t), [guess = 0]
        ddx3(t), [guess = 0]
        ddy3(t), [guess = -g]
        F3(t), [guess  = 0]
        θ3(t), [state_priority=10, guess = 0]
        dθ3(t), [guess = 0]
        ddθ3(t), [guess = 0]
    end

    eqs = [D(x1) ~ dx1
        D(dx1) ~ ddx1
        D(y1) ~ dy1
        D(dy1) ~ ddy1
        D(θ1) ~ dθ1
        D(dθ1) ~ ddθ1
        cos(θ1) ~ +x1 / L
        sin(θ1) ~ -y1 / L
        m * ddx1 ~ -F1 * cos(θ1) + F2 * cos(θ2)
        m * ddy1 ~ +F1 * sin(θ1) - F2 * sin(θ2) - m * g

        D(x2) ~ dx2
        D(dx2) ~ ddx2
        D(y2) ~ dy2
        D(dy2) ~ ddy2
        D(θ2) ~ dθ2
        D(dθ2) ~ ddθ2
        cos(θ2) ~ +(x2 - x1) / L
        sin(θ2) ~ -(y2 - y1) / L
        m * ddx2 ~ -F2 * cos(θ2) + F3 * cos(θ3)
        m * ddy2 ~ +F2 * sin(θ2) - F3 * sin(θ3) - m * g
        
        D(x3) ~ dx3
        D(dx3) ~ ddx3
        D(y3) ~ dy3
        D(dy3) ~ ddy3
        D(θ3) ~ dθ3
        D(dθ3) ~ ddθ3
        cos(θ3) ~ +(x3 - x2) / L
        sin(θ3) ~ -(y3 - y2) / L
        m * ddx3 ~ -F3 * cos(θ3)
        m * ddy3 ~ +F3 * sin(θ3) - m * g
        ]
    return ODESystem(eqs, t, vars, pars; name)
end

@mtkbuild sys = System()

initsys = ModelingToolkit.generate_initializesystem(sys) #54 equations/54 unknowns
initsys = structural_simplify(initsys) # ERROR! InvalidSystemException: The system is structurally singular!

Error & Stacktrace ⚠️

Stacktrace:
  [1] check_consistency(state::TearingState{NonlinearSystem}, orig_inputs::Set{Any})
    @ ModelingToolkit.StructuralTransformations C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\structural_transformation\utils.jl:108
  [2] _structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, dummy_derivative::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\systems\systemstructure.jl:691
  [3] _structural_simplify!
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\systems\systemstructure.jl:674 [inlined]
  [4] #structural_simplify!#1096
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\systems\systemstructure.jl:667 [inlined]
  [5] __structural_simplify(sys::NonlinearSystem, io::Nothing; simplify::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\systems\systems.jl:83
  [6] __structural_simplify
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\systems\systems.jl:64 [inlined]
  [7] structural_simplify(sys::NonlinearSystem, io::Nothing; simplify::Bool, split::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\systems\systems.jl:24
  [8] structural_simplify
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\systems\systems.jl:21 [inlined]
  [9] structural_simplify(sys::NonlinearSystem)
    @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\U2JaS\src\systems\systems.jl:21
 [10] top-level scope
    @ c:\Work\Packages\Issue953\mwe.jl:86

Environment (please complete the following information):
Using ModelingToolkit v9.24.0

Additional context
Here is the equivalent Modelica model...

model MTK
	parameter Real m = 1.0;
	parameter Real L = 0.5;
	parameter Real g = 9.81;
	Real x1(start = 0.5, fixed=true);
	Real y1(start = 0, fixed=true);
	Real dx1(start = 0);
	Real dy1(start = 0);
	Real ddx1(start = 0);
	Real ddy1(start = -9.81);
	Real F1(start = 0);
	Real theta1(start = 0);
	Real dtheta1(start = 0);
	Real ddtheta1(start = 0);
	Real x2(start = 1.0, fixed=true);
	Real y2(start = 0, fixed=true);
	Real dx2(start = 0);
	Real dy2(start = 0);
	Real ddx2(start = 0);
	Real ddy2(start = -9.81);
	Real F2(start = 0);
	Real theta2(start = 0);
	Real dtheta2(start = 0);
	Real ddtheta2(start = 0);
	Real x3(start = 1.5, fixed=true);
	Real y3(start = 0, fixed=true);
	Real dx3(start = 0);
	Real dy3(start = 0);
	Real ddx3(start = 0);
	Real ddy3(start = -9.81);
	Real F3(start = 0);
	Real theta3(start = 0);
	Real dtheta3(start = 0);
	Real ddtheta3(start = 0);
equation
	der(x1) = dx1;
	der(dx1) = ddx1;
	der(y1) = dy1;
	der(dy1) = ddy1;
	der(theta1) = dtheta1;
	der(dtheta1) = ddtheta1;
	cos(theta1) = x1 / L;
	sin(theta1) = (-y1) / L;
	m*ddx1 = cos(theta2)*F2 - cos(theta1)*F1;
	m*ddy1 = -g*m + sin(theta1)*F1 - sin(theta2)*F2;
	der(x2) = dx2;
	der(dx2) = ddx2;
	der(y2) = dy2;
	der(dy2) = ddy2;
	der(theta2) = dtheta2;
	der(dtheta2) = ddtheta2;
	cos(theta2) = (-x1 + x2) / L;
	sin(theta2) = (y1 - y2) / L;
	m*ddx2 = -cos(theta2)*F2 + F3*cos(theta3);
	m*ddy2 = -g*m - sin(theta3)*F3 + sin(theta2)*F2;
	der(x3) = dx3;
	der(dx3) = ddx3;
	der(y3) = dy3;
	der(dy3) = ddy3;
	der(theta3) = dtheta3;
	der(dtheta3) = ddtheta3;
	cos(theta3) = (x3 - x2) / L;
	sin(theta3) = (-y3 + y2) / L;
	m*ddx3 = -F3*cos(theta3);
	m*ddy3 = -g*m + sin(theta3)*F3;
end MTK;