Menu
▾
▴
[python-control-discuss] SF.net SVN: python-control:[168] trunk
|
From: <mur...@us...> - 2011-07-16 06:05:52
|
Revision: 168
http://python-control.svn.sourceforge.net/python-control/?rev=168&view=rev
Author: murrayrm
Date: 2011-07-16 06:05:46 +0000 (Sat, 16 Jul 2011)
Log Message:
-----------
* Moved margin command from freqplot.py to margins.py (new file); updated
matlab.py, init.py, etc to reflect changes
* Added PhaseCrossoverFrequencies() function; contributed by Steffen Waldherr
* Added initial unit tests for StabilityMargins(), PhaseCrossoverFrequencies()
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/__init__.py
trunk/src/freqplot.py
trunk/src/matlab.py
trunk/src/statesp.py
trunk/tests/matlab_test.py
Added Paths:
-----------
trunk/src/margins.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2011-07-13 16:03:32 UTC (rev 167)
+++ trunk/ChangeLog 2011-07-16 06:05:46 UTC (rev 168)
@@ -1,3 +1,24 @@
+2011-07-15 Richard Murray <murray@malabar.local>
+
+ * tests/matlab_test.py (TestMatlab): added unittest for margin()
+ commands (calling format only)
+
+ * src/statesp.py (StateSpace): updated comments
+
+ * tests/margin_test.py: set up unit tests for StabilityMargins() and
+ PhaseCrossoverFrequencies()
+
+ * src/__init__.py: added margins.py to __init__
+
+2011-07-14 Richard Murray <murray@malabar.local>
+
+ * src/margins.py (GainPhaseMargin): moved freqplot.MarginPlot to
+ margin.StabilityMargins (MarginPlot didn't actually plot anything)
+
+ * src/margins.py (PhaseCrossoverFrequencies): added new function to
+ compute frequencies that we cross real axis. Contributed by Steffen
+ Waldherr <wal...@is...>
+
2011-07-11 Richard Murray <murray@malabar.local>
* src/rlocus.py: added real() and imag() to list of functions
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2011-07-13 16:03:32 UTC (rev 167)
+++ trunk/src/__init__.py 2011-07-16 06:05:46 UTC (rev 168)
@@ -57,15 +57,16 @@
# Import functions from within the control system library
#! Should probably only import the exact functions we use...
-from xferfcn import *
-from statesp import *
-from freqplot import *
-from nichols import *
from bdalg import *
-from statefbk import *
from delay import *
+from freqplot import *
+from margins import *
+from mateqn import *
from modelsimp import *
+from nichols import *
from rlocus import *
-from mateqn import *
+from statefbk import *
+from statesp import *
from timeresp import ForcedResponse, InitialResponse, StepResponse, \
ImpulseResponse
+from xferfcn import *
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2011-07-13 16:03:32 UTC (rev 167)
+++ trunk/src/freqplot.py 2011-07-16 06:05:46 UTC (rev 168)
@@ -299,101 +299,6 @@
phase = np.squeeze(phase_tmp)
plt.subplot(224); plt.loglog(omega, mag);
-
-# gain and phase margins
-# contributed by Sawyer B. Fuller <mi...@ca...>
-def MarginPlot(sysdata, deg=True):
- """Calculate gain and phase margins and associated crossover frequencies
-
- Usage:
-
- gm, pm, sm, wg, wp, ws = margin(sysdata, deg=True)
-
- Parameters
- ----------
- sysdata: linsys or (mag, phase, omega) sequence
- sys : linsys
- Linear SISO system
- mag, phase, omega : sequence of array_like
- Input magnitude, phase, and frequencies (rad/sec) sequence from
- bode frequency response data
- deg=True: boolean
- If true, all input and output phases in degrees, else in radians
-
- Returns
- -------
- gm, pm, sm, wg, wp, ws: float
- Gain margin gm, phase margin pm, stability margin sm, and
- associated crossover
- frequencies wg, wp, and ws of SISO open-loop. If more than
- one crossover frequency is detected, returns the lowest corresponding
- margin.
- """
- #TODO do this precisely without the effects of discretization of frequencies?
- #TODO assumes SISO
- #TODO unit tests, margin plot
-
- if (not getattr(sysdata, '__iter__', False)):
- sys = sysdata
- mag, phase, omega = bode(sys, deg=deg, Plot=False)
- elif len(sysdata) == 3:
- mag, phase, omega = sysdata
- else:
- raise ValueError("Margin sysdata must be either a linear system or a 3-sequence of mag, phase, omega.")
-
- if deg:
- cycle = 360.
- crossover = 180.
- else:
- cycle = 2 * np.pi
- crossover = np.pi
-
- wrapped_phase = -np.mod(phase, cycle)
-
- # phase margin from minimum phase among all gain crossovers
- neg_mag_crossings_i = np.nonzero(np.diff(mag < 1) > 0)[0]
- mag_crossings_p = wrapped_phase[neg_mag_crossings_i]
- if len(neg_mag_crossings_i) == 0:
- if mag[0] < 1: # gain always less than one
- wp = np.nan
- pm = np.inf
- else: # gain always greater than one
- print "margin: no magnitude crossings found"
- wp = np.nan
- pm = np.nan
- else:
- min_mag_crossing_i = neg_mag_crossings_i[np.argmin(mag_crossings_p)]
- wp = omega[min_mag_crossing_i]
- pm = crossover + phase[min_mag_crossing_i]
- if pm < 0:
- print "warning: system unstable: negative phase margin"
-
- # gain margin from minimum gain margin among all phase crossovers
- neg_phase_crossings_i = np.nonzero(np.diff(wrapped_phase < -crossover) > 0)[0]
- neg_phase_crossings_g = mag[neg_phase_crossings_i]
- if len(neg_phase_crossings_i) == 0:
- wg = np.nan
- gm = np.inf
- else:
- min_phase_crossing_i = neg_phase_crossings_i[
- np.argmax(neg_phase_crossings_g)]
- wg = omega[min_phase_crossing_i]
- gm = abs(1/mag[min_phase_crossing_i])
- if gm < 1:
- print "warning: system unstable: gain margin < 1"
-
- # stability margin from minimum abs distance from -1 point
- if deg:
- phase_rad = phase * np.pi / 180.
- else:
- phase_rad = phase
- L = mag * np.exp(1j * phase_rad) # complex loop response to -1 pt
- min_Lplus1_i = np.argmin(np.abs(L + 1))
- sm = np.abs(L[min_Lplus1_i] + 1)
- ws = phase[min_Lplus1_i]
-
- return gm, pm, sm, wg, wp, ws
-
#
# Utility functions
#
@@ -461,4 +366,3 @@
bode = BodePlot
nyquist = NyquistPlot
gangof4 = GangOf4Plot
-margin = MarginPlot
Added: trunk/src/margins.py
===================================================================
--- trunk/src/margins.py (rev 0)
+++ trunk/src/margins.py 2011-07-16 06:05:46 UTC (rev 168)
@@ -0,0 +1,193 @@
+"""margin.py
+
+Functions for computing stability margins and related functions.
+
+Routeins in this module:
+
+margin.StabilityMargins
+margin.PhaseCrossoverFrequencies
+"""
+
+"""Copyright (c) 2011 by California Institute of Technology
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the California Institute of Technology nor
+ the names of its contributors may be used to endorse or promote
+ products derived from this software without specific prior
+ written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
+OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
+
+Author: Richard M. Murray
+Date: 14 July 2011
+
+$Id: xferfcn.py 165 2011-06-26 02:44:09Z murrayrm $
+
+"""
+
+import xferfcn
+from freqplot import bode
+import numpy as np
+
+# gain and phase margins
+# contributed by Sawyer B. Fuller <mi...@ca...>
+#! TODO - need to add unit test functions
+def StabilityMargins(sysdata, deg=True):
+ """Calculate gain, phase and stability margins and associated
+ crossover frequencies.
+
+ Usage:
+
+ gm, pm, sm, wg, wp, ws = StabilityMargins(sysdata, deg=True)
+
+ Parameters
+ ----------
+ sysdata: linsys or (mag, phase, omega) sequence
+ sys : linsys
+ Linear SISO system
+ mag, phase, omega : sequence of array_like
+ Input magnitude, phase, and frequencies (rad/sec) sequence from
+ bode frequency response data
+ deg=True: boolean
+ If true, all input and output phases in degrees, else in radians
+
+ Returns
+ -------
+ gm, pm, sm, wg, wp, ws: float
+ Gain margin gm, phase margin pm, stability margin sm, and
+ associated crossover
+ frequencies wg, wp, and ws of SISO open-loop. If more than
+ one crossover frequency is detected, returns the lowest corresponding
+ margin.
+ """
+ #TODO do this precisely without the effects of discretization of frequencies?
+ #TODO assumes SISO
+ #TODO unit tests, margin plot
+
+ if (not getattr(sysdata, '__iter__', False)):
+ sys = sysdata
+ mag, phase, omega = bode(sys, deg=deg, Plot=False)
+ elif len(sysdata) == 3:
+ mag, phase, omega = sysdata
+ else:
+ raise ValueError("Margin sysdata must be either a linear system or a 3-sequence of mag, phase, omega.")
+
+ if deg:
+ cycle = 360.
+ crossover = 180.
+ else:
+ cycle = 2 * np.pi
+ crossover = np.pi
+
+ wrapped_phase = -np.mod(phase, cycle)
+
+ # phase margin from minimum phase among all gain crossovers
+ neg_mag_crossings_i = np.nonzero(np.diff(mag < 1) > 0)[0]
+ mag_crossings_p = wrapped_phase[neg_mag_crossings_i]
+ if len(neg_mag_crossings_i) == 0:
+ if mag[0] < 1: # gain always less than one
+ wp = np.nan
+ pm = np.inf
+ else: # gain always greater than one
+ print "margin: no magnitude crossings found"
+ wp = np.nan
+ pm = np.nan
+ else:
+ min_mag_crossing_i = neg_mag_crossings_i[np.argmin(mag_crossings_p)]
+ wp = omega[min_mag_crossing_i]
+ pm = crossover + phase[min_mag_crossing_i]
+ if pm < 0:
+ print "warning: system unstable: negative phase margin"
+
+ # gain margin from minimum gain margin among all phase crossovers
+ neg_phase_crossings_i = np.nonzero(np.diff(wrapped_phase < -crossover) > 0)[0]
+ neg_phase_crossings_g = mag[neg_phase_crossings_i]
+ if len(neg_phase_crossings_i) == 0:
+ wg = np.nan
+ gm = np.inf
+ else:
+ min_phase_crossing_i = neg_phase_crossings_i[
+ np.argmax(neg_phase_crossings_g)]
+ wg = omega[min_phase_crossing_i]
+ gm = abs(1/mag[min_phase_crossing_i])
+ if gm < 1:
+ print "warning: system unstable: gain margin < 1"
+
+ # stability margin from minimum abs distance from -1 point
+ if deg:
+ phase_rad = phase * np.pi / 180.
+ else:
+ phase_rad = phase
+ L = mag * np.exp(1j * phase_rad) # complex loop response to -1 pt
+ min_Lplus1_i = np.argmin(np.abs(L + 1))
+ sm = np.abs(L[min_Lplus1_i] + 1)
+ ws = phase[min_Lplus1_i]
+
+ return gm, pm, sm, wg, wp, ws
+
+# Contributed by Steffen Waldherr <wal...@is...>
+#! TODO - need to add test functions
+def PhaseCrossoverFrequencies(sys):
+ """
+ Compute frequencies and gains at intersections with real axis
+ in Nyquist plot.
+
+ Call as:
+ omega, gain = PhaseCrossoverFrequencies()
+
+ Returns
+ -------
+ omega: 1d array of (non-negative) frequencies where Nyquist plot
+ intersects the real axis
+
+ gain: 1d array of corresponding gains
+
+ Examples
+ --------
+ >>> tf = TransferFunction([1], [1, 2, 3, 4])
+ >>> PhaseCrossoverFrequenies(tf)
+ (array([ 1.73205081, 0. ]), array([-0.5 , 0.25]))
+ """
+
+ # Convert to a transfer function
+ tf = xferfcn._convertToTransferFunction(sys)
+
+ # if not siso, fall back to (0,0) element
+ #! TODO: should add a check and warning here
+ num = tf.num[0][0]
+ den = tf.den[0][0]
+
+ # Compute frequencies that we cross over the real axis
+ numj = (1.j)**np.arange(len(num)-1,-1,-1)*num
+ denj = (-1.j)**np.arange(len(den)-1,-1,-1)*den
+ allfreq = np.roots(np.imag(np.polymul(numj,denj)))
+ realfreq = np.real(allfreq[np.isreal(allfreq)])
+ realposfreq = realfreq[realfreq >= 0.]
+
+ # using real() to avoid rounding errors and results like 1+0j
+ # it would be nice to have a vectorized version of self.evalfr here
+ gain = np.real(np.asarray([tf.evalfr(f)[0][0] for f in realposfreq]))
+
+ return realposfreq, gain
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2011-07-13 16:03:32 UTC (rev 167)
+++ trunk/src/matlab.py 2011-07-16 06:05:46 UTC (rev 168)
@@ -72,6 +72,7 @@
import ctrlutil
import freqplot
import timeresp
+import margins
from statesp import StateSpace, _rss_generate, _convertToStateSpace
from xferfcn import TransferFunction, _convertToTransferFunction
from lti import Lti #base class of StateSpace, TransferFunction
@@ -1025,7 +1026,6 @@
rlist = RootLocus(sys, klist, **keywords)
return rlist, klist
-
def margin(*args):
"""Calculate gain and phase margins and associated crossover frequencies
@@ -1060,16 +1060,15 @@
"""
if len(args) == 1:
sys = args[0]
- margins = freqplot.margin(sys)
+ margin = margins.StabilityMargins(sys)
elif len(args) == 3:
- margins = freqplot.margin(args)
+ margin = margins.StabilityMargins(args)
else:
raise ValueError("Margin needs 1 or 3 arguments; received %i."
% len(args))
- return margins[0], margins[1], margins[3], margins[4]
+ return margin[0], margin[1], margin[3], margin[4]
-
def dcgain(*args):
'''
Compute the gain of the system in steady state.
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2011-07-13 16:03:32 UTC (rev 167)
+++ trunk/src/statesp.py 2011-07-16 06:05:46 UTC (rev 168)
@@ -208,7 +208,7 @@
return StateSpace(self.A, self.B, -self.C, -self.D)
- # Addition of two transfer functions (parallel interconnection)
+ # Addition of two state space systems (parallel interconnection)
def __add__(self, other):
"""Add two LTI systems (parallel connection)."""
@@ -244,7 +244,7 @@
return self + other
- # Subtraction of two transfer functions (parallel interconnection)
+ # Subtraction of two state space systems (parallel interconnection)
def __sub__(self, other):
"""Subtract two LTI systems."""
@@ -255,7 +255,7 @@
return other + (-self)
- # Multiplication of two transfer functions (series interconnection)
+ # Multiplication of two state space systems (series interconnection)
def __mul__(self, other):
"""Multiply two LTI objects (serial connection)."""
@@ -284,7 +284,7 @@
return StateSpace(A, B, C, D)
- # Right multiplication of two transfer functions (series interconnection)
+ # Right multiplication of two state space systems (series interconnection)
# Just need to convert LH argument to a state space object
def __rmul__(self, other):
"""Right multiply two LTI objects (serial connection)."""
Modified: trunk/tests/matlab_test.py
===================================================================
--- trunk/tests/matlab_test.py 2011-07-13 16:03:32 UTC (rev 167)
+++ trunk/tests/matlab_test.py 2011-07-16 06:05:46 UTC (rev 168)
@@ -186,6 +186,13 @@
[-0.1391, 48.9776]])
yout, _t, _xout = lsim(self.mimo_ss1, u, t, x0)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
+
+ def testMargin(self):
+ #! TODO: check results to make sure they are OK
+ gm, pm, wg, wp = margin(self.siso_tf1);
+ gm, pm, wg, wp = margin(self.siso_tf2);
+ gm, pm, wg, wp = margin(self.siso_ss1);
+ gm, pm, wg, wp = margin(self.siso_ss2);
def testDcgain(self):
#Create different forms of a SISO system
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
