python-control-discuss

[python-control-discuss] SF.net SVN: python-control:[33] trunk/src

[<mur...@us...>]

Revision: 33
          http://python-control.svn.sourceforge.net/python-control/?rev=33&view=rev
Author:   murrayrm
Date:     2010-11-26 21:59:57 +0000 (Fri, 26 Nov 2010)

Log Message:
-----------
* Added Nichols plot, contributed by Allan McInnes
* Allan also split out code to determine default frequency range

Modified Paths:
--------------
    trunk/src/freqplot.py
    trunk/src/matlab.py

Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py	2010-11-26 21:56:26 UTC (rev 32)
+++ trunk/src/freqplot.py	2010-11-26 21:59:57 UTC (rev 33)
@@ -4,7 +4,7 @@
 # Date: 24 May 09
 # 
 # This file contains some standard control system plots: Bode plots,
-# Nyquist plots and pole-zero diagrams
+# Nyquist plots, Nichols plots and pole-zero diagrams
 #
 # Copyright (c) 2010 by California Institute of Technology
 # All rights reserved.
@@ -46,6 +46,13 @@
 from ctrlutil import unwrap
 from bdalg import feedback
 
+#
+# Main plotting functions
+#
+# This section of the code contains the functions for generating
+# frequency domain plots
+#
+   
 # Bode plot
 def bode(syslist, omega=None, dB=False, Hz=False):
     """Bode plot for a system
@@ -81,36 +88,10 @@
     if (not getattr(syslist, '__iter__', False)):
         syslist = (syslist,)
 
-    #
     # Select a default range if none is provided
-    #
-    # This code looks at the poles and zeros of all of the systems that
-    # we are plotting and sets the frequency range to be one decade above
-    # and below the min and max feature frequencies, rounded to the nearest
-    # integer.  It excludes poles and zeros at the origin.  If no features
-    # are found, it turns logspace(-1, 1)
-    #
     if (omega == None):
-        # Find the list of all poles and zeros in the systems
-        features = np.array(())
-        for sys in syslist:
-            # Add new features to the list
-            features = np.concatenate((features, np.abs(sys.poles)))
-            features = np.concatenate((features, np.abs(sys.zeros)))
+        omega = default_frequency_range(syslist)
 
-        # Get rid of poles and zeros at the origin
-        features = features[features != 0];
-
-        # Make sure there is at least one point in the range
-        if (features.shape[0] == 0): features = [1];
-
-        # Take the log of the features
-        features = np.log10(features)
-
-        # Set the to be an order of magnitude beyond any features
-        omega = sp.logspace(np.floor(np.min(features))-1, 
-                            np.ceil(np.max(features))+1)
-
     for sys in syslist:
         # Get the magnitude and phase of the system
         mag, phase, omega = sys.freqresp(omega)
@@ -154,7 +135,7 @@
     return (211, 212)
 
 # Nyquist plot
-def nyquist(sys, omega=None):
+def nyquist(syslist, omega=None):
     """Nyquist plot for a system
 
     Usage
@@ -165,8 +146,8 @@
 
     Parameters
     ----------
-    sys : linsys
-        Linear input/output system
+    syslist : linsys
+        List of linear input/output systems (single system is OK)
     omega : freq_range
         Range of frequencies (list or bounds) in rad/sec
 
@@ -174,27 +155,82 @@
     -------------
     None
     """
-
+    # If argument was a singleton, turn it into a list
+    if (not getattr(syslist, '__iter__', False)):
+        syslist = (syslist,)
+        
     # Select a default range if none is provided
-    #! TODO: This needs to be made more intelligent
     if (omega == None):
-        omega = sp.logspace(-2, 2);
+        omega = default_frequency_range(syslist)
 
-    # Get the magnitude and phase of the system
-    mag, phase, omega = sys.freqresp(omega)
+    for sys in syslist:
+        # Get the magnitude and phase of the system
+        mag, phase, omega = sys.freqresp(omega)
+    
+        # Compute the primary curve
+        x = sp.multiply(mag, sp.cos(phase));
+        y = sp.multiply(mag, sp.sin(phase));
+    
+        # Plot the primary curve and mirror image
+        plt.plot(x, y, '-');
+        plt.plot(x, -y, '--');
 
-    # Compute the primary curve
-    x = sp.multiply(mag, sp.cos(phase));
-    y = sp.multiply(mag, sp.sin(phase));
+    # Mark the -1 point
+    plt.plot([-1], [0], 'r+')
 
-    # Plot the primary curve and mirror image
-    plt.plot(x, y, '-');
-    plt.plot(x, -y, '--');
+# Nichols plot
+# Contributed by Allan McInnes <All...@ca...>
+#! TODO: need unit test code
+def nichols(syslist, omega=None):
+    """Nichols plot for a system
 
