diff --git a/control/mateqn.py b/control/mateqn.py index b23e42def..312efad13 100644 --- a/control/mateqn.py +++ b/control/mateqn.py @@ -14,7 +14,7 @@ import numpy as np import scipy as sp -from numpy import eye, finfo, inexact +from numpy import eye, finfo from scipy.linalg import eigvals, solve from .exception import ControlArgument, ControlDimension, ControlSlycot, \ @@ -81,7 +81,7 @@ def _warn_ill_conditioned_E(E): # -def lyap(A, Q, C=None, E=None, method=None): +def lyap(A, Q, C=None, E=None, method=None, **kwargs): """Solves the continuous-time Lyapunov equation. X = lyap(A, Q) solves @@ -117,6 +117,9 @@ def lyap(A, Q, C=None, E=None, method=None): Set the method used for computing the result. Current methods are 'slycot' and 'scipy'. If set to None (default), try 'slycot' first and then 'scipy'. + symmetric_kwargs : dict, optional + Keyword arguments passed to the SciPy symmetry/Hermitian check, + such as `atol` and `rtol`. Returns ------- @@ -143,6 +146,12 @@ def lyap(A, Q, C=None, E=None, method=None): equations", Advances in Computational Mathematics, 8:33-48, 1998. """ + + symmetric_kwargs = kwargs.pop("symmetric_kwargs", {}) + + # Make sure there were no extraneous keywords + if kwargs: + raise TypeError("unrecognized keyword(s): ", str(kwargs)) # Decide what method to use method = _slycot_or_scipy(method) if method == 'slycot': @@ -169,11 +178,11 @@ def lyap(A, Q, C=None, E=None, method=None): # Solve standard Lyapunov equation if C is None and E is None: # Check to make sure input matrices are the right shape and type - _check_shape(Q, n, n, square=True, symmetric=True, name="Q") + _check_shape(Q, n, n, square=True, symmetric=True, name="Q", symmetric_kwargs=symmetric_kwargs) if method == 'scipy': - # Solve the Lyapunov equation using SciPy - return sp.linalg.solve_continuous_lyapunov(A, -Q) + # Solve the Lyapunov equation using SciPy + return sp.linalg.solve_continuous_lyapunov(A, -Q) # Solve the Lyapunov equation by calling Slycot function sb03md with warnings.catch_warnings(): @@ -243,7 +252,7 @@ def lyap(A, Q, C=None, E=None, method=None): return X -def dlyap(A, Q, C=None, E=None, method=None): +def dlyap(A, Q, C=None, E=None, method=None, **kwargs): """Solves the discrete-time Lyapunov equation. X = dlyap(A, Q) solves @@ -279,6 +288,9 @@ def dlyap(A, Q, C=None, E=None, method=None): Set the method used for computing the result. Current methods are 'slycot' and 'scipy'. If set to None (default), try 'slycot' first and then 'scipy'. + symmetric_kwargs : dict, optional + Keyword arguments passed to the SciPy symmetry/Hermitian check, + such as `atol` and `rtol`. Returns ------- @@ -316,6 +328,11 @@ def dlyap(A, Q, C=None, E=None, method=None): equation AX + XB = C", Comm. ACM, 15(9), pp. 820-826, 1972. """ + + symmetric_kwargs = kwargs.pop("symmetric_kwargs", {}) + # Make sure there were no extraneous keywords + if kwargs: + raise TypeError("unrecognized keyword(s): ", str(kwargs)) # Decide what method to use method = _slycot_or_scipy(method) @@ -346,7 +363,7 @@ def dlyap(A, Q, C=None, E=None, method=None): # Solve standard Lyapunov equation if C is None and E is None: # Check to make sure input matrices are the right shape and type - _check_shape(Q, n, n, square=True, symmetric=True, name="Q") + _check_shape(Q, n, n, square=True, symmetric=True, name="Q", symmetric_kwargs=symmetric_kwargs) if method == 'scipy': # Solve the Lyapunov equation using SciPy @@ -450,7 +467,7 @@ def dlyap(A, Q, C=None, E=None, method=None): # def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None, - _As="A", _Bs="B", _Qs="Q", _Rs="R", _Ss="S", _Es="E"): + _As="A", _Bs="B", _Qs="Q", _Rs="R", _Ss="S", _Es="E", **kwargs): """Solves the continuous-time algebraic Riccati equation. X, L, G = care(A, B, Q, R=None) solves @@ -484,6 +501,9 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None, Set the method used for computing the result. Current methods are 'slycot' and 'scipy'. If set to None (default), try 'slycot' first and then 'scipy'. + symmetric_kwargs : dict, optional + Keyword arguments passed to the SciPy symmetry/Hermitian check, + such as `atol` and `rtol`. stabilizing : bool, optional If `method` is 'slycot', unstabilized eigenvalues will be returned in the initial elements of `L`. Not supported for 'scipy'. @@ -498,6 +518,13 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None, Gain matrix. """ + + symmetric_kwargs = kwargs.pop("symmetric_kwargs", {}) + + # Make sure there were no extraneous keywords + if kwargs: + raise TypeError("unrecognized keyword(s): ", str(kwargs)) + # Decide what method to use method = _slycot_or_scipy(method) @@ -518,8 +545,8 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None, # Check to make sure input matrices are the right shape and type _check_shape(A, n, n, square=True, name=_As) _check_shape(B, n, m, name=_Bs) - _check_shape(Q, n, n, square=True, symmetric=True, name=_Qs) - _check_shape(R, m, m, square=True, symmetric=True, name=_Rs) + _check_shape(Q, n, n, square=True, symmetric=True, name=_Qs, symmetric_kwargs=symmetric_kwargs) + _check_shape(R, m, m, square=True, symmetric=True, name=_Rs, symmetric_kwargs=symmetric_kwargs) # Solve the standard algebraic Riccati equation if S is None and E is None: @@ -606,7 +633,7 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None, return X, L, G def dare(A, B, Q, R, S=None, E=None, stabilizing=True, method=None, - _As="A", _Bs="B", _Qs="Q", _Rs="R", _Ss="S", _Es="E"): + _As="A", _Bs="B", _Qs="Q", _Rs="R", _Ss="S", _Es="E", **kwargs): """Solves the discrete-time algebraic Riccati equation. X, L, G = dare(A, B, Q, R) solves @@ -640,6 +667,9 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True, method=None, Set the method used for computing the result. Current methods are 'slycot' and 'scipy'. If set to None (default), try 'slycot' first and then 'scipy'. + symmetric_kwargs : dict, optional + Keyword arguments passed to the SciPy symmetry/Hermitian check, + such as `atol` and `rtol`. stabilizing : bool, optional If `method` is 'slycot', unstabilized eigenvalues will be returned in the initial elements of `L`. Not supported for 'scipy'. @@ -654,6 +684,13 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True, method=None, Gain matrix. """ + + symmetric_kwargs = kwargs.pop("symmetric_kwargs", {}) + + # Make sure there were no extraneous keywords + if kwargs: + raise TypeError("unrecognized keyword(s): ", str(kwargs)) + # Decide what method to use method = _slycot_or_scipy(method) @@ -674,8 +711,8 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True, method=None, # Check to make sure input matrices are the right shape and type _check_shape(A, n, n, square=True, name=_As) _check_shape(B, n, m, name=_Bs) - _check_shape(Q, n, n, square=True, symmetric=True, name=_Qs) - _check_shape(R, m, m, square=True, symmetric=True, name=_Rs) + _check_shape(Q, n, n, square=True, symmetric=True, name=_Qs, symmetric_kwargs=symmetric_kwargs) + _check_shape(R, m, m, square=True, symmetric=True, name=_Rs, symmetric_kwargs=symmetric_kwargs) if E is not None: _check_shape(E, n, n, square=True, name=_Es) if S is not None: @@ -740,7 +777,8 @@ def _slycot_or_scipy(method): # Utility function to check matrix dimensions -def _check_shape(M, n, m, square=False, symmetric=False, name="??"): +def _check_shape(M, n, m, square=False, symmetric=False, name="??", + symmetric_kwargs=None): """Check the shape and properties of a 2D array. This function can be used to check to make sure a 2D array_like has the @@ -773,7 +811,7 @@ def _check_shape(M, n, m, square=False, symmetric=False, name="??"): if (square or symmetric) and M.shape[0] != M.shape[1]: raise ControlDimension("%s must be a square matrix" % name) - if symmetric and not _is_symmetric(M): + if symmetric and not _is_symmetric(M, symmetric_kwargs=symmetric_kwargs): raise ControlArgument("%s must be a symmetric matrix" % name) if M.shape[0] != n or M.shape[1] != m: @@ -785,10 +823,13 @@ def _check_shape(M, n, m, square=False, symmetric=False, name="??"): # Utility function to check if a matrix is symmetric -def _is_symmetric(M): +def _is_symmetric(M, symmetric_kwargs=None): M = np.atleast_2d(M) - if isinstance(M[0, 0], inexact): - eps = finfo(M.dtype).eps - return ((M - M.T) < eps).all() + symmetric_kwargs = ( + symmetric_kwargs.copy() if symmetric_kwargs else {} + ) + + if np.iscomplexobj(M): + return sp.linalg.ishermitian(M, **symmetric_kwargs) else: - return (M == M.T).all() + return sp.linalg.issymmetric(M, **symmetric_kwargs) \ No newline at end of file diff --git a/control/stochsys.py b/control/stochsys.py index 22e2f738b..83909acfa 100644 --- a/control/stochsys.py +++ b/control/stochsys.py @@ -83,6 +83,9 @@ def lqe(*args, **kwargs): Set the method used for computing the result. Current methods are 'slycot' and 'scipy'. If set to None (default), try 'slycot' first and then 'scipy'. + symmetric_kwargs : dict, optional + Keyword arguments passed to the SciPy symmetry/Hermitian check, + such as `atol` and `rtol`. Returns ------- @@ -133,6 +136,8 @@ def lqe(*args, **kwargs): # Get the method to use (if specified as a keyword) method = kwargs.pop('method', None) + symmetric_kwargs = kwargs.pop('symmetric_kwargs', {}) + if kwargs: raise TypeError("unrecognized keyword(s): ", str(kwargs)) @@ -178,7 +183,7 @@ def lqe(*args, **kwargs): # Compute the result (dimension and symmetry checking done in care()) P, E, LT = care(A.T, C.T, G @ QN @ G.T, RN, method=method, - _Bs="C", _Qs="QN", _Rs="RN", _Ss="NN") + _Bs="C", _Qs="QN", _Rs="RN", _Ss="NN", symmetric_kwargs=symmetric_kwargs) return LT.T, P, E @@ -220,6 +225,9 @@ def dlqe(*args, **kwargs): Set the method used for computing the result. Current methods are 'slycot' and 'scipy'. If set to None (default), try 'slycot' first and then 'scipy'. + symmetric_kwargs : dict, optional + Keyword arguments passed to the SciPy symmetry/Hermitian check, + such as `atol` and `rtol`. Returns ------- @@ -252,6 +260,8 @@ def dlqe(*args, **kwargs): # Get the method to use (if specified as a keyword) method = kwargs.pop('method', None) + symmetric_kwargs = kwargs.pop('symmetric_kwargs', {}) + if kwargs: raise TypeError("unrecognized keyword(s): ", str(kwargs)) @@ -299,7 +309,7 @@ def dlqe(*args, **kwargs): # Compute the result (dimension and symmetry checking done in dare()) P, E, LT = dare(A.T, C.T, G @ QN @ G.T, RN, method=method, - _Bs="C", _Qs="QN", _Rs="RN", _Ss="NN") + _Bs="C", _Qs="QN", _Rs="RN", _Ss="NN", symmetric_kwargs=symmetric_kwargs) return LT.T, P, E diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py index 566b35a28..bc8e1ef65 100644 --- a/control/tests/kwargs_test.py +++ b/control/tests/kwargs_test.py @@ -99,9 +99,12 @@ def test_kwarg_search(module, prefix): @pytest.mark.parametrize( "function, nsssys, ntfsys, moreargs, kwargs", [(control.append, 2, 0, (), {}), + (control.care, 0, 0, ([[-1]], [[1]], [[1]], [[1]]), {'method': 'scipy'}), (control.combine_tf, 0, 0, ([[1, 0], [0, 1]], ), {}), + (control.dare, 0, 0, ([[0.5]], [[1]], [[1]], [[1]]), {'method': 'scipy'}), (control.dlqe, 1, 0, ([[1]], [[1]]), {}), (control.dlqr, 1, 0, ([[1, 0], [0, 1]], [[1]]), {}), + (control.dlyap, 0, 0, ([[0.5]], [[1]]), {'method': 'scipy'}), (control.drss, 0, 0, (2, 1, 1), {}), (control.feedback, 2, 0, (), {}), (control.flatsys.flatsys, 1, 0, (), {}), @@ -109,6 +112,7 @@ def test_kwarg_search(module, prefix): (control.lqe, 1, 0, ([[1]], [[1]]), {}), (control.lqr, 1, 0, ([[1, 0], [0, 1]], [[1]]), {}), (control.linearize, 1, 0, (0, 0), {}), + (control.lyap, 0, 0, ([[-1]], [[1]]), {'method': 'scipy'}), (control.negate, 1, 0, (), {}), (control.nlsys, 0, 0, (lambda t, x, u, params: np.array([0]),), {}), (control.parallel, 2, 0, (), {}), @@ -250,14 +254,17 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo): 'bode': test_response_plot_kwargs, 'bode_plot': test_response_plot_kwargs, 'LTI.bode_plot': test_response_plot_kwargs, # tested via bode_plot + 'care': test_unrecognized_kwargs, 'combine_tf': test_unrecognized_kwargs, 'create_estimator_iosystem': stochsys_test.test_estimator_errors, 'create_statefbk_iosystem': statefbk_test.TestStatefbk.test_statefbk_errors, + 'dare': test_unrecognized_kwargs, 'describing_function_plot': test_matplotlib_kwargs, 'describing_function_response': descfcn_test.test_describing_function_exceptions, 'dlqe': test_unrecognized_kwargs, 'dlqr': test_unrecognized_kwargs, + 'dlyap': test_unrecognized_kwargs, 'drss': test_unrecognized_kwargs, 'feedback': test_unrecognized_kwargs, 'find_eqpt': iosys_test.test_find_operating_point, @@ -275,6 +282,7 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo): 'linearize': test_unrecognized_kwargs, 'lqe': test_unrecognized_kwargs, 'lqr': test_unrecognized_kwargs, + 'lyap': test_unrecognized_kwargs, 'LTI.forced_response': statesp_test.test_convenience_aliases, 'LTI.impulse_response': statesp_test.test_convenience_aliases, 'LTI.initial_response': statesp_test.test_convenience_aliases, diff --git a/control/tests/mateqn_test.py b/control/tests/mateqn_test.py index 5cc34c2ab..7c604fdc0 100644 --- a/control/tests/mateqn_test.py +++ b/control/tests/mateqn_test.py @@ -39,7 +39,7 @@ import pytest from scipy.linalg import eigvals, solve -from control.mateqn import lyap, dlyap, care, dare +from control.mateqn import lyap, dlyap, care, dare, _is_symmetric from control.exception import ControlArgument, ControlDimension @@ -466,3 +466,25 @@ def test_raise(self): cdare(A, B, Qfs, R, S, E) with pytest.raises(ControlArgument): cdare(A, B, Q, Rfs, S, E) + + def test_is_symmetric_scale_aware(self): + M = np.array([ + [1e8, 1e8], + [1e8 + 1e-8, 1e8] + ]) + assert not _is_symmetric(M) + assert _is_symmetric(M,symmetric_kwargs={"rtol": 1e-12},) + + def test_is_symmetric_rejects_asymmetric(self): + M = np.array([ + [1., 2.], + [5., 1.] + ]) + assert not _is_symmetric(M) + + def test_is_symmetric_complex_hermitian(self): + M = np.array([ + [1., 2. + 1.j], + [2. - 1.j, 3.] + ]) + assert _is_symmetric(M) \ No newline at end of file diff --git a/control/tests/stochsys_test.py b/control/tests/stochsys_test.py index 20e799643..77317d918 100644 --- a/control/tests/stochsys_test.py +++ b/control/tests/stochsys_test.py @@ -34,7 +34,52 @@ def test_LQE(method): A, G, C, QN, RN = (np.array([[X]]) for X in [0., .1, 1., 10., 2.]) L, P, poles = lqe(A, G, C, QN, RN, method=method) check_LQE(L, P, poles, G, QN, RN) - +def test_lqe_symmetric_kwargs(): + A = np.array([ + [-1., 0.], + [0., -2.] + ]) + G = np.eye(2) + C = np.eye(2) + QN = np.array([ + [1., 0.2], + [0.2 + 1e-15, 1.] + ]) + RN = np.eye(2) + + # Exact symmetry check should fail + with pytest.raises(ControlArgument, match="symmetric"): + lqe(A, G, C, QN, RN, method="scipy") + + # Passing rtol should allow the nearly symmetric matrix + lqe( + A, G, C, QN, RN, + method="scipy", + symmetric_kwargs={"rtol": 1e-12}, + ) +def test_dlqe_symmetric_kwargs(): + A = np.array([ + [0.5, 0.], + [0., 0.4] + ]) + G = np.eye(2) + C = np.eye(2) + QN = np.array([ + [1., 0.2], + [0.2 + 1e-15, 1.] + ]) + RN = np.eye(2) + + # Exact symmetry check should fail + with pytest.raises(ControlArgument, match="symmetric"): + dlqe(A, G, C, QN, RN, method="scipy") + + # Passing rtol should allow the nearly symmetric matrix + dlqe( + A, G, C, QN, RN, + method="scipy", + symmetric_kwargs={"rtol": 1e-12}, + ) @pytest.mark.parametrize("cdlqe", [lqe, dlqe]) def test_lqe_call_format(cdlqe): # Create a random state space system for testing