Fixed docstrings per roryyorke plus discrete-time corrections

murrayrm · murrayrm · commit 7407265496ec · 2022-04-16T13:49:13.000-07:00
* changed return type for poles() and zeros() to complex
* updated (missing) discrete-time tests for stability warnings
* changed processing of near imaginary poles to pure imaginary poles
diff --git a/control/freqplot.py b/control/freqplot.py
@@ -644,14 +644,14 @@ def nyquist_plot(
         Linestyles for mirror image of the Nyquist curve.  The first element
         is used for unscaled portions of the Nyquist curve, the second element
         is used for portions that are scaled (using max_curve_magnitude).  If
-        `False` then omit completely.  Default linestyle (['--', '-.']) is
+        `False` then omit completely.  Default linestyle (['--', ':']) is
         determined by config.defaults['nyquist.mirror_style'].
 
     primary_style : [str, str], optional
         Linestyles for primary image of the Nyquist curve.  The first
         element is used for unscaled portions of the Nyquist curve,
         the second element is used for portions that are scaled (using
-        max_curve_magnitude). Default linestyle (['-', ':']) is
+        max_curve_magnitude).  Default linestyle (['-', '-.']) is
         determined by config.defaults['nyquist.mirror_style'].
 
     start_marker : str, optional
@@ -750,6 +750,9 @@ def _parse_linestyle(style_name, allow_false=False):
         if isinstance(style, str):
             # Only one style provided, use the default for the other
             style = [style, _nyquist_defaults['nyquist.' + style_name][1]]
+            warnings.warn(
+                "use of a single string for linestyle will be deprecated "
+                " in a future release", PendingDeprecationWarning)
         if (allow_false and style is False) or \
            (isinstance(style, list) and len(style) == 2):
             return style
@@ -765,7 +768,7 @@ def _parse_linestyle(style_name, allow_false=False):
 
     # Determine the range of frequencies to use, based on args/features
     omega, omega_range_given = _determine_omega_vector(
-        syslist, omega, omega_limits, omega_num)
+        syslist, omega, omega_limits, omega_num, feature_periphery_decades=2)
 
     # If omega was not specified explicitly, start at omega = 0
     if not omega_range_given:
@@ -790,7 +793,7 @@ def _parse_linestyle(style_name, allow_false=False):
 
         # Determine the contour used to evaluate the Nyquist curve
         if sys.isdtime(strict=True):
-            # Transform frequencies in for discrete-time systems
+            # Restrict frequencies for discrete-time systems
             nyquistfrq = math.pi / sys.dt
             if not omega_range_given:
                 # limit up to and including nyquist frequency
@@ -817,12 +820,12 @@ def _parse_linestyle(style_name, allow_false=False):
                 # because we don't need to indent for them
                 zplane_poles = sys.poles()
                 zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)]
-                splane_poles = np.log(zplane_poles)/sys.dt
+                splane_poles = np.log(zplane_poles) / sys.dt
 
                 zplane_cl_poles = sys.feedback().poles()
                 zplane_cl_poles = zplane_cl_poles[
                     ~np.isclose(abs(zplane_poles), 0.)]
-                splane_cl_poles = np.log(zplane_cl_poles)/sys.dt
+                splane_cl_poles = np.log(zplane_cl_poles) / sys.dt
 
             #
             # Check to make sure indent radius is small enough
@@ -851,8 +854,8 @@ def _parse_linestyle(style_name, allow_false=False):
             # See if we should add some frequency points near imaginary poles
             #
             for p in splane_poles:
-                # See if we need to process this pole (skip any that is on
-                # the not near or on the negative omega axis + user override)
+                # See if we need to process this pole (skip if on the negative
+                # imaginary axis or not near imaginary axis + user override)
                 if p.imag < 0 or abs(p.real) > indent_radius or \
                    omega_range_given:
                     continue
@@ -894,13 +897,13 @@ def _parse_linestyle(style_name, allow_false=False):
                         - (s - p).real
 
                     # Figure out which way to offset the contour point
-                    if p.real < 0 or (np.isclose(p.real, 0)
-                                      and indent_direction == 'right'):
+                    if p.real < 0 or (p.real == 0 and
+                                      indent_direction == 'right'):
                         # Indent to the right
                         splane_contour[i] += offset
 
-                    elif p.real > 0 or (np.isclose(p.real, 0)
-                                        and indent_direction == 'left'):
+                    elif p.real > 0 or (p.real == 0 and
+                                         indent_direction == 'left'):
                         # Indent to the left
                         splane_contour[i] -= offset
 
@@ -937,9 +940,21 @@ def _parse_linestyle(style_name, allow_false=False):
         # Nyquist criterion is actually satisfied.
         #
         if isinstance(sys, (StateSpace, TransferFunction)):
-            P = (sys.poles().real > 0).sum() if indent_direction == 'right' \
-                else (sys.poles().real >= 0).sum()
-            Z = (sys.feedback().poles().real >= 0).sum()
+            # Count the number of open/closed loop RHP poles
+            if sys.isctime():
+                if indent_direction == 'right':
+                    P = (sys.poles().real > 0).sum()
+                else:
+                    P = (sys.poles().real >= 0).sum()
+                Z = (sys.feedback().poles().real >= 0).sum()
+            else:
+                if indent_direction == 'right':
+                    P = (np.abs(sys.poles()) > 1).sum()
+                else:
+                    P = (np.abs(sys.poles()) >= 1).sum()
+                Z = (np.abs(sys.feedback().poles()) >= 1).sum()
+
+            # Check to make sure the results make sense; warn if not
             if Z != count + P and warn_encirclements:
                 warnings.warn(
                     "number of encirclements does not match Nyquist criterion;"
@@ -976,7 +991,7 @@ def _parse_linestyle(style_name, allow_false=False):
             # Find the different portions of the curve (with scaled pts marked)
             reg_mask = np.logical_or(
                 np.abs(resp) > max_curve_magnitude,
-                contour.real != 0)
+                splane_contour.real != 0)
             # reg_mask = np.logical_or(
             #     np.abs(resp.real) > max_curve_magnitude,
             #     np.abs(resp.imag) > max_curve_magnitude)
