From fe8888f94e2f126599cc28dfdec22926f0b28984 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 21 Oct 2023 22:48:58 -0700 Subject: [PATCH 1/3] add unit test illustrating issue #935 + add keyword for tf2ss --- control/statesp.py | 22 +++++++++++++++------- control/tests/statesp_test.py | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 7 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index 38dd2388d..8ec7f0fc4 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -60,7 +60,7 @@ from scipy.signal import StateSpace as signalStateSpace from warnings import warn -from .exception import ControlSlycot +from .exception import ControlSlycot, slycot_check from .frdata import FrequencyResponseData from .lti import LTI, _process_frequency_response from .iosys import InputOutputSystem, common_timebase, isdtime, \ @@ -1615,10 +1615,13 @@ def ss(*args, **kwargs): warn("state labels specified for " "non-unique state space realization") + # Allow method to be specified (eg, tf2ss) + method = kwargs.pop('method', None) + # Create a state space system from an LTI system sys = StateSpace( _convert_to_statespace( - sys, + sys, method=method, use_prefix_suffix=not sys._generic_name_check()), **kwargs) @@ -2189,7 +2192,7 @@ def _f2s(f): return s -def _convert_to_statespace(sys, use_prefix_suffix=False): +def _convert_to_statespace(sys, use_prefix_suffix=False, method=None): """Convert a system to state space form (if needed). If sys is already a state space, then it is returned. If sys is a @@ -2213,13 +2216,17 @@ def _convert_to_statespace(sys, use_prefix_suffix=False): raise ValueError("transfer function is non-proper; can't " "convert to StateSpace system") - try: + if method is None and slycot_check() or method == 'slycot': + if not slycot_check(): + raise ValueError("method='slycot' requires slycot") + from slycot import td04ad # Change the numerator and denominator arrays so that the transfer # function matrix has a common denominator. # matrices are also sized/padded to fit td04ad num, den, denorder = sys.minreal()._common_den() + num, den, denorder = sys._common_den() # transfer function to state space conversion now should work! ssout = td04ad('C', sys.ninputs, sys.noutputs, @@ -2230,9 +2237,8 @@ def _convert_to_statespace(sys, use_prefix_suffix=False): ssout[1][:states, :states], ssout[2][:states, :sys.ninputs], ssout[3][:sys.noutputs, :states], ssout[4], sys.dt) - except ImportError: - # No Slycot. Scipy tf->ss can't handle MIMO, but static - # MIMO is an easy special case we can check for here + elif method in [None, 'scipy']: + # Scipy tf->ss can't handle MIMO, but SISO is OK maxn = max(max(len(n) for n in nrow) for nrow in sys.num) maxd = max(max(len(d) for d in drow) @@ -2250,6 +2256,8 @@ def _convert_to_statespace(sys, use_prefix_suffix=False): A, B, C, D = \ sp.signal.tf2ss(squeeze(sys.num), squeeze(sys.den)) newsys = StateSpace(A, B, C, D, sys.dt) + else: + raise ValueError(f"unknown {method=}") # Copy over the signal (and system) names newsys._copy_names( diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py index 20fd62ca2..1c1a050aa 100644 --- a/control/tests/statesp_test.py +++ b/control/tests/statesp_test.py @@ -1203,3 +1203,37 @@ def test_params_warning(): sys.output(0, [0], [0], {'k': 5}) +# Check that tf2ss returns stable system (see issue #935) +@pytest.mark.parametrize("method", [ + # pytest.param(None), # use this one when SLICOT bug is sorted out + pytest.param( # remove this one when SLICOT bug is sorted out + None, marks=pytest.mark.xfail( + ct.slycot_check(), reason="tf2ss SLICOT bug")), + pytest.param( + 'slycot', marks=[ + pytest.mark.xfail( + not ct.slycot_check(), reason="slycot not installed"), + pytest.mark.xfail( # remove this one when SLICOT bug is sorted out + ct.slycot_check(), reason="tf2ss SLICOT bug")]), + pytest.param('scipy') +]) +def test_tf2ss_unstable(method): + num = np.array([ + 9.94004350e-13, 2.67602795e-11, 2.31058712e-10, 1.15119493e-09, + 5.04635153e-09, 1.34066064e-08, 2.11938725e-08, 2.39940325e-08, + 2.05897777e-08, 1.17092854e-08, 4.71236875e-09, 1.19497537e-09, + 1.90815347e-10, 1.00655454e-11, 1.47388887e-13, 8.40314881e-16, + 1.67195685e-18]) + den = np.array([ + 9.43513863e-11, 6.05312352e-08, 7.92752628e-07, 5.23764693e-06, + 1.82502556e-05, 1.24355899e-05, 8.68206174e-06, 2.73818482e-06, + 4.29133144e-07, 3.85554417e-08, 1.62631575e-09, 8.41098151e-12, + 9.85278302e-15, 4.07646645e-18, 5.55496497e-22, 3.06560494e-26, + 5.98908988e-31]) + + tf_sys = ct.tf(num, den) + ss_sys = ct.tf2ss(tf_sys, method=method) + + tf_poles = np.sort(tf_sys.poles()) + ss_poles = np.sort(ss_sys.poles()) + np.testing.assert_allclose(tf_poles, ss_poles, rtol=1e-4) From 2580d1862e5bbcb6b7d4e4ca03299332b0cf69db Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sun, 22 Oct 2023 13:13:07 -0700 Subject: [PATCH 2/3] add some documentation and notes --- control/statesp.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/control/statesp.py b/control/statesp.py index 8ec7f0fc4..eb52e848e 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -1583,6 +1583,13 @@ def ss(*args, **kwargs): -------- tf, ss2tf, tf2ss + Notes + ----- + If a transfer function is passed as the sole positional argument, the + system will be converted to state space form in the same way as calling + :func:`~control.tf2ss`. The `method` keyword can be used to select the + method for conversion. + Examples -------- Create a Linear I/O system object from matrices. @@ -1768,6 +1775,10 @@ def tf2ss(*args, **kwargs): name : string, optional System name. If unspecified, a generic name is generated with a unique integer id. + method : str, optional + 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' (SISO only). Raises ------ @@ -1784,6 +1795,13 @@ def tf2ss(*args, **kwargs): tf ss2tf + Notes + ----- + The ``slycot`` routine used to convert a transfer function into state + space form appears to have a bug and in some (rare) instances may not + return a system with the same poles as the input transfer function. + For SISO systems, setting ``method=scipy`` can be used as an alternative. + Examples -------- >>> num = [[[1., 2.], [3., 4.]], [[5., 6.], [7., 8.]]] From 31793e6d4a75faa9f81d730df5f100b4eb79ccd1 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Tue, 7 Nov 2023 22:17:59 -0800 Subject: [PATCH 3/3] updated MIMO not implemented w/ unit test, per @sawyerbfuller --- control/statesp.py | 9 +++++---- control/tests/convert_test.py | 4 ++-- control/tests/statesp_test.py | 12 ++++++++++++ 3 files changed, 19 insertions(+), 6 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index eb52e848e..e14a8358a 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -60,10 +60,10 @@ from scipy.signal import StateSpace as signalStateSpace from warnings import warn -from .exception import ControlSlycot, slycot_check +from .exception import ControlSlycot, slycot_check, ControlMIMONotImplemented from .frdata import FrequencyResponseData from .lti import LTI, _process_frequency_response -from .iosys import InputOutputSystem, common_timebase, isdtime, \ +from .iosys import InputOutputSystem, common_timebase, isdtime, issiso, \ _process_iosys_keywords, _process_dt_keyword, _process_signal_list from .nlsys import NonlinearIOSystem, InterconnectedSystem from . import config @@ -2268,8 +2268,9 @@ def _convert_to_statespace(sys, use_prefix_suffix=False, method=None): D[i, j] = sys.num[i][j][0] / sys.den[i][j][0] newsys = StateSpace([], [], [], D, sys.dt) else: - if sys.ninputs != 1 or sys.noutputs != 1: - raise TypeError("No support for MIMO without slycot") + if not issiso(sys): + raise ControlMIMONotImplemented( + "MIMO system conversion not supported without Slycot") A, B, C, D = \ sp.signal.tf2ss(squeeze(sys.num), squeeze(sys.den)) diff --git a/control/tests/convert_test.py b/control/tests/convert_test.py index 14f3133e1..7975bbe5a 100644 --- a/control/tests/convert_test.py +++ b/control/tests/convert_test.py @@ -21,7 +21,7 @@ from control import rss, ss, ss2tf, tf, tf2ss from control.statefbk import ctrb, obsv from control.freqplot import bode -from control.exception import slycot_check +from control.exception import slycot_check, ControlMIMONotImplemented from control.tests.conftest import slycotonly @@ -167,7 +167,7 @@ def testConvertMIMO(self): # Convert to state space and look for an error if (not slycot_check()): - with pytest.raises(TypeError): + with pytest.raises(ControlMIMONotImplemented): tf2ss(tsys) else: ssys = tf2ss(tsys) diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py index 1c1a050aa..59f441456 100644 --- a/control/tests/statesp_test.py +++ b/control/tests/statesp_test.py @@ -1237,3 +1237,15 @@ def test_tf2ss_unstable(method): tf_poles = np.sort(tf_sys.poles()) ss_poles = np.sort(ss_sys.poles()) np.testing.assert_allclose(tf_poles, ss_poles, rtol=1e-4) + + +def test_tf2ss_mimo(): + sys_tf = ct.tf([[[1], [1, 1, 1]]], [[[1, 1, 1], [1, 2, 1]]]) + + if ct.slycot_check(): + sys_ss = ct.ss(sys_tf) + np.testing.assert_allclose( + np.sort(sys_tf.poles()), np.sort(sys_ss.poles())) + else: + with pytest.raises(ct.ControlMIMONotImplemented): + sys_ss = ct.ss(sys_tf)