两个结构完全相同的ModelingToolkit模型——一个能求解，一个在 `solve` 时崩溃

编程语言 2026-07-09

我有两个使用 ModelingToolkit.jl 的Julia模型，在各方面都完全相同，唯一的差别是一个子组件的名称。尽管如此，structural_simplify 对每个模型产生的未知数数量不同（10 与 12），一个模型成功完成，另一个则失败。我无法理解为什么一个组件的名称会影响结构简化。

唯一的差异

Working (Model_Working.jl) — 第二个摆被命名为 pendulum1:

@named pendulum1 = Pendulum(mass=5.0, length=5.0)

Not working (Model_NotWorking.jl) — 第二个摆被命名为 pendulum2:

@named pendulum2 = Pendulum(mass=5.0, length=5.0)

完整的工作代码（`Model_Working.jl`）

using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots

t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits

const STATESELECT_NEVER   = -100
const STATESELECT_AVOID   = -10
const STATESELECT_DEFAULT = 0
const STATESELECT_PREFER  = 50
const STATESELECT_ALWAYS  = 100

_outerProduct(u, v) = reshape(u, :, 1) .* reshape(v, 1, :)
_skew(e) = [0 -e[3] e[2]; e[3] 0 -e[1]; -e[2] e[1] 0]

@connector function Frame(; name)
    @variables begin
        (r_0(t))[1:3]
        (T(t))[1:3, 1:3]
        (w(t))[1:3]
        (f(t))[1:3], [connect=Flow]
        (t_(t))[1:3], [connect=Flow]
    end
    ODESystem(Equation[], t, [collect(r_0)..., collect(T)..., collect(w)..., collect(f)..., collect(t_)...], []; name=name)
end

function Fixed(; name, r=[0,0,0])
    @parameters (r[1:3]) = r
    @named frame_b = Frame()
    eqs = [
        frame_b.r_0 ~ r,
        frame_b.T   ~ [1 0 0; 0 1 0; 0 0 1],
        frame_b.w   ~ [0,0,0],
    ]
    ODESystem(eqs, t, [], [collect(r)...]; systems=[frame_b], name=name)
end

function FixedTranslation(; name, r=[0,1,0])
    @parameters (r[1:3]) = r
    @named frame_a = Frame()
    @named frame_b = Frame()
    eqs = [
        frame_b.r_0  ~ frame_a.r_0 + transpose(frame_a.T) * r,
        frame_b.T    ~ frame_a.T,
        frame_b.w    ~ frame_a.w,
        zeros(3)     ~ frame_a.f + frame_b.f,
        zeros(3)     ~ frame_a.t_ + frame_b.t_ + cross(r, frame_b.f),
    ]
    ODESystem(eqs, t, [], [collect(r)...]; systems=[frame_a, frame_b], name=name)
end

function Body(; name, g_0=[0,0,-9.81], r_CM=[0,0,0], m=1,
               I_11=0.001, I_22=0.001, I_33=0.001, I_21=0, I_31=0, I_32=0)
    @parameters begin
        (g_0[1:3]) = g_0; (r_CM[1:3]) = r_CM; m = m
        I_11 = I_11; I_22 = I_22; I_33 = I_33
        I_21 = I_21; I_31 = I_31; I_32 = I_32
    end
    I = [I_11 I_21 I_31; I_21 I_22 I_32; I_31 I_32 I_33]
    @variables begin
        (r_0(t))[1:3]; (v_0(t))[1:3]; (a_0(t))[1:3]
        (w_a(t))[1:3]; (z_a(t))[1:3]
    end
    @named frame_a = Frame()
    eqs = [
        r_0 ~ frame_a.r_0,  v_0 ~ D(r_0),  a_0 ~ D(v_0),
        w_a ~ frame_a.w,    z_a ~ D(w_a),
        frame_a.f  ~ m*(frame_a.T*(a_0 - g_0) + cross(z_a, r_CM) + cross(w_a, cross(w_a, r_CM))),
        frame_a.t_ ~ I*z_a + cross(w_a, I*w_a) + cross(r_CM, frame_a.f),
    ]
    ODESystem(eqs, t,
        [collect(r_0)..., collect(v_0)..., collect(a_0)..., collect(w_a)..., collect(z_a)...],
        [collect(g_0)..., collect(r_CM)..., m, I_11, I_22, I_33, I_21, I_31, I_32];
        systems=[frame_a], name=name)
