Repository: python-control/python-control

Issue:
sys.dcgain() should return a real value rather than a complex value (until/unless we support complex coefficients)

Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/38101161/badge)](https://coveralls.io/builds/38101161)

Coverage increased (+0.004%) to 88.949% when pulling **f9db86a4cacff9f9bddf0c0d28c57c3042076953 on bnavigator:dcgain-real** into **22b99531589a2bbc2cdb9459c08eaf511eff03d8 on python-control:master**.


Patch:
diff --git a/control/iosys.py b/control/iosys.py
index e28de59f..7ed4c8b0 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -358,88 +358,26 @@ class InputOutputSystem(object):
         if (warning):
             warn("Parameters passed to InputOutputSystem ignored.")
 
-    def _rhs(self, t, x, u, params={}):
+    def _rhs(self, t, x, u):
         """Evaluate right hand side of a differential or difference equation.
 
         Private function used to compute the right hand side of an
-        input/output system model. Intended for fast
-        evaluation; for a more user-friendly interface
-        you may want to use :meth:`dynamics`.
+        input/output system model.
 
         """
         NotImplemented("Evaluation not implemented for system of type ",
                        type(self))
 
-    def dynamics(self, t, x, u):
-        """Compute the dynamics of a differential or difference equation.
-
-        Given time `t`, input `u` and state `x`, returns the value of the
-        right hand side of the dynamical system. If the system is continuous,
-        returns the time derivative
-
-            dx/dt = f(t, x, u)
-
-        where `f` is the system's (possibly nonlinear) dynamics function.
-        If the system is discrete-time, returns the next value of `x`:
-
-            x[t+dt] = f(t, x[t], u[t])
-
-        Where `t` is a scalar.
-
-        The inputs `x` and `u` must be of the correct length.
-
-        Parameters
-        ----------
-        t : float
-            the time at which to evaluate
-        x : array_like
-            current state
-        u : array_like
-            input
-
-        Returns
-        -------
-        dx/dt or x[t+dt] : ndarray
-        """
-        return self._rhs(t, x, u)
-
     def _out(self, t, x, u, params={}):
         """Evaluate the output of a system at a given state, input, and time
 
         Private function used to compute the output of of an input/output
-        system model given the state, input, parameters. Intended for fast
-        evaluation; for a more user-friendly interface you may want to use
-        :meth:`output`.
+        system model given the state, input, parameters, and time.
 
         """
         # If no output function was defined in subclass, return state
         return x
 
-    def output(self, t, x, u):
-        """Compute the output of the system
-
-        Given time `t`, input `u` and state `x`, returns the output of the
-        system:
-
-            y = g(t, x, u)
-
-        The inputs `x` and `u` must be of the correct length.
-
-        Parameters
-        ----------
-        t : float
-            the time at which to evaluate
-        x : array_like
-            current state
-        u : array_like
-            input
-
-        Returns
-        -------
-        y : ndarray
-        """
-        return self._out(t, x, u)
-
     def set_inputs(self, inputs, prefix='u'):
         """Set the number/names of the system inputs.
 
diff --git a/control/lti.py b/control/lti.py
index 30569863..01d04e02 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -208,6 +208,13 @@ class LTI:
         raise NotImplementedError("dcgain not implemented for %s objects" %
                                   str(self.__class__))
 
+    def _dcgain(self, warn_infinite):
+        zeroresp = self(0 if self.isctime() else 1,
+                        warn_infinite=warn_infinite)
+        if np.all(np.logical_or(np.isreal(zeroresp), np.isnan(zeroresp.imag))):
+            return zeroresp.real
+        else:
+            return zeroresp
 
 # Test to see if a system is SISO
 def issiso(sys, strict=False):
diff --git a/control/margins.py b/control/margins.py
index 0b53f26e..d22b6a0f 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -218,7 +218,7 @@ def _likely_numerical_inaccuracy(sys):
         # * z**(-p_q)
         x = [1] + [0] * (-p_q)
         p1 = np.polymul(p1, x)
-    return np.linalg.norm(p1) < 1e-4 * np.linalg.norm(p2)
+    return np.linalg.norm(p1) < 1e-3 * np.linalg.norm(p2)
 
 # Took the framework for the old function by
 # Sawyer B. Fuller <minster@uw.edu>, removed a lot of the innards
@@ -337,8 +337,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0, method='best'):
         if isinstance(sys, xferfcn.TransferFunction) and not sys.isctime():
             if _likely_numerical_inaccuracy(sys):
                 warn("stability_margins: Falling back to 'frd' method "
-                "because of chance of numerical inaccuracy in 'poly' method.",
-                stacklevel=2)
+                "because of chance of numerical inaccuracy in 'poly' method.")
                 omega_sys = freqplot._default_frequency_range(sys)
                 omega_sys = omega_sys[omega_sys < np.pi / sys.dt]
                 sys = frdata.FRD(sys, omega_sys, smooth=True)
diff --git a/control/mateqn.py b/control/mateqn.py
index 28b01d28..0b129fd9 100644
--- a/control/mateqn.py
+++ b/control/mateqn.py
@@ -5,6 +5,9 @@
 #
 # Author: Bjorn Olofsson
 
+# Python 3 compatibility (needs to go here)
+from __future__ import print_function
+
 # Copyright (c) 2011, All rights reserved.
 
 # Redistribution and use in source and binary forms, with or without
@@ -41,34 +44,6 @@ from scipy.linalg import eigvals, solve_discrete_are, solve
 from .exception import ControlSlycot, ControlArgument
 from .statesp import _ssmatrix
 
-# Make sure we have access to the right slycot routines
-try:
-    from slycot import sb03md57
-    # wrap without the deprecation warning
-    def sb03md(n, C, A, U, dico, job='X',fact='N',trana='N',ldwork=None):
-        ret = sb03md57(A, U, C, dico, job, fact, trana, ldwork)
-        return ret[2:]
-except ImportError:
-    try:
-        from slycot import sb03md
-    except ImportError:
-        sb03md = None
-
-try:
-    from slycot import sb04md
-except ImportError:
-    sb04md = None
-
-try:
-    from slycot import sb04qd
-except ImportError:
-    sb0qmd = None
-
-try:
-    from slycot import sg03ad
-except ImportError:
-    sb04ad = None
-
 __all__ = ['lyap', 'dlyap', 'dare', 'care']
 
 #
@@ -118,11 +93,16 @@ def lyap(A, Q, C=None, E=None):
     state space operations.  See :func:`~control.use_numpy_matrix`.
     """
 
-    if sb03md is None:
+    # Make sure we have access to the right slycot routines
+    try:
+        from slycot import sb03md
+    except ImportError:
         raise ControlSlycot("can't find slycot module 'sb03md'")
-    if sb04md is None:
-        raise ControlSlycot("can't find slycot module 'sb04md'")
 
+    try:
+        from slycot import sb04md
+    except ImportError:
+        raise ControlSlycot("can't find slycot module 'sb04md'")
 
     # Reshape 1-d arrays
     if len(shape(A)) == 1:
@@ -299,11 +279,19 @@ def dlyap(A, Q, C=None, E=None):
     of the same dimension. """
 
     # Make sure we have access to the right slycot routines
-    if sb03md is None:
+    try:
+        from slycot import sb03md
+    except ImportError:
         raise ControlSlycot("can't find slycot module 'sb03md'")
-    if sb04qd is None:
+
+    try:
+        from slycot import sb04qd
+    except ImportError:
         raise ControlSlycot("can't find slycot module 'sb04qd'")
-    if sg03ad is None:
+
+    try:
+        from slycot import sg03ad
+    except ImportError:
         raise ControlSlycot("can't find slycot module 'sg03ad'")
 
     # Reshape 1-d arrays
diff --git a/control/statefbk.py b/control/statefbk.py
index 0017412a..7106449a 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -47,25 +47,6 @@ from .mateqn import care
 from .statesp import _ssmatrix
 from .exception import ControlSlycot, ControlArgument, ControlDimension
 
-# Make sure we have access to the right slycot routines
-try:
-    from slycot import sb03md57
-    # wrap without the deprecation warning
-    def sb03md(n, C, A, U, dico, job='X',fact='N',trana='N',ldwork=None):
-        ret = sb03md57(A, U, C, dico, job, fact, trana, ldwork)
-        return ret[2:]
-except ImportError:
-    try:
-        from slycot import sb03md
-    except ImportError:
-        sb03md = None
-
-try:
-    from slycot import sb03od
-except ImportError:
-    sb03od = None
-
-
 __all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 'lqe',
            'acker']
 
@@ -593,7 +574,7 @@ def gram(sys, type):
         * if `type` is not 'c', 'o', 'cf' or 'of'
         * if system is unstable (sys.A has eigenvalues not in left half plane)
 
-    ControlSlycot
+    ImportError
         if slycot routine sb03md cannot be found
         if slycot routine sb03od cannot be found
 
@@ -633,7 +614,9 @@ def gram(sys, type):
     if type == 'c' or type == 'o':
         # Compute Gramian by the Slycot routine sb03md
         # make sure Slycot is installed
-        if sb03md is None:
+        try:
+            from slycot import sb03md
+        except ImportError:
             raise ControlSlycot("can't find slycot module 'sb03md'")
         if type == 'c':
             tra = 'T'
@@ -651,7 +634,9 @@ def gram(sys, type):
 
     elif type == 'cf' or type == 'of':
         # Compute cholesky factored gramian from slycot routine sb03od
-        if sb03od is None:
+        try:
+            from slycot import sb03od
+        except ImportError:
             raise ControlSlycot("can't find slycot module 'sb03od'")
         tra = 'N'
         n = sys.nstates
diff --git a/control/statesp.py b/control/statesp.py
index 758b91ed..c75e6f66 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -1216,104 +1216,27 @@ class StateSpace(LTI):
 
         .. math: G(1) = C (I - A)^{-1} B + D
 
-        Returns
-        -------
-        gain : ndarray
-            An array of shape (outputs,inputs); the array will either be the
-            zero-frequency (or DC) gain, or, if the frequency response is
-            singular, the array will be filled with (inf + nanj).
-
-        """
-        return self(0, warn_infinite=warn_infinite) if self.isctime() \
-            else self(1, warn_infinite=warn_infinite)
-
-    def dynamics(self, t, x, u=0):
-        """Compute the dynamics of the system
-
-        Given input `u` and state `x`, returns the dynamics of the state-space
-        system. If the system is continuous, returns the time derivative dx/dt
-
-            dx/dt = A x + B u
-
-        where A and B are the state-space matrices of the system. If the
-        system is discrete-time, returns the next value of `x`:
-
-            x[t+dt] = A x[t] + B u[t]
-
-        The inputs `x` and `u` must be of the correct length for the system.
-
-        The first argument `t` is ignored because :class:`StateSpace` systems
-        are time-invariant. It is included so that the dynamics can be passed
-        to most numerical integrators, such as :func:`scipy.integrate.solve_ivp`
-        and for consistency with :class:`IOSystem` systems.
-
-        Parameters
-        ----------
-        t : float (ignored)
-            time
-        x : array_like
-            current state
-        u : array_like (optional)
-            input, zero if omitted
-
-        Returns
-        -------
-        dx/dt or x[t+dt] : ndarray
-        """
-        x = np.reshape(x, (-1, 1)) # force to a column in case matrix
-        if np.size(x) != self.nstates:
-            raise ValueError("len(x) must be equal to number of states")
-        if u is 0:
-            return self.A.dot(x).reshape((-1,)) # return as row vector
-        else: # received t, x, and u, ignore t
-            u = np.reshape(u, (-1, 1)) # force to a column in case matrix
-            if np.size(u) != self.ninputs:
-                raise ValueError("len(u) must be equal to number of inputs")
-            return self.A.dot(x).reshape((-1,)) \
-                 + self.B.dot(u).reshape((-1,)) # return as row vector
-
-    def output(self, t, x, u=0):
-        """Compute the output of the system
-
-        Given input `u` and state `x`, returns the output `y` of the
-        state-space system:
-
-            y = C x + D u
-
-        where A and B are the state-space matrices of the system.
-
-        The first argument `t` is ignored because :class:`StateSpace` systems
-        are time-invariant. It is included so that the dynamics can be passed
-        to most numerical integrators, such as scipy's `integrate.solve_ivp` and
-        for consistency with :class:`IOSystem` systems.
-
-        The inputs `x` and `u` must be of the correct length for the system.
-
         Parameters
         ----------
-        t : float (ignored)
-            time
-        x : array_like
-            current state
-        u : array_like (optional)
-            input (zero if omitted)
+        warn_infinite : bool, optional
+            By default, don't issue a warning message if the zero-frequency
+            gain is infinite.  Setting `warn_infinite` to generate the warning
+            message.
 
         Returns
         -------
-        y : ndarray
+        gain : (outputs, inputs) ndarray or scalar
+            Array or scalar value for SISO systems, depending on
+            config.defaults['control.squeeze_frequency_response'].
+            The value of the array elements or the scalar is either the
+            zero-frequency (or DC) gain, or `inf`, if the frequency response
+            is singular.
+
+            For real valued systems, the empty imaginary part of the
+            complex zero-frequency response is discarded and a real array or
+            scalar is returned.
         """
-        x = np.reshape(x, (-1, 1)) # force to a column in case matrix
-        if np.size(x) != self.nstates:
-            raise ValueError("len(x) must be equal to number of states")
-
-        if u is 0:
-            return self.C.dot(x).reshape((-1,)) # return as row vector
-        else: # received t, x, and u, ignore t
-            u = np.reshape(u, (-1, 1)) # force to a column in case matrix
-            if np.size(u) != self.ninputs:
-                raise ValueError("len(u) must be equal to number of inputs")
-            return self.C.dot(x).reshape((-1,)) \
-                 + self.D.dot(u).reshape((-1,)) # return as row vector
+        return self._dcgain(warn_infinite)
 
     def _isstatic(self):
         """True if and only if the system has no dynamics, that is,
@@ -1321,7 +1244,6 @@ class StateSpace(LTI):
         return not np.any(self.A) and not np.any(self.B)
 
 
-
 # TODO: add discrete time check
 def _convert_to_statespace(sys, **kw):
     """Convert a system to state space form (if needed).
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 50e4870a..3e48c7f2 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -1070,12 +1070,19 @@ class TransferFunction(LTI):
 
         Returns
         -------
-        gain : ndarray
-            The zero-frequency gain
+        gain : (outputs, inputs) ndarray or scalar
+            Array or scalar value for SISO systems, depending on
+            config.defaults['control.squeeze_frequency_response'].
+            The value of the array elements or the scalar is either the
+            zero-frequency (or DC) gain, or `inf`, if the frequency response
+            is singular.
+
+            For real valued systems, the empty imaginary part of the
+            complex zero-frequency response is discarded and a real array or
+            scalar is returned.
 
         """
-        return self(0, warn_infinite=warn_infinite) if self.isctime() \
-            else self(1, warn_infinite=warn_infinite)
+        return self._dcgain(warn_infinite)
 
     def _isstatic(self):
          """returns True if and only if all of the numerator and denominator


Test patch:
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index 983330af..321580ba 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -433,9 +433,9 @@ def test_dcgain_consistency():
     np.testing.assert_equal(
         sys_ss(0j, warn_infinite=False), complex(np.inf, np.nan))
     np.testing.assert_equal(
-        sys_tf.dcgain(warn_infinite=False), complex(np.inf, np.nan))
+        sys_tf.dcgain(), np.inf)
     np.testing.assert_equal(
-        sys_ss.dcgain(warn_infinite=False), complex(np.inf, np.nan))
+        sys_ss.dcgain(), np.inf)
 
     # Set up transfer function with pole, zero at the origin
     sys_tf = ctrl.tf([1, 0], [1, 0])
