Merge pull request #937 from murrayrm/tf2ss_unstable_bug-21Oct2023

murrayrm · web-flow · commit 2ddc9ea647ca · 2023-11-07T22:39:13.000-08:00
diff --git 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
+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
@@ -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.
@@ -1615,10 +1622,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)
 
@@ -1765,6 +1775,10 @@ def tf2ss(*args, **kwargs):
     name : string, optional
         System name. If unspecified, a generic name <sys[id]> 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
     ------
@@ -1781,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.]]]
@@ -2189,7 +2210,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 +2234,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 +2255,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)
@@ -2244,12 +2268,15 @@ def _convert_to_statespace(sys, use_prefix_suffix=False):
                     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))
                 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/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
@@ -1203,3 +1203,49 @@ 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)
+
+
+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)
