updated Nichols chart code to match 0.3d

murrayrm · murrayrm · commit 0f31ff6fb628 · 2011-04-02T16:15:41.000Z
diff --git a/ChangeLog b/ChangeLog
@@ -1,3 +1,15 @@
+2011-04-02  Richard Murray  <murray@malabar.local>
+
+	* tests/nichols_test.py (TestStateSpace.testNgrid): updated testcode
+	to turn off grid in initial Nichols chart plot.
+
+	* src/freqplot.py: updated comments at top of file to reflect
+	nichols chart move
+
+	* src/nichols.py: transferred over changes from v0.3d
+
+	* src/matlab.py (ngrid): moved import to function
+
 2011-03-31  Richard Murray  <murray@malabar.local>
 
 	* examples/pvtol-nested.py: updated stability margin plot to use
diff --git a/src/freqplot.py b/src/freqplot.py
@@ -4,7 +4,8 @@
 # Date: 24 May 09
 # 
 # This file contains some standard control system plots: Bode plots,
-# Nyquist plots, Nichols plots and pole-zero diagrams
+# Nyquist plots and pole-zero diagrams.  The code for Nichols charts
+# is in nichols.py.
 #
 # Copyright (c) 2010 by California Institute of Technology
 # All rights reserved.
diff --git a/src/matlab.py b/src/matlab.py
@@ -75,7 +75,7 @@
 # Import MATLAB-like functions that can be used as-is
 from ctrlutil import unwrap
 from freqplot import nyquist, gangof4
-from nichols import nichols, nichols_grid
+from nichols import nichols
 from bdalg import series, parallel, negate, feedback
 from pzmap import pzmap
 from statefbk import ctrb, obsv, gram, place, lqr
@@ -753,6 +753,7 @@ def ngrid():
     =====
     ngrid()
     """
+    from nichols import nichols_grid
     nichols_grid()
 
 #
diff --git a/src/nichols.py b/src/nichols.py
@@ -45,7 +45,8 @@
 from ctrlutil import unwrap
 from freqplot import default_frequency_range
 
-def nichols(syslist, omega=None):
+# Nichols plot
+def nichols(syslist, omega=None, grid=True):
     """Nichols plot for a system
 
     Usage
@@ -60,6 +61,8 @@ def nichols(syslist, omega=None):
         List of linear input/output systems (single system is OK)
     omega : freq_range
         Range of frequencies (list or bounds) in rad/sec
+    grid : boolean, optional
+        True if the plot should include a Nichols-chart grid. Default is True.
 
     Return values
     -------------
@@ -71,7 +74,7 @@ def nichols(syslist, omega=None):
         syslist = (syslist,)
 
     # Select a default range if none is provided
-    if (omega == None):
+    if omega is None:
         omega = default_frequency_range(syslist)
 
     for sys in syslist:
@@ -94,88 +97,110 @@ def nichols(syslist, omega=None):
 
     # Mark the -180 point
     plt.plot([-180], [0], 'r+')
+
+    # Add grid
+    if grid:
+        nichols_grid()
     
 # Nichols grid
-def nichols_grid():
+#! TODO: Consider making linestyle configurable
+def nichols_grid(cl_mags=None, cl_phases=None):
     """Nichols chart grid
     
     Usage
     =====
     nichols_grid()
 
-    Plots a Nichols chart grid on the current axis.
+    Plots a Nichols chart grid on the current axis, or creates a new chart
+    if no plot already exists.
 
     Parameters
     ----------
-    None
+    cl_mags : array-like (dB), optional
+        Array of closed-loop magnitudes defining the iso-gain lines on a
+        custom Nichols chart.
+    cl_phases : array-like (degrees), optional
+        Array of closed-loop phases defining the iso-phase lines on a custom
+        Nichols chart. Must be in the range -360 < cl_phases < 0
 
     Return values
     -------------
     None    
     """
-    mag_min_default = -40.0 # dB
-    mag_step = 20.0         # dB
-
-    # Chart defaults
-    phase_min, phase_max, mag_min, mag_max = -360.0, 0.0, mag_min_default, 40.0
-
-    # Set actual chart bounds based on current plot
+    # Default chart size
+    ol_phase_min = -359.99
+    ol_phase_max = 0.0
+    ol_mag_min = -40.0
+    ol_mag_max = default_ol_mag_max = 50.0
+    
+    # Find bounds of the current dataset, if there is one. 
     if plt.gcf().gca().has_data():
-        phase_min, phase_max, mag_min, mag_max = plt.axis()
-        
+        ol_phase_min, ol_phase_max, ol_mag_min, ol_mag_max = plt.axis()
+
     # M-circle magnitudes.
