You can subscribe to this list here.
2010 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(19) |
Jul
|
Aug
|
Sep
|
Oct
|
Nov
(11) |
Dec
(5) |
---|---|---|---|---|---|---|---|---|---|---|---|---|
2011 |
Jan
|
Feb
(113) |
Mar
(12) |
Apr
(27) |
May
(2) |
Jun
(16) |
Jul
(6) |
Aug
(6) |
Sep
|
Oct
(3) |
Nov
(9) |
Dec
(2) |
2012 |
Jan
(5) |
Feb
(11) |
Mar
|
Apr
(3) |
May
|
Jun
|
Jul
|
Aug
(3) |
Sep
(7) |
Oct
(18) |
Nov
(18) |
Dec
|
2013 |
Jan
(4) |
Feb
(1) |
Mar
(3) |
Apr
(1) |
May
|
Jun
(33) |
Jul
(2) |
Aug
(5) |
Sep
|
Oct
|
Nov
|
Dec
(1) |
2014 |
Jan
(1) |
Feb
|
Mar
(8) |
Apr
|
May
(3) |
Jun
(3) |
Jul
(9) |
Aug
(5) |
Sep
(6) |
Oct
|
Nov
|
Dec
|
2015 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
(4) |
2017 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
(5) |
Nov
|
Dec
|
2018 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(1) |
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
2019 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
(6) |
Sep
|
Oct
|
Nov
(2) |
Dec
|
2020 |
Jan
(1) |
Feb
(1) |
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
2021 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
(4) |
Oct
|
Nov
|
Dec
|
From: Sawyer F. <mi...@ca...> - 2011-04-06 21:55:04
|
Hey i just downloaded the latest version of python-control and it looks like everything has gone MIMO. Pretty cool, but is there any way to way to stick to SISO, eg a subclass or flag? A little bit of a drag to constantly be appending [0][0] to every frequency response for the most common use case. I think matlab does a special-case for SISO. Anything like that in the works or somewhere i'm not finding it? S ------------------------ Sawyer B. Fuller Ph.D. Candidate, Bioengineering California Institute of Technology http://www.cds.caltech.edu/~minster/ <http://alumni.media.mit.edu/~minster/> |
From: <mur...@us...> - 2011-04-03 04:55:59
|
Revision: 152 http://python-control.svn.sourceforge.net/python-control/?rev=152&view=rev Author: murrayrm Date: 2011-04-03 04:55:51 +0000 (Sun, 03 Apr 2011) Log Message: ----------- Small fixes based on testing + re-added root locus * Implemented syntax checks for MATLAB functions * Fixed root locus plots to work with new SF, TF classes * Transfer functions with integer coeffs now convered to floats * More detailed list of changes in ChangeLog Modified Paths: -------------- trunk/ChangeLog trunk/Pending trunk/examples/secord-matlab.py trunk/src/__init__.py trunk/src/exception.py trunk/src/freqplot.py trunk/src/matlab.py trunk/src/rlocus.py trunk/src/xferfcn.py trunk/tests/matlab_test.py Modified: trunk/ChangeLog =================================================================== --- trunk/ChangeLog 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/ChangeLog 2011-04-03 04:55:51 UTC (rev 152) @@ -1,5 +1,38 @@ 2011-04-02 Richard Murray <murray@malabar.local> + * src/xferfcn.py (TransferFunction.__add__): fixed bug when adding a + transfer function to state space system; _convertToTransferFunction + was being called with input/output keywords. Picked up in unit test. + + * tests/matlab_test.py: added calls to all of the functions that are + currently implemented in the library, to make sure there are no + hidden issues. These calls do *not* test functionality, they just + make sure that MATLAB compatibility functions accept the right types + of arguments. + + * examples/secord-matlab.py: added root locus plot to list of + figures that are produced + + * src/__init__.py: added rlocus to list of modules that are imported + by control module + + * src/exception.py (ControlMIMONotImplemented): added exception for + functions that are not yet implemented for MIMO systems + + * src/xferfcn.py (TransferFunction.__init__): convert integer + numerator and denominator objects to floats to eliminate errors when + creating common denominator (also used on pole command). This fixes + and error were tf([1], [1, 2, 1]).pole() would generate an error. + + * src/freqplot.py (bode): Tweaked documentation string to remove 'h' + from mag and phase return arguments + + * src/rlocus.py (RootLocus): removed commands that set figure number + inside of RootLocus command, for consistency with other frequency + plot routines. Added Plot=True argument and updated documentation. + + * src/matlab.py: added rlocus() command (calls RootLocus) + * MANIFEST.in: Added MANIFEST.in file to include examples and tests in source distribution Modified: trunk/Pending =================================================================== --- trunk/Pending 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/Pending 2011-04-03 04:55:51 UTC (rev 152) @@ -12,7 +12,7 @@ * matlab.step() doesn't handle systems with a pole at the origin (use lsim2) * TF <-> SS transformations are buggy; see tests/convert_test.py * hsvd returns different value than MATLAB (2010a); see modelsimp_test.py - * MIMO common denominator fails unit test; see convert_test.py + * lsim doesn't work for StateSpace systems (signal.lsim2 bug??) Transfer code from Roberto Bucher's yottalab to python-control acker - pole placement using Ackermann method @@ -44,23 +44,15 @@ * tests/freqresp.py needs to be converted to unit test * Convert examples/test-{response,statefbk}.py to unit tests -TransferFunction class fixes - * evalfr is not working (num, den stored as ndarrays, not poly1ds) +Root locus plot improvements + * Make sure that scipy.signal.lti objects still work + * Update calling syntax to be consistent with other plotting commands -Block diagram algebra fixes - * Implement state space block diagram algebra - * Produce minimal realizations to avoid later errors - State space class fixes - * Convert pvtol to state space systems and rerun example * Implement pzmap for state space systems -LTI updates - * Implement control.matlab.step (with semantics similar to MATLAB) - Basic functions to be added * margin - compute gain and phase margin (no plot) - * lqr - compute optimal feedback gains (use SLICOT SB02ND.f) * lyap - solve Lyapunov equation (use SLICOT SB03MD.f) * See http://www.slicot.org/shared/libindex.html for list of functions Modified: trunk/examples/secord-matlab.py =================================================================== --- trunk/examples/secord-matlab.py 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/examples/secord-matlab.py 2011-04-03 04:55:51 UTC (rev 152) @@ -27,3 +27,7 @@ # Nyquist plot for the system figure(3) nyquist(sys, logspace(-2, 2)) + +# Root lcous plut for the system +figure(4) +rlocus(sys) Modified: trunk/src/__init__.py =================================================================== --- trunk/src/__init__.py 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/src/__init__.py 2011-04-03 04:55:51 UTC (rev 152) @@ -64,3 +64,4 @@ from statefbk import * from delay import * from modelsimp import * +from rlocus import * Modified: trunk/src/exception.py =================================================================== --- trunk/src/exception.py 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/src/exception.py 2011-04-03 04:55:51 UTC (rev 152) @@ -52,6 +52,10 @@ """Raised when arguments to a function are not correct""" pass +class ControlMIMONotImplemented(Exception): + """Function is not currently implemented for MIMO systems""" + pass + class ControlNotImplemented(Exception): """Functionality is not yet implemented""" pass Modified: trunk/src/freqplot.py =================================================================== --- trunk/src/freqplot.py 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/src/freqplot.py 2011-04-03 04:55:51 UTC (rev 152) @@ -60,7 +60,7 @@ Usage ===== - (magh, phaseh, omega) = bode(syslist, omega=None, dB=False, Hz=False, color=None, Plot=True) + (mag, phase, omega) = bode(syslist, omega=None, dB=False, Hz=False, color=None, Plot=True) Plots a Bode plot for the system over a (optional) frequency range. @@ -79,8 +79,8 @@ Return values ------------- - magh : magnitude array - phaseh : phase array + mag : magnitude array + phase : phase array omega : frequency array Notes Modified: trunk/src/matlab.py =================================================================== --- trunk/src/matlab.py 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/src/matlab.py 2011-04-03 04:55:51 UTC (rev 152) @@ -179,7 +179,7 @@ \* ss/modred - model order reduction Compensator design - rlocus - evans root locus +\* rlocus - evans root locus \* place - pole placement estim - form estimator given estimator gain reg - form regulator given state-feedback and estimator gains @@ -187,7 +187,7 @@ LQR/LQG design ss/lqg - single-step LQG design \* lqr - linear-Quadratic (LQ) state-feedback regulator -\* dlqr - discrete-time LQ state-feedback regulator + dlqr - discrete-time LQ state-feedback regulator lqry - lq regulator with output weighting lqrd - discrete LQ regulator for continuous plant ss/lqi - linear-Quadratic-Integral (LQI) controller @@ -243,8 +243,8 @@ Overloaded arithmetic operations \* + and - - add and subtract systems (parallel connection) \* \* - multiply systems (series connection) -\* / - left divide -- sys1\sys2 means inv(sys1)\*sys2 -\* / - right divide -- sys1/sys2 means sys1\*inv(sys2) + / - left divide -- sys1\sys2 means inv(sys1)\*sys2 +- \ - right divide -- sys1/sys2 means sys1\*inv(sys2) ^ - powers of a given system ' - pertransposition .' - transposition of input/output map @@ -758,12 +758,35 @@ from nichols import nichols_grid nichols_grid() +# Root locus plot +def rlocus(sys, klist = None, **keywords): + """Root locus plot + + Parameters + ---------- + sys: StateSpace or TransferFunction object + klist: optional list of gains + + Returns + ------- + rlist: list of roots for each gain + klist: list of gains used to compute roots + """ + from rlocus import RootLocus + if (klist == None): + #! TODO: update with a smart cacluation of the gains + klist = logspace(-3, 3) + + rlist = RootLocus(sys, klist, **keywords) + return rlist, klist + + # # Modifications to scipy.signal functions # # Redefine lsim to use lsim2 -def lsim(*args, **keywords): +def lsim(sys, U=None, T=None, X0=None, **keywords): """Simulate the output of a linear system Examples @@ -784,12 +807,13 @@ yout: response of the system xout: time evolution of the state vector """ - sys = args[0] + # Convert the system to an signal.lti for simulation + #! This should send a warning for MIMO systems ltiobjs = sys.returnScipySignalLti() ltiobj = ltiobjs[0][0] - - return sp.signal.lsim2(ltiobj, **keywords) + return sp.signal.lsim2(ltiobj, U, T, X0, **keywords) + #! Redefine step to use lsim2 #! Not yet implemented def step(*args, **keywords): Modified: trunk/src/rlocus.py =================================================================== --- trunk/src/rlocus.py 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/src/rlocus.py 2011-04-03 04:55:51 UTC (rev 152) @@ -37,65 +37,93 @@ # * Added BSD copyright info to file (per Ryan) # * Added code to convert (num, den) to poly1d's if they aren't already. # This allows Ryan's code to run on a standard signal.ltisys object -# or a controls.TransferFunction object. +# or a control.TransferFunction object. # * Added some comments to make sure I understand the code +# +# RMM, 2 April 2011: modified to work with new Lti structure (see ChangeLog) +# * Not tested: should still work on signal.ltisys objects # # $Id$ # Packages used by this module from scipy import * +import scipy.signal # signal processing toolbox +import pylab # plotting routines +import xferfcn # transfer function manipulation # Main function: compute a root locus diagram -def RootLocus(sys, kvect, fig=None, fignum=1, \ - clear=True, xlim=None, ylim=None, plotstr='-'): +def RootLocus(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.""" + of kvect. + Parameters + ---------- + sys : linsys + Linear input/output systems (SISO only, for now) + klist : gain_range (default = None) + List of gains to use in computing diagram + Plot : boolean (default = True) + If True, plot magnitude and phase + + Return values + ------------- + rlist : list of computed root locations + klist : list of gains + """ + # Convert numerator and denominator to polynomials if they aren't (nump, denp) = _systopoly1d(sys); - # Set up the figure - if fig is None: - import pylab - fig = pylab.figure(fignum) - if clear: - fig.clf() - ax = fig.add_subplot(111) - # Compute out the loci mymat = _RLFindRoots(sys, kvect) mymat = _RLSortRoots(sys, mymat) - # plot open loop poles - poles = array(denp.r) - ax.plot(real(poles), imag(poles), 'x') + # Create the plot + if (Plot): + ax = pylab.axes(); - # plot open loop zeros - zeros = array(nump.r) - if zeros.any(): - ax.plot(real(zeros), imag(zeros), 'o') + # plot open loop poles + poles = array(denp.r) + ax.plot(real(poles), imag(poles), 'x') - # Now plot the loci - for col in mymat.T: - ax.plot(real(col), imag(col), plotstr) + # plot open loop zeros + zeros = array(nump.r) + if zeros.any(): + ax.plot(real(zeros), imag(zeros), 'o') - # Set up plot axes and labels - if xlim: - ax.set_xlim(xlim) - if ylim: - ax.set_ylim(ylim) - ax.set_xlabel('Real') - ax.set_ylabel('Imaginary') + # Now plot the loci + for col in mymat.T: + ax.plot(real(col), imag(col), plotstr) + + # Set up plot axes and labels + if xlim: + ax.set_xlim(xlim) + if ylim: + ax.set_ylim(ylim) + ax.set_xlabel('Real') + ax.set_ylabel('Imaginary') + return mymat # Utility function to extract numerator and denominator polynomials def _systopoly1d(sys): """Extract numerator and denominator polynomails for a system""" + # Allow inputs from the signal processing toolbox + if (isinstance(sys, scipy.signal.lti)): + nump = sys.num; denp = sys.den; - # Start by extracting the numerator and denominator from system object - nump = sys.num; denp = sys.den; + else: + # Convert to a transfer function, if needed + sys = xferfcn._convertToTransferFunction(sys) + # Make sure we have a SISO system + if (sys.inputs > 1 or sys.outputs > 1): + raise ControlMIMONotImplemented() + + # Start by extracting the numerator and denominator from system object + nump = sys.num[0][0]; denp = sys.den[0][0]; + # Check to see if num, den are already polynomials; otherwise convert if (not isinstance(nump, poly1d)): nump = poly1d(nump) if (not isinstance(denp, poly1d)): denp = poly1d(denp) Modified: trunk/src/xferfcn.py =================================================================== --- trunk/src/xferfcn.py 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/src/xferfcn.py 2011-04-03 04:55:51 UTC (rev 152) @@ -129,11 +129,20 @@ for i in range(len(data)): if isinstance(data[i], (int, float, long, complex)): # Convert scalar to list of list of array. - data[i] = [[array([data[i]])]] + if (isinstance(data[i], int)): + # Convert integers to floats at this point + data[i] = [[array([data[i]], dtype=float)]] + else: + data[i] = [[array([data[i]])]] elif (isinstance(data[i], (list, tuple, ndarray)) and isinstance(data[i][0], (int, float, long, complex))): # Convert array to list of list of array. - data[i] = [[array(data[i])]] + if (isinstance(data[i][0], int)): + # Convert integers to floats at this point + #! Not sure this covers all cases correctly + data[i] = [[array(data[i], dtype=float)]] + else: + data[i] = [[array(data[i])]] elif (isinstance(data[i], list) and isinstance(data[i][0], list) and isinstance(data[i][0][0], (list, tuple, ndarray)) and @@ -142,7 +151,10 @@ # coefficient vectors to arrays, if necessary. for j in range(len(data[i])): for k in range(len(data[i][j])): - data[i][j][k] = array(data[i][j][k]) + if (isinstance(data[i][j][k], int)): + data[i][j][k] = array(data[i][j][k], dtype=float) + else: + data[i][j][k] = array(data[i][j][k]) else: # If the user passed in anything else, then it's unclear what # the meaning is. @@ -274,7 +286,9 @@ """Add two LTI objects (parallel connection).""" # Convert the second argument to a transfer function. - if not isinstance(other, TransferFunction): + if (isinstance(other, statesp.StateSpace)): + other = _convertToTransferFunction(other) + elif not isinstance(other, TransferFunction): other = _convertToTransferFunction(other, inputs=self.inputs, outputs=self.outputs) @@ -513,6 +527,7 @@ # A sorted list to keep track of cumulative poles found as we scan # self.den. poles = [] + # A 3-D list to keep track of common denominator poles not present in # the self.den[i][j]. missingpoles = [[[] for j in range(self.inputs)] for i in Modified: trunk/tests/matlab_test.py =================================================================== --- trunk/tests/matlab_test.py 2011-04-02 19:32:33 UTC (rev 151) +++ trunk/tests/matlab_test.py 2011-04-03 04:55:51 UTC (rev 152) @@ -2,47 +2,229 @@ # # matlab_test.py - test MATLAB compatibility # RMM, 30 Mar 2011 (based on TestMatlab from v0.4a) +# +# This test suite just goes through and calls all of the MATLAB +# functions using different systems and arguments to make sure that +# nothing crashes. It doesn't test actual functionality; the module +# specific unit tests will do that. import unittest import numpy as np from control.matlab import * class TestMatlab(unittest.TestCase): - def testStep(self): + def setUp(self): + """Set up some systems for testing out MATLAB functions""" A = np.matrix("1. -2.; 3. -4.") B = np.matrix("5.; 7.") C = np.matrix("6. 8.") D = np.matrix("9.") - sys = ss(A,B,C,D) + self.siso_ss1 = ss(A,B,C,D) + + # Create some transfer functions + self.siso_tf1 = tf([1], [1, 2, 1]); + self.siso_tf2 = tf([1, 1], [1, 2, 3, 1]); + + # Conversions + self.siso_tf3 = tf(self.siso_ss1); + self.siso_ss2 = ss(self.siso_tf2); + self.siso_ss3 = tf2ss(self.siso_tf3); + self.siso_tf4 = ss2tf(self.siso_ss2); + + def testParallel(self): + sys1 = parallel(self.siso_ss1, self.siso_ss2) + sys1 = parallel(self.siso_ss1, self.siso_tf2) + sys1 = parallel(self.siso_tf1, self.siso_ss2) + sys1 = parallel(1, self.siso_ss2) + sys1 = parallel(1, self.siso_tf2) + sys1 = parallel(self.siso_ss1, 1) + sys1 = parallel(self.siso_tf1, 1) + + def testSeries(self): + sys1 = series(self.siso_ss1, self.siso_ss2) + sys1 = series(self.siso_ss1, self.siso_tf2) + sys1 = series(self.siso_tf1, self.siso_ss2) + sys1 = series(1, self.siso_ss2) + sys1 = series(1, self.siso_tf2) + sys1 = series(self.siso_ss1, 1) + sys1 = series(self.siso_tf1, 1) + + def testFeedback(self): + sys1 = feedback(self.siso_ss1, self.siso_ss2) + sys1 = feedback(self.siso_ss1, self.siso_tf2) + sys1 = feedback(self.siso_tf1, self.siso_ss2) + sys1 = feedback(1, self.siso_ss2) + sys1 = feedback(1, self.siso_tf2) + sys1 = feedback(self.siso_ss1, 1) + sys1 = feedback(self.siso_tf1, 1) + + def testPoleZero(self): + pole(self.siso_ss1); + pole(self.siso_tf1); + pole(self.siso_tf2); + zero(self.siso_ss1); + zero(self.siso_tf1); + zero(self.siso_tf2); + + def testPZmap(self): + # pzmap(self.siso_ss1); not implemented + # pzmap(self.siso_ss2); not implemented + pzmap(self.siso_tf1); + pzmap(self.siso_tf2); + pzmap(self.siso_tf2, Plot=False); + + def testStep(self): + sys = self.siso_ss1 t = np.linspace(0, 1, 10) t, yout = step(sys, T=t) youttrue = np.matrix("9. 17.6457 24.7072 30.4855 35.2234 39.1165 42.3227 44.9694 47.1599 48.9776") np.testing.assert_array_almost_equal(yout, youttrue,decimal=4) def testImpulse(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") - sys = ss(A,B,C,D) + sys = self.siso_ss1 t = np.linspace(0, 1, 10) t, yout = impulse(sys, T=t) youttrue = np.matrix("86. 70.1808 57.3753 46.9975 38.5766 31.7344 26.1668 21.6292 17.9245 14.8945") np.testing.assert_array_almost_equal(yout, youttrue,decimal=4) # def testInitial(self): -# A = np.matrix("1. -2.; 3. -4.") -# B = np.matrix("5.; 7.") -# C = np.matrix("6. 8.") -# D = np.matrix("9.") -# sys = ss(A,B,C,D) + sys = self.siso_ss1 # t = np.linspace(0, 1, 10) # x0 = np.matrix(".5; 1.") # t, yout = initial(sys, T=t, X0=x0) # youttrue = np.matrix("11. 8.1494 5.9361 4.2258 2.9118 1.9092 1.1508 0.5833 0.1645 -0.1391") # np.testing.assert_array_almost_equal(yout, youttrue,decimal=4) + def testLsim(self): + T = range(1, 100) + u = np.sin(T) + lsim(self.siso_tf1, u, T) + # lsim(self.siso_ss1, u, T) # generates error?? + # lsim(self.siso_ss1, u, T, self.siso_ss1.B) + def testBode(self): + bode(self.siso_ss1) + bode(self.siso_tf1) + bode(self.siso_tf2) + (mag, phase, freq) = bode(self.siso_tf2, Plot=False) + bode(self.siso_tf1, self.siso_tf2) + w = logspace(-3, 3); + bode(self.siso_ss1, w) + bode(self.siso_ss1, self.siso_tf2, w) + bode(self.siso_ss1, '-', self.siso_tf1, 'b--', self.siso_tf2, 'k.') + + def testRlocus(self): + rlocus(self.siso_ss1) + rlocus(self.siso_tf1) + rlocus(self.siso_tf2) + rlist, klist = rlocus(self.siso_tf2, klist=[1, 10, 100], Plot=False) + + def testNyquist(self): + nyquist(self.siso_ss1) + nyquist(self.siso_tf1) + nyquist(self.siso_tf2) + w = logspace(-3, 3); + nyquist(self.siso_tf2, w) + (real, imag, freq) = nyquist(self.siso_tf2, w, Plot=False) + + def testNichols(self): + nichols(self.siso_ss1) + nichols(self.siso_tf1) + nichols(self.siso_tf2) + w = logspace(-3, 3); + nichols(self.siso_tf2, w) + nichols(self.siso_tf2, grid=False) + + def testFreqresp(self): + w = logspace(-3, 3) + freqresp(self.siso_ss1, w) + freqresp(self.siso_ss2, w) + freqresp(self.siso_ss3, w) + freqresp(self.siso_tf1, w) + freqresp(self.siso_tf2, w) + freqresp(self.siso_tf3, w) + + def testEvalfr(self): + w = 1 + evalfr(self.siso_ss1, w) + evalfr(self.siso_ss2, w) + evalfr(self.siso_ss3, w) + evalfr(self.siso_tf1, w) + evalfr(self.siso_tf2, w) + evalfr(self.siso_tf3, w) + + def testHsvd(self): + hsvd(self.siso_ss1) + hsvd(self.siso_ss2) + hsvd(self.siso_ss3) + + def testBalred(self): + balred(self.siso_ss1, 1) + balred(self.siso_ss2, 2) + balred(self.siso_ss3, [2, 2]) + + def testModred(self): + modred(self.siso_ss1, [1]) + modred(self.siso_ss2 * self.siso_ss3, [1, 2]) + modred(self.siso_ss3, [1], 'matchdc') + modred(self.siso_ss3, [1], 'truncate') + + def testPlace(self): + place(self.siso_ss1.A, self.siso_ss1.B, [-2, -2]) + + def testLQR(self): + (K, S, E) = lqr(self.siso_ss1.A, self.siso_ss1.B, np.eye(2), np.eye(1)) + (K, S, E) = lqr(self.siso_ss2.A, self.siso_ss2.B, np.eye(3), \ + np.eye(1), [[1], [1], [2]]) + + def testRss(self): + rss(1) + rss(2) + rss(2, 3, 1) + + def testDrss(self): + drss(1) + drss(2) + drss(2, 3, 1) + + def testCtrb(self): + ctrb(self.siso_ss1.A, self.siso_ss1.B) + ctrb(self.siso_ss2.A, self.siso_ss2.B) + + def testObsv(self): + obsv(self.siso_ss1.A, self.siso_ss1.C) + obsv(self.siso_ss2.A, self.siso_ss2.C) + + def testGram(self): + gram(self.siso_ss1, 'c') + gram(self.siso_ss2, 'c') + gram(self.siso_ss1, 'o') + gram(self.siso_ss2, 'o') + + def testPade(self): + pade(1, 1) + pade(1, 2) + pade(5, 4) + + def testOpers(self): + self.siso_ss1 + self.siso_ss2 + self.siso_tf1 + self.siso_tf2 + self.siso_ss1 + self.siso_tf2 + self.siso_tf1 + self.siso_ss2 + self.siso_ss1 * self.siso_ss2 + self.siso_tf1 * self.siso_tf2 + self.siso_ss1 * self.siso_tf2 + self.siso_tf1 * self.siso_ss2 + # self.siso_ss1 / self.siso_ss2 not implemented yet + # self.siso_tf1 / self.siso_tf2 + # self.siso_ss1 / self.siso_tf2 + # self.siso_tf1 / self.siso_ss2 + + def testUnwrap(self): + phase = np.array(range(1, 100)) / 10.; + wrapped = phase % (2 * np.pi) + unwrapped = unwrap(wrapped) + def suite(): return unittest.TestLoader().loadTestsFromTestCase(TestMatlab) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-04-02 19:32:40
|
Revision: 151 http://python-control.svn.sourceforge.net/python-control/?rev=151&view=rev Author: murrayrm Date: 2011-04-02 19:32:33 +0000 (Sat, 02 Apr 2011) Log Message: ----------- added examples/ and tests/ to source distribution; updated version number to 0.4c Modified Paths: -------------- trunk/ChangeLog trunk/README trunk/setup.py Added Paths: ----------- trunk/MANIFEST.in Modified: trunk/ChangeLog =================================================================== --- trunk/ChangeLog 2011-04-02 18:34:21 UTC (rev 150) +++ trunk/ChangeLog 2011-04-02 19:32:33 UTC (rev 151) @@ -1,5 +1,10 @@ 2011-04-02 Richard Murray <murray@malabar.local> + * MANIFEST.in: Added MANIFEST.in file to include examples and tests + in source distribution + + * README: Updated to include information on how to run unit tests. + * setup.py: updated version number to 0.4c ---- control-0.4b released ----- Added: trunk/MANIFEST.in =================================================================== --- trunk/MANIFEST.in (rev 0) +++ trunk/MANIFEST.in 2011-04-02 19:32:33 UTC (rev 151) @@ -0,0 +1,2 @@ +include examples/*.py +include tests/*.py Modified: trunk/README =================================================================== --- trunk/README 2011-04-02 18:34:21 UTC (rev 150) +++ trunk/README 2011-04-02 19:32:33 UTC (rev 151) @@ -16,3 +16,10 @@ examples/secord-matlab.py (using ipython -pylab). It should generate a step response, Bode plot and Nyquist plot for a simple second order linear system. + +You can also run a set of unit tests to make sure that everything is +working correctly. After installation, run + + python tests/test_all.py + +from the source distribution directory. Modified: trunk/setup.py =================================================================== --- trunk/setup.py 2011-04-02 18:34:21 UTC (rev 150) +++ trunk/setup.py 2011-04-02 19:32:33 UTC (rev 151) @@ -2,13 +2,13 @@ from distutils.core import setup -setup(name='control', - version='0.4c', - description='Python Control Systems Library', - author='Richard Murray', - author_email='mu...@cd...', - url='http://python-control.sourceforge.net', - requires='scipy', - packages=['control'], - package_dir = {'control':'src'}, +setup(name = 'control', + version = '0.4c', + description = 'Python Control Systems Library', + author = 'Richard Murray', + author_email = 'mu...@cd...', + url = 'http://python-control.sourceforge.net', + requires = ['scipy', 'matplotlib'], + package_dir = {'control' : 'src'}, + packages = ['control'], ) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: Richard M. <mu...@cd...> - 2011-04-02 18:58:55
|
freqplot.py is probably the place to put it. Make sure to update to the latest version of trunk (may require a fresh checkout since lots of things moved around). It would also be useful to put in a unit test for the code (see tests/ subdirectory for examples). -richard On 2 Apr 2011, at 11:56 , Sawyer Fuller wrote: > hey, i've got a nice margin command i've written up. (it doesn't plot anything yet) > > where do you think it goes? I'm thinking freqplot > > ------------------------ > Sawyer B. Fuller > Ph.D. Candidate, Bioengineering > California Institute of Technology > http://www.cds.caltech.edu/~minster/ > > > ------------------------------------------------------------------------------ > Create and publish websites with WebMatrix > Use the most popular FREE web apps or write code yourself; > WebMatrix provides all the features you need to develop and > publish your website. http://p.sf.net/sfu/ms-webmatrix-sf > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
From: Sawyer F. <mi...@ca...> - 2011-04-02 18:56:20
|
hey, i've got a nice margin command i've written up. (it doesn't plot anything yet) where do you think it goes? I'm thinking freqplot ------------------------ Sawyer B. Fuller Ph.D. Candidate, Bioengineering California Institute of Technology http://www.cds.caltech.edu/~minster/ <http://alumni.media.mit.edu/~minster/> |
From: <mur...@us...> - 2011-04-02 18:34:27
|
Revision: 150 http://python-control.svn.sourceforge.net/python-control/?rev=150&view=rev Author: murrayrm Date: 2011-04-02 18:34:21 +0000 (Sat, 02 Apr 2011) Log Message: ----------- tag of v0.4b + updated trunk to 0.4c Modified Paths: -------------- trunk/ChangeLog trunk/setup.py Added Paths: ----------- tags/control-0.4b/ Modified: trunk/ChangeLog =================================================================== --- trunk/ChangeLog 2011-04-02 18:24:51 UTC (rev 149) +++ trunk/ChangeLog 2011-04-02 18:34:21 UTC (rev 150) @@ -1,5 +1,11 @@ 2011-04-02 Richard Murray <murray@malabar.local> + * setup.py: updated version number to 0.4c + +---- control-0.4b released ----- + +2011-04-02 Richard Murray <murray@malabar.local> + * src/__init__.py: removed import of tests module (moved to tests/) * src/matlab.py: Added hsvd, balred, modred to list of functions Modified: trunk/setup.py =================================================================== --- trunk/setup.py 2011-04-02 18:24:51 UTC (rev 149) +++ trunk/setup.py 2011-04-02 18:34:21 UTC (rev 150) @@ -3,7 +3,7 @@ from distutils.core import setup setup(name='control', - version='0.4b', + version='0.4c', description='Python Control Systems Library', author='Richard Murray', author_email='mu...@cd...', This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-04-02 18:24:57
|
Revision: 149 http://python-control.svn.sourceforge.net/python-control/?rev=149&view=rev Author: murrayrm Date: 2011-04-02 18:24:51 +0000 (Sat, 02 Apr 2011) Log Message: ----------- making v0.4b new trunk Added Paths: ----------- trunk/ Removed Paths: ------------- branches/control-0.4b/ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-04-02 18:24:32
|
Revision: 148 http://python-control.svn.sourceforge.net/python-control/?rev=148&view=rev Author: murrayrm Date: 2011-04-02 18:24:26 +0000 (Sat, 02 Apr 2011) Log Message: ----------- moving trunk to branch (in case we need to update) Added Paths: ----------- branches/control-0.3/ Removed Paths: ------------- trunk/ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-04-02 18:23:22
|
Revision: 147 http://python-control.svn.sourceforge.net/python-control/?rev=147&view=rev Author: murrayrm Date: 2011-04-02 18:23:16 +0000 (Sat, 02 Apr 2011) Log Message: ----------- moving v0.4a branch to tags (Princeton changes) Added Paths: ----------- tags/control-0.4a/ Removed Paths: ------------- branches/control-0.4a/ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-04-02 18:13:30
|
Revision: 146 http://python-control.svn.sourceforge.net/python-control/?rev=146&view=rev Author: murrayrm Date: 2011-04-02 18:13:23 +0000 (Sat, 02 Apr 2011) Log Message: ----------- Changes in preparation for using v0.4 as trunk: * Added missing MATLAB compatible functions for module reduction * Turned off print statements in unit tests (turn back on using verbose flag) * Removed bug warnings from conversion routines; posted open bug to Pending * Updated examples/pvtol-nested-ss to reflect interface changes * See ChangeLog for a detailed list of changes Modified Paths: -------------- branches/control-0.4b/ChangeLog branches/control-0.4b/Pending branches/control-0.4b/examples/pvtol-nested-ss.py branches/control-0.4b/src/__init__.py branches/control-0.4b/src/matlab.py branches/control-0.4b/src/modelsimp.py branches/control-0.4b/src/statesp.py branches/control-0.4b/src/xferfcn.py branches/control-0.4b/tests/convert_test.py branches/control-0.4b/tests/modelsimp_test.py branches/control-0.4b/tests/slycot_convert_test.py branches/control-0.4b/tests/test_all.py Modified: branches/control-0.4b/ChangeLog =================================================================== --- branches/control-0.4b/ChangeLog 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/ChangeLog 2011-04-02 18:13:23 UTC (rev 146) @@ -1,5 +1,31 @@ 2011-04-02 Richard Murray <murray@malabar.local> + * src/__init__.py: removed import of tests module (moved to tests/) + + * src/matlab.py: Added hsvd, balred, modred to list of functions + that are imported for use as is. Updated documentation string to + indicate that these are implemented, along with a few other + functions (zero, lqr) that weren't properly listed. + + * src/modelsimp.py (balred): Removed extraneous print statements + (modred): Set method to be 'matchdc' by default (to match MATLAB) + + * src/__init__.py: added missing import of modelsimp functions + + * tests/slycot_convert_test.py (TestSlycot.testTF): turned off print + statements in unit test to make it easier to see results. Use + verbose=True to turn back on. + + * tests/convert_test.py (TestConvert.testConvert): got rid of print + statements in unittest; clutters the output so that you can't see + the errors clearly. Use verbose=True to turn back on. + + * src/statesp.py (_convertToStateSpace): removed "buggy" print + statements + + * src/xferfcn.py (_convertToTransferFunction): removed "buggy" print + statements + * tests/nichols_test.py (TestStateSpace.testNgrid): updated testcode to turn off grid in initial Nichols chart plot. Modified: branches/control-0.4b/Pending =================================================================== --- branches/control-0.4b/Pending 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/Pending 2011-04-02 18:13:23 UTC (rev 146) @@ -10,9 +10,9 @@ OPEN BUGS * matlab.step() doesn't handle systems with a pole at the origin (use lsim2) - * tests/convert_test.py not working yet - * tests/test_all.py should report on failed tests - * tests/freqresp.py needs to be converted to unit test + * TF <-> SS transformations are buggy; see tests/convert_test.py + * hsvd returns different value than MATLAB (2010a); see modelsimp_test.py + * MIMO common denominator fails unit test; see convert_test.py Transfer code from Roberto Bucher's yottalab to python-control acker - pole placement using Ackermann method @@ -40,6 +40,9 @@ * Put together unit tests for all functions (after deciding on framework) * Figure out how to import 'figure' command properly (version issue?) * Figure out source of BadCoefficients warning messages (pvtol-lqr and others) + * tests/test_all.py should report on failed tests + * tests/freqresp.py needs to be converted to unit test + * Convert examples/test-{response,statefbk}.py to unit tests TransferFunction class fixes * evalfr is not working (num, den stored as ndarrays, not poly1ds) Modified: branches/control-0.4b/examples/pvtol-nested-ss.py =================================================================== --- branches/control-0.4b/examples/pvtol-nested-ss.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/examples/pvtol-nested-ss.py 2011-04-02 18:13:23 UTC (rev 146) @@ -95,21 +95,21 @@ # (gm, pm, wgc, wpc) = margin(L); figure(6); clf; subplot(221); -(magh, phaseh) = bode(L); +bode(L, logspace(-4, 3)); # Add crossover line -subplot(magh); hold(True); -loglog([10^-4, 10^3], [1, 1], 'k-') +subplot(211); hold(True); +loglog([1e-4, 1e3], [1, 1], 'k-') # Replot phase starting at -90 degrees -bode(L, logspace(-4, 3)); (mag, phase, w) = freqresp(L, logspace(-4, 3)); phase = phase - 360; -subplot(phaseh); -semilogx([10^-4, 10^3], [-180, -180], 'k-') + +subplot(212); +semilogx([1e-4, 1e3], [-180, -180], 'k-') hold(True); semilogx(w, np.squeeze(phase), 'b-') -axis([10^-4, 10^3, -360, 0]); +axis([1e-4, 1e3, -360, 0]); xlabel('Frequency [deg]'); ylabel('Phase [deg]'); # set(gca, 'YTick', [-360, -270, -180, -90, 0]); # set(gca, 'XTick', [10^-4, 10^-2, 1, 100]); @@ -152,8 +152,8 @@ #TODO: PZmap for statespace systems has not yet been implemented. figure(10); clf(); -#(P, Z) = pzmap(T, Plot=True) -#print "Closed loop poles and zeros: ", P, Z +# (P, Z) = pzmap(T, Plot=True) +# print "Closed loop poles and zeros: ", P, Z # Gang of Four figure(11); clf(); Modified: branches/control-0.4b/src/__init__.py =================================================================== --- branches/control-0.4b/src/__init__.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/src/__init__.py 2011-04-02 18:13:23 UTC (rev 146) @@ -63,5 +63,4 @@ from bdalg import * from statefbk import * from delay import * - -from test import * +from modelsimp import * Modified: branches/control-0.4b/src/matlab.py =================================================================== --- branches/control-0.4b/src/matlab.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/src/matlab.py 2011-04-02 18:13:23 UTC (rev 146) @@ -80,6 +80,7 @@ from pzmap import pzmap from statefbk import ctrb, obsv, gram, place, lqr from delay import pade +from modelsimp import hsvd, balred, modred __doc__ += """ The control.matlab module defines functions that are roughly the @@ -139,9 +140,9 @@ lti/bandwidth - system bandwidth lti/norm - h2 and Hinfinity norms of LTI models \* lti/pole - system poles - lti/zero - system (transmission) zeros +\* lti/zero - system (transmission) zeros lti/order - model order (number of states) -\* pzmap - pole-zero map +\* pzmap - pole-zero map (TF only) lti/iopzmap - input/output pole-zero map damp - natural frequency and damping of system poles esort - sort continuous poles by real part @@ -173,19 +174,20 @@ Model simplification minreal - minimal realization and pole/zero cancellation ss/sminreal - structurally minimal realization (state space) - lti/hsvd - hankel singular values (state contributions) - lti/balred - reduced-order approximations of LTI models - ss/modred - model order reduction +\* lti/hsvd - hankel singular values (state contributions) +\* lti/balred - reduced-order approximations of LTI models +\* ss/modred - model order reduction Compensator design rlocus - evans root locus - place - pole placement +\* place - pole placement estim - form estimator given estimator gain reg - form regulator given state-feedback and estimator gains LQR/LQG design ss/lqg - single-step LQG design - lqr, dlqr - linear-Quadratic (LQ) state-feedback regulator +\* lqr - linear-Quadratic (LQ) state-feedback regulator +\* dlqr - discrete-time LQ state-feedback regulator lqry - lq regulator with output weighting lqrd - discrete LQ regulator for continuous plant ss/lqi - linear-Quadratic-Integral (LQI) controller Modified: branches/control-0.4b/src/modelsimp.py =================================================================== --- branches/control-0.4b/src/modelsimp.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/src/modelsimp.py 2011-04-02 18:13:23 UTC (rev 146) @@ -84,7 +84,7 @@ # Return the Hankel singular values return hsv -def modred(sys,ELIM,method): +def modred(sys,ELIM,method='matchdc'): """Model reduction of sys by eliminating the states in ELIM using a given method Parameters @@ -207,8 +207,8 @@ #Check system is stable D,V = np.linalg.eig(sys.A) - print D.shape - print D + # print D.shape + # print D for e in D: if e.real >= 0: raise ValueError, "Oops, the system is unstable!" Modified: branches/control-0.4b/src/statesp.py =================================================================== --- branches/control-0.4b/src/statesp.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/src/statesp.py 2011-04-02 18:13:23 UTC (rev 146) @@ -469,8 +469,6 @@ # Repeat the common denominator along the rows. den = array([den for i in range(sys.outputs)]) # TODO: transfer function to state space conversion is still buggy! - print "Warning: transfer function to state space conversion by td04ad \ -is still buggy! Advise converting state space sys back to tf to verify the transformation was correct." #print num #print shape(num) ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=0.0) Modified: branches/control-0.4b/src/xferfcn.py =================================================================== --- branches/control-0.4b/src/xferfcn.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/src/xferfcn.py 2011-04-02 18:13:23 UTC (rev 146) @@ -713,8 +713,6 @@ # Use Slycot to make the transformation. TODO: this is still somewhat # buggy! - print "Warning: state space to transfer function conversion by tb04ad \ -is still buggy!" tfout = tb04ad(sys.states, sys.inputs, sys.outputs, sys.A, sys.B, sys.C, sys.D,tol1=0.0) @@ -727,8 +725,8 @@ num[i][j] = list(tfout[6][i, j, :]) # Each transfer function matrix row has a common denominator. den[i][j] = list(tfout[5][i, :]) - print num - print den + # print num + # print den return TransferFunction(num, den) elif isinstance(sys, (int, long, float, complex)): if "inputs" in kw: Modified: branches/control-0.4b/tests/convert_test.py =================================================================== --- branches/control-0.4b/tests/convert_test.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/tests/convert_test.py 2011-04-02 18:13:23 UTC (rev 146) @@ -1,6 +1,6 @@ #!/usr/bin/env python -"""TestConvert.py +"""convert_test.py Test state space and transfer function conversion. @@ -40,7 +40,7 @@ print "sys%i:\n" % ind print sys - def testConvert(self): + def testConvert(self, verbose=0): """Test state space to transfer function conversion.""" #Currently it only tests that a TF->SS->TF generates an unchanged TF @@ -52,16 +52,20 @@ #start with a random SS system and transform to TF #then back to SS, check that the matrices are the same. ssOriginal = matlab.rss(states, inputs, outputs) - self.printSys(ssOriginal, 1) + if (verbose): + self.printSys(ssOriginal, 1) tfOriginal = matlab.tf(ssOriginal) - self.printSys(tfOriginal, 2) + if (verbose): + self.printSys(tfOriginal, 2) ssTransformed = matlab.ss(tfOriginal) - self.printSys(ssTransformed, 3) + if (verbose): + self.printSys(ssTransformed, 3) tfTransformed = matlab.tf(ssTransformed) - self.printSys(tfTransformed, 4) + if (verbose): + self.printSys(tfTransformed, 4) for inputNum in range(inputs): for outputNum in range(outputs): @@ -81,8 +85,6 @@ #[mag,phase,freq]=bode(sys) #it doesn't seem to...... #This should be added. - - def suite(): return unittest.TestLoader().loadTestsFromTestCase(TestConvert) Modified: branches/control-0.4b/tests/modelsimp_test.py =================================================================== --- branches/control-0.4b/tests/modelsimp_test.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/tests/modelsimp_test.py 2011-04-02 18:13:23 UTC (rev 146) @@ -16,7 +16,7 @@ D = np.matrix("9.") sys = ss(A,B,C,D) hsv = hsvd(sys) - hsvtrue = np.matrix("24.42686 0.5731395") + hsvtrue = np.matrix("24.42686 0.5731395") # from MATLAB np.testing.assert_array_almost_equal(hsv, hsvtrue) def testMarkov(self): Modified: branches/control-0.4b/tests/slycot_convert_test.py =================================================================== --- branches/control-0.4b/tests/slycot_convert_test.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/tests/slycot_convert_test.py 2011-04-02 18:13:23 UTC (rev 146) @@ -26,24 +26,23 @@ self.maxI = 1 self.maxO = 1 - def testTF(self): + def testTF(self, verbose=False): """ Directly tests the functions tb04ad and td04ad through direct comparison of transfer function coefficients. - Similar to TestConvert, but tests at a lower level. + Similar to convert_test, but tests at a lower level. """ for states in range(1, self.maxStates): for inputs in range(1, self.maxI+1): for outputs in range(1, self.maxO+1): for testNum in range(self.numTests): - ssOriginal = matlab.rss(states, inputs, outputs) + if (verbose): + print '====== Original SS ==========' + print ssOriginal + print 'states=',states + print 'inputs=',inputs + print 'outputs=',outputs - print '====== Original SS ==========' - print ssOriginal - print 'states=',states - print 'inputs=',inputs - print 'outputs=',outputs - tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\ tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad(states,inputs,outputs,\ ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0) @@ -55,15 +54,16 @@ tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad(ssTransformed_nr,\ inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,ssTransformed_D,tol1=0.0) #print 'size(Trans_A)=',ssTransformed_A.shape - print '===== Transformed SS ==========' - print matlab.ss(ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D) - #print 'Trans_nr=',ssTransformed_nr - #print 'tfOrig_index=',tfOriginal_index - #print 'tfOrig_ucoeff=',tfOriginal_ucoeff - #print 'tfOrig_dcoeff=',tfOriginal_dcoeff - #print 'tfTrans_index=',tfTransformed_index - #print 'tfTrans_ucoeff=',tfTransformed_ucoeff - #print 'tfTrans_dcoeff=',tfTransformed_dcoeff + if (verbose): + print '===== Transformed SS ==========' + print matlab.ss(ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D) + # print 'Trans_nr=',ssTransformed_nr + # print 'tfOrig_index=',tfOriginal_index + # print 'tfOrig_ucoeff=',tfOriginal_ucoeff + # print 'tfOrig_dcoeff=',tfOriginal_dcoeff + # print 'tfTrans_index=',tfTransformed_index + # print 'tfTrans_ucoeff=',tfTransformed_ucoeff + # print 'tfTrans_dcoeff=',tfTransformed_dcoeff #Compare the TF directly, must match #numerators np.testing.assert_array_almost_equal(tfOriginal_ucoeff,tfTransformed_ucoeff,decimal=3) Modified: branches/control-0.4b/tests/test_all.py =================================================================== --- branches/control-0.4b/tests/test_all.py 2011-04-02 16:15:41 UTC (rev 145) +++ branches/control-0.4b/tests/test_all.py 2011-04-02 18:13:23 UTC (rev 146) @@ -7,7 +7,7 @@ import re # regular expressions import os # operating system commands -def test_all(verbosity=2): +def test_all(verbosity=0): """ Runs all tests written for python-control. """ try: # autodiscovery (python 2.7+) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-04-02 16:15:47
|
Revision: 145 http://python-control.svn.sourceforge.net/python-control/?rev=145&view=rev Author: murrayrm Date: 2011-04-02 16:15:41 +0000 (Sat, 02 Apr 2011) Log Message: ----------- updated Nichols chart code to match 0.3d Modified Paths: -------------- branches/control-0.4b/ChangeLog branches/control-0.4b/src/freqplot.py branches/control-0.4b/src/matlab.py branches/control-0.4b/src/nichols.py branches/control-0.4b/tests/nichols_test.py Modified: branches/control-0.4b/ChangeLog =================================================================== --- branches/control-0.4b/ChangeLog 2011-04-02 15:52:38 UTC (rev 144) +++ branches/control-0.4b/ChangeLog 2011-04-02 16:15:41 UTC (rev 145) @@ -1,3 +1,15 @@ +2011-04-02 Richard Murray <murray@malabar.local> + + * tests/nichols_test.py (TestStateSpace.testNgrid): updated testcode + to turn off grid in initial Nichols chart plot. + + * src/freqplot.py: updated comments at top of file to reflect + nichols chart move + + * src/nichols.py: transferred over changes from v0.3d + + * src/matlab.py (ngrid): moved import to function + 2011-03-31 Richard Murray <murray@malabar.local> * examples/pvtol-nested.py: updated stability margin plot to use Modified: branches/control-0.4b/src/freqplot.py =================================================================== --- branches/control-0.4b/src/freqplot.py 2011-04-02 15:52:38 UTC (rev 144) +++ branches/control-0.4b/src/freqplot.py 2011-04-02 16:15:41 UTC (rev 145) @@ -4,7 +4,8 @@ # Date: 24 May 09 # # This file contains some standard control system plots: Bode plots, -# Nyquist plots, Nichols plots and pole-zero diagrams +# Nyquist plots and pole-zero diagrams. The code for Nichols charts +# is in nichols.py. # # Copyright (c) 2010 by California Institute of Technology # All rights reserved. Modified: branches/control-0.4b/src/matlab.py =================================================================== --- branches/control-0.4b/src/matlab.py 2011-04-02 15:52:38 UTC (rev 144) +++ branches/control-0.4b/src/matlab.py 2011-04-02 16:15:41 UTC (rev 145) @@ -75,7 +75,7 @@ # Import MATLAB-like functions that can be used as-is from ctrlutil import unwrap from freqplot import nyquist, gangof4 -from nichols import nichols, nichols_grid +from nichols import nichols from bdalg import series, parallel, negate, feedback from pzmap import pzmap from statefbk import ctrb, obsv, gram, place, lqr @@ -753,6 +753,7 @@ ===== ngrid() """ + from nichols import nichols_grid nichols_grid() # Modified: branches/control-0.4b/src/nichols.py =================================================================== --- branches/control-0.4b/src/nichols.py 2011-04-02 15:52:38 UTC (rev 144) +++ branches/control-0.4b/src/nichols.py 2011-04-02 16:15:41 UTC (rev 145) @@ -45,7 +45,8 @@ from ctrlutil import unwrap from freqplot import default_frequency_range -def nichols(syslist, omega=None): +# Nichols plot +def nichols(syslist, omega=None, grid=True): """Nichols plot for a system Usage @@ -60,6 +61,8 @@ List of linear input/output systems (single system is OK) omega : freq_range Range of frequencies (list or bounds) in rad/sec + grid : boolean, optional + True if the plot should include a Nichols-chart grid. Default is True. Return values ------------- @@ -71,7 +74,7 @@ syslist = (syslist,) # Select a default range if none is provided - if (omega == None): + if omega is None: omega = default_frequency_range(syslist) for sys in syslist: @@ -94,88 +97,110 @@ # Mark the -180 point plt.plot([-180], [0], 'r+') + + # Add grid + if grid: + nichols_grid() # Nichols grid -def nichols_grid(): +#! TODO: Consider making linestyle configurable +def nichols_grid(cl_mags=None, cl_phases=None): """Nichols chart grid Usage ===== nichols_grid() - Plots a Nichols chart grid on the current axis. + Plots a Nichols chart grid on the current axis, or creates a new chart + if no plot already exists. Parameters ---------- - None + cl_mags : array-like (dB), optional + Array of closed-loop magnitudes defining the iso-gain lines on a + custom Nichols chart. + cl_phases : array-like (degrees), optional + Array of closed-loop phases defining the iso-phase lines on a custom + Nichols chart. Must be in the range -360 < cl_phases < 0 Return values ------------- None """ - mag_min_default = -40.0 # dB - mag_step = 20.0 # dB + # Default chart size + ol_phase_min = -359.99 + ol_phase_max = 0.0 + ol_mag_min = -40.0 + ol_mag_max = default_ol_mag_max = 50.0 + + # Find bounds of the current dataset, if there is one. + if plt.gcf().gca().has_data(): + ol_phase_min, ol_phase_max, ol_mag_min, ol_mag_max = plt.axis() - # Chart defaults - phase_min, phase_max, mag_min, mag_max = -360.0, 0.0, mag_min_default, 40.0 - - # Set actual chart bounds based on current plot - if plt.gcf().gca().has_data(): - phase_min, phase_max, mag_min, mag_max = plt.axis() - # M-circle magnitudes. - # The "fixed" set are always generated, since this guarantees a recognizable - # Nichols chart grid. - mags_fixed = np.array([-40.0, -20.0, -12.0, -6.0, -3.0, -1.0, -0.5, 0.0, - 0.25, 0.5, 1.0, 3.0, 6.0, 12.0]) - - if mag_min < mag_min_default: - # Outside the "fixed" set of magnitudes, the generated M-circles - # are extended in steps of 'mag_step' dB to cover anything made - # visible by the range of the existing plot - mags_adjust = np.arange(mag_step*np.floor(mag_min/mag_step), - mag_min_default, mag_step) - mags = np.concatenate((mags_adjust, mags_fixed)) - else: - mags = mags_fixed + if cl_mags is None: + # Default chart magnitudes + # The key set of magnitudes are always generated, since this + # guarantees a recognizable Nichols chart grid. + key_cl_mags = np.array([-40.0, -20.0, -12.0, -6.0, -3.0, -1.0, -0.5, 0.0, + 0.25, 0.5, 1.0, 3.0, 6.0, 12.0]) + # Extend the range of magnitudes if necessary. The extended arange + # will end up empty if no extension is required. Assumes that closed-loop + # magnitudes are approximately aligned with open-loop magnitudes beyond + # the value of np.min(key_cl_mags) + cl_mag_step = -20.0 # dB + extended_cl_mags = np.arange(np.min(key_cl_mags), + ol_mag_min + cl_mag_step, cl_mag_step) + cl_mags = np.concatenate((extended_cl_mags, key_cl_mags)) # N-circle phases (should be in the range -360 to 0) - phases = np.array([-0.25, -10.0, -20.0, -30.0, -45.0, -60.0, -90.0, - -120.0, -150.0, -180.0, -210.0, -240.0, -270.0, - -310.0, -325.0, -340.0, -350.0, -359.75]) + if cl_phases is None: + # Choose a reasonable set of default phases (denser if the open-loop + # data is restricted to a relatively small range of phases). + key_cl_phases = np.array([-0.25, -45.0, -90.0, -180.0, -270.0, -325.0, -359.75]) + if np.abs(ol_phase_max - ol_phase_min) < 90.0: + other_cl_phases = np.arange(-10.0, -360.0, -10.0) + else: + other_cl_phases = np.arange(-10.0, -360.0, -20.0) + cl_phases = np.concatenate((key_cl_phases, other_cl_phases)) + else: + assert ((-360.0 < np.min(cl_phases)) and (np.max(cl_phases) < 0.0)) # Find the M-contours - m = m_circles(mags, phase_min=np.min(phases), phase_max=np.max(phases)) + m = m_circles(cl_mags, phase_min=np.min(cl_phases), phase_max=np.max(cl_phases)) m_mag = 20*sp.log10(np.abs(m)) m_phase = sp.mod(sp.degrees(sp.angle(m)), -360.0) # Unwrap # Find the N-contours - n = n_circles(phases, mag_min=np.min(mags), mag_max=np.max(mags)) + n = n_circles(cl_phases, mag_min=np.min(cl_mags), mag_max=np.max(cl_mags)) n_mag = 20*sp.log10(np.abs(n)) n_phase = sp.mod(sp.degrees(sp.angle(n)), -360.0) # Unwrap # Plot the contours behind other plot elements. # The "phase offset" is used to produce copies of the chart that cover # the entire range of the plotted data, starting from a base chart computed - # over the range -360 < phase < 0 (see above). Given the range + # over the range -360 < phase < 0. Given the range # the base chart is computed over, the phase offset should be 0 - # for -360 < phase_min < 0. - phase_offset_min = 360.0*np.ceil(phase_min/360.0) - phase_offset_max = 360.0*np.ceil(phase_max/360.0) + 360.0 + # for -360 < ol_phase_min < 0. + phase_offset_min = 360.0*np.ceil(ol_phase_min/360.0) + phase_offset_max = 360.0*np.ceil(ol_phase_max/360.0) + 360.0 phase_offsets = np.arange(phase_offset_min, phase_offset_max, 360.0) + for phase_offset in phase_offsets: + # Draw M and N contours plt.plot(m_phase + phase_offset, m_mag, color='gray', - linestyle='dashed', zorder=0) + linestyle='dotted', zorder=0) plt.plot(n_phase + phase_offset, n_mag, color='gray', - linestyle='dashed', zorder=0) + linestyle='dotted', zorder=0) - # Add magnitude labels - for x, y, m in zip(m_phase[:][-1], m_mag[:][-1], mags): - align = 'right' if m < 0.0 else 'left' - plt.text(x, y, str(m) + ' dB', size='small', ha=align) + # Add magnitude labels + for x, y, m in zip(m_phase[:][-1] + phase_offset, m_mag[:][-1], cl_mags): + align = 'right' if m < 0.0 else 'left' + plt.text(x, y, str(m) + ' dB', size='small', ha=align, color='gray') - # Make sure axes conform to any pre-existing plot. - plt.axis([phase_min, phase_max, mag_min, mag_max]) + # Fit axes to generated chart + plt.axis([phase_offset_min - 360.0, phase_offset_max - 360.0, + np.min(cl_mags), np.max([ol_mag_max, default_ol_mag_max])]) # # Utility functions @@ -185,38 +210,44 @@ # # Compute contours of a closed-loop transfer function -def closed_loop_contours(Hmag, Hphase): - """Contours of the function H = G/(1+G). +def closed_loop_contours(Gcl_mags, Gcl_phases): + """Contours of the function Gcl = Gol/(1+Gol), where + Gol is an open-loop transfer function, and Gcl is a corresponding + closed-loop transfer function. Usage ===== - contours = closed_loop_contours(mags, phases) + contours = closed_loop_contours(Gcl_mags, Gcl_phases) Parameters ---------- - mags : array-like - Meshgrid array of magnitudes of the contours - phases : array-like - Meshgrid array of phases in radians of the contours + Gcl_mags : array-like + Array of magnitudes of the contours + Gcl_phases : array-like + Array of phases in radians of the contours Return values ------------- contours : complex array Array of complex numbers corresponding to the contours. """ - # Compute the contours in H-space - H = Hmag*sp.exp(1.j*Hphase) + # Compute the contours in Gcl-space. Since we're given closed-loop + # magnitudes and phases, this is just a case of converting them into + # a complex number. + Gcl = Gcl_mags*sp.exp(1.j*Gcl_phases) - # Invert H = G/(1+G) to get an expression for the contours in G-space - return H/(1.0 - H) + # Invert Gcl = Gol/(1+Gol) to map the contours into the open-loop space + return Gcl/(1.0 - Gcl) # M-circle def m_circles(mags, phase_min=-359.75, phase_max=-0.25): - """Constant-magnitude contours of the function H = G/(1+G). + """Constant-magnitude contours of the function Gcl = Gol/(1+Gol), where + Gol is an open-loop transfer function, and Gcl is a corresponding + closed-loop transfer function. Usage ===== - contours = m_circles(mags) + contours = m_circles(mags, phase_min, phase_max) Parameters ---------- @@ -234,17 +265,19 @@ """ # Convert magnitudes and phase range into a grid suitable for # building contours - phases = sp.radians(sp.linspace(phase_min, phase_max, 500)) - Hmag, Hphase = sp.meshgrid(10.0**(mags/20.0), phases) - return closed_loop_contours(Hmag, Hphase) + phases = sp.radians(sp.linspace(phase_min, phase_max, 2000)) + Gcl_mags, Gcl_phases = sp.meshgrid(10.0**(mags/20.0), phases) + return closed_loop_contours(Gcl_mags, Gcl_phases) # N-circle def n_circles(phases, mag_min=-40.0, mag_max=12.0): - """Constant-phase contours of the function H = G/(1+G). + """Constant-phase contours of the function Gcl = Gol/(1+Gol), where + Gol is an open-loop transfer function, and Gcl is a corresponding + closed-loop transfer function. Usage ===== - contour = n_circles(angles) + contours = n_circles(phases, mag_min, mag_max) Parameters ---------- @@ -263,5 +296,5 @@ # Convert phases and magnitude range into a grid suitable for # building contours mags = sp.linspace(10**(mag_min/20.0), 10**(mag_max/20.0), 2000) - Hphase, Hmag = sp.meshgrid(sp.radians(phases), mags) - return closed_loop_contours(Hmag, Hphase) + Gcl_phases, Gcl_mags = sp.meshgrid(sp.radians(phases), mags) + return closed_loop_contours(Gcl_mags, Gcl_phases) Modified: branches/control-0.4b/tests/nichols_test.py =================================================================== --- branches/control-0.4b/tests/nichols_test.py 2011-04-02 15:52:38 UTC (rev 144) +++ branches/control-0.4b/tests/nichols_test.py 2011-04-02 16:15:41 UTC (rev 145) @@ -20,13 +20,13 @@ self.sys = StateSpace(A, B, C, D) - def testNichols(self): + def testNicholsPlain(self): """Generate a Nichols plot.""" nichols(self.sys) def testNgrid(self): """Generate a Nichols plot.""" - nichols(self.sys) + nichols(self.sys, grid=False) ngrid() def suite(): This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-04-02 15:52:46
|
Revision: 144 http://python-control.svn.sourceforge.net/python-control/?rev=144&view=rev Author: murrayrm Date: 2011-04-02 15:52:38 +0000 (Sat, 02 Apr 2011) Log Message: ----------- * Moved nichols() code to separate file, added unit test; add'l updates pending * Moved slycot imports inside functions that need them * Updated pvtol-nested example to work with v0.4 * See ChangeLog for detailed listing of changes Modified Paths: -------------- branches/control-0.4b/ChangeLog branches/control-0.4b/examples/pvtol-nested.py branches/control-0.4b/src/__init__.py branches/control-0.4b/src/freqplot.py branches/control-0.4b/src/matlab.py branches/control-0.4b/src/statesp.py branches/control-0.4b/src/xferfcn.py Added Paths: ----------- branches/control-0.4b/src/nichols.py branches/control-0.4b/tests/nichols_test.py Modified: branches/control-0.4b/ChangeLog =================================================================== --- branches/control-0.4b/ChangeLog 2011-04-01 04:46:30 UTC (rev 143) +++ branches/control-0.4b/ChangeLog 2011-04-02 15:52:38 UTC (rev 144) @@ -1,3 +1,27 @@ +2011-03-31 Richard Murray <murray@malabar.local> + + * examples/pvtol-nested.py: updated stability margin plot to use + proper calling format for bode(). + + * src/statesp.py (_convertToStateSpace): moved slycot import + to the location where it is actually needed (allows running some + commands without slycot installed) + + * src/xferfcn.py (_convertToTransferFunction): moved slycot import + to the location where it is actually needed (allows running some + commands without slycot installed) + + * src/nichols.py: new file for Nichols plot routines; move + nichols(), nichols_grid(), closed_loop_contours(), m_circles(), + n_circles() + + * src/__init__.py, src/freqresp.py, src/matlab.py: updated to match + new file structure for Nichols charts + + * src/nichols.py (nichols): updated processing of freqresp to take + into account the fact that return arguments are now a matrix of + results (even for a SISO system) + 2011-03-30 Richard Murray <murray@malabar.local> * tests/: added top level subdirectory, to be used for unit tests. Modified: branches/control-0.4b/examples/pvtol-nested.py =================================================================== --- branches/control-0.4b/examples/pvtol-nested.py 2011-04-01 04:46:30 UTC (rev 143) +++ branches/control-0.4b/examples/pvtol-nested.py 2011-04-02 15:52:38 UTC (rev 144) @@ -84,22 +84,22 @@ # (gm, pm, wgc, wpc) = margin(L); #! TODO: this figure has something wrong; axis limits mismatch -figure(6); clf; subplot(221); -(magh, phaseh) = bode(L); +figure(6); clf; +bode(L); # Add crossover line -subplot(magh); hold(True); -loglog([10^-4, 10^3], [1, 1], 'k-') +subplot(211); hold(True); +loglog([1e-4, 1e3], [1, 1], 'k-') # Replot phase starting at -90 degrees bode(L, logspace(-4, 3)); (mag, phase, w) = freqresp(L, logspace(-4, 3)); phase = phase - 360; -subplot(phaseh); -semilogx([10^-4, 10^3], [-180, -180], 'k-') +subplot(212); +semilogx([1e-4, 1e3], [-180, -180], 'k-') hold(True); semilogx(w, np.squeeze(phase), 'b-') -axis([10^-4, 10^3, -360, 0]); +axis([1e-4, 1e3, -360, 0]); xlabel('Frequency [deg]'); ylabel('Phase [deg]'); # set(gca, 'YTick', [-360, -270, -180, -90, 0]); # set(gca, 'XTick', [10^-4, 10^-2, 1, 100]); Modified: branches/control-0.4b/src/__init__.py =================================================================== --- branches/control-0.4b/src/__init__.py 2011-04-01 04:46:30 UTC (rev 143) +++ branches/control-0.4b/src/__init__.py 2011-04-02 15:52:38 UTC (rev 144) @@ -59,6 +59,7 @@ from xferfcn import * from statesp import * from freqplot import * +from nichols import * from bdalg import * from statefbk import * from delay import * Modified: branches/control-0.4b/src/freqplot.py =================================================================== --- branches/control-0.4b/src/freqplot.py 2011-04-01 04:46:30 UTC (rev 143) +++ branches/control-0.4b/src/freqplot.py 2011-04-02 15:52:38 UTC (rev 144) @@ -212,139 +212,7 @@ plt.plot([-1], [0], 'r+') return x, y, omega -# Nichols plot -# Contributed by Allan McInnes <All...@ca...> -#! TODO: need unit test code -def nichols(syslist, omega=None): - """Nichols plot for a system - Usage - ===== - magh = nichols(sys, omega=None) - - Plots a Nichols plot for the system over a (optional) frequency range. - - Parameters - ---------- - syslist : linsys - List of linear input/output systems (single system is OK) - omega : freq_range - Range of frequencies (list or bounds) in rad/sec - - Return values - ------------- - None - """ - - # If argument was a singleton, turn it into a list - if (not getattr(syslist, '__iter__', False)): - syslist = (syslist,) - - # Select a default range if none is provided - if (omega == None): - omega = default_frequency_range(syslist) - - for sys in syslist: - # Get the magnitude and phase of the system - mag, phase, omega = sys.freqresp(omega) - - # Convert to Nichols-plot format (phase in degrees, - # and magnitude in dB) - x = unwrap(sp.degrees(phase), 360) - y = 20*sp.log10(mag) - - # Generate the plot - plt.plot(x, y) - - plt.xlabel('Phase (deg)') - plt.ylabel('Magnitude (dB)') - plt.title('Nichols Plot') - - # Mark the -180 point - plt.plot([-180], [0], 'r+') - -# Nichols grid -def nichols_grid(): - """Nichols chart grid - - Usage - ===== - nichols_grid() - - Plots a Nichols chart grid on the current axis. - - Parameters - ---------- - None - - Return values - ------------- - None - """ - mag_min_default = -40.0 # dB - mag_step = 20.0 # dB - - # Chart defaults - phase_min, phase_max, mag_min, mag_max = -360.0, 0.0, mag_min_default, 40.0 - - # Set actual chart bounds based on current plot - if plt.gcf().gca().has_data(): - phase_min, phase_max, mag_min, mag_max = plt.axis() - - # M-circle magnitudes. - # The "fixed" set are always generated, since this guarantees a recognizable - # Nichols chart grid. - mags_fixed = np.array([-40.0, -20.0, -12.0, -6.0, -3.0, -1.0, -0.5, 0.0, - 0.25, 0.5, 1.0, 3.0, 6.0, 12.0]) - - if mag_min < mag_min_default: - # Outside the "fixed" set of magnitudes, the generated M-circles - # are extended in steps of 'mag_step' dB to cover anything made - # visible by the range of the existing plot - mags_adjust = np.arange(mag_step*np.floor(mag_min/mag_step), - mag_min_default, mag_step) - mags = np.concatenate((mags_adjust, mags_fixed)) - else: - mags = mags_fixed - - # N-circle phases (should be in the range -360 to 0) - phases = np.array([-0.25, -10.0, -20.0, -30.0, -45.0, -60.0, -90.0, - -120.0, -150.0, -180.0, -210.0, -240.0, -270.0, - -310.0, -325.0, -340.0, -350.0, -359.75]) - - # Find the M-contours - m = m_circles(mags, phase_min=np.min(phases), phase_max=np.max(phases)) - m_mag = 20*sp.log10(np.abs(m)) - m_phase = sp.mod(sp.degrees(sp.angle(m)), -360.0) # Unwrap - - # Find the N-contours - n = n_circles(phases, mag_min=np.min(mags), mag_max=np.max(mags)) - n_mag = 20*sp.log10(np.abs(n)) - n_phase = sp.mod(sp.degrees(sp.angle(n)), -360.0) # Unwrap - - # Plot the contours behind other plot elements. - # The "phase offset" is used to produce copies of the chart that cover - # the entire range of the plotted data, starting from a base chart computed - # over the range -360 < phase < 0 (see above). Given the range - # the base chart is computed over, the phase offset should be 0 - # for -360 < phase_min < 0. - phase_offset_min = 360.0*np.ceil(phase_min/360.0) - phase_offset_max = 360.0*np.ceil(phase_max/360.0) + 360.0 - phase_offsets = np.arange(phase_offset_min, phase_offset_max, 360.0) - for phase_offset in phase_offsets: - plt.plot(m_phase + phase_offset, m_mag, color='gray', - linestyle='dashed', zorder=0) - plt.plot(n_phase + phase_offset, n_mag, color='gray', - linestyle='dashed', zorder=0) - - # Add magnitude labels - for x, y, m in zip(m_phase[:][-1], m_mag[:][-1], mags): - align = 'right' if m < 0.0 else 'left' - plt.text(x, y, str(m) + ' dB', size='small', ha=align) - - # Make sure axes conform to any pre-existing plot. - plt.axis([phase_min, phase_max, mag_min, mag_max]) - # Gang of Four #! TODO: think about how (and whether) to handle lists of systems def gangof4(P, C, omega=None): @@ -462,84 +330,3 @@ return omega -# Compute contours of a closed-loop transfer function -def closed_loop_contours(Hmag, Hphase): - """Contours of the function H = G/(1+G). - - Usage - ===== - contours = closed_loop_contours(mags, phases) - - Parameters - ---------- - mags : array-like - Meshgrid array of magnitudes of the contours - phases : array-like - Meshgrid array of phases in radians of the contours - - Return values - ------------- - contours : complex array - Array of complex numbers corresponding to the contours. - """ - # Compute the contours in H-space - H = Hmag*sp.exp(1.j*Hphase) - - # Invert H = G/(1+G) to get an expression for the contours in G-space - return H/(1.0 - H) - -# M-circle -def m_circles(mags, phase_min=-359.75, phase_max=-0.25): - """Constant-magnitude contours of the function H = G/(1+G). - - Usage - ===== - contours = m_circles(mags) - - Parameters - ---------- - mags : array-like - Array of magnitudes in dB of the M-circles - phase_min : degrees - Minimum phase in degrees of the N-circles - phase_max : degrees - Maximum phase in degrees of the N-circles - - Return values - ------------- - contours : complex array - Array of complex numbers corresponding to the contours. - """ - # Convert magnitudes and phase range into a grid suitable for - # building contours - phases = sp.radians(sp.linspace(phase_min, phase_max, 500)) - Hmag, Hphase = sp.meshgrid(10.0**(mags/20.0), phases) - return closed_loop_contours(Hmag, Hphase) - -# N-circle -def n_circles(phases, mag_min=-40.0, mag_max=12.0): - """Constant-phase contours of the function H = G/(1+G). - - Usage - ===== - contour = n_circles(angles) - - Parameters - ---------- - phases : array-like - Array of phases in degrees of the N-circles - mag_min : dB - Minimum magnitude in dB of the N-circles - mag_max : dB - Maximum magnitude in dB of the N-circles - - Return values - ------------- - contours : complex array - Array of complex numbers corresponding to the contours. - """ - # Convert phases and magnitude range into a grid suitable for - # building contours - mags = sp.linspace(10**(mag_min/20.0), 10**(mag_max/20.0), 2000) - Hphase, Hmag = sp.meshgrid(sp.radians(phases), mags) - return closed_loop_contours(Hmag, Hphase) Modified: branches/control-0.4b/src/matlab.py =================================================================== --- branches/control-0.4b/src/matlab.py 2011-04-01 04:46:30 UTC (rev 143) +++ branches/control-0.4b/src/matlab.py 2011-04-02 15:52:38 UTC (rev 144) @@ -74,7 +74,8 @@ # Import MATLAB-like functions that can be used as-is from ctrlutil import unwrap -from freqplot import nyquist, nichols, gangof4 +from freqplot import nyquist, gangof4 +from nichols import nichols, nichols_grid from bdalg import series, parallel, negate, feedback from pzmap import pzmap from statefbk import ctrb, obsv, gram, place, lqr @@ -752,7 +753,7 @@ ===== ngrid() """ - freqplot.nichols_grid() + nichols_grid() # # Modifications to scipy.signal functions Added: branches/control-0.4b/src/nichols.py =================================================================== --- branches/control-0.4b/src/nichols.py (rev 0) +++ branches/control-0.4b/src/nichols.py 2011-04-02 15:52:38 UTC (rev 144) @@ -0,0 +1,267 @@ +# nichols.py - Nichols plot +# +# Contributed by Allan McInnes <All...@ca...> +# +# This file contains some standard control system plots: Bode plots, +# Nyquist plots, Nichols plots and pole-zero diagrams +# +# Copyright (c) 2010 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. +# +# $Id: freqplot.py 139 2011-03-30 16:19:59Z murrayrm $ + +import matplotlib.pyplot as plt +import scipy as sp +import numpy as np +from ctrlutil import unwrap +from freqplot import default_frequency_range + +def nichols(syslist, omega=None): + """Nichols plot for a system + + Usage + ===== + magh = nichols(sys, omega=None) + + Plots a Nichols plot for the system over a (optional) frequency range. + + Parameters + ---------- + syslist : linsys + List of linear input/output systems (single system is OK) + omega : freq_range + Range of frequencies (list or bounds) in rad/sec + + Return values + ------------- + None + """ + + # If argument was a singleton, turn it into a list + if (not getattr(syslist, '__iter__', False)): + syslist = (syslist,) + + # Select a default range if none is provided + if (omega == None): + omega = default_frequency_range(syslist) + + for sys in syslist: + # Get the magnitude and phase of the system + mag_tmp, phase_tmp, omega = sys.freqresp(omega) + mag = np.squeeze(mag_tmp) + phase = np.squeeze(phase_tmp) + + # Convert to Nichols-plot format (phase in degrees, + # and magnitude in dB) + x = unwrap(sp.degrees(phase), 360) + y = 20*sp.log10(mag) + + # Generate the plot + plt.plot(x, y) + + plt.xlabel('Phase (deg)') + plt.ylabel('Magnitude (dB)') + plt.title('Nichols Plot') + + # Mark the -180 point + plt.plot([-180], [0], 'r+') + +# Nichols grid +def nichols_grid(): + """Nichols chart grid + + Usage + ===== + nichols_grid() + + Plots a Nichols chart grid on the current axis. + + Parameters + ---------- + None + + Return values + ------------- + None + """ + mag_min_default = -40.0 # dB + mag_step = 20.0 # dB + + # Chart defaults + phase_min, phase_max, mag_min, mag_max = -360.0, 0.0, mag_min_default, 40.0 + + # Set actual chart bounds based on current plot + if plt.gcf().gca().has_data(): + phase_min, phase_max, mag_min, mag_max = plt.axis() + + # M-circle magnitudes. + # The "fixed" set are always generated, since this guarantees a recognizable + # Nichols chart grid. + mags_fixed = np.array([-40.0, -20.0, -12.0, -6.0, -3.0, -1.0, -0.5, 0.0, + 0.25, 0.5, 1.0, 3.0, 6.0, 12.0]) + + if mag_min < mag_min_default: + # Outside the "fixed" set of magnitudes, the generated M-circles + # are extended in steps of 'mag_step' dB to cover anything made + # visible by the range of the existing plot + mags_adjust = np.arange(mag_step*np.floor(mag_min/mag_step), + mag_min_default, mag_step) + mags = np.concatenate((mags_adjust, mags_fixed)) + else: + mags = mags_fixed + + # N-circle phases (should be in the range -360 to 0) + phases = np.array([-0.25, -10.0, -20.0, -30.0, -45.0, -60.0, -90.0, + -120.0, -150.0, -180.0, -210.0, -240.0, -270.0, + -310.0, -325.0, -340.0, -350.0, -359.75]) + + # Find the M-contours + m = m_circles(mags, phase_min=np.min(phases), phase_max=np.max(phases)) + m_mag = 20*sp.log10(np.abs(m)) + m_phase = sp.mod(sp.degrees(sp.angle(m)), -360.0) # Unwrap + + # Find the N-contours + n = n_circles(phases, mag_min=np.min(mags), mag_max=np.max(mags)) + n_mag = 20*sp.log10(np.abs(n)) + n_phase = sp.mod(sp.degrees(sp.angle(n)), -360.0) # Unwrap + + # Plot the contours behind other plot elements. + # The "phase offset" is used to produce copies of the chart that cover + # the entire range of the plotted data, starting from a base chart computed + # over the range -360 < phase < 0 (see above). Given the range + # the base chart is computed over, the phase offset should be 0 + # for -360 < phase_min < 0. + phase_offset_min = 360.0*np.ceil(phase_min/360.0) + phase_offset_max = 360.0*np.ceil(phase_max/360.0) + 360.0 + phase_offsets = np.arange(phase_offset_min, phase_offset_max, 360.0) + for phase_offset in phase_offsets: + plt.plot(m_phase + phase_offset, m_mag, color='gray', + linestyle='dashed', zorder=0) + plt.plot(n_phase + phase_offset, n_mag, color='gray', + linestyle='dashed', zorder=0) + + # Add magnitude labels + for x, y, m in zip(m_phase[:][-1], m_mag[:][-1], mags): + align = 'right' if m < 0.0 else 'left' + plt.text(x, y, str(m) + ' dB', size='small', ha=align) + + # Make sure axes conform to any pre-existing plot. + plt.axis([phase_min, phase_max, mag_min, mag_max]) + +# +# Utility functions +# +# This section of the code contains some utility functions for +# generating Nichols plots +# + +# Compute contours of a closed-loop transfer function +def closed_loop_contours(Hmag, Hphase): + """Contours of the function H = G/(1+G). + + Usage + ===== + contours = closed_loop_contours(mags, phases) + + Parameters + ---------- + mags : array-like + Meshgrid array of magnitudes of the contours + phases : array-like + Meshgrid array of phases in radians of the contours + + Return values + ------------- + contours : complex array + Array of complex numbers corresponding to the contours. + """ + # Compute the contours in H-space + H = Hmag*sp.exp(1.j*Hphase) + + # Invert H = G/(1+G) to get an expression for the contours in G-space + return H/(1.0 - H) + +# M-circle +def m_circles(mags, phase_min=-359.75, phase_max=-0.25): + """Constant-magnitude contours of the function H = G/(1+G). + + Usage + ===== + contours = m_circles(mags) + + Parameters + ---------- + mags : array-like + Array of magnitudes in dB of the M-circles + phase_min : degrees + Minimum phase in degrees of the N-circles + phase_max : degrees + Maximum phase in degrees of the N-circles + + Return values + ------------- + contours : complex array + Array of complex numbers corresponding to the contours. + """ + # Convert magnitudes and phase range into a grid suitable for + # building contours + phases = sp.radians(sp.linspace(phase_min, phase_max, 500)) + Hmag, Hphase = sp.meshgrid(10.0**(mags/20.0), phases) + return closed_loop_contours(Hmag, Hphase) + +# N-circle +def n_circles(phases, mag_min=-40.0, mag_max=12.0): + """Constant-phase contours of the function H = G/(1+G). + + Usage + ===== + contour = n_circles(angles) + + Parameters + ---------- + phases : array-like + Array of phases in degrees of the N-circles + mag_min : dB + Minimum magnitude in dB of the N-circles + mag_max : dB + Maximum magnitude in dB of the N-circles + + Return values + ------------- + contours : complex array + Array of complex numbers corresponding to the contours. + """ + # Convert phases and magnitude range into a grid suitable for + # building contours + mags = sp.linspace(10**(mag_min/20.0), 10**(mag_max/20.0), 2000) + Hphase, Hmag = sp.meshgrid(sp.radians(phases), mags) + return closed_loop_contours(Hmag, Hphase) Modified: branches/control-0.4b/src/statesp.py =================================================================== --- branches/control-0.4b/src/statesp.py 2011-04-01 04:46:30 UTC (rev 143) +++ branches/control-0.4b/src/statesp.py 2011-04-02 15:52:38 UTC (rev 144) @@ -78,7 +78,6 @@ from numpy.linalg import inv, det, solve from numpy.linalg.linalg import LinAlgError from scipy.signal import lti -from slycot import td04ad from lti import Lti import xferfcn @@ -457,6 +456,7 @@ # Already a state space system; just return it return sys elif isinstance(sys, xferfcn.TransferFunction): + from slycot import td04ad if len(kw): raise TypeError("If sys is a TransferFunction, _convertToStateSpace \ cannot take keywords.") Modified: branches/control-0.4b/src/xferfcn.py =================================================================== --- branches/control-0.4b/src/xferfcn.py 2011-04-01 04:46:30 UTC (rev 143) +++ branches/control-0.4b/src/xferfcn.py 2011-04-02 15:52:38 UTC (rev 144) @@ -79,7 +79,6 @@ polyadd, polymul, polyval, roots, sort, sqrt, zeros from scipy.signal import lti from copy import deepcopy -from slycot import tb04ad from lti import Lti import statesp @@ -707,6 +706,7 @@ return sys elif isinstance(sys, statesp.StateSpace): + from slycot import tb04ad if len(kw): raise TypeError("If sys is a StateSpace, _convertToTransferFunction \ cannot take keywords.") Added: branches/control-0.4b/tests/nichols_test.py =================================================================== --- branches/control-0.4b/tests/nichols_test.py (rev 0) +++ branches/control-0.4b/tests/nichols_test.py 2011-04-02 15:52:38 UTC (rev 144) @@ -0,0 +1,37 @@ +#!/usr/bin/env python +# +# nichols_test.py - test Nichols plot +# RMM, 31 Mar 2011 + +import unittest +import numpy as np +from control.matlab import * + +class TestStateSpace(unittest.TestCase): + """Tests for the Nichols plots.""" + + def setUp(self): + """Set up a system to test operations on.""" + + A = [[-3., 4., 2.], [-1., -3., 0.], [2., 5., 3.]] + B = [[1.], [-3.], [-2.]] + C = [[4., 2., -3.]] + D = [[0.]] + + self.sys = StateSpace(A, B, C, D) + + def testNichols(self): + """Generate a Nichols plot.""" + nichols(self.sys) + + def testNgrid(self): + """Generate a Nichols plot.""" + nichols(self.sys) + ngrid() + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestStateSpace) + + +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-04-01 04:46:36
|
Revision: 143 http://python-control.svn.sourceforge.net/python-control/?rev=143&view=rev Author: murrayrm Date: 2011-04-01 04:46:30 +0000 (Fri, 01 Apr 2011) Log Message: ----------- copied Sawyer's unwrap change to 0.4b Modified Paths: -------------- branches/control-0.4b/src/ctrlutil.py Property Changed: ---------------- branches/control-0.4b/src/ctrlutil.py Modified: branches/control-0.4b/src/ctrlutil.py =================================================================== --- branches/control-0.4b/src/ctrlutil.py 2011-03-31 16:43:19 UTC (rev 142) +++ branches/control-0.4b/src/ctrlutil.py 2011-04-01 04:46:30 UTC (rev 143) @@ -73,9 +73,9 @@ for k in range(len(angle)): # See if we need to account for angle wrapping if (angle[k] - last > period/2): - wrap = -period + wrap -= period elif (last - angle[k] > period/2): - wrap = period + wrap += period # Update the last value we have sene last = angle[k] Property changes on: branches/control-0.4b/src/ctrlutil.py ___________________________________________________________________ Added: svn:mergeinfo + /trunk/src/ctrlutil.py:140-142 This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bli...@us...> - 2011-03-31 16:43:25
|
Revision: 142 http://python-control.svn.sourceforge.net/python-control/?rev=142&view=rev Author: blinkminster Date: 2011-03-31 16:43:19 +0000 (Thu, 31 Mar 2011) Log Message: ----------- updated control.unwrap to unwrap phase changes greater than 2pi Modified Paths: -------------- trunk/src/ctrlutil.py Modified: trunk/src/ctrlutil.py =================================================================== --- trunk/src/ctrlutil.py 2011-03-30 19:41:58 UTC (rev 141) +++ trunk/src/ctrlutil.py 2011-03-31 16:43:19 UTC (rev 142) @@ -71,9 +71,9 @@ for k in range(len(angle)): # See if we need to account for angle wrapping if (angle[k] - last > period/2): - wrap = -period + wrap -= period elif (last - angle[k] > period/2): - wrap = period + wrap += period # Update the last value we have sene last = angle[k] This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-03-30 19:42:06
|
Revision: 141 http://python-control.svn.sourceforge.net/python-control/?rev=141&view=rev Author: murrayrm Date: 2011-03-30 19:41:58 +0000 (Wed, 30 Mar 2011) Log Message: ----------- moved unit tests into a separate directory; use 'python tests/test_all.py' to run all tests Modified Paths: -------------- branches/control-0.4b/ChangeLog branches/control-0.4b/Pending branches/control-0.4b/setup.py Added Paths: ----------- branches/control-0.4b/tests/ branches/control-0.4b/tests/bdalg_test.py branches/control-0.4b/tests/convert_test.py branches/control-0.4b/tests/freqresp.py branches/control-0.4b/tests/matlab_test.py branches/control-0.4b/tests/modelsimp_test.py branches/control-0.4b/tests/slycot_convert_test.py branches/control-0.4b/tests/statefbk_test.py branches/control-0.4b/tests/statesp_test.py branches/control-0.4b/tests/test_all.py branches/control-0.4b/tests/xferfcn_test.py Removed Paths: ------------- branches/control-0.4b/src/TestBDAlg.py branches/control-0.4b/src/TestConvert.py branches/control-0.4b/src/TestFreqRsp.py branches/control-0.4b/src/TestMatlab.py branches/control-0.4b/src/TestModelsimp.py branches/control-0.4b/src/TestSlycot.py branches/control-0.4b/src/TestStateSp.py branches/control-0.4b/src/TestStatefbk.py branches/control-0.4b/src/TestXferFcn.py Modified: branches/control-0.4b/ChangeLog =================================================================== --- branches/control-0.4b/ChangeLog 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/ChangeLog 2011-03-30 19:41:58 UTC (rev 141) @@ -1,3 +1,15 @@ +2011-03-30 Richard Murray <murray@malabar.local> + + * tests/: added top level subdirectory, to be used for unit tests. + The idea in putting the code here is that you can do 'setup.py test' + during installation to make sure everything is working correctly. + The test code would normally *not* be callable from the installed + module. + + * tests/*_test.py: moved from src/Test*.py + + * setup.py: updated version number. + 2011-02-13 Richard Murray <murray@sumatra.local> * src/*.py: added svn:keywords Id properly Modified: branches/control-0.4b/Pending =================================================================== --- branches/control-0.4b/Pending 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/Pending 2011-03-30 19:41:58 UTC (rev 141) @@ -9,7 +9,10 @@ to be implemented. OPEN BUGS - * step() doesn't handle systems with a pole at the origin (use lsim2) + * matlab.step() doesn't handle systems with a pole at the origin (use lsim2) + * tests/convert_test.py not working yet + * tests/test_all.py should report on failed tests + * tests/freqresp.py needs to be converted to unit test Transfer code from Roberto Bucher's yottalab to python-control acker - pole placement using Ackermann method Modified: branches/control-0.4b/setup.py =================================================================== --- branches/control-0.4b/setup.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/setup.py 2011-03-30 19:41:58 UTC (rev 141) @@ -3,12 +3,12 @@ from distutils.core import setup setup(name='control', - version='0.4a', + version='0.4b', description='Python Control Systems Library', author='Richard Murray', author_email='mu...@cd...', url='http://python-control.sourceforge.net', requires='scipy', + packages=['control'], package_dir = {'control':'src'}, - packages=['control'], ) Deleted: branches/control-0.4b/src/TestBDAlg.py =================================================================== --- branches/control-0.4b/src/TestBDAlg.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestBDAlg.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,166 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -from xferfcn import TransferFunction -from statesp import StateSpace -from bdalg import feedback -import unittest - -class TestFeedback(unittest.TestCase): - """These are tests for the feedback function in bdalg.py. Currently, some - of the tests are not implemented, or are not working properly. TODO: these - need to be fixed.""" - - def setUp(self): - """This contains some random LTI systems and scalars for testing.""" - - # Two random SISO systems. - self.sys1 = TransferFunction([1, 2], [1, 2, 3]) - self.sys2 = StateSpace([[1., 4.], [3., 2.]], [[1.], [-4.]], - [[1., 0.]], [[0.]]) - # Two random scalars. - self.x1 = 2.5 - self.x2 = -3. - - def testScalarScalar(self): - """Scalar system with scalar feedback block.""" - - ans1 = feedback(self.x1, self.x2) - ans2 = feedback(self.x1, self.x2, 1.) - - self.assertAlmostEqual(ans1.num[0][0][0] / ans1.den[0][0][0], - -2.5 / 6.5) - self.assertAlmostEqual(ans2.num[0][0][0] / ans2.den[0][0][0], 2.5 / 8.5) - - def testScalarSS(self): - """Scalar system with state space feedback block.""" - - ans1 = feedback(self.x1, self.sys2) - ans2 = feedback(self.x1, self.sys2, 1.) - - np.testing.assert_array_almost_equal(ans1.A, [[-1.5, 4.], [13., 2.]]) - np.testing.assert_array_almost_equal(ans1.B, [[2.5], [-10.]]) - np.testing.assert_array_almost_equal(ans1.C, [[-2.5, 0.]]) - np.testing.assert_array_almost_equal(ans1.D, [[2.5]]) - np.testing.assert_array_almost_equal(ans2.A, [[3.5, 4.], [-7., 2.]]) - np.testing.assert_array_almost_equal(ans2.B, [[2.5], [-10.]]) - np.testing.assert_array_almost_equal(ans2.C, [[2.5, 0.]]) - np.testing.assert_array_almost_equal(ans2.D, [[2.5]]) - - def testScalarTF(self): - """Scalar system with transfer function feedback block.""" - - ans1 = feedback(self.x1, self.sys1) - ans2 = feedback(self.x1, self.sys1, 1.) - - np.testing.assert_array_almost_equal(ans1.num, [[[2.5, 5., 7.5]]]) - np.testing.assert_array_almost_equal(ans1.den, [[[1., 4.5, 8.]]]) - np.testing.assert_array_almost_equal(ans2.num, [[[2.5, 5., 7.5]]]) - np.testing.assert_array_almost_equal(ans2.den, [[[1., -0.5, -2.]]]) - - def testSSScalar(self): - """State space system with scalar feedback block.""" - - ans1 = feedback(self.sys2, self.x1) - ans2 = feedback(self.sys2, self.x1, 1.) - - np.testing.assert_array_almost_equal(ans1.A, [[-1.5, 4.], [13., 2.]]) - np.testing.assert_array_almost_equal(ans1.B, [[1.], [-4.]]) - np.testing.assert_array_almost_equal(ans1.C, [[1., 0.]]) - np.testing.assert_array_almost_equal(ans1.D, [[0.]]) - np.testing.assert_array_almost_equal(ans2.A, [[3.5, 4.], [-7., 2.]]) - np.testing.assert_array_almost_equal(ans2.B, [[1.], [-4.]]) - np.testing.assert_array_almost_equal(ans2.C, [[1., 0.]]) - np.testing.assert_array_almost_equal(ans2.D, [[0.]]) - - def testSSSS1(self): - """State space system with state space feedback block.""" - - ans1 = feedback(self.sys2, self.sys2) - ans2 = feedback(self.sys2, self.sys2, 1.) - - np.testing.assert_array_almost_equal(ans1.A, [[1., 4., -1., 0.], - [3., 2., 4., 0.], [1., 0., 1., 4.], [-4., 0., 3., 2]]) - np.testing.assert_array_almost_equal(ans1.B, [[1.], [-4.], [0.], [0.]]) - np.testing.assert_array_almost_equal(ans1.C, [[1., 0., 0., 0.]]) - np.testing.assert_array_almost_equal(ans1.D, [[0.]]) - np.testing.assert_array_almost_equal(ans2.A, [[1., 4., 1., 0.], - [3., 2., -4., 0.], [1., 0., 1., 4.], [-4., 0., 3., 2.]]) - np.testing.assert_array_almost_equal(ans2.B, [[1.], [-4.], [0.], [0.]]) - np.testing.assert_array_almost_equal(ans2.C, [[1., 0., 0., 0.]]) - np.testing.assert_array_almost_equal(ans2.D, [[0.]]) - - def testSSSS2(self): - """State space system with state space feedback block, including a - direct feedthrough term.""" - - sys3 = StateSpace([[-1., 4.], [2., -3]], [[2.], [3.]], [[-3., 1.]], - [[-2.]]) - sys4 = StateSpace([[-3., -2.], [1., 4.]], [[-2.], [-6.]], [[2., -3.]], - [[3.]]) - - ans1 = feedback(sys3, sys4) - ans2 = feedback(sys3, sys4, 1.) - - np.testing.assert_array_almost_equal(ans1.A, - [[-4.6, 5.2, 0.8, -1.2], [-3.4, -1.2, 1.2, -1.8], - [-1.2, 0.4, -1.4, -4.4], [-3.6, 1.2, 5.8, -3.2]]) - np.testing.assert_array_almost_equal(ans1.B, - [[-0.4], [-0.6], [-0.8], [-2.4]]) - np.testing.assert_array_almost_equal(ans1.C, [[0.6, -0.2, -0.8, 1.2]]) - np.testing.assert_array_almost_equal(ans1.D, [[0.4]]) - np.testing.assert_array_almost_equal(ans2.A, - [[-3.57142857142857, 4.85714285714286, 0.571428571428571, - -0.857142857142857], - [-1.85714285714286, -1.71428571428571, 0.857142857142857, - -1.28571428571429], - [0.857142857142857, -0.285714285714286, -1.85714285714286, - -3.71428571428571], - [2.57142857142857, -0.857142857142857, 4.42857142857143, - -1.14285714285714]]) - np.testing.assert_array_almost_equal(ans2.B, [[0.285714285714286], - [0.428571428571429], [0.571428571428571], [1.71428571428571]]) - np.testing.assert_array_almost_equal(ans2.C, [[-0.428571428571429, - 0.142857142857143, -0.571428571428571, 0.857142857142857]]) - np.testing.assert_array_almost_equal(ans2.D, [[-0.285714285714286]]) - - - def testSSTF(self): - """State space system with transfer function feedback block.""" - - # This functionality is not implemented yet. - pass - - def testTFScalar(self): - """Transfer function system with scalar feedback block.""" - - ans1 = feedback(self.sys1, self.x1) - ans2 = feedback(self.sys1, self.x1, 1.) - - np.testing.assert_array_almost_equal(ans1.num, [[[1., 2.]]]) - np.testing.assert_array_almost_equal(ans1.den, [[[1., 4.5, 8.]]]) - np.testing.assert_array_almost_equal(ans2.num, [[[1., 2.]]]) - np.testing.assert_array_almost_equal(ans2.den, [[[1., -0.5, -2.]]]) - - def testTFSS(self): - """Transfer function system with state space feedback block.""" - - # This functionality is not implemented yet. - pass - - def testTFTF(self): - """Transfer function system with transfer function feedback block.""" - - ans1 = feedback(self.sys1, self.sys1) - ans2 = feedback(self.sys1, self.sys1, 1.) - - np.testing.assert_array_almost_equal(ans1.num, [[[1., 4., 7., 6.]]]) - np.testing.assert_array_almost_equal(ans1.den, - [[[1., 4., 11., 16., 13.]]]) - np.testing.assert_array_almost_equal(ans2.num, [[[1., 4., 7., 6.]]]) - np.testing.assert_array_almost_equal(ans2.den, [[[1., 4., 9., 8., 5.]]]) -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestFeedback) - -if __name__ == "__main__": - unittest.main() Deleted: branches/control-0.4b/src/TestConvert.py =================================================================== --- branches/control-0.4b/src/TestConvert.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestConvert.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,91 +0,0 @@ -#!/usr/bin/env python - -"""TestConvert.py - -Test state space and transfer function conversion. - -Currently, this unit test script is not complete. It converts several random -state spaces back and forth between state space and transfer function -representations. Ideally, it should be able to assert that the conversion -outputs are correct. This is not yet implemented. - -Also, the conversion seems to enter an infinite loop once in a while. The cause -of this is unknown. - -""" - -import numpy as np -import matlab -import unittest - -class TestConvert(unittest.TestCase): - """Test state space and transfer function conversions.""" - - def setUp(self): - """Set up testing parameters.""" - - # Number of times to run each of the randomized tests. - self.numTests = 1 #almost guarantees failure - # Maximum number of states to test + 1 - self.maxStates = 20 - # Maximum number of inputs and outputs to test + 1 - self.maxIO = 20 - # Set to True to print systems to the output. - self.debug = True - - def printSys(self, sys, ind): - """Print system to the standard output.""" - - if self.debug: - print "sys%i:\n" % ind - print sys - - def testConvert(self): - """Test state space to transfer function conversion.""" - #Currently it only tests that a TF->SS->TF generates an unchanged TF - - #print __doc__ - - for states in range(1, self.maxStates): - for inputs in range(1, self.maxIO): - for outputs in range(1, self.maxIO): - #start with a random SS system and transform to TF - #then back to SS, check that the matrices are the same. - ssOriginal = matlab.rss(states, inputs, outputs) - self.printSys(ssOriginal, 1) - - tfOriginal = matlab.tf(ssOriginal) - self.printSys(tfOriginal, 2) - - ssTransformed = matlab.ss(tfOriginal) - self.printSys(ssTransformed, 3) - - tfTransformed = matlab.tf(ssTransformed) - self.printSys(tfTransformed, 4) - - for inputNum in range(inputs): - for outputNum in range(outputs): - np.testing.assert_array_almost_equal(\ - tfOriginal.num[outputNum][inputNum],\ - tfTransformed.num[outputNum][inputNum]) - - np.testing.assert_array_almost_equal(\ - tfOriginal.den[outputNum][inputNum],\ - tfTransformed.den[outputNum][inputNum]) - - #To test the ss systems is harder because they aren't the same - #realization. This could be done with checking that they have the - #same freq response with bode, but apparently it doesn't work - #the way it should right now: - ## Bode should work like this: - #[mag,phase,freq]=bode(sys) - #it doesn't seem to...... - #This should be added. - - - -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestConvert) - -if __name__ == "__main__": - unittest.main() Deleted: branches/control-0.4b/src/TestFreqRsp.py =================================================================== --- branches/control-0.4b/src/TestFreqRsp.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestFreqRsp.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,63 +0,0 @@ -#!/usr/bin/env python - -### MUST BE CONVERTED TO A UNIT TEST!!! - - -# Script to test frequency response and frequency response plots like bode, nyquist and gang of 4. -# Especially need to ensure that nothing SISO is broken and that MIMO at least handles exceptions and has some default to SISO in place. - - -import unittest -from statesp import StateSpace -from matlab import ss, tf, bode -import numpy as np -import matplotlib.pyplot as plt - -# SISO -plt.close('all') - -A = np.matrix('1,1;0,1') -B = np.matrix('0;1') -C = np.matrix('1,0') -D = 0 -sys = StateSpace(A,B,C,D) -#or try -#sys = ss(A,B,C,D) - -# test frequency response -omega = np.linspace(10e-2,10e2,1000) -frq=sys.freqresp(omega) - -# MIMO -B = np.matrix('1,0;0,1') -D = np.matrix('0,0') -sysMIMO = ss(A,B,C,D) -frqMIMO = sysMIMO.freqresp(omega) - -plt.figure(1) -bode(sys) - -systf = tf(sys) -tfMIMO = tf(sysMIMO) - -print systf.pole() -#print tfMIMO.pole() # - should throw not implemented exception -#print tfMIMO.zero() # - should throw not implemented exception - -plt.figure(2) -bode(systf) - -#bode(sysMIMO) # - should throw not implemented exception -#bode(tfMIMO) # - should throw not implemented exception - -#plt.figure(3) -#plt.semilogx(omega,20*np.log10(np.squeeze(frq[0]))) - -#plt.figure(4) -#bode(sysMIMO,omega) - - -def suite(): - pass - #Uncomment this once it is a real unittest - #return unittest.TestLoader().loadTestsFromTestCase(TestFreqRsp) Deleted: branches/control-0.4b/src/TestMatlab.py =================================================================== --- branches/control-0.4b/src/TestMatlab.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestMatlab.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,47 +0,0 @@ -#!/usr/bin/env python - -from matlab import * -import numpy as np -import unittest - -class TestMatlab(unittest.TestCase): - def testStep(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") - sys = ss(A,B,C,D) - t = np.linspace(0, 1, 10) - t, yout = step(sys, T=t) - youttrue = np.matrix("9. 17.6457 24.7072 30.4855 35.2234 39.1165 42.3227 44.9694 47.1599 48.9776") - np.testing.assert_array_almost_equal(yout, youttrue,decimal=4) - - def testImpulse(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") - sys = ss(A,B,C,D) - t = np.linspace(0, 1, 10) - t, yout = impulse(sys, T=t) - youttrue = np.matrix("86. 70.1808 57.3753 46.9975 38.5766 31.7344 26.1668 21.6292 17.9245 14.8945") - np.testing.assert_array_almost_equal(yout, youttrue,decimal=4) - -# def testInitial(self): -# A = np.matrix("1. -2.; 3. -4.") -# B = np.matrix("5.; 7.") -# C = np.matrix("6. 8.") -# D = np.matrix("9.") -# sys = ss(A,B,C,D) -# t = np.linspace(0, 1, 10) -# x0 = np.matrix(".5; 1.") -# t, yout = initial(sys, T=t, X0=x0) -# youttrue = np.matrix("11. 8.1494 5.9361 4.2258 2.9118 1.9092 1.1508 0.5833 0.1645 -0.1391") -# np.testing.assert_array_almost_equal(yout, youttrue,decimal=4) - - -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestMatlab) - -if __name__ == '__main__': - unittest.main() Deleted: branches/control-0.4b/src/TestModelsimp.py =================================================================== --- branches/control-0.4b/src/TestModelsimp.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestModelsimp.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,96 +0,0 @@ -#!/usr/bin/env python - -from modelsimp import * -from matlab import * -import numpy as np -import unittest - -class TestModelsimp(unittest.TestCase): - def testHSVD(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") - sys = ss(A,B,C,D) - hsv = hsvd(sys) - hsvtrue = np.matrix("24.42686 0.5731395") - np.testing.assert_array_almost_equal(hsv, hsvtrue) - - def testMarkov(self): - U = np.matrix("1.; 1.; 1.; 1.; 1.") - Y = U - M = 3 - H = markov(Y,U,M) - Htrue = np.matrix("1.; 0.; 0.") - np.testing.assert_array_almost_equal( H, Htrue ) - - def testModredMatchDC(self): - #balanced realization computed in matlab for the transfer function: - # num = [1 11 45 32], den = [1 15 60 200 60] - A = np.matrix('-1.958, -1.194, 1.824, -1.464; \ - -1.194, -0.8344, 2.563, -1.351; \ - -1.824, -2.563, -1.124, 2.704; \ - -1.464, -1.351, -2.704, -11.08') - B = np.matrix('-0.9057; -0.4068; -0.3263; -0.3474') - C = np.matrix('-0.9057, -0.4068, 0.3263, -0.3474') - D = np.matrix('0.') - sys = ss(A,B,C,D) - rsys = modred(sys,[2, 3],'matchdc') - Artrue = np.matrix('-4.431, -4.552; -4.552, -5.361') - Brtrue = np.matrix('-1.362; -1.031') - Crtrue = np.matrix('-1.362, -1.031') - Drtrue = np.matrix('-0.08384') - np.testing.assert_array_almost_equal(rsys.A, Artrue,decimal=3) - np.testing.assert_array_almost_equal(rsys.B, Brtrue,decimal=3) - np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=3) - np.testing.assert_array_almost_equal(rsys.D, Drtrue,decimal=2) - - def testModredTruncate(self): - #balanced realization computed in matlab for the transfer function: - # num = [1 11 45 32], den = [1 15 60 200 60] - A = np.matrix('-1.958, -1.194, 1.824, -1.464; \ - -1.194, -0.8344, 2.563, -1.351; \ - -1.824, -2.563, -1.124, 2.704; \ - -1.464, -1.351, -2.704, -11.08') - B = np.matrix('-0.9057; -0.4068; -0.3263; -0.3474') - C = np.matrix('-0.9057, -0.4068, 0.3263, -0.3474') - D = np.matrix('0.') - sys = ss(A,B,C,D) - rsys = modred(sys,[2, 3],'truncate') - Artrue = np.matrix('-1.958, -1.194; -1.194, -0.8344') - Brtrue = np.matrix('-0.9057; -0.4068') - Crtrue = np.matrix('-0.9057, -0.4068') - Drtrue = np.matrix('0.') - np.testing.assert_array_almost_equal(rsys.A, Artrue) - np.testing.assert_array_almost_equal(rsys.B, Brtrue) - np.testing.assert_array_almost_equal(rsys.C, Crtrue) - np.testing.assert_array_almost_equal(rsys.D, Drtrue) - - def testBalredTruncate(self): - #controlable canonical realization computed in matlab for the transfer function: - # num = [1 11 45 32], den = [1 15 60 200 60] - A = np.matrix('-15., -7.5, -6.25, -1.875; \ - 8., 0., 0., 0.; \ - 0., 4., 0., 0.; \ - 0., 0., 1., 0.') - B = np.matrix('2.; 0.; 0.; 0.') - C = np.matrix('0.5, 0.6875, 0.7031, 0.5') - D = np.matrix('0.') - sys = ss(A,B,C,D) - orders = 2 - rsys = balred(sys,orders,method='truncate') - Artrue = np.matrix('-1.958, -1.194; -1.194, -0.8344') - Brtrue = np.matrix('0.9057; 0.4068') - Crtrue = np.matrix('0.9057, 0.4068') - Drtrue = np.matrix('0.') - np.testing.assert_array_almost_equal(rsys.A, Artrue,decimal=2) - np.testing.assert_array_almost_equal(rsys.B, Brtrue,decimal=4) - np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=4) - np.testing.assert_array_almost_equal(rsys.D, Drtrue,decimal=4) - -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestModelsimp) - - -if __name__ == '__main__': - unittest.main() Deleted: branches/control-0.4b/src/TestSlycot.py =================================================================== --- branches/control-0.4b/src/TestSlycot.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestSlycot.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,124 +0,0 @@ -#!/usr/bin/env python - - -import numpy as np -from slycot import tb04ad, td04ad -import matlab -import unittest - - -class TestSlycot(unittest.TestCase): - """TestSlycot compares transfer function and state space conversions for - various numbers of inputs,outputs and states. - 1. Usually passes for SISO systems of any state dim, occasonally, there will be a dimension mismatch if the original randomly generated ss system is not minimal because td04ad returns a minimal system. - - 2. For small systems with many inputs, n<<m, the tests fail because td04ad returns a minimal ss system which has fewer states than the original system. It is typical for systems with many more inputs than states to have extraneous states. - - 3. For systems with larger dimensions, n~>5 and with 2 or more outputs the conversion to statespace (td04ad) intermittently results in an equivalent realization of higher order than the original tf order. We think this has to do with minimum realization tolerances in the Fortran. The algorithm doesn't recognize that two denominators are identical and so it creates a system with nearly duplicate eigenvalues and double the state dimension. This should not be a problem in the python-control usage because the common_den() method finds repeated roots within a tolerance that we specify. - - Matlab: Matlab seems to force its statespace system output to have order less than or equal to the order of denominators provided, avoiding the problem of very large state dimension we describe in 3. It does however, still have similar problems with pole/zero cancellation such as we encounter in 2, where a statespace system may have fewer states than the original order of transfer function. - """ - def setUp(self): - """Define some test parameters.""" - self.numTests = 5 - self.maxStates = 10 - self.maxI = 1 - self.maxO = 1 - - def testTF(self): - """ Directly tests the functions tb04ad and td04ad through direct comparison of transfer function coefficients. - Similar to TestConvert, but tests at a lower level. - """ - for states in range(1, self.maxStates): - for inputs in range(1, self.maxI+1): - for outputs in range(1, self.maxO+1): - for testNum in range(self.numTests): - - ssOriginal = matlab.rss(states, inputs, outputs) - - print '====== Original SS ==========' - print ssOriginal - print 'states=',states - print 'inputs=',inputs - print 'outputs=',outputs - - - tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\ - tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad(states,inputs,outputs,\ - ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0) - - ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\ - = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=0.0) - - tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\ - tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad(ssTransformed_nr,\ - inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,ssTransformed_D,tol1=0.0) - #print 'size(Trans_A)=',ssTransformed_A.shape - print '===== Transformed SS ==========' - print matlab.ss(ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D) - #print 'Trans_nr=',ssTransformed_nr - #print 'tfOrig_index=',tfOriginal_index - #print 'tfOrig_ucoeff=',tfOriginal_ucoeff - #print 'tfOrig_dcoeff=',tfOriginal_dcoeff - #print 'tfTrans_index=',tfTransformed_index - #print 'tfTrans_ucoeff=',tfTransformed_ucoeff - #print 'tfTrans_dcoeff=',tfTransformed_dcoeff - #Compare the TF directly, must match - #numerators - np.testing.assert_array_almost_equal(tfOriginal_ucoeff,tfTransformed_ucoeff,decimal=3) - #denominators - np.testing.assert_array_almost_equal(tfOriginal_dcoeff,tfTransformed_dcoeff,decimal=3) - - def testFreqResp(self): - """Compare the bode reponses of the SS systems and TF systems to the original SS - They generally are different realizations but have same freq resp. - Currently this test may only be applied to SISO systems. - """ - for states in range(1,self.maxStates): - for testNum in range(self.numTests): - for inputs in range(1,1): - for outputs in range(1,1): - ssOriginal = matlab.rss(states, inputs, outputs) - - tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\ - tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad(states,inputs,outputs,\ - ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0) - - ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\ - = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=0.0) - - tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\ - tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad(\ - ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\ - ssTransformed_D,tol1=0.0) - - numTransformed = np.array(tfTransformed_ucoeff) - denTransformed = np.array(tfTransformed_dcoeff) - numOriginal = np.array(tfOriginal_ucoeff) - denOriginal = np.array(tfOriginal_dcoeff) - - ssTransformed = matlab.ss(ssTransformed_A,ssTransformed_B,ssTransformed_C,ssTransformed_D) - for inputNum in range(inputs): - for outputNum in range(outputs): - [ssOriginalMag,ssOriginalPhase,freq] = matlab.bode(ssOriginal,Plot=False) - [tfOriginalMag,tfOriginalPhase,freq] = matlab.bode(matlab.tf(numOriginal[outputNum][inputNum],denOriginal[outputNum]),Plot=False) - [ssTransformedMag,ssTransformedPhase,freq] = matlab.bode(ssTransformed,freq,Plot=False) - [tfTransformedMag,tfTransformedPhase,freq] = matlab.bode(matlab.tf(numTransformed[outputNum][inputNum],denTransformed[outputNum]),freq,Plot=False) - #print 'numOrig=',numOriginal[outputNum][inputNum] - #print 'denOrig=',denOriginal[outputNum] - #print 'numTrans=',numTransformed[outputNum][inputNum] - #print 'denTrans=',denTransformed[outputNum] - np.testing.assert_array_almost_equal(ssOriginalMag,tfOriginalMag,decimal=3) - np.testing.assert_array_almost_equal(ssOriginalPhase,tfOriginalPhase,decimal=3) - np.testing.assert_array_almost_equal(ssOriginalMag,ssTransformedMag,decimal=3) - np.testing.assert_array_almost_equal(ssOriginalPhase,ssTransformedPhase,decimal=3) - np.testing.assert_array_almost_equal(tfOriginalMag,tfTransformedMag,decimal=3) - np.testing.assert_array_almost_equal(tfOriginalPhase,tfTransformedPhase,decimal=2) - -#These are here for once the above is made into a unittest. -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestSlycot) - -if __name__=='__main__': - unittest.main() - Deleted: branches/control-0.4b/src/TestStateSp.py =================================================================== --- branches/control-0.4b/src/TestStateSp.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestStateSp.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,209 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import unittest -import matlab -from statesp import StateSpace - -class TestStateSpace(unittest.TestCase): - """Tests for the StateSpace class.""" - - def setUp(self): - """Set up a MIMO system to test operations on.""" - - A = [[-3., 4., 2.], [-1., -3., 0.], [2., 5., 3.]] - B = [[1., 4.], [-3., -3.], [-2., 1.]] - C = [[4., 2., -3.], [1., 4., 3.]] - D = [[-2., 4.], [0., 1.]] - - a = [[4., 1.], [2., -3]] - b = [[5., 2.], [-3., -3.]] - c = [[2., -4], [0., 1.]] - d = [[3., 2.], [1., -1.]] - - self.sys1 = StateSpace(A, B, C, D) - self.sys2 = StateSpace(a, b, c, d) - - def testPole(self): - """Evaluate the poles of a MIMO system.""" - - p = self.sys1.pole() - - np.testing.assert_array_almost_equal(p, [3.34747678408874, - -3.17373839204437 + 1.47492908003839j, - -3.17373839204437 - 1.47492908003839j]) - - def testZero(self): - """Evaluate the zeros of a SISO system.""" - - sys = StateSpace(self.sys1.A, [[3.], [-2.], [4.]], [[-1., 3., 2.]], [[-4.]]) - z = sys.zero() - - np.testing.assert_array_almost_equal(z, [4.26864638637134, - -3.75932319318567 + 1.10087776649554j, - -3.75932319318567 - 1.10087776649554j]) - - def testAdd(self): - """Add two MIMO systems.""" - - A = [[-3., 4., 2., 0., 0.], [-1., -3., 0., 0., 0.], - [2., 5., 3., 0., 0.], [0., 0., 0., 4., 1.], [0., 0., 0., 2., -3.]] - B = [[1., 4.], [-3., -3.], [-2., 1.], [5., 2.], [-3., -3.]] - C = [[4., 2., -3., 2., -4.], [1., 4., 3., 0., 1.]] - D = [[1., 6.], [1., 0.]] - - sys = self.sys1 + self.sys2 - - np.testing.assert_array_almost_equal(sys.A, A) - np.testing.assert_array_almost_equal(sys.B, B) - np.testing.assert_array_almost_equal(sys.C, C) - np.testing.assert_array_almost_equal(sys.D, D) - - def testSub(self): - """Subtract two MIMO systems.""" - - A = [[-3., 4., 2., 0., 0.], [-1., -3., 0., 0., 0.], - [2., 5., 3., 0., 0.], [0., 0., 0., 4., 1.], [0., 0., 0., 2., -3.]] - B = [[1., 4.], [-3., -3.], [-2., 1.], [5., 2.], [-3., -3.]] - C = [[4., 2., -3., -2., 4.], [1., 4., 3., 0., -1.]] - D = [[-5., 2.], [-1., 2.]] - - sys = self.sys1 - self.sys2 - - np.testing.assert_array_almost_equal(sys.A, A) - np.testing.assert_array_almost_equal(sys.B, B) - np.testing.assert_array_almost_equal(sys.C, C) - np.testing.assert_array_almost_equal(sys.D, D) - - def testMul(self): - """Multiply two MIMO systems.""" - - A = [[4., 1., 0., 0., 0.], [2., -3., 0., 0., 0.], [2., 0., -3., 4., 2.], - [-6., 9., -1., -3., 0.], [-4., 9., 2., 5., 3.]] - B = [[5., 2.], [-3., -3.], [7., -2.], [-12., -3.], [-5., -5.]] - C = [[-4., 12., 4., 2., -3.], [0., 1., 1., 4., 3.]] - D = [[-2., -8.], [1., -1.]] - - sys = self.sys1 * self.sys2 - - np.testing.assert_array_almost_equal(sys.A, A) - np.testing.assert_array_almost_equal(sys.B, B) - np.testing.assert_array_almost_equal(sys.C, C) - np.testing.assert_array_almost_equal(sys.D, D) - - def testEvalFr(self): - """Evaluate the frequency response at one frequency.""" - - A = [[-2, 0.5], [0.5, -0.3]] - B = [[0.3, -1.3], [0.1, 0.]] - C = [[0., 0.1], [-0.3, -0.2]] - D = [[0., -0.8], [-0.3, 0.]] - sys = StateSpace(A, B, C, D) - - resp = [[4.37636761487965e-05 - 0.0152297592997812j, - -0.792603938730853 + 0.0261706783369803j], - [-0.331544857768052 + 0.0576105032822757j, - 0.128919037199125 - 0.143824945295405j]] - - np.testing.assert_almost_equal(sys.evalfr(1.), resp) - - def testFreqResp(self): - """Evaluate the frequency response at multiple frequencies.""" - - A = [[-2, 0.5], [0.5, -0.3]] - B = [[0.3, -1.3], [0.1, 0.]] - C = [[0., 0.1], [-0.3, -0.2]] - D = [[0., -0.8], [-0.3, 0.]] - sys = StateSpace(A, B, C, D) - - truemag = [[[0.0852992637230322, 0.00103596611395218], - [0.935374692849736, 0.799380720864549]], - [[0.55656854563842, 0.301542699860857], - [0.609178071542849, 0.0382108097985257]]] - truephase = [[[-0.566195599644593, -1.68063565332582], - [3.0465958317514, 3.14141384339534]], - [[2.90457947657161, 3.10601268291914], - [-0.438157380501337, -1.40720969147217]]] - trueomega = [0.1, 10.] - - mag, phase, omega = sys.freqresp(trueomega) - - np.testing.assert_almost_equal(mag, truemag) - np.testing.assert_almost_equal(phase, truephase) - np.testing.assert_equal(omega, trueomega) - -class TestRss(unittest.TestCase): - """These are tests for the proper functionality of statesp.rss.""" - - def setUp(self): - # Number of times to run each of the randomized tests. - self.numTests = 100 - # Maxmimum number of states to test + 1 - self.maxStates = 10 - # Maximum number of inputs and outputs to test + 1 - self.maxIO = 5 - - def testShape(self): - """Test that rss outputs have the right state, input, and output - size.""" - - for states in range(1, self.maxStates): - for inputs in range(1, self.maxIO): - for outputs in range(1, self.maxIO): - sys = matlab.rss(states, inputs, outputs) - self.assertEqual(sys.states, states) - self.assertEqual(sys.inputs, inputs) - self.assertEqual(sys.outputs, outputs) - - def testPole(self): - """Test that the poles of rss outputs have a negative real part.""" - - for states in range(1, self.maxStates): - for inputs in range(1, self.maxIO): - for outputs in range(1, self.maxIO): - sys = matlab.rss(states, inputs, outputs) - p = sys.pole() - for z in p: - self.assertTrue(z.real < 0) - -class TestDrss(unittest.TestCase): - """These are tests for the proper functionality of statesp.drss.""" - - def setUp(self): - # Number of times to run each of the randomized tests. - self.numTests = 100 - # Maxmimum number of states to test + 1 - self.maxStates = 10 - # Maximum number of inputs and outputs to test + 1 - self.maxIO = 5 - - def testShape(self): - """Test that drss outputs have the right state, input, and output - size.""" - - for states in range(1, self.maxStates): - for inputs in range(1, self.maxIO): - for outputs in range(1, self.maxIO): - sys = matlab.drss(states, inputs, outputs) - self.assertEqual(sys.states, states) - self.assertEqual(sys.inputs, inputs) - self.assertEqual(sys.outputs, outputs) - - def testPole(self): - """Test that the poles of drss outputs have less than unit magnitude.""" - - for states in range(1, self.maxStates): - for inputs in range(1, self.maxIO): - for outputs in range(1, self.maxIO): - sys = matlab.drss(states, inputs, outputs) - p = sys.pole() - for z in p: - self.assertTrue(abs(z) < 1) - - -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestStateSpace) - - -if __name__ == "__main__": - unittest.main() Deleted: branches/control-0.4b/src/TestStatefbk.py =================================================================== --- branches/control-0.4b/src/TestStatefbk.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestStatefbk.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,88 +0,0 @@ -#!/usr/bin/env python - -from statefbk import ctrb, obsv, place, lqr, gram -from matlab import * -import numpy as np -import unittest - -class TestStatefbk(unittest.TestCase): - def testCtrbSISO(self): - A = np.matrix("1. 2.; 3. 4.") - B = np.matrix("5.; 7.") - Wctrue = np.matrix("5. 19.; 7. 43.") - Wc = ctrb(A,B) - np.testing.assert_array_almost_equal(Wc, Wctrue) - - def testCtrbMIMO(self): - A = np.matrix("1. 2.; 3. 4.") - B = np.matrix("5. 6.; 7. 8.") - Wctrue = np.matrix("5. 6. 19. 22.; 7. 8. 43. 50.") - Wc = ctrb(A,B) - np.testing.assert_array_almost_equal(Wc, Wctrue) - - def testObsvSISO(self): - A = np.matrix("1. 2.; 3. 4.") - C = np.matrix("5. 7.") - Wotrue = np.matrix("5. 7.; 26. 38.") - Wo = obsv(A,C) - np.testing.assert_array_almost_equal(Wo, Wotrue) - - def testObsvMIMO(self): - A = np.matrix("1. 2.; 3. 4.") - C = np.matrix("5. 6.; 7. 8.") - Wotrue = np.matrix("5. 6.; 7. 8.; 23. 34.; 31. 46.") - Wo = obsv(A,C) - np.testing.assert_array_almost_equal(Wo, Wotrue) - - def testCtrbObsvDuality(self): - A = np.matrix("1.2 -2.3; 3.4 -4.5") - B = np.matrix("5.8 6.9; 8. 9.1") - Wc = ctrb(A,B); - A = np.transpose(A) - C = np.transpose(B) - Wo = np.transpose(obsv(A,C)); - np.testing.assert_array_almost_equal(Wc,Wo) - - def testGramWc(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5. 6.; 7. 8.") - C = np.matrix("4. 5.; 6. 7.") - D = np.matrix("13. 14.; 15. 16.") - sys = ss(A, B, C, D) - Wctrue = np.matrix("18.5 24.5; 24.5 32.5") - Wc = gram(sys,'c') - np.testing.assert_array_almost_equal(Wc, Wctrue) - - def testGramWo(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5. 6.; 7. 8.") - C = np.matrix("4. 5.; 6. 7.") - D = np.matrix("13. 14.; 15. 16.") - sys = ss(A, B, C, D) - Wotrue = np.matrix("257.5 -94.5; -94.5 56.5") - Wo = gram(sys,'o') - np.testing.assert_array_almost_equal(Wo, Wotrue) - - def testGramWo2(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") - sys = ss(A,B,C,D) - Wotrue = np.matrix("198. -72.; -72. 44.") - Wo = gram(sys,'o') - np.testing.assert_array_almost_equal(Wo, Wotrue) - - def testGramsys(self): - num =[1.] - den = [1., 1., 1.] - sys = tf(num,den) - self.assertRaises(ValueError, gram, sys, 'o') - self.assertRaises(ValueError, gram, sys, 'c') - -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestStatefbk) - - -if __name__ == '__main__': - unittest.main() Deleted: branches/control-0.4b/src/TestXferFcn.py =================================================================== --- branches/control-0.4b/src/TestXferFcn.py 2011-03-30 16:27:27 UTC (rev 140) +++ branches/control-0.4b/src/TestXferFcn.py 2011-03-30 19:41:58 UTC (rev 141) @@ -1,437 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -from statesp import StateSpace -from xferfcn import TransferFunction, _convertToTransferFunction -import unittest - -class TestXferFcn(unittest.TestCase): - """These are tests for functionality and correct reporting of the transfer - function class. Throughout these tests, we will give different input - formats to the xTranferFunction constructor, to try to break it. These - tests have been verified in MATLAB.""" - - # Tests for raising exceptions. - - def testBadInputType(self): - """Give the constructor invalid input types.""" - - self.assertRaises(TypeError, TransferFunction, [[0., 1.], [2., 3.]], - [[5., 2.], [3., 0.]]) - - def testInconsistentDimension(self): - """Give the constructor a numerator and denominator of different - sizes.""" - - self.assertRaises(ValueError, TransferFunction, [[[1.]]], - [[[1.], [2., 3.]]]) - self.assertRaises(ValueError, TransferFunction, [[[1.]]], - [[[1.]], [[2., 3.]]]) - self.assertRaises(ValueError, TransferFunction, [[[1.]]], - [[[1.], [1., 2.]], [[5., 2.], [2., 3.]]]) - - def testInconsistentColumns(self): - """Give the constructor inputs that do not have the same number of - columns in each row.""" - - self.assertRaises(ValueError, TransferFunction, 1., - [[[1.]], [[2.], [3.]]]) - self.assertRaises(ValueError, TransferFunction, [[[1.]], [[2.], [3.]]], - 1.) - - def testZeroDenominator(self): - """Give the constructor a transfer function with a zero denominator.""" - - self.assertRaises(ValueError, TransferFunction, 1., 0.) - self.assertRaises(ValueError, TransferFunction, - [[[1.], [2., 3.]], [[-1., 4.], [3., 2.]]], - [[[1., 0.], [0.]], [[0., 0.], [2.]]]) - - def testAddInconsistentDimension(self): - """Add two transfer function matrices of different sizes.""" - - sys1 = TransferFunction([[[1., 2.]]], [[[4., 5.]]]) - sys2 = TransferFunction([[[4., 3.]], [[1., 2.]]], - [[[1., 6.]], [[2., 4.]]]) - self.assertRaises(ValueError, sys1.__add__, sys2) - self.assertRaises(ValueError, sys1.__sub__, sys2) - self.assertRaises(ValueError, sys1.__radd__, sys2) - self.assertRaises(ValueError, sys1.__rsub__, sys2) - - def testMulInconsistentDimension(self): - """Multiply two transfer function matrices of incompatible sizes.""" - - sys1 = TransferFunction([[[1., 2.], [4., 5.]], [[2., 5.], [4., 3.]]], - [[[6., 2.], [4., 1.]], [[6., 7.], [2., 4.]]]) - sys2 = TransferFunction([[[1.]], [[2.]], [[3.]]], - [[[4.]], [[5.]], [[6.]]]) - self.assertRaises(ValueError, sys1.__mul__, sys2) - self.assertRaises(ValueError, sys2.__mul__, sys1) - self.assertRaises(ValueError, sys1.__rmul__, sys2) - self.assertRaises(ValueError, sys2.__rmul__, sys1) - - # Tests for TransferFunction._truncatecoeff - - def testTruncateCoeff1(self): - """Remove extraneous zeros in polynomial representations.""" - - sys1 = TransferFunction([0., 0., 1., 2.], [[[0., 0., 0., 3., 2., 1.]]]) - - np.testing.assert_array_equal(sys1.num, [[[1., 2.]]]) - np.testing.assert_array_equal(sys1.den, [[[3., 2., 1.]]]) - - def testTruncateCoeff2(self): - """Remove extraneous zeros in polynomial representations.""" - - sys1 = TransferFunction([0., 0., 0.], 1.) - - np.testing.assert_array_equal(sys1.num, [[[0.]]]) - np.testing.assert_array_equal(sys1.den, [[[1.]]]) - - # Tests for TransferFunction.__neg__ - - def testNegScalar(self): - """Negate a direct feedthrough system.""" - - sys1 = TransferFunction(2., np.array([-3])) - sys2 = - sys1 - - np.testing.assert_array_equal(sys2.num, [[[-2.]]]) - np.testing.assert_array_equal(sys2.den, [[[-3.]]]) - - def testNegSISO(self): - """Negate a SISO system.""" - - sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1.]) - sys2 = - sys1 - - np.testing.assert_array_equal(sys2.num, [[[-1., -3., -5.]]]) - np.testing.assert_array_equal(sys2.den, [[[1., 6., 2., -1.]]]) - - def testNegMIMO(self): - """Negate a MIMO system.""" - - num1 = [[[1., 2.], [0., 3.], [2., -1.]], - [[1.], [4., 0.], [1., -4., 3.]]] - num3 = [[[-1., -2.], [0., -3.], [-2., 1.]], - [[-1.], [-4., 0.], [-1., 4., -3.]]] - den1 = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]], - [[3., 0., .0], [2., -1., -1.], [1.]]] - - sys1 = TransferFunction(num1, den1) - sys2 = - sys1 - sys3 = TransferFunction(num3, den1) - - for i in range(sys3.outputs): - for j in range(sys3.inputs): - np.testing.assert_array_equal(sys2.num[i][j], sys3.num[i][j]) - np.testing.assert_array_equal(sys2.den[i][j], sys3.den[i][j]) - - # Tests for TransferFunction.__add__ - - def testAddScalar(self): - """Add two direct feedthrough systems.""" - - sys1 = TransferFunction(1., [[[1.]]]) - sys2 = TransferFunction(np.array([2.]), [1.]) - sys3 = sys1 + sys2 - - np.testing.assert_array_equal(sys3.num, 3.) - np.testing.assert_array_equal(sys3.den, 1.) - - def testAddSISO(self): - """Add two SISO systems.""" - - sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1]) - sys2 = TransferFunction([[np.array([-1., 3.])]], [[[1., 0., -1.]]]) - sys3 = sys1 + sys2 - - # If sys3.num is [[[0., 20., 4., -8.]]], then this is wrong! - np.testing.assert_array_equal(sys3.num, [[[20., 4., -8]]]) - np.testing.assert_array_equal(sys3.den, [[[1., 6., 1., -7., -2., 1.]]]) - - def testAddMIMO(self): - """Add two MIMO systems.""" - - num1 = [[[1., 2.], [0., 3.], [2., -1.]], - [[1.], [4., 0.], [1., -4., 3.]]] - den1 = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]], - [[3., 0., .0], [2., -1., -1.], [1.]]] - num2 = [[[0., 0., -1], [2.], [-1., -1.]], - [[1., 2.], [-1., -2.], [4.]]] - den2 = [[[-1.], [1., 2., 3.], [-1., -1.]], - [[-4., -3., 2.], [0., 1.], [1., 0.]]] - num3 = [[[3., -3., -6], [5., 6., 9.], [-4., -2., 2]], - [[3., 2., -3., 2], [-2., -3., 7., 2.], [1., -4., 3., 4]]] - den3 = [[[3., -2., -4.], [1., 2., 3., 0., 0.], [-2., -1., 1.]], - [[-12., -9., 6., 0., 0.], [2., -1., -1.], [1., 0.]]] - - sys1 = TransferFunction(num1, den1) - sys2 = TransferFunction(num2, den2) - sys3 = sys1 + sys2 - - for i in range(sys3.outputs): - for j in range(sys3.inputs): - np.testing.assert_array_equal(sys3.num[i][j], num3[i][j]) - np.testing.assert_array_equal(sys3.den[i][j], den3[i][j]) - - # Tests for TransferFunction.__sub__ - - def testSubScalar(self): - """Add two direct feedthrough systems.""" - - sys1 = TransferFunction(1., [[[1.]]]) - sys2 = TransferFunction(np.array([2.]), [1.]) - sys3 = sys1 - sys2 - - np.testing.assert_array_equal(sys3.num, -1.) - np.testing.assert_array_equal(sys3.den, 1.) - - def testSubSISO(self): - """Add two SISO systems.""" - - sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1]) - sys2 = TransferFunction([[np.array([-1., 3.])]], [[[1., 0., -1.]]]) - sys3 = sys1 - sys2 - sys4 = sys2 - sys1 - - np.testing.assert_array_equal(sys3.num, [[[2., 6., -12., -10., -2.]]]) - np.testing.assert_array_equal(sys3.den, [[[1., 6., 1., -7., -2., 1.]]]) - np.testing.assert_array_equal(sys4.num, [[[-2., -6., 12., 10., 2.]]]) - np.testing.assert_array_equal(sys4.den, [[[1., 6., 1., -7., -2., 1.]]]) - - def testSubMIMO(self): - """Add two MIMO systems.""" - - num1 = [[[1., 2.], [0., 3.], [2., -1.]], - [[1.], [4., 0.], [1., -4., 3.]]] - den1 = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]], - [[3., 0., .0], [2., -1., -1.], [1.]]] - num2 = [[[0., 0., -1], [2.], [-1., -1.]], - [[1., 2.], [-1., -2.], [4.]]] - den2 = [[[-1.], [1., 2., 3.], [-1., -1.]], - [[-4., -3., 2.], [0., 1.], [1., 0.]]] - num3 = [[[-3., 1., 2.], [1., 6., 9.], [0.]], - [[-3., -10., -3., 2], [2., 3., 1., -2], [1., -4., 3., -4]]] - den3 = [[[3., -2., -4], [1., 2., 3., 0., 0.], [1]], - [[-12., -9., 6., 0., 0.], [2., -1., -1], [1., 0.]]] - - sys1 = TransferFunction(num1, den1) - sys2 = TransferFunction(num2, den2) - sys3 = sys1 - sys2 - - for i in range(sys3.outputs): - for j in range(sys3.inputs): - np.testing.assert_array_equal(sys3.num[i][j], num3[i][j]) - np.testing.assert_array_equal(sys3.den[i][j], den3[i][j]) - - # Tests for TransferFunction.__mul__ - - def testMulScalar(self): - """Multiply two direct feedthrough systems.""" - - sys1 = TransferFunction(2., [1.]) - sys2 = TransferFunction(1., 4.) - sys3 = sys1 * sys2 - sys4 = sys1 * sys2 - - np.testing.assert_array_equal(sys3.num, [[[2.]]]) - np.testing.assert_array_equal(sys3.den, [[[4.]]]) - np.testing.assert_array_equal(sys3.num, sys4.num) - np.testing.assert_array_equal(sys3.den, sys4.den) - - def testMulSISO(self): - """Multiply two SISO systems.""" - - sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1]) - sys2 = TransferFunction([[[-1., 3.]]], [[[1., 0., -1.]]]) - sys3 = sys1 * sys2 - sys4 = sys2 * sys1 - - np.testing.assert_array_equal(sys3.num, [[[-1., 0., 4., 15.]]]) - np.testing.assert_array_equal(sys3.den, [[[1., 6., 1., -7., -2., 1.]]]) - np.testing.assert_array_equal(sys3.num, sys4.num) - np.testing.assert_array_equal(sys3.den, sys4.den) - - def testMulMIMO(self): - """Multiply two MIMO systems.""" - - num1 = [[[1., 2.], [0., 3.], [2., -1.]], - [[1.], [4., 0.], [1., -4., 3.]]] - den1 = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]], - [[3., 0., .0], [2., -1., -1.], [1.]]] - num2 = [[[0., 1., 2.]], - [[1., -5.]], - [[-2., 1., 4.]]] - den2 = [[[1., 0., 0., 0.]], - [[-2., 1., 3.]], - [[4., -1., -1., 0.]]] - num3 = [[[-24., 52., -14., 245., -490., -115., 467., -95., -56., 12., - 0., 0., 0.]], - [[24., -132., 138., 345., -768., -106., 510., 41., -79., -69., - -23., 17., 6., 0.]]] - den3 = [[[48., -92., -84., 183., 44., -97., -2., 12., 0., 0., 0., 0., - 0., 0.]], - [[-48., 60., 84., -81., -45., 21., 9., 0., 0., 0., 0., 0., 0.]]] - - sys1 = TransferFunction(num1, den1) - sys2 = TransferFunction(num2, den2) - sys3 = sys1 * sys2 - - for i in range(sys3.outputs): - for j in range(sys3.inputs): - np.testing.assert_array_equal(sys3.num[i][j], num3[i][j]) - np.testing.assert_array_equal(sys3.den[i][j], den3[i][j]) - - # Tests for TransferFunction.__div__ - - def testDivScalar(self): - """Divide two direct feedthrough systems.""" - - sys1 = TransferFunction(np.array([3.]), -4.) - sys2 = TransferFunction(5., 2.) - sys3 = sys1 / sys2 - - np.testing.assert_array_equal(sys3.num, [[[6.]]]) - np.testing.assert_array_equal(sys3.den, [[[-20.]]]) - - def testDivSISO(self): - """Divide two SISO systems.""" - - sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1]) - sys2 = TransferFunction([[[-1., 3.]]], [[[1., 0., -1.]]]) - sys3 = sys1 / sys2 - sys4 = sys2 / sys1 - - np.testing.assert_array_equal(sys3.num, [[[1., 3., 4., -3., -5.]]]) - np.testing.assert_array_equal(sys3.den, [[[-1., -3., 16., 7., -3.]]]) - np.testing.assert_array_equal(sys4.num, sys3.den) - np.testing.assert_array_equal(sys4.den, sys3.num) - - # Tests for TransferFunction.evalfr. - - def testEvalFrSISO(self): - """Evaluate the frequency response of a SISO system at one frequency.""" - - sys = TransferFunction([1., 3., 5], [1., 6., 2., -1]) - - np.testing.assert_array_almost_equal(sys.evalfr(1.), - np.array([[-0.5 - 0.5j]])) - np.testing.assert_array_almost_equal(sys.evalfr(32.), - np.array([[0.00281959302585077 - 0.030628473607392j]])) - - def testEvalFrMIMO(self): - """Evaluate the frequency response of a MIMO system at one frequency.""" - - num = [[[1., 2.], [0., 3.], [2., -1.]], - [[1.], [4., 0.], [1., -4., 3.]]] - den = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]], - [[3., 0., .0], [2., -1., -1.], [1.]]] - sys = TransferFunction(num, den) - resp = [[0.147058823529412 + 0.0882352941176471j, -0.75, 1.], - [-0.083333333333333, -0.188235294117647 - 0.847058823529412j, - -1. - 8.j]] - - np.testing.assert_array_almost_equal(sys.evalfr(2.), resp) - - # Tests for TransferFunction.freqresp. - - def testFreqRespSISO(self): - """Evaluate the magnitude and phase of a SISO system at multiple - frequencies.""" - - sys = TransferFunction([1., 3., 5], [1., 6., 2., -1]) - - truemag = [[[4.63507337473906, 0.707106781186548, 0.0866592803995351]]] - truephase = [[[-2.89596891081488, -2.35619449019234, - -1.32655885133871]]] - trueomega = [0.1, 1., 10.] - - mag, phase, omega = sys.freqresp(trueomega) - - np.testing.assert_array_almost_equal(mag, truemag) - np.testing.assert_array_almost_equal(phase, truephase) - np.testing.assert_array_almost_equal(omega, trueomega) - - def testFreqRespMIMO(self): - """Evaluate the magnitude and phase of a MIMO system at multiple - frequencies.""" - - num = [[[1., 2.], [0., 3.], [2., -1.]], - [[1.], [4., 0.], [1., -4., 3.]]] - den = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]], - [[3., 0., .0], [2., -1., -1.], [1.]]] - sys = TransferFunction(num, den) - - trueomega = [0.1, 1., 10.] - truemag = [[[0.496287094505259, 0.307147558416976, 0.0334738176210382], - [300., 3., 0.03], [1., 1., 1.]], - [[33.3333333333333, 0.333333333333333, 0.00333333333333333], - [0.390285696125482, 1.26491106406735, 0.198759144198533], - [3.01663720059274, 4.47213595499958, 104.92378186093]]] - truephase = [[[3.7128711165168e-4, 0.185347949995695, 1.30770596539255], - [-np.pi, -np.pi, -np.pi], [0., 0., 0.]], - [[-np.pi, -np.pi, -np.pi], - [-1.66852323415362, -1.89254688119154, -1.62050658356412], - [-0.132989648369409, -1.1071487177940, -2.7504672066207]]] - - mag, phase, omega = sys.freqresp(trueomega) - - np.testing.assert_array_almost_equal(mag, truemag) - np.testing.assert_array_almost_equal(phase, truephase) - np.testing.assert_array_equal(omega, trueomega) - - # Tests for TransferFunction.pole and TransferFunction.zero. - - def testPoleMIMO(self): - """Test for correct MIMO poles.""" - - sys = TransferFunction([[[1.], [1.]], [[1.], [1.]]], - [[[1., 2.], [1., 3.]], [[1., 4., 4.], [1., 9., 14.]]]) - p = sys.pole() - - np.testing.assert_array_almost_equal(p, [-7., -3., -2., -2.]) - - # Tests for TransferFunction.feedback. - - def testFeedbackSISO(self): - """Test for correct SISO transfer function feedback.""" - - sys1 = TransferFunction([-1., 4.], [1., 3., 5.]) - sys2 = TransferFunction([2., 3., 0.], [1., -3., 4., 0]) - - sys3 = sys1.feedback(sys2) - sys4 = sys1.feedback(sys2, 1) - - np.testing.assert_array_equal(sys3.num, [[[-1., 7., -16., 16., 0.]]]) - np.testing.assert_array_equal(sys3.den, [[[1., 0., -2., 2., 32., 0.]]]) - np.testing.assert_array_equal(sys4.num, [[[-1., 7., -16., 16., 0.]]]) - np.testing.assert_array_equal(sys4.den, [[[1., 0., 2., -8., 8., 0.]]]) - - def testConvertToTransferFunctio... [truncated message content] |
From: <mur...@us...> - 2011-03-30 16:27:33
|
Revision: 140 http://python-control.svn.sourceforge.net/python-control/?rev=140&view=rev Author: murrayrm Date: 2011-03-30 16:27:27 +0000 (Wed, 30 Mar 2011) Log Message: ----------- fixed bug in specifying axis limits; numbers were in wrong form Modified Paths: -------------- trunk/examples/pvtol-nested.py Modified: trunk/examples/pvtol-nested.py =================================================================== --- trunk/examples/pvtol-nested.py 2011-03-30 16:19:59 UTC (rev 139) +++ trunk/examples/pvtol-nested.py 2011-03-30 16:27:27 UTC (rev 140) @@ -88,17 +88,17 @@ # Add crossover line subplot(magh); hold(True); -loglog([10^-4, 10^3], [1, 1], 'k-') +loglog([1e-4, 1e3], [1, 1], 'k-') # Replot phase starting at -90 degrees bode(L, logspace(-4, 3)); (mag, phase, w) = freqresp(L, logspace(-4, 3)); phase = phase - 360; subplot(phaseh); -semilogx([10^-4, 10^3], [-180, -180], 'k-') +semilogx([10^-4, 1e3], [-180, -180], 'k-') hold(True); semilogx(w, phase, 'b-') -axis([10^-4, 10^3, -360, 0]); +axis([1e-4, 10^3, -360, 0]); xlabel('Frequency [deg]'); ylabel('Phase [deg]'); # set(gca, 'YTick', [-360, -270, -180, -90, 0]); # set(gca, 'XTick', [10^-4, 10^-2, 1, 100]); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mur...@us...> - 2011-03-30 16:20:05
|
Revision: 139 http://python-control.svn.sourceforge.net/python-control/?rev=139&view=rev Author: murrayrm Date: 2011-03-30 16:19:59 +0000 (Wed, 30 Mar 2011) Log Message: ----------- creating new revision of version 0.4b, in preparation for making 0.4 the default (trunk) Added Paths: ----------- branches/control-0.4b/ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: Richard M. <mu...@cd...> - 2011-03-24 09:05:00
|
If you are going to make minor changes, I think it is fine to just do them in trunk. The main guideline is to make sure that you don't commit changes to trunk that break the test cases. As we convert to v0.4a (which I hope to make the default in a few weeks, if I can find the time), we'll start having easier ways to run the test suite. If you are going to do something substantial (ala the Princeton group changes), then you should create a branch (eg, "svn cp trunk branches/murray-24Mar11") and work there. When everything is working, we can merge the changes in the branch with trunk. Right now the density of work is sparse enough that I don't think we need to do too much to keep branches separate, etc, but if someone has suggestions on simple things to do that might improve our productivity, just send comments to the list. -richard On 21 Mar 2011, at 18:11 , Sawyer Fuller wrote: > also, any procedures you can point me to regarding etiquette? i haven't ever developed in a team before. Is there a branch i should commit to? is there a post somewhere that addresses how to contribute? > > ------------------------ > Sawyer B. Fuller > California Institute of Technology > http://www.cds.caltech.edu/~minster/ > > > > > On Mon, Mar 21, 2011 at 6:07 PM, Sawyer Fuller <mi...@ca...> wrote: > cool > > it's blinkminster > > ------------------------ > Sawyer B. Fuller > California Institute of Technology > http://www.cds.caltech.edu/~minster/ > > > > > On Mon, Mar 21, 2011 at 2:53 AM, Richard Murray <mu...@cd...> wrote: > What's your sourceforge id? I can add it to the list of people who can commit. > > -richard > > On 20 Mar 2011, at 18:23 , Sawyer Fuller wrote: > > > hey, i've a quick fix to ctrlutil.unwrap to make it successfully unwrap phase changes of greater than 4pi: > > > > in this section: > > > > for k in range(len(angle)): > > # See if we need to account for angle wrapping > > if (angle[k] - last > period/2): > > wrap = -period > > elif (last - angle[k] > period/2): > > wrap = period > > > > i've replaced setting the wrap value with incrementing it: > > > > for k in range(len(angle)): > > # See if we need to account for angle wrapping > > if (angle[k] - last > period/2): > > wrap -= period > > elif (last - angle[k] > period/2): > > wrap += period > > > > also, what's the procedure for getting commit access to the repository? > > > > ------------------------ > > Sawyer B. Fuller > > California Institute of Technology > > http://www.cds.caltech.edu/~minster/ > > > > > > ------------------------------------------------------------------------------ > > Colocation vs. Managed Hosting > > A question and answer guide to determining the best fit > > for your organization - today and in the future. > > http://p.sf.net/sfu/internap-sfd2d_______________________________________________ > > python-control-discuss mailing list > > pyt...@li... > > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > > |
From: Sawyer F. <mi...@ca...> - 2011-03-22 01:12:00
|
also, any procedures you can point me to regarding etiquette? i haven't ever developed in a team before. Is there a branch i should commit to? is there a post somewhere that addresses how to contribute? ------------------------ Sawyer B. Fuller California Institute of Technology http://www.cds.caltech.edu/~minster/ <http://alumni.media.mit.edu/~minster/> On Mon, Mar 21, 2011 at 6:07 PM, Sawyer Fuller <mi...@ca...> wrote: > cool > > it's blinkminster > > ------------------------ > Sawyer B. Fuller > California Institute of Technology > http://www.cds.caltech.edu/~minster/<http://alumni.media.mit.edu/~minster/> > > > > > On Mon, Mar 21, 2011 at 2:53 AM, Richard Murray <mu...@cd...>wrote: > >> What's your sourceforge id? I can add it to the list of people who can >> commit. >> >> -richard >> >> On 20 Mar 2011, at 18:23 , Sawyer Fuller wrote: >> >> > hey, i've a quick fix to ctrlutil.unwrap to make it successfully unwrap >> phase changes of greater than 4pi: >> > >> > in this section: >> > >> > for k in range(len(angle)): >> > # See if we need to account for angle wrapping >> > if (angle[k] - last > period/2): >> > wrap = -period >> > elif (last - angle[k] > period/2): >> > wrap = period >> > >> > i've replaced setting the wrap value with incrementing it: >> > >> > for k in range(len(angle)): >> > # See if we need to account for angle wrapping >> > if (angle[k] - last > period/2): >> > wrap -= period >> > elif (last - angle[k] > period/2): >> > wrap += period >> > >> > also, what's the procedure for getting commit access to the repository? >> > >> > ------------------------ >> > Sawyer B. Fuller >> > California Institute of Technology >> > http://www.cds.caltech.edu/~minster/ >> > >> > >> > >> ------------------------------------------------------------------------------ >> > Colocation vs. Managed Hosting >> > A question and answer guide to determining the best fit >> > for your organization - today and in the future. >> > >> http://p.sf.net/sfu/internap-sfd2d_______________________________________________ >> > python-control-discuss mailing list >> > pyt...@li... >> > https://lists.sourceforge.net/lists/listinfo/python-control-discuss >> >> > |
From: Sawyer F. <mi...@ca...> - 2011-03-21 01:23:50
|
hey, i've a quick fix to ctrlutil.unwrap to make it successfully unwrap phase changes of greater than 4pi: in this section: for k in range(len(angle)): # See if we need to account for angle wrapping if (angle[k] - last > period/2): wrap = -period elif (last - angle[k] > period/2): wrap = period i've replaced setting the wrap value with incrementing it: for k in range(len(angle)): # See if we need to account for angle wrapping if (angle[k] - last > period/2): wrap -= period elif (last - angle[k] > period/2): wrap += period also, what's the procedure for getting commit access to the repository? ------------------------ Sawyer B. Fuller California Institute of Technology http://www.cds.caltech.edu/~minster/ <http://alumni.media.mit.edu/~minster/> |
From: Lauren P. <lpa...@Pr...> - 2011-03-11 23:36:22
|
Hi all, All of the additional Slycot routines needed for python-control branch 0.4a, (currently at revision 138) are now available through the Slycot repository on github. Enrico has tagged the python-control compatible version for us as v0.1.0 since the Slycot library is under frequent development and may change in the future. One item of note: The SLICOT library is down, with the following message at www.slicot.org: "Dear prospective SLICOT user, due to inconsistencies or misunderstanding of the SLICOT licensing policies we currently have to stop access to the sources and to re-organize the licensing. Please visit us later again. Thank you. In urgent cases please contact gf...@sy...." I'm not sure if this is cause for concern but wanted to bring it to your attention. Cheers, Lauren |
From: <lpa...@us...> - 2011-03-11 19:36:50
|
Revision: 138 http://python-control.svn.sourceforge.net/python-control/?rev=138&view=rev Author: lpadilla Date: 2011-03-11 19:36:44 +0000 (Fri, 11 Mar 2011) Log Message: ----------- Changed signature of tb04ad to match the routine now available in Slycot version on GitHub tagged as v0.1.0. Modified Paths: -------------- branches/control-0.4a/src/TestSlycot.py branches/control-0.4a/src/xferfcn.py Modified: branches/control-0.4a/src/TestSlycot.py =================================================================== --- branches/control-0.4a/src/TestSlycot.py 2011-02-19 03:15:18 UTC (rev 137) +++ branches/control-0.4a/src/TestSlycot.py 2011-03-11 19:36:44 UTC (rev 138) @@ -44,16 +44,15 @@ tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\ - tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad('R',states,inputs,outputs,\ + tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad(states,inputs,outputs,\ ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0) ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\ = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=0.0) tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\ - tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad('R',\ - ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\ - ssTransformed_D,tol1=0.0) + tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad(ssTransformed_nr,\ + inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,ssTransformed_D,tol1=0.0) #print 'size(Trans_A)=',ssTransformed_A.shape print '===== Transformed SS ==========' print matlab.ss(ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D) @@ -82,14 +81,14 @@ ssOriginal = matlab.rss(states, inputs, outputs) tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\ - tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad('R',states,inputs,outputs,\ + tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad(states,inputs,outputs,\ ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0) ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\ = td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=0.0) tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\ - tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad('R',\ + tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad(\ ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\ ssTransformed_D,tol1=0.0) Modified: branches/control-0.4a/src/xferfcn.py =================================================================== --- branches/control-0.4a/src/xferfcn.py 2011-02-19 03:15:18 UTC (rev 137) +++ branches/control-0.4a/src/xferfcn.py 2011-03-11 19:36:44 UTC (rev 138) @@ -715,7 +715,7 @@ # buggy! print "Warning: state space to transfer function conversion by tb04ad \ is still buggy!" - tfout = tb04ad('R',sys.states, sys.inputs, sys.outputs, sys.A, sys.B, sys.C, + tfout = tb04ad(sys.states, sys.inputs, sys.outputs, sys.A, sys.B, sys.C, sys.D,tol1=0.0) # Preallocate outputs. This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: Lauren P. <lpa...@pr...> - 2011-03-06 14:19:43
|
Hi Richard and others, So far 3 out of the 5 five additional routines are in now in Slycot, available with latest download from Github. Still waiting on tb04ad and td04ad. Will keep you informed. Lauren On Mar 5, 2011, at 3:37 PM, Richard Murray wrote: All of the functions have internal documentation that you can access with import control.matlab help control.matlab help control.<function> There are also some fairly detailed examples on the project wiki: http://sourceforge.net/apps/mediawiki/python-control/index.php?title=Main_Page The new version (2.0a) will have more complete documentation. Just waiting for slycot to have the right wrappers so that the installation will work without lots of manual intervention (Lauren: any updates on this?). -richard On 5 Mar 2011, at 12:32 , JACOBSON CLAS wrote: > Installation works fine for the control-python library. Is there documentation for the library somewhere? > > Clas A. Jacobson > 138 Metcalf Road > Tolland, CT 06084 > > 860 875 7577 (home) > jac...@sb... > ------------------------------------------------------------------------------ > What You Don't Know About Data Connectivity CAN Hurt You > This paper provides an overview of data connectivity, details > its effect on application quality, and explores various alternative > solutions. http://p.sf.net/sfu/progress-d2d_______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss ------------------------------------------------------------------------------ What You Don't Know About Data Connectivity CAN Hurt You This paper provides an overview of data connectivity, details its effect on application quality, and explores various alternative solutions. http://p.sf.net/sfu/progress-d2d _______________________________________________ python-control-discuss mailing list pyt...@li... https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
From: Richard M. <mu...@cd...> - 2011-03-05 20:37:48
|
All of the functions have internal documentation that you can access with import control.matlab help control.matlab help control.<function> There are also some fairly detailed examples on the project wiki: http://sourceforge.net/apps/mediawiki/python-control/index.php?title=Main_Page The new version (2.0a) will have more complete documentation. Just waiting for slycot to have the right wrappers so that the installation will work without lots of manual intervention (Lauren: any updates on this?). -richard On 5 Mar 2011, at 12:32 , JACOBSON CLAS wrote: > Installation works fine for the control-python library. Is there documentation for the library somewhere? > > Clas A. Jacobson > 138 Metcalf Road > Tolland, CT 06084 > > 860 875 7577 (home) > jac...@sb... > ------------------------------------------------------------------------------ > What You Don't Know About Data Connectivity CAN Hurt You > This paper provides an overview of data connectivity, details > its effect on application quality, and explores various alternative > solutions. http://p.sf.net/sfu/progress-d2d_______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
From: JACOBSON C. <jac...@sb...> - 2011-03-05 20:32:08
|
Installation works fine for the control-python library. Is there documentation for the library somewhere? Clas A. Jacobson 138 Metcalf Road Tolland, CT 06084 860 875 7577 (home) jac...@sb... |