steadystate() cannot solve a system whose equilibrium is set by a loop rather than by any single block.

Minimal case, a PID around a nonlinear plant:

ref, err = Constant(1.0), Adder("+-")
ctrl = PID(Kp=1.0, Ki=2.0, Kd=0.0, f_max=50)
plant = DynamicalSystem(
    func_dyn=lambda x, u, t: -x**2 + u,
    func_alg=lambda x, u, t: x,
    initial_value=1.0,
    jac_dyn=lambda x, u, t: -2*x,
    )

Sim = Simulation(
    blocks=[ref, err, ctrl, plant, Scope()],
    connections=[
        Connection(ref, err[0]), Connection(plant, err[1]),
        Connection(err, ctrl), Connection(ctrl, plant),
        Connection(plant, Scope()),
        ],
    dt=0.01,
    )

Sim.steadystate()   # RuntimeError: did not converge

The equilibrium exists and is unique: the integral action forces err = 0, so x = 1, u = x**2 = 1, and the PID integrator state is 0.5. The solve stalls at x = 0.618, err = 0.382, with the PID integrator never leaving 0.

Why

_solve is a Gauss-Seidel sweep: _update(t), then each block calls its own solve(t, dt) on its own equations with its own Jacobian. Nothing in that loop can express a condition that couples one block's state to another block's residual.

The PID's integrator has dx/dt = u, so its own residual carries no information about its own state. It is determined entirely by the surrounding loop driving u to zero, which the block-local solve cannot see. The Newton matrix is singular in exactly that direction (#241), so after that fix the state receives no correction there, and a state that never moves also contributes a zero column to the Anderson history, which then cannot move it either.

Substituting a plain fixed-point step x <- g on that direction was tried while fixing #241 and is worse: on a closed loop it injects the full residual every sweep and the iteration diverges (residual 1e11) instead of stalling.

Not a regression

#241 changed the failure from LinAlgError: SVD did not converge in Linear Least Squares, raised from an unrelated place, into a truthful non-convergence report with per-block residuals. The underlying inability to solve this class of system predates it.

Direction

The fix is not in the block-local solve. It needs an outer stage over the coupling: a small global solve over the few quantities the loop determines, with the existing per-block steadystate as the inner iteration. That is the same shape periodic_steady_state already uses (outer Anderson on the period map, inner transient run), and the same structure an inverse solve needs.

Worth noting that this is also what pushed the trim() proposal in #237 towards a separate solver stack. A single mechanism would cover both.