Copied over latest changes from v0.3d:

murrayrm · murrayrm · commit 1d0afdf7decd · 2011-02-13T16:49:58.000Z
* Added in latest nichols() code from Allan McInnes (including 13 Feb commits)
* Turned svn:keyword Id back on for source code files
* Set version number to 0.4a
diff --git a/ChangeLog b/ChangeLog
@@ -1,3 +1,18 @@
+2011-02-13  Richard Murray  <murray@sumatra.local>
+
+	* src/*.py: added svn:keywords Id properly
+
+	* src/matlab.py (ngrid): added ngrid() from v0.3d
+
+	* src/freqplot.py (nichols_grid, closed_loop_contours, m_circles,
+	n_circles): copied over changes from Allan McInnes in v0.3d; ngrid()
+	functiality + split out some of the nichols chart code into separate
+	functions
+
+2011-02-12  Richard Murray  <murray@sumatra.local>
+
+	* setup.py: updated version number to 0.4a
+
 2010-11-05  Richard Murray  <murray@sumatra.local>
 
 	* external/yottalab.py: New file containing Roberto Bucher's control
diff --git a/examples/bdalg-matlab.py b/examples/bdalg-matlab.py
@@ -0,0 +1,17 @@
+# bdalg-matlab.py - demonstrate some MATLAB commands for block diagram altebra
+# RMM, 29 May 09
+
+from control.matlab import *    # MATLAB-like functions
+
+# System matrices
+A1 = [[0, 1.], [-4, -1]]
+B1 = [[0], [1.]]
+C1 = [[1., 0]]
+sys1ss = ss(A1, B1, C1, 0)
+sys1tf = ss2tf(sys1ss)
+
+sys2tf = tf([1, 0.5], [1, 5]);
+sys2ss = tf2ss(sys2tf);
+
+# Series composition
+series1 = sys1ss + sys2ss;
diff --git a/examples/test-response.py b/examples/test-response.py
@@ -0,0 +1,18 @@
+# test-response.py - Unit tests for system response functions
+# RMM, 11 Sep 2010
+
+from matplotlib.pyplot import * # Grab MATLAB plotting functions
+from control.matlab import *    # Load the controls systems library
+from scipy import arange        # function to create range of numbers
+
+# Create several systems for testing
+sys1 = tf([1], [1, 2, 1])
+sys2 = tf([1, 1], [1, 1, 0])
+
+# Generate step responses
+(T1a, y1a) = step(sys1)
+(T1b, y1b) = step(sys1, T = arange(0, 10, 0.1))
+(T1c, y1c) = step(sys1, X0 = [1, 0])
+(T2a, y2a) = step(sys2, T = arange(0, 10, 0.1))
+
+plot(T1a, y1a, T1b, y1b, T1c, y1c, T2a, y2a)
diff --git a/examples/test-statefbk.py b/examples/test-statefbk.py
@@ -0,0 +1,27 @@
+# test-statefbk.py - Unit tests for state feedback code
+# RMM, 6 Sep 2010
+
+import numpy as np              # Numerical library
+from scipy import *             # Load the scipy functions
+from control.matlab import *    # Load the controls systems library
+
+# Parameters defining the system
+m = 250.0			# system mass
+k = 40.0			# spring constant
+b = 60.0			# damping constant
+
+# System matrices
+A = matrix([[1, -1, 1.], [1, -k/m, -b/m], [1, 1, 1]])
+B = matrix([[0], [1/m], [1]])
+C = matrix([[1., 0, 1.]])
+sys = ss(A, B, C, 0);
+
+# Controllability
+Wc = ctrb(A, B)
+print "Wc = ", Wc
+
+# Observability
+Wo = obsv(A, C)
+print "Wo = ", Wo
+
+
diff --git a/setup.py b/setup.py
@@ -3,7 +3,7 @@
 from distutils.core import setup
 
 setup(name='control',
-      version='0.3c',
+      version='0.4a',
       description='Python Control Systems Library',
       author='Richard Murray',
       author_email='murray@cds.caltech.edu',
diff --git a/src/__init__.py b/src/__init__.py
@@ -37,7 +37,7 @@
 # OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 # SUCH DAMAGE.
 # 
-# $Id: __init__.py 29 2010-11-06 13:03:55Z murrayrm $
+# $Id$
 
 """Control System Library
 
diff --git a/src/bdalg.py b/src/bdalg.py
@@ -47,7 +47,7 @@
 Date: 24 May 09
 Revised: Kevin K. Chen, Dec 10
 
-$Id: bdalg.py 17 2010-05-29 23:50:52Z murrayrm $
+$Id$
 
 """
 
diff --git a/src/ctrlutil.py b/src/ctrlutil.py
@@ -38,7 +38,7 @@
 # OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 # SUCH DAMAGE.
 # 
-# $Id: ctrlutil.py 21 2010-06-06 17:29:42Z murrayrm $
+# $Id$
 
 # Packages that we need access to
 import scipy as sp
diff --git a/src/delay.py b/src/delay.py
@@ -38,7 +38,7 @@
 # OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 # SUCH DAMAGE.
 # 
-# $Id: pade.py 17 2010-05-29 23:50:52Z murrayrm $
+# $Id$
 
 def pade(T, n=1):
     """ 
diff --git a/src/freqplot.py b/src/freqplot.py
@@ -38,7 +38,7 @@
 # OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 # SUCH DAMAGE.
 # 
-# $Id: freqplot.py 33 2010-11-26 21:59:57Z murrayrm $
+# $Id$
 
 import matplotlib.pyplot as plt
 import scipy as sp
@@ -212,16 +212,15 @@ def nyquist(syslist, omega=None, Plot=True):
                 plt.plot([-1], [0], 'r+')
 
         return x, y, omega
-
 # Nichols plot
 # Contributed by Allan McInnes <Allan.McInnes@canterbury.ac.nz>
 #! TODO: need unit test code
-def nichols(syslist, omega=None, Plot=True):
+def nichols(syslist, omega=None):
     """Nichols plot for a system
 
     Usage
     =====
-    mag, phase, freq = nichols(sys, omega=None, Plot=True)
+    magh = nichols(sys, omega=None)
 
     Plots a Nichols plot for the system over a (optional) frequency range.
 
@@ -231,15 +230,10 @@ def nichols(syslist, omega=None, Plot=True):
         List of linear input/output systems (single system is OK)
     omega : freq_range
         Range of frequencies (list or bounds) in rad/sec
-    Plot : boolean
-        if True, plot the Nichols frequency response 
 
     Return values
     -------------
-    mag : array of magnitudes
-    phase : array of phases
-    freq : array of frequencies
-
+    None
     """
 
     # If argument was a singleton, turn it into a list
@@ -250,34 +244,106 @@ def nichols(syslist, omega=None, Plot=True):
     if (omega == None):
         omega = default_frequency_range(syslist)
 
-
     for sys in syslist:
-        if (sys.inputs > 1 or sys.outputs > 1):
-            #TODO: Add MIMO nichols plots. 
-            raise NotImplementedError("Nichols is currently only implemented for SISO systems.")
-        else:
-            # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega = sys.freqresp(omega)
-            mag = np.squeeze(mag_tmp)
-            phase = np.squeeze(phase_tmp)
- 
-            # Convert to Nichols-plot format (phase in degrees, 
-            # and magnitude in dB)
-            x = unwrap(sp.degrees(phase), 360)
-            y = 20*sp.log10(mag)
+        # Get the magnitude and phase of the system
+        mag, phase, omega = sys.freqresp(omega)
     
-            if (Plot):
-                # Generate the plot
-                plt.plot(x, y)
+        # Convert to Nichols-plot format (phase in degrees, 
+        # and magnitude in dB)
+        x = unwrap(sp.degrees(phase), 360)
+        y = 20*sp.log10(mag)
+    
+        # Generate the plot
+        plt.plot(x, y)
         
-                plt.xlabel('Phase (deg)')
-                plt.ylabel('Magnitude (dB)')
-                plt.title('Nichols Plot')
+    plt.xlabel('Phase (deg)')
+    plt.ylabel('Magnitude (dB)')
+    plt.title('Nichols Plot')
 
-                # Mark the -180 point
-                plt.plot([-180], [0], 'r+')
+    # Mark the -180 point
+    plt.plot([-180], [0], 'r+')
     
-            return mag, phase, omega
+# Nichols grid
+def nichols_grid():
+    """Nichols chart grid
+    
+    Usage
+    =====
+    nichols_grid()
+
+    Plots a Nichols chart grid on the current axis.
+
+    Parameters
+    ----------
+    None
+
+    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
+    if plt.gcf().gca().has_data():
+        phase_min, phase_max, mag_min, 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
+                     
+    # 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])
+
+    # Find the M-contours
+    m = m_circles(mags, phase_min=np.min(phases), phase_max=np.max(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_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 
+    # 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
+    phase_offsets = np.arange(phase_offset_min, phase_offset_max, 360.0)
+    for phase_offset in phase_offsets:
+        plt.plot(m_phase + phase_offset, m_mag, color='gray',
+                 linestyle='dashed', zorder=0)
+        plt.plot(n_phase + phase_offset, n_mag, color='gray',
+                 linestyle='dashed', 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)
+
+    # Make sure axes conform to any pre-existing plot.
+    plt.axis([phase_min, phase_max, mag_min, mag_max])
 
 # Gang of Four
 #! TODO: think about how (and whether) to handle lists of systems
@@ -395,3 +461,85 @@ def default_frequency_range(syslist):
                         np.ceil(np.max(features))+1)   
                         
     return omega
+
+# Compute contours of a closed-loop transfer function
+def closed_loop_contours(Hmag, Hphase):
+    """Contours of the function H = G/(1+G).
+
+    Usage
+    =====
+    contours = closed_loop_contours(mags, phases)
+
+    Parameters
+    ----------
+    mags : array-like
+        Meshgrid array of magnitudes of the contours
+    phases : array-like
+        Meshgrid 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)
+
+    # Invert H = G/(1+G) to get an expression for the contours in G-space
+    return H/(1.0 - H)
+
+# M-circle
+def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
+    """Constant-magnitude contours of the function H = G/(1+G).
+
+    Usage
+    =====
+    contours = m_circles(mags)
+
+    Parameters
+    ----------
+    mags : array-like
+        Array of magnitudes in dB of the M-circles
+    phase_min : degrees
+        Minimum phase in degrees of the N-circles
+    phase_max : degrees
+        Maximum phase in degrees of the N-circles
+
+    Return values
+    -------------
+    contours : complex array
+        Array of complex numbers corresponding to the contours.
+    """
+    # 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)
+
+# N-circle
+def n_circles(phases, mag_min=-40.0, mag_max=12.0):
+    """Constant-phase contours of the function H = G/(1+G).
+
+    Usage
+    ===== 
+    contour = n_circles(angles)
+
+    Parameters
+    ----------
+    phases : array-like
+        Array of phases in degrees of the N-circles
+    mag_min : dB
+        Minimum magnitude in dB of the N-circles
+    mag_max : dB
+        Maximum magnitude in dB of the N-circles
+
+    Return values
+    -------------
+    contours : complex array
+        Array of complex numbers corresponding to the contours.
+    """
+    # 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)
diff --git a/src/matlab.py b/src/matlab.py
diff --git a/src/rlocus.py b/src/rlocus.py
diff --git a/src/statesp.py b/src/statesp.py
diff --git a/src/xferfcn.py b/src/xferfcn.py

Original file line number	Diff line number	Diff line change
`@@ -37,7 +37,7 @@`
`37`	`37`	`# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF`
`38`	`38`	`# SUCH DAMAGE.`
`39`	`39`	`#`
`40`		`-# $Id: __init__.py 29 2010-11-06 13:03:55Z murrayrm $`
	`40`	`+# $Id$`
`41`	`41`
`42`	`42`	`"""Control System Library`
`43`	`43`
Original file line number	Diff line number	Diff line change
`@@ -38,7 +38,7 @@`
`38`	`38`	`# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF`
`39`	`39`	`# SUCH DAMAGE.`
`40`	`40`	`#`
`41`		`-# $Id: ctrlutil.py 21 2010-06-06 17:29:42Z murrayrm $`
	`41`	`+# $Id$`
`42`	`42`
`43`	`43`	`# Packages that we need access to`
`44`	`44`	`import scipy as sp`