-    # The "fixed" set are always generated, since this guarantees a recognizable
-    # Nichols chart grid.
-    mags_fixed = np.array([-40.0, -20.0, -12.0, -6.0, -3.0, -1.0, -0.5, 0.0,
-                                 0.25, 0.5, 1.0, 3.0, 6.0, 12.0])
-
-    if mag_min < mag_min_default:
-        # Outside the "fixed" set of magnitudes, the generated M-circles
-        # are extended in steps of 'mag_step' dB to cover anything made
-        # visible by the range of the existing plot
-        mags_adjust = np.arange(mag_step*np.floor(mag_min/mag_step),
-                                mag_min_default, mag_step)
-        mags = np.concatenate((mags_adjust, mags_fixed))
-    else:
-        mags = mags_fixed
+    if cl_mags is None:
+        # Default chart magnitudes
+        # The key set of magnitudes are always generated, since this
+        # guarantees a recognizable Nichols chart grid.
+        key_cl_mags = np.array([-40.0, -20.0, -12.0, -6.0, -3.0, -1.0, -0.5, 0.0,
+                                     0.25, 0.5, 1.0, 3.0, 6.0, 12.0])
+        # Extend the range of magnitudes if necessary. The extended arange
+        # will end up empty if no extension is required. Assumes that closed-loop
+        # magnitudes are approximately aligned with open-loop magnitudes beyond
+        # the value of np.min(key_cl_mags)
+        cl_mag_step = -20.0 # dB
+        extended_cl_mags = np.arange(np.min(key_cl_mags),
+                                     ol_mag_min + cl_mag_step, cl_mag_step)
+        cl_mags = np.concatenate((extended_cl_mags, key_cl_mags))
                      
     # N-circle phases (should be in the range -360 to 0)
-    phases = np.array([-0.25, -10.0, -20.0, -30.0, -45.0, -60.0, -90.0,
-                       -120.0, -150.0, -180.0, -210.0, -240.0, -270.0,
-                       -310.0, -325.0, -340.0, -350.0, -359.75])
+    if cl_phases is None:
+        # Choose a reasonable set of default phases (denser if the open-loop
+        # data is restricted to a relatively small range of phases).
+        key_cl_phases = np.array([-0.25, -45.0, -90.0, -180.0, -270.0, -325.0, -359.75])
+        if np.abs(ol_phase_max - ol_phase_min) < 90.0:
+            other_cl_phases = np.arange(-10.0, -360.0, -10.0)
+        else:
+            other_cl_phases = np.arange(-10.0, -360.0, -20.0)
+        cl_phases = np.concatenate((key_cl_phases, other_cl_phases))
+    else:
+        assert ((-360.0 < np.min(cl_phases)) and (np.max(cl_phases) < 0.0))
 
     # Find the M-contours
-    m = m_circles(mags, phase_min=np.min(phases), phase_max=np.max(phases))
+    m = m_circles(cl_mags, phase_min=np.min(cl_phases), phase_max=np.max(cl_phases))
     m_mag = 20*sp.log10(np.abs(m))
     m_phase = sp.mod(sp.degrees(sp.angle(m)), -360.0) # Unwrap
 
     # Find the N-contours
-    n = n_circles(phases, mag_min=np.min(mags), mag_max=np.max(mags))
+    n = n_circles(cl_phases, mag_min=np.min(cl_mags), mag_max=np.max(cl_mags))
     n_mag = 20*sp.log10(np.abs(n))
     n_phase = sp.mod(sp.degrees(sp.angle(n)), -360.0) # Unwrap
 
     # Plot the contours behind other plot elements.
     # The "phase offset" is used to produce copies of the chart that cover
     # the entire range of the plotted data, starting from a base chart computed
-    # over the range -360 < phase < 0 (see above). Given the range 
+    # over the range -360 < phase < 0. Given the range 
     # the base chart is computed over, the phase offset should be 0
-    # for -360 < phase_min < 0.
-    phase_offset_min = 360.0*np.ceil(phase_min/360.0)
-    phase_offset_max = 360.0*np.ceil(phase_max/360.0) + 360.0
+    # for -360 < ol_phase_min < 0.
+    phase_offset_min = 360.0*np.ceil(ol_phase_min/360.0)
+    phase_offset_max = 360.0*np.ceil(ol_phase_max/360.0) + 360.0
     phase_offsets = np.arange(phase_offset_min, phase_offset_max, 360.0)
+
     for phase_offset in phase_offsets:
+        # Draw M and N contours
         plt.plot(m_phase + phase_offset, m_mag, color='gray',
-                 linestyle='dashed', zorder=0)
+                 linestyle='dotted', zorder=0)
         plt.plot(n_phase + phase_offset, n_mag, color='gray',
-                 linestyle='dashed', zorder=0)
+                 linestyle='dotted', zorder=0)
 
-    # Add magnitude labels
-    for x, y, m in zip(m_phase[:][-1], m_mag[:][-1], mags):
-        align = 'right' if m < 0.0 else 'left'
-        plt.text(x, y, str(m) + ' dB', size='small', ha=align)
+        # Add magnitude labels
+        for x, y, m in zip(m_phase[:][-1] + phase_offset, m_mag[:][-1], cl_mags):
+            align = 'right' if m < 0.0 else 'left'
+            plt.text(x, y, str(m) + ' dB', size='small', ha=align, color='gray')
 
-    # Make sure axes conform to any pre-existing plot.
-    plt.axis([phase_min, phase_max, mag_min, mag_max])
+    # Fit axes to generated chart
+    plt.axis([phase_offset_min - 360.0, phase_offset_max - 360.0,
+              np.min(cl_mags), np.max([ol_mag_max, default_ol_mag_max])])
 
 #
 # Utility functions