-    # Mark the -1 point
-    plt.plot([-1], [0], '+k')
+    Usage
+    =====
+    magh = nichols(sys, omega=None)
 
+    Plots a Nichols plot for the system over a (optional) frequency range.
+
+    Parameters
+    ----------
+    syslist : linsys
+        List of linear input/output systems (single system is OK)
+    omega : freq_range
+        Range of frequencies (list or bounds) in rad/sec
+
+    Return values
+    -------------
+    None
+    """
+
+    # If argument was a singleton, turn it into a list
+    if (not getattr(syslist, '__iter__', False)):
+        syslist = (syslist,)
+
+    # Select a default range if none is provided
+    if (omega == None):
+        omega = default_frequency_range(syslist)
+
+    for sys in syslist:
+        # Get the magnitude and phase of the system
+        mag, phase, omega = sys.freqresp(omega)
+    
+        # 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')
+
+    # Mark the -180 point
+    plt.plot([-180], [0], 'r+')
+    
 # Gang of Four
+#! TODO: think about how (and whether) to handle lists of systems
 def gangof4(P, C, omega=None):
     """Plot the "Gang of 4" transfer functions for a system
 
@@ -220,7 +256,7 @@
     # Select a default range if none is provided
     #! TODO: This needs to be made more intelligent
     if (omega == None):
-        omega = sp.logspace(-2, 2);
+        omega = default_frequency_range((P,C))
 
     # Compute the senstivity functions
     L = P*C;
@@ -240,3 +276,60 @@
 
     mag, phase, omega = S.freqresp(omega);
     plt.subplot(224); plt.loglog(omega, mag);
+
+#
+# Utility functions
+#
+# This section of the code contains some utility functions for
+# generating frequency domain plots
+#
+   
+# Compute reasonable defaults for axes
+def default_frequency_range(syslist):
+    """Compute a reasonable default frequency range for frequency
+    domain plots.
+
+    Usage
+    =====
+    omega = default_frequency_range(syslist)
+
+    Finds a reasonable default frequency range by examining the features
+    (poles and zeros) of the systems in syslist.
+
+    Parameters
+    ----------
+    syslist : linsys
+        List of linear input/output systems (single system is OK)
+
+    Return values
+    -------------
+    omega : freq_range
+        Range of frequencies in rad/sec
+    """
+    # This code looks at the poles and zeros of all of the systems that
+    # we are plotting and sets the frequency range to be one decade above
+    # and below the min and max feature frequencies, rounded to the nearest
+    # integer.  It excludes poles and zeros at the origin.  If no features
+    # are found, it turns logspace(-1, 1)
+    
+    # Find the list of all poles and zeros in the systems
+    features = np.array(())
+    for sys in syslist:
+        # Add new features to the list
+        features = np.concatenate((features, np.abs(sys.poles)))
+        features = np.concatenate((features, np.abs(sys.zeros)))
+
+    # Get rid of poles and zeros at the origin
+    features = features[features != 0];
+
+    # Make sure there is at least one point in the range
+    if (features.shape[0] == 0): features = [1];
+
+    # Take the log of the features
+    features = np.log10(features)
+
+    # Set the range to be an order of magnitude beyond any features
+    omega = sp.logspace(np.floor(np.min(features))-1, 
+                        np.ceil(np.max(features))+1)   
+                        
+    return omega

Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py	2010-11-26 21:56:26 UTC (rev 32)
+++ trunk/src/matlab.py	2010-11-26 21:59:57 UTC (rev 33)
@@ -66,7 +66,7 @@
 
 # Import MATLAB-like functions that can be used as-is
 from ctrlutil import unwrap
-from freqplot import nyquist, gangof4
+from freqplot import nyquist, nichols, gangof4
 from bdalg import series, parallel, negate, feedback
 from pzmap import pzmap
 from statefbk import ctrb, obsv, place, lqr
@@ -155,7 +155,7 @@
    lti/bodemag    - Bode magnitude diagram only
    sigma          - singular value frequency plot
 *  nyquist        - Nyquist plot
-   nichols        - Nichols plot
+*  nichols        - Nichols plot
    margin         - gain and phase margins
    lti/allmargin  - all crossover frequencies and related gain/phase margins
    lti/freqresp   - frequency response over a frequency grid


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[python-control-discuss] SF.net SVN: python-control:[128] trunk/src

[<mur...@us...>]

Revision: 128
          http://python-control.svn.sourceforge.net/python-control/?rev=128&view=rev
