Repository: python-control/python-control
Pull request #997 (2024-04-24 11:35:54)

Add extra information to nyquist plots
Adds the possibility to add the unit circle (to get gain crossover frequency),  Ms-circle (sensitivity) and Mt-circle (complementary sensitivity) to nyquist plots similar to the approach taken in [JuliaControl](https://github.com/JuliaControl/ControlSystems.jl/blob/master/lib/ControlSystemsBase/src/plotting.jl)

![image](https://github.com/python-control/python-control/assets/74975295/66b36818-6166-41e2-893c-3e2a3928942f)


murrayrm: It would be nice to have unit tests, even if they call the code to get coverage (better would be to check to make sure the lines are in the right place).

If the system is discrete time, are these still the right circles to draw?

Any interaction with axis scaling?  Eg, if all of the action is very near the origin, might this reset the axes in a way that makes it harder to see the Nyquist plot?  If so, is there a better way to do that (eg, fix the axes and let the circle disappear at the edges of the axes)?

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

coverage: 94.405% (-0.1%) from 94.503%
when pulling **ba9dbf0e93a8a8b2985d1ce5267ac3167952bf03 on wueestry:nyquist_improvements**
into **e4a03e8419231b040d846ef089babd029e18168c on python-control:main**.

wueestry: > It would be nice to have unit tests, even if they call the code to get coverage (better would be to check to make sure the lines are in the right place).
> 
> If the system is discrete time, are these still the right circles to draw?
> 
> Any interaction with axis scaling? Eg, if all of the action is very near the origin, might this reset the axes in a way that makes it harder to see the Nyquist plot? If so, is there a better way to do that (eg, fix the axes and let the circle disappear at the edges of the axes)?

This works as well for discrete systems. And tests have now been added.
Not sure yet if fixing the axis is the best idea, as it would result in issues when plotting a second system in the same plot afterwards

Patch:
diff --git a/control/freqplot.py b/control/freqplot.py
index a63ef20d..8917690b 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -1060,7 +1060,9 @@ _nyquist_defaults = {
     'nyquist.max_curve_magnitude': 20,          # clip large values
     'nyquist.max_curve_offset': 0.02,           # offset of primary/mirror
     'nyquist.start_marker': 'o',                # marker at start of curve
-    'nyquist.start_marker_size': 4,             # size of the maker
+    'nyquist.start_marker_size': 4,             # size of the marker
+    'nyquist.circle_style':                     # style for unit circles
+      {'color': 'black', 'linestyle': 'dashed', 'linewidth': 1}
 }
 
 
@@ -1477,8 +1479,8 @@ def nyquist_response(
 
 def nyquist_plot(
         data, omega=None, plot=None, label_freq=0, color=None, label=None,
-        return_contour=None, title=None, legend_loc='upper right',
-        ax=None, **kwargs):
+        return_contour=None, title=None, legend_loc='upper right', ax=None,
+        unit_circle=False, mt_circles=None, ms_circles=None, **kwargs):
     """Nyquist plot for a system.
 
     Generates a Nyquist plot for the system over a (optional) frequency
@@ -1501,7 +1503,13 @@ def nyquist_plot(
         ``omega_limits``.
     color : string, optional
         Used to specify the color of the line and arrowhead.
-
+    unit_circle : bool, optional
+        If ``True``, display the unit circle, to read gain crossover frequency.
+    mt_circles : array_like, optional
+        Draw circles corresponding to the given magnitudes of sensitivity.
+    ms_circles : array_like, optional
+        Draw circles corresponding to the given magnitudes of complementary
+        sensitivity.
     **kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional
         Additional keywords (passed to `matplotlib`)
 
@@ -1856,7 +1864,48 @@ def nyquist_plot(
         # Mark the -1 point
         plt.plot([-1], [0], 'r+')
 
-        # Label the frequencies of the points
+        #
+        # Draw circles for gain crossover and sensitivity functions
+        #
+        theta = np.linspace(0, 2*np.pi, 100)
+        cos = np.cos(theta)
+        sin = np.sin(theta)
+        label_pos = 15
+
+        # Display the unit circle, to read gain crossover frequency
+        if unit_circle:
+            plt.plot(cos, sin, **config.defaults['nyquist.circle_style'])
+        
+        # Draw circles for given magnitudes of sensitivity
+        if ms_circles is not None:
+            for ms in ms_circles:
+                pos_x = -1 + (1/ms)*cos
+                pos_y = (1/ms)*sin
+                plt.plot(
+                    pos_x, pos_y, **config.defaults['nyquist.circle_style'])
+                plt.text(pos_x[label_pos], pos_y[label_pos], ms)
+
+        # Draw circles for given magnitudes of complementary sensitivity
+        if mt_circles is not None:
+            for mt in mt_circles:
+                if mt != 1:
+                    ct = -mt**2/(mt**2-1)  # Mt center
+                    rt = mt/(mt**2-1)  # Mt radius
+                    pos_x = ct+rt*cos
+                    pos_y = rt*sin
+                    plt.plot(
+                        pos_x, pos_y,
+                        **config.defaults['nyquist.circle_style'])
+                    plt.text(pos_x[label_pos], pos_y[label_pos], mt)
+                else:
+                    _, _, ymin, ymax = plt.axis()
+                    pos_y = np.linspace(ymin, ymax, 100)
+                    plt.vlines(
+                        -0.5, ymin=ymin, ymax=ymax,
+                        **config.defaults['nyquist.circle_style'])
+                    plt.text(-0.5, pos_y[label_pos], 1)
+
+        # Label the frequencies of the points on the Nyquist curve
         if label_freq:
             ind = slice(None, None, label_freq)
             omega_sys = np.imag(splane_contour[np.real(splane_contour) == 0])


Test patch:
diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py
index 8354932d..af950535 100644
--- a/control/tests/nyquist_test.py
+++ b/control/tests/nyquist_test.py
@@ -214,6 +214,22 @@ def test_nyquist_arrows(arrows):
     assert _Z(sys) == response.count + _P(sys)
 
 
+def test_sensitivity_circles():
+    A = np.array([
+        [-3.56355873, -1.22980795, -1.5626527 , -0.4626829],
+        [-8.52361371, -3.60331459, -3.71574266, -0.43839201],
+        [-2.50458726, -0.72361335, -1.77795489, -0.4038419],
+        [-0.281183  ,  0.23391825,  0.19096003, -0.9771515]])
+    B = np.array([[-0.], [-1.42827213], [ 0.76806551], [-1.07987454]])
+    C = np.array([[-0.,  0.35557249,  0.35941791, -0.]])
+    D = np.array([[0]])
+    sys1 = ct.ss(A, B, C, D)
+    sys2 = ct.ss(A, B, C, D, dt=0.1)
+    plt.figure()
+    ct.nyquist_plot(sys1, unit_circle=True, mt_circles=[0.9,1,1.1,1.2], ms_circles=[0.9,1,1.1,1.2])
+    ct.nyquist_plot(sys2, unit_circle=True, mt_circles=[0.9,1,1.1,1.2], ms_circles=[0.9,1,1.1,1.2])
+
+
 def test_nyquist_encirclements():
     # Example 14.14: effect of friction in a cart-pendulum system
     s = ct.tf('s')
@@ -518,6 +534,9 @@ if __name__ == "__main__":
     test_nyquist_arrows(3)
     test_nyquist_arrows([0.1, 0.5, 0.9])
 
+    print("Test sensitivity circles")
+    test_sensitivity_circles()
+
     print("Stability checks")
     test_nyquist_encirclements()