Menu
▾
▴
python-control-discuss — Discussion list for Python Control Systems Library
You can subscribe to this list here.
| 2010 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(19) |
Jul
|
Aug
|
Sep
|
Oct
|
Nov
(11) |
Dec
(5) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2011 |
Jan
|
Feb
(113) |
Mar
(12) |
Apr
(27) |
May
(2) |
Jun
(16) |
Jul
(6) |
Aug
(6) |
Sep
|
Oct
(3) |
Nov
(9) |
Dec
(2) |
| 2012 |
Jan
(5) |
Feb
(11) |
Mar
|
Apr
(3) |
May
|
Jun
|
Jul
|
Aug
(3) |
Sep
(7) |
Oct
(18) |
Nov
(18) |
Dec
|
| 2013 |
Jan
(4) |
Feb
(1) |
Mar
(3) |
Apr
(1) |
May
|
Jun
(33) |
Jul
(2) |
Aug
(5) |
Sep
|
Oct
|
Nov
|
Dec
(1) |
| 2014 |
Jan
(1) |
Feb
|
Mar
(8) |
Apr
|
May
(3) |
Jun
(3) |
Jul
(9) |
Aug
(5) |
Sep
(6) |
Oct
|
Nov
|
Dec
|
| 2015 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
(4) |
| 2017 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
(5) |
Nov
|
Dec
|
| 2018 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(1) |
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2019 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
(6) |
Sep
|
Oct
|
Nov
(2) |
Dec
|
| 2020 |
Jan
(1) |
Feb
(1) |
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2021 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
(4) |
Oct
|
Nov
|
Dec
|
|
From: <all...@us...> - 2011-02-19 03:15:24
|
Revision: 137
http://python-control.svn.sourceforge.net/python-control/?rev=137&view=rev
Author: allanmcinnes
Date: 2011-02-19 03:15:18 +0000 (Sat, 19 Feb 2011)
Log Message:
-----------
Bugfix freqplot.default_frequency_range(), which failed to deal properly with single systems. Also add an option to control the number of points. Set default to a larger number of points than the logspace default, so that we're more likely to get a smooth curve on bode, nyquist, and nichols plots.
Modified Paths:
--------------
trunk/src/freqplot.py
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2011-02-19 03:02:07 UTC (rev 136)
+++ trunk/src/freqplot.py 2011-02-19 03:15:18 UTC (rev 137)
@@ -192,10 +192,10 @@
Parameters
----------
- cl_mags : array-like (dB)
+ cl_mags : array-like (dB), optional
Array of closed-loop magnitudes defining a custom set of
M-circle iso-gain lines.
- cl_phases : array-like (degrees)
+ cl_phases : array-like (degrees), optional
Array of closed-loop phases defining a custom set of
N-circle iso-phase lines. Must be in the range -180.0 < cl_phases < 180.0
@@ -273,8 +273,8 @@
List of linear input/output systems (single system is OK)
omega : freq_range
Range of frequencies (list or bounds) in rad/sec
- grid : boolean
- True if the plot should include a Nichols-chart grid
+ grid : boolean, optional
+ True if the plot should include a Nichols-chart grid. Default is True.
Return values
-------------
@@ -326,10 +326,10 @@
Parameters
----------
- cl_mags : array-like (dB)
+ cl_mags : array-like (dB), optional
Array of closed-loop magnitudes defining the iso-gain lines on a
custom Nichols chart.
- cl_phases : array-like (degrees)
+ cl_phases : array-like (degrees), optional
Array of closed-loop phases defining the iso-phase lines on a custom
Nichols chart. Must be in the range -360 < cl_phases < 0
@@ -468,7 +468,7 @@
#
# Compute reasonable defaults for axes
-def default_frequency_range(syslist):
+def default_frequency_range(syslist, num=1000):
"""Compute a reasonable default frequency range for frequency
domain plots.
@@ -483,6 +483,8 @@
----------
syslist : linsys
List of linear input/output systems (single system is OK)
+ num : integer, optional
+ Number of samples to generate. Default is 50.
Return values
-------------
@@ -495,6 +497,10 @@
# integer. It excludes poles and zeros at the origin. If no features
# are found, it turns logspace(-1, 1)
+ # If argument was a singleton, turn it into a list
+ if (not getattr(syslist, '__iter__', False)):
+ syslist = (syslist,)
+
# Find the list of all poles and zeros in the systems
features = np.array(())
for sys in syslist:
@@ -513,7 +519,7 @@
# 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)
+ np.ceil(np.max(features))+1, num)
return omega
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <all...@us...> - 2011-02-19 03:02:13
|
Revision: 136
http://python-control.svn.sourceforge.net/python-control/?rev=136&view=rev
Author: allanmcinnes
Date: 2011-02-19 03:02:07 +0000 (Sat, 19 Feb 2011)
Log Message:
-----------
Add a Nyquist grid (aka Hall chart) showing M-circles and N-circles on the complex plane. Tweak Nichols grid to eliminate unnecessary use of kwargs.
Modified Paths:
--------------
trunk/src/freqplot.py
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2011-02-14 07:49:05 UTC (rev 135)
+++ trunk/src/freqplot.py 2011-02-19 03:02:07 UTC (rev 136)
@@ -178,6 +178,83 @@
# Mark the -1 point
plt.plot([-1], [0], 'r+')
+# Nyquist grid
+#! TODO: Consider making linestyle configurable
+def nyquist_grid(cl_mags=None, cl_phases=None):
+ """Nyquist plot grid of M-circles and N-circles (aka "Hall chart")
+
+ Usage
+ =====
+ nyquist_grid()
+
+ Plots a grid of M-circles and N-circles on the current axis, or
+ creates a default grid if no plot already exists.
+
+ Parameters
+ ----------
+ cl_mags : array-like (dB)
+ Array of closed-loop magnitudes defining a custom set of
+ M-circle iso-gain lines.
+ cl_phases : array-like (degrees)
+ Array of closed-loop phases defining a custom set of
+ N-circle iso-phase lines. Must be in the range -180.0 < cl_phases < 180.0
+
+ Return values
+ -------------
+ None
+ """
+ # Default chart size
+ re_min = -4.0
+ re_max = 3.0
+ im_min = -2.0
+ im_max = 2.0
+
+ # Find bounds of the current dataset, if there is one.
+ if plt.gcf().gca().has_data():
+ re_min, re_max, im_min, im_max = plt.axis()
+
+ # M-circle magnitudes.
+ if cl_mags is None:
+ cl_mags = np.array([-20.0, -10.0, -6.0, -4.0, -2.0, 0.0,
+ 2.0, 4.0, 6.0, 10.0, 20.0])
+
+ # N-circle phases (should be in the range -180.0 to 180.0)
+ if cl_phases is None:
+ cl_phases = np.array([-90.0, -60.0, -45.0, -30.0, -15.0,
+ 15.0, 30.0, 45.0, 60.0, 90.0])
+ else:
+ assert ((-180.0 < np.min(cl_phases)) and (np.max(cl_phases) < 180.0))
+
+ # Find the M-contours and N-contours
+ m = m_circles(cl_mags, phase_min=0.0, phase_max=359.99)
+ n = n_circles(cl_phases, mag_min=-40.0, mag_max=40.0)
+
+ # Draw contours
+ plt.plot(np.real(m), np.imag(m), color='gray', linestyle='dotted', zorder=0)
+ plt.plot(np.real(n), np.imag(n), color='gray', linestyle='dotted', zorder=0)
+
+ # Add magnitude labels
+ for i in range(0, len(cl_mags)):
+ if not cl_mags[i] == 0.0:
+ mag = 10.0**(cl_mags[i]/20.0)
+ x = -mag**2.0/(mag**2.0 - 1.0) # Center of the M-circle
+ y = np.abs(mag/(mag**2.0 - 1.0)) # Maximum point
+ else:
+ x, y = -0.5, im_max
+ plt.text(x, y, str(cl_mags[i]) + ' dB', size='small', color='gray')
+
+ # Add phase labels
+ for i in range(0, len(cl_phases)):
+ y = np.sign(cl_phases[i])*np.max(np.abs(np.imag(n)[:,i]))
+ p = str(cl_phases[i])
+ plt.text(-0.5, y, p + '$^\circ$', size='small', color='gray')
+
+ # Fit axes to original plot
+ plt.axis([re_min, re_max, im_min, im_max])
+
+# Make an alias
+hall_grid = nyquist_grid
+
# Nichols plot
# Contributed by Allan McInnes <All...@ca...>
#! TODO: need unit test code
@@ -209,7 +286,7 @@
syslist = (syslist,)
# Select a default range if none is provided
- if (omega == None):
+ if omega is None:
omega = default_frequency_range(syslist)
for sys in syslist:
@@ -236,7 +313,8 @@
nichols_grid()
# Nichols grid
-def nichols_grid(**kwargs):
+#! TODO: Consider making linestyle configurable
+def nichols_grid(cl_mags=None, cl_phases=None):
"""Nichols chart grid
Usage
@@ -270,10 +348,7 @@
ol_phase_min, ol_phase_max, ol_mag_min, ol_mag_max = plt.axis()
# M-circle magnitudes.
- if kwargs.has_key('cl_mags'):
- # Custom chart
- cl_mags = kwargs['cl_mags']
- else:
+ if cl_mags is None:
# Default chart magnitudes
# The key set of magnitudes are always generated, since this
# guarantees a recognizable Nichols chart grid.
@@ -289,11 +364,7 @@
cl_mags = np.concatenate((extended_cl_mags, key_cl_mags))
# N-circle phases (should be in the range -360 to 0)
- if kwargs.has_key('cl_phases'):
- # Custom chart
- cl_phases = kwargs['cl_phases']
- assert ((-360.0 < np.min(cl_phases)) and (np.max(cl_phases) < 0.0))
- else:
+ if cl_phases is None:
# Choose a reasonable set of default phases (denser if the open-loop
# data is restricted to a relatively small range of phases).
key_cl_phases = np.array([-0.25, -45.0, -90.0, -180.0, -270.0, -325.0, -359.75])
@@ -302,6 +373,8 @@
else:
other_cl_phases = np.arange(-10.0, -360.0, -20.0)
cl_phases = np.concatenate((key_cl_phases, other_cl_phases))
+ else:
+ assert ((-360.0 < np.min(cl_phases)) and (np.max(cl_phases) < 0.0))
# Find the M-contours
m = m_circles(cl_mags, phase_min=np.min(cl_phases), phase_max=np.max(cl_phases))
@@ -326,14 +399,14 @@
for phase_offset in phase_offsets:
# Draw M and N contours
plt.plot(m_phase + phase_offset, m_mag, color='gray',
- linestyle='dashed', zorder=0)
+ linestyle='dotted', zorder=0)
plt.plot(n_phase + phase_offset, n_mag, color='gray',
- linestyle='dashed', zorder=0)
+ linestyle='dotted', zorder=0)
# Add magnitude labels
for x, y, m in zip(m_phase[:][-1] + phase_offset, m_mag[:][-1], cl_mags):
align = 'right' if m < 0.0 else 'left'
- plt.text(x, y, str(m) + ' dB', size='small', ha=align)
+ plt.text(x, y, str(m) + ' dB', size='small', ha=align, color='gray')
# Fit axes to generated chart
plt.axis([phase_offset_min - 360.0, phase_offset_max - 360.0,
@@ -500,7 +573,7 @@
"""
# Convert magnitudes and phase range into a grid suitable for
# building contours
- phases = sp.radians(sp.linspace(phase_min, phase_max, 500))
+ phases = sp.radians(sp.linspace(phase_min, phase_max, 2000))
Gcl_mags, Gcl_phases = sp.meshgrid(10.0**(mags/20.0), phases)
return closed_loop_contours(Gcl_mags, Gcl_phases)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <all...@us...> - 2011-02-14 07:49:12
|
Revision: 135
http://python-control.svn.sourceforge.net/python-control/?rev=135&view=rev
Author: allanmcinnes
Date: 2011-02-14 07:49:05 +0000 (Mon, 14 Feb 2011)
Log Message:
-----------
Add the ability to define custom Nichols charts by passing magnitude and/or phase arrays to nichols_grid. Also some tweaks to the chart grid generation and axis configuration, and further clean-up of documentation and naming.
Modified Paths:
--------------
trunk/src/freqplot.py
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2011-02-13 16:49:58 UTC (rev 134)
+++ trunk/src/freqplot.py 2011-02-14 07:49:05 UTC (rev 135)
@@ -181,7 +181,7 @@
# Nichols plot
# Contributed by Allan McInnes <All...@ca...>
#! TODO: need unit test code
-def nichols(syslist, omega=None):
+def nichols(syslist, omega=None, grid=True):
"""Nichols plot for a system
Usage
@@ -196,6 +196,8 @@
List of linear input/output systems (single system is OK)
omega : freq_range
Range of frequencies (list or bounds) in rad/sec
+ grid : boolean
+ True if the plot should include a Nichols-chart grid
Return values
-------------
@@ -228,88 +230,114 @@
# Mark the -180 point
plt.plot([-180], [0], 'r+')
+
+ # Add grid
+ if grid:
+ nichols_grid()
# Nichols grid
-def nichols_grid():
+def nichols_grid(**kwargs):
"""Nichols chart grid
Usage
=====
nichols_grid()
- Plots a Nichols chart grid on the current axis.
+ Plots a Nichols chart grid on the current axis, or creates a new chart
+ if no plot already exists.
Parameters
----------
- None
+ cl_mags : array-like (dB)
+ Array of closed-loop magnitudes defining the iso-gain lines on a
+ custom Nichols chart.
+ cl_phases : array-like (degrees)
+ Array of closed-loop phases defining the iso-phase lines on a custom
+ Nichols chart. Must be in the range -360 < cl_phases < 0
Return values
-------------
None
"""
- mag_min_default = -40.0 # dB
- mag_step = 20.0 # dB
+ # Default chart size
+ ol_phase_min = -359.99
+ ol_phase_max = 0.0
+ ol_mag_min = -40.0
+ ol_mag_max = default_ol_mag_max = 50.0
+
+ # Find bounds of the current dataset, if there is one.
+ if plt.gcf().gca().has_data():
+ ol_phase_min, ol_phase_max, ol_mag_min, ol_mag_max = plt.axis()
- # 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))
+ if kwargs.has_key('cl_mags'):
+ # Custom chart
+ cl_mags = kwargs['cl_mags']
else:
- mags = mags_fixed
+ # Default chart magnitudes
+ # The key set of magnitudes are always generated, since this
+ # guarantees a recognizable Nichols chart grid.
+ key_cl_mags = 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])
+ # Extend the range of magnitudes if necessary. The extended arange
+ # will end up empty if no extension is required. Assumes that closed-loop
+ # magnitudes are approximately aligned with open-loop magnitudes beyond
+ # the value of np.min(key_cl_mags)
+ cl_mag_step = -20.0 # dB
+ extended_cl_mags = np.arange(np.min(key_cl_mags),
+ ol_mag_min + cl_mag_step, cl_mag_step)
+ cl_mags = np.concatenate((extended_cl_mags, key_cl_mags))
# 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])
+ if kwargs.has_key('cl_phases'):
+ # Custom chart
+ cl_phases = kwargs['cl_phases']
+ assert ((-360.0 < np.min(cl_phases)) and (np.max(cl_phases) < 0.0))
+ else:
+ # Choose a reasonable set of default phases (denser if the open-loop
+ # data is restricted to a relatively small range of phases).
+ key_cl_phases = np.array([-0.25, -45.0, -90.0, -180.0, -270.0, -325.0, -359.75])
+ if np.abs(ol_phase_max - ol_phase_min) < 90.0:
+ other_cl_phases = np.arange(-10.0, -360.0, -10.0)
+ else:
+ other_cl_phases = np.arange(-10.0, -360.0, -20.0)
+ cl_phases = np.concatenate((key_cl_phases, other_cl_phases))
# Find the M-contours
- m = m_circles(mags, phase_min=np.min(phases), phase_max=np.max(phases))
+ m = m_circles(cl_mags, phase_min=np.min(cl_phases), phase_max=np.max(cl_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 = n_circles(cl_phases, mag_min=np.min(cl_mags), mag_max=np.max(cl_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
+ # over the range -360 < phase < 0. 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
+ # for -360 < ol_phase_min < 0.
+ phase_offset_min = 360.0*np.ceil(ol_phase_min/360.0)
+ phase_offset_max = 360.0*np.ceil(ol_phase_max/360.0) + 360.0
phase_offsets = np.arange(phase_offset_min, phase_offset_max, 360.0)
+
for phase_offset in phase_offsets:
+ # Draw M and N contours
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)
+ # Add magnitude labels
+ for x, y, m in zip(m_phase[:][-1] + phase_offset, m_mag[:][-1], cl_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])
+ # Fit axes to generated chart
+ plt.axis([phase_offset_min - 360.0, phase_offset_max - 360.0,
+ np.min(cl_mags), np.max([ol_mag_max, default_ol_mag_max])])
# Gang of Four
#! TODO: think about how (and whether) to handle lists of systems
@@ -417,38 +445,44 @@
return omega
# Compute contours of a closed-loop transfer function
-def closed_loop_contours(Hmag, Hphase):
- """Contours of the function H = G/(1+G).
+def closed_loop_contours(Gcl_mags, Gcl_phases):
+ """Contours of the function Gcl = Gol/(1+Gol), where
+ Gol is an open-loop transfer function, and Gcl is a corresponding
+ closed-loop transfer function.
Usage
=====
- contours = closed_loop_contours(mags, phases)
+ contours = closed_loop_contours(Gcl_mags, Gcl_phases)
Parameters
----------
- mags : array-like
- Meshgrid array of magnitudes of the contours
- phases : array-like
- Meshgrid array of phases in radians of the contours
+ Gcl_mags : array-like
+ Array of magnitudes of the contours
+ Gcl_phases : array-like
+ 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)
+ # Compute the contours in Gcl-space. Since we're given closed-loop
+ # magnitudes and phases, this is just a case of converting them into
+ # a complex number.
+ Gcl = Gcl_mags*sp.exp(1.j*Gcl_phases)
- # Invert H = G/(1+G) to get an expression for the contours in G-space
- return H/(1.0 - H)
+ # Invert Gcl = Gol/(1+Gol) to map the contours into the open-loop space
+ return Gcl/(1.0 - Gcl)
# M-circle
def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
- """Constant-magnitude contours of the function H = G/(1+G).
+ """Constant-magnitude contours of the function Gcl = Gol/(1+Gol), where
+ Gol is an open-loop transfer function, and Gcl is a corresponding
+ closed-loop transfer function.
Usage
=====
- contours = m_circles(mags)
+ contours = m_circles(mags, phase_min, phase_max)
Parameters
----------
@@ -467,16 +501,18 @@
# 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)
+ Gcl_mags, Gcl_phases = sp.meshgrid(10.0**(mags/20.0), phases)
+ return closed_loop_contours(Gcl_mags, Gcl_phases)
# N-circle
def n_circles(phases, mag_min=-40.0, mag_max=12.0):
- """Constant-phase contours of the function H = G/(1+G).
+ """Constant-phase contours of the function Gcl = Gol/(1+Gol), where
+ Gol is an open-loop transfer function, and Gcl is a corresponding
+ closed-loop transfer function.
Usage
=====
- contour = n_circles(angles)
+ contours = n_circles(phases, mag_min, mag_max)
Parameters
----------
@@ -495,5 +531,5 @@
# 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)
+ Gcl_phases, Gcl_mags = sp.meshgrid(sp.radians(phases), mags)
+ return closed_loop_contours(Gcl_mags, Gcl_phases)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2011-02-13 16:50:05
|
Revision: 134
http://python-control.svn.sourceforge.net/python-control/?rev=134&view=rev
Author: murrayrm
Date: 2011-02-13 16:49:58 +0000 (Sun, 13 Feb 2011)
Log Message:
-----------
Copied over latest changes from v0.3d:
* 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
Modified Paths:
--------------
branches/control-0.4a/ChangeLog
branches/control-0.4a/setup.py
branches/control-0.4a/src/__init__.py
branches/control-0.4a/src/bdalg.py
branches/control-0.4a/src/ctrlutil.py
branches/control-0.4a/src/delay.py
branches/control-0.4a/src/freqplot.py
branches/control-0.4a/src/matlab.py
branches/control-0.4a/src/rlocus.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Added Paths:
-----------
branches/control-0.4a/examples/bdalg-matlab.py
branches/control-0.4a/examples/test-response.py
branches/control-0.4a/examples/test-statefbk.py
Property Changed:
----------------
branches/control-0.4a/src/__init__.py
branches/control-0.4a/src/bdalg.py
branches/control-0.4a/src/ctrlutil.py
branches/control-0.4a/src/delay.py
branches/control-0.4a/src/exception.py
branches/control-0.4a/src/freqplot.py
branches/control-0.4a/src/lti.py
branches/control-0.4a/src/matlab.py
branches/control-0.4a/src/modelsimp.py
branches/control-0.4a/src/pzmap.py
branches/control-0.4a/src/rlocus.py
branches/control-0.4a/src/robust.py
branches/control-0.4a/src/statefbk.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/test.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/ChangeLog
===================================================================
--- branches/control-0.4a/ChangeLog 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/ChangeLog 2011-02-13 16:49:58 UTC (rev 134)
@@ -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
Copied: branches/control-0.4a/examples/bdalg-matlab.py (from rev 132, trunk/examples/bdalg-matlab.py)
===================================================================
--- branches/control-0.4a/examples/bdalg-matlab.py (rev 0)
+++ branches/control-0.4a/examples/bdalg-matlab.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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;
Copied: branches/control-0.4a/examples/test-response.py (from rev 132, trunk/examples/test-response.py)
===================================================================
--- branches/control-0.4a/examples/test-response.py (rev 0)
+++ branches/control-0.4a/examples/test-response.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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)
Copied: branches/control-0.4a/examples/test-statefbk.py (from rev 132, trunk/examples/test-statefbk.py)
===================================================================
--- branches/control-0.4a/examples/test-statefbk.py (rev 0)
+++ branches/control-0.4a/examples/test-statefbk.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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
+
+
Modified: branches/control-0.4a/setup.py
===================================================================
--- branches/control-0.4a/setup.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/setup.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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='mu...@cd...',
Modified: branches/control-0.4a/src/__init__.py
===================================================================
--- branches/control-0.4a/src/__init__.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/__init__.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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
Property changes on: branches/control-0.4a/src/__init__.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: branches/control-0.4a/src/bdalg.py
===================================================================
--- branches/control-0.4a/src/bdalg.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/bdalg.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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$
"""
Property changes on: branches/control-0.4a/src/bdalg.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: branches/control-0.4a/src/ctrlutil.py
===================================================================
--- branches/control-0.4a/src/ctrlutil.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/ctrlutil.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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
Property changes on: branches/control-0.4a/src/ctrlutil.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: branches/control-0.4a/src/delay.py
===================================================================
--- branches/control-0.4a/src/delay.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/delay.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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: branches/control-0.4a/src/delay.py
___________________________________________________________________
Added: svn:keywords
+ Id
Property changes on: branches/control-0.4a/src/exception.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: branches/control-0.4a/src/freqplot.py
===================================================================
--- branches/control-0.4a/src/freqplot.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/freqplot.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -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 @@
plt.plot([-1], [0], 'r+')
return x, y, omega
-
# Nichols plot
# Contributed by Allan McInnes <All...@ca...>
#! 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 @@
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,35 +244,107 @@
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
def gangof4(P, C, omega=None):
@@ -395,3 +461,85 @@
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)
Property changes on: branches/control-0.4a/src/freqplot.py
___________________________________________________________________
Added: svn:keywords
+ Id
Property changes on: branches/control-0.4a/src/lti.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: branches/control-0.4a/src/matlab.py
===================================================================
--- branches/control-0.4a/src/matlab.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/matlab.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -51,7 +51,7 @@
Date: 29 May 09
Revised: Kevin K. Chen, Dec 10
-$Id: matlab.py 33 2010-11-26 21:59:57Z murrayrm $
+$Id$
"""
@@ -744,6 +744,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
#
Property changes on: branches/control-0.4a/src/matlab.py
___________________________________________________________________
Added: svn:keywords
+ Id
Property changes on: branches/control-0.4a/src/modelsimp.py
___________________________________________________________________
Added: svn:keywords
+ Id
Property changes on: branches/control-0.4a/src/pzmap.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: branches/control-0.4a/src/rlocus.py
===================================================================
--- branches/control-0.4a/src/rlocus.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/rlocus.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -40,7 +40,7 @@
# or a controls.TransferFunction object.
# * Added some comments to make sure I understand the code
#
-# $Id: rlocus.py 29 2010-11-06 13:03:55Z murrayrm $
+# $Id$
# Packages used by this module
from scipy import *
Property changes on: branches/control-0.4a/src/rlocus.py
___________________________________________________________________
Added: svn:keywords
+ Id
Property changes on: branches/control-0.4a/src/robust.py
___________________________________________________________________
Added: svn:keywords
+ Id
Property changes on: branches/control-0.4a/src/statefbk.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/statesp.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -68,7 +68,7 @@
Date: 24 May 09
Revised: Kevin K. Chen, Dec 10
-$Id: statepy 21 2010-06-06 17:29:42Z murrayrm $
+$Id$
"""
Property changes on: branches/control-0.4a/src/statesp.py
___________________________________________________________________
Added: svn:keywords
+ Id
Property changes on: branches/control-0.4a/src/test.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-13 13:21:29 UTC (rev 133)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-13 16:49:58 UTC (rev 134)
@@ -70,7 +70,7 @@
Date: 24 May 09
Revised: Kevin K. Chewn, Dec 10
-$Id: xferfcn.py 21 2010-06-06 17:29:42Z murrayrm $
+$Id$
"""
Property changes on: branches/control-0.4a/src/xferfcn.py
___________________________________________________________________
Added: svn:keywords
+ Id
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <all...@us...> - 2011-02-13 13:21:35
|
Revision: 133
http://python-control.svn.sourceforge.net/python-control/?rev=133&view=rev
Author: allanmcinnes
Date: 2011-02-13 13:21:29 +0000 (Sun, 13 Feb 2011)
Log Message:
-----------
Split out common code from M-circle and N-circle utility functions. Make variable naming through nichols_grid and utility functions more consistent, and provide more explanation of what's going on.
Modified Paths:
--------------
trunk/src/freqplot.py
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2011-02-13 03:35:38 UTC (rev 132)
+++ trunk/src/freqplot.py 2011-02-13 13:21:29 UTC (rev 133)
@@ -247,55 +247,69 @@
-------------
None
"""
- gmin_default = -40.0 # dB
- gain_step = 20.0 # dB
+ 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():
- pmin, pmax, gmin, gmax = plt.axis()
+ 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:
- 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))
+ 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
- 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
+ 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
- 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
+ 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
- for phase_offset in np.arange(phases_low_end, phases_high_end, 360.0):
- plt.plot(phase_m + phase_offset, mag_m, color='gray',
+ # 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(phase_n + phase_offset, mag_n, color='gray',
+ plt.plot(n_phase + phase_offset, n_mag, color='gray',
linestyle='dashed', zorder=0)
# Add magnitude labels
- for x, y, m in zip(phase_m[:][-1], mag_m[:][-1], mags):
+ 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([pmin, pmax, gmin, gmax])
+ plt.axis([phase_min, phase_max, mag_min, mag_max])
# Gang of Four
#! TODO: think about how (and whether) to handle lists of systems
@@ -402,54 +416,84 @@
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, pmin=-359.75, pmax=-0.25):
+def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
"""Constant-magnitude contours of the function H = G/(1+G).
Usage
=====
- contour = m_circle(mags)
+ 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
-------------
- contour : complex array
- Array of complex numbers corresponding to the contour.
+ contours : complex array
+ Array of complex numbers corresponding to the contours.
"""
- # Compute the contours in H-space
- phases = sp.radians(sp.linspace(pmin, pmax, 500))
+ # 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)
- 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)
+ return closed_loop_contours(Hmag, Hphase)
# N-circle
-def n_circles(phases, gmin=-40.0, gmax=12.0):
+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_circle(angles)
+ contour = n_circles(angles)
Parameters
----------
phases : array-like
- Array of phases in degrees of the N-circle
+ 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
-------------
- contour : complex array
+ contours : 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)
+ # 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)
- 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)
+ return closed_loop_contours(Hmag, Hphase)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2011-02-13 03:35:44
|
Revision: 132
http://python-control.svn.sourceforge.net/python-control/?rev=132&view=rev
Author: murrayrm
Date: 2011-02-13 03:35:38 +0000 (Sun, 13 Feb 2011)
Log Message:
-----------
tagging version 0.3d
Added Paths:
-----------
tags/control-0.3d/
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2011-02-13 03:32:19
|
Revision: 131
http://python-control.svn.sourceforge.net/python-control/?rev=131&view=rev
Author: murrayrm
Date: 2011-02-13 03:32:11 +0000 (Sun, 13 Feb 2011)
Log Message:
-----------
* Adding a few unit tests (will probably get redone in v0.4a)
* Updated list of changes
* Incremented version number to 0.3d
This is intended to be final set of changes before moving over to v0.4a
Modified Paths:
--------------
trunk/ChangeLog
trunk/setup.py
Added Paths:
-----------
trunk/examples/bdalg-matlab.py
trunk/examples/test-statefbk.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2011-02-13 03:25:52 UTC (rev 130)
+++ trunk/ChangeLog 2011-02-13 03:32:11 UTC (rev 131)
@@ -1,3 +1,16 @@
+2011-02-12 Richard Murray <murray@sumatra.local>
+
+ * setup.py: Updated version number to 0.3d, in preparation for release
+
+ * src/statefbk.py (lqr): Updated sb02md calling signature to match
+ latest slycot version (Lauren Padilla)
+
+ * src/freqplot.py (nichols_grid): new function from Allan McInnes
+ <all...@ca...> to generate a Nichols chart for a
+ given plot (equivalent to ngrid in MATLAB).
+
+ * src/matlab.py (ngrid): MATLAB compatible ngrid() command
+
2010-11-05 Richard Murray <murray@sumatra.local>
* external/yottalab.py: New file containing Roberto Bucher's control
Added: trunk/examples/bdalg-matlab.py
===================================================================
--- trunk/examples/bdalg-matlab.py (rev 0)
+++ trunk/examples/bdalg-matlab.py 2011-02-13 03:32:11 UTC (rev 131)
@@ -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;
Added: trunk/examples/test-statefbk.py
===================================================================
--- trunk/examples/test-statefbk.py (rev 0)
+++ trunk/examples/test-statefbk.py 2011-02-13 03:32:11 UTC (rev 131)
@@ -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
+
+
Modified: trunk/setup.py
===================================================================
--- trunk/setup.py 2011-02-13 03:25:52 UTC (rev 130)
+++ trunk/setup.py 2011-02-13 03:32:11 UTC (rev 131)
@@ -3,7 +3,7 @@
from distutils.core import setup
setup(name='control',
- version='0.3c',
+ version='0.3d',
description='Python Control Systems Library',
author='Richard Murray',
author_email='mu...@cd...',
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2011-02-13 03:25:58
|
Revision: 130
http://python-control.svn.sourceforge.net/python-control/?rev=130&view=rev
Author: murrayrm
Date: 2011-02-13 03:25:52 +0000 (Sun, 13 Feb 2011)
Log Message:
-----------
missing test file (was referred to in ChangeLog, but handn't been committed
Added Paths:
-----------
trunk/examples/test-response.py
Added: trunk/examples/test-response.py
===================================================================
--- trunk/examples/test-response.py (rev 0)
+++ trunk/examples/test-response.py 2011-02-13 03:25:52 UTC (rev 130)
@@ -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)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2011-02-13 03:17:51
|
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.
|
|
From: <mur...@us...> - 2011-02-13 03:17:16
|
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.
|
|
From: Lauren P. <lpa...@pr...> - 2011-02-10 18:47:28
|
Hi all, If you would like to try out the latest python-control version in branches/control-0.4a, at revision 127, you will need to replace some of your Slycot files with those of the same name in this attachment. This is just a temporary patch until our additions are officially incorporated into the Slycot library on Github. I will send an email as soon as that happens. To use the patch, just do the following: 1. Verify you have a working copy of Slycot and can re- build and install it. 2. Replace the following files with ours of the same name: in <your_slycot_dir>/slycot, replace: __init__.py, analysis.py, synthesis.py, transform.py in <your_slycot_dir>/slycot/src, replace: analysis.pyf, synthesis.pyf, transform.pyf 3. touch _wrapper.pyf 4. Build and install Slycot again Best wishes, Lauren |
|
From: <lpa...@us...> - 2011-02-10 18:37:19
|
Revision: 127
http://python-control.svn.sourceforge.net/python-control/?rev=127&view=rev
Author: lpadilla
Date: 2011-02-10 18:37:09 +0000 (Thu, 10 Feb 2011)
Log Message:
-----------
A few changes to unittests and usage of slycot routines.
Modified Paths:
--------------
branches/control-0.4a/src/TestSlycot.py
branches/control-0.4a/src/modelsimp.py
branches/control-0.4a/src/statefbk.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestSlycot.py
===================================================================
--- branches/control-0.4a/src/TestSlycot.py 2011-02-10 17:32:53 UTC (rev 126)
+++ branches/control-0.4a/src/TestSlycot.py 2011-02-10 18:37:09 UTC (rev 127)
@@ -8,19 +8,30 @@
class TestSlycot(unittest.TestCase):
+ """TestSlycot compares transfer function and state space conversions for
+ various numbers of inputs,outputs and states.
+ 1. Usually passes for SISO systems of any state dim, occasonally, there will be a dimension mismatch if the original randomly generated ss system is not minimal because td04ad returns a minimal system.
+
+ 2. For small systems with many inputs, n<<m, the tests fail because td04ad returns a minimal ss system which has fewer states than the original system. It is typical for systems with many more inputs than states to have extraneous states.
+
+ 3. For systems with larger dimensions, n~>5 and with 2 or more outputs the conversion to statespace (td04ad) intermittently results in an equivalent realization of higher order than the original tf order. We think this has to do with minimum realization tolerances in the Fortran. The algorithm doesn't recognize that two denominators are identical and so it creates a system with nearly duplicate eigenvalues and double the state dimension. This should not be a problem in the python-control usage because the common_den() method finds repeated roots within a tolerance that we specify.
+
+ Matlab: Matlab seems to force its statespace system output to have order less than or equal to the order of denominators provided, avoiding the problem of very large state dimension we describe in 3. It does however, still have similar problems with pole/zero cancellation such as we encounter in 2, where a statespace system may have fewer states than the original order of transfer function.
+ """
def setUp(self):
"""Define some test parameters."""
- self.numTests = 1
- self.maxStates = 10 #seems to fail rarely with 4, sometimes with 5, frequently with 6. Could it be a problem with the subsystems?
- self.maxIO = 10
-
+ self.numTests = 5
+ self.maxStates = 10
+ self.maxI = 1
+ self.maxO = 1
+
def testTF(self):
""" Directly tests the functions tb04ad and td04ad through direct comparison of transfer function coefficients.
Similar to TestConvert, but tests at a lower level.
"""
for states in range(1, self.maxStates):
- for inputs in range(1, self.maxIO+1):
- for outputs in range(1, self.maxIO+1):
+ for inputs in range(1, self.maxI+1):
+ for outputs in range(1, self.maxO+1):
for testNum in range(self.numTests):
ssOriginal = matlab.rss(states, inputs, outputs)
@@ -34,23 +45,25 @@
tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad('R',states,inputs,outputs,\
- ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=1e-10)
+ ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0)
ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\
- = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=1e-8)
+ = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=0.0)
tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\
tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad('R',\
ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\
- ssTransformed_D,tol1=1e-10)
- print 'size(Trans_A)=',ssTransformed_A.shape
- print 'Trans_nr=',ssTransformed_nr
- print 'tfOrig_index=',tfOriginal_index
- print 'tfOrig_ucoeff=',tfOriginal_ucoeff
- print 'tfOrig_dcoeff=',tfOriginal_dcoeff
- print 'tfTrans_index=',tfTransformed_index
- print 'tfTrans_ucoeff=',tfTransformed_ucoeff
- print 'tfTrans_dcoeff=',tfTransformed_dcoeff
+ ssTransformed_D,tol1=0.0)
+ #print 'size(Trans_A)=',ssTransformed_A.shape
+ print '===== Transformed SS =========='
+ print matlab.ss(ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D)
+ #print 'Trans_nr=',ssTransformed_nr
+ #print 'tfOrig_index=',tfOriginal_index
+ #print 'tfOrig_ucoeff=',tfOriginal_ucoeff
+ #print 'tfOrig_dcoeff=',tfOriginal_dcoeff
+ #print 'tfTrans_index=',tfTransformed_index
+ #print 'tfTrans_ucoeff=',tfTransformed_ucoeff
+ #print 'tfTrans_dcoeff=',tfTransformed_dcoeff
#Compare the TF directly, must match
#numerators
np.testing.assert_array_almost_equal(tfOriginal_ucoeff,tfTransformed_ucoeff,decimal=3)
@@ -61,24 +74,24 @@
"""Compare the bode reponses of the SS systems and TF systems to the original SS
They generally are different realizations but have same freq resp.
Currently this test may only be applied to SISO systems.
-
+ """
for states in range(1,self.maxStates):
for testNum in range(self.numTests):
- for inputs in range(1,self.maxIO+1):
- for outputs in range(1,self.maxIO+1):
+ for inputs in range(1,1):
+ for outputs in range(1,1):
ssOriginal = matlab.rss(states, inputs, outputs)
tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad('R',states,inputs,outputs,\
- ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=1e-10)
+ ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0)
ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\
- = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=1e-8)
+ = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=0.0)
tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\
tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad('R',\
ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\
- ssTransformed_D,tol1=1e-10)
+ ssTransformed_D,tol1=0.0)
numTransformed = np.array(tfTransformed_ucoeff)
denTransformed = np.array(tfTransformed_dcoeff)
@@ -88,21 +101,21 @@
ssTransformed = matlab.ss(ssTransformed_A,ssTransformed_B,ssTransformed_C,ssTransformed_D)
for inputNum in range(inputs):
for outputNum in range(outputs):
- #[ssOriginalMag,ssOriginalPhase,freq] = matlab.bode(ssOriginal,Plot=False)
+ [ssOriginalMag,ssOriginalPhase,freq] = matlab.bode(ssOriginal,Plot=False)
[tfOriginalMag,tfOriginalPhase,freq] = matlab.bode(matlab.tf(numOriginal[outputNum][inputNum],denOriginal[outputNum]),Plot=False)
- #[ssTransformedMag,ssTransformedPhase,freq] = matlab.bode(ssTransformed,freq,Plot=False)
+ [ssTransformedMag,ssTransformedPhase,freq] = matlab.bode(ssTransformed,freq,Plot=False)
[tfTransformedMag,tfTransformedPhase,freq] = matlab.bode(matlab.tf(numTransformed[outputNum][inputNum],denTransformed[outputNum]),freq,Plot=False)
- print 'numOrig=',numOriginal[outputNum][inputNum]
- print 'denOrig=',denOriginal[outputNum]
- print 'numTrans=',numTransformed[outputNum][inputNum]
- print 'denTrans=',denTransformed[outputNum]
- #np.testing.assert_array_almost_equal(ssOriginalMag,tfOriginalMag,decimal=3)
- #np.testing.assert_array_almost_equal(ssOriginalPhase,tfOriginalPhase,decimal=3)
- #np.testing.assert_array_almost_equal(ssOriginalMag,ssTransformedMag,decimal=3)
- #np.testing.assert_array_almost_equal(ssOriginalPhase,ssTransformedPhase,decimal=3)
- #np.testing.assert_array_almost_equal(tfOriginalMag,tfTransformedMag,decimal=3)
+ #print 'numOrig=',numOriginal[outputNum][inputNum]
+ #print 'denOrig=',denOriginal[outputNum]
+ #print 'numTrans=',numTransformed[outputNum][inputNum]
+ #print 'denTrans=',denTransformed[outputNum]
+ np.testing.assert_array_almost_equal(ssOriginalMag,tfOriginalMag,decimal=3)
+ np.testing.assert_array_almost_equal(ssOriginalPhase,tfOriginalPhase,decimal=3)
+ np.testing.assert_array_almost_equal(ssOriginalMag,ssTransformedMag,decimal=3)
+ np.testing.assert_array_almost_equal(ssOriginalPhase,ssTransformedPhase,decimal=3)
+ np.testing.assert_array_almost_equal(tfOriginalMag,tfTransformedMag,decimal=3)
np.testing.assert_array_almost_equal(tfOriginalPhase,tfTransformedPhase,decimal=2)
- """
+
#These are here for once the above is made into a unittest.
def suite():
return unittest.TestLoader().loadTestsFromTestCase(TestSlycot)
Modified: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py 2011-02-10 17:32:53 UTC (rev 126)
+++ branches/control-0.4a/src/modelsimp.py 2011-02-10 18:37:09 UTC (rev 127)
@@ -222,19 +222,10 @@
raise ControlSlycot("can't find slycot subroutine ab09ad")
job = 'B' # balanced (B) or not (N)
equil = 'N' # scale (S) or not (N)
- ordsel = 'F' # fixed truncation level (F) or find the truncation level given tol (A)
n = np.size(sys.A,0)
m = np.size(sys.B,1)
p = np.size(sys.C,0)
- nr = orders
- tol = 0.
- out = ab09ad(dico,job,equil,ordsel, n, m, p, nr, sys.A, sys.B, sys.C,tol)
- Ar = out[0][0:nr,0:nr]
- Br = out[1][0:nr,0:m]
- Cr = out[2][0:p,0:nr]
- hsv = out[3]
- iwarn = out[4]
- info = out[5]
+ Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=orders,tol=0.0)
rsys = StateSpace(Ar, Br, Cr, sys.D)
else:
Modified: branches/control-0.4a/src/statefbk.py
===================================================================
--- branches/control-0.4a/src/statefbk.py 2011-02-10 17:32:53 UTC (rev 126)
+++ branches/control-0.4a/src/statefbk.py 2011-02-10 18:37:09 UTC (rev 127)
@@ -177,7 +177,7 @@
sb02mt(nstates, ninputs, B, R, A, Q, N, jobl='N');
# Call the SLICOT function
- X,rcond,w,S,U = sb02md(nstates, A_b, G, Q_b, 'C')
+ X,rcond,w,S,U,A_inv = sb02md(nstates, A_b, G, Q_b, 'C')
# Now compute the return value
K = np.dot(np.linalg.inv(R), (np.dot(B.T, X) + N.T));
@@ -288,10 +288,10 @@
if e.real >= 0:
raise ValueError, "Oops, the system is unstable!"
if type=='c':
- trana = 'T'
+ tra = 'T'
C = -np.dot(sys.B,sys.B.transpose())
elif type=='o':
- trana = 'N'
+ tra = 'N'
C = -np.dot(sys.C.transpose(),sys.C)
else:
raise ValueError, "Oops, neither observable, nor controllable!"
@@ -305,7 +305,7 @@
n = sys.states
U = np.zeros((n,n))
A = sys.A
- out = sb03md(n, C, A, U, dico, 'X', 'N', trana)
- gram = out[0]
+ X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
+ gram = X
return gram
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-10 17:32:53 UTC (rev 126)
+++ branches/control-0.4a/src/statesp.py 2011-02-10 18:37:09 UTC (rev 127)
@@ -163,7 +163,6 @@
useless = []
# Search for useless states.
- tol = 1e-16
for i in range(self.states):
if (all(self.A[i, :] == zeros((1, self.states))) and
all(self.B[i, :] == zeros((1, self.inputs)))):
@@ -474,7 +473,7 @@
is still buggy! Advise converting state space sys back to tf to verify the transformation was correct."
#print num
#print shape(num)
- ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=1e-8)
+ ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=0.0)
states = ssout[0]
return StateSpace(ssout[1][:states, :states],
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-10 17:32:53 UTC (rev 126)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-10 18:37:09 UTC (rev 127)
@@ -212,7 +212,6 @@
"""
# Beware: this is a shallow copy. This should be okay.
- tol = 1e-16
data = [self.num, self.den]
for p in range(len(data)):
for i in range(self.outputs):
@@ -717,7 +716,7 @@
print "Warning: state space to transfer function conversion by tb04ad \
is still buggy!"
tfout = tb04ad('R',sys.states, sys.inputs, sys.outputs, sys.A, sys.B, sys.C,
- sys.D,tol1=1e-10)
+ sys.D,tol1=0.0)
# Preallocate outputs.
num = [[[] for j in range(sys.inputs)] for i in range(sys.outputs)]
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <lpa...@us...> - 2011-02-10 17:32:59
|
Revision: 126
http://python-control.svn.sourceforge.net/python-control/?rev=126&view=rev
Author: lpadilla
Date: 2011-02-10 17:32:53 +0000 (Thu, 10 Feb 2011)
Log Message:
-----------
Revised list of return values for sb02md to be compatible with Slycot commit c98ba6f30899ff2078e0.
Modified Paths:
--------------
trunk/src/statefbk.py
Modified: trunk/src/statefbk.py
===================================================================
--- trunk/src/statefbk.py 2011-02-09 16:08:17 UTC (rev 125)
+++ trunk/src/statefbk.py 2011-02-10 17:32:53 UTC (rev 126)
@@ -176,7 +176,7 @@
sb02mt(nstates, ninputs, B, R, A, Q, N, jobl='N');
# Call the SLICOT function
- X,rcond,w,S,U = sb02md(nstates, A_b, G, Q_b, 'C')
+ X,rcond,w,S,U,A_inv = sb02md(nstates, A_b, G, Q_b, 'C')
# Now compute the return value
K = np.dot(np.linalg.inv(R), (np.dot(B.T, X) + N.T));
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Richard M. <mu...@cd...> - 2011-02-10 15:07:47
|
Hi Riccardo. Thanks for the heads up. I'll try to have a look at this over the weekend. Python-control developers: if anyone has time to go in and do a quick fix on either trunk or v0.4a, go for it! -richard On 10 Feb 2011, at 6:26 , Riccardo Gori wrote: > Dear Professor Murray, > we are currently using python-control for our university work but we are experiencing problems with the last version of slycot: after this commit > > https://github.com/avventi/Slycot/commit/c98ba6f30899ff2078e0e8524657993c12cc4bb6 > > the API of sb02md changed (the number of outputs increased by one) and calling lqr may cause the "too many value to unpack" python error in statefbk.py file. > The fix should be straightforward for that line, but I don't know if in the rest of the python control codebase are there functions whom relies on the old API. > > Thanks for this great library, > Riccardo Gori |
|
From: <kk...@us...> - 2011-02-09 16:08:28
|
Revision: 125
http://python-control.svn.sourceforge.net/python-control/?rev=125&view=rev
Author: kkchen
Date: 2011-02-09 16:08:17 +0000 (Wed, 09 Feb 2011)
Log Message:
-----------
Removed build/* and some other stuff from version control.
Removed Paths:
-------------
branches/control-0.4a/MANIFEST.in
branches/control-0.4a/build/
Property Changed:
----------------
branches/control-0.4a/
branches/control-0.4a/doc/
Property changes on: branches/control-0.4a
___________________________________________________________________
Modified: svn:ignore
- build/
+ MANIFEST*
build
dist
Deleted: branches/control-0.4a/MANIFEST.in
===================================================================
--- branches/control-0.4a/MANIFEST.in 2011-02-08 23:04:46 UTC (rev 124)
+++ branches/control-0.4a/MANIFEST.in 2011-02-09 16:08:17 UTC (rev 125)
@@ -1,5 +0,0 @@
-include README
-include setup.py
-include src/*.py
-include examples/README examples/*.py
-prune examples/*-test.py
Property changes on: branches/control-0.4a/doc
___________________________________________________________________
Modified: svn:ignore
- _build/
_templates/
_static/
+ _build
_templates
_static
*.pdf
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Richard M. <mu...@cd...> - 2011-02-09 00:56:51
|
When I set up the subversion repository, I decided that it would be useful for developers to see commits as they came through. The Princeton group is sending through all of their commits from the past several months (so that we can keep track of this history), so there are ~100 message in a few hours. Sorry about that. You can all safely delete all of the messages for the individual commits. I will look through these and flag anything that I think people want to know about. Next time we do a bunch of commits in a row, I'll flip over the subversion e-mail alerts so that they don't go to everyone. (Kevin and crew: let me know if you have a lot more commits coming through and I'll do that just to spare everyone the many e-mails.) -richard On 8 Feb 2011, at 14:48 , Scott C. Livingston wrote: > Is the mailing list set to track commits, or is it just me? > > ~Scott > > On 8 February 2011 14:17, <kk...@us...> wrote: >> Revision: 99 >> http://python-control.svn.sourceforge.net/python-control/?rev=99&view=rev >> Author: kkchen >> Date: 2011-02-08 22:17:55 +0000 (Tue, 08 Feb 2011) >> >> Log Message: >> ----------- >> Implemented deep copy of A in gram. This is a temporary solution. >> >> Kevin K. Chen <kk...@pr...> >> >> Modified Paths: >> -------------- >> branches/control-0.4a/src/statefbk.py >> >> Modified: branches/control-0.4a/src/statefbk.py >> =================================================================== >> --- branches/control-0.4a/src/statefbk.py 2011-02-08 22:17:51 UTC (rev 98) >> +++ branches/control-0.4a/src/statefbk.py 2011-02-08 22:17:55 UTC (rev 99) >> @@ -266,6 +266,9 @@ >> >>> Wo = gram(sys,'o') >> >> """ >> + >> + from copy import deepcopy >> + >> #Check for ss system object, need a utility for this? >> >> #TODO: Check for continous or discrete, only continuous supported right now >> @@ -299,7 +302,8 @@ >> raise ControlSlycot("can't find slycot module 'sb03md'") >> n = sys.states >> U = np.zeros((n,n)) >> - out = sb03md(n, C, sys.A, U, dico, 'X', 'N', trana) >> + A = deepcopy(sys.A) >> + out = sb03md(n, C, A, U, dico, 'X', 'N', trana) >> gram = out[0] >> return gram >> >> >> >> This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. >> >> ------------------------------------------------------------------------------ >> The ultimate all-in-one performance toolkit: Intel(R) Parallel Studio XE: >> Pinpoint memory and threading errors before they happen. >> Find and fix more than 250 security defects in the development cycle. >> Locate bottlenecks in serial and parallel code that limit performance. >> http://p.sf.net/sfu/intel-dev2devfeb >> _______________________________________________ >> python-control-discuss mailing list >> pyt...@li... >> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >> >> > > > > -- > http://scottman.net > 865-964-5384 > > ------------------------------------------------------------------------------ > The ultimate all-in-one performance toolkit: Intel(R) Parallel Studio XE: > Pinpoint memory and threading errors before they happen. > Find and fix more than 250 security defects in the development cycle. > Locate bottlenecks in serial and parallel code that limit performance. > http://p.sf.net/sfu/intel-dev2devfeb > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: <kk...@us...> - 2011-02-08 23:04:52
|
Revision: 124
http://python-control.svn.sourceforge.net/python-control/?rev=124&view=rev
Author: kkchen
Date: 2011-02-08 23:04:46 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Removed .hgignore; set some svn:ignore properties.
Removed Paths:
-------------
.hgignore
Property Changed:
----------------
/
branches/control-0.4a/
branches/control-0.4a/build/
branches/control-0.4a/doc/
branches/control-0.4a/src/
Property changes on:
___________________________________________________________________
Added: svn:ignore
+ .DS_Store
Deleted: .hgignore
===================================================================
--- .hgignore 2011-02-08 22:19:58 UTC (rev 123)
+++ .hgignore 2011-02-08 23:04:46 UTC (rev 124)
@@ -1,17 +0,0 @@
-
-# Automatically generated by `hgimportsvn`
-syntax:glob
-.svn
-.hgsvn
-
-# These lines are suggested according to the svn:ignore property
-# Feel free to enable them by uncommenting them
-syntax:glob
-
-.DS_Store
-.*.swp
-build/*
-doc/_build/*
-doc/_templates/*
-doc/_static/*
-*.pyc
Property changes on: branches/control-0.4a
___________________________________________________________________
Added: svn:ignore
+ build/
Property changes on: branches/control-0.4a/build
___________________________________________________________________
Added: svn:ignore
+ *
Property changes on: branches/control-0.4a/doc
___________________________________________________________________
Added: svn:ignore
+ _build/
_templates/
_static/
Property changes on: branches/control-0.4a/src
___________________________________________________________________
Added: svn:ignore
+ *.pyc
.*.swp
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Scott C. L. <sli...@ca...> - 2011-02-08 22:49:36
|
Is the mailing list set to track commits, or is it just me? ~Scott On 8 February 2011 14:17, <kk...@us...> wrote: > Revision: 99 > http://python-control.svn.sourceforge.net/python-control/?rev=99&view=rev > Author: kkchen > Date: 2011-02-08 22:17:55 +0000 (Tue, 08 Feb 2011) > > Log Message: > ----------- > Implemented deep copy of A in gram. This is a temporary solution. > > Kevin K. Chen <kk...@pr...> > > Modified Paths: > -------------- > branches/control-0.4a/src/statefbk.py > > Modified: branches/control-0.4a/src/statefbk.py > =================================================================== > --- branches/control-0.4a/src/statefbk.py 2011-02-08 22:17:51 UTC (rev 98) > +++ branches/control-0.4a/src/statefbk.py 2011-02-08 22:17:55 UTC (rev 99) > @@ -266,6 +266,9 @@ > >>> Wo = gram(sys,'o') > > """ > + > + from copy import deepcopy > + > #Check for ss system object, need a utility for this? > > #TODO: Check for continous or discrete, only continuous supported right now > @@ -299,7 +302,8 @@ > raise ControlSlycot("can't find slycot module 'sb03md'") > n = sys.states > U = np.zeros((n,n)) > - out = sb03md(n, C, sys.A, U, dico, 'X', 'N', trana) > + A = deepcopy(sys.A) > + out = sb03md(n, C, A, U, dico, 'X', 'N', trana) > gram = out[0] > return gram > > > > This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. > > ------------------------------------------------------------------------------ > The ultimate all-in-one performance toolkit: Intel(R) Parallel Studio XE: > Pinpoint memory and threading errors before they happen. > Find and fix more than 250 security defects in the development cycle. > Locate bottlenecks in serial and parallel code that limit performance. > http://p.sf.net/sfu/intel-dev2devfeb > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > -- http://scottman.net 865-964-5384 |
|
From: <kk...@us...> - 2011-02-08 22:20:04
|
Revision: 123
http://python-control.svn.sourceforge.net/python-control/?rev=123&view=rev
Author: kkchen
Date: 2011-02-08 22:19:58 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Separated tests in TestSlycot.py
Lauren Padilla <lpa...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestSlycot.py
Modified: branches/control-0.4a/src/TestSlycot.py
===================================================================
--- branches/control-0.4a/src/TestSlycot.py 2011-02-08 22:19:52 UTC (rev 122)
+++ branches/control-0.4a/src/TestSlycot.py 2011-02-08 22:19:58 UTC (rev 123)
@@ -1,42 +1,112 @@
#!/usr/bin/env python
-#### THIS MUST BE MADE INTO A UNITTEST TO BE PART OF THE TESTING FUNCTIONS!!!!
-
import numpy as np
from slycot import tb04ad, td04ad
import matlab
+import unittest
-numTests = 1
-maxStates = 3
-maxIO = 3
-for states in range(1, maxStates):
- for inputs in range(1, maxIO):
- for outputs in range(1, maxIO):
- sys1 = matlab.rss(states, inputs, outputs)
- print "sys1"
- print sys1
+class TestSlycot(unittest.TestCase):
+ def setUp(self):
+ """Define some test parameters."""
+ self.numTests = 1
+ self.maxStates = 10 #seems to fail rarely with 4, sometimes with 5, frequently with 6. Could it be a problem with the subsystems?
+ self.maxIO = 10
+
+ def testTF(self):
+ """ Directly tests the functions tb04ad and td04ad through direct comparison of transfer function coefficients.
+ Similar to TestConvert, but tests at a lower level.
+ """
+ for states in range(1, self.maxStates):
+ for inputs in range(1, self.maxIO+1):
+ for outputs in range(1, self.maxIO+1):
+ for testNum in range(self.numTests):
+
+ ssOriginal = matlab.rss(states, inputs, outputs)
+
+ print '====== Original SS =========='
+ print ssOriginal
+ print 'states=',states
+ print 'inputs=',inputs
+ print 'outputs=',outputs
+
+
+ tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
+ tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad('R',states,inputs,outputs,\
+ ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=1e-10)
+
+ ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\
+ = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=1e-8)
+
+ tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\
+ tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad('R',\
+ ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\
+ ssTransformed_D,tol1=1e-10)
+ print 'size(Trans_A)=',ssTransformed_A.shape
+ print 'Trans_nr=',ssTransformed_nr
+ print 'tfOrig_index=',tfOriginal_index
+ print 'tfOrig_ucoeff=',tfOriginal_ucoeff
+ print 'tfOrig_dcoeff=',tfOriginal_dcoeff
+ print 'tfTrans_index=',tfTransformed_index
+ print 'tfTrans_ucoeff=',tfTransformed_ucoeff
+ print 'tfTrans_dcoeff=',tfTransformed_dcoeff
+ #Compare the TF directly, must match
+ #numerators
+ np.testing.assert_array_almost_equal(tfOriginal_ucoeff,tfTransformed_ucoeff,decimal=3)
+ #denominators
+ np.testing.assert_array_almost_equal(tfOriginal_dcoeff,tfTransformed_dcoeff,decimal=3)
+
+ def testFreqResp(self):
+ """Compare the bode reponses of the SS systems and TF systems to the original SS
+ They generally are different realizations but have same freq resp.
+ Currently this test may only be applied to SISO systems.
+
+ for states in range(1,self.maxStates):
+ for testNum in range(self.numTests):
+ for inputs in range(1,self.maxIO+1):
+ for outputs in range(1,self.maxIO+1):
+ ssOriginal = matlab.rss(states, inputs, outputs)
+
+ tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
+ tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad('R',states,inputs,outputs,\
+ ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=1e-10)
+
+ ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\
+ = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=1e-8)
+
+ tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\
+ tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad('R',\
+ ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\
+ ssTransformed_D,tol1=1e-10)
- sys2 = tb04ad(states,inputs,outputs,sys1.A,sys1.B,sys1.C,sys1.D,outputs,outputs,inputs)
- print "sys2"
- print sys2
-
- ldwork = 100*max(1,states+max(states,max(3*inputs,3*outputs)))
- dwork = np.zeros(ldwork)
- sys3 = td04ad(inputs,outputs,sys2[4],sys2[5],sys2[6])
- #sys3 = td04ad(inputs,outputs,sys2[4],sys2[5],sys2[6],ldwork)
- print "sys3"
- print sys3
-
- sys4 = tb04ad(states,inputs,outputs,sys3[1][0:states,0:states],sys3[2][0:states,0:inputs],sys3[3][0:outputs,0:states],sys3[4],outputs,outputs,inputs)
- print "sys4"
- print sys4
-
+ numTransformed = np.array(tfTransformed_ucoeff)
+ denTransformed = np.array(tfTransformed_dcoeff)
+ numOriginal = np.array(tfOriginal_ucoeff)
+ denOriginal = np.array(tfOriginal_dcoeff)
+
+ ssTransformed = matlab.ss(ssTransformed_A,ssTransformed_B,ssTransformed_C,ssTransformed_D)
+ for inputNum in range(inputs):
+ for outputNum in range(outputs):
+ #[ssOriginalMag,ssOriginalPhase,freq] = matlab.bode(ssOriginal,Plot=False)
+ [tfOriginalMag,tfOriginalPhase,freq] = matlab.bode(matlab.tf(numOriginal[outputNum][inputNum],denOriginal[outputNum]),Plot=False)
+ #[ssTransformedMag,ssTransformedPhase,freq] = matlab.bode(ssTransformed,freq,Plot=False)
+ [tfTransformedMag,tfTransformedPhase,freq] = matlab.bode(matlab.tf(numTransformed[outputNum][inputNum],denTransformed[outputNum]),freq,Plot=False)
+ print 'numOrig=',numOriginal[outputNum][inputNum]
+ print 'denOrig=',denOriginal[outputNum]
+ print 'numTrans=',numTransformed[outputNum][inputNum]
+ print 'denTrans=',denTransformed[outputNum]
+ #np.testing.assert_array_almost_equal(ssOriginalMag,tfOriginalMag,decimal=3)
+ #np.testing.assert_array_almost_equal(ssOriginalPhase,tfOriginalPhase,decimal=3)
+ #np.testing.assert_array_almost_equal(ssOriginalMag,ssTransformedMag,decimal=3)
+ #np.testing.assert_array_almost_equal(ssOriginalPhase,ssTransformedPhase,decimal=3)
+ #np.testing.assert_array_almost_equal(tfOriginalMag,tfTransformedMag,decimal=3)
+ np.testing.assert_array_almost_equal(tfOriginalPhase,tfTransformedPhase,decimal=2)
+ """
#These are here for once the above is made into a unittest.
def suite():
- return unittest.TestLoader().loadTestsFromTestCase(TestSlycot)
+ return unittest.TestLoader().loadTestsFromTestCase(TestSlycot)
if __name__=='__main__':
- unittest.main()
+ unittest.main()
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:19:58
|
Revision: 122
http://python-control.svn.sourceforge.net/python-control/?rev=122&view=rev
Author: kkchen
Date: 2011-02-08 22:19:52 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Modified tb04ad and td04ad function calls to match changes to the slycot wrappers.
Lauren Padilla <lpa...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestConvert.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestConvert.py
===================================================================
--- branches/control-0.4a/src/TestConvert.py 2011-02-08 22:19:48 UTC (rev 121)
+++ branches/control-0.4a/src/TestConvert.py 2011-02-08 22:19:52 UTC (rev 122)
@@ -25,13 +25,13 @@
"""Set up testing parameters."""
# Number of times to run each of the randomized tests.
- self.numTests = 10 #almost guarantees failure
+ self.numTests = 1 #almost guarantees failure
# Maximum number of states to test + 1
self.maxStates = 20
# Maximum number of inputs and outputs to test + 1
self.maxIO = 20
# Set to True to print systems to the output.
- self.debug = False
+ self.debug = True
def printSys(self, sys, ind):
"""Print system to the standard output."""
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:19:48 UTC (rev 121)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:19:52 UTC (rev 122)
@@ -163,6 +163,7 @@
useless = []
# Search for useless states.
+ tol = 1e-16
for i in range(self.states):
if (all(self.A[i, :] == zeros((1, self.states))) and
all(self.B[i, :] == zeros((1, self.inputs)))):
@@ -473,7 +474,7 @@
is still buggy! Advise converting state space sys back to tf to verify the transformation was correct."
#print num
#print shape(num)
- ssout = td04ad(sys.inputs, sys.outputs, index, den, num)
+ ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=1e-8)
states = ssout[0]
return StateSpace(ssout[1][:states, :states],
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:19:48 UTC (rev 121)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:19:52 UTC (rev 122)
@@ -212,6 +212,7 @@
"""
# Beware: this is a shallow copy. This should be okay.
+ tol = 1e-16
data = [self.num, self.den]
for p in range(len(data)):
for i in range(self.outputs):
@@ -219,7 +220,7 @@
# Find the first nontrivial coefficient.
nonzero = None
for k in range(data[p][i][j].size):
- if data[p][i][j][k]:
+ if (data[p][i][j][k]):
nonzero = k
break
@@ -594,6 +595,8 @@
# Modify the numerators so that they each take the common denominator.
num = deepcopy(self.num)
+ if isinstance(den,float):
+ den = array([den])
for i in range(self.outputs):
for j in range(self.inputs):
# The common denominator has leading coefficient 1. Scale out
@@ -713,8 +716,8 @@
# buggy!
print "Warning: state space to transfer function conversion by tb04ad \
is still buggy!"
- tfout = tb04ad(sys.states, sys.inputs, sys.outputs, sys.A, sys.B, sys.C,
- sys.D, sys.outputs, sys.outputs, sys.inputs)
+ tfout = tb04ad('R',sys.states, sys.inputs, sys.outputs, sys.A, sys.B, sys.C,
+ sys.D,tol1=1e-10)
# Preallocate outputs.
num = [[[] for j in range(sys.inputs)] for i in range(sys.outputs)]
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:19:54
|
Revision: 121
http://python-control.svn.sourceforge.net/python-control/?rev=121&view=rev
Author: kkchen
Date: 2011-02-08 22:19:48 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added a keyword for plotting and now returning mag/phase/freq etc. for nyquist and nichols.
Lauren Padilla <lpa...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/freqplot.py
Modified: branches/control-0.4a/src/freqplot.py
===================================================================
--- branches/control-0.4a/src/freqplot.py 2011-02-08 22:19:44 UTC (rev 120)
+++ branches/control-0.4a/src/freqplot.py 2011-02-08 22:19:48 UTC (rev 121)
@@ -153,12 +153,12 @@
return mag, phase, omega
# Nyquist plot
-def nyquist(syslist, omega=None):
+def nyquist(syslist, omega=None, Plot=True):
"""Nyquist plot for a system
Usage
=====
- magh = nyquist(sys, omega=None)
+ real, imag, freq = nyquist(sys, omega=None, Plot=True)
Plots a Nyquist plot for the system over a (optional) frequency range.
@@ -168,10 +168,14 @@
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 magnitude
Return values
-------------
- None
+ real : real part of the frequency response array
+ imag : imaginary part of the frequency response array
+ freq : frequencies
"""
# If argument was a singleton, turn it into a list
if (not getattr(syslist, '__iter__', False)):
@@ -199,22 +203,25 @@
# 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, '--');
- # Mark the -1 point
- plt.plot([-1], [0], 'r+')
+ if (Plot):
+ # Plot the primary curve and mirror image
+ plt.plot(x, y, '-');
+ plt.plot(x, -y, '--');
+ # Mark the -1 point
+ plt.plot([-1], [0], 'r+')
+
+ return x, y, omega
+
# Nichols plot
# Contributed by Allan McInnes <All...@ca...>
#! TODO: need unit test code
-def nichols(syslist, omega=None):
+def nichols(syslist, omega=None, Plot=True):
"""Nichols plot for a system
Usage
=====
- magh = nichols(sys, omega=None)
+ mag, phase, freq = nichols(sys, omega=None, Plot=True)
Plots a Nichols plot for the system over a (optional) frequency range.
@@ -224,10 +231,15 @@
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
-------------
- None
+ mag : array of magnitudes
+ phase : array of phases
+ freq : array of frequencies
+
"""
# If argument was a singleton, turn it into a list
@@ -254,16 +266,19 @@
x = unwrap(sp.degrees(phase), 360)
y = 20*sp.log10(mag)
- # Generate the plot
- plt.plot(x, y)
+ if (Plot):
+ # 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
+
# Gang of Four
#! TODO: think about how (and whether) to handle lists of systems
def gangof4(P, C, omega=None):
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:19:50
|
Revision: 120
http://python-control.svn.sourceforge.net/python-control/?rev=120&view=rev
Author: kkchen
Date: 2011-02-08 22:19:44 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Created copy constructors for StateSpace and TransferFunction.
The copy() methods have been removed, since copy constructors should be used.
Also, calling ss on a StateSpace object and tf on a TransferFunction object now
return deep copies.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/matlab.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/matlab.py
===================================================================
--- branches/control-0.4a/src/matlab.py 2011-02-08 22:19:39 UTC (rev 119)
+++ branches/control-0.4a/src/matlab.py 2011-02-08 22:19:44 UTC (rev 120)
@@ -59,6 +59,7 @@
import scipy as sp # SciPy library (used all over)
import numpy as np # NumPy library
import scipy.signal as signal # Signal processing library
+from copy import deepcopy
# Import MATLAB-like functions that are defined in other packages
from scipy.signal import zpk2ss, ss2zpk, tf2zpk, zpk2tf
@@ -305,7 +306,7 @@
elif len(args) == 1:
sys = args[0]
if isinstance(sys, StateSpace):
- return sys
+ return deepcopy(sys)
elif isinstance(sys, TransferFunction):
return tf2ss(sys)
else:
@@ -366,7 +367,7 @@
if isinstance(sys, StateSpace):
return ss2tf(sys)
elif isinstance(sys, TransferFunction):
- return sys
+ return deepcopy(sys)
else:
raise TypeError("tf(sys): sys must be a StateSpace or \
TransferFunction object. It is %s." % type(sys))
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:19:39 UTC (rev 119)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:19:44 UTC (rev 120)
@@ -78,7 +78,6 @@
from numpy.linalg import inv, det, solve
from numpy.linalg.linalg import LinAlgError
from scipy.signal import lti
-from copy import deepcopy
from slycot import td04ad
from lti import Lti
import xferfcn
@@ -96,9 +95,30 @@
"""
- def __init__(self, A=0, B=0, C=0, D=1):
- """Construct a state space object. The default is unit static gain."""
+ def __init__(self, *args):
+ """Construct a state space object.
+ The default constructor is StateSpace(A, B, C, D), where A, B, C, D are
+ matrices or equivalent objects. To call the copy constructor, call
+ StateSpace(sys), where sys is a StateSpace object.
+
+ """
+
+ if len(args) == 4:
+ # The user provided A, B, C, and D matrices.
+ (A, B, C, D) = args
+ elif len(args) == 1:
+ # Use the copy constructor.
+ if not isinstance(args[0], StateSpace):
+ raise TypeError("The one-argument constructor can only take in \
+a StateSpace object. Recived %s." % type(args[0]))
+ A = args[0].A
+ B = args[0].B
+ C = args[0].C
+ D = args[0].D
+ else:
+ raise ValueError("Needs 1 or 4 arguments; received %i." % len(args))
+
# Here we're going to convert inputs to matrices, if the user gave a
# non-matrix type.
matrices = [A, B, C, D]
@@ -171,11 +191,6 @@
self.inputs = self.B.shape[1]
self.outputs = self.C.shape[0]
- def copy(self):
- """Return a deep copy of the instance."""
-
- return deepcopy(self)
-
def __str__(self):
"""String representation of the state space."""
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:19:39 UTC (rev 119)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:19:44 UTC (rev 120)
@@ -101,9 +101,29 @@
"""
- def __init__(self, num=1, den=1):
- """Construct a transfer function. The default is unit static gain."""
+ def __init__(self, *args):
+ """Construct a transfer function.
+
+ The default constructor is TransferFunction(num, den), where num and den
+ are lists of lists of arrays containing polynomial coefficients. To
+ call the copy constructor, call TransferFunction(sys), where sys is a
+ TransferFunction object.
+ """
+
+ if len(args) == 2:
+ # The user provided a numerator and a denominator.
+ (num, den) = args
+ elif len(args) == 1:
+ # Use the copy constructor.
+ if not isinstance(args[0], TransferFunction):
+ raise TypeError("The one-argument constructor can only take in \
+a TransferFunction object. Received %s." % type(args[0]))
+ num = args[0].num
+ den = args[0].den
+ else:
+ raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
+
# Make num and den into lists of lists of arrays, if necessary. Beware:
# this is a shallow copy! This should be okay, but be careful.
data = [num, den]
@@ -211,11 +231,6 @@
data[p][i][j] = data[p][i][j][nonzero:]
[self.num, self.den] = data
- def copy(self):
- """Return a deep copy of the instance."""
-
- return deepcopy(self)
-
def __str__(self):
"""String representation of the transfer function."""
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:19:45
|
Revision: 119
http://python-control.svn.sourceforge.net/python-control/?rev=119&view=rev
Author: kkchen
Date: 2011-02-08 22:19:39 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added rss-balred.py to /examples and bode now returns mag, phase, freq, with the option to supress plotting.
Lauren Padilla <lpa...@pr...>
Modified Paths:
--------------
branches/control-0.4a/examples/secord-matlab.py
branches/control-0.4a/src/freqplot.py
branches/control-0.4a/src/modelsimp.py
branches/control-0.4a/src/xferfcn.py
Added Paths:
-----------
branches/control-0.4a/examples/rss-balred.py
Added: branches/control-0.4a/examples/rss-balred.py
===================================================================
--- branches/control-0.4a/examples/rss-balred.py (rev 0)
+++ branches/control-0.4a/examples/rss-balred.py 2011-02-08 22:19:39 UTC (rev 119)
@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+import numpy as np
+import control.modelsimp as msimp
+import control.matlab as mt
+from control.statesp import StateSpace
+import matplotlib.pyplot as plt
+
+plt.close('all')
+
+#controlable canonical realization computed in matlab for the transfer function:
+# num = [1 11 45 32], den = [1 15 60 200 60]
+A = np.matrix('-15., -7.5, -6.25, -1.875; \
+8., 0., 0., 0.; \
+0., 4., 0., 0.; \
+0., 0., 1., 0.')
+B = np.matrix('2.; 0.; 0.; 0.')
+C = np.matrix('0.5, 0.6875, 0.7031, 0.5')
+D = np.matrix('0.')
+
+# The full system
+fsys = StateSpace(A,B,C,D)
+# The reduced system, truncating the order by 1
+ord = 3
+rsys = msimp.balred(fsys,ord, method = 'truncate')
+
+# Comparison of the step responses of the full and reduced systems
+plt.figure(1)
+t, y = mt.step(fsys)
+tr, yr = mt.step(rsys)
+plt.plot(t.T, y.T)
+plt.hold(True)
+plt.plot(tr.T, yr.T)
+
+# Repeat balanced reduction, now with 100-dimensional random state space
+sysrand = mt.rss(100, 1, 1)
+rsysrand = msimp.balred(sysrand,10,method ='truncate')
+
+# Comparison of the impulse responses of the full and reduced random systems
+plt.figure(2)
+trand, yrand = mt.impulse(sysrand)
+trandr, yrandr = mt.impulse(rsysrand)
+plt.plot(trand.T, yrand.T,trandr.T, yrandr.T)
+
Property changes on: branches/control-0.4a/examples/rss-balred.py
___________________________________________________________________
Added: svn:executable
+ *
Modified: branches/control-0.4a/examples/secord-matlab.py
===================================================================
--- branches/control-0.4a/examples/secord-matlab.py 2011-02-08 22:19:33 UTC (rev 118)
+++ branches/control-0.4a/examples/secord-matlab.py 2011-02-08 22:19:39 UTC (rev 119)
@@ -22,7 +22,7 @@
# Bode plot for the system
figure(2)
-bode(sys, logspace(-2, 2))
+mag,phase,om = bode(sys, logspace(-2, 2),Plot=True)
# Nyquist plot for the system
figure(3)
Modified: branches/control-0.4a/src/freqplot.py
===================================================================
--- branches/control-0.4a/src/freqplot.py 2011-02-08 22:19:33 UTC (rev 118)
+++ branches/control-0.4a/src/freqplot.py 2011-02-08 22:19:39 UTC (rev 119)
@@ -54,12 +54,12 @@
#
# Bode plot
-def bode(syslist, omega=None, dB=False, Hz=False, color=None):
+def bode(syslist, omega=None, dB=False, Hz=False, color=None, Plot=True):
"""Bode plot for a system
Usage
=====
- (magh, phaseh) = bode(syslist, omega=None, dB=False, Hz=False)
+ (magh, phaseh, omega) = bode(syslist, omega=None, dB=False, Hz=False, color=None, Plot=True)
Plots a Bode plot for the system over a (optional) frequency range.
@@ -73,16 +73,18 @@
If True, plot result in dB
Hz : boolean
If True, plot frequency in Hz (omega must be provided in rad/sec)
+ Plot : boolean
+ If True, plot magnitude and phase
Return values
-------------
- magh : graphics handle to magnitude plot (for rescaling, etc)
- phaseh : graphics handle to phase plot
+ magh : magnitude array
+ phaseh : phase array
+ omega : frequency array
Notes
-----
- 1. Use (mag, phase, freq) = sys.freqresp(freq) to generate the
- frequency response for a system.
+ 1. Alternatively, may 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)):
@@ -108,46 +110,47 @@
# 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:
- if color==None:
- plt.semilogx(omega, mag)
+ if (Plot):
+ # 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)")
else:
- plt.semilogx(omega, mag, color=color)
- plt.ylabel("Magnitude (dB)")
- else:
- if color==None:
- plt.loglog(omega, mag)
- else:
- plt.loglog(omega, mag, color=color)
- plt.ylabel("Magnitude")
+ if color==None:
+ plt.loglog(omega, mag)
+ else:
+ plt.loglog(omega, mag, color=color)
+ plt.ylabel("Magnitude")
- # Add a grid to the plot
- plt.grid(True)
- plt.grid(True, which='minor')
- plt.hold(True);
+ # Add a grid to the plot
+ plt.grid(True)
+ plt.grid(True, which='minor')
+ plt.hold(True);
- # Phase plot
- plt.subplot(212);
- if color==None:
- plt.semilogx(omega, phase)
- else:
- plt.semilogx(omega, phase, color=color)
- plt.hold(True)
+ # Phase plot
+ plt.subplot(212);
+ if color==None:
+ plt.semilogx(omega, phase)
+ else:
+ plt.semilogx(omega, phase, color=color)
+ plt.hold(True)
- # Add a grid to the plot
- plt.grid(True)
- plt.grid(True, which='minor')
- plt.ylabel("Phase (deg)")
+ # 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)")
+ # Label the frequency axis
+ if Hz:
+ plt.xlabel("Frequency (Hz)")
+ else:
+ plt.xlabel("Frequency (rad/sec)")
- return (211, 212)
+ return mag, phase, omega
# Nyquist plot
def nyquist(syslist, omega=None):
Modified: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py 2011-02-08 22:19:33 UTC (rev 118)
+++ branches/control-0.4a/src/modelsimp.py 2011-02-08 22:19:39 UTC (rev 119)
@@ -207,6 +207,8 @@
#Check system is stable
D,V = np.linalg.eig(sys.A)
+ print D.shape
+ print D
for e in D:
if e.real >= 0:
raise ValueError, "Oops, the system is unstable!"
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:19:33 UTC (rev 118)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:19:39 UTC (rev 119)
@@ -710,7 +710,8 @@
num[i][j] = list(tfout[6][i, j, :])
# Each transfer function matrix row has a common denominator.
den[i][j] = list(tfout[5][i, :])
-
+ print num
+ print den
return TransferFunction(num, den)
elif isinstance(sys, (int, long, float, complex)):
if "inputs" in kw:
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:19:39
|
Revision: 118
http://python-control.svn.sourceforge.net/python-control/?rev=118&view=rev
Author: kkchen
Date: 2011-02-08 22:19:33 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Back to manually adding tests for python <2.7
My method for finding existing tests required you to be in the build
directory, it didn't work from any directory once installed. It now DOES
work from any directory. The cost though is tests must be manually added to a
list in test.py 'testModules'.
bb...@ra...
Modified Paths:
--------------
branches/control-0.4a/src/test.py
Modified: branches/control-0.4a/src/test.py
===================================================================
--- branches/control-0.4a/src/test.py 2011-02-08 22:19:28 UTC (rev 117)
+++ branches/control-0.4a/src/test.py 2011-02-08 22:19:33 UTC (rev 118)
@@ -75,7 +75,9 @@
print 'Completed tests in',mod
except:
#find the modules ourselves without unittest discovery
- testModules = findTests()
+ #testModules = findTests()
+ testModules =['TestBDAlg','TestConvert','TestFreqRsp','TestMatlab','TestModelsimp',\
+ 'TestSlycot','TestStatefbk','TestStateSp','TestXferFcn']
print 'Tests may be incomplete, will attempt to run tests in modules:'
for m in testModules:
print m
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:19:34
|
Revision: 117
http://python-control.svn.sourceforge.net/python-control/?rev=117&view=rev
Author: kkchen
Date: 2011-02-08 22:19:28 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
minor changes to test.py
bb...@ra...
Modified Paths:
--------------
branches/control-0.4a/src/test.py
Modified: branches/control-0.4a/src/test.py
===================================================================
--- branches/control-0.4a/src/test.py 2011-02-08 22:19:23 UTC (rev 116)
+++ branches/control-0.4a/src/test.py 2011-02-08 22:19:28 UTC (rev 117)
@@ -56,7 +56,7 @@
###############################################################################
-def tests():
+def tests(verbosity=2):
import unittest
try: #auto test discovery is only implemented in python 2.7+
start_dir='./' #change to a tests directory eventually.
@@ -74,19 +74,13 @@
t.run(tests)
print 'Completed tests in',mod
except:
- #If can't do auto discovery, for now it is hard-coded. This is not ideal for
- #when new tests are added or existing ones are reorganized/renamed.
-
- #remove all of the print commands once tests are debugged and converted to
- #unittests.
-
- print 'Tests may be incomplete'
- t=unittest.TextTestRunner()
-
+ #find the modules ourselves without unittest discovery
testModules = findTests()
+ print 'Tests may be incomplete, will attempt to run tests in modules:'
+ for m in testModules:
+ print m
- suite = unittest.TestSuite()
- suiteList=[]
+ suiteList=[] #list of unittest.TestSuite objects
for mod in testModules:
exec('import '+mod+' as currentModule')
#After tests have been debugged and made into unittests, remove
@@ -101,8 +95,8 @@
print 'The test module '+mod+' doesnt have '+\
'a proper suite() function that returns a unittest.TestSuite object'+\
' Please fix this!'
- alltests = unittest.TestSuite(suiteList)
- t.run(unittest.TestSuite(alltests))
+ t=unittest.TextTestRunner(verbosity=verbosity)
+ t.run(unittest.TestSuite(unittest.TestSuite(suiteList)))
###############################################################################
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
7 messages has been excluded from this view by a project administrator.