Author:   murrayrm
Date:     2011-02-13 03:17:10 +0000 (Sun, 13 Feb 2011)

Log Message:
-----------
added new ngrid() command from Allan McInnes

Modified Paths:
--------------
    trunk/src/freqplot.py
    trunk/src/matlab.py

Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py	2011-02-10 18:37:09 UTC (rev 127)
+++ trunk/src/freqplot.py	2011-02-13 03:17:10 UTC (rev 128)
@@ -229,6 +229,74 @@
     # Mark the -180 point
     plt.plot([-180], [0], 'r+')
     
+# 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    
+    """
+    gmin_default = -40.0 # dB
+    gain_step = 20.0     # dB
+
+    if plt.gcf().gca().has_data():
+        pmin, pmax, gmin, gmax = plt.axis()
+    else:
+        pmin, pmax = -360.0, 0.0
+        gmin, gmax = gmin_default, 40.0
+
+    # Determine the bounds of the chart
+    mags_low_end = np.min([gmin_default, gain_step*np.floor(gmin/gain_step)])
+    phases_low_end = 360.0*np.ceil(pmin/360.0)
+    phases_high_end = 360.0*(1.0 + np.ceil(pmax/360.0))
+
+    # M-circle magnitudes - we adjust the lower end of the range to match
+    # any existing data
+    mags_fixed = np.array([-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])
+    mags_adjustable = np.arange(mags_low_end, np.min(mags_fixed), gain_step)
+    mags = np.concatenate((mags_adjustable, 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
+    mcs = m_circles(mags, pmin=np.min(phases), pmax=np.max(phases))
+    mag_m = 20*sp.log10(np.abs(mcs))
+    phase_m = sp.mod(sp.degrees(sp.angle(mcs)), -360.0) # Unwrap
+
+    # Find the N-contours
+    ncs = n_circles(phases, gmin=np.min(mags), gmax=np.max(mags))
+    mag_n = 20*sp.log10(np.abs(ncs))
+    phase_n = sp.mod(sp.degrees(sp.angle(ncs)), -360.0) # Unwrap
+
+    # Plot the contours
+    for phase_offset in np.arange(phases_low_end, phases_high_end, 360.0):
+        plt.plot(phase_m + phase_offset, mag_m, color='gray',
+                 linestyle='dashed', zorder=0)
+        plt.plot(phase_n + phase_offset, mag_n, color='gray',
+                 linestyle='dashed', zorder=0)
+
+    # Add magnitude labels
+    for x, y, m in zip(phase_m[:][-1], mag_m[:][-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([pmin, pmax, gmin, gmax])
+
 # Gang of Four
 #! TODO: think about how (and whether) to handle lists of systems
 def gangof4(P, C, omega=None):
@@ -333,3 +401,55 @@
                         np.ceil(np.max(features))+1)   
                         
     return omega
+
+# M-circle
+def m_circles(mags, pmin=-359.75, pmax=-0.25):
+    """Constant-magnitude contours of the function H = G/(1+G).
+
+    Usage
+    =====
+    contour = m_circle(mags)
+
+    Parameters
+    ----------
+    mags : array-like
+        Array of magnitudes in dB of the M-circles
+
+    Return values
+    -------------
+    contour : complex array
+        Array of complex numbers corresponding to the contour.
+    """
+    # Compute the contours in H-space
+    phases = sp.radians(sp.linspace(pmin, pmax, 500))
+    Hmag, Hphase = sp.meshgrid(10.0**(mags/20.0), phases)
+    H = Hmag*sp.exp(1.j*Hphase)
+    
+    # Invert H = G/(1+G) to get an expression for the contour in G-space
+    return H/(1.0 - H)
+
+# N-circle
+def n_circles(phases, gmin=-40.0, gmax=12.0):
+    """Constant-phase contours of the function H = G/(1+G).
+
+    Usage
+    ===== 
+    contour = n_circle(angles)
+
+    Parameters
+    ----------
+    phases : array-like
+        Array of phases in degrees of the N-circle
+
+    Return values
+    -------------
+    contour : complex array
+        Array of complex numbers corresponding to the contours.
+    """
+    # Compute the contours in H-space
+    mags = sp.linspace(10**(gmin/20.0), 10**(gmax/20.0), 2000)
+    Hphase, Hmag = sp.meshgrid(sp.radians(phases), mags)
+    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)

Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py	2011-02-10 18:37:09 UTC (rev 127)
+++ trunk/src/matlab.py	2011-02-13 03:17:10 UTC (rev 128)
@@ -371,6 +371,16 @@
     # Call the bode command
     return freqplot.bode(syslist, omega, **keywords)
 
+# Nichols chart grid
+def ngrid():
+    """Nichols chart grid.
+
+    Usage
+    =====
+    ngrid()
+    """
+    freqplot.nichols_grid()
+
 #
 # Modifications to scipy.signal functions
 #


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[python-control-discuss] SF.net SVN: python-control:[129] trunk/src

[<mur...@us...>]

Revision: 129
          http://python-control.svn.sourceforge.net/python-control/?rev=129&view=rev
Author:   murrayrm
Date:     2011-02-13 03:17:45 +0000 (Sun, 13 Feb 2011)

Log Message:
-----------
added svn:keyword ID

Modified Paths:
--------------
    trunk/src/delay.py

Property Changed:
----------------
    trunk/src/delay.py
    trunk/src/exception.py
    trunk/src/statefbk.py

Modified: trunk/src/delay.py
===================================================================
--- trunk/src/delay.py	2011-02-13 03:17:10 UTC (rev 128)
+++ trunk/src/delay.py	2011-02-13 03:17:45 UTC (rev 129)
@@ -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):
     """ 


