Menu
▾
▴
python-control-discuss
|
From: <mur...@us...> - 2011-06-26 02:44:16
|
Revision: 165
http://python-control.svn.sourceforge.net/python-control/?rev=165&view=rev
Author: murrayrm
Date: 2011-06-26 02:44:09 +0000 (Sun, 26 Jun 2011)
Log Message:
-----------
* Made fixes to documentation (intro, matlab)
* Got rid of unit test error in hsvd (convert slycot args to numpy array)
* Tracked down possible errors in convert_test (still in progress)
See ChangeLog for more detailed list of code that was modified
Modified Paths:
--------------
trunk/ChangeLog
trunk/doc/intro.rst
trunk/src/matlab.py
trunk/src/statefbk.py
trunk/src/timeresp.py
trunk/src/xferfcn.py
trunk/tests/convert_test.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2011-06-22 19:26:01 UTC (rev 164)
+++ trunk/ChangeLog 2011-06-26 02:44:09 UTC (rev 165)
@@ -1,3 +1,27 @@
+2011-06-25 Richard Murray <murray@malabar.local>
+
+ * src/xferfcn.py (TransferFunction._common_den): changed tolerance
+ for detecting complex valued poles to a user-settable parameter,
+ with default value 1e-8. This was an attempt to fix errors in the
+ convert_test.py unittest script (conversion routine was
+ misclassifying some poles as imaginary when they weren't).
+
+ * src/xferfcn.py (_convertToTransferFunction): converted arguments
+ to tb04ad to numpy arrays; fixes a unit test error in convert_test.py.
+
+ * src/statefbk.py (gram): convert system matrix passed to sb03md to
+ numpy array; this fixes a unit test error in modelsimp_test.py.
+
+ * src/matlab.py (impulse): got rid of X0 argument for impulse
+ response (not implemented in MATLAB).
+
+ * doc/intro.rst: added some quick start information
+
+ * src/matlab.py: added documentation for step, impulse, initial, lsim
+
+ * src/timeresp.py: fixed some MATLAB specific function names in
+ function doc strings
+
2011-06-22 Richard Murray <murray@malabar.local>
* doc/intro.rst: fixed some small types
Modified: trunk/doc/intro.rst
===================================================================
--- trunk/doc/intro.rst 2011-06-22 19:26:01 UTC (rev 164)
+++ trunk/doc/intro.rst 2011-06-26 02:44:09 UTC (rev 165)
@@ -8,9 +8,10 @@
that implement common operations for the analysis and design of
feedback control systems. The initial goal is to implement all of the
functionality required to work through the examples in the textbook
-Feedback Systems by \xC5str\xF6m and Murray. A MATLAB compatibility package
-(control.matlab) is available that provides functions corresponding to
-the commands available in the MATLAB Control Systems Toolbox.
+Feedback Systems by Astrom and Murray. A MATLAB compatibility package
+(control.matlab) is available that provides many of the common
+functions corresponding to commands available in the MATLAB Control
+Systems Toolbox.
In addition to the documentation here, there is a project wiki that
contains some additional information about how to use the package
@@ -26,3 +27,19 @@
* Transfer functions are only implemented for SISO systems (due to
limitations in the underlying signals.lti class); use state space
representations for MIMO systems.
+
+Getting Started
+---------------
+1. Download the latest release from http:sf.net/projects/python-control/files.
+2. Untar the source code in a temporary directory and run 'python setup.py
+ install' to build and install the code
+3. To see if things are working correctly, run ipython -pylab and run the
+ script 'examples/secord-matlab.py'. This should generate a set response,
+ Bode plot and Nyquist plot for a simple second order system.
+4. To see the commands that are available, run the following commands in
+ ipython::
+ >>> import control
+ >>> ?control.matlab
+5. If you want to have a MATLAB-like environment for running the control
+ toolbox, use::
+ >>> from control.matlab import *
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2011-06-22 19:26:01 UTC (rev 164)
+++ trunk/src/matlab.py 2011-06-26 02:44:09 UTC (rev 165)
@@ -1128,21 +1128,217 @@
# Call corresponding functions in timeresp, with arguments transposed
def step(sys, T=None, X0=0., input=0, output=0, **keywords):
+ '''
+ Step response of a linear system
+
+ If the system has multiple inputs or outputs (MIMO), one input and one
+ output have to be selected for the simulation. The parameters `input`
+ and `output` do this. All other inputs are set to 0, all other outputs
+ are ignored.
+
+ Parameters
+ ----------
+ sys: StateSpace, or TransferFunction
+ LTI system to simulate
+
+ T: array-like object, optional
+ Time vector (argument is autocomputed if not given)
+
+ X0: array-like or number, optional
+ Initial condition (default = 0)
+
+ Numbers are converted to constant arrays with the correct shape.
+
+ input: int
+ Index of the input that will be used in this simulation.
+
+ output: int
+ Index of the output that will be used in this simulation.
+
+ **keywords:
+ Additional keyword arguments control the solution algorithm for the
+ differential equations. These arguments are passed on to the function
+ :func:`control.ForcedResponse`, which in turn passes them on to
+ :func:`scipy.integrate.odeint`. See the documentation for
+ :func:`scipy.integrate.odeint` for information about these
+ arguments.
+
+ Returns
+ -------
+ T: array
+ Time values of the output
+
+ yout: array
+ Response of the system
+
+ See Also
+ --------
+ lsim, initial, impulse
+
+ Examples
+ --------
+ >>> T, yout = step(sys, T, X0)
+ '''
T, yout = timeresp.StepResponse(sys, T, X0, input, output,
transpose = True, **keywords)
return T, yout
-def impulse(sys, T=None, X0=0., input=0, output=0, **keywords):
- T, yout = timeresp.ImpulseResponse(sys, T, X0, input, output,
+def impulse(sys, T=None, input=0, output=0, **keywords):
+ '''
+ Impulse response of a linear system
+
+ If the system has multiple inputs or outputs (MIMO), one input and
+ one output must be selected for the simulation. The parameters
+ `input` and `output` do this. All other inputs are set to 0, all
+ other outputs are ignored.
+
+ Parameters
+ ----------
+ sys: StateSpace, TransferFunction
+ LTI system to simulate
+
+ T: array-like object, optional
+ Time vector (argument is autocomputed if not given)
+
+ input: int
+ Index of the input that will be used in this simulation.
+
+ output: int
+ Index of the output that will be used in this simulation.
+
+ **keywords:
+ Additional keyword arguments control the solution algorithm for the
+ differential equations. These arguments are passed on to the function
+ :func:`lsim`, which in turn passes them on to
+ :func:`scipy.integrate.odeint`. See the documentation for
+ :func:`scipy.integrate.odeint` for information about these
+ arguments.
+
+ Returns
+ -------
+ T: array
+ Time values of the output
+ yout: array
+ Response of the system
+
+ See Also
+ --------
+ lsim, step, initial
+
+ Examples
+ --------
+ >>> T, yout = impulse(sys, T)
+ '''
+ T, yout = timeresp.ImpulseResponse(sys, T, 0, input, output,
transpose = True, **keywords)
return T, yout
def initial(sys, T=None, X0=0., input=0, output=0, **keywords):
+ '''
+ Initial condition response of a linear system
+
+ If the system has multiple inputs or outputs (MIMO), one input and one
+ output have to be selected for the simulation. The parameters `input`
+ and `output` do this. All other inputs are set to 0, all other outputs
+ are ignored.
+
+ Parameters
+ ----------
+ sys: StateSpace, or TransferFunction
+ LTI system to simulate
+
+ T: array-like object, optional
+ Time vector (argument is autocomputed if not given)
+
+ X0: array-like object or number, optional
+ Initial condition (default = 0)
+
+ Numbers are converted to constant arrays with the correct shape.
+
+ input: int
+ Index of the input that will be used in this simulation.
+
+ output: int
+ Index of the output that will be used in this simulation.
+
+ **keywords:
+ Additional keyword arguments control the solution algorithm for the
+ differential equations. These arguments are passed on to the function
+ :func:`lsim`, which in turn passes them on to
+ :func:`scipy.integrate.odeint`. See the documentation for
+ :func:`scipy.integrate.odeint` for information about these
+ arguments.
+
+
+ Returns
+ -------
+ T: array
+ Time values of the output
+ yout: array
+ Response of the system
+
+ See Also
+ --------
+ lsim, step, impulse
+
+ Examples
+ --------
+ >>> T, yout = initial(sys, T, X0)
+ '''
T, yout = timeresp.InitialResponse(sys, T, X0, input, output,
transpose = True, **keywords)
return T, yout
def lsim(sys, U=0., T=None, X0=0., **keywords):
+ '''
+ Simulate the output of a linear system.
+
+ As a convenience for parameters `U`, `X0`:
+ Numbers (scalars) are converted to constant arrays with the correct shape.
+ The correct shape is inferred from arguments `sys` and `T`.
+
+ Parameters
+ ----------
+ sys: Lti (StateSpace, or TransferFunction)
+ LTI system to simulate
+
+ U: array-like or number, optional
+ Input array giving input at each time `T` (default = 0).
+
+ If `U` is ``None`` or ``0``, a special algorithm is used. This special
+ algorithm is faster than the general algorithm, which is used otherwise.
+
+ T: array-like
+ Time steps at which the input is defined, numbers must be (strictly
+ monotonic) increasing.
+
+ X0: array-like or number, optional
+ Initial condition (default = 0).
+
+ **keywords:
+ Additional keyword arguments control the solution algorithm for the
+ differential equations. These arguments are passed on to the function
+ :func:`scipy.integrate.odeint`. See the documentation for
+ :func:`scipy.integrate.odeint` for information about these
+ arguments.
+
+ Returns
+ -------
+ T: array
+ Time values of the output.
+ yout: array
+ Response of the system.
+ xout: array
+ Time evolution of the state vector.
+
+ See Also
+ --------
+ step, initial, impulse
+
+ Examples
+ --------
+ >>> T, yout, xout = lsim(sys, U, T, X0)
+ '''
T, yout, xout = timeresp.ForcedResponse(sys, T, U, X0,
transpose = True, **keywords)
return T, yout, xout
Modified: trunk/src/statefbk.py
===================================================================
--- trunk/src/statefbk.py 2011-06-22 19:26:01 UTC (rev 164)
+++ trunk/src/statefbk.py 2011-06-26 02:44:09 UTC (rev 165)
@@ -328,7 +328,7 @@
raise ControlSlycot("can't find slycot module 'sb03md'")
n = sys.states
U = np.zeros((n,n))
- A = sys.A
+ A = np.array(sys.A) # convert to NumPy array for slycot
X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
gram = X
return gram
Modified: trunk/src/timeresp.py
===================================================================
--- trunk/src/timeresp.py 2011-06-22 19:26:01 UTC (rev 164)
+++ trunk/src/timeresp.py 2011-06-26 02:44:09 UTC (rev 165)
@@ -419,11 +419,11 @@
See Also
--------
- lsim, initial, impulse
+ ForcedResponse, InitialResponse, ImpulseResponse
Examples
--------
- >>> T, yout = step(sys, T, X0)
+ >>> T, yout = StepResponse(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
sys = _mimo2siso(sys, input, output, warn_conversion=True)
@@ -491,11 +491,11 @@
See Also
--------
- lsim, step, impulse
+ ForcedResponse, ImpulseResponse, StepResponse
Examples
--------
- >>> T, yout = InitialResponsesys, T, X0)
+ >>> T, yout = InitialResponse(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
sys = _mimo2siso(sys, input, output, warn_conversion=True)
@@ -564,7 +564,7 @@
See Also
--------
- lsim, step, initial
+ ForcedReponse, InitialResponse, StepResponse
Examples
--------
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2011-06-22 19:26:01 UTC (rev 164)
+++ trunk/src/xferfcn.py 2011-06-26 02:44:09 UTC (rev 165)
@@ -527,22 +527,46 @@
return out
- def _common_den(self):
- """Compute MIMO common denominator; return it and an adjusted numerator.
+ def _common_den(self, imag_tol=None):
+ """
+ Compute MIMO common denominator; return it and an adjusted numerator.
+ This function computes the single denominator containing all
+ the poles of sys.den, and reports it as the array d. The
+ output numerator array n is modified to use the common
+ denominator; the coefficient arrays are also padded with zeros
+ to be the same size as d. n is an sys.outputs by sys.inputs
+ by len(d) array.
+
+ Parameters
+ ----------
+ imag_tol: float
+ Threshold for the imaginary part of a root to use in detecting
+ complex poles
+
+ Returns
+ -------
+ num: array
+ Multi-dimensional array of numerator coefficients. num[i][j]
+ gives the numerator coefficient array for the ith input and jth
+ output
+
+ den: array
+ Array of coefficients for common denominator polynomial
+
+ Examples
+ --------
>>> n, d = sys._common_den()
-
- computes the single denominator containing all the poles of sys.den, and
- reports it as the array d. The output numerator array n is modified to
- use the common denominator; the coefficient arrays are also padded with
- zeros to be the same size as d. n is an sys.outputs-by-sys.inputs-by-
- len(d) array.
-
+
"""
# Machine precision for floats.
eps = finfo(float).eps
+ # Decide on the tolerance to use in deciding of a pole is complex
+ if (imag_tol == None):
+ imag_tol = 1e-8 #! TODO: figure out the right number to use
+
# A sorted list to keep track of cumulative poles found as we scan
# self.den.
poles = []
@@ -615,8 +639,9 @@
# To prevent buildup of imaginary part error, handle complex
# pole pairs together.
quad = polymul([1., -poles[n]], [1., -poles[n+1]])
- assert all(quad.imag < eps), "The quadratic has a nontrivial \
-imaginary part: %g" % quad.imag.max()
+ assert all(quad.imag < 10 * eps), \
+ "The quadratic has a nontrivial imaginary part: %g" \
+ % quad.imag.max()
quad = quad.real
den = polymul(den, quad)
@@ -747,10 +772,10 @@
raise TypeError("If sys is a StateSpace, _convertToTransferFunction \
cannot take keywords.")
- # Use Slycot to make the transformation. TODO: this is still somewhat
- # buggy!
- tfout = tb04ad(sys.states, sys.inputs, sys.outputs, sys.A, sys.B, sys.C,
- sys.D,tol1=0.0)
+ # Use Slycot to make the transformation
+ # Make sure to convert system matrices to numpy arrays
+ tfout = tb04ad(sys.states, sys.inputs, sys.outputs, array(sys.A),
+ array(sys.B), array(sys.C), array(sys.D), tol1=0.0)
# Preallocate outputs.
num = [[[] for j in range(sys.inputs)] for i in range(sys.outputs)]
Modified: trunk/tests/convert_test.py
===================================================================
--- trunk/tests/convert_test.py 2011-06-22 19:26:01 UTC (rev 164)
+++ trunk/tests/convert_test.py 2011-06-26 02:44:09 UTC (rev 165)
@@ -29,9 +29,9 @@
# Maximum number of states to test + 1
self.maxStates = 20
# Maximum number of inputs and outputs to test + 1
- self.maxIO = 20
+ self.maxIO = 10
# Set to True to print systems to the output.
- self.debug = True
+ self.debug = False
def printSys(self, sys, ind):
"""Print system to the standard output."""
@@ -40,9 +40,10 @@
print "sys%i:\n" % ind
print sys
- def testConvert(self, verbose=0):
+ def testConvert(self):
"""Test state space to transfer function conversion."""
#Currently it only tests that a TF->SS->TF generates an unchanged TF
+ verbose = self.debug
#print __doc__
@@ -70,12 +71,14 @@
for inputNum in range(inputs):
for outputNum in range(outputs):
np.testing.assert_array_almost_equal(\
- tfOriginal.num[outputNum][inputNum],\
- tfTransformed.num[outputNum][inputNum])
+ tfOriginal.num[outputNum][inputNum], \
+ tfTransformed.num[outputNum][inputNum], \
+ err_msg='numerator mismatch')
np.testing.assert_array_almost_equal(\
- tfOriginal.den[outputNum][inputNum],\
- tfTransformed.den[outputNum][inputNum])
+ tfOriginal.den[outputNum][inputNum], \
+ tfTransformed.den[outputNum][inputNum],
+ err_msg='denominator mismatch')
#To test the ss systems is harder because they aren't the same
#realization. This could be done with checking that they have the
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2011-07-25 05:00:48
|
Revision: 169
http://python-control.svn.sourceforge.net/python-control/?rev=169&view=rev
Author: murrayrm
Date: 2011-07-25 05:00:41 +0000 (Mon, 25 Jul 2011)
Log Message:
-----------
* Added first cut at phase portrait command (PhasePlot) from MATLAB
code; needs more work to pythonize the call signatures
* Added simple unit tests for margin command
See ChangeLog for more detailed list of changes
Modified Paths:
--------------
trunk/ChangeLog
Added Paths:
-----------
trunk/examples/README
trunk/examples/phaseplots.py
trunk/src/phaseplot.py
trunk/tests/margin_test.py
trunk/tests/phaseplot_test.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2011-07-16 06:05:46 UTC (rev 168)
+++ trunk/ChangeLog 2011-07-25 05:00:41 UTC (rev 169)
@@ -1,3 +1,17 @@
+2011-07-24 Richard Murray <murray@malabar.local>
+
+ * tests/margin_test.py: added simple unit tests for margin functions
+ (initial versions just call functions; some comparisons missing)
+
+ * examples/README: added missing README file
+
+ * examples/phaseplots.py: FBS examples for phaseplot
+
+ * tests/phaseplot_test.py: unit tests for phaseplot
+
+ * src/phaseplot.py: initial cut at phase portrait function, built
+ from amphaseplot (Feeback Systems [FBS], Astrom and Murray, 2008)
+
2011-07-15 Richard Murray <murray@malabar.local>
* tests/matlab_test.py (TestMatlab): added unittest for margin()
Added: trunk/examples/README
===================================================================
--- trunk/examples/README (rev 0)
+++ trunk/examples/README 2011-07-25 05:00:41 UTC (rev 169)
@@ -0,0 +1,4 @@
+This directory contains worked examples using the python-control
+library. Each example should work by running 'ipython -pylab' and
+then running the given file (make sure to have the python-control
+module in your path).
Added: trunk/examples/phaseplots.py
===================================================================
--- trunk/examples/phaseplots.py (rev 0)
+++ trunk/examples/phaseplots.py 2011-07-25 05:00:41 UTC (rev 169)
@@ -0,0 +1,120 @@
+# phaseplots.py - examples of phase portraits
+# RMM, 24 July 2011
+#
+# This file contains examples of phase portraits pulled from "Feedback
+# Systems" by Astrom and Murray (Princeton University Press, 2008).
+
+import numpy as np
+import matplotlib.pyplot as mpl
+from control.phaseplot import PhasePlot
+from numpy import pi
+
+# Clear out any figures that are present
+mpl.close('all')
+
+#
+# Inverted pendulum
+#
+
+# Define the ODEs for a damped (inverted) pendulum
+def invpend_ode(x, t, m=1., l=1., b=0.2, g=1):
+ return (x[1], -b/m*x[1] + (g*l/m) * np.sin(x[0]))
+
+# Set up the figure the way we want it to look
+mpl.figure(); mpl.clf();
+mpl.axis([-2*pi, 2*pi, -2.1, 2.1]);
+mpl.title('Inverted pendlum')
+
+# Outer trajectories
+PhasePlot(invpend_ode,
+ 'logtime', (3, 0.7), None,
+ [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
+ [-1, 2.1], [4.2, 2.1], [5, 2.1],
+ [2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
+ [1, -2.1], [-4.2, -2.1], [-5, -2.1] ],
+ np.linspace(0, 40, 200))
+
+# Separatrices
+mpl.hold(True);
+PhasePlot(invpend_ode, 'auto', 0, None, [[-2.3056, 2.1], [2.3056, -2.1]], 6)
+mpl.show();
+
+#
+# Systems of ODEs: damped oscillator example (simulation + phase portrait)
+#
+
+def oscillator_ode(x, t, m=1., b=1, k=1):
+ return (x[1], -k/m*x[0] - b/m*x[1])
+
+# Generate a vector plot for the damped oscillator
+mpl.figure(); mpl.clf();
+PhasePlot(oscillator_ode, [-1, 1, 10], [-1, 1, 10], 0.15);
+mpl.hold(True); mpl.plot([0], [0], '.');
+# a=gca; set(a,'FontSize',20); set(a,'DataAspectRatio',[1,1,1]);
+mpl.xlabel('x1'); mpl.ylabel('x2');
+
+# Generate a phase plot for the damped oscillator
+mpl.figure(); mpl.clf();
+mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1, 1, 1]);
+PhasePlot(oscillator_ode,
+ 'timepts', [0.25, 0.8, 2, 3], None, [
+ [-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
+ [1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
+ ], np.linspace(0, 8, 80))
+mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
+# set(gca,'DataAspectRatio',[1,1,1]);
+mpl.xlabel('x1'); mpl.ylabel('x2');
+
+mpl.show()
+
+#
+# Stability definitions
+#
+# This set of plots illustrates the various types of equilibrium points.
+#
+
+# Saddle point vector field
+def saddle_ode(x, t):
+ return (x[0] - 3*x[1], -3*x[0] + x[1]);
+
+# Asy stable
+m = 1; b = 1; k = 1; # default values
+mpl.figure(); mpl.clf();
+mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
+PhasePlot(oscillator_ode, 'timepts', [0.3, 1, 2, 3], None,
+ [[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1],
+ [1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1],
+ [-1.3,-1]],
+ np.linspace(0, 10, 100), parms = (m, b, k));
+mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
+# set(gca,'FontSize', 16);
+mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+
+# Saddle
+mpl.figure(); mpl.clf();
+mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
+PhasePlot(saddle_ode, 'timepts', [0.2, 0.5, 0.8], None,
+ [ [-1, -1], [1, 1],
+ [-1, -0.95], [-1, -0.9], [-1, -0.8], [-1, -0.6], [-1, -0.4], [-1, -0.2],
+ [-0.95, -1], [-0.9, -1], [-0.8, -1], [-0.6, -1], [-0.4, -1], [-0.2, -1],
+ [1, 0.95], [1, 0.9], [1, 0.8], [1, 0.6], [1, 0.4], [1, 0.2],
+ [0.95, 1], [0.9, 1], [0.8, 1], [0.6, 1], [0.4, 1], [0.2, 1],
+ [-0.5, -0.45], [-0.45, -0.5], [0.5, 0.45], [0.45, 0.5],
+ [-0.04, 0.04], [0.04, -0.04] ], np.linspace(0, 2, 20));
+mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
+# set(gca,'FontSize', 16);
+mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+
+# Stable isL
+m = 1; b = 0; k = 1; # zero damping
+mpl.figure(); mpl.clf();
+mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
+PhasePlot(oscillator_ode, 'timepts',
+ [pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6, 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi], None,
+ [ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ],
+ np.linspace(0, 20, 200), parms = (m, b, k));
+mpl.hold(True); mpl.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3);
+# set(gca,'FontSize', 16);
+mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+
+mpl.show()
Added: trunk/src/phaseplot.py
===================================================================
--- trunk/src/phaseplot.py (rev 0)
+++ trunk/src/phaseplot.py 2011-07-25 05:00:41 UTC (rev 169)
@@ -0,0 +1,265 @@
+# phaseplot.py - generate 2D phase portraits
+#
+# Author: Richard M. Murray
+# Date: Fall 2002 (MATLAB version), based on a version by Kristi Morgansen
+#
+# 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. The name of the author may not be used to endorse or promote products
+# derived from this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 THE AUTHOR 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.
+
+import numpy as np
+import matplotlib.pyplot as mpl
+from matplotlib.mlab import frange, find
+from exception import ControlNotImplemented
+from scipy.integrate import odeint
+
+def PhasePlot(odefun, xlimv, ylimv, scale=1, xinit=None, T=None, parms=(),
+ verbose=True):
+ """
+ PhasePlot Phase plot for 2D dynamical systems
+
+ PhasePlot(F, X1range, X2range, scale) produces a quiver plot for
+ the function F. X1range and X2range have the form [min, max, num]
+ and specify the axis limits for the plot, along with the number of
+ subdivisions in each axis, used to produce an quiver plot for the
+ vector field. The vector field is scaled by a factor 'scale'
+ (default = 1).
+
+ The function F should be the same for as used for scipy.integrate.
+ Namely, it should be a function of the form dxdt = F(x, t) that
+ accepts a state x of dimension 2 and returns a derivative dxdt of
+ dimension 2.
+
+ PhasePlot(F, X1range, X2range, scale, Xinit) produces a phase
+ plot for the function F, consisting of the quiver plot plus stream
+ lines. The streamlines start at the initial conditions listed in
+ Xinit, which should be a matrix whose rows give the desired inital
+ conditions for x1 and x2. X1range or X2range is 'auto', the arrows
+ are produced based on the stream lines. If 'scale' is negative,
+ dots are placed at the base of the arrows. If 'scale' is zero, no
+ dots are produced.
+
+ PhasePlot(F, X1range, X2range, scale, boxgrid(X1range2, X2range2))
+ produces a phase plot with stream lines generated at the edges of
+ the rectangle defined by X1range2, X2range2. These ranges are in
+ the same form as X1range, X2range.
+
+ PhasePlot(F, X1range, X2range, scale, Xinit, T) produces a phase
+ plot where the streamlines are simluted for time T (default = 50).
+
+ PhasePlot(F, X1range, X2range, scale, Xinit, T, P1, P2, ...)
+ passes additional parameters to the function F, in the same way as
+ ODE45.
+
+ Instead of drawing arrows on a grid, arrows can also be drawn on
+ streamlines by usinge the X1range and X2range arguments as follows:
+
+ X1range X2range
+ ------- -------
+ 'auto' N Draw N arrows using equally space time points
+ 'logtime' [N, lam] Draw N arrows using exponential time constant lam
+ 'timepts' [t1, t2, ...] Draw arrows at the list times
+ """
+
+ #
+ # Parameters defining the plot
+ #
+ # The constants below define the parameters that control how the
+ # plot looks. These can be modified to customize the look of the
+ # phase plot.
+ #
+
+ #! TODO: convert to keywords
+ #! TODO: eliminate old parameters that aren't used
+ # PP color = ['m', 'c', 'r', 'g', 'b', 'k', 'y'];
+ PP_stream_color = ('b'); # color array for streamlines
+ PP_stream_linewidth = 1; # line width for stream lines
+
+ PP_arrow_linewidth = 1; # line width for arrows (quiver)
+ PP_arrow_markersize = 10; # size of arrow base marker
+
+ #
+ # Figure out ranges for phase plot (argument processing)
+ #
+ auto = 0; logtime = 0; timepts = 0; Narrows = 0;
+ if (isinstance(xlimv, str) and xlimv == 'auto'):
+ auto = 1;
+ Narrows = ylimv;
+ if (verbose):
+ print 'Using auto arrows\n';
+
+ elif (isinstance(xlimv, str) and xlimv == 'logtime'):
+ logtime = 1;
+ Narrows = ylimv[0];
+ timefactor = ylimv[1];
+ if (verbose):
+ print 'Using logtime arrows\n';
+
+ elif (isinstance(xlimv, str) and xlimv == 'timepts'):
+ timepts = 1;
+ Narrows = len(ylimv);
+
+ else:
+ # Figure out the set of points for the quiver plot
+ (x1, x2) = np.meshgrid(
+ frange(xlimv[0], xlimv[1], float(xlimv[1]-xlimv[0])/xlimv[2]),
+ frange(ylimv[0], ylimv[1], float(ylimv[1]-ylimv[0])/ylimv[2]));
+
+ if ((not auto) and (not logtime) and (not timepts)):
+ # Now calculate the vector field at those points
+ (nr,nc) = x1.shape;
+ dx = np.empty((nr, nc, 2))
+ for i in range(nr):
+ for j in range(nc):
+ dx[i, j, :] = np.squeeze(odefun((x1[i,j], x2[i,j]), 0, *parms))
+
+ # Plot the quiver plot
+ #! TODO: figure out arguments to make arrows show up correctly
+ if (scale == None):
+ mpl.quiver(x1, x2, dx[:,:,1], dx[:,:,2], angles='xy')
+ elif (scale != 0):
+ #! TODO: optimize parameters for arrows
+ #! TODO: figure out arguments to make arrows show up correctly
+ xy = mpl.quiver(x1, x2, dx[:,:,0]*np.abs(scale),
+ dx[:,:,1]*np.abs(scale), angles='xy')
+ # set(xy, 'LineWidth', PP_arrow_linewidth, 'Color', 'b');
+
+ #! TODO: Tweak the shape of the plot
+ # a=gca; set(a,'DataAspectRatio',[1,1,1]);
+ # set(a,'XLim',xlimv(1:2)); set(a,'YLim',ylimv(1:2));
+ mpl.xlabel('x1'); mpl.ylabel('x2');
+
+ # See if we should also generate the streamlines
+ if (xinit == None or len(xinit) == 0):
+ return
+
+ # Convert initial conditions to a numpy array
+ xinit = np.array(xinit);
+ (nr, nc) = np.shape(xinit);
+
+ # Generate some empty matrices to keep arrow information
+ x1 = np.empty((nr, Narrows)); x2 = np.empty((nr, Narrows));
+ dx = np.empty((nr, Narrows, 2))
+
+ # See if we were passed a simulation time
+ if (T == None):
+ T = 50
+
+ # Parse the time we were passed
+ TSPAN = T;
+ if (isinstance(T, (int, float))):
+ TSPAN = np.linspace(0, T, 100);
+
+ # Figure out the limits for the plot
+ if (scale == None):
+ # Assume that the current axis are set as we want them
+ alim = mpl.axis();
+ xmin = alim[0]; xmax = alim[1];
+ ymin = alim[2]; ymax = alim[3];
+ else:
+ # Use the maximum extent of all trajectories
+ xmin = np.min(xinit[:,0]); xmax = np.max(xinit[:,0]);
+ ymin = np.min(xinit[:,1]); ymax = np.max(xinit[:,1]);
+
+ # Generate the streamlines for each initial condition
+ for i in range(nr):
+ state = odeint(odefun, xinit[i], TSPAN, args=parms);
+ time = TSPAN
+ mpl.hold(True);
+ mpl.plot(state[:,0], state[:,1])
+ #! TODO: add back in colors for stream lines
+ # PP_stream_color(np.mod(i-1, len(PP_stream_color))+1));
+ # set(h[i], 'LineWidth', PP_stream_linewidth);
+
+ # Plot arrows if quiver parameters were 'auto'
+ if (auto or logtime or timepts):
+ # Compute the locations of the arrows
+ #! TODO: check this logic to make sure it works in python
+ for j in range(Narrows):
+
+ # Figure out starting index; headless arrows start at 0
+ k = -1 if scale == None else 0;
+
+ # Figure out what time index to use for the next point
+ if (auto):
+ # Use a linear scaling based on ODE time vector
+ tind = np.floor((len(time)/Narrows) * (j-k)) + k;
+ elif (logtime):
+ # Use an exponential time vector
+ # MATLAB: tind = find(time < (j-k) / lambda, 1, 'last');
+ tarr = find(time < (j-k) / timefactor);
+ tind = tarr[-1] if len(tarr) else 0;
+ elif (timepts):
+ # Use specified time points
+ # MATLAB: tind = find(time < ylimv[j], 1, 'last');
+ tarr = find(time < ylimv[j]);
+ tind = tarr[-1] if len(tarr) else 0;
+
+ # For tailless arrows, skip the first point
+ if (tind == 0 and scale == None):
+ continue;
+
+ # Figure out the arrow at this point on the curve
+ x1[i,j] = state[tind, 0];
+ x2[i,j] = state[tind, 1];
+
+ # Skip arrows outside of initial condition box
+ if (scale != None or
+ (x1[i,j] <= xmax and x1[i,j] >= xmin and
+ x2[i,j] <= ymax and x2[i,j] >= ymin)):
+ v = odefun((x1[i,j], x2[i,j]), 0, *parms)
+ dx[i, j, 0] = v[0]; dx[i, j, 1] = v[1];
+ else:
+ dx[i, j, 0] = 0; dx[i, j, 1] = 0;
+
+ # Set the plot shape before plotting arrows to avoid warping
+ # a=gca;
+ # if (scale != None):
+ # set(a,'DataAspectRatio', [1,1,1]);
+ # if (xmin != xmax and ymin != ymax):
+ # mpl.axis([xmin, xmax, ymin, ymax]);
+ # set(a, 'Box', 'on');
+
+ # Plot arrows on the streamlines
+ if (scale == None and Narrows > 0):
+ # Use a tailless arrow
+ #! TODO: figure out arguments to make arrows show up correctly
+ mpl.quiver(x1, x2, dx[:,:,0], dx[:,:,1], angles='xy')
+ elif (scale != 0 and Narrows > 0):
+ #! TODO: figure out arguments to make arrows show up correctly
+ xy = mpl.quiver(x1, x2, dx[:,:,0]*abs(scale), dx[:,:,1]*abs(scale),
+ angles='xy')
+ # set(xy, 'LineWidth', PP_arrow_linewidth);
+ # set(xy, 'AutoScale', 'off');
+ # set(xy, 'AutoScaleFactor', 0);
+
+ if (scale < 0):
+ bp = mpl.plot(x1, x2, 'b.'); # add dots at base
+ # set(bp, 'MarkerSize', PP_arrow_markersize);
+
+ return;
Added: trunk/tests/margin_test.py
===================================================================
--- trunk/tests/margin_test.py (rev 0)
+++ trunk/tests/margin_test.py 2011-07-25 05:00:41 UTC (rev 169)
@@ -0,0 +1,47 @@
+#!/usr/bin/env python
+#
+# margin_test.py - test suit for stability margin commands
+# RMM, 15 Jul 2011
+
+import unittest
+import numpy as np
+from control.xferfcn import TransferFunction
+from control.statesp import StateSpace
+from control.margin import *
+
+class TestMargin(unittest.TestCase):
+ """These are tests for the margin commands in margin.py."""
+
+ def setUp(self):
+ self.sys1 = TransferFunction([1, 2], [1, 2, 3])
+ self.sys2 = TransferFunction([1], [1, 2, 3, 4])
+ self.sys3 = StateSpace([[1., 4.], [3., 2.]], [[1.], [-4.]],
+ [[1., 0.]], [[0.]])
+
+ def testGainPhaseMargin(self):
+ gm, pm, sm, wg, wp, ws = StabilityMargins(self.sys1);
+ gm, pm, sm, wg, wp, ws = StabilityMargins(self.sys2);
+ gm, pm, sm, wg, wp, ws = StabilityMargins(self.sys3);
+
+ def testPhaseCrossoverFrequencies(self):
+ omega, gain = PhaseCrossoverFrequencies(self.sys2)
+ np.testing.assert_array_almost_equal(omega, [1.73205, 0.])
+ np.testing.assert_array_almost_equal(gain, [-0.5, 0.25])
+
+ tf = TransferFunction([1],[1,1])
+ omega, gain = PhaseCrossoverFrequencies(tf)
+ np.testing.assert_array_almost_equal(omega, [0.])
+ np.testing.assert_array_almost_equal(gain, [1.])
+
+ # testing MIMO, only (0,0) element is considered
+ tf = TransferFunction([[[1],[2]],[[3],[4]]],
+ [[[1, 2, 3, 4],[1,1]],[[1,1],[1,1]]])
+ omega, gain = PhaseCrossoverFrequencies(tf)
+ np.testing.assert_array_almost_equal(omega, [1.73205081, 0.])
+ np.testing.assert_array_almost_equal(gain, [-0.5, 0.25])
+
+def suite():
+ return unittest.TestLoader().loadTestsFromTestCase(TestMargin)
+
+if __name__ == "__main__":
+ unittest.main()
Added: trunk/tests/phaseplot_test.py
===================================================================
--- trunk/tests/phaseplot_test.py (rev 0)
+++ trunk/tests/phaseplot_test.py 2011-07-25 05:00:41 UTC (rev 169)
@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+#
+# phaseplot_test.py - test phase plot functions
+# RMM, 17 24 2011 (based on TestMatlab from v0.4c)
+#
+# This test suite calls various phaseplot functions. Since the plots
+# themselves can't be verified, this is mainly here to make sure all
+# of the function arguments are handled correctly. If you run an
+# individual test by itself and then type show(), it should pop open
+# the figures so that you can check them visually.
+
+import unittest
+import numpy as np
+import scipy as sp
+import matplotlib.pyplot as mpl
+from control.phaseplot import *
+from numpy import pi
+
+class TestPhasePlot(unittest.TestCase):
+ def setUp(self):
+ pass;
+
+ def testInvPendNoSims(self):
+ PhasePlot(self.invpend_ode, (-6,6,10), (-6,6,10));
+
+ def testInvPendSims(self):
+ PhasePlot(self.invpend_ode, (-6,6,10), (-6,6,10),
+ xinit = ([1,1], [-1,1]));
+
+ def testInvPendTimePoints(self):
+ PhasePlot(self.invpend_ode, (-6,6,10), (-6,6,10),
+ xinit = ([1,1], [-1,1]), T=np.linspace(0,5,100));
+
+ def testInvPendLogtime(self):
+ PhasePlot(self.invpend_ode,
+ 'logtime', (3, 0.7), None,
+ [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
+ [-1, 2.1], [4.2, 2.1], [5, 2.1],
+ [2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
+ [1, -2.1], [-4.2, -2.1], [-5, -2.1] ],
+ np.linspace(0, 40, 200), verbose=False)
+
+ def testInvPendAuto(self):
+ PhasePlot(self.invpend_ode, 'auto', 0, None,
+ [[-2.3056, 2.1], [2.3056, -2.1]], 6, verbose=False)
+
+ def testOscillatorParams(self):
+ m = 1; b = 1; k = 1; # default values
+ PhasePlot(self.oscillator_ode, 'timepts', [0.3, 1, 2, 3], None,
+ [[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1],
+ [1,1], [1.3,1], [1,-1], [0.3,-1], [0,-1], [-0.25,-1],
+ [-0.5,-1], [-0.7,-1], [-1,-1], [-1.3,-1]],
+ np.linspace(0, 10, 100), parms = (m, b, k));
+
+ # Sample dynamical systems - inverted pendulum
+ def invpend_ode(self, x, t, m=1., l=1., b=0, g=9.8):
+ import numpy as np
+ return (x[1], -b/m*x[1] + (g*l/m) * np.sin(x[0]))
+
+ # Sample dynamical systems - oscillator
+ def oscillator_ode(self, x, t, m=1., b=1, k=1, extra=None):
+ return (x[1], -k/m*x[0] - b/m*x[1])
+
+def suite():
+ return unittest.TestLoader().loadTestsFromTestCase(TestPhasePlot)
+
+if __name__ == '__main__':
+ unittest.main()
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2011-07-26 02:41:10
|
Revision: 170
http://python-control.svn.sourceforge.net/python-control/?rev=170&view=rev
Author: murrayrm
Date: 2011-07-26 02:41:03 +0000 (Tue, 26 Jul 2011)
Log Message:
-----------
* Updated phaseplot module with new (more pythonic) call signature
* Added boxplot function to generate initial conditions around a box
* New genetic switch example (from FBS), demonstrating phase plots
Modified Paths:
--------------
trunk/ChangeLog
trunk/examples/phaseplots.py
trunk/src/__init__.py
trunk/src/phaseplot.py
Added Paths:
-----------
trunk/examples/genswitch.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2011-07-25 05:00:41 UTC (rev 169)
+++ trunk/ChangeLog 2011-07-26 02:41:03 UTC (rev 170)
@@ -1,3 +1,17 @@
+2011-07-25 Richard Murray <murray@malabar.local>
+
+ * examples/phaseplots.py: updated calls to PhasePlot to use new
+ argument structure
+
+ * src/phaseplot.py (PhasePlot): Updated call signature to be
+ more pythonic and fixed up documentation.
+
+ * examples/genswitch.py (genswitch): added new example showing
+ PhasePlot functionality
+
+ * src/phaseplot.py (boxgrid): added function to compute initial
+ conditions around the edges of a box
+
2011-07-24 Richard Murray <murray@malabar.local>
* tests/margin_test.py: added simple unit tests for margin functions
Added: trunk/examples/genswitch.py
===================================================================
--- trunk/examples/genswitch.py (rev 0)
+++ trunk/examples/genswitch.py 2011-07-26 02:41:03 UTC (rev 170)
@@ -0,0 +1,81 @@
+# genswitch_plot.m - run the Collins genetic switch model
+# RMM, 24 Jan 07
+#
+# This file contains an example from FBS of a simple dynamical model
+# of a genetic switch. Plots time traces and a phase portrait using
+# the python-control library.
+
+import numpy as np
+import matplotlib.pyplot as mpl
+from scipy.integrate import odeint
+from matplotlib.mlab import frange
+from control import PhasePlot, boxgrid
+
+# Simple model of a genetic switch
+#
+# This function implements the basic model of the genetic switch
+# Parameters taken from Gardner, Cantor and Collins, Nature, 2000
+def genswitch(y, t, mu=4, n=2):
+ return (mu / (1 + y[1]**n) - y[0], mu / (1 + y[0]**n) - y[1])
+
+# Run a simulation from an initial condition
+tim1 = np.linspace(0, 10, 100)
+sol1 = odeint(genswitch, [1, 5], tim1)
+
+# Extract the equlibirum points
+mu = 4; n = 2; # switch parameters
+eqpt = np.empty(3);
+eqpt[0] = sol1[0,-1]
+eqpt[1] = sol1[1,-1]
+eqpt[2] = 0; # fzero(@(x) mu/(1+x^2) - x, 2);
+
+# Run another simulation showing switching behavior
+tim2 = np.linspace(11, 25, 100);
+sol2 = odeint(genswitch, sol1[-1,:] + [2, -2], tim2)
+
+# First plot out the curves that define the equilibria
+u = frange(0, 4.5, 0.1)
+f = np.divide(mu, (1 + u**n)) # mu / (1 + u^n), elementwise
+
+mpl.figure(1); mpl.clf();
+mpl.axis([0, 5, 0, 5]); # box on;
+mpl.plot(u, f, '-', f, u, '--') # 'LineWidth', AM_data_linewidth);
+mpl.legend('z1, f(z1)', 'z2, f(z2)') # legend(lgh, 'boxoff');
+mpl.plot([0, 3], [0, 3], 'k-') # 'LineWidth', AM_ref_linewidth);
+mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.',
+ eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3);
+mpl.xlabel('z1, f(z2)');
+mpl.ylabel('z2, f(z1)');
+
+# Time traces
+mpl.figure(3); mpl.clf(); # subplot(221);
+mpl.plot(tim1, sol1[:,0], 'b-', tim1, sol1[:,1], 'g--');
+# set(pl, 'LineWidth', AM_data_linewidth);
+mpl.hold(True);
+mpl.plot([tim1[-1], tim1[-1]+1],
+ [sol1[-1,0], sol2[0,1]], 'ko:',
+ [tim1[-1], tim1[-1]+1], [sol1[-1,1], sol2[0,0]], 'ko:');
+# set(pl, 'LineWidth', AM_data_linewidth, 'MarkerSize', AM_data_markersize);
+mpl.plot(tim2, sol2[:,0], 'b-', tim2, sol2[:,1], 'g--');
+# set(pl, 'LineWidth', AM_data_linewidth);
+mpl.axis([0, 25, 0, 5]);
+
+mpl.xlabel('Time {\itt} [scaled]');
+mpl.ylabel('Protein concentrations [scaled]');
+mpl.legend('z1 (A)', 'z2 (B)') # 'Orientation', 'horizontal');
+# legend(legh, 'boxoff');
+
+# Phase portrait
+mpl.figure(2); mpl.clf(); # subplot(221);
+mpl.axis([0, 5, 0, 5]); # set(gca, 'DataAspectRatio', [1, 1, 1]);
+PhasePlot(genswitch, X0 = boxgrid([0, 5, 6], [0, 5, 6]), T = 10,
+ timepts = [0.2, 0.6, 1.2])
+
+# Add the stable equilibrium points
+mpl.hold(True);
+mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.',
+ eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3);
+
+mpl.xlabel('Protein A [scaled]');
+mpl.ylabel('Protein B [scaled]'); # 'Rotation', 90);
+
Modified: trunk/examples/phaseplots.py
===================================================================
--- trunk/examples/phaseplots.py 2011-07-25 05:00:41 UTC (rev 169)
+++ trunk/examples/phaseplots.py 2011-07-26 02:41:03 UTC (rev 170)
@@ -27,16 +27,16 @@
# Outer trajectories
PhasePlot(invpend_ode,
- 'logtime', (3, 0.7), None,
- [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
- [-1, 2.1], [4.2, 2.1], [5, 2.1],
- [2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
- [1, -2.1], [-4.2, -2.1], [-5, -2.1] ],
- np.linspace(0, 40, 200))
+ X0 = [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
+ [-1, 2.1], [4.2, 2.1], [5, 2.1],
+ [2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
+ [1, -2.1], [-4.2, -2.1], [-5, -2.1] ],
+ T = np.linspace(0, 40, 200),
+ logtime = (3, 0.7) )
# Separatrices
mpl.hold(True);
-PhasePlot(invpend_ode, 'auto', 0, None, [[-2.3056, 2.1], [2.3056, -2.1]], 6)
+PhasePlot(invpend_ode, X0 = [[-2.3056, 2.1], [2.3056, -2.1]], T=6, lingrid=0)
mpl.show();
#
@@ -57,10 +57,10 @@
mpl.figure(); mpl.clf();
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1, 1, 1]);
PhasePlot(oscillator_ode,
- 'timepts', [0.25, 0.8, 2, 3], None, [
- [-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
- [1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
- ], np.linspace(0, 8, 80))
+ X0 = [
+ [-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
+ [1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
+ ], T = np.linspace(0, 8, 80), timepts = [0.25, 0.8, 2, 3])
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
# set(gca,'DataAspectRatio',[1,1,1]);
mpl.xlabel('x1'); mpl.ylabel('x2');
@@ -81,11 +81,13 @@
m = 1; b = 1; k = 1; # default values
mpl.figure(); mpl.clf();
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
-PhasePlot(oscillator_ode, 'timepts', [0.3, 1, 2, 3], None,
- [[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1],
- [1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1],
- [-1.3,-1]],
- np.linspace(0, 10, 100), parms = (m, b, k));
+PhasePlot(oscillator_ode,
+ X0 = [
+ [-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1],
+ [1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1],
+ [-1.3,-1]
+ ], T = np.linspace(0, 10, 100),
+ timepts = [0.3, 1, 2, 3], parms = (m, b, k));
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
# set(gca,'FontSize', 16);
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
@@ -93,14 +95,14 @@
# Saddle
mpl.figure(); mpl.clf();
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
-PhasePlot(saddle_ode, 'timepts', [0.2, 0.5, 0.8], None,
+PhasePlot(saddle_ode, scale = 2, timepts = [0.2, 0.5, 0.8], X0 =
[ [-1, -1], [1, 1],
[-1, -0.95], [-1, -0.9], [-1, -0.8], [-1, -0.6], [-1, -0.4], [-1, -0.2],
[-0.95, -1], [-0.9, -1], [-0.8, -1], [-0.6, -1], [-0.4, -1], [-0.2, -1],
[1, 0.95], [1, 0.9], [1, 0.8], [1, 0.6], [1, 0.4], [1, 0.2],
[0.95, 1], [0.9, 1], [0.8, 1], [0.6, 1], [0.4, 1], [0.2, 1],
[-0.5, -0.45], [-0.45, -0.5], [0.5, 0.45], [0.45, 0.5],
- [-0.04, 0.04], [0.04, -0.04] ], np.linspace(0, 2, 20));
+ [-0.04, 0.04], [0.04, -0.04] ], T = np.linspace(0, 2, 20));
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
# set(gca,'FontSize', 16);
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
@@ -109,10 +111,10 @@
m = 1; b = 0; k = 1; # zero damping
mpl.figure(); mpl.clf();
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
-PhasePlot(oscillator_ode, 'timepts',
- [pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6, 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi], None,
- [ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ],
- np.linspace(0, 20, 200), parms = (m, b, k));
+PhasePlot(oscillator_ode, timepts =
+ [pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6, 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi],
+ X0 = [ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ],
+ T = np.linspace(0, 20, 200), parms = (m, b, k));
mpl.hold(True); mpl.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3);
# set(gca,'FontSize', 16);
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2011-07-25 05:00:41 UTC (rev 169)
+++ trunk/src/__init__.py 2011-07-26 02:41:03 UTC (rev 170)
@@ -64,6 +64,7 @@
from mateqn import *
from modelsimp import *
from nichols import *
+from phaseplot import PhasePlot, boxgrid
from rlocus import *
from statefbk import *
from statesp import *
Modified: trunk/src/phaseplot.py
===================================================================
--- trunk/src/phaseplot.py 2011-07-25 05:00:41 UTC (rev 169)
+++ trunk/src/phaseplot.py 2011-07-26 02:41:03 UTC (rev 170)
@@ -1,7 +1,8 @@
# phaseplot.py - generate 2D phase portraits
#
# Author: Richard M. Murray
-# Date: Fall 2002 (MATLAB version), based on a version by Kristi Morgansen
+# Date: 24 July 2011, converted from MATLAB version (2002); based on
+# a version by Kristi Morgansen
#
# Copyright (c) 2011 by California Institute of Technology
# All rights reserved.
@@ -38,99 +39,107 @@
from exception import ControlNotImplemented
from scipy.integrate import odeint
-def PhasePlot(odefun, xlimv, ylimv, scale=1, xinit=None, T=None, parms=(),
- verbose=True):
+def PhasePlot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
+ lingrid=None, lintime=None, logtime=None, timepts=None,
+ parms=(), verbose=True):
"""
- PhasePlot Phase plot for 2D dynamical systems
+ Phase plot for 2D dynamical systems
- PhasePlot(F, X1range, X2range, scale) produces a quiver plot for
- the function F. X1range and X2range have the form [min, max, num]
- and specify the axis limits for the plot, along with the number of
- subdivisions in each axis, used to produce an quiver plot for the
- vector field. The vector field is scaled by a factor 'scale'
- (default = 1).
+ Produces a vector field or stream line plot for a planar system.
- The function F should be the same for as used for scipy.integrate.
- Namely, it should be a function of the form dxdt = F(x, t) that
- accepts a state x of dimension 2 and returns a derivative dxdt of
- dimension 2.
+ Call signatures:
+ PhasePlot(func, X, Y, ...) - display vector field on meshgrid
+ PhasePlot(func, X, Y, scale, ...) - scale arrows
+ PhasePlot(func. X0=(...), T=Tmax, ...) - display stream lines
+ PhasePlot(func, X, Y, X0=[...], T=Tmax, ...) - plot both
+ PhasePlot(func, X0=[...], T=Tmax, lingrid=N, ...) - plot both
+ PhasePlot(func, X0=[...], lintime=N, ...) - stream lines with arrows
- PhasePlot(F, X1range, X2range, scale, Xinit) produces a phase
- plot for the function F, consisting of the quiver plot plus stream
- lines. The streamlines start at the initial conditions listed in
- Xinit, which should be a matrix whose rows give the desired inital
- conditions for x1 and x2. X1range or X2range is 'auto', the arrows
- are produced based on the stream lines. If 'scale' is negative,
- dots are placed at the base of the arrows. If 'scale' is zero, no
- dots are produced.
-
- PhasePlot(F, X1range, X2range, scale, boxgrid(X1range2, X2range2))
- produces a phase plot with stream lines generated at the edges of
- the rectangle defined by X1range2, X2range2. These ranges are in
- the same form as X1range, X2range.
+ Parameters
+ ----------
+ func : callable(x, t, ...)
+ Computes the time derivative of y (compatible with odeint).
+ The function should be the same for as used for
+ scipy.integrate. Namely, it should be a function of the form
+ dxdt = F(x, t) that accepts a state x of dimension 2 and
+ returns a derivative dx/dt of dimension 2.
- PhasePlot(F, X1range, X2range, scale, Xinit, T) produces a phase
- plot where the streamlines are simluted for time T (default = 50).
+ X, Y: ndarray, optional
+ Two 1-D arrays representing x and y coordinates of a grid.
+ These arguments are passed to meshgrid and generate the lists
+ of points at which the vector field is plotted. If absent (or
+ None), the vector field is not plotted.
- PhasePlot(F, X1range, X2range, scale, Xinit, T, P1, P2, ...)
- passes additional parameters to the function F, in the same way as
- ODE45.
+ scale: float, optional
+ Scale size of arrows; default = 1
- Instead of drawing arrows on a grid, arrows can also be drawn on
- streamlines by usinge the X1range and X2range arguments as follows:
+ X0: ndarray of initial conditions, optional
+ List of initial conditions from which streamlines are plotted.
+ Each initial condition should be a pair of numbers.
- X1range X2range
- ------- -------
- 'auto' N Draw N arrows using equally space time points
- 'logtime' [N, lam] Draw N arrows using exponential time constant lam
- 'timepts' [t1, t2, ...] Draw arrows at the list times
- """
+ T: array-like or number, optional
+ Length of time to run simulations that generate streamlines.
+ If a single number, the same simulation time is used for all
+ initial conditions. Otherwise, should be a list of length
+ len(X0) that gives the simulation time for each initial
+ condition. Default value = 50.
- #
- # Parameters defining the plot
- #
- # The constants below define the parameters that control how the
- # plot looks. These can be modified to customize the look of the
- # phase plot.
- #
+ lingrid = N or (N, M): integer or 2-tuple of integers, optional
+ If X0 is given and X, Y are missing, a grid of arrows is
+ produced using the limits of the initial conditions, with N
+ grid points in each dimension or N grid points in x and M grid
+ points in y.
- #! TODO: convert to keywords
- #! TODO: eliminate old parameters that aren't used
- # PP color = ['m', 'c', 'r', 'g', 'b', 'k', 'y'];
- PP_stream_color = ('b'); # color array for streamlines
- PP_stream_linewidth = 1; # line width for stream lines
+ lintime = N: integer, optional
+ Draw N arrows using equally space time points
- PP_arrow_linewidth = 1; # line width for arrows (quiver)
- PP_arrow_markersize = 10; # size of arrow base marker
+ logtime = (N, lambda): (integer, float), optional
+ Draw N arrows using exponential time constant lambda
+ timepts = [t1, t2, ...]: array-like, optional
+ Draw arrows at the given list times
+
+ parms: tuple, optional
+ List of parameters to pass to vector field: func(x, t, *parms)
+
+ See also
+ --------
+ boxgrid(X, Y): construct box-shaped grid of initial conditions
+
+ Examples
+ --------
+ """
+
#
# Figure out ranges for phase plot (argument processing)
#
- auto = 0; logtime = 0; timepts = 0; Narrows = 0;
- if (isinstance(xlimv, str) and xlimv == 'auto'):
- auto = 1;
- Narrows = ylimv;
+ #! TODO: need to add error checking to arguments
+ autoFlag = False; logtimeFlag = False; timeptsFlag = False; Narrows = 0;
+ if (lingrid != None):
+ autoFlag = True;
+ Narrows = lingrid;
if (verbose):
print 'Using auto arrows\n';
- elif (isinstance(xlimv, str) and xlimv == 'logtime'):
- logtime = 1;
- Narrows = ylimv[0];
- timefactor = ylimv[1];
+ elif (logtime != None):
+ logtimeFlag = True;
+ Narrows = logtime[0];
+ timefactor = logtime[1];
if (verbose):
print 'Using logtime arrows\n';
- elif (isinstance(xlimv, str) and xlimv == 'timepts'):
- timepts = 1;
- Narrows = len(ylimv);
+ elif (timepts != None):
+ timeptsFlag = True;
+ Narrows = len(timepts);
else:
# Figure out the set of points for the quiver plot
+ #! TODO: Add sanity checks
(x1, x2) = np.meshgrid(
- frange(xlimv[0], xlimv[1], float(xlimv[1]-xlimv[0])/xlimv[2]),
- frange(ylimv[0], ylimv[1], float(ylimv[1]-ylimv[0])/ylimv[2]));
+ frange(X[0], X[1], float(X[1]-X[0])/X[2]),
+ frange(Y[0], Y[1], float(Y[1]-Y[0])/Y[2]));
- if ((not auto) and (not logtime) and (not timepts)):
+ if ((not autoFlag) and (not logtimeFlag) and (not timeptsFlag)):
# Now calculate the vector field at those points
(nr,nc) = x1.shape;
dx = np.empty((nr, nc, 2))
@@ -151,16 +160,16 @@
#! TODO: Tweak the shape of the plot
# a=gca; set(a,'DataAspectRatio',[1,1,1]);
- # set(a,'XLim',xlimv(1:2)); set(a,'YLim',ylimv(1:2));
+ # set(a,'XLim',X(1:2)); set(a,'YLim',Y(1:2));
mpl.xlabel('x1'); mpl.ylabel('x2');
# See if we should also generate the streamlines
- if (xinit == None or len(xinit) == 0):
+ if (X0 == None or len(X0) == 0):
return
# Convert initial conditions to a numpy array
- xinit = np.array(xinit);
- (nr, nc) = np.shape(xinit);
+ X0 = np.array(X0);
+ (nr, nc) = np.shape(X0);
# Generate some empty matrices to keep arrow information
x1 = np.empty((nr, Narrows)); x2 = np.empty((nr, Narrows));
@@ -183,12 +192,12 @@
ymin = alim[2]; ymax = alim[3];
else:
# Use the maximum extent of all trajectories
- xmin = np.min(xinit[:,0]); xmax = np.max(xinit[:,0]);
- ymin = np.min(xinit[:,1]); ymax = np.max(xinit[:,1]);
+ xmin = np.min(X0[:,0]); xmax = np.max(X0[:,0]);
+ ymin = np.min(X0[:,1]); ymax = np.max(X0[:,1]);
# Generate the streamlines for each initial condition
for i in range(nr):
- state = odeint(odefun, xinit[i], TSPAN, args=parms);
+ state = odeint(odefun, X0[i], TSPAN, args=parms);
time = TSPAN
mpl.hold(True);
mpl.plot(state[:,0], state[:,1])
@@ -197,7 +206,7 @@
# set(h[i], 'LineWidth', PP_stream_linewidth);
# Plot arrows if quiver parameters were 'auto'
- if (auto or logtime or timepts):
+ if (autoFlag or logtimeFlag or timeptsFlag):
# Compute the locations of the arrows
#! TODO: check this logic to make sure it works in python
for j in range(Narrows):
@@ -206,18 +215,18 @@
k = -1 if scale == None else 0;
# Figure out what time index to use for the next point
- if (auto):
+ if (autoFlag):
# Use a linear scaling based on ODE time vector
tind = np.floor((len(time)/Narrows) * (j-k)) + k;
- elif (logtime):
+ elif (logtimeFlag):
# Use an exponential time vector
# MATLAB: tind = find(time < (j-k) / lambda, 1, 'last');
tarr = find(time < (j-k) / timefactor);
tind = tarr[-1] if len(tarr) else 0;
- elif (timepts):
+ elif (timeptsFlag):
# Use specified time points
- # MATLAB: tind = find(time < ylimv[j], 1, 'last');
- tarr = find(time < ylimv[j]);
+ # MATLAB: tind = find(time < Y[j], 1, 'last');
+ tarr = find(time < timepts[j]);
tind = tarr[-1] if len(tarr) else 0;
# For tailless arrows, skip the first point
@@ -263,3 +272,20 @@
# set(bp, 'MarkerSize', PP_arrow_markersize);
return;
+
+# Utility function for generating initial conditions around a box
+def boxgrid(xlimp, ylimp):
+ """BOXGRID generate list of points on edge of box
+
+ list = BOXGRID([xmin xmax xnum], [ymin ymax ynum]) generates a
+ list of points that correspond to a uniform grid at the end of the
+ box defined by the corners [xmin ymin] and [xmax ymax].
+ """
+
+ sx10 = frange(xlimp[0], xlimp[1], float(xlimp[1]-xlimp[0])/xlimp[2])
+ sy10 = frange(ylimp[0], ylimp[1], float(ylimp[1]-ylimp[0])/ylimp[2])
+
+ sx1 = np.hstack((0, sx10, 0*sy10+sx10[0], sx10, 0*sy10+sx10[-1]))
+ sx2 = np.hstack((0, 0*sx10+sy10[0], sy10, 0*sx10+sy10[-1], sy10))
+
+ return np.transpose( np.vstack((sx1, sx2)) )
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2011-08-07 15:13:05
|
Revision: 172
http://python-control.svn.sourceforge.net/python-control/?rev=172&view=rev
Author: murrayrm
Date: 2011-08-07 15:12:58 +0000 (Sun, 07 Aug 2011)
Log Message:
-----------
* Changed function names to be more consistent with python coding
standards (PEP 8): StepReponse is now step_response, etc. No
changes in control.matlab interface.
* Changed version number (for next release) to 0.5a due to changes in
call signatures.
Modified Paths:
--------------
trunk/ChangeLog
trunk/examples/genswitch.py
trunk/examples/phaseplots.py
trunk/examples/pvtol-lqr.py
trunk/examples/pvtol-nested-ss.py
trunk/examples/pvtol-nested.py
trunk/examples/secord-matlab.py
trunk/examples/test-response.py
trunk/setup.py
trunk/src/__init__.py
trunk/src/freqplot.py
trunk/src/margins.py
trunk/src/matlab.py
trunk/src/nichols.py
trunk/src/phaseplot.py
trunk/src/rlocus.py
trunk/src/timeresp.py
trunk/tests/margin_test.py
trunk/tests/phaseplot_test.py
trunk/tests/test_all.py
trunk/tests/timeresp_test.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/ChangeLog 2011-08-07 15:12:58 UTC (rev 172)
@@ -1,3 +1,37 @@
+2011-08-07 Richard Murray <murray@malabar.local>
+
+ * examples/secord-matlab.py, examples/pvtol-nested-ss.py,
+ examples/test-response.py, examples/pvtol-nested.py: fixed small bug
+ in order of output arguments for step command
+
+2011-08-06 Richard Murray <murray@malabar.local>
+
+ * src/matlab.py (ngrid): copy documentation from nichols_grid
+
+ * src/__init__.py: changed import commands to import specific
+ functions rather than '*' (better modularity)
+
+ * src/freqplot.py: default function names are now bode_plot,
+ nyquist_plot and gangof4_plot (still with aliases to non-"_plot"
+ versions)
+
+ * src/nichols.py (nichols_plot): updated nichols to nichols_plot for
+ consistency with other python-control plotting functions. Set up
+ alias for original name
+
+ * src/margins.py: StabilityMargins, PhaseCrossoverFrequencies ->
+ stability_margins, phase_crossover_frequencies
+
+ * src/phaseplot.py: changed PhasePlot and boxgrid to phase_plot,
+ box_grid
+
+ * src/timeresp.py: changed ForcedReponse, InitialReponse,
+ ImpulseReponse and StepResponse to forced_response,
+ initial_response, impulse_response and step_response.
+
+ * src/rlocus.py: changed RootLocus to root_locus for better
+ compatability with PEP 8. Also updated unit tests and examples.
+
2011-07-25 Richard Murray <murray@malabar.local>
* tests/phaseplot_test.py: updated unit tests to use new call
Modified: trunk/examples/genswitch.py
===================================================================
--- trunk/examples/genswitch.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/examples/genswitch.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -9,7 +9,7 @@
import matplotlib.pyplot as mpl
from scipy.integrate import odeint
from matplotlib.mlab import frange
-from control import PhasePlot, boxgrid
+from control import phase_plot, box_grid
# Simple model of a genetic switch
#
@@ -68,7 +68,7 @@
# Phase portrait
mpl.figure(2); mpl.clf(); # subplot(221);
mpl.axis([0, 5, 0, 5]); # set(gca, 'DataAspectRatio', [1, 1, 1]);
-PhasePlot(genswitch, X0 = boxgrid([0, 5, 6], [0, 5, 6]), T = 10,
+phase_plot(genswitch, X0 = box_grid([0, 5, 6], [0, 5, 6]), T = 10,
timepts = [0.2, 0.6, 1.2])
# Add the stable equilibrium points
Modified: trunk/examples/phaseplots.py
===================================================================
--- trunk/examples/phaseplots.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/examples/phaseplots.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -6,7 +6,7 @@
import numpy as np
import matplotlib.pyplot as mpl
-from control.phaseplot import PhasePlot
+from control.phaseplot import phase_plot
from numpy import pi
# Clear out any figures that are present
@@ -26,7 +26,7 @@
mpl.title('Inverted pendlum')
# Outer trajectories
-PhasePlot(invpend_ode,
+phase_plot(invpend_ode,
X0 = [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
[-1, 2.1], [4.2, 2.1], [5, 2.1],
[2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
@@ -36,7 +36,7 @@
# Separatrices
mpl.hold(True);
-PhasePlot(invpend_ode, X0 = [[-2.3056, 2.1], [2.3056, -2.1]], T=6, lingrid=0)
+phase_plot(invpend_ode, X0 = [[-2.3056, 2.1], [2.3056, -2.1]], T=6, lingrid=0)
mpl.show();
#
@@ -48,7 +48,7 @@
# Generate a vector plot for the damped oscillator
mpl.figure(); mpl.clf();
-PhasePlot(oscillator_ode, [-1, 1, 10], [-1, 1, 10], 0.15);
+phase_plot(oscillator_ode, [-1, 1, 10], [-1, 1, 10], 0.15);
mpl.hold(True); mpl.plot([0], [0], '.');
# a=gca; set(a,'FontSize',20); set(a,'DataAspectRatio',[1,1,1]);
mpl.xlabel('x1'); mpl.ylabel('x2');
@@ -56,7 +56,7 @@
# Generate a phase plot for the damped oscillator
mpl.figure(); mpl.clf();
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1, 1, 1]);
-PhasePlot(oscillator_ode,
+phase_plot(oscillator_ode,
X0 = [
[-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
[1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
@@ -81,7 +81,7 @@
m = 1; b = 1; k = 1; # default values
mpl.figure(); mpl.clf();
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
-PhasePlot(oscillator_ode,
+phase_plot(oscillator_ode,
X0 = [
[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1],
[1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1],
@@ -95,7 +95,7 @@
# Saddle
mpl.figure(); mpl.clf();
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
-PhasePlot(saddle_ode, scale = 2, timepts = [0.2, 0.5, 0.8], X0 =
+phase_plot(saddle_ode, scale = 2, timepts = [0.2, 0.5, 0.8], X0 =
[ [-1, -1], [1, 1],
[-1, -0.95], [-1, -0.9], [-1, -0.8], [-1, -0.6], [-1, -0.4], [-1, -0.2],
[-0.95, -1], [-0.9, -1], [-0.8, -1], [-0.6, -1], [-0.4, -1], [-0.2, -1],
@@ -111,7 +111,7 @@
m = 1; b = 0; k = 1; # zero damping
mpl.figure(); mpl.clf();
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
-PhasePlot(oscillator_ode, timepts =
+phase_plot(oscillator_ode, timepts =
[pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6, 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi],
X0 = [ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ],
T = np.linspace(0, 20, 200), parms = (m, b, k));
Modified: trunk/examples/pvtol-lqr.py
===================================================================
--- trunk/examples/pvtol-lqr.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/examples/pvtol-lqr.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -111,11 +111,11 @@
# Step response for the first input
H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx);
-(Tx, Yx) = step(H1ax, T=linspace(0,10,100));
+(Yx, Tx) = step(H1ax, T=linspace(0,10,100));
# Step response for the second input
H1ay = ss(Ay - By*K1a[1,alt], By*K1a[1,alt]*yd[alt,:], Cy, Dy);
-(Ty, Yy) = step(H1ay, T=linspace(0,10,100));
+(Yy, Ty) = step(H1ay, T=linspace(0,10,100));
subplot(221); title("Identity weights")
# plot(T, Y[:,1, 1], '-', T, Y[:,2, 2], '--'); hold(True);
@@ -136,9 +136,9 @@
Qu1c = (200**2)*diag([1, 1]); (K1c, X, E) = lqr(A, B, Qx1, Qu1c);
H1cx = ss(Ax - Bx*K1c[0,lat], Bx*K1c[0,lat]*xd[lat,:],Cx, Dx);
-[T1, Y1] = step(H1ax, T=linspace(0,10,100));
-[T2, Y2] = step(H1bx, T=linspace(0,10,100));
-[T3, Y3] = step(H1cx, T=linspace(0,10,100));
+[Y1, T1] = step(H1ax, T=linspace(0,10,100));
+[Y2, T2] = step(H1bx, T=linspace(0,10,100));
+[Y3, T3] = step(H1cx, T=linspace(0,10,100));
subplot(222); title("Effect of input weights")
plot(T1.T, Y1.T, 'b-'); hold(True);
@@ -160,8 +160,8 @@
H2y = ss(Ay - By*K2[1,alt], By*K2[1,alt]*yd[alt,:], Cy, Dy);
subplot(223); title("Output weighting")
-[T2x, Y2x] = step(H2x, T=linspace(0,10,100));
-[T2y, Y2y] = step(H2y, T=linspace(0,10,100));
+[Y2x, T2x] = step(H2x, T=linspace(0,10,100));
+[Y2y, T2y] = step(H2y, T=linspace(0,10,100));
plot(T2x.T, Y2x.T, T2y.T, Y2y.T)
ylabel('position');
xlabel('time'); ylabel('position');
@@ -183,8 +183,8 @@
H3y = ss(Ay - By*K3[1,alt], By*K3[1,alt]*yd[alt,:], Cy, Dy);
subplot(224)
# step(H3x, H3y, 10);
-[T3x, Y3x] = step(H3x, T=linspace(0,10,100));
-[T3y, Y3y] = step(H3y, T=linspace(0,10,100));
+[Y3x, T3x] = step(H3x, T=linspace(0,10,100));
+[Y3y, T3y] = step(H3y, T=linspace(0,10,100));
plot(T3x.T, Y3x.T, T3y.T, Y3y.T)
title("Physically motivated weights")
xlabel('time');
Modified: trunk/examples/pvtol-nested-ss.py
===================================================================
--- trunk/examples/pvtol-nested-ss.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/examples/pvtol-nested-ss.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -144,10 +144,10 @@
# 'EdgeColor', color, 'FaceColor', color);
figure(9);
-(Tvec, Yvec) = step(T, linspace(1, 20));
+(Yvec, Tvec) = step(T, linspace(1, 20));
plot(Tvec.T, Yvec.T); hold(True);
-(Tvec, Yvec) = step(Co*S, linspace(1, 20));
+(Yvec, Tvec) = step(Co*S, linspace(1, 20));
plot(Tvec.T, Yvec.T);
#TODO: PZmap for statespace systems has not yet been implemented.
Modified: trunk/examples/pvtol-nested.py
===================================================================
--- trunk/examples/pvtol-nested.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/examples/pvtol-nested.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -134,10 +134,10 @@
# 'EdgeColor', color, 'FaceColor', color);
figure(9);
-(Tvec, Yvec) = step(T, linspace(1, 20));
+(Yvec, Tvec) = step(T, linspace(0, 20));
plot(Tvec.T, Yvec.T); hold(True);
-(Tvec, Yvec) = step(Co*S, linspace(1, 20));
+(Yvec, Tvec) = step(Co*S, linspace(0, 20));
plot(Tvec.T, Yvec.T);
figure(10); clf();
Modified: trunk/examples/secord-matlab.py
===================================================================
--- trunk/examples/secord-matlab.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/examples/secord-matlab.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -17,7 +17,7 @@
# Step response for the system
figure(1)
-T, yout = step(sys)
+yout, T = step(sys)
plot(T.T, yout.T)
# Bode plot for the system
Modified: trunk/examples/test-response.py
===================================================================
--- trunk/examples/test-response.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/examples/test-response.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -10,9 +10,9 @@
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))
+(y1a, T1a) = step(sys1)
+(y1b, T1b) = step(sys1, T = arange(0, 10, 0.1))
+(y1c, T1c) = step(sys1, X0 = [1, 0])
+(y2a, T2a) = step(sys2, T = arange(0, 10, 0.1))
plot(T1a, y1a, T1b, y1b, T1c, y1c, T2a, y2a)
Modified: trunk/setup.py
===================================================================
--- trunk/setup.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/setup.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -4,7 +4,7 @@
from setuptools import setup
setup(name = 'control',
- version = '0.4d',
+ version = '0.5a',
description = 'Python Control Systems Library',
author = 'Richard Murray',
author_email = 'mu...@cd...',
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/src/__init__.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -57,17 +57,18 @@
# Import functions from within the control system library
#! Should probably only import the exact functions we use...
-from bdalg import *
-from delay import *
-from freqplot import *
-from margins import *
-from mateqn import *
-from modelsimp import *
-from nichols import *
-from phaseplot import PhasePlot, boxgrid
-from rlocus import *
-from statefbk import *
-from statesp import *
-from timeresp import ForcedResponse, InitialResponse, StepResponse, \
- ImpulseResponse
-from xferfcn import *
+from bdalg import series, parallel, negate, feedback
+from delay import pade
+from freqplot import bode_plot, nyquist_plot, gangof4_plot
+from freqplot import bode, nyquist, gangof4
+from margins import stability_margins, phase_crossover_frequencies
+from mateqn import lyap, dlyap, care, dare
+from modelsimp import hsvd, modred, balred, era, markov
+from nichols import nichols_plot, nichols
+from phaseplot import phase_plot, box_grid
+from rlocus import root_locus
+from statefbk import place, lqr, ctrb, obsv, gram
+from statesp import StateSpace
+from timeresp import forced_response, initial_response, step_response, \
+ impulse_response
+from xferfcn import TransferFunction
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/src/freqplot.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -55,7 +55,7 @@
#
# Bode plot
-def BodePlot(syslist, omega=None, dB=False, Hz=False, deg=True,
+def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True,
color=None, Plot=True):
"""Bode plot for a system
@@ -178,7 +178,7 @@
return mags, phases, omegas
# Nyquist plot
-def NyquistPlot(syslist, omega=None, Plot=True):
+def nyquist_plot(syslist, omega=None, Plot=True):
"""Nyquist plot for a system
Plots a Nyquist plot for the system over a (optional) frequency range.
@@ -245,7 +245,7 @@
# Gang of Four
#! TODO: think about how (and whether) to handle lists of systems
-def GangOf4Plot(P, C, omega=None):
+def gangof4_plot(P, C, omega=None):
"""Plot the "Gang of 4" transfer functions for a system
Generates a 2x2 plot showing the "Gang of 4" sensitivity functions
@@ -363,6 +363,6 @@
return omega
# Function aliases
-bode = BodePlot
-nyquist = NyquistPlot
-gangof4 = GangOf4Plot
+bode = bode_plot
+nyquist = nyquist_plot
+gangof4 = gangof4_plot
Modified: trunk/src/margins.py
===================================================================
--- trunk/src/margins.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/src/margins.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -4,8 +4,8 @@
Routeins in this module:
-margin.StabilityMargins
-margin.PhaseCrossoverFrequencies
+margin.stability_margins
+margin.phase_crossover_frequencies
"""
"""Copyright (c) 2011 by California Institute of Technology
@@ -54,13 +54,13 @@
# gain and phase margins
# contributed by Sawyer B. Fuller <mi...@ca...>
#! TODO - need to add unit test functions
-def StabilityMargins(sysdata, deg=True):
+def stability_margins(sysdata, deg=True):
"""Calculate gain, phase and stability margins and associated
crossover frequencies.
Usage:
- gm, pm, sm, wg, wp, ws = StabilityMargins(sysdata, deg=True)
+ gm, pm, sm, wg, wp, ws = stability_margins(sysdata, deg=True)
Parameters
----------
@@ -149,13 +149,13 @@
# Contributed by Steffen Waldherr <wal...@is...>
#! TODO - need to add test functions
-def PhaseCrossoverFrequencies(sys):
+def phase_crossover_frequencies(sys):
"""
Compute frequencies and gains at intersections with real axis
in Nyquist plot.
Call as:
- omega, gain = PhaseCrossoverFrequencies()
+ omega, gain = phase_crossover_frequencies()
Returns
-------
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/src/matlab.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -62,6 +62,7 @@
# Libraries that we make use of
import scipy as sp # SciPy library (used all over)
import numpy as np # NumPy library
+import re # regular expressions
from copy import deepcopy
# Import MATLAB-like functions that are defined in other packages
@@ -990,15 +991,10 @@
return freqplot.bode(syslist, omega, **keywords)
# Nichols chart grid
+from nichols import nichols_grid
def ngrid():
- """Nichols chart grid.
-
- Examples
- --------
- >>> ngrid()
- """
- from nichols import nichols_grid
nichols_grid()
+ngrid.__doc__ = re.sub('nichols_grid', 'ngrid', nichols_grid.__doc__)
# Root locus plot
def rlocus(sys, klist = None, **keywords):
@@ -1018,12 +1014,12 @@
klist:
list of gains used to compute roots
"""
- from rlocus import RootLocus
+ from rlocus import root_locus
if (klist == None):
#! TODO: update with a smart cacluation of the gains
klist = logspace(-3, 3)
- rlist = RootLocus(sys, klist, **keywords)
+ rlist = root_locus(sys, klist, **keywords)
return rlist, klist
def margin(*args):
@@ -1060,9 +1056,9 @@
"""
if len(args) == 1:
sys = args[0]
- margin = margins.StabilityMargins(sys)
+ margin = margins.stability_margins(sys)
elif len(args) == 3:
- margin = margins.StabilityMargins(args)
+ margin = margins.stability_margins(args)
else:
raise ValueError("Margin needs 1 or 3 arguments; received %i."
% len(args))
@@ -1157,7 +1153,7 @@
**keywords:
Additional keyword arguments control the solution algorithm for the
differential equations. These arguments are passed on to the function
- :func:`control.ForcedResponse`, which in turn passes them on to
+ :func:`control.forced_response`, which in turn passes them on to
:func:`scipy.integrate.odeint`. See the documentation for
:func:`scipy.integrate.odeint` for information about these
arguments.
@@ -1178,7 +1174,7 @@
--------
>>> T, yout = step(sys, T, X0)
'''
- T, yout = timeresp.StepResponse(sys, T, X0, input, output,
+ T, yout = timeresp.step_response(sys, T, X0, input, output,
transpose = True, **keywords)
return yout, T
@@ -1228,7 +1224,7 @@
--------
>>> T, yout = impulse(sys, T)
'''
- T, yout = timeresp.ImpulseResponse(sys, T, 0, input, output,
+ T, yout = timeresp.impulse_response(sys, T, 0, input, output,
transpose = True, **keywords)
return yout, T
@@ -1284,7 +1280,7 @@
--------
>>> T, yout = initial(sys, T, X0)
'''
- T, yout = timeresp.InitialResponse(sys, T, X0, input, output,
+ T, yout = timeresp.initial_response(sys, T, X0, input, output,
transpose = True, **keywords)
return yout, T
@@ -1338,6 +1334,6 @@
--------
>>> T, yout, xout = lsim(sys, U, T, X0)
'''
- T, yout, xout = timeresp.ForcedResponse(sys, T, U, X0,
+ T, yout, xout = timeresp.forced_response(sys, T, U, X0,
transpose = True, **keywords)
return yout, T, xout
Modified: trunk/src/nichols.py
===================================================================
--- trunk/src/nichols.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/src/nichols.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -46,7 +46,7 @@
from freqplot import default_frequency_range
# Nichols plot
-def nichols(syslist, omega=None, grid=True):
+def nichols_plot(syslist, omega=None, grid=True):
"""Nichols plot for a system
Plots a Nichols plot for the system over a (optional) frequency range.
@@ -294,3 +294,6 @@
mags = sp.linspace(10**(mag_min/20.0), 10**(mag_max/20.0), 2000)
Gcl_phases, Gcl_mags = sp.meshgrid(sp.radians(phases), mags)
return closed_loop_contours(Gcl_mags, Gcl_phases)
+
+# Function aliases
+nichols = nichols_plot
Modified: trunk/src/phaseplot.py
===================================================================
--- trunk/src/phaseplot.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/src/phaseplot.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -39,7 +39,7 @@
from exception import ControlNotImplemented
from scipy.integrate import odeint
-def PhasePlot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
+def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
lingrid=None, lintime=None, logtime=None, timepts=None,
parms=(), verbose=True):
"""
@@ -48,12 +48,12 @@
Produces a vector field or stream line plot for a planar system.
Call signatures:
- PhasePlot(func, X, Y, ...) - display vector field on meshgrid
- PhasePlot(func, X, Y, scale, ...) - scale arrows
- PhasePlot(func. X0=(...), T=Tmax, ...) - display stream lines
- PhasePlot(func, X, Y, X0=[...], T=Tmax, ...) - plot both
- PhasePlot(func, X0=[...], T=Tmax, lingrid=N, ...) - plot both
- PhasePlot(func, X0=[...], lintime=N, ...) - stream lines with arrows
+ phase_plot(func, X, Y, ...) - display vector field on meshgrid
+ phase_plot(func, X, Y, scale, ...) - scale arrows
+ phase_plot(func. X0=(...), T=Tmax, ...) - display stream lines
+ phase_plot(func, X, Y, X0=[...], T=Tmax, ...) - plot both
+ phase_plot(func, X0=[...], T=Tmax, lingrid=N, ...) - plot both
+ phase_plot(func, X0=[...], lintime=N, ...) - stream lines with arrows
Parameters
----------
@@ -104,7 +104,7 @@
See also
--------
- boxgrid(X, Y): construct box-shaped grid of initial conditions
+ box_grid(X, Y): construct box-shaped grid of initial conditions
Examples
--------
@@ -274,10 +274,10 @@
return;
# Utility function for generating initial conditions around a box
-def boxgrid(xlimp, ylimp):
- """BOXGRID generate list of points on edge of box
+def box_grid(xlimp, ylimp):
+ """box_grid generate list of points on edge of box
- list = BOXGRID([xmin xmax xnum], [ymin ymax ynum]) generates a
+ list = box_grid([xmin xmax xnum], [ymin ymax ynum]) generates a
list of points that correspond to a uniform grid at the end of the
box defined by the corners [xmin ymin] and [xmax ymax].
"""
Modified: trunk/src/rlocus.py
===================================================================
--- trunk/src/rlocus.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/src/rlocus.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -52,7 +52,7 @@
import xferfcn # transfer function manipulation
# Main function: compute a root locus diagram
-def RootLocus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True):
+def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True):
"""Calculate the root locus by finding the roots of 1+k*TF(s)
where TF is self.num(s)/self.den(s) and each k is an element
of kvect.
Modified: trunk/src/timeresp.py
===================================================================
--- trunk/src/timeresp.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/src/timeresp.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -53,8 +53,8 @@
This is a convention for function arguments and return values that
represent time series: sequences of values that change over time. It
is used throughout the library, for example in the functions
-:func:`ForcedResponse`, :func:`StepResponse`, :func:`ImpulseResponse`,
-and :func:`InitialResponse`.
+:func:`forced_response`, :func:`step_response`, :func:`impulse_response`,
+and :func:`initial_response`.
.. note::
This convention is different from the convention used in the library
@@ -235,7 +235,7 @@
return out_array
# Forced response of a linear system
-def ForcedResponse(sys, T=None, U=0., X0=0., transpose=False, **keywords):
+def forced_response(sys, T=None, U=0., X0=0., transpose=False, **keywords):
"""Simulate the output of a linear system.
As a convenience for parameters `U`, `X0`:
@@ -285,11 +285,11 @@
See Also
--------
- StepResponse, InitialResponse, ImpulseResponse
+ step_response, initial_response, impulse_response
Examples
--------
- >>> T, yout, xout = ForcedResponse(sys, T, u, X0)
+ >>> T, yout, xout = forced_response(sys, T, u, X0)
"""
if not isinstance(sys, Lti):
raise TypeError('Parameter ``sys``: must be a ``Lti`` object. '
@@ -365,7 +365,7 @@
return T, yout, xout
-def StepResponse(sys, T=None, X0=0., input=0, output=0, \
+def step_response(sys, T=None, X0=0., input=0, output=0, \
transpose = False, **keywords):
#pylint: disable=W0622
"""Step response of a linear system
@@ -419,11 +419,11 @@
See Also
--------
- ForcedResponse, InitialResponse, ImpulseResponse
+ forced_response, initial_response, impulse_response
Examples
--------
- >>> T, yout = StepResponse(sys, T, X0)
+ >>> T, yout = step_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
sys = _mimo2siso(sys, input, output, warn_conversion=True)
@@ -431,13 +431,13 @@
T = _default_response_times(sys.A, 100)
U = np.ones_like(T)
- T, yout, _xout = ForcedResponse(sys, T, U, X0,
+ T, yout, _xout = forced_response(sys, T, U, X0,
transpose=transpose, **keywords)
return T, yout
-def InitialResponse(sys, T=None, X0=0., input=0, output=0, transpose=False,
+def initial_response(sys, T=None, X0=0., input=0, output=0, transpose=False,
**keywords):
#pylint: disable=W0622
"""Initial condition response of a linear system
@@ -491,26 +491,26 @@
See Also
--------
- ForcedResponse, ImpulseResponse, StepResponse
+ forced_response, impulse_response, step_response
Examples
--------
- >>> T, yout = InitialResponse(sys, T, X0)
+ >>> T, yout = initial_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
sys = _mimo2siso(sys, input, output, warn_conversion=True)
- #Create time and input vectors; checking is done in ForcedResponse(...)
- #The initial vector X0 is created in ForcedResponse(...) if necessary
+ #Create time and input vectors; checking is done in forced_response(...)
+ #The initial vector X0 is created in forced_response(...) if necessary
if T is None:
T = _default_response_times(sys.A, 100)
U = np.zeros_like(T)
- T, yout, _xout = ForcedResponse(sys, T, U, X0, transpose=transpose,
+ T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
**keywords)
return T, yout
-def ImpulseResponse(sys, T=None, X0=0., input=0, output=0,
+def impulse_response(sys, T=None, X0=0., input=0, output=0,
transpose=False, **keywords):
#pylint: disable=W0622
"""Impulse response of a linear system
@@ -564,11 +564,11 @@
See Also
--------
- ForcedReponse, InitialResponse, StepResponse
+ ForcedReponse, initial_response, step_response
Examples
--------
- >>> T, yout = ImpulseResponse(sys, T, X0)
+ >>> T, yout = impulse_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
sys = _mimo2siso(sys, input, output, warn_conversion=True)
@@ -597,7 +597,7 @@
T = _default_response_times(sys.A, 100)
U = np.zeros_like(T)
- T, yout, _xout = ForcedResponse(sys, T, U, new_X0, \
+ T, yout, _xout = forced_response(sys, T, U, new_X0, \
transpose=transpose, **keywords)
return T, yout
Modified: trunk/tests/margin_test.py
===================================================================
--- trunk/tests/margin_test.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/tests/margin_test.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -7,7 +7,7 @@
import numpy as np
from control.xferfcn import TransferFunction
from control.statesp import StateSpace
-from control.margin import *
+from control.margins import *
class TestMargin(unittest.TestCase):
"""These are tests for the margin commands in margin.py."""
@@ -18,25 +18,25 @@
self.sys3 = StateSpace([[1., 4.], [3., 2.]], [[1.], [-4.]],
[[1., 0.]], [[0.]])
- def testGainPhaseMargin(self):
- gm, pm, sm, wg, wp, ws = StabilityMargins(self.sys1);
- gm, pm, sm, wg, wp, ws = StabilityMargins(self.sys2);
- gm, pm, sm, wg, wp, ws = StabilityMargins(self.sys3);
+ def test_stability_margins(self):
+ gm, pm, sm, wg, wp, ws = stability_margins(self.sys1);
+ gm, pm, sm, wg, wp, ws = stability_margins(self.sys2);
+ gm, pm, sm, wg, wp, ws = stability_margins(self.sys3);
- def testPhaseCrossoverFrequencies(self):
- omega, gain = PhaseCrossoverFrequencies(self.sys2)
+ def test_phase_crossover_frequencies(self):
+ omega, gain = phase_crossover_frequencies(self.sys2)
np.testing.assert_array_almost_equal(omega, [1.73205, 0.])
np.testing.assert_array_almost_equal(gain, [-0.5, 0.25])
tf = TransferFunction([1],[1,1])
- omega, gain = PhaseCrossoverFrequencies(tf)
+ omega, gain = phase_crossover_frequencies(tf)
np.testing.assert_array_almost_equal(omega, [0.])
np.testing.assert_array_almost_equal(gain, [1.])
# testing MIMO, only (0,0) element is considered
tf = TransferFunction([[[1],[2]],[[3],[4]]],
[[[1, 2, 3, 4],[1,1]],[[1,1],[1,1]]])
- omega, gain = PhaseCrossoverFrequencies(tf)
+ omega, gain = phase_crossover_frequencies(tf)
np.testing.assert_array_almost_equal(omega, [1.73205081, 0.])
np.testing.assert_array_almost_equal(gain, [-0.5, 0.25])
Modified: trunk/tests/phaseplot_test.py
===================================================================
--- trunk/tests/phaseplot_test.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/tests/phaseplot_test.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -21,18 +21,18 @@
pass;
def testInvPendNoSims(self):
- PhasePlot(self.invpend_ode, (-6,6,10), (-6,6,10));
+ phase_plot(self.invpend_ode, (-6,6,10), (-6,6,10));
def testInvPendSims(self):
- PhasePlot(self.invpend_ode, (-6,6,10), (-6,6,10),
+ phase_plot(self.invpend_ode, (-6,6,10), (-6,6,10),
X0 = ([1,1], [-1,1]));
def testInvPendTimePoints(self):
- PhasePlot(self.invpend_ode, (-6,6,10), (-6,6,10),
+ phase_plot(self.invpend_ode, (-6,6,10), (-6,6,10),
X0 = ([1,1], [-1,1]), T=np.linspace(0,5,100));
def testInvPendLogtime(self):
- PhasePlot(self.invpend_ode, X0 =
+ phase_plot(self.invpend_ode, X0 =
[ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
[-1, 2.1], [4.2, 2.1], [5, 2.1],
[2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
@@ -42,12 +42,12 @@
verbose=False)
def testInvPendAuto(self):
- PhasePlot(self.invpend_ode, lingrid = 0, X0=
+ phase_plot(self.invpend_ode, lingrid = 0, X0=
[[-2.3056, 2.1], [2.3056, -2.1]], T=6, verbose=False)
def testOscillatorParams(self):
m = 1; b = 1; k = 1; # default values
- PhasePlot(self.oscillator_ode, timepts = [0.3, 1, 2, 3], X0 =
+ phase_plot(self.oscillator_ode, timepts = [0.3, 1, 2, 3], X0 =
[[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1],
[1,1], [1.3,1], [1,-1], [0.3,-1], [0,-1], [-0.25,-1],
[-0.5,-1], [-0.7,-1], [-1,-1], [-1.3,-1]],
Modified: trunk/tests/test_all.py
===================================================================
--- trunk/tests/test_all.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/tests/test_all.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -1,6 +1,6 @@
#!/usr/bin/env python
#
-# test_al.py - test suit for python-control
+# test_all.py - test suit for python-control
# RMM, 30 Mar 2011
import unittest # unit test module
Modified: trunk/tests/timeresp_test.py
===================================================================
--- trunk/tests/timeresp_test.py 2011-07-26 05:44:01 UTC (rev 171)
+++ trunk/tests/timeresp_test.py 2011-08-07 15:12:58 UTC (rev 172)
@@ -43,7 +43,7 @@
"0. 9. ")
self.mimo_ss1 = StateSpace(A, B, C, D)
- def testStepResponse(self):
+ def test_step_response(self):
#Test SISO system
sys = self.siso_ss1
t = np.linspace(0, 1, 10)
@@ -51,64 +51,64 @@
42.3227, 44.9694, 47.1599, 48.9776])
# SISO call
- tout, yout = StepResponse(sys, T=t)
+ tout, yout = step_response(sys, T=t)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
np.testing.assert_array_almost_equal(tout, t)
# Play with arguments
- tout, yout = StepResponse(sys, T=t, X0=0)
+ tout, yout = step_response(sys, T=t, X0=0)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
np.testing.assert_array_almost_equal(tout, t)
X0 = np.array([0, 0]);
- tout, yout = StepResponse(sys, T=t, X0=X0)
+ tout, yout = step_response(sys, T=t, X0=X0)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
np.testing.assert_array_almost_equal(tout, t)
# Test MIMO system, which contains ``siso_ss1`` twice
sys = self.mimo_ss1
- _t, y_00 = StepResponse(sys, T=t, input=0, output=0)
- _t, y_11 = StepResponse(sys, T=t, input=1, output=1)
+ _t, y_00 = step_response(sys, T=t, input=0, output=0)
+ _t, y_11 = step_response(sys, T=t, input=1, output=1)
np.testing.assert_array_almost_equal(y_00, youttrue, decimal=4)
np.testing.assert_array_almost_equal(y_11, youttrue, decimal=4)
- def testImpulseResponse(self):
+ def test_impulse_response(self):
#Test SISO system
sys = self.siso_ss1
t = np.linspace(0, 1, 10)
youttrue = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766, 31.7344,
26.1668, 21.6292, 17.9245, 14.8945])
- tout, yout = ImpulseResponse(sys, T=t)
+ tout, yout = impulse_response(sys, T=t)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
np.testing.assert_array_almost_equal(tout, t)
#Test MIMO system, which contains ``siso_ss1`` twice
sys = self.mimo_ss1
- _t, y_00 = ImpulseResponse(sys, T=t, input=0, output=0)
- _t, y_11 = ImpulseResponse(sys, T=t, input=1, output=1)
+ _t, y_00 = impulse_response(sys, T=t, input=0, output=0)
+ _t, y_11 = impulse_response(sys, T=t, input=1, output=1)
np.testing.assert_array_almost_equal(y_00, youttrue, decimal=4)
np.testing.assert_array_almost_equal(y_11, youttrue, decimal=4)
- def testInitialResponse(self):
+ def test_initial_response(self):
#Test SISO system
sys = self.siso_ss1
t = np.linspace(0, 1, 10)
x0 = np.array([[0.5], [1]]);
youttrue = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118, 1.9092,
1.1508, 0.5833, 0.1645, -0.1391])
- tout, yout = InitialResponse(sys, T=t, X0=x0)
+ tout, yout = initial_response(sys, T=t, X0=x0)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
np.testing.assert_array_almost_equal(tout, t)
#Test MIMO system, which contains ``siso_ss1`` twice
sys = self.mimo_ss1
x0 = np.matrix(".5; 1.; .5; 1.")
- _t, y_00 = InitialResponse(sys, T=t, X0=x0, input=0, output=0)
- _t, y_11 = InitialResponse(sys, T=t, X0=x0, input=1, output=1)
+ _t, y_00 = initial_response(sys, T=t, X0=x0, input=0, output=0)
+ _t, y_11 = initial_response(sys, T=t, X0=x0, input=1, output=1)
np.testing.assert_array_almost_equal(y_00, youttrue, decimal=4)
np.testing.assert_array_almost_equal(y_11, youttrue, decimal=4)
- def testForcedResponse(self):
+ def test_forced_response(self):
t = np.linspace(0, 1, 10)
#compute step response - test with state space, and transfer function
@@ -116,10 +116,10 @@
u = np.array([1., 1, 1, 1, 1, 1, 1, 1, 1, 1])
youttrue = np.array([9., 17.6457, 24.7072, 30.4855, 35.2234, 39.1165,
42.3227, 44.9694, 47.1599, 48.9776])
- tout, yout, _xout = ForcedResponse(self.siso_ss1, t, u)
+ tout, yout, _xout = forced_response(self.siso_ss1, t, u)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
np.testing.assert_array_almost_equal(tout, t)
- _t, yout, _xout = ForcedResponse(self.siso_tf2, t, u)
+ _t, yout, _xout = forced_response(self.siso_tf2, t, u)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
#test with initial value and special algorithm for ``U=0``
@@ -127,7 +127,7 @@
x0 = np.matrix(".5; 1.")
youttrue = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118, 1.9092,
1.1508, 0.5833, 0.1645, -0.1391])
- _t, yout, _xout = ForcedResponse(self.siso_ss1, t, u, x0)
+ _t, yout, _xout = forced_response(self.siso_ss1, t, u, x0)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
#Test MIMO system, which contains ``siso_ss1`` twice
@@ -139,7 +139,7 @@
1.1508, 0.5833, 0.1645, -0.1391],
[9., 17.6457, 24.7072, 30.4855, 35.2234, 39.1165,
42.3227, 44.9694, 47.1599, 48.9776]])
- _t, yout, _xout = ForcedResponse(self.mimo_ss1, t, u, x0)
+ _t, yout, _xout = forced_response(self.mimo_ss1, t, u, x0)
np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
def suite():
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X