change infinite gain value to inf (from nan/inf)

murrayrm · murrayrm · commit b67eb7d3c80d · 2021-02-10T21:44:39.000-08:00
diff --git a/control/statesp.py b/control/statesp.py
@@ -668,7 +668,7 @@ def __rdiv__(self, other):
         raise NotImplementedError(
             "StateSpace.__rdiv__ is not implemented yet.")
 
-    def __call__(self, x, squeeze=None):
+    def __call__(self, x, squeeze=None, warn_infinite=True):
         """Evaluate system's transfer function at complex frequency.
 
         Returns the complex frequency response `sys(x)` where `x` is `s` for
@@ -689,6 +689,8 @@ def __call__(self, x, squeeze=None):
             keep all indices (output, input and, if omega is array_like,
             frequency) even if the system is SISO. The default value can be
             set using config.defaults['control.squeeze_frequency_response'].
+        warn_infinite : bool, optional
+            If set to `False`, don't warn if frequency response is infinite.
 
         Returns
         -------
@@ -702,7 +704,7 @@ def __call__(self, x, squeeze=None):
 
         """
         # Use Slycot if available
-        out = self.horner(x)
+        out = self.horner(x, warn_infinite=warn_infinite)
         return _process_frequency_response(self, x, out, squeeze=squeeze)
 
     def slycot_laub(self, x):
@@ -758,7 +760,7 @@ def slycot_laub(self, x):
             out[:, :, kk+1] = result[0] + self.D
         return out
 
-    def horner(self, x):
+    def horner(self, x, warn_infinite=True):
         """Evaluate system's transfer function at complex frequency
         using Laub's or Horner's method.
 
@@ -795,14 +797,22 @@ def horner(self, x):
                 raise ValueError("input list must be 1D")
 
             # Preallocate
-            out = empty((self.noutputs, self.ninputs, len(x_arr)), dtype=complex)
+            out = empty((self.noutputs, self.ninputs, len(x_arr)),
+                        dtype=complex)
 
             #TODO: can this be vectorized?
             for idx, x_idx in enumerate(x_arr):
-                out[:,:,idx] = \
-                    np.dot(self.C,
+                try:
+                    out[:,:,idx] = np.dot(
+                        self.C,
                         solve(x_idx * eye(self.nstates) - self.A, self.B)) \
-                    + self.D
+                        + self.D
+                except LinAlgError:
+                    if warn_infinite:
+                        warn("frequency response is not finite",
+                             RuntimeWarning)
+                    # TODO: check for nan cases
+                    out[:,:,idx] = np.inf
         return out
 
     def freqresp(self, omega):
@@ -1200,7 +1210,7 @@ def dcgain(self):
         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 np.nan.
+            response is singular, the array will be filled with np.inf.
         """
         try:
             if self.isctime():
@@ -1210,7 +1220,7 @@ def dcgain(self):
                 gain = np.squeeze(self.horner(1))
         except LinAlgError:
             # eigenvalue at DC
-            gain = np.tile(np.nan, (self.noutputs, self.ninputs))
+            gain = np.tile(np.inf, (self.noutputs, self.ninputs))
         return np.squeeze(gain)
 
     def _isstatic(self):
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
@@ -498,7 +498,7 @@ def test_dc_gain_cont(self):
         np.testing.assert_allclose(sys2.dcgain(), expected)
 
         sys3 = StateSpace(0., 1., 1., 0.)
-        np.testing.assert_equal(sys3.dcgain(), 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."""
@@ -516,7 +516,7 @@ def test_dc_gain_discr(self):
 
         # summer
         sys = StateSpace(1, 1, 1, 0, True)
-        np.testing.assert_equal(sys.dcgain(), np.nan)
+        np.testing.assert_equal(sys.dcgain(), np.inf)
 
     @pytest.mark.parametrize("outputs", range(1, 6))
     @pytest.mark.parametrize("inputs", range(1, 6))
@@ -539,7 +539,7 @@ def test_dc_gain_integrator(self, outputs, inputs, dt):
         c = np.eye(max(outputs, states))[:outputs, :states]
         d = np.zeros((outputs, inputs))
         sys = StateSpace(a, b, c, d, dt)
-        dc = np.squeeze(np.full_like(d, np.nan))
+        dc = np.squeeze(np.full_like(d, np.inf))
         np.testing.assert_array_equal(dc, sys.dcgain())
 
     def test_scalar_static_gain(self):
diff --git a/control/xferfcn.py b/control/xferfcn.py
@@ -234,7 +234,7 @@ def __init__(self, *args, **kwargs):
                     dt = config.defaults['control.default_dt']
         self.dt = dt
 
-    def __call__(self, x, squeeze=None):
+    def __call__(self, x, squeeze=None, warn_infinite=True):
         """Evaluate system's transfer function at complex frequencies.
 
         Returns the complex frequency response `sys(x)` where `x` is `s` for
@@ -262,6 +262,8 @@ def __call__(self, x, squeeze=None):
             If True and the system is single-input single-output (SISO),
             return a 1D array rather than a 3D array.  Default value (True)
             set by config.defaults['control.squeeze_frequency_response'].
+        warn_infinite : bool, optional
+            If set to `False`, turn off divide by zero warning.
 
         Returns
         -------
@@ -274,10 +276,10 @@ def __call__(self, x, squeeze=None):
             then single-dimensional axes are removed.
 
         """
-        out = self.horner(x)
+        out = self.horner(x, warn_infinite=warn_infinite)
         return _process_frequency_response(self, x, out, squeeze=squeeze)
 
-    def horner(self, x):
+    def horner(self, x, warn_infinite=True):
         """Evaluate system's transfer function at complex frequency
         using Horner's method.
 
@@ -307,8 +309,9 @@ def horner(self, x):
         out = empty((self.noutputs, self.ninputs, len(x_arr)), dtype=complex)
         for i in range(self.noutputs):
             for j in range(self.ninputs):
-                out[i][j] = (polyval(self.num[i][j], x) /
-                             polyval(self.den[i][j], x))
+                with np.errstate(divide='warn' if warn_infinite else 'ignore'):
+                    out[i][j] = (polyval(self.num[i][j], x) /
+                                 polyval(self.den[i][j], x))
         return out
 
     def _truncatecoeff(self):