Property changes on: trunk/src/delay.py
___________________________________________________________________
Added: svn:keywords
   + Id


Property changes on: trunk/src/exception.py
___________________________________________________________________
Added: svn:keywords
   + Id


Property changes on: trunk/src/statefbk.py
___________________________________________________________________
Added: svn:keywords
   + Id


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[python-control-discuss] SF.net SVN: python-control:[153] trunk/src

[<bli...@us...>]

Revision: 153
          http://python-control.svn.sourceforge.net/python-control/?rev=153&view=rev
Author:   blinkminster
Date:     2011-04-07 20:09:44 +0000 (Thu, 07 Apr 2011)

Log Message:
-----------
added margin command and associated matlab-like wrapper (does not plot yet). bugfixes to bode and default_frequency_range.

Modified Paths:
--------------
    trunk/src/freqplot.py
    trunk/src/matlab.py

Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py	2011-04-03 04:55:51 UTC (rev 152)
+++ trunk/src/freqplot.py	2011-04-07 20:09:44 UTC (rev 153)
@@ -55,12 +55,13 @@
 #
    
 # Bode plot
-def bode(syslist, omega=None, dB=False, Hz=False, color=None, Plot=True):
+def bode(syslist, omega=None, dB=False, Hz=False, deg=True, 
+        color=None, Plot=True):
     """Bode plot for a system
 
     Usage
     =====
-    (mag, phase, omega) = bode(syslist, omega=None, dB=False, Hz=False, color=None, Plot=True)
+    (mag, phase, omega) = bode(syslist, omega=None, dB=False, Hz=False, color=None, deg=True, Plot=True)
 
     Plots a Bode plot for the system over a (optional) frequency range.
 
@@ -74,23 +75,30 @@
         If True, plot result in dB
     Hz : boolean
         If True, plot frequency in Hz (omega must be provided in rad/sec)
+    color : matplotlib color
+        Color of line in bode plot
+    deg : boolean
+        If True, return phase in degrees (else radians)
     Plot : boolean
         If True, plot magnitude and phase
 
     Return values
     -------------
-    mag : magnitude array
-    phase : phase array
-    omega : frequency array
+    mag : magnitude array (list if len(syslist) > 1)
+    phase : phase array (list if len(syslist) > 1)
+    omega : frequency array (list if len(syslist) > 1)
 
     Notes
     -----
-    1. Alternatively, may use (mag, phase, freq) = sys.freqresp(freq) to generate the frequency response for a system.
+    1. Alternatively, you may use the lower-level method 
+    (mag, phase, freq) = sys.freqresp(freq) to generate the frequency 
+    response for a system, but it returns a MIMO response.
     """
     # If argument was a singleton, turn it into a list
     if (not getattr(syslist, '__iter__', False)):
         syslist = (syslist,)
 
+    mags, phases, omegas = [], [], []
     for sys in syslist:
         if (sys.inputs > 1 or sys.outputs > 1):
             #TODO: Add MIMO bode plots. 
@@ -104,10 +112,14 @@
             mag_tmp, phase_tmp, omega = sys.freqresp(omega)
             mag = np.squeeze(mag_tmp)
             phase = np.squeeze(phase_tmp)
+            phase = unwrap(phase)
             if Hz: omega = omega/(2*sp.pi)
             if dB: mag = 20*sp.log10(mag)
-            phase = unwrap(phase*180/sp.pi, 360)
-
+            if deg: phase = phase * 180 / sp.pi
+            
+            mags.append(mag)
+            phases.append(phase)
+            omegas.append(omega)
             # Get the dimensions of the current axis, which we will divide up
             #! TODO: Not current implemented; just use subplot for now
 