@@ -448,7 +448,7 @@ def test_dcgain_consistency():
     np.testing.assert_equal(
         sys_tf(0j, warn_infinite=False), complex(np.nan, np.nan))
     np.testing.assert_equal(
-        sys_tf.dcgain(warn_infinite=False), complex(np.nan, np.nan))
+        sys_tf.dcgain(), np.nan)
 
     # Set up state space version
     sys_ss = ctrl.tf2ss(ctrl.tf([1, 0], [1, 1])) * \
@@ -462,7 +462,7 @@ def test_dcgain_consistency():
         np.testing.assert_equal(
             sys_ss(0j, warn_infinite=False), complex(np.nan, np.nan))
         np.testing.assert_equal(
-            sys_ss.dcgain(warn_infinite=False), complex(np.nan, np.nan))
+            sys_ss.dcgain(), np.nan)
     elif 0 in sys_ss.pole():
         # Pole at the origin, but zero elsewhere => should get (inf + nanj)
         np.testing.assert_equal(
@@ -470,7 +470,7 @@ def test_dcgain_consistency():
         np.testing.assert_equal(
             sys_ss(0j, warn_infinite=False), complex(np.inf, np.nan))
         np.testing.assert_equal(
-            sys_ss.dcgain(warn_infinite=False), complex(np.inf, np.nan))
+            sys_ss.dcgain(), np.inf)
     else:
         # Near pole/zero cancellation => nothing sensible to check
         pass
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
index a1246103..afeb71ef 100644
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -336,41 +336,41 @@ def test_zmore_stability_margins(tsys_zmore):
 @pytest.mark.parametrize(
     'cnum, cden, dt,'
     'ref,'
-    'rtol, poly_is_inaccurate',
+    'rtol',
     [( # gh-465
       [2], [1, 3, 2, 0], 1e-2,
-      [ 2.955761, 32.398492,  0.429535,  1.403725,  0.749367,  0.923898],
-      1e-5, True),
+      [2.9558, 32.390, 0.43584, 1.4037, 0.74951, 0.97079],
+      2e-3), # the gradient of the function reduces numerical precision
      ( # 2/(s+1)**3
       [2], [1, 3, 3, 1], .1,
       [3.4927, 65.4212, 0.5763, 1.6283, 0.76625, 1.2019],
-      1e-4, True),
-     ( # gh-523 a
+      1e-4),
+     ( # gh-523
       [1.1 * 4 * np.pi**2], [1, 2 * 0.2 * 2 * np.pi,  4 * np.pi**2], .05,
       [2.3842, 18.161, 0.26953, 11.712, 8.7478, 9.1504],
-      1e-4, False),
-     ( # gh-523 b
-       # H1 = w1**2 / (z**2 + 2*zt*w1 * z + w1**2)
-       # H2 = w2**2 / (z**2 + 2*zt*w2 * z + w2**2)
-       # H = H1 * H2
-       # w1 = 1, w2 = 100, zt = 0.5
-      [5e4], [1., 101., 10101., 10100., 10000.], 1e-3,
-      [18.8766, 26.3564, 0.406841, 9.76358, 2.32933, 2.55986],
-      1e-5, True),
+      1e-4),
      ])
