Skip to content

Add unit test illustrating issue #935 + add method keyword for tf2ss #937

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Nov 8, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
47 changes: 37 additions & 10 deletions control/statesp.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
------
Expand All @@ -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.]]]
Expand Down Expand Up @@ -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
Expand All @@ -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,
Expand All @@ -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)
Expand All @@ -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 = \
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change for the two lines above:

if not issiso(sys):
    raise ControlMIMONotImplemented("MIMO system conversion not supported without Slycot")

sp.signal.tf2ss(squeeze(sys.num), squeeze(sys.den))
newsys = StateSpace(A, B, C, D, sys.dt)
else:
raise ValueError(f"unknown {method=}")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

First instantiation of self-documenting F-strings in python-control? : )


# Copy over the signal (and system) names
newsys._copy_names(
Expand Down
4 changes: 2 additions & 2 deletions control/tests/convert_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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)
Expand Down
46 changes: 46 additions & 0 deletions control/tests/statesp_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)