|
From: <mur...@us...> - 2010-11-26 22:00:03
|
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.
|