end

@connector function Flange(; name)
    @variables begin; phi(t); tau(t), [connect=Flow]; end
    ODESystem(Equation[], t, [phi, tau], []; name=name)
end

axisAngleToT(e, angle) =
    _outerProduct(e, e) + (Matrix(I, 3, 3) - _outerProduct(e, e))*cos(angle) - _skew(e)*sin(angle)

function Rotational_Fixed(; name, phi0=0)
    @parameters phi0 = phi0
    @named flange = Flange()
    ODESystem([flange.phi ~ phi0], t, [], [phi0]; systems=[flange], name=name)
end

function Revolute(; name, n=[0,0,1], stateSelect=STATESELECT_PREFER)
    e = n/sqrt(dot(n, n))
    @parameters (n[1:3]) = n
    @variables begin
        phi(t), [state_priority = stateSelect];  w(t), [state_priority = stateSelect]
        a(t);  tau(t);  angle(t);  (T_rel(t))[1:3, 1:3]
    end
    @named axis = Flange(); @named support = Flange()
    @named frame_a = Frame(); @named frame_b = Frame()
    @named fixed = Rotational_Fixed()
    eqs = [
        angle ~ phi,  w ~ D(phi),  a ~ D(w),
        frame_b.r_0 ~ frame_a.r_0,
        T_rel       ~ axisAngleToT(e, angle),
        frame_b.T   ~ frame_a.T * T_rel,
        frame_b.w   ~ frame_b.T*transpose(frame_a.T)*(frame_a.w + w * e),
        frame_a.f   ~ -frame_a.T*transpose(frame_b.T)*frame_b.f,
        frame_a.t_  ~ frame_a.T*transpose(frame_b.T)*frame_b.t_,
        tau ~ -dot(frame_b.t_, e),
        phi ~ axis.phi,  tau ~ axis.tau,
        connect(fixed.flange, support),
    ]
    ODESystem(eqs, t, [phi, w, a, tau, angle, collect(T_rel)...], [collect(n)...];
              systems=[axis, support, frame_a, frame_b, fixed], name=name)
end

function Pendulum(; name, mass=1, length=1)
    @named frame = Frame()
    @named revolute         = Revolute(n=[1,0,0])
    @named fixedTranslation = FixedTranslation(r=[0,length,0])
    @named body             = Body(m=mass)
    eqs = [
        connect(frame, revolute.frame_a),
        connect(revolute.frame_b, fixedTranslation.frame_a),
        connect(fixedTranslation.frame_b, body.frame_a),
    ]
    ODESystem(eqs, t, [], []; systems=[frame, revolute, fixedTranslation, body], name=name)
end

function MyModel(; name)
    @named fixed     = Fixed(r=[0,0,0])
    @named pendulum  = Pendulum(mass=2.5, length=1.0)
    @named pendulum1 = Pendulum(mass=5.0, length=5.0)   # <-- named pendulum1
    eqs = [
        connect(fixed.frame_b, pendulum.frame),
        connect(pendulum1.frame, fixed.frame_b),
    ]
    ODESystem(eqs, t, [], []; systems=[fixed, pendulum, pendulum1], name=name)
end

@named _sys_raw = MyModel()
_sys = structural_simplify(_sys_raw)
println("Structural simplify done. n_unknowns = ", length(unknowns(_sys)))

_u0      = [v => 0.0 for v in unknowns(_sys)]
_tspan   = (0.0, 10.0)
_guesses = Dict{Any,Float64}(v => 0.0 for v in unknowns(_sys))
_prob    = ODEProblem(_sys, _u0, _tspan, []; guesses=_guesses, build_initializeprob=false)
_sol     = solve(_prob, Rodas5(); reltol=1e-6, abstol=1e-8)

完整的不可工作代码（`Model_NotWorking.jl`）

这段代码是完全相同，除了 MyModel 内的这两行：

@named pendulum2 = Pendulum(mass=5.0, length=5.0)   # <-- named pendulum2
# ...
ODESystem(eqs, t, [], []; systems=[fixed, pendulum, pendulum2], name=name)

为什么将 pendulum1 → pendulum2 重命名会把 structural_simplify 的未知数从10个变成12个？