@@ -134,10 +146,14 @@
 
                 # Phase plot
                 plt.subplot(212);
+                if deg: 
+                    phase_deg = phase
+                else:
+                    phase_deg = phase * 180 / sp.pi
                 if color==None:
-                    plt.semilogx(omega, phase)
+                    plt.semilogx(omega, phase_deg)
                 else:
-                    plt.semilogx(omega, phase, color=color)
+                    plt.semilogx(omega, phase_deg, color=color)
                 plt.hold(True)
 
                 # Add a grid to the plot
@@ -151,7 +167,10 @@
                 else:
                     plt.xlabel("Frequency (rad/sec)")
 
-    return mag, phase, omega
+    if len(syslist) == 1:
+        return mags[0], phases[0], omegas[0]
+    else:
+        return mags, phases, omegas
 
 # Nyquist plot
 def nyquist(syslist, omega=None, Plot=True):
@@ -274,6 +293,101 @@
         phase = np.squeeze(phase_tmp)
         plt.subplot(224); plt.loglog(omega, mag);
 
+
+# gain and phase margins
+# contributed by Sawyer B. Fuller <mi...@ca...>
+def margin(sysdata, deg=True):
+    """Calculate gain and phase margins and associated crossover frequencies
+
+    Usage:
+    
+    gm, pm, sm, wg, wp, ws = margin(sysdata, deg=True)
+    
+    Parameters
+    ----------
+    sysdata: linsys or (mag, phase, omega) sequence 
+        sys : linsys
+            Linear SISO system
+        mag, phase, omega : sequence of array_like
+            Input magnitude, phase, and frequencies (rad/sec) sequence from 
+            bode frequency response data 
+    deg=True: boolean  
+        If true, all input and output phases in degrees, else in radians
+        
+    Returns
+    -------
+    gm, pm, sm, wg, wp, ws: float
+        Gain margin gm, phase margin pm, stability margin sm, and 
+        associated crossover
+        frequencies wg, wp, and ws of SISO open-loop. If more than
+        one crossover frequency is detected, returns the lowest corresponding
+        margin. 
+    """
+    #TODO do this precisely without the effects of discretization of frequencies?
+    #TODO assumes SISO
+    #TODO unit tests, margin plot
+
+    if (not getattr(sysdata, '__iter__', False)):
+        sys = sysdata
+        mag, phase, omega = bode(sys, deg=deg, Plot=False)
+    elif len(sysdata) == 3:
+        mag, phase, omega = sysdata
+    else: 
+        raise ValueError("Margin sysdata must be either a linear system or a 3-sequence of mag, phase, omega.")
+        
+    if deg:
+        cycle = 360. 
+        crossover = 180. 
+    else:
+        cycle = 2 * np.pi
+        crossover = np.pi
+        
+    wrapped_phase = -np.mod(phase, cycle)
+    
+    # phase margin from minimum phase among all gain crossovers
+    neg_mag_crossings_i = np.nonzero(np.diff(mag < 1) > 0)[0]
+    mag_crossings_p = wrapped_phase[neg_mag_crossings_i]
+    if len(neg_mag_crossings_i) == 0:
+        if mag[0] < 1: # gain always less than one
+            wp = np.nan
+            pm = np.inf
+        else: # gain always greater than one
+            print "margin: no magnitude crossings found"
+            wp = np.nan
+            pm = np.nan
+    else:
+        min_mag_crossing_i = neg_mag_crossings_i[np.argmin(mag_crossings_p)]
+        wp = omega[min_mag_crossing_i]
+        pm = crossover + phase[min_mag_crossing_i] 
+        if pm < 0:
+            print "warning: system unstable: negative phase margin"
+    
+    # gain margin from minimum gain margin among all phase crossovers
+    neg_phase_crossings_i = np.nonzero(np.diff(wrapped_phase < -crossover) > 0)[0]
+    neg_phase_crossings_g = mag[neg_phase_crossings_i]
+    if len(neg_phase_crossings_i) == 0:
+        wg = np.nan
+        gm = np.inf
+    else:
+        min_phase_crossing_i = neg_phase_crossings_i[
+            np.argmax(neg_phase_crossings_g)]
+        wg = omega[min_phase_crossing_i]
+        gm = abs(1/mag[min_phase_crossing_i])
+        if gm < 1: 
+            print "warning: system unstable: gain margin < 1"    
+
+    # stability margin from minimum abs distance from -1 point
+    if deg:
+        phase_rad = phase * np.pi / 180.
+    else:
+        phase_rad = phase
+    L = mag * np.exp(1j * phase_rad) # complex loop response to -1 pt
+    min_Lplus1_i = np.argmin(np.abs(L + 1))
+    sm = np.abs(L[min_Lplus1_i] + 1)
+    ws = phase[min_Lplus1_i]
+
+    return gm, pm, sm, wg, wp, ws 
+
 #
 # Utility functions
 #
