Menu
▾
▴
python-control-discuss
|
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-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-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 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.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X