Menu
▾
▴
python-control-discuss — Discussion list for Python Control Systems Library
You can subscribe to this list here.
| 2010 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(19) |
Jul
|
Aug
|
Sep
|
Oct
|
Nov
(11) |
Dec
(5) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2011 |
Jan
|
Feb
(113) |
Mar
(12) |
Apr
(27) |
May
(2) |
Jun
(16) |
Jul
(6) |
Aug
(6) |
Sep
|
Oct
(3) |
Nov
(9) |
Dec
(2) |
| 2012 |
Jan
(5) |
Feb
(11) |
Mar
|
Apr
(3) |
May
|
Jun
|
Jul
|
Aug
(3) |
Sep
(7) |
Oct
(18) |
Nov
(18) |
Dec
|
| 2013 |
Jan
(4) |
Feb
(1) |
Mar
(3) |
Apr
(1) |
May
|
Jun
(33) |
Jul
(2) |
Aug
(5) |
Sep
|
Oct
|
Nov
|
Dec
(1) |
| 2014 |
Jan
(1) |
Feb
|
Mar
(8) |
Apr
|
May
(3) |
Jun
(3) |
Jul
(9) |
Aug
(5) |
Sep
(6) |
Oct
|
Nov
|
Dec
|
| 2015 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
(4) |
| 2017 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
(5) |
Nov
|
Dec
|
| 2018 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(1) |
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2019 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
(6) |
Sep
|
Oct
|
Nov
(2) |
Dec
|
| 2020 |
Jan
(1) |
Feb
(1) |
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2021 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
(4) |
Oct
|
Nov
|
Dec
|
|
From: <kk...@us...> - 2011-02-08 22:13:15
|
Revision: 41
http://python-control.svn.sourceforge.net/python-control/?rev=41&view=rev
Author: kkchen
Date: 2011-02-08 22:13:09 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added rss unit tests.
Kevin K. Chen <kk...@pr...>
Added Paths:
-----------
branches/control-0.4a/src/TestStateSp.py
Added: branches/control-0.4a/src/TestStateSp.py
===================================================================
--- branches/control-0.4a/src/TestStateSp.py (rev 0)
+++ branches/control-0.4a/src/TestStateSp.py 2011-02-08 22:13:09 UTC (rev 41)
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+import numpy as np
+import statesp as ss
+import unittest
+
+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
+
+ def testShape(self):
+ """Test that rss outputs have the right state, input, and output
+ size."""
+
+ for states in range(1, 10):
+ for inputs in range(1, 5):
+ for outputs in range(1, 5):
+ sys = ss.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, 10):
+ for inputs in range(1, 5):
+ for outputs in range(1, 5):
+ sys = ss.rss(states, inputs, outputs)
+ p = sys.poles()
+ for z in p:
+ self.assertTrue(z.real < 0)
+
+if __name__ == "__main__":
+ unittest.main()
\ No newline at end of file
Property changes on: branches/control-0.4a/src/TestStateSp.py
___________________________________________________________________
Added: svn:executable
+ *
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:13:10
|
Revision: 40
http://python-control.svn.sourceforge.net/python-control/?rev=40&view=rev
Author: kkchen
Date: 2011-02-08 22:13:04 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added hsvd and markov functionality and test script TestModelsimp.py.
Steven Brunton <sbr...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestStatefbk.py
branches/control-0.4a/src/modelsimp.py
Added Paths:
-----------
branches/control-0.4a/src/TestModelsimp.py
Added: branches/control-0.4a/src/TestModelsimp.py
===================================================================
--- branches/control-0.4a/src/TestModelsimp.py (rev 0)
+++ branches/control-0.4a/src/TestModelsimp.py 2011-02-08 22:13:04 UTC (rev 40)
@@ -0,0 +1,28 @@
+#!/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 )
+
+if __name__ == '__main__':
+ unittest.main()
Modified: branches/control-0.4a/src/TestStatefbk.py
===================================================================
--- branches/control-0.4a/src/TestStatefbk.py 2011-02-08 22:12:57 UTC (rev 39)
+++ branches/control-0.4a/src/TestStatefbk.py 2011-02-08 22:13:04 UTC (rev 40)
@@ -48,8 +48,7 @@
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)
- sys = 1.
+ 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)
@@ -60,7 +59,6 @@
C = np.matrix("4. 5.; 6. 7.")
D = np.matrix("13. 14.; 15. 16.")
sys = ss(A, B, C, D)
- sys = 1.
Wotrue = np.matrix("257.5 -94.5; -94.5 56.5")
Wo = gram(sys,'o')
np.testing.assert_array_almost_equal(Wo, Wotrue)
Modified: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py 2011-02-08 22:12:57 UTC (rev 39)
+++ branches/control-0.4a/src/modelsimp.py 2011-02-08 22:13:04 UTC (rev 40)
@@ -53,6 +53,8 @@
=====
H = hsvd(sys)
+ The Hankel singular values are the singular values of the Hankel operator. In practice, we compute the square root of the eigenvalues of the matrix formed by taking the product of the observability and controllability gramians. There are other (more efficient) methods based on solving the Lyapunov equation in a particular way (more details soon).
+
Inputs
------
sys : a state space system
@@ -63,15 +65,16 @@
"""
- # Make sure that SLICOT is installed
- try:
- from slycot import sb01bd
- except ImportError:
- raise ControlSlycot("can't find slycot module 'sb01bd'")
-
- H = 1.
+ Wc = gram(sys,'c')
+ Wo = gram(sys,'o')
+
+ WoWc = np.dot(Wo, Wc)
+ w, v = LA.eig(WoWc)
+
+ hsv = np.sqrt(w)
+ hsv = np.matrix(hsv)
# Return the Hankel singular values
- return H
+ return hsv
def era(YY,m,n,nin,nout,r):
"""Calculate an ERA model of order r based on the impulse-response data YY
@@ -100,6 +103,7 @@
Usage
=====
H = markov(Y,U,M)
+ Currently only works for SISO
Inputs
------
@@ -113,3 +117,23 @@
"""
+ # Convert input parameters to matrices (if they aren't already)
+ Ymat = np.mat(Y)
+ Umat = np.mat(U)
+ n = np.size(U)
+
+ # Construct a matrix of control inputs to invert
+ UU = Umat
+ for i in range(1, M-1):
+ newCol = np.vstack((0, UU[0:n-1,i-2]))
+ UU = np.hstack((UU, newCol))
+ Ulast = np.vstack((0, UU[0:n-1,M-2]))
+ for i in range(n-1,0,-1):
+ Ulast[i] = np.sum(Ulast[0:i-1])
+ UU = np.hstack((UU, Ulast))
+
+ # Invert and solve for Markov parameters
+ H = UU.I
+ H = np.dot(H, Y)
+
+ return H
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:13:04
|
Revision: 39
http://python-control.svn.sourceforge.net/python-control/?rev=39&view=rev
Author: kkchen
Date: 2011-02-08 22:12:57 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Removed scipy.signal.lti from statesp.py; added rss.
The StateSpace class now uses seven data members: A, B, C, D, states, inputs,
and outputs. scipy.signal.lti was removed because it does not support MIMO
functionality.
Also, a working rss has been added to statesp.py. It still has yet to be
tested.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/statesp.py
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:12:53 UTC (rev 38)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:12:57 UTC (rev 39)
@@ -55,14 +55,32 @@
# of the functions that already existing in that package to be used
# directly.
#
-class StateSpace(signal.lti):
+class StateSpace:
"""The StateSpace class is used to represent linear input/output systems.
"""
# Initialization
- def __init__(self, *args, **keywords):
- # First initialize the parent object
- signal.lti.__init__(self, *args, **keywords)
+ def __init__(self, A, B, C, D):
+ self.A = A
+ self.B = B
+ self.C = C
+ self.D = D
+ self.states = A.shape[0]
+ self.inputs = B.shape[1]
+ self.outputs = C.shape[0]
+
+ # Check that the matrix sizes are consistent.
+ if self.states != A.shape[1]:
+ raise ValueError("A must be square.")
+ if self.states != B.shape[0]:
+ raise ValueError("B must have the same row size as A.")
+ if self.states != C.shape[1]:
+ raise ValueError("C must have the same column size as A.")
+ if self.inputs != D.shape[1]:
+ raise ValueError("D must have the same column size as B.")
+ if self.outputs != D.shape[0]:
+ raise ValueError("D must have the same row size as C.")
+
# Style to use for printing
def __str__(self):
str = "A = " + self.A.__str__() + "\n\n"
@@ -76,7 +94,7 @@
"""Compute the response of a system to a list of frequencies"""
# Generate and save a transfer function matrix
#! TODO: This is currently limited to SISO systems
- nout, nin = self.D.shape
+ #nout, nin = self.D.shape
# Compute the denominator from the A matrix
den = sp.poly1d(sp.poly(self.A))
@@ -99,7 +117,9 @@
return None
# Compute poles and zeros
- def poles(self): return sp.roots(sp.poly(self.A))
+ def poles(self):
+ return sp.roots(sp.poly(self.A))
+
def zeros(self):
den = sp.poly1d(sp.poly(self.A))
@@ -248,3 +268,83 @@
else:
raise TypeError("can't convert given type to StateSpace system")
+
+def rss(states=1, inputs=1, outputs=1):
+ """Create a stable random state space object."""
+
+ import numpy
+ from numpy.random import rand, randn
+
+ # Make some poles for A. Preallocate a complex array.
+ poles = numpy.zeros(states) + numpy.zeros(states) * 0.j
+ i = 0
+ while i < states - 1:
+ if rand() < 0.05 and i != 0:
+ # Small chance of copying poles, if we're not at the first element.
+ if poles[i-1].imag == 0:
+ # Copy previous real pole.
+ poles[i] = poles[i-1]
+ i += 1
+ else:
+ # Copy previous complex conjugate pair of poles.
+ poles[i:i+2] = poles[i-2:i]
+ i += 2
+ elif rand() < 0.6:
+ # Real pole.
+ poles[i] = -sp.exp(randn()) + 0.j
+ i += 1
+ else:
+ # Complex conjugate pair of poles.
+ poles[i] = complex(-sp.exp(randn()), sp.exp(randn()))
+ poles[i+1] = complex(poles[i].real, -poles[i].imag)
+ i += 2
+ # When we reach this point, we either have one or zero poles left to fill.
+ # Put a real pole if there is one space left.
+ if i == states - 1:
+ poles[i] = -sp.exp(randn()) + 0.j
+
+ # Now put the poles in A as real blocks on the diagonal.
+ A = numpy.zeros((states, states))
+ i = 0
+ while i < states:
+ if poles[i].imag == 0:
+ A[i, i] = poles[i].real
+ i += 1
+ else:
+ A[i, i] = A[i+1, i+1] = poles[i].real
+ A[i, i+1] = poles[i].imag
+ A[i+1, i] = -poles[i].imag
+ i += 2
+
+ # Finally, apply a transformation so that A is not block-diagonal.
+ while True:
+ T = randn(states, states)
+ try:
+ A = numpy.dot(numpy.linalg.solve(T, A), T) # A = T \ A * T
+ break
+ except numpy.linalg.linalg.LinAlgError:
+ # In the unlikely event that T is rank-deficient, iterate again.
+ pass
+
+ # Make the remaining matrices.
+ B = randn(states, inputs)
+ C = randn(outputs, states)
+ D = randn(outputs, inputs)
+
+ # Make masks to zero out some of the elements.
+ while True:
+ B_mask = rand(states, inputs) < 0.8
+ if sp.any(B_mask): # Retry if we get all zeros.
+ break
+ while True:
+ C_mask = rand(outputs, states) < 0.8
+ if sp.any(C_mask): # Retry if we get all zeros.
+ break
+ D_mask = rand(outputs, inputs) < 0.3
+
+ # Apply masks.
+ B = B * B_mask
+ C = C * C_mask
+ D = D * D_mask
+
+ return StateSpace(A, B, C, D)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:12:59
|
Revision: 38
http://python-control.svn.sourceforge.net/python-control/?rev=38&view=rev
Author: kkchen
Date: 2011-02-08 22:12:53 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added modelsimp.py with interfaces for hsvd, era, and markov functions. Also, modified TestStatefbk.py so that we import numpy as np, as in other files.
Steven Brunton <sbr...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestStatefbk.py
Added Paths:
-----------
branches/control-0.4a/src/modelsimp.py
Modified: branches/control-0.4a/src/TestStatefbk.py
===================================================================
--- branches/control-0.4a/src/TestStatefbk.py 2011-02-08 22:11:21 UTC (rev 37)
+++ branches/control-0.4a/src/TestStatefbk.py 2011-02-08 22:12:53 UTC (rev 38)
@@ -2,68 +2,68 @@
from statefbk import ctrb, obsv, place, lqr, gram
from matlab import *
-import numpy as N
+import numpy as np
import unittest
class TestStatefbk(unittest.TestCase):
def testCtrbSISO(self):
- A = N.matrix("1. 2.; 3. 4.")
- B = N.matrix("5.; 7.")
- Wctrue = N.matrix("5. 19.; 7. 43.")
+ A = np.matrix("1. 2.; 3. 4.")
+ B = np.matrix("5.; 7.")
+ Wctrue = np.matrix("5. 19.; 7. 43.")
Wc = ctrb(A,B)
- N.testing.assert_array_almost_equal(Wc, Wctrue)
+ np.testing.assert_array_almost_equal(Wc, Wctrue)
def testCtrbMIMO(self):
- A = N.matrix("1. 2.; 3. 4.")
- B = N.matrix("5. 6.; 7. 8.")
- Wctrue = N.matrix("5. 6. 19. 22.; 7. 8. 43. 50.")
+ 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)
- N.testing.assert_array_almost_equal(Wc, Wctrue)
+ np.testing.assert_array_almost_equal(Wc, Wctrue)
def testObsvSISO(self):
- A = N.matrix("1. 2.; 3. 4.")
- C = N.matrix("5. 7.")
- Wotrue = N.matrix("5. 7.; 26. 38.")
+ A = np.matrix("1. 2.; 3. 4.")
+ C = np.matrix("5. 7.")
+ Wotrue = np.matrix("5. 7.; 26. 38.")
Wo = obsv(A,C)
- N.testing.assert_array_almost_equal(Wo, Wotrue)
+ np.testing.assert_array_almost_equal(Wo, Wotrue)
def testObsvMIMO(self):
- A = N.matrix("1. 2.; 3. 4.")
- C = N.matrix("5. 6.; 7. 8.")
- Wotrue = N.matrix("5. 6.; 7. 8.; 23. 34.; 31. 46.")
+ 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)
- N.testing.assert_array_almost_equal(Wo, Wotrue)
+ np.testing.assert_array_almost_equal(Wo, Wotrue)
def testCtrbObsvDuality(self):
- A = N.matrix("1.2 -2.3; 3.4 -4.5")
- B = N.matrix("5.8 6.9; 8. 9.1")
+ 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 = N.transpose(A)
- C = N.transpose(B)
- Wo = N.transpose(obsv(A,C));
- N.testing.assert_array_almost_equal(Wc,Wo)
+ 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 = N.matrix("1. -2.; 3. -4.")
- B = N.matrix("5. 6.; 7. 8.")
- C = N.matrix("4. 5.; 6. 7.")
- D = N.matrix("13. 14.; 15. 16.")
+ 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)
sys = 1.
- Wctrue = N.matrix("18.5 24.5; 24.5 32.5")
+ Wctrue = np.matrix("18.5 24.5; 24.5 32.5")
Wc = gram(sys,'c')
- N.testing.assert_array_almost_equal(Wc, Wctrue)
+ np.testing.assert_array_almost_equal(Wc, Wctrue)
def testGramWo(self):
- A = N.matrix("1. -2.; 3. -4.")
- B = N.matrix("5. 6.; 7. 8.")
- C = N.matrix("4. 5.; 6. 7.")
- D = N.matrix("13. 14.; 15. 16.")
+ 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)
sys = 1.
- Wotrue = N.matrix("257.5 -94.5; -94.5 56.5")
+ Wotrue = np.matrix("257.5 -94.5; -94.5 56.5")
Wo = gram(sys,'o')
- N.testing.assert_array_almost_equal(Wo, Wotrue)
+ np.testing.assert_array_almost_equal(Wo, Wotrue)
if __name__ == '__main__':
unittest.main()
Added: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py (rev 0)
+++ branches/control-0.4a/src/modelsimp.py 2011-02-08 22:12:53 UTC (rev 38)
@@ -0,0 +1,115 @@
+# modelsimp.py - tools for model simplification
+#
+# Author: Steve Brunton, Kevin Chen, Lauren Padilla
+# Date: 30 Nov 2010
+#
+# This file contains routines for obtaining reduced order models
+#
+# 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$
+
+# External packages and modules
+import numpy as np
+import ctrlutil
+from control.exception import *
+
+# Hankel Singular Value Decomposition
+# The following returns the Hankel singular values, which are singular values of the matrix formed by multiplying the controllability and observability grammians
+def hsvd(sys):
+ """Calculate the Hankel singular values
+
+ Usage
+ =====
+ H = hsvd(sys)
+
+ Inputs
+ ------
+ sys : a state space system
+
+ Outputs
+ -------
+ H : a list of Hankel singular values
+
+ """
+
+ # Make sure that SLICOT is installed
+ try:
+ from slycot import sb01bd
+ except ImportError:
+ raise ControlSlycot("can't find slycot module 'sb01bd'")
+
+ H = 1.
+ # Return the Hankel singular values
+ return H
+
+def era(YY,m,n,nin,nout,r):
+ """Calculate an ERA model of order r based on the impulse-response data YY
+
+ Usage
+ =====
+ sys = era(YY,m,n,nin,nout,r)
+
+ Inputs
+ ------
+ YY : nout x nin dimensional impulse-response data
+ m : number of rows in Hankel matrix
+ n : number of columns in Hankel matrix
+ nin : number of input variables
+ nout : number of output variables
+ r : order of model
+
+ Outputs
+ -------
+ sys : a reduced order model sys=ss(Ar,Br,Cr,Dr)
+
+ """
+def markov(Y,U,M):
+ """Calculate the first M Markov parameters [D CB CAB ...] from input U, output Y
+
+ Usage
+ =====
+ H = markov(Y,U,M)
+
+ Inputs
+ ------
+ Y : output data
+ U : input data
+ M : number of Markov parameters to output
+
+ Outputs
+ -------
+ H : first M Markov parameters
+
+ """
+
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:11:28
|
Revision: 37
http://python-control.svn.sourceforge.net/python-control/?rev=37&view=rev
Author: kkchen
Date: 2011-02-08 22:11:21 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added shell interface for gram along with tests.
Steven Brunton <sbr...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestStatefbk.py
branches/control-0.4a/src/statefbk.py
Modified: branches/control-0.4a/src/TestStatefbk.py
===================================================================
--- branches/control-0.4a/src/TestStatefbk.py 2011-02-08 22:03:57 UTC (rev 36)
+++ branches/control-0.4a/src/TestStatefbk.py 2011-02-08 22:11:21 UTC (rev 37)
@@ -1,6 +1,7 @@
#!/usr/bin/env python
-from statefbk import ctrb, obsv, place, lqr
+from statefbk import ctrb, obsv, place, lqr, gram
+from matlab import *
import numpy as N
import unittest
@@ -42,5 +43,27 @@
Wo = N.transpose(obsv(A,C));
N.testing.assert_array_almost_equal(Wc,Wo)
+ def testGramWc(self):
+ A = N.matrix("1. -2.; 3. -4.")
+ B = N.matrix("5. 6.; 7. 8.")
+ C = N.matrix("4. 5.; 6. 7.")
+ D = N.matrix("13. 14.; 15. 16.")
+ # sys = ss(A, B, C, D)
+ sys = 1.
+ Wctrue = N.matrix("18.5 24.5; 24.5 32.5")
+ Wc = gram(sys,'c')
+ N.testing.assert_array_almost_equal(Wc, Wctrue)
+
+ def testGramWo(self):
+ A = N.matrix("1. -2.; 3. -4.")
+ B = N.matrix("5. 6.; 7. 8.")
+ C = N.matrix("4. 5.; 6. 7.")
+ D = N.matrix("13. 14.; 15. 16.")
+ sys = ss(A, B, C, D)
+ sys = 1.
+ Wotrue = N.matrix("257.5 -94.5; -94.5 56.5")
+ Wo = gram(sys,'o')
+ N.testing.assert_array_almost_equal(Wo, Wotrue)
+
if __name__ == '__main__':
unittest.main()
Modified: branches/control-0.4a/src/statefbk.py
===================================================================
--- branches/control-0.4a/src/statefbk.py 2011-02-08 22:03:57 UTC (rev 36)
+++ branches/control-0.4a/src/statefbk.py 2011-02-08 22:11:21 UTC (rev 37)
@@ -190,7 +190,7 @@
Usage
=====
- Wc = ctrb(A, B)
+ C = ctrb(A, B)
Inputs
------
@@ -198,7 +198,7 @@
Outputs
-------
- Wc: Controllability matrix
+ C: Controllability matrix
"""
# Convert input parameters to matrices (if they aren't already)
@@ -216,7 +216,7 @@
Usage
=====
- Wc = obsv(A, C)
+ O = obsv(A, C)
Inputs
------
@@ -224,7 +224,7 @@
Outputs
-------
- Wc: Observability matrix
+ O: Observability matrix
"""
# Convert input parameters to matrices (if they aren't already)
@@ -237,3 +237,22 @@
for i in range(1, n):
obsv = np.vstack((obsv, cmat*amat**i))
return obsv
+
+def gram(sys,type):
+ """Gramian
+
+ Usage
+ =====
+ Wc = gram(sys,'c')
+ Wo = gram(sys,'o')
+ """
+
+ if type=='c':
+ print "controllable"
+ elif type=='o':
+ print "observable"
+ else:
+ raise ValueError, "Oops, neither observable, nor controllable!"
+
+ gram = 0.
+ return gram
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 22:04:03
|
Revision: 36
http://python-control.svn.sourceforge.net/python-control/?rev=36&view=rev
Author: kkchen
Date: 2011-02-08 22:03:57 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added test that ctrb and obsv are dual. Also removed accidental extra .hgignore that I shouldn't have added.
Steven Brunton <sbr...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestStatefbk.py
Modified: branches/control-0.4a/src/TestStatefbk.py
===================================================================
--- branches/control-0.4a/src/TestStatefbk.py 2011-02-08 20:57:52 UTC (rev 35)
+++ branches/control-0.4a/src/TestStatefbk.py 2011-02-08 22:03:57 UTC (rev 36)
@@ -32,6 +32,15 @@
Wotrue = N.matrix("5. 6.; 7. 8.; 23. 34.; 31. 46.")
Wo = obsv(A,C)
N.testing.assert_array_almost_equal(Wo, Wotrue)
+
+ def testCtrbObsvDuality(self):
+ A = N.matrix("1.2 -2.3; 3.4 -4.5")
+ B = N.matrix("5.8 6.9; 8. 9.1")
+ Wc = ctrb(A,B);
+ A = N.transpose(A)
+ C = N.transpose(B)
+ Wo = N.transpose(obsv(A,C));
+ N.testing.assert_array_almost_equal(Wc,Wo)
if __name__ == '__main__':
unittest.main()
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 20:57:58
|
Revision: 35
http://python-control.svn.sourceforge.net/python-control/?rev=35&view=rev
Author: kkchen
Date: 2011-02-08 20:57:52 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added TestStatefbk.py, which currently tests obsv and ctrb. Also, obsv and ctrb had their algorithms backwards, so I fixed this.
Modified Paths:
--------------
branches/control-0.4a/src/statefbk.py
Added Paths:
-----------
branches/control-0.4a/src/TestStatefbk.py
Added: branches/control-0.4a/src/TestStatefbk.py
===================================================================
--- branches/control-0.4a/src/TestStatefbk.py (rev 0)
+++ branches/control-0.4a/src/TestStatefbk.py 2011-02-08 20:57:52 UTC (rev 35)
@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+
+from statefbk import ctrb, obsv, place, lqr
+import numpy as N
+import unittest
+
+class TestStatefbk(unittest.TestCase):
+ def testCtrbSISO(self):
+ A = N.matrix("1. 2.; 3. 4.")
+ B = N.matrix("5.; 7.")
+ Wctrue = N.matrix("5. 19.; 7. 43.")
+ Wc = ctrb(A,B)
+ N.testing.assert_array_almost_equal(Wc, Wctrue)
+
+ def testCtrbMIMO(self):
+ A = N.matrix("1. 2.; 3. 4.")
+ B = N.matrix("5. 6.; 7. 8.")
+ Wctrue = N.matrix("5. 6. 19. 22.; 7. 8. 43. 50.")
+ Wc = ctrb(A,B)
+ N.testing.assert_array_almost_equal(Wc, Wctrue)
+
+ def testObsvSISO(self):
+ A = N.matrix("1. 2.; 3. 4.")
+ C = N.matrix("5. 7.")
+ Wotrue = N.matrix("5. 7.; 26. 38.")
+ Wo = obsv(A,C)
+ N.testing.assert_array_almost_equal(Wo, Wotrue)
+
+ def testObsvMIMO(self):
+ A = N.matrix("1. 2.; 3. 4.")
+ C = N.matrix("5. 6.; 7. 8.")
+ Wotrue = N.matrix("5. 6.; 7. 8.; 23. 34.; 31. 46.")
+ Wo = obsv(A,C)
+ N.testing.assert_array_almost_equal(Wo, Wotrue)
+
+if __name__ == '__main__':
+ unittest.main()
Modified: branches/control-0.4a/src/statefbk.py
===================================================================
--- branches/control-0.4a/src/statefbk.py 2011-02-08 18:04:02 UTC (rev 34)
+++ branches/control-0.4a/src/statefbk.py 2011-02-08 20:57:52 UTC (rev 35)
@@ -205,11 +205,10 @@
amat = np.mat(A)
bmat = np.mat(B)
n = np.shape(amat)[0]
-
# Construct the controllability matrix
ctrb = bmat
for i in range(1, n):
- ctrb = np.vstack((ctrb, amat**i*bmat))
+ ctrb = np.hstack((ctrb, amat**i*bmat))
return ctrb
def obsv(A, C):
@@ -236,5 +235,5 @@
# Construct the controllability matrix
obsv = cmat
for i in range(1, n):
- obsv = np.hstack((obsv, cmat*amat**i))
+ obsv = np.vstack((obsv, cmat*amat**i))
return obsv
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <kk...@us...> - 2011-02-08 18:04:11
|
Revision: 34
http://python-control.svn.sourceforge.net/python-control/?rev=34&view=rev
Author: kkchen
Date: 2011-02-08 18:04:02 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Imported SourceForce repo to hg. Created control-0.4a branch.
Added Paths:
-----------
.hgignore
branches/control-0.4a/
branches/control-0.4a/ChangeLog
branches/control-0.4a/MANIFEST.in
branches/control-0.4a/Pending
branches/control-0.4a/README
branches/control-0.4a/build/
branches/control-0.4a/build/lib/
branches/control-0.4a/build/lib/control/
branches/control-0.4a/build/lib/control/TestBDAlg.py
branches/control-0.4a/build/lib/control/TestConvert.py
branches/control-0.4a/build/lib/control/TestFreqRsp.py
branches/control-0.4a/build/lib/control/TestMatlab.py
branches/control-0.4a/build/lib/control/TestModelsimp.py
branches/control-0.4a/build/lib/control/TestSlycot.py
branches/control-0.4a/build/lib/control/TestStateSp.py
branches/control-0.4a/build/lib/control/TestStatefbk.py
branches/control-0.4a/build/lib/control/TestXferFcn.py
branches/control-0.4a/build/lib/control/__init__.py
branches/control-0.4a/build/lib/control/bdalg.py
branches/control-0.4a/build/lib/control/ctrlutil.py
branches/control-0.4a/build/lib/control/delay.py
branches/control-0.4a/build/lib/control/exception.py
branches/control-0.4a/build/lib/control/freqplot.py
branches/control-0.4a/build/lib/control/lti.py
branches/control-0.4a/build/lib/control/matlab.py
branches/control-0.4a/build/lib/control/modelsimp.py
branches/control-0.4a/build/lib/control/pzmap.py
branches/control-0.4a/build/lib/control/rlocus.py
branches/control-0.4a/build/lib/control/robust.py
branches/control-0.4a/build/lib/control/statefbk.py
branches/control-0.4a/build/lib/control/statesp.py
branches/control-0.4a/build/lib/control/test.py
branches/control-0.4a/build/lib/control/xferfcn.py
branches/control-0.4a/doc/
branches/control-0.4a/doc/control.tex
branches/control-0.4a/examples/
branches/control-0.4a/examples/pvtol-lqr.py
branches/control-0.4a/examples/pvtol-nested-ss.py
branches/control-0.4a/examples/pvtol-nested.py
branches/control-0.4a/examples/secord-matlab.py
branches/control-0.4a/examples/slicot-test.py
branches/control-0.4a/examples/type2_type3.py
branches/control-0.4a/external/
branches/control-0.4a/external/controls.py
branches/control-0.4a/external/yottalab.py
branches/control-0.4a/setup.py
branches/control-0.4a/src/
branches/control-0.4a/src/__init__.py
branches/control-0.4a/src/bdalg.py
branches/control-0.4a/src/ctrlutil.py
branches/control-0.4a/src/delay.py
branches/control-0.4a/src/exception.py
branches/control-0.4a/src/freqplot.py
branches/control-0.4a/src/matlab.py
branches/control-0.4a/src/pzmap.py
branches/control-0.4a/src/rlocus.py
branches/control-0.4a/src/statefbk.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Added: .hgignore
===================================================================
--- .hgignore (rev 0)
+++ .hgignore 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,17 @@
+
+# Automatically generated by `hgimportsvn`
+syntax:glob
+.svn
+.hgsvn
+
+# These lines are suggested according to the svn:ignore property
+# Feel free to enable them by uncommenting them
+syntax:glob
+
+.DS_Store
+.*.swp
+build/*
+doc/_build/*
+doc/_templates/*
+doc/_static/*
+*.pyc
Added: branches/control-0.4a/ChangeLog
===================================================================
--- branches/control-0.4a/ChangeLog (rev 0)
+++ branches/control-0.4a/ChangeLog 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,144 @@
+2010-11-05 Richard Murray <murray@sumatra.local>
+
+ * external/yottalab.py: New file containing Roberto Bucher's control
+ library functions. OK to start pulling these into the main library,
+ with attribution, but note that they use modifications of the
+ default library => some rewrites will be needed.
+
+2010-09-11 Richard Murray <murray@sumatra.local>
+
+ * src/matlab.py (step): Added local step response function that uses
+ lsim2() instead of signal.step (which can't handle integrators).
+ This function may not be needed when new scipy step2() function is
+ available.
+ (impulse): Added local impulse response function that sets the
+ initial condition based on the input matrix and then uses the
+ lsim2() function to compute the response.
+
+ * examples/test-response.py: Added test script for making sure that
+ time repsonse functions are working as desired
+
+ * src/matlab.py (lsim): Added local version of lsim that calls
+ signal.lsim2 (actual ODE integrator)
+
+2010-09-06 Richard Murray <murray@sumatra.local>
+
+ * src/statefbk.py (ctrb): new function for testing controllability
+ * src/statefbk.py (obsv): new function for testing observabiilty
+
+2010-09-02 Richard Murray <murray@sumatra.local>
+
+ * src/statefbk.py (place): Use np.size() instead of len() for
+ finding length of placed_eigs for better compatability with
+ different python versions [courtesy of Roberto Bucher]
+
+ * src/delay.py (pade): New file for delay-based computations +
+ initial implementation of pade() [courtesy Sawyer Fuller]
+
+2010-06-17 Richard Murray <murray@sumatra.local>
+
+ * src/rlocus.py: changed num, den to nump, denp for clarity
+ * src/rlocus.py: new file with Ryan Krauss's root locus code
+
+2010-06-06 Richard Murray <murray@sumatra.local>
+
+ * examples/pvtol-lqr.py: Added example to test out LQR routines
+
+ * src/matlab.py (bode): created a wrapper that allows MATLAB style
+ arguments for bode (eg, bode(sys1, sys2))
+
+ * src/ctrlutil.py (issys): added new function to check if an object
+ is a system (state space or transfer function). Will generalize
+ this latter to look for other compatible classes.
+
+ * src/freqplot.py (bode): Compute frequency range of bode plot based
+ on poles and zeros
+ (bode): Allow bode plot to be passed a list (or tuple) as the first
+ argument, in which case multiple bode plots are generated
+
+ * src/statesp.py (StateSpace.zeros): new function to compute zeros
+ for a state space system
+ (StateSpace): defined new functions to compute poles of a state
+ space system
+
+ * src/xferfcn.py (TransferFunction): defined new functions to
+ compute poles and zeros of a transfer function.
+
+2010-05-31 Richard Murray <murray@sumatra.local>
+
+ * src/exception.py (ControlNotImplemented): added new exception, to
+ be used for functions that are not yet implemented
+
+ * src/statefbk.py (lqr): added lqr function (using slycot). Still
+ needs to be verified to make sure calculations are correct.
+
+ * ChangeLog: converted to standard GNU formation (old style below)
+ * setup.py: updated package number to v0.3, changed URL to
+ sourceforge wiki
+
+------------------
+31 May 2010, RMM: added place() function using slycot
+ * New module: statefbk - functions to design state feedback controllers
+ * Uses Enrico Avventi slycot wrappers (http://github.com/avventi/Slycot)
+ * Also added some exception types: ControlSlycot and ControlDimension
+ * Added new example to test slycot interface (directly)
+
+29 May 2010, RMM: updated function documentation
+ * Added __doc__ strings for all current functions
+ * Added __doc__ string to matlab module, listing control toolbox functions
+
+22 May 2010, RMM: tweaked comments and released v0.3a
+ * Changed copyright information on modified files to 2010
+ * Updated "to do" comments to use "#! TODO:" as prefix
+
+11 Feb 2010, GR: implemented and tested state space feedback
+15 Jan 2010, GR: added new example, improved bode
+
+4 Jan 2010, GR: updated bode plots
+ * made bode plot more like matlab
+ * added options for plotting in dB, Hz
+
+27 Dec 2009, GR: important bug fix: feedback TFs were being divided by two
+
+10 Oct 09, RMM: reset matplotlib import in secord-matlab
+ * Using 'from matplotlib import *' causes error with figures
+ * On my other computer, got error when trying to import 'matplotlib.pyplot'
+ * Need to sort out versions and figure out proper import structure
+
+13 Sep 09, RMM: added basic state space functionality
+ * Updated StateSpace routines to allow BD algebra with constants
+ * Updated pvtol-nested example to try to use state space representation
+ (not completely working yet)
+
+12 Sep 09, RMM: code restructuring for transfer functions
+ * Implemented feedback() method in bldalg; partially working
+ for mixture of TF, SS and numbers (not thoroughly tested yet)
+ * New feedback method for TransferFunctions
+ * Updated gangof4 to use new feedback function
+
+9 Sep 09, RMM: updated pzmap to generate a plot (turn off with Plot=False)
+
+8 Sep 09, RMM: rewrite of xferfcn to handle type casting better
+ * Appropriate functions now call convertToTransferFunction
+ * Restricted transfer function to SISO only
+
+7 Sep 09, RMM: additional fixes
+ * Implemented block diagram operations for SISO transfer functions
+ * Modified frequency response functions for transfer functions
+ * Added rudimentary pole/zero computations
+ * Added comments on things that need to be fixed (search for !)
+
+5 Sep 09, RMM: updates to get standard examples working
+ * Copied and converted pvtol_nested.py from AM08, Chapter 11
+ * Updated freqresp to use poly1d for computing values
+ * Added outputs to bode plot - return mag and phase subplot handles
+
+2009-05-24 Richard Murray <murray@kona-2.local>
+ * ./ Initial creation of package files and ChangeLog
+ * Using simpler text format since emacs python mode doesn't support
+ ChangeLog entries anyway
+
+
+Local Variables:
+mode:text
+End:
Added: branches/control-0.4a/MANIFEST.in
===================================================================
--- branches/control-0.4a/MANIFEST.in (rev 0)
+++ branches/control-0.4a/MANIFEST.in 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,5 @@
+include README
+include setup.py
+include src/*.py
+include examples/README examples/*.py
+prune examples/*-test.py
Added: branches/control-0.4a/Pending
===================================================================
--- branches/control-0.4a/Pending (rev 0)
+++ branches/control-0.4a/Pending 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,65 @@
+List of Pending changes for control-python
+RMM, 5 Sep 09
+
+This file contains brief notes on features that need to be added to
+the python control library. Mainly intended to keep track of "bigger
+picture" things that need to be done.
+
+--> See src/matlab.py for a list of MATLAB functions that eventually need
+ to be implemented.
+
+OPEN BUGS
+ * step() doesn't handle systems with a pole at the origin (use lsim2)
+
+Transfer code from Roberto Bucher's yottalab to python-control
+ acker - pole placement using Ackermann method
+ c2d - contimous to discrete time conversion
+ full_obs - full order observer
+ red_obs - reduced order observer
+ comp_form - state feedback controller+observer in compact form
+ comp_form_i - state feedback controller+observer+integ in compact form
+ dsimul - simulate discrete time systems
+ dstep - step response (plot) of discrete time systems
+ dimpulse - imoulse response (plot) of discrete time systems
+ bb_step - step response (plot) of continous time systems
+ sysctr - system+controller+observer+feedback
+ care - Solve Riccati equation for contimous time systems
+ dare - Solve Riccati equation for discrete time systems
+ dlqr - discrete linear quadratic regulator
+ minreal - minimal state space representation
+
+Transfer code from Ryan Krauss's control.py to python-control
+ * phase margin computations (as part of margin command)
+ * step reponse
+ * c2d, c2d_tustin (compare to Bucher version first)
+
+Examples and test cases
+ * 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)
+
+TransferFunction class fixes
+ * evalfr is not working (num, den stored as ndarrays, not poly1ds)
+
+Block diagram algebra fixes
+ * Implement state space block diagram algebra
+ * Produce minimal realizations to avoid later errors
+
+State space class fixes
+ * Convert pvtol to state space systems and rerun example
+ * Implement pzmap for state space systems
+
+LTI updates
+ * Implement control.matlab.step (with semantics similar to MATLAB)
+
+Basic functions to be added
+ * margin - compute gain and phase margin (no plot)
+ * lqr - compute optimal feedback gains (use SLICOT SB02ND.f)
+ * lyap - solve Lyapunov equation (use SLICOT SB03MD.f)
+ * See http://www.slicot.org/shared/libindex.html for list of functions
+
+----
+Instructions for building python package
+ * python setup.py build
+ * python setup.py install
+ * python setup.py sdist
Added: branches/control-0.4a/README
===================================================================
--- branches/control-0.4a/README (rev 0)
+++ branches/control-0.4a/README 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,18 @@
+Python Control System Library
+RMM, 23 May 09
+
+This directory contains the source code for the Python Control Systems
+Library (python-control). This library is still under development,
+but is intended to serve as a wrapper for standard control system
+algorithms in the python programming environment.
+
+Installation instructions
+-------------------------
+Standard python package installation:
+
+ python setup.py install
+
+To see if things are working, you can run the script
+examples/secord-matlab.py (using ipython -pylab). It should generate
+a step response, Bode plot and Nyquist plot for a simple second order
+linear system.
Added: branches/control-0.4a/build/lib/control/TestBDAlg.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestBDAlg.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestBDAlg.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,166 @@
+#!/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()
Added: branches/control-0.4a/build/lib/control/TestConvert.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestConvert.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestConvert.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,91 @@
+#!/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()
Added: branches/control-0.4a/build/lib/control/TestFreqRsp.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestFreqRsp.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestFreqRsp.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,63 @@
+#!/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)
Added: branches/control-0.4a/build/lib/control/TestMatlab.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestMatlab.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestMatlab.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,47 @@
+#!/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()
Added: branches/control-0.4a/build/lib/control/TestModelsimp.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestModelsimp.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestModelsimp.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,96 @@
+#!/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()
Added: branches/control-0.4a/build/lib/control/TestSlycot.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestSlycot.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestSlycot.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+
+#### THIS MUST BE MADE INTO A UNITTEST TO BE PART OF THE TESTING FUNCTIONS!!!!
+
+
+import numpy as np
+from slycot import tb04ad, td04ad
+import matlab
+
+numTests = 1
+maxStates = 3
+maxIO = 3
+
+for states in range(1, maxStates):
+ for inputs in range(1, maxIO):
+ for outputs in range(1, maxIO):
+ sys1 = matlab.rss(states, inputs, outputs)
+ print "sys1"
+ print sys1
+
+ sys2 = tb04ad(states,inputs,outputs,sys1.A,sys1.B,sys1.C,sys1.D,outputs,outputs,inputs)
+ print "sys2"
+ print sys2
+
+ ldwork = 100*max(1,states+max(states,max(3*inputs,3*outputs)))
+ dwork = np.zeros(ldwork)
+ sys3 = td04ad(inputs,outputs,sys2[4],sys2[5],sys2[6])
+ #sys3 = td04ad(inputs,outputs,sys2[4],sys2[5],sys2[6],ldwork)
+ print "sys3"
+ print sys3
+
+ sys4 = tb04ad(states,inputs,outputs,sys3[1][0:states,0:states],sys3[2][0:states,0:inputs],sys3[3][0:outputs,0:states],sys3[4],outputs,outputs,inputs)
+ print "sys4"
+ print sys4
+
+#These are here for once the above is made into a unittest.
+def suite():
+ return unittest.TestLoader().loadTestsFromTestCase(TestSlycot)
+
+if __name__=='__main__':
+ unittest.main()
+
Added: branches/control-0.4a/build/lib/control/TestStateSp.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestStateSp.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestStateSp.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,209 @@
+#!/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()
Added: branches/control-0.4a/build/lib/control/TestStatefbk.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestStatefbk.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestStatefbk.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,88 @@
+#!/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()
Added: branches/control-0.4a/build/lib/control/TestXferFcn.py
===================================================================
--- branches/control-0.4a/build/lib/control/TestXferFcn.py (rev 0)
+++ branches/control-0.4a/build/lib/control/TestXferFcn.py 2011-02-08 18:04:02 UTC (rev 34)
@@ -0,0 +1,437 @@
+#!/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.],
+...
[truncated message content] |
|
From: Richard M. <mu...@cd...> - 2011-02-08 16:44:26
|
This sounds good. I'll try to spend some time looking at the new branch and then quickly switching to control-0.4a as the working version (trunk). -richard On 8 Feb 2011, at 8:39 , Kevin Chen wrote: > Developers, > > We are going to push our changes to a branch called control-0.4a (not control-0.3d, sorry! I hope this is okay. It took a while to get this branch working.). This is my first time using hgsvn, so I hope it works okay. After the push, the trunk should remain unchanged from 0.3c, and the new code should appear in branches/control-0.4a. > > - the Princeton group > > On Tue, Feb 8, 2011 at 11:32, Lauren Padilla <lpa...@pr...> wrote: > > > Begin forwarded message: > > From: Richard Murray <mu...@cd...> > Date: February 6, 2011 5:47:40 PM EST > To: pyt...@li... > Subject: Re: [python-control-discuss] Major updates from the Princeton group > Reply-To: pyt...@li... > > Brandt, Steve, Kevin, and Lauren: > > These look like really great changes and will definitely push the python-control library forward! What I suggest you do is create a copy of the current trunk as a branch (branches/control-0.3d) and then post all of your changes. This will let those of us working on the current trunk use the branch of needed, but hopefully switch to the new structure as soon as possible. > > I suggest that we updated the revision number after your changes, so we will create tags/control-0.4a based on your updated trunk. If you can create a page on the wiki that has the details of your revision (basically the e-mail below), that should help us keep track of what changed. > > As soon as I have a chance, I'll go through and start updating old functions to fit the new documentation on unit test structure, just to make sure I understand how to do this. > > Thanks for all of the hard work! I hope you all continue to contribute changes now and then, despite your class being done. > > -richard > > On 3 Feb 2011, at 14:07 , Lauren Padilla wrote: > > > Richard and developers, > > > > We are writing this note to discuss the major changes that we, the Princeton group, have made to our internal python-control repository. We have more or less completed our contribution to the project, and would like to push our changes to the SVN repository. > > > > == Unfinished business == > > > > - The ss2tf and tf2ss conversions are carried out in SLICOT, wrapped in Slycot, then wrapped in statesp._convertToStateSpace() and xferfcn._convertToTransferFunction(). (As for the user interface, these are then wrapped by tf2ss() and ss2tf(), and finally wrapped by ss() and tf().) The conversions seem to be buggy sometimes. This can be seen by running TestConvert.py, which makes a random state space (sys1), converts it to a transfer function (sys2), back to a state space (sys3), and back to a transfer function (sys4). Sometimes the outputs are correct; other times they are not; and yet other times, our machines enter an infinite loop. We did not have time to solve this issue completely, so it remains open. > > - hsvd() works, but there are better ways to compute the Hankel singular values. Better algorithms, however, may be substantial and difficult to program. Also, we did include MATLAB's baldata output. balred() can return the balancing data. > > - balred() does not yet support MatchDC. > > - h2syn() and hinfsyn() are implemented, but not yet tested. > > - We used Sphinx to create our python-control documentation, as previously discussed on this board. For the most part, however, we only added Sphinx-style documentation to the functions we worked on directly. Some of the docstrings in the library should be updated to this format, if it's desirable to continue using Sphinx. > > - Plotting MIMO systems has not been implemented. In the future, we envision the following common interface for most plotting routines (step, bode, impulse, ...): > > 1. return (i,j,k) array where i = # inputs, j = # outputs, k = # frequencies or times > > 2. plot=True optional argument to actually plot instead of just returning values > > 3. squeeze output so that for SISO system, return is not an (i,j,k), but rather a (k,). > > - It may be desirable to have the TransferFunction class represent SISO numerators and denominators as numpy arrays, instead of lists of lists of arrays. This might be nontrivial to implement. > > > > == LTI Class structure == > > > > As previously discussed on this board, we decided to remove the use of scipy.signal.lti as a parent class. Its restriction to SISO and SIMO systems makes it impractical to use. Therefore, we now have our own structure. The parent class, "Lti," currently only has two fields and nothing else, but this can be expanded in the future. "StateSpace" and "TransferFunction" are now child classes of "Lti," and almost have full MIMO functionality. > > > > In addition to implementing this new class structure, the classes themselves have been developed quite a bit. For instance, both StateSpace.__init__() and TransferFunction.__init__() now have guards to make sure the inputs are valid. In addition, many (but not all) functions have been expanded to support MIMO. For instance, almost all of the operator overloads in TransferFunction were completely rewritten to support MIMO; however, __div__() and __rdiv__() are not yet implemented for MIMO. See the source code for more details. Functions not completed raise NotImplementedError. > > > > == Slycot == > > > > We have wrapped five additional SLICOT functions into Slycot. These Slycot wrappers are used in h2syn(), hinfsyn(), ss2tf(), tf2ss(), and balred(). We are working with Enrico Avventi to merge our changes to his Slycot repository on github. We anticipate some iteration with Enrico before a fully integrated Slycot library is available for download. Until then, we will email a patch and instructions for use by python-control developers and users. > > > > == Unit tests == > > > > The unit tests are packaged, and can be run together using > >>>> import control > >>>> control.tests() > > > > There are a total of nine unit test files (Test*.py). Most of these should be self-explanatory. However, these may not be: > > - TestConvert.py: This file is not finished. Ideally, it would convert random state spaces back and forth between state space and transfer function form, and check that the conversions are correct. However, the conversions are sometimes buggy (see Unfinished Business below). > > - TestFreqResp.py: This is not a unit test, but can be used to check for errors and exceptions in freqresp.py and matlab.bode(). > > - TestSlycot.py: Also not a unit test. This checks for correct functionality of the ss2tf and tf2ss Slycot wrappers tb04ad and td04ad. > > > > == Functions implemented == > > > > This is a brief list of functions we either created or revised significantly. > > > > - StateSpace, TransferFunction > > - As stated above, these two classes were heavily revamped. In addition to these revisions, we added: > > - __init__(): The initializer can now also act as a copy constructor. StateSpace(ss_obj) and TransferFunction(tf_obj) perform deep copies. I believe that Python does not support assignment (=) overloading. Generally speaking, one should never use > >>>> sys1 = sys2 > > since this will be a shallow copy. > > - returnScipySignalLti: scipy.signal.lti is still useful for its few methods. For a MIMO LTI system, this function returns a list of a list of SISO scipy.signal.lti objects--one object for each element of the MIMO transfer function. We use this in lsim, step, initial, and impulse. > > - We also added some "private" methods to these two classes. (E.g., TransferFunction._truncatecoeff(), TransferFunction._common_den(), etc.) In general, these methods are for bookkeeping and do not affect the interface. > > - Random state space generation > > - rss() > > - drss() > > - the generating function _rss_generate() > > - Controllability and observability functions > > - gram() > > - obsv() > > - ctrb() > > - Model reduction tools > > - modred() > > - hsvd() > > - balred() > > - markov() > > - Synthesis tools; these have not been tested > > - h2syn() > > - hinfsyn() > > - Plotting tools > > - matlab.bode() > > - freqresp.bode() > > - default_frequency_range() > > - step() > > - impulse() > > > > == Contact == > > > > If you need to contact us in the future, email us at > > - bb...@pr... > > - sbr...@pr... > > - kk...@pr... > > - lpa...@pr... > > > > Happy coding, > > Brandt, Steve, Kevin, and Lauren > > ------------------------------------------------------------------------------ > > The modern datacenter depends on network connectivity to access resources > > and provide services. The best practices for maximizing a physical server's > > connectivity to a physical network are well understood - see how these > > rules translate into the virtual world? > > http://p.sf.net/sfu/oracle-sfdevnlfb > > _______________________________________________ > > python-control-discuss mailing list > > pyt...@li... > > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > > ------------------------------------------------------------------------------ > The modern datacenter depends on network connectivity to access resources > and provide services. The best practices for maximizing a physical server's > connectivity to a physical network are well understood - see how these > rules translate into the virtual world? > http://p.sf.net/sfu/oracle-sfdevnlfb > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > > ------------------------------------------------------------------------------ > The ultimate all-in-one performance toolkit: Intel(R) Parallel Studio XE: > Pinpoint memory and threading errors before they happen. > Find and fix more than 250 security defects in the development cycle. > Locate bottlenecks in serial and parallel code that limit performance. > http://p.sf.net/sfu/intel-dev2devfeb_______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Kevin C. <kk...@pr...> - 2011-02-08 16:39:49
|
Developers, We are going to push our changes to a branch called control-0.4a (not control-0.3d, sorry! I hope this is okay. It took a while to get this branch working.). This is my first time using hgsvn, so I hope it works okay. After the push, the trunk should remain unchanged from 0.3c, and the new code should appear in branches/control-0.4a. - the Princeton group On Tue, Feb 8, 2011 at 11:32, Lauren Padilla <lpa...@pr...> wrote: > > > Begin forwarded message: > > From: Richard Murray <mu...@cd...> > Date: February 6, 2011 5:47:40 PM EST > To: pyt...@li... > Subject: Re: [python-control-discuss] Major updates from the Princeton > group > Reply-To: pyt...@li... > > Brandt, Steve, Kevin, and Lauren: > > These look like really great changes and will definitely push the > python-control library forward! What I suggest you do is create a copy of > the current trunk as a branch (branches/control-0.3d) and then post all of > your changes. This will let those of us working on the current trunk use > the branch of needed, but hopefully switch to the new structure as soon as > possible. > > I suggest that we updated the revision number after your changes, so we > will create tags/control-0.4a based on your updated trunk. If you can > create a page on the wiki that has the details of your revision (basically > the e-mail below), that should help us keep track of what changed. > > As soon as I have a chance, I'll go through and start updating old > functions to fit the new documentation on unit test structure, just to make > sure I understand how to do this. > > Thanks for all of the hard work! I hope you all continue to contribute > changes now and then, despite your class being done. > > -richard > > On 3 Feb 2011, at 14:07 , Lauren Padilla wrote: > > > Richard and developers, > > > > We are writing this note to discuss the major changes that we, the > Princeton group, have made to our internal python-control repository. We > have more or less completed our contribution to the project, and would like > to push our changes to the SVN repository. > > > > == Unfinished business == > > > > - The ss2tf and tf2ss conversions are carried out in SLICOT, wrapped in > Slycot, then wrapped in statesp._convertToStateSpace() and > xferfcn._convertToTransferFunction(). (As for the user interface, these are > then wrapped by tf2ss() and ss2tf(), and finally wrapped by ss() and tf().) > The conversions seem to be buggy sometimes. This can be seen by running > TestConvert.py, which makes a random state space (sys1), converts it to a > transfer function (sys2), back to a state space (sys3), and back to a > transfer function (sys4). Sometimes the outputs are correct; other times > they are not; and yet other times, our machines enter an infinite loop. We > did not have time to solve this issue completely, so it remains open. > > - hsvd() works, but there are better ways to compute the Hankel singular > values. Better algorithms, however, may be substantial and difficult to > program. Also, we did include MATLAB's baldata output. balred() can return > the balancing data. > > - balred() does not yet support MatchDC. > > - h2syn() and hinfsyn() are implemented, but not yet tested. > > - We used Sphinx to create our python-control documentation, as > previously discussed on this board. For the most part, however, we only > added Sphinx-style documentation to the functions we worked on directly. > Some of the docstrings in the library should be updated to this format, if > it's desirable to continue using Sphinx. > > - Plotting MIMO systems has not been implemented. In the future, we > envision the following common interface for most plotting routines (step, > bode, impulse, ...): > > 1. return (i,j,k) array where i = # inputs, j = # outputs, k = # > frequencies or times > > 2. plot=True optional argument to actually plot instead of just > returning values > > 3. squeeze output so that for SISO system, return is not an (i,j,k), > but rather a (k,). > > - It may be desirable to have the TransferFunction class represent SISO > numerators and denominators as numpy arrays, instead of lists of lists of > arrays. This might be nontrivial to implement. > > > > == LTI Class structure == > > > > As previously discussed on this board, we decided to remove the use of > scipy.signal.lti as a parent class. Its restriction to SISO and SIMO > systems makes it impractical to use. Therefore, we now have our own > structure. The parent class, "Lti," currently only has two fields and > nothing else, but this can be expanded in the future. "StateSpace" and > "TransferFunction" are now child classes of "Lti," and almost have full MIMO > functionality. > > > > In addition to implementing this new class structure, the classes > themselves have been developed quite a bit. For instance, both > StateSpace.__init__() and TransferFunction.__init__() now have guards to > make sure the inputs are valid. In addition, many (but not all) functions > have been expanded to support MIMO. For instance, almost all of the > operator overloads in TransferFunction were completely rewritten to support > MIMO; however, __div__() and __rdiv__() are not yet implemented for MIMO. > See the source code for more details. Functions not completed raise > NotImplementedError. > > > > == Slycot == > > > > We have wrapped five additional SLICOT functions into Slycot. These > Slycot wrappers are used in h2syn(), hinfsyn(), ss2tf(), tf2ss(), and > balred(). We are working with Enrico Avventi to merge our changes to his > Slycot repository on github. We anticipate some iteration with Enrico before > a fully integrated Slycot library is available for download. Until then, we > will email a patch and instructions for use by python-control developers and > users. > > > > == Unit tests == > > > > The unit tests are packaged, and can be run together using > >>>> import control > >>>> control.tests() > > > > There are a total of nine unit test files (Test*.py). Most of these > should be self-explanatory. However, these may not be: > > - TestConvert.py: This file is not finished. Ideally, it would convert > random state spaces back and forth between state space and transfer function > form, and check that the conversions are correct. However, the conversions > are sometimes buggy (see Unfinished Business below). > > - TestFreqResp.py: This is not a unit test, but can be used to check > for errors and exceptions in freqresp.py and matlab.bode(). > > - TestSlycot.py: Also not a unit test. This checks for correct > functionality of the ss2tf and tf2ss Slycot wrappers tb04ad and td04ad. > > > > == Functions implemented == > > > > This is a brief list of functions we either created or revised > significantly. > > > > - StateSpace, TransferFunction > > - As stated above, these two classes were heavily revamped. In > addition to these revisions, we added: > > - __init__(): The initializer can now also act as a copy constructor. > StateSpace(ss_obj) and TransferFunction(tf_obj) perform deep copies. I > believe that Python does not support assignment (=) overloading. Generally > speaking, one should never use > >>>> sys1 = sys2 > > since this will be a shallow copy. > > - returnScipySignalLti: scipy.signal.lti is still useful for its few > methods. For a MIMO LTI system, this function returns a list of a list of > SISO scipy.signal.lti objects--one object for each element of the MIMO > transfer function. We use this in lsim, step, initial, and impulse. > > - We also added some "private" methods to these two classes. (E.g., > TransferFunction._truncatecoeff(), TransferFunction._common_den(), etc.) In > general, these methods are for bookkeeping and do not affect the interface. > > - Random state space generation > > - rss() > > - drss() > > - the generating function _rss_generate() > > - Controllability and observability functions > > - gram() > > - obsv() > > - ctrb() > > - Model reduction tools > > - modred() > > - hsvd() > > - balred() > > - markov() > > - Synthesis tools; these have not been tested > > - h2syn() > > - hinfsyn() > > - Plotting tools > > - matlab.bode() > > - freqresp.bode() > > - default_frequency_range() > > - step() > > - impulse() > > > > == Contact == > > > > If you need to contact us in the future, email us at > > - bb...@pr... > > - sbr...@pr... > > - kk...@pr... > > - lpa...@pr... > > > > Happy coding, > > Brandt, Steve, Kevin, and Lauren > > > ------------------------------------------------------------------------------ > > The modern datacenter depends on network connectivity to access resources > > and provide services. The best practices for maximizing a physical > server's > > connectivity to a physical network are well understood - see how these > > rules translate into the virtual world? > > http://p.sf.net/sfu/oracle-sfdevnlfb > > _______________________________________________ > > python-control-discuss mailing list > > pyt...@li... > > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > > > ------------------------------------------------------------------------------ > The modern datacenter depends on network connectivity to access resources > and provide services. The best practices for maximizing a physical server's > connectivity to a physical network are well understood - see how these > rules translate into the virtual world? > http://p.sf.net/sfu/oracle-sfdevnlfb > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > |
|
From: Richard M. <mu...@cd...> - 2011-02-06 22:47:49
|
Brandt, Steve, Kevin, and Lauren: These look like really great changes and will definitely push the python-control library forward! What I suggest you do is create a copy of the current trunk as a branch (branches/control-0.3d) and then post all of your changes. This will let those of us working on the current trunk use the branch of needed, but hopefully switch to the new structure as soon as possible. I suggest that we updated the revision number after your changes, so we will create tags/control-0.4a based on your updated trunk. If you can create a page on the wiki that has the details of your revision (basically the e-mail below), that should help us keep track of what changed. As soon as I have a chance, I'll go through and start updating old functions to fit the new documentation on unit test structure, just to make sure I understand how to do this. Thanks for all of the hard work! I hope you all continue to contribute changes now and then, despite your class being done. -richard On 3 Feb 2011, at 14:07 , Lauren Padilla wrote: > Richard and developers, > > We are writing this note to discuss the major changes that we, the Princeton group, have made to our internal python-control repository. We have more or less completed our contribution to the project, and would like to push our changes to the SVN repository. > > == Unfinished business == > > - The ss2tf and tf2ss conversions are carried out in SLICOT, wrapped in Slycot, then wrapped in statesp._convertToStateSpace() and xferfcn._convertToTransferFunction(). (As for the user interface, these are then wrapped by tf2ss() and ss2tf(), and finally wrapped by ss() and tf().) The conversions seem to be buggy sometimes. This can be seen by running TestConvert.py, which makes a random state space (sys1), converts it to a transfer function (sys2), back to a state space (sys3), and back to a transfer function (sys4). Sometimes the outputs are correct; other times they are not; and yet other times, our machines enter an infinite loop. We did not have time to solve this issue completely, so it remains open. > - hsvd() works, but there are better ways to compute the Hankel singular values. Better algorithms, however, may be substantial and difficult to program. Also, we did include MATLAB's baldata output. balred() can return the balancing data. > - balred() does not yet support MatchDC. > - h2syn() and hinfsyn() are implemented, but not yet tested. > - We used Sphinx to create our python-control documentation, as previously discussed on this board. For the most part, however, we only added Sphinx-style documentation to the functions we worked on directly. Some of the docstrings in the library should be updated to this format, if it's desirable to continue using Sphinx. > - Plotting MIMO systems has not been implemented. In the future, we envision the following common interface for most plotting routines (step, bode, impulse, ...): > 1. return (i,j,k) array where i = # inputs, j = # outputs, k = # frequencies or times > 2. plot=True optional argument to actually plot instead of just returning values > 3. squeeze output so that for SISO system, return is not an (i,j,k), but rather a (k,). > - It may be desirable to have the TransferFunction class represent SISO numerators and denominators as numpy arrays, instead of lists of lists of arrays. This might be nontrivial to implement. > > == LTI Class structure == > > As previously discussed on this board, we decided to remove the use of scipy.signal.lti as a parent class. Its restriction to SISO and SIMO systems makes it impractical to use. Therefore, we now have our own structure. The parent class, "Lti," currently only has two fields and nothing else, but this can be expanded in the future. "StateSpace" and "TransferFunction" are now child classes of "Lti," and almost have full MIMO functionality. > > In addition to implementing this new class structure, the classes themselves have been developed quite a bit. For instance, both StateSpace.__init__() and TransferFunction.__init__() now have guards to make sure the inputs are valid. In addition, many (but not all) functions have been expanded to support MIMO. For instance, almost all of the operator overloads in TransferFunction were completely rewritten to support MIMO; however, __div__() and __rdiv__() are not yet implemented for MIMO. See the source code for more details. Functions not completed raise NotImplementedError. > > == Slycot == > > We have wrapped five additional SLICOT functions into Slycot. These Slycot wrappers are used in h2syn(), hinfsyn(), ss2tf(), tf2ss(), and balred(). We are working with Enrico Avventi to merge our changes to his Slycot repository on github. We anticipate some iteration with Enrico before a fully integrated Slycot library is available for download. Until then, we will email a patch and instructions for use by python-control developers and users. > > == Unit tests == > > The unit tests are packaged, and can be run together using >>>> import control >>>> control.tests() > > There are a total of nine unit test files (Test*.py). Most of these should be self-explanatory. However, these may not be: > - TestConvert.py: This file is not finished. Ideally, it would convert random state spaces back and forth between state space and transfer function form, and check that the conversions are correct. However, the conversions are sometimes buggy (see Unfinished Business below). > - TestFreqResp.py: This is not a unit test, but can be used to check for errors and exceptions in freqresp.py and matlab.bode(). > - TestSlycot.py: Also not a unit test. This checks for correct functionality of the ss2tf and tf2ss Slycot wrappers tb04ad and td04ad. > > == Functions implemented == > > This is a brief list of functions we either created or revised significantly. > > - StateSpace, TransferFunction > - As stated above, these two classes were heavily revamped. In addition to these revisions, we added: > - __init__(): The initializer can now also act as a copy constructor. StateSpace(ss_obj) and TransferFunction(tf_obj) perform deep copies. I believe that Python does not support assignment (=) overloading. Generally speaking, one should never use >>>> sys1 = sys2 > since this will be a shallow copy. > - returnScipySignalLti: scipy.signal.lti is still useful for its few methods. For a MIMO LTI system, this function returns a list of a list of SISO scipy.signal.lti objects--one object for each element of the MIMO transfer function. We use this in lsim, step, initial, and impulse. > - We also added some "private" methods to these two classes. (E.g., TransferFunction._truncatecoeff(), TransferFunction._common_den(), etc.) In general, these methods are for bookkeeping and do not affect the interface. > - Random state space generation > - rss() > - drss() > - the generating function _rss_generate() > - Controllability and observability functions > - gram() > - obsv() > - ctrb() > - Model reduction tools > - modred() > - hsvd() > - balred() > - markov() > - Synthesis tools; these have not been tested > - h2syn() > - hinfsyn() > - Plotting tools > - matlab.bode() > - freqresp.bode() > - default_frequency_range() > - step() > - impulse() > > == Contact == > > If you need to contact us in the future, email us at > - bb...@pr... > - sbr...@pr... > - kk...@pr... > - lpa...@pr... > > Happy coding, > Brandt, Steve, Kevin, and Lauren > ------------------------------------------------------------------------------ > The modern datacenter depends on network connectivity to access resources > and provide services. The best practices for maximizing a physical server's > connectivity to a physical network are well understood - see how these > rules translate into the virtual world? > http://p.sf.net/sfu/oracle-sfdevnlfb > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Scott C. L. <sli...@ca...> - 2011-02-04 00:43:52
|
Hi Lauren and collaborators, Are you waiting for approval before committing? If so, can I see your work in the meantime, e.g. as a snapshot of your internal repo? Thanks, Scott Livingston On 3 February 2011 14:07, Lauren Padilla <lpa...@pr...> wrote: > Richard and developers, > > We are writing this note to discuss the major changes that we, the Princeton group, have made to our internal python-control repository. We have more or less completed our contribution to the project, and would like to push our changes to the SVN repository. > > == Unfinished business == > > - The ss2tf and tf2ss conversions are carried out in SLICOT, wrapped in Slycot, then wrapped in statesp._convertToStateSpace() and xferfcn._convertToTransferFunction(). (As for the user interface, these are then wrapped by tf2ss() and ss2tf(), and finally wrapped by ss() and tf().) The conversions seem to be buggy sometimes. This can be seen by running TestConvert.py, which makes a random state space (sys1), converts it to a transfer function (sys2), back to a state space (sys3), and back to a transfer function (sys4). Sometimes the outputs are correct; other times they are not; and yet other times, our machines enter an infinite loop. We did not have time to solve this issue completely, so it remains open. > - hsvd() works, but there are better ways to compute the Hankel singular values. Better algorithms, however, may be substantial and difficult to program. Also, we did include MATLAB's baldata output. balred() can return the balancing data. > - balred() does not yet support MatchDC. > - h2syn() and hinfsyn() are implemented, but not yet tested. > - We used Sphinx to create our python-control documentation, as previously discussed on this board. For the most part, however, we only added Sphinx-style documentation to the functions we worked on directly. Some of the docstrings in the library should be updated to this format, if it's desirable to continue using Sphinx. > - Plotting MIMO systems has not been implemented. In the future, we envision the following common interface for most plotting routines (step, bode, impulse, ...): > 1. return (i,j,k) array where i = # inputs, j = # outputs, k = # frequencies or times > 2. plot=True optional argument to actually plot instead of just returning values > 3. squeeze output so that for SISO system, return is not an (i,j,k), but rather a (k,). > - It may be desirable to have the TransferFunction class represent SISO numerators and denominators as numpy arrays, instead of lists of lists of arrays. This might be nontrivial to implement. > > == LTI Class structure == > > As previously discussed on this board, we decided to remove the use of scipy.signal.lti as a parent class. Its restriction to SISO and SIMO systems makes it impractical to use. Therefore, we now have our own structure. The parent class, "Lti," currently only has two fields and nothing else, but this can be expanded in the future. "StateSpace" and "TransferFunction" are now child classes of "Lti," and almost have full MIMO functionality. > > In addition to implementing this new class structure, the classes themselves have been developed quite a bit. For instance, both StateSpace.__init__() and TransferFunction.__init__() now have guards to make sure the inputs are valid. In addition, many (but not all) functions have been expanded to support MIMO. For instance, almost all of the operator overloads in TransferFunction were completely rewritten to support MIMO; however, __div__() and __rdiv__() are not yet implemented for MIMO. See the source code for more details. Functions not completed raise NotImplementedError. > > == Slycot == > > We have wrapped five additional SLICOT functions into Slycot. These Slycot wrappers are used in h2syn(), hinfsyn(), ss2tf(), tf2ss(), and balred(). We are working with Enrico Avventi to merge our changes to his Slycot repository on github. We anticipate some iteration with Enrico before a fully integrated Slycot library is available for download. Until then, we will email a patch and instructions for use by python-control developers and users. > > == Unit tests == > > The unit tests are packaged, and can be run together using > >>> import control > >>> control.tests() > > There are a total of nine unit test files (Test*.py). Most of these should be self-explanatory. However, these may not be: > - TestConvert.py: This file is not finished. Ideally, it would convert random state spaces back and forth between state space and transfer function form, and check that the conversions are correct. However, the conversions are sometimes buggy (see Unfinished Business below). > - TestFreqResp.py: This is not a unit test, but can be used to check for errors and exceptions in freqresp.py and matlab.bode(). > - TestSlycot.py: Also not a unit test. This checks for correct functionality of the ss2tf and tf2ss Slycot wrappers tb04ad and td04ad. > > == Functions implemented == > > This is a brief list of functions we either created or revised significantly. > > - StateSpace, TransferFunction > - As stated above, these two classes were heavily revamped. In addition to these revisions, we added: > - __init__(): The initializer can now also act as a copy constructor. StateSpace(ss_obj) and TransferFunction(tf_obj) perform deep copies. I believe that Python does not support assignment (=) overloading. Generally speaking, one should never use > >>> sys1 = sys2 > since this will be a shallow copy. > - returnScipySignalLti: scipy.signal.lti is still useful for its few methods. For a MIMO LTI system, this function returns a list of a list of SISO scipy.signal.lti objects--one object for each element of the MIMO transfer function. We use this in lsim, step, initial, and impulse. > - We also added some "private" methods to these two classes. (E.g., TransferFunction._truncatecoeff(), TransferFunction._common_den(), etc.) In general, these methods are for bookkeeping and do not affect the interface. > - Random state space generation > - rss() > - drss() > - the generating function _rss_generate() > - Controllability and observability functions > - gram() > - obsv() > - ctrb() > - Model reduction tools > - modred() > - hsvd() > - balred() > - markov() > - Synthesis tools; these have not been tested > - h2syn() > - hinfsyn() > - Plotting tools > - matlab.bode() > - freqresp.bode() > - default_frequency_range() > - step() > - impulse() > > == Contact == > > If you need to contact us in the future, email us at > - bb...@pr... > - sbr...@pr... > - kk...@pr... > - lpa...@pr... > > Happy coding, > Brandt, Steve, Kevin, and Lauren > ------------------------------------------------------------------------------ > The modern datacenter depends on network connectivity to access resources > and provide services. The best practices for maximizing a physical server's > connectivity to a physical network are well understood - see how these > rules translate into the virtual world? > http://p.sf.net/sfu/oracle-sfdevnlfb > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > -- http://scottman.net 865-964-5384 |
|
From: Lauren P. <lpa...@pr...> - 2011-02-03 22:07:35
|
Richard and developers,
We are writing this note to discuss the major changes that we, the Princeton group, have made to our internal python-control repository. We have more or less completed our contribution to the project, and would like to push our changes to the SVN repository.
== Unfinished business ==
- The ss2tf and tf2ss conversions are carried out in SLICOT, wrapped in Slycot, then wrapped in statesp._convertToStateSpace() and xferfcn._convertToTransferFunction(). (As for the user interface, these are then wrapped by tf2ss() and ss2tf(), and finally wrapped by ss() and tf().) The conversions seem to be buggy sometimes. This can be seen by running TestConvert.py, which makes a random state space (sys1), converts it to a transfer function (sys2), back to a state space (sys3), and back to a transfer function (sys4). Sometimes the outputs are correct; other times they are not; and yet other times, our machines enter an infinite loop. We did not have time to solve this issue completely, so it remains open.
- hsvd() works, but there are better ways to compute the Hankel singular values. Better algorithms, however, may be substantial and difficult to program. Also, we did include MATLAB's baldata output. balred() can return the balancing data.
- balred() does not yet support MatchDC.
- h2syn() and hinfsyn() are implemented, but not yet tested.
- We used Sphinx to create our python-control documentation, as previously discussed on this board. For the most part, however, we only added Sphinx-style documentation to the functions we worked on directly. Some of the docstrings in the library should be updated to this format, if it's desirable to continue using Sphinx.
- Plotting MIMO systems has not been implemented. In the future, we envision the following common interface for most plotting routines (step, bode, impulse, ...):
1. return (i,j,k) array where i = # inputs, j = # outputs, k = # frequencies or times
2. plot=True optional argument to actually plot instead of just returning values
3. squeeze output so that for SISO system, return is not an (i,j,k), but rather a (k,).
- It may be desirable to have the TransferFunction class represent SISO numerators and denominators as numpy arrays, instead of lists of lists of arrays. This might be nontrivial to implement.
== LTI Class structure ==
As previously discussed on this board, we decided to remove the use of scipy.signal.lti as a parent class. Its restriction to SISO and SIMO systems makes it impractical to use. Therefore, we now have our own structure. The parent class, "Lti," currently only has two fields and nothing else, but this can be expanded in the future. "StateSpace" and "TransferFunction" are now child classes of "Lti," and almost have full MIMO functionality.
In addition to implementing this new class structure, the classes themselves have been developed quite a bit. For instance, both StateSpace.__init__() and TransferFunction.__init__() now have guards to make sure the inputs are valid. In addition, many (but not all) functions have been expanded to support MIMO. For instance, almost all of the operator overloads in TransferFunction were completely rewritten to support MIMO; however, __div__() and __rdiv__() are not yet implemented for MIMO. See the source code for more details. Functions not completed raise NotImplementedError.
== Slycot ==
We have wrapped five additional SLICOT functions into Slycot. These Slycot wrappers are used in h2syn(), hinfsyn(), ss2tf(), tf2ss(), and balred(). We are working with Enrico Avventi to merge our changes to his Slycot repository on github. We anticipate some iteration with Enrico before a fully integrated Slycot library is available for download. Until then, we will email a patch and instructions for use by python-control developers and users.
== Unit tests ==
The unit tests are packaged, and can be run together using
>>> import control
>>> control.tests()
There are a total of nine unit test files (Test*.py). Most of these should be self-explanatory. However, these may not be:
- TestConvert.py: This file is not finished. Ideally, it would convert random state spaces back and forth between state space and transfer function form, and check that the conversions are correct. However, the conversions are sometimes buggy (see Unfinished Business below).
- TestFreqResp.py: This is not a unit test, but can be used to check for errors and exceptions in freqresp.py and matlab.bode().
- TestSlycot.py: Also not a unit test. This checks for correct functionality of the ss2tf and tf2ss Slycot wrappers tb04ad and td04ad.
== Functions implemented ==
This is a brief list of functions we either created or revised significantly.
- StateSpace, TransferFunction
- As stated above, these two classes were heavily revamped. In addition to these revisions, we added:
- __init__(): The initializer can now also act as a copy constructor. StateSpace(ss_obj) and TransferFunction(tf_obj) perform deep copies. I believe that Python does not support assignment (=) overloading. Generally speaking, one should never use
>>> sys1 = sys2
since this will be a shallow copy.
- returnScipySignalLti: scipy.signal.lti is still useful for its few methods. For a MIMO LTI system, this function returns a list of a list of SISO scipy.signal.lti objects--one object for each element of the MIMO transfer function. We use this in lsim, step, initial, and impulse.
- We also added some "private" methods to these two classes. (E.g., TransferFunction._truncatecoeff(), TransferFunction._common_den(), etc.) In general, these methods are for bookkeeping and do not affect the interface.
- Random state space generation
- rss()
- drss()
- the generating function _rss_generate()
- Controllability and observability functions
- gram()
- obsv()
- ctrb()
- Model reduction tools
- modred()
- hsvd()
- balred()
- markov()
- Synthesis tools; these have not been tested
- h2syn()
- hinfsyn()
- Plotting tools
- matlab.bode()
- freqresp.bode()
- default_frequency_range()
- step()
- impulse()
== Contact ==
If you need to contact us in the future, email us at
- bb...@pr...
- sbr...@pr...
- kk...@pr...
- lpa...@pr...
Happy coding,
Brandt, Steve, Kevin, and Lauren
|
|
From: Richard M. <mu...@cd...> - 2010-12-18 03:33:41
|
A couple of comments on your suggested changes (hopefully others will chime in as well): >> 1. MIMO - The most important item is a new LTI base class to support MIMO (and StateSpace and TransferFunction derived classes). This is discussed in the document lti_overhaul.pdf, and we have a skeleton of the python code, skeleton.py. The basic approach sounds fine, but I would suggest having the class named control.lti since I don't think python is going to get confused between that class and signal.lti. To some extent all of this is hidden from the end user, so probably doesn't matter much (just aesthetically I never like new classes that use '2' as the modifier). Other possibilities would be MIMO, LtiMIMO, etc. To Roberto's point, I think we should also think through how we want to do discrete time systems. Comments (and alternatives) on the class structure I outline would be great. I would eventually like to allow for time delayed systems as well (similar to MATLAB). >> 2. Documentation - Python, NumPy, and SciPy use Sphinx to generate the documentation found at their respective websites. For purposes of consistency, we will use the same. An example is provided in html.zip Sounds right to me. >> 3. New Functions and Unit Tests - We have implemented the gram, rss, and hsvd functions. We have also implemented unit tests for these functions as well as ctrb and obsv using the Python unittest framework. Next we are planning on wrapping SLICOT routines to add balred and modred functions. Good. We need to go through and start implementing unit tests for all of the other basic functions, so we can make sure that we have fairly comprehensive testing built into the package. >> 4. SVN Merge - After getting your feedback, we would like to implement the new MIMO Lti2 class. Once we have all of the wrinkles ironed out, we intend to make a merge back into the SVN repository. Until then, we are doing internal version control using Mercurial. Let me know what you think about the class name, but I think it is OK to commit when you are ready. If you are going to break some functions due to your changes, we should create a branch and update things there, so that trunk is always a working copy. -richard |
|
From: Richard M. <mu...@cd...> - 2010-12-17 19:32:49
|
What I discussed with the Princeton group when they visited me in Pasadena was actually defining a separate class for discrete time systems, rather than having a field that indicated the type of system. This seems cleaner and with the late binding features in python should let existing functions that make sense in discrete time work without change. So we might have a class structure that looks like
* Lti2
* StateSpace
* DTStateSpace
* CTStateSpace
* TDStateSpace (time-delayed systems)
* TransferFunction
* DTTransferFunction
* CTTransferFunction
Of course, many functions would just work on the StateSpace and TransferFunction classes directly (or even LTI2, for things that are really generic).
More comments on the original proposals later in the day.
-richard
On 17 Dec 2010, at 6:58 , Roberto Bucher wrote:
> Please don't forget to implement a field to handle the sampling time for
> managing continous time and discrete time systems (for example Tsamp=0->
> continous time system, Tsamp > 0 =Ts -> discrete time system)
>
> Regards
>
> Roberto
>
>
> On Friday 17 December 2010 13:35:33 Richard Murray wrote:
>> Steve, Kevin and Lauren sent me this update a while back, but it looks like
>> my forward to the list didn't go through. Sending again, so that I (and
>> others) can respond to it.
>>
>> -richard
>>
>> Begin forwarded message:
>>> From: Steve Brunton <sbr...@Pr...>
>>> Date: 8 December 2010 8:08:41 PST
>>> To: Richard Murray <mu...@cd...>
>>> Cc: Kevin Chen <kk...@Pr...>, Lauren Padilla
>>> <lpa...@Pr...>, Steve Brunton <sbr...@Pr...>, Brandt
>>> Belson <bb...@Pr...> Subject: python-control update
>>>
>>> Hi Richard,
>>>
>>> We wanted to give you an update on our work with python-control.
>>>
>>> 1. MIMO - The most important item is a new LTI base class to support
>>> MIMO (and StateSpace and TransferFunction derived classes). This is
>>> discussed in the document lti_overhaul.pdf, and we have a skeleton of
>>> the python code, skeleton.py.
>>>
>>> 2. Documentation - Python, NumPy, and SciPy use Sphinx to generate the
>>> documentation found at their respective websites. For purposes of
>>> consistency, we will use the same. An example is provided in html.zip
>>>
>>> 3. New Functions and Unit Tests - We have implemented the gram, rss,
>>> and hsvd functions. We have also implemented unit tests for these
>>> functions as well as ctrb and obsv using the Python unittest framework.
>>> Next we are planning on wrapping SLICOT routines to add balred and
>>> modred functions.
>>>
>>> 4. SVN Merge - After getting your feedback, we would like to implement
>>> the new MIMO Lti2 class. Once we have all of the wrinkles ironed out,
>>> we intend to make a merge back into the SVN repository. Until then, we
>>> are doing internal version control using Mercurial.
>>>
>>> Thanks, and we look forward to getting your feedback.
>>>
>>> Steve, Kevin, and Lauren
>
> --
> -----------------------------------------------------------------------------
> Great spirits have always encountered violent opposition
> from mediocre minds (A. Einstein)
> ----------------------------------------------------------------------------
> University of Applied Sciences of Southern Switzerland
> Dept. Innovative Technologies
> CH-6928 Lugano-Manno
> http://web.dti.supsi.ch/~bucher
|
|
From: Roberto B. <rob...@su...> - 2010-12-17 14:58:53
|
Please don't forget to implement a field to handle the sampling time for managing continous time and discrete time systems (for example Tsamp=0-> continous time system, Tsamp > 0 =Ts -> discrete time system) Regards Roberto On Friday 17 December 2010 13:35:33 Richard Murray wrote: > Steve, Kevin and Lauren sent me this update a while back, but it looks like > my forward to the list didn't go through. Sending again, so that I (and > others) can respond to it. > > -richard > > Begin forwarded message: > > From: Steve Brunton <sbr...@Pr...> > > Date: 8 December 2010 8:08:41 PST > > To: Richard Murray <mu...@cd...> > > Cc: Kevin Chen <kk...@Pr...>, Lauren Padilla > > <lpa...@Pr...>, Steve Brunton <sbr...@Pr...>, Brandt > > Belson <bb...@Pr...> Subject: python-control update > > > > Hi Richard, > > > > We wanted to give you an update on our work with python-control. > > > > 1. MIMO - The most important item is a new LTI base class to support > > MIMO (and StateSpace and TransferFunction derived classes). This is > > discussed in the document lti_overhaul.pdf, and we have a skeleton of > > the python code, skeleton.py. > > > > 2. Documentation - Python, NumPy, and SciPy use Sphinx to generate the > > documentation found at their respective websites. For purposes of > > consistency, we will use the same. An example is provided in html.zip > > > > 3. New Functions and Unit Tests - We have implemented the gram, rss, > > and hsvd functions. We have also implemented unit tests for these > > functions as well as ctrb and obsv using the Python unittest framework. > > Next we are planning on wrapping SLICOT routines to add balred and > > modred functions. > > > > 4. SVN Merge - After getting your feedback, we would like to implement > > the new MIMO Lti2 class. Once we have all of the wrinkles ironed out, > > we intend to make a merge back into the SVN repository. Until then, we > > are doing internal version control using Mercurial. > > > > Thanks, and we look forward to getting your feedback. > > > > Steve, Kevin, and Lauren -- ----------------------------------------------------------------------------- Great spirits have always encountered violent opposition from mediocre minds (A. Einstein) ---------------------------------------------------------------------------- University of Applied Sciences of Southern Switzerland Dept. Innovative Technologies CH-6928 Lugano-Manno http://web.dti.supsi.ch/~bucher |
|
From: Richard M. <mu...@cd...> - 2010-12-17 12:35:45
|
Steve, Kevin and Lauren sent me this update a while back, but it looks like my forward to the list didn't go through. Sending again, so that I (and others) can respond to it. -richard Begin forwarded message: > From: Steve Brunton <sbr...@Pr...> > Date: 8 December 2010 8:08:41 PST > To: Richard Murray <mu...@cd...> > Cc: Kevin Chen <kk...@Pr...>, Lauren Padilla <lpa...@Pr...>, Steve Brunton <sbr...@Pr...>, Brandt Belson <bb...@Pr...> > Subject: python-control update > > Hi Richard, > > We wanted to give you an update on our work with python-control. > > 1. MIMO - The most important item is a new LTI base class to support MIMO (and StateSpace and TransferFunction derived classes). This is discussed in the document lti_overhaul.pdf, and we have a skeleton of the python code, skeleton.py. > > 2. Documentation - Python, NumPy, and SciPy use Sphinx to generate the documentation found at their respective websites. For purposes of consistency, we will use the same. An example is provided in html.zip > > 3. New Functions and Unit Tests - We have implemented the gram, rss, and hsvd functions. We have also implemented unit tests for these functions as well as ctrb and obsv using the Python unittest framework. Next we are planning on wrapping SLICOT routines to add balred and modred functions. > > 4. SVN Merge - After getting your feedback, we would like to implement the new MIMO Lti2 class. Once we have all of the wrinkles ironed out, we intend to make a merge back into the SVN repository. Until then, we are doing internal version control using Mercurial. > > Thanks, and we look forward to getting your feedback. > > Steve, Kevin, and Lauren > |
|
From: Steve B. <sbr...@Pr...> - 2010-12-08 18:14:41
|
Hi Richard, We wanted to give you an update on our work with python-control. 1. MIMO - The most important item is a new LTI base class to support MIMO (and StateSpace and TransferFunction derived classes). This is discussed in the document lti_overhaul.pdf, and we have a skeleton of the python code, skeleton.py. 2. Documentation - Python, NumPy, and SciPy use Sphinx to generate the documentation found at their respective websites. For purposes of consistency, we will use the same. An example is provided in html.zip 3. New Functions and Unit Tests - We have implemented the gram, rss, and hsvd functions. We have also implemented unit tests for these functions as well as ctrb and obsv using the Python unittest framework. Next we are planning on wrapping SLICOT routines to add balred and modred functions. 4. SVN Merge - After getting your feedback, we would like to implement the new MIMO Lti2 class. Once we have all of the wrinkles ironed out, we intend to make a merge back into the SVN repository. Until then, we are doing internal version control using Mercurial. Thanks, and we look forward to getting your feedback. Steve, Kevin, and Lauren |
|
From: <mur...@us...> - 2010-11-26 22:00:03
|
Revision: 33
http://python-control.svn.sourceforge.net/python-control/?rev=33&view=rev
Author: murrayrm
Date: 2010-11-26 21:59:57 +0000 (Fri, 26 Nov 2010)
Log Message:
-----------
* Added Nichols plot, contributed by Allan McInnes
* Allan also split out code to determine default frequency range
Modified Paths:
--------------
trunk/src/freqplot.py
trunk/src/matlab.py
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2010-11-26 21:56:26 UTC (rev 32)
+++ trunk/src/freqplot.py 2010-11-26 21:59:57 UTC (rev 33)
@@ -4,7 +4,7 @@
# Date: 24 May 09
#
# This file contains some standard control system plots: Bode plots,
-# Nyquist plots and pole-zero diagrams
+# Nyquist plots, Nichols plots and pole-zero diagrams
#
# Copyright (c) 2010 by California Institute of Technology
# All rights reserved.
@@ -46,6 +46,13 @@
from ctrlutil import unwrap
from bdalg import feedback
+#
+# Main plotting functions
+#
+# This section of the code contains the functions for generating
+# frequency domain plots
+#
+
# Bode plot
def bode(syslist, omega=None, dB=False, Hz=False):
"""Bode plot for a system
@@ -81,36 +88,10 @@
if (not getattr(syslist, '__iter__', False)):
syslist = (syslist,)
- #
# Select a default range if none is provided
- #
- # This code looks at the poles and zeros of all of the systems that
- # we are plotting and sets the frequency range to be one decade above
- # and below the min and max feature frequencies, rounded to the nearest
- # integer. It excludes poles and zeros at the origin. If no features
- # are found, it turns logspace(-1, 1)
- #
if (omega == None):
- # Find the list of all poles and zeros in the systems
- features = np.array(())
- for sys in syslist:
- # Add new features to the list
- features = np.concatenate((features, np.abs(sys.poles)))
- features = np.concatenate((features, np.abs(sys.zeros)))
+ omega = default_frequency_range(syslist)
- # Get rid of poles and zeros at the origin
- features = features[features != 0];
-
- # Make sure there is at least one point in the range
- if (features.shape[0] == 0): features = [1];
-
- # Take the log of the features
- features = np.log10(features)
-
- # Set the to be an order of magnitude beyond any features
- omega = sp.logspace(np.floor(np.min(features))-1,
- np.ceil(np.max(features))+1)
-
for sys in syslist:
# Get the magnitude and phase of the system
mag, phase, omega = sys.freqresp(omega)
@@ -154,7 +135,7 @@
return (211, 212)
# Nyquist plot
-def nyquist(sys, omega=None):
+def nyquist(syslist, omega=None):
"""Nyquist plot for a system
Usage
@@ -165,8 +146,8 @@
Parameters
----------
- sys : linsys
- Linear input/output system
+ syslist : linsys
+ List of linear input/output systems (single system is OK)
omega : freq_range
Range of frequencies (list or bounds) in rad/sec
@@ -174,27 +155,82 @@
-------------
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
- #! TODO: This needs to be made more intelligent
if (omega == None):
- omega = sp.logspace(-2, 2);
+ omega = default_frequency_range(syslist)
- # Get the magnitude and phase of the system
- mag, phase, omega = sys.freqresp(omega)
+ for sys in syslist:
+ # Get the magnitude and phase of the system
+ mag, phase, omega = sys.freqresp(omega)
+
+ # Compute the primary curve
+ x = sp.multiply(mag, sp.cos(phase));
+ y = sp.multiply(mag, sp.sin(phase));
+
+ # Plot the primary curve and mirror image
+ plt.plot(x, y, '-');
+ plt.plot(x, -y, '--');
- # Compute the primary curve
- x = sp.multiply(mag, sp.cos(phase));
- y = sp.multiply(mag, sp.sin(phase));
+ # Mark the -1 point
+ plt.plot([-1], [0], 'r+')
- # Plot the primary curve and mirror image
- plt.plot(x, y, '-');
- plt.plot(x, -y, '--');
+# Nichols plot
+# Contributed by Allan McInnes <All...@ca...>
+#! TODO: need unit test code
+def nichols(syslist, omega=None):
+ """Nichols plot for a system
- # Mark the -1 point
- plt.plot([-1], [0], '+k')
+ 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+')
+
# Gang of Four
+#! TODO: think about how (and whether) to handle lists of systems
def gangof4(P, C, omega=None):
"""Plot the "Gang of 4" transfer functions for a system
@@ -220,7 +256,7 @@
# Select a default range if none is provided
#! TODO: This needs to be made more intelligent
if (omega == None):
- omega = sp.logspace(-2, 2);
+ omega = default_frequency_range((P,C))
# Compute the senstivity functions
L = P*C;
@@ -240,3 +276,60 @@
mag, phase, omega = S.freqresp(omega);
plt.subplot(224); plt.loglog(omega, mag);
+
+#
+# Utility functions
+#
+# This section of the code contains some utility functions for
+# generating frequency domain plots
+#
+
+# Compute reasonable defaults for axes
+def default_frequency_range(syslist):
+ """Compute a reasonable default frequency range for frequency
+ domain plots.
+
+ Usage
+ =====
+ omega = default_frequency_range(syslist)
+
+ Finds a reasonable default frequency range by examining the features
+ (poles and zeros) of the systems in syslist.
+
+ Parameters
+ ----------
+ syslist : linsys
+ List of linear input/output systems (single system is OK)
+
+ Return values
+ -------------
+ omega : freq_range
+ Range of frequencies in rad/sec
+ """
+ # This code looks at the poles and zeros of all of the systems that
+ # we are plotting and sets the frequency range to be one decade above
+ # and below the min and max feature frequencies, rounded to the nearest
+ # integer. It excludes poles and zeros at the origin. If no features
+ # are found, it turns logspace(-1, 1)
+
+ # Find the list of all poles and zeros in the systems
+ features = np.array(())
+ for sys in syslist:
+ # Add new features to the list
+ features = np.concatenate((features, np.abs(sys.poles)))
+ features = np.concatenate((features, np.abs(sys.zeros)))
+
+ # Get rid of poles and zeros at the origin
+ features = features[features != 0];
+
+ # Make sure there is at least one point in the range
+ if (features.shape[0] == 0): features = [1];
+
+ # Take the log of the features
+ features = np.log10(features)
+
+ # Set the range to be an order of magnitude beyond any features
+ omega = sp.logspace(np.floor(np.min(features))-1,
+ np.ceil(np.max(features))+1)
+
+ return omega
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2010-11-26 21:56:26 UTC (rev 32)
+++ trunk/src/matlab.py 2010-11-26 21:59:57 UTC (rev 33)
@@ -66,7 +66,7 @@
# Import MATLAB-like functions that can be used as-is
from ctrlutil import unwrap
-from freqplot import nyquist, gangof4
+from freqplot import nyquist, nichols, gangof4
from bdalg import series, parallel, negate, feedback
from pzmap import pzmap
from statefbk import ctrb, obsv, place, lqr
@@ -155,7 +155,7 @@
lti/bodemag - Bode magnitude diagram only
sigma - singular value frequency plot
* nyquist - Nyquist plot
- nichols - Nichols plot
+* nichols - Nichols plot
margin - gain and phase margins
lti/allmargin - all crossover frequencies and related gain/phase margins
lti/freqresp - frequency response over a frequency grid
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2010-11-26 21:56:32
|
Revision: 32
http://python-control.svn.sourceforge.net/python-control/?rev=32&view=rev
Author: murrayrm
Date: 2010-11-26 21:56:26 +0000 (Fri, 26 Nov 2010)
Log Message:
-----------
found minor bug in resetting axis limits for Bode plots
Modified Paths:
--------------
trunk/examples/pvtol-nested.py
Modified: trunk/examples/pvtol-nested.py
===================================================================
--- trunk/examples/pvtol-nested.py 2010-11-26 21:55:41 UTC (rev 31)
+++ trunk/examples/pvtol-nested.py 2010-11-26 21:56:26 UTC (rev 32)
@@ -82,6 +82,7 @@
#! Not yet implemented
# (gm, pm, wgc, wpc) = margin(L);
+#! TODO: this figure has something wrong; axis limits mismatch
figure(6); clf; subplot(221);
(magh, phaseh) = bode(L);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2010-11-26 21:55:47
|
Revision: 31
http://python-control.svn.sourceforge.net/python-control/?rev=31&view=rev
Author: murrayrm
Date: 2010-11-26 21:55:41 +0000 (Fri, 26 Nov 2010)
Log Message:
-----------
added header comment
Modified Paths:
--------------
trunk/external/yottalab.py
Modified: trunk/external/yottalab.py
===================================================================
--- trunk/external/yottalab.py 2010-11-06 16:26:19 UTC (rev 30)
+++ trunk/external/yottalab.py 2010-11-26 21:55:41 UTC (rev 31)
@@ -1,3 +1,5 @@
+# yottalab.py - Roberto Bucher's yottalab library
+
"""
This is a procedural interface to the yttalab library
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Ryan K. <rk...@si...> - 2010-11-06 17:06:55
|
I guess I don't feel bad then. I have one project (rst2beamer) where I do everything in git, but the other developers use svn. So, I have two different directories, one for the svn version and one for the git version. I use unison to sync the two folders if I make a change in git or they make one in svn. This is a fairly ugly hack, but I haven't seen a better, more reliable solution. There are supposed to be built in svn capabilities in git, but I can't seem to make them work. -- Ryan Krauss, Ph.D. Assistant Professor Mechanical Engineering Southern Illinois University Edwardsville On Sat, Nov 6, 2010 at 11:44 AM, Richard Murray <mu...@cd...> wrote: > I *wish* I had something to make this easy. Since your code is nicely contained in a single file, I just grabbed the contents from github and replaced the body of the original file. Works fine in this case, but would not work if there were dozens of files in a package. > > -richard > > On 6 Nov 2010, at 9:41 , Ryan Krauss wrote: > >> Just for my information, what are you using to make git and svn play >> nice with one another. You seem to have made efficient work of >> pulling in my code from git. I know a little git and very little svn. >> >> -- >> Ryan Krauss, Ph.D. >> Assistant Professor >> Mechanical Engineering >> Southern Illinois University Edwardsville >> >> >> >> On Sat, Nov 6, 2010 at 10:25 AM, Richard Murray <mu...@cd...> wrote: >>> Thanks, Ryan. I'll pull the latest version of your package into the sourceforge repository (so that developers can access it quickly) and update the python-control step response to use your code. >>> >>> -richard >>> >>> On 6 Nov 2010, at 8:22 , Ryan Krauss wrote: >>> >>>> FYI, I saw some mention of lsim2 in a recent email. My code includes >>>> a step response function that first tries lsim then uses lsim2 if >>>> there is an exception. The exceptions usually are caused by pure >>>> integrators. I think that defaulting to lsim2 will make the >>>> simulation slower. >>>> >>>> I have fixed a few bugs in my code. There may even be a bug fix in >>>> the digital conversion code that hasn't made it into my main branch. >>>> The most recent versions of my code are on github and on my course >>>> website for this semester: >>>> >>>> https://github.com/ryanGT/controls >>>> http://www.cs.siue.edu/~rkrauss/mechatronics/2010/python_controls/python.html >>>> >>>> -- >>>> Ryan Krauss, Ph.D. >>>> Assistant Professor >>>> Mechanical Engineering >>>> Southern Illinois University Edwardsville >>>> >>>> >>>> >>>> On Sat, Nov 6, 2010 at 8:45 AM, Richard Murray <mu...@cd...> wrote: >>>>> A couple of updates on status of the the python-control package, related to the commit earlier today: >>>>> >>>>> * Ryan Krauss has sent me his control.py package that he uses for his class at Southern Illinois University Edwardsville. He has implemented functions for continuous to discrete time conversion and root locus plus. Ryan's code is now in trunk/external/control.py and Ryan has graciously given us permission to incorporate functions into the python-control package. >>>>> >>>>> * Roberto Bucher has sent me his yottalab.py package that has a number of functions for state space computations in continuous and discrete time. This code is in trunk/external/yottalab.py. Roberto has also give us permission to incoporate functions into the python-control library. >>>>> >>>>> * I've added a pade() command developed by Sawyer Fuller at Caltech, using the description in Golub and Van Loan. This code is in trunk/src/delay.py. >>>>> >>>>> * A group of students at Princeton (Lauren Padilla, Kevin Chen and Steve Brunton) are working on a class project this term and will be helping to implement some additional functions in the library. I've encouraged them to use this list to send out thoughts and ideas about their plans. >>>>> >>>>> * In integrating some of the functions from Ryan and Roberto, I've made various changes to things and fixed a couple of bugs. See trunk/ChangeLog for details. >>>>> >>>>> That's it for now. Comments and feedback is always welcome. >>>>> >>>>> -richard >>>>> >>>>> >>>>> >>>>> ------------------------------------------------------------------------------ >>>>> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >>>>> David G. Thomson, author of the best-selling book "Blueprint to a >>>>> Billion" shares his insights and actions to help propel your >>>>> business during the next growth cycle. Listen Now! >>>>> http://p.sf.net/sfu/SAP-dev2dev >>>>> _______________________________________________ >>>>> python-control-discuss mailing list >>>>> pyt...@li... >>>>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >>>>> >>>> >>>> ------------------------------------------------------------------------------ >>>> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >>>> David G. Thomson, author of the best-selling book "Blueprint to a >>>> Billion" shares his insights and actions to help propel your >>>> business during the next growth cycle. Listen Now! >>>> http://p.sf.net/sfu/SAP-dev2dev >>>> _______________________________________________ >>>> python-control-discuss mailing list >>>> pyt...@li... >>>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >>> >>> >>> ------------------------------------------------------------------------------ >>> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >>> David G. Thomson, author of the best-selling book "Blueprint to a >>> Billion" shares his insights and actions to help propel your >>> business during the next growth cycle. Listen Now! >>> http://p.sf.net/sfu/SAP-dev2dev >>> _______________________________________________ >>> python-control-discuss mailing list >>> pyt...@li... >>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >>> >> >> ------------------------------------------------------------------------------ >> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >> David G. Thomson, author of the best-selling book "Blueprint to a >> Billion" shares his insights and actions to help propel your >> business during the next growth cycle. Listen Now! >> http://p.sf.net/sfu/SAP-dev2dev >> _______________________________________________ >> python-control-discuss mailing list >> pyt...@li... >> https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > > ------------------------------------------------------------------------------ > The Next 800 Companies to Lead America's Growth: New Video Whitepaper > David G. Thomson, author of the best-selling book "Blueprint to a > Billion" shares his insights and actions to help propel your > business during the next growth cycle. Listen Now! > http://p.sf.net/sfu/SAP-dev2dev > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > |
|
From: Richard M. <mu...@cd...> - 2010-11-06 16:45:01
|
I *wish* I had something to make this easy. Since your code is nicely contained in a single file, I just grabbed the contents from github and replaced the body of the original file. Works fine in this case, but would not work if there were dozens of files in a package. -richard On 6 Nov 2010, at 9:41 , Ryan Krauss wrote: > Just for my information, what are you using to make git and svn play > nice with one another. You seem to have made efficient work of > pulling in my code from git. I know a little git and very little svn. > > -- > Ryan Krauss, Ph.D. > Assistant Professor > Mechanical Engineering > Southern Illinois University Edwardsville > > > > On Sat, Nov 6, 2010 at 10:25 AM, Richard Murray <mu...@cd...> wrote: >> Thanks, Ryan. I'll pull the latest version of your package into the sourceforge repository (so that developers can access it quickly) and update the python-control step response to use your code. >> >> -richard >> >> On 6 Nov 2010, at 8:22 , Ryan Krauss wrote: >> >>> FYI, I saw some mention of lsim2 in a recent email. My code includes >>> a step response function that first tries lsim then uses lsim2 if >>> there is an exception. The exceptions usually are caused by pure >>> integrators. I think that defaulting to lsim2 will make the >>> simulation slower. >>> >>> I have fixed a few bugs in my code. There may even be a bug fix in >>> the digital conversion code that hasn't made it into my main branch. >>> The most recent versions of my code are on github and on my course >>> website for this semester: >>> >>> https://github.com/ryanGT/controls >>> http://www.cs.siue.edu/~rkrauss/mechatronics/2010/python_controls/python.html >>> >>> -- >>> Ryan Krauss, Ph.D. >>> Assistant Professor >>> Mechanical Engineering >>> Southern Illinois University Edwardsville >>> >>> >>> >>> On Sat, Nov 6, 2010 at 8:45 AM, Richard Murray <mu...@cd...> wrote: >>>> A couple of updates on status of the the python-control package, related to the commit earlier today: >>>> >>>> * Ryan Krauss has sent me his control.py package that he uses for his class at Southern Illinois University Edwardsville. He has implemented functions for continuous to discrete time conversion and root locus plus. Ryan's code is now in trunk/external/control.py and Ryan has graciously given us permission to incorporate functions into the python-control package. >>>> >>>> * Roberto Bucher has sent me his yottalab.py package that has a number of functions for state space computations in continuous and discrete time. This code is in trunk/external/yottalab.py. Roberto has also give us permission to incoporate functions into the python-control library. >>>> >>>> * I've added a pade() command developed by Sawyer Fuller at Caltech, using the description in Golub and Van Loan. This code is in trunk/src/delay.py. >>>> >>>> * A group of students at Princeton (Lauren Padilla, Kevin Chen and Steve Brunton) are working on a class project this term and will be helping to implement some additional functions in the library. I've encouraged them to use this list to send out thoughts and ideas about their plans. >>>> >>>> * In integrating some of the functions from Ryan and Roberto, I've made various changes to things and fixed a couple of bugs. See trunk/ChangeLog for details. >>>> >>>> That's it for now. Comments and feedback is always welcome. >>>> >>>> -richard >>>> >>>> >>>> >>>> ------------------------------------------------------------------------------ >>>> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >>>> David G. Thomson, author of the best-selling book "Blueprint to a >>>> Billion" shares his insights and actions to help propel your >>>> business during the next growth cycle. Listen Now! >>>> http://p.sf.net/sfu/SAP-dev2dev >>>> _______________________________________________ >>>> python-control-discuss mailing list >>>> pyt...@li... >>>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >>>> >>> >>> ------------------------------------------------------------------------------ >>> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >>> David G. Thomson, author of the best-selling book "Blueprint to a >>> Billion" shares his insights and actions to help propel your >>> business during the next growth cycle. Listen Now! >>> http://p.sf.net/sfu/SAP-dev2dev >>> _______________________________________________ >>> python-control-discuss mailing list >>> pyt...@li... >>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >> >> >> ------------------------------------------------------------------------------ >> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >> David G. Thomson, author of the best-selling book "Blueprint to a >> Billion" shares his insights and actions to help propel your >> business during the next growth cycle. Listen Now! >> http://p.sf.net/sfu/SAP-dev2dev >> _______________________________________________ >> python-control-discuss mailing list >> pyt...@li... >> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >> > > ------------------------------------------------------------------------------ > The Next 800 Companies to Lead America's Growth: New Video Whitepaper > David G. Thomson, author of the best-selling book "Blueprint to a > Billion" shares his insights and actions to help propel your > business during the next growth cycle. Listen Now! > http://p.sf.net/sfu/SAP-dev2dev > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Ryan K. <rk...@si...> - 2010-11-06 16:42:04
|
Just for my information, what are you using to make git and svn play nice with one another. You seem to have made efficient work of pulling in my code from git. I know a little git and very little svn. -- Ryan Krauss, Ph.D. Assistant Professor Mechanical Engineering Southern Illinois University Edwardsville On Sat, Nov 6, 2010 at 10:25 AM, Richard Murray <mu...@cd...> wrote: > Thanks, Ryan. I'll pull the latest version of your package into the sourceforge repository (so that developers can access it quickly) and update the python-control step response to use your code. > > -richard > > On 6 Nov 2010, at 8:22 , Ryan Krauss wrote: > >> FYI, I saw some mention of lsim2 in a recent email. My code includes >> a step response function that first tries lsim then uses lsim2 if >> there is an exception. The exceptions usually are caused by pure >> integrators. I think that defaulting to lsim2 will make the >> simulation slower. >> >> I have fixed a few bugs in my code. There may even be a bug fix in >> the digital conversion code that hasn't made it into my main branch. >> The most recent versions of my code are on github and on my course >> website for this semester: >> >> https://github.com/ryanGT/controls >> http://www.cs.siue.edu/~rkrauss/mechatronics/2010/python_controls/python.html >> >> -- >> Ryan Krauss, Ph.D. >> Assistant Professor >> Mechanical Engineering >> Southern Illinois University Edwardsville >> >> >> >> On Sat, Nov 6, 2010 at 8:45 AM, Richard Murray <mu...@cd...> wrote: >>> A couple of updates on status of the the python-control package, related to the commit earlier today: >>> >>> * Ryan Krauss has sent me his control.py package that he uses for his class at Southern Illinois University Edwardsville. He has implemented functions for continuous to discrete time conversion and root locus plus. Ryan's code is now in trunk/external/control.py and Ryan has graciously given us permission to incorporate functions into the python-control package. >>> >>> * Roberto Bucher has sent me his yottalab.py package that has a number of functions for state space computations in continuous and discrete time. This code is in trunk/external/yottalab.py. Roberto has also give us permission to incoporate functions into the python-control library. >>> >>> * I've added a pade() command developed by Sawyer Fuller at Caltech, using the description in Golub and Van Loan. This code is in trunk/src/delay.py. >>> >>> * A group of students at Princeton (Lauren Padilla, Kevin Chen and Steve Brunton) are working on a class project this term and will be helping to implement some additional functions in the library. I've encouraged them to use this list to send out thoughts and ideas about their plans. >>> >>> * In integrating some of the functions from Ryan and Roberto, I've made various changes to things and fixed a couple of bugs. See trunk/ChangeLog for details. >>> >>> That's it for now. Comments and feedback is always welcome. >>> >>> -richard >>> >>> >>> >>> ------------------------------------------------------------------------------ >>> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >>> David G. Thomson, author of the best-selling book "Blueprint to a >>> Billion" shares his insights and actions to help propel your >>> business during the next growth cycle. Listen Now! >>> http://p.sf.net/sfu/SAP-dev2dev >>> _______________________________________________ >>> python-control-discuss mailing list >>> pyt...@li... >>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >>> >> >> ------------------------------------------------------------------------------ >> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >> David G. Thomson, author of the best-selling book "Blueprint to a >> Billion" shares his insights and actions to help propel your >> business during the next growth cycle. Listen Now! >> http://p.sf.net/sfu/SAP-dev2dev >> _______________________________________________ >> python-control-discuss mailing list >> pyt...@li... >> https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > > ------------------------------------------------------------------------------ > The Next 800 Companies to Lead America's Growth: New Video Whitepaper > David G. Thomson, author of the best-selling book "Blueprint to a > Billion" shares his insights and actions to help propel your > business during the next growth cycle. Listen Now! > http://p.sf.net/sfu/SAP-dev2dev > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > |
|
From: <mur...@us...> - 2010-11-06 16:26:26
|
Revision: 30
http://python-control.svn.sourceforge.net/python-control/?rev=30&view=rev
Author: murrayrm
Date: 2010-11-06 16:26:19 +0000 (Sat, 06 Nov 2010)
Log Message:
-----------
Latest version of Ryan Krauss's control package. Significant changes
(from git commit log):
2010-10-26
added simplify method to cancel poles and zeros at the same location
2010-10-08
Commented out one line for mpl compatability (new version of mpl)
2010-09-10
Added a simplify method to the TF class to cancel poles and zeros that are with a certain tolerance of one another
2010-09-01
Bug fix update for scipy 0.8 and its changes to scipy.io
2010-06-07
version 1.1.3 - I don't remember the changes (I am a bad git user).
2010-04-19
updated to version 1.1.2 mainly getting rid of annoying error message
2010-03-31
Merge branch 'master' of github.com:ryanGT/controls
- debugged leading zeros in numerators adding to Carten's patch
- debugged TF.lsim and lsim2; set-up build_controls_py.py to make exe
2009-09-09
Small bug fixes
2009-07-04
Add ButterWorth filter Compensator
Modified Paths:
--------------
trunk/external/controls.py
Modified: trunk/external/controls.py
===================================================================
--- trunk/external/controls.py 2010-11-06 13:03:55 UTC (rev 29)
+++ trunk/external/controls.py 2010-11-06 16:26:19 UTC (rev 30)
@@ -11,10 +11,13 @@
from math import atan2, log10
from scipy import *
+from scipy import signal
+from scipy import interpolate, integrate
from scipy.linalg import inv as inverse
from scipy.optimize import newton, fmin, fminbound
-from scipy.io import read_array, save, loadmat, write_array
+#from scipy.io import read_array, save, loadmat, write_array
from scipy import signal
+from numpy.linalg import LinAlgError
from IPython.Debugger import Pdb
@@ -22,7 +25,7 @@
from matplotlib.ticker import LogFormatterMathtext
-version = '1.1.0'
+version = '1.1.4'
class MyFormatter(LogFormatterMathtext):
def __call__(self, x, pos=None):
@@ -126,7 +129,10 @@
def prependzeros(num, den):
nd = len(den)
- nn = len(num)
+ if isscalar(num):
+ nn = 1
+ else:
+ nn = len(num)
if nn < nd:
zvect = zeros(nd-nn)
numout = r_[zvect, num]
@@ -388,26 +394,42 @@
def __setattr__(self, attr, val):
realizable = False
if hasattr(self, 'den') and hasattr(self, 'num'):
- realizable = _realizable(self.num, self.den)
+ realizable = _realizable(self.num, self.den)
if realizable:
- signal.lti.__setattr__(self, attr, val)
+ signal.lti.__setattr__(self, attr, val)
else:
- self.__dict__[attr] = val
+ self.__dict__[attr] = val
- def __init__(self, num, den, dt=0.01, maxt=5.0, myvar='s'):
+ def __init__(self, num, den, dt=0.01, maxt=5.0, myvar='s', label='G'):
"""num and den are either scalar constants or lists that are
passed to scipy.poly1d to create a list of coefficients."""
#print('in TransferFunction.__init__, dt=%s' % dt)
if _realizable(num, den):
- signal.lti.__init__(self, num, den)
+ num = atleast_1d(num)
+ den = atleast_1d(den)
+ start_num_ind = nonzero(num)[0][0]
+ start_den_ind = nonzero(den)[0][0]
+ num_ = num[start_num_ind:]
+ den_ = den[start_den_ind:]
+ signal.lti.__init__(self, num_, den_)
+ else:
+ z, p, k = signal.tf2zpk(num, den)
+ self.gain = k
self.num = poly1d(num)
self.den = poly1d(den)
self.dt = dt
self.myvar = myvar
self.maxt = maxt
+ self.label = label
+ def print_poles(self, label=None):
+ if label is None:
+ label = self.label
+ print(label +' poles =' + str(self.poles))
+
+
def __repr__(self, labelstr='controls.TransferFunction'):
nstr=str(self.num)#.strip()
dstr=str(self.den)#.strip()
@@ -502,6 +524,16 @@
return self.dt
+ def simplify(self, rtol=1e-5, atol=1e-10):
+ """Return a new TransferFunction object with poles and zeros
+ that nearly cancel (within real or absolutie tolerance rtol
+ and atol) removed."""
+ gain = self.gain
+ new_num, new_den = polyfactor(self.num, self.den, prepend=False)
+ newtf = self.__class__(new_num*gain, new_den)
+ return newtf
+
+
def ToLatex(self, eps=1e-12, fmt='%0.5g', ds=True):
mynum = self.num
myden = self.den
@@ -590,7 +622,7 @@
def rlocfind(self, poleloc):
self.poleloc = poleloc
kinit,pinit = _k_poles(self,poleloc)
- k = optimize.fmin(self.opt,[kinit])[0]
+ k = fmin(self.opt,[kinit])[0]
poles = self._RLFindRoots([k])
poles = self._RLSortRoots(poles)
return k, poles
@@ -646,26 +678,102 @@
self.num = self.num/const
self.den = self.den/const
- def lsim(self, u, t, interp=0, returnall=False, X0=None):
+ def lsim(self, u, t, interp=0, returnall=False, X0=None, hmax=None):
"""Find the response of the TransferFunction to the input u
with time vector t. Uses signal.lsim.
return y the response of the system."""
+ try:
+ out = signal.lsim(self, u, t, interp=interp, X0=X0)
+ except LinAlgError:
+ #if the system has a pure integrator, lsim won't work.
+ #Call lsim2.
+ out = self.lsim2(u, t, X0=X0, returnall=True, hmax=hmax)
+ #override returnall because it is handled below
if returnall:#most users will just want the system output y,
#but some will need the (t, y, x) tuple that
#signal.lsim returns
- return signal.lsim(self, u, t, interp=interp, X0=X0)
+ return out
else:
- return signal.lsim(self, u, t, interp=interp, X0=X0)[1]
+ return out[1]
- def lsim2(self, u, t, returnall=False, X0=None):
- #tempsys=signal.lti(self.num,self.den)
- if returnall:
- return signal.lsim2(self, u, t, X0=X0)
- else:
- return signal.lsim2(self, u, t, X0=X0)[1]
+## def lsim2(self, u, t, returnall=False, X0=None):
+## #tempsys=signal.lti(self.num,self.den)
+## if returnall:
+## return signal.lsim2(self, u, t, X0=X0)
+## else:
+## return signal.lsim2(self, u, t, X0=X0)[1]
+ def lsim2(self, U, T, X0=None, returnall=False, hmax=None):
+ """Simulate output of a continuous-time linear system, using ODE solver.
+ Inputs:
+
+ system -- an instance of the LTI class or a tuple describing the
+ system. The following gives the number of elements in
+ the tuple and the interpretation.
+ 2 (num, den)
+ 3 (zeros, poles, gain)
+ 4 (A, B, C, D)
+ U -- an input array describing the input at each time T
+ (linear interpolation is assumed between given times).
+ If there are multiple inputs, then each column of the
+ rank-2 array represents an input.
+ T -- the time steps at which the input is defined and at which
+ the output is desired.
+ X0 -- (optional, default=0) the initial conditions on the state vector.
+
+ Outputs: (T, yout, xout)
+
+ T -- the time values for the output.
+ yout -- the response of the system.
+ xout -- the time-evolution of the state-vector.
+ """
+ # system is an lti system or a sequence
+ # with 2 (num, den)
+ # 3 (zeros, poles, gain)
+ # 4 (A, B, C, D)
+ # describing the system
+ # U is an input vector at times T
+ # if system describes multiple outputs
+ # then U can be a rank-2 array with the number of columns
+ # being the number of inputs
+
+ # rather than use lsim, use direct integration and matrix-exponential.
+ if hmax is None:
+ hmax = T[1]-T[0]
+ U = atleast_1d(U)
+ T = atleast_1d(T)
+ if len(U.shape) == 1:
+ U = U.reshape((U.shape[0],1))
+ sU = U.shape
+ if len(T.shape) != 1:
+ raise ValueError, "T must be a rank-1 array."
+ if sU[0] != len(T):
+ raise ValueError, "U must have the same number of rows as elements in T."
+ if sU[1] != self.inputs:
+ raise ValueError, "System does not define that many inputs."
+
+ if X0 is None:
+ X0 = zeros(self.B.shape[0],self.A.dtype)
+
+ # for each output point directly integrate assume zero-order hold
+ # or linear interpolation.
+
+ ufunc = interpolate.interp1d(T, U, kind='linear', axis=0, \
+ bounds_error=False)
+
+ def fprime(x, t, self, ufunc):
+ return dot(self.A,x) + squeeze(dot(self.B,nan_to_num(ufunc([t]))))
+
+ xout = integrate.odeint(fprime, X0, T, args=(self, ufunc), hmax=hmax)
+ yout = dot(self.C,transpose(xout)) + dot(self.D,transpose(U))
+ if returnall:
+ return T, squeeze(transpose(yout)), xout
+ else:
+ return squeeze(transpose(yout))
+
+
def residue(self, tol=1e-3, verbose=0):
"""from scipy.signal.residue:
@@ -748,7 +856,8 @@
def FreqResp(self, f, fignum=1, fig=None, clear=True, \
- grid=True, legend=None, legloc=1, legsub=1, **kwargs):
+ grid=True, legend=None, legloc=1, legsub=1, \
+ use_rad=False, **kwargs):
"""Compute the frequency response of the transfer function
using the frequency vector f, returning a complex vector.
@@ -765,10 +874,14 @@
if testvect.all():
s=f#then you really sent me s and not f
else:
- s=2.0j*pi*f
+ if use_rad:
+ s = 1.0j*f
+ else:
+ s=2.0j*pi*f
self.comp = self.num(s)/self.den(s)
self.dBmag = 20*log10(abs(self.comp))
- self.phase = angle(self.comp,1)
+ rphase = unwrap(angle(self.comp))
+ self.phase = rphase*180.0/pi
if fig is None:
if fignum is not None:
@@ -778,19 +891,19 @@
if fig is not None:
if clear:
fig.clf()
+ ax1 = fig.add_subplot(2,1,1)
+ ax2 = fig.add_subplot(2,1,2, sharex=ax1)
+ else:
+ ax1 = fig.axes[0]
+ ax2 = fig.axes[1]
if fig is not None:
- myargs=['linetype','colors','linewidth']
+ myargs=['linetype','linewidth']
subkwargs={}
for key in myargs:
if kwargs.has_key(key):
subkwargs[key]=kwargs[key]
- if clear:
- fig.clf()
- ax1 = fig.add_subplot(2,1,1)
- #if clear:
- # ax1.cla()
- myind=ax1._get_lines.count
+ #myind=ax1._get_lines.count
mylines=_PlotMag(f, self, axis=ax1, **subkwargs)
ax1.set_ylabel('Mag. Ratio (dB)')
ax1.xaxis.set_major_formatter(MyFormatter())
@@ -798,12 +911,12 @@
ax1.grid(1)
if legend is not None and legsub==1:
ax1.legend(legend, legloc)
- ax2 = fig.add_subplot(2,1,2, sharex=ax1)
- #if clear:
- # ax2.cla()
mylines=_PlotPhase(f, self, axis=ax2, **subkwargs)
ax2.set_ylabel('Phase (deg.)')
- ax2.set_xlabel('Freq. (Hz)')
+ if use_rad:
+ ax2.set_xlabel('$\\omega$ (rad./sec.)')
+ else:
+ ax2.set_xlabel('Freq. (Hz)')
ax2.xaxis.set_major_formatter(MyFormatter())
if grid:
ax2.grid(1)
@@ -1142,8 +1255,8 @@
yvect = shift(yvect, curY)
return Ydig*dv
+TF = TransferFunction
-
class Input(TransferFunction):
def __repr__(self):
return TransferFunction.__repr__(self, labelstr='controls.Input')
@@ -1259,8 +1372,17 @@
else:
return vin
+class ButterworthFilter(Compensator):
+ def __init__(self,fc,mag=1.0):
+ """Create a compensator that is a second order Butterworth
+ filter. fc is the corner frequency in Hz and mag is the low
+ frequency magnitude so that the transfer function will be
+ mag*wn**2/(s**2+2*z*wn*s+wn**2) where z=1/sqrt(2) and
+ wn=2.0*pi*fc."""
+ z=1.0/sqrt(2.0)
+ wn=2.0*pi*fc
+ Compensator.__init__(self,mag*wn**2,[1.0,2.0*z*wn,wn**2])
-
class Closed_Loop_System_with_Sat(object):
def __init__(self, plant_tf, Kp, sat):
self.plant_tf = plant_tf
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
7 messages has been excluded from this view by a project administrator.