-@pytest.mark.filterwarnings("error")
-def test_stability_margins_discrete(cnum, cden, dt,
-                                    ref,
-                                    rtol, poly_is_inaccurate):
+def test_stability_margins_discrete(cnum, cden, dt, ref, rtol):
     """Test stability_margins with discrete TF input"""
     tf = TransferFunction(cnum, cden).sample(dt)
-    if poly_is_inaccurate:
-        with pytest.warns(UserWarning, match="numerical inaccuracy in 'poly'"):
-            out = stability_margins(tf)
-        # cover the explicit frd branch and make sure it yields the same
-        # results as the fallback mechanism
-        out_frd = stability_margins(tf, method='frd')
-        assert_allclose(out, out_frd)
-    else:
-        out = stability_margins(tf)
+    out = stability_margins(tf, method='poly')
     assert_allclose(out, ref, rtol=rtol)
+
+
+def test_stability_margins_methods():
+    # the following system gives slightly inaccurate result for DT systems
+    # because of numerical issues
+    omegan = 1
+    zeta = 0.5
+    resonance = TransferFunction(omegan**2, [1, 2*zeta*omegan, omegan**2])
+    omegan2 = 100
+    resonance2 = TransferFunction(omegan2**2, [1, 2*zeta*omegan2, omegan2**2])
+    sys = 5 * resonance * resonance2
+    sysd = sys.sample(0.001, 'zoh')
+    """Test stability_margins() function with different methods"""
+    out = stability_margins(sysd, method='best')
+    # confirm getting reasonable results using FRD method
+    assert_allclose(
+        (18.876634740386308, 26.356358386241055, 0.40684127995261044,
+         9.763585494645046, 2.3293357226374805, 2.55985695034263),
+        stability_margins(sysd, method='frd'), rtol=1e-5)
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 1eec5ead..6a750900 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -8,7 +8,6 @@ BG,  26 Jul 2020 merge statesp_array_test.py differences into statesp_test.py
 """
 
 import numpy as np
-from numpy.testing import assert_array_almost_equal
 import pytest
 import operator
 from numpy.linalg import solve
@@ -48,17 +47,6 @@ class TestStateSpace:
         """3-states square system (2 inputs x 2 outputs)"""
         return StateSpace(*sys322ABCD)
 
-    @pytest.fixture
-    def sys121(self):
-        """2 state, 1 input, 1 output (siso)  system"""
-        A121 = [[4., 1.],
-                [2., -3]]
-        B121 = [[5.],
-                [-3.]]
-        C121 = [[2., -4]]
-        D121 = [[3.]]
-        return StateSpace(A121, B121, C121, D121)
-
     @pytest.fixture
     def sys222(self):
         """2-states square system (2 inputs x 2 outputs)"""
@@ -510,7 +498,7 @@ class TestStateSpace:
         np.testing.assert_allclose(sys2.dcgain(), expected)
 
         sys3 = StateSpace(0., 1., 1., 0.)
-        np.testing.assert_equal(sys3.dcgain(), complex(np.inf, np.nan))
+        np.testing.assert_equal(sys3.dcgain(), np.inf)
 
     def test_dc_gain_discr(self):
         """Test DC gain for discrete-time state-space systems."""
@@ -528,7 +516,7 @@ class TestStateSpace:
 
         # summer
         sys = StateSpace(1, 1, 1, 0, True)
-        np.testing.assert_equal(sys.dcgain(), complex(np.inf, np.nan))
+        np.testing.assert_equal(sys.dcgain(), np.inf)
 
     @pytest.mark.parametrize("outputs", range(1, 6))
     @pytest.mark.parametrize("inputs", range(1, 6))
@@ -551,7 +539,7 @@ class TestStateSpace:
         c = np.eye(max(outputs, states))[:outputs, :states]
         d = np.zeros((outputs, inputs))
         sys = StateSpace(a, b, c, d, dt)
-        dc = np.full_like(d, complex(np.inf, np.nan), dtype=complex)
+        dc = np.full_like(d, np.inf, dtype=float)
         if sys.issiso():
             dc = dc.squeeze()
 
@@ -763,70 +751,6 @@ class TestStateSpace:
             np.squeeze(sys322.horner(1.j)),
             mag[:, :, 0] * np.exp(1.j * phase[:, :, 0]))
 
-    @pytest.mark.parametrize('x',
-        [[1, 1], [[1], [1]], np.atleast_2d([1,1]).T])
-    @pytest.mark.parametrize('u', [0, 1, np.atleast_1d(2)])
-    def test_dynamics_and_output_siso(self, x, u, sys121):
-        assert_array_almost_equal(
-            sys121.dynamics(0, x, u),
-            sys121.A.dot(x).reshape((-1,)) + sys121.B.dot(u).reshape((-1,)))
-        assert_array_almost_equal(
-            sys121.output(0, x, u),
-            sys121.C.dot(x).reshape((-1,)) + sys121.D.dot(u).reshape((-1,)))
-        assert_array_almost_equal(
-            sys121.dynamics(0, x),
-            sys121.A.dot(x).reshape((-1,)))
-        assert_array_almost_equal(
-            sys121.output(0, x),
-            sys121.C.dot(x).reshape((-1,)))
-
-    # too few and too many states and inputs
-    @pytest.mark.parametrize('x', [0, 1, [], [1, 2, 3], np.atleast_1d(2)])
-    def test_error_x_dynamics_and_output_siso(self, x, sys121):
-        with pytest.raises(ValueError):
-            sys121.dynamics(0, x)
-        with pytest.raises(ValueError):
-            sys121.output(0, x)
-    @pytest.mark.parametrize('u', [[1, 1], np.atleast_1d((2, 2))])
-    def test_error_u_dynamics_output_siso(self, u, sys121):
-        with pytest.raises(ValueError):
-            sys121.dynamics(0, 1, u)
-        with pytest.raises(ValueError):
-            sys121.output(0, 1, u)
-
-    @pytest.mark.parametrize('x',
-        [[1, 1], [[1], [1]], np.atleast_2d([1,1]).T])
-    @pytest.mark.parametrize('u',
-        [[1, 1], [[1], [1]], np.atleast_2d([1,1]).T])
-    def test_dynamics_and_output_mimo(self, x, u, sys222):
-        assert_array_almost_equal(
-            sys222.dynamics(0, x, u),
-            sys222.A.dot(x).reshape((-1,)) + sys222.B.dot(u).reshape((-1,)))
-        assert_array_almost_equal(
-            sys222.output(0, x, u),
-            sys222.C.dot(x).reshape((-1,)) + sys222.D.dot(u).reshape((-1,)))
-        assert_array_almost_equal(
-            sys222.dynamics(0, x),
-            sys222.A.dot(x).reshape((-1,)))
-        assert_array_almost_equal(
-            sys222.output(0, x),
-            sys222.C.dot(x).reshape((-1,)))
-
-    # too few and too many states and inputs
-    @pytest.mark.parametrize('x', [0, 1, [1, 1, 1]])
-    def test_error_x_dynamics_mimo(self, x, sys222):
-        with pytest.raises(ValueError):
-            sys222.dynamics(0, x)
-        with pytest.raises(ValueError):
-            sys222.output(0, x)
-    @pytest.mark.parametrize('u', [1, [1, 1, 1]])
-    def test_error_u_dynamics_mimo(self, u, sys222):
-        with pytest.raises(ValueError):
-            sys222.dynamics(0, (1, 1), u)
-        with pytest.raises(ValueError):
-            sys222.output(0, (1, 1), u)
-
-
 class TestRss:
     """These are tests for the proper functionality of statesp.rss."""
 
@@ -1029,4 +953,3 @@ def test_xferfcn_ndarray_precedence(op, tf, arr):
     ss = ct.tf2ss(tf)
     result = op(arr, ss)
     assert isinstance(result, ct.StateSpace)
-
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index b892655e..06e7fc9d 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -807,7 +807,7 @@ class TestXferFcn:
         np.testing.assert_equal(sys2.dcgain(), 2)
 
         sys3 = TransferFunction(6, [1, 0])
-        np.testing.assert_equal(sys3.dcgain(), complex(np.inf, np.nan))
+        np.testing.assert_equal(sys3.dcgain(), np.inf)
 
         num = [[[15], [21], [33]], [[10], [14], [22]]]
         den = [[[1, 3], [2, 3], [3, 3]], [[1, 5], [2, 7], [3, 11]]]
@@ -827,13 +827,13 @@ class TestXferFcn:
 
         # differencer
         sys = TransferFunction(1, [1, -1], True)
-        np.testing.assert_equal(sys.dcgain(), complex(np.inf, np.nan))
+        np.testing.assert_equal(sys.dcgain(), np.inf)
 
         # differencer, with warning
         sys = TransferFunction(1, [1, -1], True)
         with pytest.warns(RuntimeWarning, match="divide by zero"):
             np.testing.assert_equal(
-                sys.dcgain(warn_infinite=True), complex(np.inf, np.nan))
+                sys.dcgain(warn_infinite=True), np.inf)
 
         # summer
         sys = TransferFunction([1, -1], [1], True)