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

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

OK, I've now set it up so that SVN commit messages go to the discuss list, so we should all see them.

-richard

On 6 Jun 2010, at 14:11 , Tony Vanelli wrote:

> Doesn't look like I received it. I also checked my spam folder just in  
> case.
> 
> -tv
> 
> 
> On Jun 6, 2010, at 10:31, Richard Murray <mu...@cd...> wrote:
> 
>> 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.
>> 
>> 
>> --- 
>> --- 
>> --- 
>> ---------------------------------------------------------------------
>> ThinkGeek and WIRED's GeekDad team up for the Ultimate
>> GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the
>> lucky parental unit.  See the prize list and enter to win:
>> http://p.sf.net/sfu/thinkgeek-promo
>> _______________________________________________
>> python-control-discuss mailing list
>> pyt...@li...
>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
> 
> ------------------------------------------------------------------------------
> ThinkGeek and WIRED's GeekDad team up for the Ultimate 
> GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the 
> lucky parental unit.  See the prize list and enter to win: 
> http://p.sf.net/sfu/thinkgeek-promo
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss