[python-control-discuss] Fwd: SF.net SVN: python-control:[21] trunk

[Richard M. <mu...@cd...>]

Tony and Erik:

Can you let me know if you got the note below from SourceForge.  I just posted some changes to the package (needed them for an example I'm putting together) and want to make sure that everyone is getting copies of the changes that got posted.  If not, I'll sort out what we need to change to make sure you get these sorts of e-mails.

Thanks.

-richard

Begin forwarded message:

> From: mur...@us...
> Date: 6 June 2010 10:29:43 PDT
> To: mu...@cd...
> Subject: SF.net SVN: python-control:[21] trunk
> 
> Revision: 21
>          http://python-control.svn.sourceforge.net/python-control/?rev=21&view=rev
> Author:   murrayrm
> Date:     2010-06-06 17:29:42 +0000 (Sun, 06 Jun 2010)
> 
> Log Message:
> -----------
> Added changes to allow bode() to plot frequency resp for multiple systems
>  * freqplot.bode now allows a list as the first argument
>  * Added matlab.bode wrapper for MATLAB syntax: bode(sys1, ..., sysN)
> 
> Added automatic frequency range determination in bode
>  * Computes range of frequencies based on pole/zero locations
>  * Eventually need to split this out as a ctrlutil function for use elsewhere
> 
> Modified Paths:
> --------------
>    trunk/ChangeLog
>    trunk/src/ctrlutil.py
>    trunk/src/freqplot.py
>    trunk/src/matlab.py
>    trunk/src/pzmap.py
>    trunk/src/statesp.py
>    trunk/src/xferfcn.py
> 
> Modified: trunk/ChangeLog
> ===================================================================
> --- trunk/ChangeLog	2010-06-01 02:08:20 UTC (rev 20)
> +++ trunk/ChangeLog	2010-06-06 17:29:42 UTC (rev 21)
> @@ -1,3 +1,25 @@
> +2010-06-06  Richard Murray  <murray@sumatra.local>
> +
> +	* src/matlab.py (bode): created a wrapper that allows MATLAB style
> +	arguments for bode (eg, bode(sys1, sys2))
> +
> +	* src/ctrlutil.py (issys): added new function to check if an object
> +	is a system (state space or transfer function).  Will generalize
> +	this latter to look for other compatible classes.
> +
> +	* src/freqplot.py (bode): Compute frequency range of bode plot based
> +	on poles and zeros
> +	(bode): Allow bode plot to be passed a list (or tuple) as the first
> +	argument, in which case multiple bode plots are generated
> +
> +	* src/statesp.py (StateSpace.zeros): new function to compute zeros
> +	for a state space system
> +	(StateSpace): defined new functions to compute poles of a state
> +	space system
> +
> +	* src/xferfcn.py (TransferFunction): defined new functions to
> +	compute poles and zeros of a transfer function.
> +
> 2010-05-31  Richard Murray  <murray@sumatra.local>
> 
> 	* src/exception.py (ControlNotImplemented): added new exception, to
> 
> Modified: trunk/src/ctrlutil.py
> ===================================================================
> --- trunk/src/ctrlutil.py	2010-06-01 02:08:20 UTC (rev 20)
> +++ trunk/src/ctrlutil.py	2010-06-06 17:29:42 UTC (rev 21)
> @@ -40,6 +40,12 @@
> # 
> # $Id$
> 
> +# Packages that we need access to
> +import scipy as sp
> +import statesp
> +import xferfcn
> +
> +# Specific functions that we use
> from scipy import pi
> 
> # Utility function to unwrap an angle measurement
> @@ -78,3 +84,12 @@
> 
>     # return the updated list
>     return angle
> +
> +# Determine if an object is a system
> +def issys(object):
> +    # Check for a member of one of the classes that we define here
> +    if (isinstance(object, (statesp.StateSpace, xferfcn.TransferFunction))): 
> +        return True
> +    
> +    # Didn't find anything that matched
> +    return False
> 
> Modified: trunk/src/freqplot.py
> ===================================================================
> --- trunk/src/freqplot.py	2010-06-01 02:08:20 UTC (rev 20)
> +++ trunk/src/freqplot.py	2010-06-06 17:29:42 UTC (rev 21)
> @@ -42,11 +42,12 @@
> 
> import matplotlib.pyplot as plt
> import scipy as sp
> +import numpy as np
> from ctrlutil import unwrap
> from bdalg import feedback
> 
> # Bode plot
> -def bode(sys, omega=None, dB=False, Hz=False):
> +def bode(syslist, omega=None, dB=False, Hz=False):
>     """Bode plot for a system
> 
>     Usage
> @@ -57,8 +58,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
>     dB : boolean
> @@ -76,45 +77,80 @@
>     1. Use (mag, phase, freq) = sys.freqresp(freq) to generate the 
>        frequency response for a system.
>     """
> +    # 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
> +    #
> +    # 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):
> -        omega = sp.logspace(-2, 2);
> +        # 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 the magnitude and phase of the system
> -    mag, phase, omega = sys.freqresp(omega)
> -    if Hz:
> -        omega = omega/(2*sp.pi)
> -    if dB:
> -        mag = 20*sp.log10(mag)
> -    phase = unwrap(phase*180/sp.pi, 360)
> +        # Get rid of poles and zeros at the origin
> +        features = features[features != 0];
> 
> -    # Get the dimensions of the current axis, which we will divide up
> -    #! TODO: Not current implemented; just use subplot for now
> +        # Make sure there is at least one point in the range
> +        if (features.shape[0] == 0): features = [1];
> 
> -    # Magnitude plot
> -    plt.subplot(211);
> -    if dB:
> -        plt.semilogx(omega, mag)
> -    else:
> -        plt.loglog(omega, mag)
> -    plt.grid(True)
> -    plt.grid(True, which='minor')
> -    if dB:
> -        plt.ylabel("Magnitude (dB)")
> -    else:
> -        plt.ylabel("Magnitude")
> +        # Take the log of the features
> +        features = np.log10(features)
> 
> -    # Phase plot
> -    plt.subplot(212);
> -    plt.semilogx(omega, phase)
> -    plt.grid(True)
> -    plt.grid(True, which='minor')
> -    plt.ylabel("Phase (deg)")
> -    if Hz:
> -        plt.xlabel("Frequency (Hz)")
> -    else:
> -        plt.xlabel("Frequency (rad/sec)")
> +        # 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)
> +        if Hz: omega = omega/(2*sp.pi)
> +        if dB: mag = 20*sp.log10(mag)
> +        phase = unwrap(phase*180/sp.pi, 360)
> +
> +        # Get the dimensions of the current axis, which we will divide up
> +        #! TODO: Not current implemented; just use subplot for now
> +
> +        # Magnitude plot
> +        plt.subplot(211); 
> +        if dB:
> +            plt.semilogx(omega, mag)
> +            plt.ylabel("Magnitude (dB)")
> +        else:
> +            plt.loglog(omega, mag)
> +            plt.ylabel("Magnitude")
> +
> +        # Add a grid to the plot
> +        plt.grid(True)
> +        plt.grid(True, which='minor')
> +        plt.hold(True);
> +
> +        # Phase plot
> +        plt.subplot(212);
> +        plt.semilogx(omega, phase)
> +        plt.hold(True)
> +
> +        # Add a grid to the plot
> +        plt.grid(True)
> +        plt.grid(True, which='minor')
> +        plt.ylabel("Phase (deg)")
> +
> +        # Label the frequency axis
> +        if Hz:
> +            plt.xlabel("Frequency (Hz)")
> +        else:
> +            plt.xlabel("Frequency (rad/sec)")
> +
>     return (211, 212)
> 
> # Nyquist plot
> 
> Modified: trunk/src/matlab.py
> ===================================================================
> --- trunk/src/matlab.py	2010-06-01 02:08:20 UTC (rev 20)
> +++ trunk/src/matlab.py	2010-06-06 17:29:42 UTC (rev 21)
> @@ -50,6 +50,7 @@
> 
> # Libraries that we make use of 
> import scipy as sp              # SciPy library (used all over)
> +import numpy as np              # NumPy library
> import scipy.signal as signal   # Signal processing library
> 
> # Import MATLAB-like functions that are defined in other packages
> @@ -57,11 +58,16 @@
> from scipy.signal import lsim, impulse, step
> from scipy import linspace, logspace
> 
> -# Import MATLAB-like functions that belong elsewhere in python-control
> -from ctrlutil import unwrap
> -from freqplot import bode, nyquist, gangof4
> +# Control system library
> +import ctrlutil
> +import freqplot
> from statesp import StateSpace
> from xferfcn import TransferFunction
> +from exception import *
> +
> +# Import MATLAB-like functions that can be used as-is
> +from ctrlutil import unwrap
> +from freqplot import nyquist, gangof4
> from bdalg import series, parallel, negate, feedback
> from pzmap import pzmap
> from statefbk import place, lqr
> @@ -304,3 +310,63 @@
> def freqresp(H, omega): 
>     """Return the frequency response for an object H at frequency omega"""
>     return H.freqresp(omega)
> +
> +# Bode plots
> +def bode(*args, **keywords):
> +    """Bode plot of the frequency response
> +
> +    Usage
> +    =====
> +    bode(sys)
> +    bode(sys, w)
> +    bode(sys1, sys2, ..., sysN)
> +    bode(sys1, sys2, ..., sysN, w)
> +    bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN')
> +    """
> +
> +    # If the first argument is a list, then assume python-control calling format
> +    if (getattr(args[0], '__iter__', False)):
> +        return freqplot.bode(*args, **keywords)
> +
> +    # Otherwise, run through the arguments and collect up arguments
> +    syslist = []; plotstyle=[]; omega=None;
> +    i = 0; 
> +    while i < len(args):
> +        # Check to see if this is a system of some sort
> +        if (ctrlutil.issys(args[i])): 
> +            # Append the system to our list of systems
> +            syslist.append(args[i])
> +            i += 1
> +
> +            # See if the next object is a plotsytle (string)
> +            if (i < len(args) and isinstance(args[i], str)):
> +                plotstyle.append(args[i])
> +                i += 1
> +
> +            # Go on to the next argument
> +            continue
> +
> +        # See if this is a frequency list
> +        elif (isinstance(args[i], (list, np.ndarray))):
> +            omega = args[i]
> +            i += 1
> +            break
> +
> +        else:
> +            raise ControlArgument("unrecognized argument type")
> +
> +    # Check to make sure that we processed all arguments
> +    if (i < len(args)):
> +        raise ControlArgument("not all arguments processed")
> +
> +    # Check to make sure we got the same number of plotstyles as systems
> +    if (len(plotstyle) != 0 and len(syslist) != len(plotstyle)):
> +        raise ControlArgument("number of systems and plotstyles should be equal")
> +
> +    # Warn about unimplemented plotstyles
> +    #! TODO: remove this when plot styles are implemented in bode()
> +    if (len(plotstyle) != 0):
> +        print("Warning (matabl.bode): plot styles not implemented");
> +
> +    # Call the bode command
> +    return freqplot.bode(syslist, omega, **keywords)
> 
> Modified: trunk/src/pzmap.py
> ===================================================================
> --- trunk/src/pzmap.py	2010-06-01 02:08:20 UTC (rev 20)
> +++ trunk/src/pzmap.py	2010-06-06 17:29:42 UTC (rev 21)
> @@ -45,6 +45,7 @@
> import xferfcn
> 
> # Compute poles and zeros for a system
> +#! TODO: extend this to handle state space systems
> def pzmap(sys, Plot=True):
>     """Plot a pole/zero map for a transfer function"""
>     if (isinstance(sys, xferfcn.TransferFunction)):
> 
> Modified: trunk/src/statesp.py
> ===================================================================
> --- trunk/src/statesp.py	2010-06-01 02:08:20 UTC (rev 20)
> +++ trunk/src/statesp.py	2010-06-06 17:29:42 UTC (rev 21)
> @@ -98,6 +98,18 @@
>         #! TODO: Not implemented
>         return None
> 
> +    # Compute poles and zeros
> +    def poles(self): return sp.roots(sp.poly(self.A))
> +    def zeros(self): 
> +        den = sp.poly1d(sp.poly(self.A))
> +
> +        # Compute the numerator based on zeros
> +        #! TODO: This is currently limited to SISO systems
> +        num = sp.poly1d(\
> +            sp.poly(self.A - sp.dot(self.B, self.C)) + (self.D[0] - 1) * den)
> +
> +        return (sp.roots(num))
> +
>     # Negation of a system
>     def __neg__(self):
>         """Negate a state space system"""
> 
> Modified: trunk/src/xferfcn.py
> ===================================================================
> --- trunk/src/xferfcn.py	2010-06-01 02:08:20 UTC (rev 20)
> +++ trunk/src/xferfcn.py	2010-06-06 17:29:42 UTC (rev 21)
> @@ -190,6 +190,10 @@
> 
>         return mag, phase, omega
> 
> +    # Compute poles and zeros
> +    def poles(self): return sp.roots(self.den)
> +    def zeros(self): return sp.roots(self.num)
> +
>     # Feedback around a transfer function
>     def feedback(sys1, sys2, sign=-1): 
>         """Feedback interconnection between two transfer functions"""
> 
> 
> This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.