@@ -311,6 +425,10 @@
     
     # Find the list of all poles and zeros in the systems
     features = np.array(())
+    
+    # detect if single sys passed by checking if it is sequence-like
+    if (not getattr(syslist, '__iter__', False)):
+        syslist = (syslist,)
     for sys in syslist:
         # Add new features to the list
         features = np.concatenate((features, np.abs(sys.pole())))

Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py	2011-04-03 04:55:51 UTC (rev 152)
+++ trunk/src/matlab.py	2011-04-07 20:09:44 UTC (rev 153)
@@ -781,6 +781,39 @@
     return rlist, klist
     
 
+def margin(*args):
+    """Calculate gain and phase margins and associated crossover frequencies
+
+    Usage:
+    gm, pm, wg, wp = margin(sys)
+    gm, pm, wg, wp = margin(mag,phase,w)
+    
+    Parameters
+    ----------
+    sys : linsys
+        Linear SISO system
+    mag, phase, w : array_like
+        Input magnitude, phase (in deg.), and frequencies (rad/sec) from bode
+        frequency response data
+
+    Returns
+    -------
+    gm, pm, wg, wp : float
+        Gain margin gm, phase margin pm (in deg), and associated crossover
+        frequencies wg and wp (in rad/sec) of SISO open-loop. If more than
+        one crossover frequency is detected, returns the lowest corresponding
+        margin. 
+    """
+    if len(args) == 1:
+        sys = args[0]
+        margins = freqplot.margin(sys)
+    elif len(args) == 3:
+        margins = freqplot.margin(args)
+    else: 
+        raise ValueError("Margin needs 1 or 3 arguments; received %i." 
+            % len(args))
+            
+    return margins[0], margins[1], margins[3], margins[4]
 #
 # Modifications to scipy.signal functions
 #


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[python-control-discuss] SF.net SVN: python-control:[180] trunk/src

[<mur...@us...>]

Revision: 180
          http://python-control.svn.sourceforge.net/python-control/?rev=180&view=rev
Author:   murrayrm
Date:     2012-01-08 02:57:46 +0000 (Sun, 08 Jan 2012)
Log Message:
-----------
small updates to documetnation + keywords property

Modified Paths:
--------------
    trunk/src/bdalg.py
    trunk/src/margins.py

Property Changed:
----------------
    trunk/src/margins.py

Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py	2012-01-08 02:56:24 UTC (rev 179)
+++ trunk/src/bdalg.py	2012-01-08 02:57:46 UTC (rev 180)
@@ -59,7 +59,7 @@
     """Return the series connection sys2 * sys1 for --> sys1 --> sys2 -->.
 
     Parameters
-    ---------
+    ----------
     sys1: scalar, StateSpace, or TransferFunction
     sys2: scalar, StateSpace, or TransferFunction
 
@@ -96,7 +96,7 @@
     Return the parallel connection sys1 + sys2.
 
     Parameters
-    ---------
+    ----------
     sys1: scalar, StateSpace, or TransferFunction
     sys2: scalar, StateSpace, or TransferFunction
 
@@ -121,7 +121,7 @@
     `sys1` is a scalar, then the output type is the type of `sys2`.
 
     Examples
-    -------
+    --------
     >>> sys3 = parallel(sys1, sys2) # Same as sys3 = sys1 + sys2.
 
     """

Modified: trunk/src/margins.py
===================================================================
--- trunk/src/margins.py	2012-01-08 02:56:24 UTC (rev 179)
+++ trunk/src/margins.py	2012-01-08 02:57:46 UTC (rev 180)
@@ -43,7 +43,7 @@
 Author: Richard M. Murray
 Date: 14 July 2011
 
-$Id: xferfcn.py 165 2011-06-26 02:44:09Z murrayrm $
+$Id$
 
 """
 
@@ -58,8 +58,8 @@
     """Calculate gain, phase and stability margins and associated
     crossover frequencies.
 
-    Usage:
-    
+    Usage
+    -----
     gm, pm, sm, wg, wp, ws = stability_margins(sysdata, deg=True)
     
     Parameters


Property changes on: trunk/src/margins.py
___________________________________________________________________
Added: svn:keywords
   + Id

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[python-control-discuss] SF.net SVN: python-control:[183] trunk/src