@@ -1508,7 +1523,7 @@ def singular_values_plot(syslist, omega=None,
 
 # Determine the frequency range to be used
 def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num,
-                            Hz=None):
+                            Hz=None, feature_periphery_decades=None):
     """Determine the frequency range for a frequency-domain plot
     according to a standard logic.
 
@@ -1554,9 +1569,9 @@ def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num,
         if omega_limits is None:
             omega_range_given = False
             # Select a default range if none is provided
-            omega_out = _default_frequency_range(syslist,
-                                                 number_of_samples=omega_num,
-                                                 Hz=Hz)
+            omega_out = _default_frequency_range(
+                syslist, number_of_samples=omega_num, Hz=Hz,
+                feature_periphery_decades=feature_periphery_decades)
         else:
             omega_limits = np.asarray(omega_limits)
             if len(omega_limits) != 2:
@@ -1640,7 +1655,7 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
 
                 features_ = np.concatenate((sys.poles(), sys.zeros()))
                 # Get rid of poles and zeros on the real axis (imag==0)
-                # * origin and real < 0
+               # * origin and real < 0
                 # * at 1.: would result in omega=0. (logaritmic plot!)
                 toreplace = np.isclose(features_.imag, 0.0) & (
                                     (features_.real <= 0.) |
diff --git a/control/statesp.py b/control/statesp.py
@@ -987,7 +987,8 @@ def freqresp(self, omega):
     def poles(self):
         """Compute the poles of a state space system."""
 
-        return eigvals(self.A) if self.nstates else np.array([])
+        return eigvals(self.A).astype(complex) if self.nstates \
+            else np.array([])
 
     def zeros(self):
         """Compute the zeros of a state space system."""
@@ -1006,8 +1007,9 @@ def zeros(self):
             if nu == 0:
                 return np.array([])
             else:
+                # Use SciPy generalized eigenvalue fucntion
                 return sp.linalg.eigvals(out[8][0:nu, 0:nu],
-                                         out[9][0:nu, 0:nu])
+                                         out[9][0:nu, 0:nu]).astype(complex)
 
         except ImportError:  # Slycot unavailable. Fall back to scipy.
             if self.C.shape[0] != self.D.shape[1]:
@@ -1031,7 +1033,7 @@ def zeros(self):
                                            (0, self.B.shape[1])), "constant")
             return np.array([x for x in sp.linalg.eigvals(L, M,
                                                           overwrite_a=True)
-                             if not isinf(x)])
+                             if not isinf(x)], dtype=complex)
 
     # Feedback around a state space system
     def feedback(self, other=1, sign=-1):
diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py
@@ -350,7 +350,8 @@ def test_linestyle_checks():
 
     # If only one line style is given use, the default value for the other
     # TODO: for now, just make sure the signature works; no correct check yet
-    ct.nyquist_plot(sys, primary_style=':', mirror_style='-.')
+    with pytest.warns(PendingDeprecationWarning, match="single string"):
+        ct.nyquist_plot(sys, primary_style=':', mirror_style='-.')
 
 @pytest.mark.usefixtures("editsdefaults")
 def test_nyquist_legacy():
@@ -363,13 +364,17 @@ def test_nyquist_legacy():
     with pytest.warns(UserWarning, match="indented contour may miss"):
         count = ct.nyquist_plot(sys)
 
-
+def test_discrete_nyquist():
+    # Make sure we can handle discrete time systems with negative poles
+    sys = ct.tf(1, [1, -0.1], dt=1) * ct.tf(1, [1, 0.1], dt=1)
+    ct.nyquist_plot(sys)
+    
 if __name__ == "__main__":
     #
     # Interactive mode: generate plots for manual viewing
     #
-    # Running this script in python (or better ipython) will show a collection of
-    # figures that should all look OK on the screeen.
+    # Running this script in python (or better ipython) will show a
+    # collection of figures that should all look OK on the screeen.
     #
 
     # In interactive mode, turn on ipython interactive graphics
@@ -411,3 +416,13 @@ def test_nyquist_legacy():
               np.array2string(sys.poles(), precision=2, separator=','))
     count = ct.nyquist_plot(sys)
     assert _Z(sys) == count + _P(sys)
+
+    print("Discrete time systems")
+    sys = ct.c2d(sys, 0.01)
+    plt.figure()
+    plt.title("Discrete-time; poles: %s" %
+              np.array2string(sys.poles(), precision=2, separator=','))
+    count = ct.nyquist_plot(sys)
+
+    
+
diff --git a/control/xferfcn.py b/control/xferfcn.py
@@ -798,7 +798,7 @@ def poles(self):
         rts = []
         for d, o in zip(den, denorder):
             rts.extend(roots(d[:o + 1]))
-        return np.array(rts)
+        return np.array(rts).astype(complex)
 
     def zeros(self):
         """Compute the zeros of a transfer function."""
@@ -808,7 +808,7 @@ def zeros(self):
                 "for SISO systems.")
         else:
             # for now, just give zeros of a SISO tf
-            return roots(self.num[0][0])
+            return roots(self.num[0][0]).astype(complex)
 
     def feedback(self, other=1, sign=-1):
         """Feedback interconnection between two LTI objects."""