@@ -185,38 +210,44 @@ def nichols_grid():
 #
    
 # Compute contours of a closed-loop transfer function
-def closed_loop_contours(Hmag, Hphase):
-    """Contours of the function H = G/(1+G).
+def closed_loop_contours(Gcl_mags, Gcl_phases):
+    """Contours of the function Gcl = Gol/(1+Gol), where
+    Gol is an open-loop transfer function, and Gcl is a corresponding
+    closed-loop transfer function.
 
     Usage
     =====
-    contours = closed_loop_contours(mags, phases)
+    contours = closed_loop_contours(Gcl_mags, Gcl_phases)
 
     Parameters
     ----------
-    mags : array-like
-        Meshgrid array of magnitudes of the contours
-    phases : array-like
-        Meshgrid array of phases in radians of the contours
+    Gcl_mags : array-like
+        Array of magnitudes of the contours
+    Gcl_phases : array-like
+        Array of phases in radians of the contours
 
     Return values
     -------------
     contours : complex array
         Array of complex numbers corresponding to the contours.
     """
-    # Compute the contours in H-space
-    H = Hmag*sp.exp(1.j*Hphase)
+    # Compute the contours in Gcl-space. Since we're given closed-loop
+    # magnitudes and phases, this is just a case of converting them into
+    # a complex number.
+    Gcl = Gcl_mags*sp.exp(1.j*Gcl_phases)
 
-    # Invert H = G/(1+G) to get an expression for the contours in G-space
-    return H/(1.0 - H)
+    # Invert Gcl = Gol/(1+Gol) to map the contours into the open-loop space
+    return Gcl/(1.0 - Gcl)
 
 # M-circle
 def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
-    """Constant-magnitude contours of the function H = G/(1+G).
+    """Constant-magnitude contours of the function Gcl = Gol/(1+Gol), where
+    Gol is an open-loop transfer function, and Gcl is a corresponding
+    closed-loop transfer function.
 
     Usage
     =====
-    contours = m_circles(mags)
+    contours = m_circles(mags, phase_min, phase_max)
 
     Parameters
     ----------
@@ -234,17 +265,19 @@ def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
     """
     # Convert magnitudes and phase range into a grid suitable for
     # building contours
-    phases = sp.radians(sp.linspace(phase_min, phase_max, 500))
-    Hmag, Hphase = sp.meshgrid(10.0**(mags/20.0), phases)
-    return closed_loop_contours(Hmag, Hphase)
+    phases = sp.radians(sp.linspace(phase_min, phase_max, 2000))
+    Gcl_mags, Gcl_phases = sp.meshgrid(10.0**(mags/20.0), phases)
+    return closed_loop_contours(Gcl_mags, Gcl_phases)
 
 # N-circle
 def n_circles(phases, mag_min=-40.0, mag_max=12.0):
-    """Constant-phase contours of the function H = G/(1+G).
+    """Constant-phase contours of the function Gcl = Gol/(1+Gol), where
+    Gol is an open-loop transfer function, and Gcl is a corresponding
+    closed-loop transfer function.
 
     Usage
     ===== 
-    contour = n_circles(angles)
+    contours = n_circles(phases, mag_min, mag_max)
 
     Parameters
     ----------
@@ -263,5 +296,5 @@ def n_circles(phases, mag_min=-40.0, mag_max=12.0):
     # Convert phases and magnitude range into a grid suitable for
     # building contours    
     mags = sp.linspace(10**(mag_min/20.0), 10**(mag_max/20.0), 2000)
-    Hphase, Hmag = sp.meshgrid(sp.radians(phases), mags)
-    return closed_loop_contours(Hmag, Hphase)
+    Gcl_phases, Gcl_mags = sp.meshgrid(sp.radians(phases), mags)
+    return closed_loop_contours(Gcl_mags, Gcl_phases)
diff --git a/tests/nichols_test.py b/tests/nichols_test.py
@@ -20,13 +20,13 @@ def setUp(self):
         
         self.sys = StateSpace(A, B, C, D) 
 
-    def testNichols(self):
+    def testNicholsPlain(self):
         """Generate a Nichols plot."""
         nichols(self.sys)
 
     def testNgrid(self):
         """Generate a Nichols plot."""
-        nichols(self.sys)
+        nichols(self.sys, grid=False)
         ngrid()
 
 def suite():

Original file line number	Diff line number	Diff line change
`@@ -4,7 +4,8 @@`
`4`	`4`	`# Date: 24 May 09`
`5`	`5`	`#`
`6`	`6`	`# This file contains some standard control system plots: Bode plots,`
`7`		`-# Nyquist plots, Nichols plots and pole-zero diagrams`
	`7`	`+# Nyquist plots and pole-zero diagrams. The code for Nichols charts`
	`8`	`+# is in nichols.py.`
`8`	`9`	`#`
`9`	`10`	`# Copyright (c) 2010 by California Institute of Technology`
`10`	`11`	`# All rights reserved.`