Inline linear SCC emits an inconsistent runtime system (rank 14/15, ‖Ax−b‖/‖b‖ ≈ 1) at a consistent state — integration goes Unstable

[baggepinnen]

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 b is mis-accounted. __reduce_linear_system!'s handling of eliminated equations (reworked in fix: handle eliminated equations in inline linear SCCs appropriately #94) and the constants/aliases topological 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).

[baggepinnen]

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\b emitted with b ∉ range(A) at runtime), but the fix needs to address the splitting, not the assembly:

1. solve runtime-rank-deficient block families jointly (the sparse joint-solve direction of Inline linear SCCs: 296×296 dense runtime solve for an ~18-unknown core; non-rank-tolerant; symbolic shrinkage intractable — solve sparse+rank-tolerant instead #95 handles this naturally), or
2. consult maybe_zeros at 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.

[baggepinnen]

Closing the loop with the physical identification. The runtime probe is now committed as test/probe_halfcar_blocks.jl in JuliaComputing/MultibodyComponents.jl#67: it harvests every emitted getindex(\(A, b), i) from equations/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 of A/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 are wheel_rotation.R_relˍt row 3, the contact-frame derivatives (auxˍt, e_long_0ˍt, delta_0ˍt), and wheel_rotation.phiˍt. The right null vector is pure phiˍt: for the spin revolute with n = [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. The b ∉ 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-bar r2.R_relˍt/phiˍt chain. 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 = 100 on the four-bar revolute r2 and on wheel_rotation (effective with #103), plus angular_state = false wheels. 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 to Success (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_zeros on 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.