[<mur...@us...>]

Revision: 183
          http://python-control.svn.sourceforge.net/python-control/?rev=183&view=rev
Author:   murrayrm
Date:     2012-08-25 16:11:20 +0000 (Sat, 25 Aug 2012)
Log Message:
-----------
improvements to nyquist and bode plots

Modified Paths:
--------------
    trunk/src/__init__.py
    trunk/src/freqplot.py

Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py	2012-02-15 06:25:15 UTC (rev 182)
+++ trunk/src/__init__.py	2012-08-25 16:11:20 UTC (rev 183)
@@ -72,3 +72,11 @@
 from timeresp import forced_response, initial_response, step_response, \
     impulse_response
 from xferfcn import TransferFunction
+
+# Import some of the more common (and benign) MATLAB shortcuts
+from matlab import ss, tf, ss2tf, tf2ss, drss
+from matlab import pole, zero, evalfr, freqresp, dcgain
+from matlab import nichols, rlocus, margin
+        # bode and nyquist come directly from freqplot.py
+from matlab import step, impulse, initial, lsim
+

Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py	2012-02-15 06:25:15 UTC (rev 182)
+++ trunk/src/freqplot.py	2012-08-25 16:11:20 UTC (rev 183)
@@ -56,7 +56,7 @@
    
 # Bode plot
 def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True, 
-        color=None, Plot=True):
+        Plot=True, *args, **kwargs):
     """Bode plot for a system
 
     Plots a Bode plot for the system over a (optional) frequency range.
@@ -71,12 +71,12 @@
         If True, plot result in dB
     Hz : boolean
         If True, plot frequency in Hz (omega must be provided in rad/sec)
-    color : matplotlib color
-        Color of line in bode plot
     deg : boolean
         If True, return phase in degrees (else radians)
     Plot : boolean
         If True, plot magnitude and phase
+    *args, **kwargs: 
+        Additional options to matplotlib (color, linestyle, etc)
 
     Returns
     -------
@@ -95,7 +95,6 @@
 
     Examples
     --------
-    >>> from matlab import ss
     >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
     >>> mag, phase, omega = bode(sys)
     """
@@ -132,45 +131,28 @@
                 # Magnitude plot
                 plt.subplot(211); 
                 if dB:
-                    if color==None:
-                        plt.semilogx(omega, mag)
-                    else:
-                        plt.semilogx(omega, mag, color=color)
-                    plt.ylabel("Magnitude (dB)")
+                    plt.semilogx(omega, mag, *args, **kwargs)
                 else:
-                    if color==None:
-                        plt.loglog(omega, mag)
-                    else: 
-                        plt.loglog(omega, mag, color=color) 
-                    plt.ylabel("Magnitude")
+                    plt.loglog(omega, mag, *args, **kwargs)
+                plt.hold(True);
 
-                # Add a grid to the plot
+                # Add a grid to the plot + labeling
                 plt.grid(True)
                 plt.grid(True, which='minor')
-                plt.hold(True);
+                plt.ylabel("Magnitude (dB)" if dB else "Magnitude")
 
                 # Phase plot
                 plt.subplot(212);
-                if deg: 
-                    phase_deg = phase
-                else:
-                    phase_deg = phase * 180 / sp.pi
-                if color==None:
-                    plt.semilogx(omega, phase_deg)
-                else:
-                    plt.semilogx(omega, phase_deg, color=color)
-                plt.hold(True)
+                plt.semilogx(omega, phase, *args, **kwargs)
+                plt.hold(True);
 
-                # Add a grid to the plot
+                # Add a grid to the plot + labeling
                 plt.grid(True)
                 plt.grid(True, which='minor')
-                plt.ylabel("Phase (deg)")
+                plt.ylabel("Phase (deg)" if deg else "Phase (rad)")
 
                 # Label the frequency axis
-                if Hz:
-                    plt.xlabel("Frequency (Hz)")
-                else:
-                    plt.xlabel("Frequency (rad/sec)")
+                plt.xlabel("Frequency (Hz)" if Hz else "Frequency (rad/sec)")
 
     if len(syslist) == 1:
         return mags[0], phases[0], omegas[0]
@@ -178,7 +160,8 @@
         return mags, phases, omegas
 
 # Nyquist plot
