Repository navigation
Inline linear SCC emits an inconsistent runtime system (rank 14/15, ‖Ax−b‖/‖b‖ ≈ 1) at a consistent state — integration goes Unstable #98
Description
Summary
The inline-linear-SCC pass emits runtime A \ b terms under the implicit contract that the linear system is solvable at runtime. On MultibodyComponents.examples.suspension.HalfCar (22 states, emitted blocks [13, 15, 12, 13, 15, 12, 295]), two of the small blocks are rank-deficient with b outside the range of A — measured at a fully consistent initialization point (initialization residual ~4e-10, found with LM + pivoted QR):
# min-norm pivoted-QR solve of the emitted (A, b) at runtime:
n=15 rank=14 relres=0.945 # ‖A x − b‖ / ‖b‖ ≈ 1: b ⟂ range(A)
n=13 rank=12 relres=0.952 # same block family, later during integration
n=13 rank=12 relres=0.0 # rank-deficient but CONSISTENT — harmless
n=295 rank=295 relres=1.4e-14 # the big block (#95) is full-rank and fine here
With the stock INLINE_LINEAR_SCC_OP = (\) this is a SingularException or, via the DyadCompilerPasses LinearSolve path (JuliaComputing/DyadCompilerPasses.jl#41), garbage that silently feeds the ODE RHS. Either way the torn variables of that block do not satisfy their equations, so the dynamics are wrong; we observe immediate Unstable (dt → 1e-21) on integration from a consistent state, with a smooth RHS and a clean FD Jacobian (max|J| ≈ 1.7e5, constant in the FD step size).
Why this looks like a construction bug rather than a model property
- The state is on the constraint manifold to ~1e-10, so a correctly-assembled linearization of consistent equations should be consistent to roundoff, not to O(1).
- A 100% relative residual on a 15×15 block means one redundant row's RHS contradicts the others — the shape you'd get if an eliminated row's contribution to
bis mis-accounted.__reduce_linear_system!'s handling of eliminated equations (reworked in fix: handle eliminated equations in inline linear SCCs appropriately #94) and theconstants/aliasestopological update would be my first suspects. - The sibling 13×13 block with the same rank deficiency is consistent (
relres = 0), so rank deficiency per se (static reaction indeterminacy) is expected and fine; the inconsistency is the anomaly.
Reproduction
HalfCar from MultibodyComponents.jl with maybe_zeros set on the joint-axis parameter vectors (16 arrays), optimize=:basic, and a rank-revealing override of the runtime linear solve to measure rank/relres (probe script available on request — it is the same setup as the measurements posted on #97 and #95). Versions: StateSelection main @ 3d10ed7 (also with #97 merged — identical numbers), ModelingToolkitTearing v1.14.1, MTK master v11.x.
Related: #95 (block size), #97 (measured there), JuliaComputing/DyadCompilerPasses.jl#41 (what the runtime does when the solve fails).
Activity
baggepinnen commentedon Jun 11, 2026
Update — root-cause correction based on testing #100's diagnostics on HalfCar (details in this comment):
The inconsistency is not a block-construction bug. With #100's (fixed) self-check, every emitted block — including the 396 → 295 reduction — is exact and consistent at generic symbol values (relres ≤ 1e-13), and the suspected buried-variable repair never fires.
The actual mechanism is cross-block indeterminacy at the degenerate parameter point: the SCC family is split into sequentially-solved blocks (13, 15, 12, 295 here). At the true parameter values (joint-axis components exactly 0) the upstream 13×13 block drops to rank 12 with a 1-dimensional nullspace; whatever solution is chosen for it (min-norm gauge, or garbage LU in production) is then incompatible with any solution of the downstream 15×15 block —
relres = 0.945— even though the union of the equations is satisfiable (the state satisfies the model to 4e-10, and the corner whose upstream block stays full-rank has a consistent downstream block, relres 2e-16).So the title condition stands (
A\bemitted withb ∉ range(A)at runtime), but the fix needs to address the splitting, not the assembly:maybe_zerosat emission time: the parameters whose zeros cause the rank drop are exactly the markable ones, so the pass can refuse to split SCC families whose runtime rank can drop.feat: jointly solve coupled rank-deficible inline-linear SCC families (…
feat: peel-on-failure for jointly-solved families; robust self-check …
rework family grouping: set-dependent linearity, greedy convex growth…
baggepinnen commentedon Jun 12, 2026
Closing the loop with the physical identification. The runtime probe is now committed as
test/probe_halfcar_blocks.jlin JuliaComputing/MultibodyComponents.jl#67: it harvests every emittedgetindex(\(A, b), i)fromequations/observed, evaluates each(A, b)at the real parameter/init point through the compiled observed chain (getu), and prints SVD rank, min-norm relres, and the null directions ofA/Aᵀmapped back to named variables.The rank-deficient blocks are not reaction-force blocks — they are velocity-level dummy-derivative blocks over rotation-matrix derivative chains, and the deficiency is parameter-degenerate: full rank at generic probe points (which is why the compile-time self-check sees nothing), singular at the actual axis-aligned joint axes. (Block sizes shifted slightly on the current dev stack —
[12,12,13,13,15,15,~296]— but the families are recognizably the same.)n=15(both suspensions; rank 14, relres ≈ 0.93–0.95): torn variables arewheel_rotation.R_relˍtrow 3, the contact-frame derivatives (auxˍt,e_long_0ˍt,delta_0ˍt), andwheel_rotation.phiˍt. The right null vector is purephiˍt: for the spin revolute withn = [0,0,±1], row 3 of the Rodrigues matrix has no φ-dependence, so the φ̇ column is identically zero — symbolically incident, numerically zero at the real axis. Theb ∉ range(A)part is gauge contamination from the upstream singular solve, consistent with the 2e16a5f analysis on the feat: jointly solve coupled rank-deficible inline-linear SCC families (#98) #100 branch.n=13(rank 12, relres = 0, mirrored suspension only — the right side's same-size block is full rank, sv ratio ~0.01): the four-barr2.R_relˍt/phiˍtchain. The two sides got different tear pivot sets (the #4612-style order sensitivity), and the mirrored one is singular at the real parameters.So the emitted systems are genuinely rank-deficient, but the indeterminacy is "which rotation-matrix derivative parametrizes the joint rate", not static reactions. The model-level fix is the textbook one — make the joint coordinates the states:
statePriority = 100on the four-bar revoluter2and onwheel_rotation(effective with #103), plusangular_state = falsewheels. With that on the #103+#99+#100 stack: every emitted block is full rank at the real point (the two former n=15 families merge into two full-rank n=46 blocks, relres ≈ 1e-14), the LinearSolve solver-failure flood disappears, and HalfCar integrates(0, 2)s toSuccess(FullCar(0, 1)s likewise). Fixed in JuliaComputing/MultibodyComponents.jl#67 (65febda).What remains compiler-side for this issue, I think, is whether the inline-linear pass should detect parameter-degenerate pivots (e.g. honoring
maybe_zeroson the axis vectors — modulo the SciML/ModelingToolkit.jl#4613 key-matching trap) or solve rank-tolerantly (#95). With the model fixed, neither is load-bearing for the suspension models anymore.test: directly exercise the buried-SCC-var repair (#98)
refactor: demote buried-b repair to a self-check assertion (#98)