add unit test illustrating issue #935 + add keyword for tf2ss

murrayrm · murrayrm · commit fe8888f94e2f · 2023-10-22T07:10:32.000-07:00
diff --git 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
@@ -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)
