GeezerVison is showing github.com — `interconnect` reports "algebraic loop detected" for MIMO systems with feedthrough that have no (or a well-posed) algebraic loop · Issue #1257 · python-control/python-control
interconnect reports "algebraic loop detected" for MIMO systems with feedthrough that have no (or a well-posed) algebraic loop #1257
DescriptionSummary
ct.interconnect raises RuntimeError: algebraic loop detected in two situations where the interconnection is fine:
- 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.
- 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
-
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.
-
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.
-
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). MetadataMetadataDevelopmentNo branches or pull requestsIssue actions
|