-def nyquist_plot(syslist, omega=None, Plot=True):
+def nyquist_plot(syslist, omega=None, Plot=True, color='b', 
+                 labelFreq=0, *args, **kwargs):
     """Nyquist plot for a system
 
     Plots a Nyquist plot for the system over a (optional) frequency range.
@@ -190,7 +173,11 @@
     omega : freq_range
         Range of frequencies (list or bounds) in rad/sec
     Plot : boolean
-        if True, plot magnitude
+        If True, plot magnitude
+    labelFreq : int
+        Label every nth frequency on the plot
+    *args, **kwargs: 
+        Additional options to matplotlib (color, linestyle, etc)
 
     Returns
     -------
@@ -203,7 +190,6 @@
 
     Examples
     --------
-    >>> from matlab import ss
     >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
     >>> real, imag, freq = nyquist(sys)
     """
@@ -214,12 +200,14 @@
     # Select a default range if none is provided
     if (omega == None):
         omega = default_frequency_range(syslist)
+
     # Interpolate between wmin and wmax if a tuple or list are provided
     elif (isinstance(omega,list) | isinstance(omega,tuple)):
         # Only accept tuple or list of length 2
         if (len(omega) != 2):
             raise ValueError("Supported frequency arguments are (wmin,wmax) tuple or list, or frequency vector. ")
-        omega = np.logspace(np.log10(omega[0]),np.log10(omega[1]),num=50,endpoint=True,base=10.0)
+        omega = np.logspace(np.log10(omega[0]), np.log10(omega[1]),
+                            num=50, endpoint=True, base=10.0)
     for sys in syslist:
         if (sys.inputs > 1 or sys.outputs > 1):
             #TODO: Add MIMO nyquist plots. 
@@ -236,11 +224,33 @@
 
             if (Plot):    
                 # Plot the primary curve and mirror image
-                plt.plot(x, y, '-');
-                plt.plot(x, -y, '--');
+                plt.plot(x, y, '-', color=color, *args, **kwargs);
+                plt.plot(x, -y, '--', color=color, *args, **kwargs);
                 # Mark the -1 point
                 plt.plot([-1], [0], 'r+')
 
+            # Label the frequencies of the points
+            if (labelFreq):
+                for xpt, ypt, omegapt in zip(x, y, omega)[::labelFreq]:
+                    # Convert to Hz
+                    f = omegapt/(2*sp.pi)       
+
+                    # Factor out multiples of 1000 and limit the
+                    # result to the range [-8, 8].
+                    pow1000 = max(min(get_pow1000(f),8),-8) 
+
+                     # Get the SI prefix.
+                    prefix = gen_prefix(pow1000)
+                    
+                    # Apply the text. (Use a space before the text to
+                    # prevent overlap with the data.)
+                    #
+                    # np.round() is used because 0.99... appears
+                    # instead of 1.0, and this would otherwise be
+                    # truncated to 0.
+                    plt.text(xpt, ypt, 
+                             ' ' + str(int(np.round(f/1000**pow1000, 0))) + 
+                             ' ' + prefix + 'Hz') 
         return x, y, omega
 
 # Gang of Four
@@ -362,6 +372,49 @@
                         
     return omega
 
+#
+# KLD 5/23/11: Two functions to create nice looking labels
+#
+def get_pow1000(num):
+    '''Determine the exponent for which the significand of a number is within the
+    range [1, 1000).
+    '''
+    # Based on algorithm from http://www.mail-archive.com/mat...@li.../msg14433.html, accessed 2010/11/7
+    # by Jason Heeris 2009/11/18
+    from decimal import Decimal
+    from math import floor
+    dnum = Decimal(str(num))
+    if dnum == 0:
+        return 0
+    elif dnum < 0:
+        dnum = -dnum
+    return int(floor(dnum.log10()/3))
+
+def gen_prefix(pow1000):
+    '''Return the SI prefix for a power of 1000.
+    '''
+    # Prefixes according to Table 5 of [BIPM 2006] (excluding hecto,
+    # deca, deci, and centi).
+    if pow1000 < -8 or pow1000 > 8:
+        raise ValueError("Value is out of the range covered by the SI prefixes.")
+    return ['Y', # yotta (10^24)
+            'Z', # zetta (10^21)
+            'E', # exa (10^18)
+            'P', # peta (10^15)
+            'T', # tera (10^12)
+            'G', # giga (10^9)
+            'M', # mega (10^6)
+            'k', # kilo (10^3)
+            '', # (10^0)
+            'm', # milli (10^-3)
+            r'$\mu$', # micro (10^-6)
+            'n', # nano (10^-9)
+            'p', # pico (10^-12)
+            'f', # femto (10^-15)
+            'a', # atto (10^-18)
+            'z', # zepto (10^-21)
+            'y'][8 - pow1000] # yocto (10^-24)
+
 # Function aliases
 bode = bode_plot
 nyquist = nyquist_plot

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.