Summary

ct.interconnect raises RuntimeError: algebraic loop detected in two situations where the interconnection is fine:

1. False positive (no algebraic loop at all). A feedthrough path that visits the same MIMO block more than once needs more passes than _compute_static_io allows.
2. Well-posed algebraic loop. I - D M is invertible, so the closed loop is well defined. ct.feedback and Simulink both handle it, but the fixed-point iteration never resolves it.

I hit (1) with decentralized PI control of a 2x2 plant.

Reproducer (python-control 0.10.2)

import numpy as np
import control as ct

# 2x2 plant with direct feedthrough only in g11: D = [[1, 0], [0, 0]]
G = ct.combine_tf([[ct.tf([1, 1], [1, 2]), ct.tf(1, [5, 9])],
                   [ct.tf(4, [1, 2]), ct.tf([1, 1], [1, 6, 9])]])
pi = ct.tf([1, 1], [1, 0])                   # PI controller, D = 1

# Case 1: no algebraic loop (D_G @ D_K is nilpotent), but interconnect fails
plant = ct.ss(G, inputs=['u[0]', 'u[1]'], outputs=['y[0]', 'y[1]'],
              name='plant')
err = ct.summing_junction(inputs=['r', '-y'], outputs='e', dimension=2,
                          name='err')
pi1 = ct.ss(pi, inputs='e[0]', outputs='u[1]', name='pi1')  # off-diagonal
pi2 = ct.ss(pi, inputs='e[1]', outputs='u[0]', name='pi2')  # pairing
try:
    ct.interconnect([err, pi1, pi2, plant], inplist=['r'], outlist=['y'])
except RuntimeError as e:
    print('case 1:', e)

K = ct.ss(ct.combine_tf([[0, pi], [pi, 0]]))
print('case 1 feedback poles:',
      np.sort_complex(ct.feedback(ct.ss(G) * K, np.eye(2)).poles()))

# Case 2: genuine but well-posed algebraic loop (static decoupler T)
T = np.array([[4/3, -2/27], [-6, 4/3]])
plantT = ct.ss(G * T, inputs=['u[0]', 'u[1]'], outputs=['y[0]', 'y[1]'],
               name='plant')
DL = ct.ss(G * T).D @ K.D
print('case 2 det(I + D_G D_K) =', np.linalg.det(np.eye(2) + DL))
try:
    ct.interconnect([err, pi1, pi2, plantT], inplist=['r'], outlist=['y'])
except RuntimeError as e:
    print('case 2:', e)
print('case 2 feedback poles:',
      np.sort_complex(ct.feedback(ct.ss(G * T) * K, np.eye(2)).poles()))

Output (poles rounded):

case 1: algebraic loop detected
case 1 feedback poles: [-3.942-1.829j -3.942+1.829j -2.366 -1.903 -0.764 -0.083]
case 2 det(I + D_G D_K) = 0.9259259259259259
case 2: algebraic loop detected
case 2 feedback poles: [-3.000 -3.000 -2.194 -2.000 -1.894 -1.800 -1.460-4.365j -1.460+4.365j -0.811 -0.109]

Both closed loops are well defined and stable (ct.feedback gives the poles above).

Analysis

Case 1. InterconnectedSystem._compute_static_io iterates at most len(self.syslist) + 1 times. Each pass moves a signal one block further. That bound assumes a feedthrough path visits each subsystem at most once, which is not true for MIMO blocks. Here the path is

y[1] -> err -> e[1] -> pi2 -> u[0] -> plant -> y[0] -> err -> e[0] -> pi1 -> u[1] -> plant

It visits err and plant twice: 6 propagation passes plus 1 to confirm convergence, but only 5 are allowed. The diagram is free of algebraic loops (D_G @ D_K = [[0, 1], [0, 0]] is nilpotent).

Splitting the plant into four SISO blocks plus summing junctions makes case 1 pass, but only because the extra blocks raise the pass limit.

Case 2. y[0] depends directly on u[1] = pi1(r[0] - y[0]), so this is a genuine algebraic loop, but det(I + D_G D_K) != 0. Plain fixed-point iteration cannot resolve it in finitely many passes with an exact-equality stopping test. It also may not converge at all when the loop gain is larger than one.

Possible fixes

1. Minimal: change the bound (nlsys.py#L957). Without an algebraic loop, the longest feedthrough path is bounded by the number of scalar internal signals, not the number of subsystems:

   cycle_count = ninputs + 1        # instead of len(self.syslist) + 1

   Monkey-patching this into 0.10.2 makes case 1 work, and the resulting A has the same eigenvalues as ct.feedback. Updating ulist within a pass (Gauss-Seidel style) would also reduce the number of passes.

2. Structural detection. Build the signal dependency graph from the feedthrough structure (nonzero pattern of D for linear blocks; all-to-all, or user-declared, for nonlinear blocks). Then evaluate the blocks in topological order in a single pass. If cycles exist, report the signals in each strongly connected component, which is far more informative than the current error message.

3. All-linear interconnections: solve exactly. At the moment LinearICSystem obtains the closed loop by numerically linearizing the interconnected system (self.linearize(0, 0), finite differences). With block-diagonal A, B, C, D, internal connection map M and input map B_in:

   y = (I - D M)^{-1} (C x + D B_in r)
   A_cl = A + B M (I - D M)^{-1} C
   

   with analogous expressions for B_cl, C_cl, D_cl. This is exact, avoids finite-difference error, and also handles well-posed algebraic loops (case 2), as feedback does. An error is needed only when I - D M is singular, i.e. the interconnection is ill-posed. For nonlinear systems with genuine loops, a Newton solve on the loop variables (as in Simulink's algebraic loop solver) would be the analogue.

Related: #461 (different cause: a non-deterministic output function).