GeezerVison is showing github.com — Add scipy fallbacks for generalized Lyapunov and discrete Sylvester e… · python-control/python-control@4608123
File treeExpand file treeCollapse file tree Open diff view settings Expand file treeCollapse file tree Open diff view settings Collapse file Copy file name to clipboardExpand all lines: control/mateqn.py +129-8Lines changed: 129 additions & 8 deletions | Original file line number | Diff line number | Diff line change |
|---|
@@ -47,15 +47,33 @@ def sb03md(n, C, A, U, dico, job='X', fact='N', trana='N', ldwork=None): | 47 | 47 | try: | 48 | 48 | from slycot import sb04qd | 49 | 49 | except ImportError: | 50 | | - sb0qmd = None | | 50 | + sb04qd = None | 51 | 51 |
| 52 | 52 | try: | 53 | 53 | from slycot import sg03ad | 54 | 54 | except ImportError: | 55 | | - sb04ad = None | | 55 | + sg03ad = None | 56 | 56 |
| 57 | 57 | __all__ = ['lyap', 'dlyap', 'dare', 'care'] | 58 | 58 |
| | 59 | + | | 60 | +def _warn_ill_conditioned_E(E): | | 61 | + """Warn that the inv(E) congruence transform will lose accuracy. | | 62 | +
| | 63 | + The scipy generalized-Lyapunov fallback reduces the problem to a | | 64 | + standard Lyapunov equation by inverting E, so a poorly conditioned E | | 65 | + costs accuracy (continuous and discrete paths alike, regardless of | | 66 | + whether the underlying scipy solve happens to warn). SLICOT sg03ad | | 67 | + (method='slycot') avoids inverting E and is preferable in that case. | | 68 | + """ | | 69 | + condE = np.linalg.cond(E) | | 70 | + if condE > 1.0 / np.sqrt(finfo(float).eps): | | 71 | + warnings.warn( | | 72 | + f"E is ill-conditioned (cond(E) = {condE:.2g}); the " | | 73 | + "method='scipy' generalized Lyapunov solution may have reduced " | | 74 | + "accuracy. Use method='slycot' (SLICOT sg03ad) for a more " | | 75 | + "robust solution.", UserWarning, stacklevel=3) | | 76 | + | 59 | 77 | # | 60 | 78 | # Lyapunov equation solvers lyap and dlyap | 61 | 79 | # | @@ -103,6 +121,21 @@ def lyap(A, Q, C=None, E=None, method=None): | 103 | 121 | X : 2D array | 104 | 122 | Solution to the Lyapunov or Sylvester equation. | 105 | 123 |
| | 124 | + Notes | | 125 | + ----- | | 126 | + For the generalized Lyapunov equation, method='slycot' uses the | | 127 | + SLICOT routine SG03AD, based on the generalized Schur method of | | 128 | + Penzl [1]_, which also handles singular E. With method='scipy', the | | 129 | + equation is transformed to a standard Lyapunov equation by inverting | | 130 | + E, which requires E to be nonsingular and loses accuracy when E is | | 131 | + ill-conditioned (a UserWarning is then issued); method='slycot' does | | 132 | + not invert E and is preferable in that case. | | 133 | +
| | 134 | + References | | 135 | + ---------- | | 136 | + .. [1] Penzl, T., "Numerical solution of generalized Lyapunov | | 137 | + equations", Advances in Computational Mathematics, 8:33-48, 1998. | | 138 | +
| 106 | 139 | """ | 107 | 140 | # Decide what method to use | 108 | 141 | method = _slycot_or_scipy(method) | @@ -162,8 +195,23 @@ def lyap(A, Q, C=None, E=None, method=None): | 162 | 195 | _check_shape(E, n, n, square=True, name="E") | 163 | 196 |
| 164 | 197 | if method == 'scipy': | 165 | | - raise ControlArgument( | 166 | | - "method='scipy' not valid for generalized Lyapunov equation") | | 198 | + # Transform to a standard Lyapunov equation by multiplying | | 199 | + # from the left by inv(E) and from the right by inv(E).T: | | 200 | + # | | 201 | + # (E^-1 A) X + X (E^-1 A)^T + E^-1 Q E^-T = 0 | | 202 | + # | | 203 | + # This requires E to be nonsingular; the SLICOT routine | | 204 | + # SG03AD used by method='slycot' (based on the generalized | | 205 | + # Schur method of Penzl (1998)) also handles singular E. | | 206 | + try: | | 207 | + At = solve(E, A) | | 208 | + Qt = solve(E, solve(E, Q).T).T | | 209 | + except np.linalg.LinAlgError: | | 210 | + raise ControlArgument( | | 211 | + "method='scipy' requires E to be nonsingular; " | | 212 | + "use method='slycot' (SLICOT sg03ad) for singular E") | | 213 | + _warn_ill_conditioned_E(E) | | 214 | + return sp.linalg.solve_continuous_lyapunov(At, -Qt) | 167 | 215 |
| 168 | 216 | # Make sure we have access to the write Slycot routine | 169 | 217 | try: | @@ -229,6 +277,32 @@ def dlyap(A, Q, C=None, E=None, method=None): | 229 | 277 | X : 2D array (or matrix) | 230 | 278 | Solution to the Lyapunov or Sylvester equation. | 231 | 279 |
| | 280 | + Notes | | 281 | + ----- | | 282 | + For the generalized Lyapunov equation, method='slycot' uses the | | 283 | + SLICOT routine SG03AD, based on the generalized Schur method of | | 284 | + Penzl [1]_, which also handles singular E. With method='scipy', the | | 285 | + equation is transformed to a standard Lyapunov equation by inverting | | 286 | + E, which requires E to be nonsingular and loses accuracy when E is | | 287 | + ill-conditioned (a UserWarning is then issued); method='slycot' does | | 288 | + not invert E and is preferable in that case. | | 289 | +
| | 290 | + For the Sylvester equation, method='slycot' uses the | | 291 | + Hessenberg-Schur method of the SLICOT routine SB04QD [2]_ and | | 292 | + method='scipy' uses the Bartels-Stewart method [3]_; both reduce the | | 293 | + coefficient matrices to (Hessenberg-)Schur form and solve the result | | 294 | + by back-substitution, with O(n^3 + m^3) cost. | | 295 | +
| | 296 | + References | | 297 | + ---------- | | 298 | + .. [1] Penzl, T., "Numerical solution of generalized Lyapunov | | 299 | + equations", Advances in Computational Mathematics, 8:33-48, 1998. | | 300 | + .. [2] Golub, G.H., Nash, S., and Van Loan, C., "A Hessenberg-Schur | | 301 | + method for the problem AX + XB = C", IEEE Trans. Automatic | | 302 | + Control, AC-24, pp. 909-913, 1979. | | 303 | + .. [3] Bartels, R.H. and Stewart, G.W., "Solution of the matrix | | 304 | + equation AX + XB = C", Comm. ACM, 15(9), pp. 820-826, 1972. | | 305 | +
| 232 | 306 | """ | 233 | 307 | # Decide what method to use | 234 | 308 | method = _slycot_or_scipy(method) | @@ -279,8 +353,40 @@ def dlyap(A, Q, C=None, E=None, method=None): | 279 | 353 | _check_shape(C, n, m, name="C") | 280 | 354 |
| 281 | 355 | if method == 'scipy': | 282 | | - raise ControlArgument( | 283 | | - "method='scipy' not valid for Sylvester equation") | | 356 | + # Solve the discrete-time Sylvester equation | | 357 | + # | | 358 | + # A X Q^T - X + C = 0 | | 359 | + # | | 360 | + # by the Bartels-Stewart method, matching the complexity of | | 361 | + # the Hessenberg-Schur algorithm of the SLICOT routine | | 362 | + # SB04QD used by method='slycot' (Golub, Nash, and Van | | 363 | + # Loan, 1979): with complex Schur forms A = U Ta U^H and | | 364 | + # Q^T = V Tq V^H and Y = U^H X V, the transformed equation | | 365 | + # Ta Y Tq - Y + U^H C V = 0 is solved column by column, | | 366 | + # each column requiring one triangular solve. O(n^3 + m^3) | | 367 | + # flops overall. | | 368 | + Ta, U = sp.linalg.schur(A, output='complex') | | 369 | + Tq, V = sp.linalg.schur(Q.T, output='complex') | | 370 | + Ct = U.conj().T @ C @ V | | 371 | + # Solvability requires lam_A * lam_Q != 1 for all pairs of | | 372 | + # eigenvalues (the diagonals of the triangular factors) | | 373 | + if np.min(np.abs(np.outer(np.diag(Tq), np.diag(Ta)) - 1.)) \ | | 374 | + < finfo(float).eps * max( | | 375 | + 1., np.abs(np.diag(Ta)).max() | | 376 | + * np.abs(np.diag(Tq)).max()): | | 377 | + raise ControlArgument( | | 378 | + "A and Q have a pair of eigenvalues whose product " | | 379 | + "is (almost) equal to 1; the discrete-time " | | 380 | + "Sylvester equation is singular") | | 381 | + Y = np.empty((n, m), dtype=complex) | | 382 | + TaY = np.empty((n, m), dtype=complex) # running Ta @ Y | | 383 | + In = np.eye(n) | | 384 | + for k in range(m): | | 385 | + rhs = -Ct[:, k] - TaY[:, :k] @ Tq[:k, k] | | 386 | + Y[:, k] = sp.linalg.solve_triangular( | | 387 | + Tq[k, k] * Ta - In, rhs) | | 388 | + TaY[:, k] = Ta @ Y[:, k] | | 389 | + return np.real(U @ Y @ V.conj().T) | 284 | 390 |
| 285 | 391 | # Solve the Sylvester equation by calling Slycot function sb04qd | 286 | 392 | X = sb04qd(n, m, -A, Q.T, C) | @@ -292,8 +398,23 @@ def dlyap(A, Q, C=None, E=None, method=None): | 292 | 398 | _check_shape(E, n, n, square=True, name="E") | 293 | 399 |
| 294 | 400 | if method == 'scipy': | 295 | | - raise ControlArgument( | 296 | | - "method='scipy' not valid for generalized Lyapunov equation") | | 401 | + # Transform to a standard Lyapunov equation by multiplying | | 402 | + # from the left by inv(E) and from the right by inv(E).T: | | 403 | + # | | 404 | + # (E^-1 A) X (E^-1 A)^T - X + E^-1 Q E^-T = 0 | | 405 | + # | | 406 | + # This requires E to be nonsingular; the SLICOT routine | | 407 | + # SG03AD used by method='slycot' (based on the generalized | | 408 | + # Schur method of Penzl (1998)) also handles singular E. | | 409 | + try: | | 410 | + At = solve(E, A) | | 411 | + Qt = solve(E, solve(E, Q).T).T | | 412 | + except np.linalg.LinAlgError: | | 413 | + raise ControlArgument( | | 414 | + "method='scipy' requires E to be nonsingular; " | | 415 | + "use method='slycot' (SLICOT sg03ad) for singular E") | | 416 | + _warn_ill_conditioned_E(E) | | 417 | + return sp.linalg.solve_discrete_lyapunov(At, Qt) | 297 | 418 |
| 298 | 419 | # Solve the generalized Lyapunov equation by calling Slycot | 299 | 420 | # function sg03ad | |
Collapse file Copy file name to clipboardExpand all lines: control/tests/mateqn_test.py +56-20Lines changed: 56 additions & 20 deletions | Original file line number | Diff line number | Diff line change |
|---|
@@ -90,19 +90,22 @@ def test_lyap_sylvester(self, method): | 90 | 90 | X_scipy = lyap(A, B, C, method='scipy') | 91 | 91 | assert_array_almost_equal(X_scipy, X) | 92 | 92 |
| 93 | | - @pytest.mark.slycot | 94 | | - def test_lyap_g(self): | | 93 | + @pytest.mark.parametrize('method', | | 94 | + ['scipy', | | 95 | + pytest.param('slycot', marks=pytest.mark.slycot)]) | | 96 | + def test_lyap_g(self, method): | 95 | 97 | A = array([[-1, 2], [-3, -4]]) | 96 | 98 | Q = array([[3, 1], [1, 1]]) | 97 | 99 | E = array([[1, 2], [2, 1]]) | 98 | | - X = lyap(A, Q, None, E) | | 100 | + X = lyap(A, Q, None, E, method=method) | 99 | 101 | # print("The solution obtained is ", X) | 100 | 102 | assert_array_almost_equal(A @ X @ E.T + E @ X @ A.T + Q, | 101 | 103 | zeros((2,2))) | 102 | 104 |
| 103 | | - # Make sure that trying to solve with SciPy generates an error | 104 | | - with pytest.raises(ControlArgument, match="'scipy' not valid"): | 105 | | - X = lyap(A, Q, None, E, method='scipy') | | 105 | + # Compare methods | | 106 | + if method == 'slycot': | | 107 | + X_scipy = lyap(A, Q, None, E, method='scipy') | | 108 | + assert_array_almost_equal(X_scipy, X) | 106 | 109 |
| 107 | 110 | @pytest.mark.parametrize('method', | 108 | 111 | ['scipy', | @@ -125,39 +128,73 @@ def test_dlyap(self, method): | 125 | 128 | X_scipy = dlyap(A,Q, method='scipy') | 126 | 129 | assert_array_almost_equal(X_scipy, X) | 127 | 130 |
| 128 | | - @pytest.mark.slycot | 129 | | - def test_dlyap_g(self): | | 131 | + @pytest.mark.parametrize('method', | | 132 | + ['scipy', | | 133 | + pytest.param('slycot', marks=pytest.mark.slycot)]) | | 134 | + def test_dlyap_g(self, method): | 130 | 135 | A = array([[-0.6, 0],[-0.1, -0.4]]) | 131 | 136 | Q = array([[3, 1],[1, 1]]) | 132 | 137 | E = array([[1, 1],[2, 1]]) | 133 | | - X = dlyap(A, Q, None, E) | | 138 | + X = dlyap(A, Q, None, E, method=method) | 134 | 139 | # print("The solution obtained is ", X) | 135 | 140 | assert_array_almost_equal(A @ X @ A.T - E @ X @ E.T + Q, | 136 | 141 | zeros((2,2))) | 137 | 142 |
| 138 | | - # Make sure that trying to solve with SciPy generates an error | 139 | | - with pytest.raises(ControlArgument, match="'scipy' not valid"): | 140 | | - X = dlyap(A, Q, None, E, method='scipy') | | 143 | + # Compare methods | | 144 | + if method == 'slycot': | | 145 | + X_scipy = dlyap(A, Q, None, E, method='scipy') | | 146 | + assert_array_almost_equal(X_scipy, X) | 141 | 147 |
| 142 | | - @pytest.mark.slycot | 143 | | - def test_dlyap_sylvester(self): | | 148 | + @pytest.mark.parametrize('method', | | 149 | + ['scipy', | | 150 | + pytest.param('slycot', marks=pytest.mark.slycot)]) | | 151 | + def test_dlyap_sylvester(self, method): | 144 | 152 | A = 5 | 145 | 153 | B = array([[4, 3], [4, 3]]) | 146 | 154 | C = array([2, 1]) | 147 | | - X = dlyap(A,B,C) | | 155 | + X = dlyap(A, B, C, method=method) | 148 | 156 | # print("The solution obtained is ", X) | 149 | 157 | assert_array_almost_equal(A * X @ B.T - X + C, zeros((1,2))) | 150 | 158 |
| 151 | 159 | A = array([[2, 1], [1, 2]]) | 152 | 160 | B = array([[1, 2], [0.5, 0.1]]) | 153 | 161 | C = array([[1, 0], [0, 1]]) | 154 | | - X = dlyap(A, B, C) | | 162 | + X = dlyap(A, B, C, method=method) | 155 | 163 | # print("The solution obtained is ", X) | 156 | 164 | assert_array_almost_equal(A @ X @ B.T - X + C, zeros((2,2))) | 157 | 165 |
| 158 | | - # Make sure that trying to solve with SciPy generates an error | 159 | | - with pytest.raises(ControlArgument, match="'scipy' not valid"): | 160 | | - X = dlyap(A, B, C, method='scipy') | | 166 | + # Compare methods | | 167 | + if method == 'slycot': | | 168 | + X_scipy = dlyap(A, B, C, method='scipy') | | 169 | + assert_array_almost_equal(X_scipy, X) | | 170 | + | | 171 | + @pytest.mark.parametrize("cdlyap", [lyap, dlyap]) | | 172 | + def test_lyap_g_singular_E_scipy(self, cdlyap): | | 173 | + """Generalized Lyapunov with singular E raises on scipy path""" | | 174 | + A = array([[-1, 2], [-3, -4]]) | | 175 | + Q = array([[3, 1], [1, 1]]) | | 176 | + E = array([[1, 2], [2, 4]]) # singular | | 177 | + with pytest.raises(ControlArgument, match="E to be nonsingular"): | | 178 | + cdlyap(A, Q, None, E, method='scipy') | | 179 | + | | 180 | + @pytest.mark.parametrize("cdlyap", [lyap, dlyap]) | | 181 | + def test_lyap_g_ill_conditioned_E_scipy(self, cdlyap): | | 182 | + """Generalized Lyapunov with ill-conditioned E warns on scipy path""" | | 183 | + A = array([[-1, 2], [-3, -4]]) | | 184 | + Q = array([[3, 1], [1, 1]]) | | 185 | + E = array([[1, 1], [1, 1 + 1e-12]]) # nonsingular, cond ~ 4e12 | | 186 | + with pytest.warns(UserWarning, match="ill-conditioned"): | | 187 | + cdlyap(A, Q, None, E, method='scipy') | | 188 | + | | 189 | + def test_dlyap_sylvester_singular_scipy(self): | | 190 | + """Discrete Sylvester with an eigenvalue product of 1 is singular""" | | 191 | + # eig(A) = {2, 3}, eig(Q) = {0.5, 0.7}; the product 2 * 0.5 = 1 | | 192 | + # makes the discrete-time Sylvester operator singular | | 193 | + A = array([[2., 0.], [0., 3.]]) | | 194 | + Q = array([[0.5, 0.], [0., 0.7]]) | | 195 | + C = array([[1., 0.], [0., 1.]]) | | 196 | + with pytest.raises(ControlArgument, match="singular"): | | 197 | + dlyap(A, Q, C, method='scipy') | 161 | 198 |
| 162 | 199 | @pytest.mark.parametrize('method', | 163 | 200 | ['scipy', | @@ -274,7 +311,6 @@ def test_dare(self, method): | 274 | 311 | lam = eigvals(A - B @ G) | 275 | 312 | assert_array_less(abs(lam), 1.0) | 276 | 313 |
| 277 | | - @pytest.mark.slycot | 278 | 314 | def test_dare_compare(self): | 279 | 315 | A = np.array([[-0.6, 0], [-0.1, -0.4]]) | 280 | 316 | Q = np.array([[2, 1], [1, 0]]) | |
|