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:15:13
|
Revision: 66
http://python-control.svn.sourceforge.net/python-control/?rev=66&view=rev
Author: kkchen
Date: 2011-02-08 22:15:06 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Old TransferFunction class removed from xferfcn.py. The "developer version"
xTransferFunction is now renamed TransferFunction.
Related fixes in bdalg.py, TestXferFcn.py, and TestBDAlg.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestBDAlg.py
branches/control-0.4a/src/TestXferFcn.py
branches/control-0.4a/src/bdalg.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestBDAlg.py
===================================================================
--- branches/control-0.4a/src/TestBDAlg.py 2011-02-08 22:15:01 UTC (rev 65)
+++ branches/control-0.4a/src/TestBDAlg.py 2011-02-08 22:15:06 UTC (rev 66)
@@ -1,7 +1,7 @@
#!/usr/bin/env python
import numpy as np
-from xferfcn import xTransferFunction
+from xferfcn import TransferFunction
from statesp import StateSpace
from bdalg import feedback
import unittest
@@ -15,7 +15,7 @@
"""This contains some random LTI systems and scalars for testing."""
# Two random SISO systems.
- self.sys1 = xTransferFunction([1, 2], [1, 2, 3])
+ self.sys1 = TransferFunction([1, 2], [1, 2, 3])
self.sys2 = StateSpace([[1., 4.], [3., 2.]], [[1.], [-4.]],
[[1., 0.]], [[0.]])
# Two random scalars.
Modified: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:15:01 UTC (rev 65)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:15:06 UTC (rev 66)
@@ -1,7 +1,7 @@
#!/usr/bin/env python
import numpy as np
-from xferfcn import xTransferFunction
+from xferfcn import TransferFunction
import unittest
class TestXferFcn(unittest.TestCase):
@@ -15,42 +15,42 @@
def testBadInputType(self):
"""Give the constructor invalid input types."""
- self.assertRaises(TypeError, xTransferFunction, [[0., 1.], [2., 3.]],
+ 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, xTransferFunction, [[[1.]]],
+ self.assertRaises(ValueError, TransferFunction, [[[1.]]],
[[[1.], [2., 3.]]])
- self.assertRaises(ValueError, xTransferFunction, [[[1.]]],
+ self.assertRaises(ValueError, TransferFunction, [[[1.]]],
[[[1.]], [[2., 3.]]])
- self.assertRaises(ValueError, xTransferFunction, [[[1.]]],
+ 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, xTransferFunction, 1.,
+ self.assertRaises(ValueError, TransferFunction, 1.,
[[[1.]], [[2.], [3.]]])
- self.assertRaises(ValueError, xTransferFunction, [[[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, xTransferFunction, 1., 0.)
- self.assertRaises(ValueError, xTransferFunction,
+ 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 = xTransferFunction([[[1., 2.]]], [[[4., 5.]]])
- sys2 = xTransferFunction([[[4., 3.]], [[1., 2.]]],
+ 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)
@@ -60,21 +60,21 @@
def testMulInconsistentDimension(self):
"""Multiply two transfer function matrices of incompatible sizes."""
- sys1 = xTransferFunction([[[1., 2.], [4., 5.]], [[2., 5.], [4., 3.]]],
+ sys1 = TransferFunction([[[1., 2.], [4., 5.]], [[2., 5.], [4., 3.]]],
[[[6., 2.], [4., 1.]], [[6., 7.], [2., 4.]]])
- sys2 = xTransferFunction([[[1.]], [[2.]], [[3.]]],
+ 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 xTransferFunction._truncatecoeff
+ # Tests for TransferFunction._truncatecoeff
def testTruncateCoeff1(self):
"""Remove extraneous zeros in polynomial representations."""
- sys1 = xTransferFunction([0., 0., 1., 2.], [[[0., 0., 0., 3., 2., 1.]]])
+ 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.]]])
@@ -82,17 +82,17 @@
def testTruncateCoeff2(self):
"""Remove extraneous zeros in polynomial representations."""
- sys1 = xTransferFunction([0., 0., 0.], 1.)
+ sys1 = TransferFunction([0., 0., 0.], 1.)
np.testing.assert_array_equal(sys1.num, [[[0.]]])
np.testing.assert_array_equal(sys1.den, [[[1.]]])
- # Tests for xTransferFunction.__neg__
+ # Tests for TransferFunction.__neg__
def testNegScalar(self):
"""Negate a direct feedthrough system."""
- sys1 = xTransferFunction(2., np.array([-3]))
+ sys1 = TransferFunction(2., np.array([-3]))
sys2 = - sys1
np.testing.assert_array_equal(sys2.num, [[[-2.]]])
@@ -101,7 +101,7 @@
def testNegSISO(self):
"""Negate a SISO system."""
- sys1 = xTransferFunction([1., 3., 5], [1., 6., 2., -1.])
+ sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1.])
sys2 = - sys1
np.testing.assert_array_equal(sys2.num, [[[-1., -3., -5.]]])
@@ -117,22 +117,22 @@
den1 = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]],
[[3., 0., .0], [2., -1., -1.], [1.]]]
- sys1 = xTransferFunction(num1, den1)
+ sys1 = TransferFunction(num1, den1)
sys2 = - sys1
- sys3 = xTransferFunction(num3, den1)
+ 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 xTransferFunction.__add__
+ # Tests for TransferFunction.__add__
def testAddScalar(self):
"""Add two direct feedthrough systems."""
- sys1 = xTransferFunction(1., [[[1.]]])
- sys2 = xTransferFunction(np.array([2.]), [1.])
+ sys1 = TransferFunction(1., [[[1.]]])
+ sys2 = TransferFunction(np.array([2.]), [1.])
sys3 = sys1 + sys2
np.testing.assert_array_equal(sys3.num, 3.)
@@ -141,8 +141,8 @@
def testAddSISO(self):
"""Add two SISO systems."""
- sys1 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
- sys2 = xTransferFunction([[np.array([-1., 3.])]], [[[1., 0., -1.]]])
+ 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!
@@ -165,8 +165,8 @@
den3 = [[[3., -2., -4.], [1., 2., 3., 0., 0.], [-2., -1., 1.]],
[[-12., -9., 6., 0., 0.], [2., -1., -1.], [1., 0.]]]
- sys1 = xTransferFunction(num1, den1)
- sys2 = xTransferFunction(num2, den2)
+ sys1 = TransferFunction(num1, den1)
+ sys2 = TransferFunction(num2, den2)
sys3 = sys1 + sys2
for i in range(sys3.outputs):
@@ -174,13 +174,13 @@
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 xTransferFunction.__sub__
+ # Tests for TransferFunction.__sub__
def testSubScalar(self):
"""Add two direct feedthrough systems."""
- sys1 = xTransferFunction(1., [[[1.]]])
- sys2 = xTransferFunction(np.array([2.]), [1.])
+ sys1 = TransferFunction(1., [[[1.]]])
+ sys2 = TransferFunction(np.array([2.]), [1.])
sys3 = sys1 - sys2
np.testing.assert_array_equal(sys3.num, -1.)
@@ -189,8 +189,8 @@
def testSubSISO(self):
"""Add two SISO systems."""
- sys1 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
- sys2 = xTransferFunction([[np.array([-1., 3.])]], [[[1., 0., -1.]]])
+ sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys2 = TransferFunction([[np.array([-1., 3.])]], [[[1., 0., -1.]]])
sys3 = sys1 - sys2
sys4 = sys2 - sys1
@@ -215,8 +215,8 @@
den3 = [[[3., -2., -4], [1., 2., 3., 0., 0.], [1]],
[[-12., -9., 6., 0., 0.], [2., -1., -1], [1., 0.]]]
- sys1 = xTransferFunction(num1, den1)
- sys2 = xTransferFunction(num2, den2)
+ sys1 = TransferFunction(num1, den1)
+ sys2 = TransferFunction(num2, den2)
sys3 = sys1 - sys2
for i in range(sys3.outputs):
@@ -224,13 +224,13 @@
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 xTransferFunction.__mul__
+ # Tests for TransferFunction.__mul__
def testMulScalar(self):
"""Multiply two direct feedthrough systems."""
- sys1 = xTransferFunction(2., [1.])
- sys2 = xTransferFunction(1., 4.)
+ sys1 = TransferFunction(2., [1.])
+ sys2 = TransferFunction(1., 4.)
sys3 = sys1 * sys2
sys4 = sys1 * sys2
@@ -242,8 +242,8 @@
def testMulSISO(self):
"""Multiply two SISO systems."""
- sys1 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
- sys2 = xTransferFunction([[[-1., 3.]]], [[[1., 0., -1.]]])
+ sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys2 = TransferFunction([[[-1., 3.]]], [[[1., 0., -1.]]])
sys3 = sys1 * sys2
sys4 = sys2 * sys1
@@ -273,8 +273,8 @@
0., 0.]],
[[-48., 60., 84., -81., -45., 21., 9., 0., 0., 0., 0., 0., 0.]]]
- sys1 = xTransferFunction(num1, den1)
- sys2 = xTransferFunction(num2, den2)
+ sys1 = TransferFunction(num1, den1)
+ sys2 = TransferFunction(num2, den2)
sys3 = sys1 * sys2
for i in range(sys3.outputs):
@@ -282,13 +282,13 @@
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 xTransferFunction.__div__
+ # Tests for TransferFunction.__div__
def testDivScalar(self):
"""Divide two direct feedthrough systems."""
- sys1 = xTransferFunction(np.array([3.]), -4.)
- sys2 = xTransferFunction(5., 2.)
+ sys1 = TransferFunction(np.array([3.]), -4.)
+ sys2 = TransferFunction(5., 2.)
sys3 = sys1 / sys2
np.testing.assert_array_equal(sys3.num, [[[6.]]])
@@ -297,8 +297,8 @@
def testDivSISO(self):
"""Divide two SISO systems."""
- sys1 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
- sys2 = xTransferFunction([[[-1., 3.]]], [[[1., 0., -1.]]])
+ sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys2 = TransferFunction([[[-1., 3.]]], [[[1., 0., -1.]]])
sys3 = sys1 / sys2
sys4 = sys2 / sys1
@@ -307,12 +307,12 @@
np.testing.assert_array_equal(sys4.num, sys3.den)
np.testing.assert_array_equal(sys4.den, sys3.num)
- # Tests for xTransferFunction.evalfr.
+ # Tests for TransferFunction.evalfr.
def testEvalFrSISO(self):
"""Evaluate the frequency response of a SISO system at one frequency."""
- sys = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
np.testing.assert_array_almost_equal(sys.evalfr(1.),
np.array([[-0.5 - 0.5j]]))
@@ -326,20 +326,20 @@
[[1.], [4., 0.], [1., -4., 3.]]]
den = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]],
[[3., 0., .0], [2., -1., -1.], [1.]]]
- sys = xTransferFunction(num, den)
+ sys = TransferFunction(num, den)
resp = [[0.147058823529412 + 0.0882352941176471j, -0.75, 1.],
[-0.083333333333333, -0.188235294117647 - 0.847058823529412j,
-1. - 8.j]]
np.testing.assert_array_almost_equal(sys.evalfr(2.), resp)
- # Tests for xTransferFunction.freqresp.
+ # Tests for TransferFunction.freqresp.
def testFreqRespSISO(self):
"""Evaluate the magnitude and phase of a SISO system at multiple
frequencies."""
- sys = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
truemag = [[[4.63507337473906, 0.707106781186548, 0.0866592803995351]]]
truephase = [[[-2.89596891081488, -2.35619449019234,
@@ -360,7 +360,7 @@
[[1.], [4., 0.], [1., -4., 3.]]]
den = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]],
[[3., 0., .0], [2., -1., -1.], [1.]]]
- sys = xTransferFunction(num, den)
+ sys = TransferFunction(num, den)
trueomega = [0.1, 1., 10.]
truemag = [[[0.496287094505259, 0.307147558416976, 0.0334738176210382],
@@ -380,12 +380,12 @@
np.testing.assert_array_almost_equal(phase, truephase)
np.testing.assert_array_equal(omega, trueomega)
- # Tests for xTransferFunction.feedback.
+ # Tests for TransferFunction.feedback.
def testFeedbackSISO(self):
- sys1 = xTransferFunction([-1., 4.], [1., 3., 5.])
- sys2 = xTransferFunction([2., 3., 0.], [1., -3., 4., 0])
+ sys1 = TransferFunction([-1., 4.], [1., 3., 5.])
+ sys2 = TransferFunction([2., 3., 0.], [1., -3., 4., 0])
sys3 = sys1.feedback(sys2)
sys4 = sys1.feedback(sys2, 1)
Modified: branches/control-0.4a/src/bdalg.py
===================================================================
--- branches/control-0.4a/src/bdalg.py 2011-02-08 22:15:01 UTC (rev 65)
+++ branches/control-0.4a/src/bdalg.py 2011-02-08 22:15:06 UTC (rev 66)
@@ -47,6 +47,21 @@
import xferfcn as tf
import statesp as ss
+def series(sys1, sys2):
+ """Return the series connection sys1 * sys2 for --> sys2 --> sys1 -->."""
+
+ return sys1 * sys2
+
+def parallel(sys1, sys2):
+ """Return the parallel connection sys1 + sys2."""
+
+ return sys1 + sys2
+
+def negate(sys):
+ """Return the negative of a system."""
+
+ return -sys;
+
def feedback(sys1, sys2, sign=-1):
"""Feedback interconnection between two I/O systems.
@@ -79,11 +94,11 @@
xferfcn.feedback."""
# Check for correct input types.
- if not isinstance(sys1, (int, long, float, complex, tf.xTransferFunction,
+ if not isinstance(sys1, (int, long, float, complex, tf.TransferFunction,
ss.StateSpace)):
raise TypeError("sys1 must be a TransferFunction or StateSpace object, \
or a scalar.")
- if not isinstance(sys2, (int, long, float, complex, tf.xTransferFunction,
+ if not isinstance(sys2, (int, long, float, complex, tf.TransferFunction,
ss.StateSpace)):
raise TypeError("sys2 must be a TransferFunction or StateSpace object, \
or a scalar.")
@@ -91,7 +106,7 @@
# If sys1 is a scalar, convert it to the appropriate LTI type so that we can
# its feedback member function.
if isinstance(sys1, (int, long, float, complex)):
- if isinstance(sys2, tf.xTransferFunction):
+ if isinstance(sys2, tf.TransferFunction):
sys1 = tf.convertToTransferFunction(sys1)
elif isinstance(sys2, ss.StateSpace):
sys1 = ss.convertToStateSpace(sys1)
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:15:01 UTC (rev 65)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:15:06 UTC (rev 66)
@@ -54,7 +54,7 @@
import statesp
from lti2 import Lti2
-class xTransferFunction(Lti2):
+class TransferFunction(Lti2):
"""The TransferFunction class is derived from the Lti2 parent class. The
main data members are 'num' and 'den', which are 2-D lists of arrays
containing MIMO numerator and denominator coefficients. For example,
@@ -210,13 +210,13 @@
for j in range(self.inputs):
num[i][j] *= -1
- return xTransferFunction(num, self.den)
+ return TransferFunction(num, self.den)
def __add__(self, other):
"""Add two transfer functions (parallel connection)."""
# Convert the second argument to a transfer function.
- if not isinstance(other, xTransferFunction):
+ if not isinstance(other, TransferFunction):
other = convertToTransferFunction(other, self.inputs, self.outputs)
# Check that the input-output sizes are consistent.
@@ -236,7 +236,7 @@
num[i][j], den[i][j] = _addSISO(self.num[i][j], self.den[i][j],
other.num[i][j], other.den[i][j])
- return xTransferFunction(num, den)
+ return TransferFunction(num, den)
def __radd__(self, other):
"""Add two transfer functions (parallel connection)."""
@@ -257,7 +257,7 @@
"""Multiply two transfer functions (serial connection)."""
# Convert the second argument to a transfer function.
- if not isinstance(other, xTransferFunction):
+ if not isinstance(other, TransferFunction):
other = convertToTransferFunction(other, self.inputs, self.inputs)
# Check that the input-output sizes are consistent.
@@ -285,7 +285,7 @@
num[i][j], den[i][j] = _addSISO(num[i][j], den[i][j],
num_summand[k], den_summand[k])
- return xTransferFunction(num, den)
+ return TransferFunction(num, den)
def __rmul__(self, other):
"""Multiply two transfer functions (serial connection)."""
@@ -298,17 +298,17 @@
if self.inputs > 1 or self.outputs > 1 or \
other.inputs > 1 or other.outputs > 1:
- raise NotImplementedError("xTransferFunction.__div__ is currently \
+ raise NotImplementedError("TransferFunction.__div__ is currently \
implemented only for SISO systems.")
# Convert the second argument to a transfer function.
- if not isinstance(other, xTransferFunction):
+ if not isinstance(other, TransferFunction):
other = convertToTransferFunction(other, 1, 1)
num = sp.polymul(self.num[0][0], other.den[0][0])
den = sp.polymul(self.den[0][0], other.num[0][0])
- return xTransferFunction(num, den)
+ return TransferFunction(num, den)
# TODO: Division of MIMO transfer function objects is quite difficult.
def __rdiv__(self, other):
@@ -316,7 +316,7 @@
if self.inputs > 1 or self.outputs > 1 or \
other.inputs > 1 or other.outputs > 1:
- raise NotImplementedError("xTransferFunction.__rdiv__ is currently \
+ raise NotImplementedError("TransferFunction.__rdiv__ is currently \
implemented only for SISO systems.")
return other / self
@@ -371,7 +371,7 @@
if self.inputs > 1 or self.outputs > 1 or \
other.inputs > 1 or other.outputs > 1:
- raise NotImplementedError("xTransferFunction.feedback is currently \
+ raise NotImplementedError("TransferFunction.feedback is currently \
only implemented for SISO functions.")
num1 = self.num[0][0]
@@ -382,185 +382,13 @@
num = sp.polymul(num1, den2)
den = sp.polyadd(sp.polymul(den2, den1), -sign * sp.polymul(num2, num1))
- return xTransferFunction(num, den)
+ return TransferFunction(num, den)
# For MIMO or SISO systems, the analytic expression is
# self / (1 - sign * other * self)
# But this does not work correctly because the state size will be too
# large.
-# This is the old TransferFunction class. It will be superceded by the
-# xTransferFunction class (which will be renamed TransferFunction) when it is
-# completed.
-class TransferFunction(signal.lti):
- """The TransferFunction class is used to represent linear
- input/output systems via its transfer function.
- """
- # Constructor
- def __init__(self, *args, **keywords):
- # First initialize the parent object
- signal.lti.__init__(self, *args, **keywords)
-
- # Make sure that this is only a SISO function
- if (self.inputs != 1 or self.outputs != 1):
- raise NotImplementedError("MIMO transfer functions not supported")
-
- # Now add a few more attributes
- self.variable = 's'
-
- # Style to use for printing (similar to MATLAB)
- def __str__(self):
- labstr = ""
- outstr = ""
- for i in range(self.inputs):
- for j in range(self.outputs):
- # Create a label for the transfer function and extract
- # numerator polynomial (depends on number of inputs/outputs)
- if (self.inputs > 1 and self.outputs > 1):
- labstr = "H[] = "; lablen = 7;
- numstr = _tfpolyToString(self.num[i,j], self.variable);
- elif (self.inputs > 1):
- labstr = "H[] = "; lablen = 7;
- numstr = _tfpolyToString(self.num[i], self.variable);
- elif (self.outputs > 1):
- labstr = "H[] = "; lablen = 7;
- numstr = _tfpolyToString(self.num[j], self.variable);
- else:
- labstr = ""; lablen = 0;
- numstr = _tfpolyToString(self.num, self.variable);
-
- # Convert the (common) denominator polynomials to strings
- denstr = _tfpolyToString(self.den, self.variable);
-
- # Figure out the length of the separating line
- dashcount = max(len(numstr), len(denstr))
- dashes = labstr + '-' * dashcount
-
- # Center the numerator or denominator
- if (len(numstr) < dashcount):
- numstr = ' ' * \
- int(round((dashcount - len(numstr))/2) + lablen) + \
- numstr
- if (len(denstr) < dashcount):
- denstr = ' ' * \
- int(round((dashcount - len(denstr))/2) + lablen) + \
- denstr
-
- outstr += "\n" + numstr + "\n" + dashes + "\n" + denstr + "\n"
- return outstr
-
- # Negation of a transfer function
- def __neg__(self):
- """Negate a transfer function"""
- return TransferFunction(-self.num, self.den)
-
- # Subtraction (use addition)
- def __sub__(self, other):
- """Subtract two transfer functions"""
- return self + (-other)
-
- def __rsub__(self, other):
- """Subtract two transfer functions"""
- return other + (-self)
-
- # Addition of two transfer functions (parallel interconnection)
- def __add__(self, sys):
- """Add two transfer functions (parallel connection)"""
- # Convert the second argument to a transfer function
- other = convertToTransferFunction(sys)
-
- # Compute the numerator and denominator of the sum
- den = sp.polymul(self.den, other.den)
- num = sp.polyadd(sp.polymul(self.num, other.den), \
- sp.polymul(other.num, self.den))
-
- return TransferFunction(num, den)
-
- # Reverse addition - just switch the order
- def __radd__(self, other):
- """Add two transfer functions (parallel connection)"""
- return self + other;
-
- # Multiplication of two transfer functions (series interconnection)
- def __mul__(self, sys):
- """Multiply two transfer functions (serial connection)"""
- # Make sure we have a transfer function (or convert to one)
- other = convertToTransferFunction(sys)
-
- # Compute the product of the transfer functions
- num = sp.polymul(self.num, other.num)
- den = sp.polymul(self.den, other.den)
- return TransferFunction(num, den)
-
- # Reverse multiplication - switch order (works for SISO)
- def __rmul__(self, other):
- """Multiply two transfer functions (serial connection)"""
- return self * other
-
- # Division between transfer functions
- def __div__(self, sys):
- """Divide two transfer functions"""
- other = convertToTransferFunction(sys);
- return TransferFunction(sp.polymul(self.num, other.den),
- sp.polymul(self.den, other.num));
-
- # Reverse division
- def __rdiv__(self, sys):
- """Divide two transfer functions"""
- other = convertToTransferFunction(sys);
- return TransferFunction(sp.polymul(other.num, self.den),
- sp.polymul(other.den, self.num));
-
- # Method for evaluating a transfer function at one frequency
- def evalfr(self, freq):
- """Evaluate a transfer function at a single frequency"""
- return sp.polyval(self.num, freq*1j) / sp.polyval(self.den, freq*1j)
-
- # Method for generating the frequency response of the system
- def freqresp(self, omega):
- """Evaluate a transfer function at a list of frequencies"""
- # Convert numerator and denomintator to 1D polynomials
- num = sp.poly1d(self.num)
- den = sp.poly1d(self.den)
-
- # Generate the frequency response at each frequency
- fresp = map(lambda w: num(w*1j) / den(w*1j), omega)
-
- mag = sp.sqrt(sp.multiply(fresp, sp.conjugate(fresp)));
- phase = sp.angle(fresp)
-
- return mag, phase, omega
-
- # Compute poles and zeros
- def poles(self): return sp.roots(self.den)
- def zeros(self): return sp.roots(self.num)
-
- # Feedback around a transfer function
- def feedback(sys1, sys2, sign=-1):
- """Feedback interconnection between two transfer functions"""
- # Get the numerator and denominator of the first system
- if (isinstance(sys1, (int, long, float, complex))):
- num1 = sys1; den1 = 1;
- elif (isinstance(sys1, TransferFunction)):
- num1 = sys1.num; den1 = sys1.den;
- else:
- raise TypeError
-
- # Get the numerator and denominator of the second system
- if (isinstance(sys2, (int, long, float, complex))):
- num2 = sys2; den2 = 1;
- elif (isinstance(sys2, TransferFunction)):
- num2 = sys2.num; den2 = sys2.den;
- else:
- raise TypeError
-
- # Compute sys1/(1 - sign*sys1*sys2)
- num = sp.polymul(num1, den2);
- den = sp.polysub(sp.polymul(den1, den2), sign * sp.polymul(num1, num2))
-
- # Return the result as a transfer function
- return TransferFunction(num, den)
-
# Utility function to convert a transfer function polynomial to a string
# Borrowed from poly1d library
def _tfpolyToString(coeffs, var='s'):
@@ -622,7 +450,7 @@
def convertToTransferFunction(sys, inputs=1, outputs=1):
"""Convert a system to transfer function form (if needed.)"""
- if isinstance(sys, xTransferFunction):
+ if isinstance(sys, TransferFunction):
return sys
elif isinstance(sys, statesp.StateSpace):
raise NotImplementedError("State space to transfer function conversion \
@@ -634,6 +462,6 @@
for i in range(min(inputs, outputs)):
num[i][i] = [sys]
- return xTransferFunction(num, den)
+ return TransferFunction(num, den)
else:
- raise TypeError("Can't convert given type to xTransferFunction system.")
+ raise TypeError("Can't convert given type to TransferFunction system.")
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:15:09
|
Revision: 65
http://python-control.svn.sourceforge.net/python-control/?rev=65&view=rev
Author: kkchen
Date: 2011-02-08 22:15:01 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Major changes regarding feedback.
- StateSpace.feedback rewritten in statesp.py, since the old one was wrong.
- bdalg.py rewritten to be simpler and to match the new feedback routines.
- TestBDAlg.py added to test bdalg.feedback.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/bdalg.py
branches/control-0.4a/src/statesp.py
Added Paths:
-----------
branches/control-0.4a/src/TestBDAlg.py
Added: branches/control-0.4a/src/TestBDAlg.py
===================================================================
--- branches/control-0.4a/src/TestBDAlg.py (rev 0)
+++ branches/control-0.4a/src/TestBDAlg.py 2011-02-08 22:15:01 UTC (rev 65)
@@ -0,0 +1,166 @@
+#!/usr/bin/env python
+
+import numpy as np
+from xferfcn import xTransferFunction
+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 = xTransferFunction([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.)
+
+ # This one doesn't work yet. The feedback system has too many states.
+ 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.)
+
+ # This one doesn't work yet. The feedback system has too many states.
+ 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.]]])
+
+if __name__ == "__main__":
+ unittest.main()
Property changes on: branches/control-0.4a/src/TestBDAlg.py
___________________________________________________________________
Added: svn:executable
+ *
Modified: branches/control-0.4a/src/bdalg.py
===================================================================
--- branches/control-0.4a/src/bdalg.py 2011-02-08 22:14:56 UTC (rev 64)
+++ branches/control-0.4a/src/bdalg.py 2011-02-08 22:15:01 UTC (rev 65)
@@ -47,86 +47,56 @@
import xferfcn as tf
import statesp as ss
-# Standard interconnections (implemented by objects)
-def series(sys1, sys2): return sys1 * sys2
-def parallel(sys1, sys2): return sys1 + sys2
-def negate(sys): return -sys;
+def feedback(sys1, sys2, sign=-1):
+ """Feedback interconnection between two I/O systems.
-# Feedback interconnection between systems
-#! TODO: This should be optimized for better performance
-#! TODO: Needs to be updated to work for state space systems
-def feedback(sys1, sys2, **keywords):
- """Feedback interconnection between two I/O systems
-
Usage
=====
- sys = feedback(sys1, sys2, **keywords)
+ sys = feedback(sys1, sys2)
+ sys = feedback(sys1, sys2, sign)
Compute the system corresponding to a feedback interconnection between
- sys1 and sys2.
+ sys1 and sys2. When sign is not specified, it assumes a value of -1
+ (negative feedback). A sign of 1 indicates positive feedback.
Parameters
----------
sys1, sys2: linsys
Linear input/output systems
+ sign: scalar
+ Feedback sign.
Return values
-------------
sys: linsys
- Keywords
- --------
- sign: float
- Sign of the interconnection (default = -1)
-
Notes
-----
- 1. This function calls calls xferfcn.feedback if the arguments are
- all either scalars or SISO transfer functions. If one or both
- of the arguments are state space systems, then the remaining
- arguments are converted to state space, as needed, and the
- statesp.feedback function is used instead. Finally, if none of
- that works, then try using sys1.feedback."""
- # Grab keyword arguments
- signparm = keywords.pop("sign", -1);
+ 1. This function calls calls xferfcn.feedback if sys1 is a TransferFunction
+ object and statesp.feedback if sys1 is a StateSpace object. If sys1 is a
+ scalar, then it is converted to sys2's type, and the corresponding
+ feedback function is used. If sys1 and sys2 are both scalars, then use
+ xferfcn.feedback."""
+
+ # Check for correct input types.
+ if not isinstance(sys1, (int, long, float, complex, tf.xTransferFunction,
+ ss.StateSpace)):
+ raise TypeError("sys1 must be a TransferFunction or StateSpace object, \
+or a scalar.")
+ if not isinstance(sys2, (int, long, float, complex, tf.xTransferFunction,
+ ss.StateSpace)):
+ raise TypeError("sys2 must be a TransferFunction or StateSpace object, \
+or a scalar.")
- #
- # Sort out which set of functions to call
- #
- # The main cases we are interested in are those where we use a
- # constant for one of the arguments to the function, in which case
- # we should use the other argument to figure out what type of
- # object we are acting on. Then call the feedback function for
- # that object.
- #
+ # If sys1 is a scalar, convert it to the appropriate LTI type so that we can
+ # its feedback member function.
+ if isinstance(sys1, (int, long, float, complex)):
+ if isinstance(sys2, tf.xTransferFunction):
+ sys1 = tf.convertToTransferFunction(sys1)
+ elif isinstance(sys2, ss.StateSpace):
+ sys1 = ss.convertToStateSpace(sys1)
+ else: # sys2 is a scalar.
+ sys1 = tf.convertToTransferFunction(sys1)
+ sys2 = tf.convertToTransferFunction(sys2)
- if (isinstance(sys1, tf.TransferFunction) and
- (isinstance(sys2, tf.TransferFunction) or
- isinstance(sys2, (int, long, float, complex)))):
- # Use transfer function feedback function
- return sys1.feedback(sys2, sign=signparm)
-
- elif (isinstance(sys2, tf.TransferFunction) and
- isinstance(sys1, (int, long, float, complex))):
- # Convert sys1 to a transfer function and then perform operation
- sys = tf.convertToTransferFunction(sys1);
- return sys.feedback(sys2, sign=signparm)
-
- elif (isinstance(sys1, ss.StateSpace) and
- (isinstance(sys2, ss.StateSpace) or
- isinstance(sys2, (int, long, float, complex)))):
- # Use state space feedback function
- return sys1.feedback(sys2, sign=signparm)
-
- elif (isinstance(sys2, ss.StateSpace) and
- isinstance(sys1, (int, long, float, complex))):
- # Convert sys1 to state space system and then perform operation
- sys = ss.convertToStateSpace(sys1);
- return sys.feedback(sys2, sign=signparm)
-
- else:
- # Assume that the first system has the right member function
- return sys1.feedback(sys2, sign=signparm)
-
- raise TypeError("can't operate on give types")
-
+ return sys1.feedback(sys2, sign)
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:14:56 UTC (rev 64)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:15:01 UTC (rev 65)
@@ -243,27 +243,36 @@
raise ValueError, "State space systems don't have compatible \
inputs/outputs for feedback."
- # note that if there is an algebraic loop then this matrix inversion
- # won't work
- # (I-D1 D2) or (I-D2 D1) will be singular
- # the easiest way to get this is to have D1 = I, D2 = I
- #! TODO: trap this error and report algebraic loop
- #! TODO: use determinant instead of inverse??
- from scipy.linalg import inv
+ from numpy.linalg import inv, det
from numpy import eye
- E21 = inv(eye(self.outputs)+sign*self.D*other.D)
- E12 = inv(eye(self.inputs)+sign*other.D*self.D)
+
+ A1 = self.A
+ B1 = self.B
+ C1 = self.C
+ D1 = self.D
+ A2 = other.A
+ B2 = other.B
+ C2 = other.C
+ D2 = other.D
- A = concatenate((
- concatenate(( self.A-sign*self.B*E12*other.D*self.C,
- -sign*self.B*E12*other.C ),axis=1),
- concatenate(( other.B*E21*self.C,
- other.A-sign*other.B*E21*self.D*other.C ),axis=1),),
- axis=0)
- B = concatenate( (self.B*E12, other.B*E21*self.D), axis=0 )
- C = concatenate( (E21*self.C, -sign*E21*self.D*other.C), axis=1 )
- D = E21*self.D
+ F = eye(self.inputs) - sign * D2 * D1
+ if abs(det(F)) < 1.e-6:
+ raise ValueError("I - sign * D2 * D1 is singular.")
+ E = inv(F)
+ T1 = eye(self.outputs) + sign * D1 * E * D2
+ T2 = eye(self.inputs) + sign * E * D2 * D1
+
+ A = concatenate(
+ (concatenate(
+ (A1 + sign * B1 * E * D2 * C1, sign * B1 * E * C2), axis=1),
+ concatenate(
+ (B2 * T1 * C1, A2 + sign * B2 * D1 * E * C2), axis=1)),
+ axis=0)
+ B = concatenate((B1 * T2, B2 * D1 * T2), axis=0)
+ C = concatenate((T1 * C1, sign * D1 * E * C2), axis=1)
+ D = D1 * T2
+
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:15:02
|
Revision: 64
http://python-control.svn.sourceforge.net/python-control/?rev=64&view=rev
Author: kkchen
Date: 2011-02-08 22:14:56 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Modified balred routine in modelsimp.py module to use real slycot routine ab09ad. User must have new slycot files installed for balred to work. testBalredTruncate unittest now passes.
Lauren Padilla <lpa...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestModelsimp.py
branches/control-0.4a/src/modelsimp.py
Modified: branches/control-0.4a/src/TestModelsimp.py
===================================================================
--- branches/control-0.4a/src/TestModelsimp.py 2011-02-08 22:14:51 UTC (rev 63)
+++ branches/control-0.4a/src/TestModelsimp.py 2011-02-08 22:14:56 UTC (rev 64)
@@ -101,14 +101,14 @@
sys = ss(A,B,C,D)
orders = 2
rsys = balred(sys,orders,'elimination','truncate')
- Artrue = np.matrix('-1.95, 0.6203; 2.314, -0.8432')
- Brtrue = np.matrix('0.7221; -0.6306')
- Crtrue = np.matrix('1.132, -0.2667')
+ 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(sys.A, Artrue)
- np.testing.assert_array_almost_equal(sys.B, Brtrue)
- np.testing.assert_array_almost_equal(sys.C, Crtrue)
- np.testing.assert_array_almost_equal(sys.D, Drtrue)
+ 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)
Modified: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py 2011-02-08 22:14:51 UTC (rev 63)
+++ branches/control-0.4a/src/modelsimp.py 2011-02-08 22:14:56 UTC (rev 64)
@@ -44,6 +44,7 @@
import ctrlutil
from control.exception import *
from statefbk import *
+from statesp import StateSpace
# 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
@@ -167,22 +168,34 @@
raise ValueError, "Oops, the system is unstable!"
if method=='matchdc':
- print "matchdc"
+ raise ValueError, "MatchDC not yet supported!"
elif method=='truncate':
- print "truncate"
+ try:
+ from slycot import ab09ad
+ except ImportError:
+ raise ControlSlycot("can't find slycot subroutine ab09ad")
+ job = 'B' # balanced (B) or not (N)
+ equil = 'N' # scale (S) or not (N)
+ ordsel = 'F' # fixed truncation level (F) or find the truncation level given tol (A)
+ n = np.size(sys.A,0)
+ m = np.size(sys.B,1)
+ p = np.size(sys.C,0)
+ nr = orders
+ tol = 0.
+ out = ab09ad(dico,job,equil,ordsel, n, m, p, nr, sys.A, sys.B, sys.C,tol)
+ Ar = out[0][0:nr,0:nr]
+ Br = out[1][0:nr,0:m]
+ Cr = out[2][0:p,0:nr]
+ hsv = out[3]
+ iwarn = out[4]
+ info = out[5]
+
+ rsys = StateSpace(Ar, Br, Cr, sys.D)
else:
raise ValueError, "Oops, method is not supported!"
- #Compute rhs using the Slycot routine XXXXXX
- #make sure Slycot is installed
- #try:
- # from slycot import XXXXXX
- #except ImportError:
- # raise ControlSlycot("can't find slycot module 'XXXXXX'")
- rsys = 0.
return rsys
-
def era(YY,m,n,nin,nout,r):
"""Calculate an ERA model of order r based on the impulse-response data YY
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:14:57
|
Revision: 63
http://python-control.svn.sourceforge.net/python-control/?rev=63&view=rev
Author: kkchen
Date: 2011-02-08 22:14:51 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Implemented xTransferFunction.feedback in xferfcn.py for SISO.
Modified bdalg.py for the new feedback functions.
Also modified xTransferFunction constructor to change the denominator to 1
whenever the numerator is 0.
Various bug fixes, e.g. in convertToTransferFunction.
Added feedback test in TestXferFcn.py.
Added NotImplementedError in StateSpace.zeros for MIMO.
Added NotImplementedError for SS -> TF and TF -> SS conversions.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestXferFcn.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:46 UTC (rev 62)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:51 UTC (rev 63)
@@ -7,7 +7,8 @@
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."""
+ formats to the xTranferFunction constructor, to try to break it. These
+ tests have been verified in MATLAB."""
# Tests for raising exceptions.
@@ -211,7 +212,7 @@
[[-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.], [-2., -1., 1.]],
+ den3 = [[[3., -2., -4], [1., 2., 3., 0., 0.], [1]],
[[-12., -9., 6., 0., 0.], [2., -1., -1], [1., 0.]]]
sys1 = xTransferFunction(num1, den1)
@@ -379,6 +380,20 @@
np.testing.assert_array_almost_equal(phase, truephase)
np.testing.assert_array_equal(omega, trueomega)
-
+ # Tests for xTransferFunction.feedback.
+
+ def testFeedbackSISO(self):
+
+ sys1 = xTransferFunction([-1., 4.], [1., 3., 5.])
+ sys2 = xTransferFunction([2., 3., 0.], [1., -3., 4., 0])
+
+ sys3 = sys1.feedback(sys2)
+ sys4 = sys1.feedback(sys2, 1)
+
+ np.testing.assert_array_equal(sys3.num, [[[-1., 7., -16., 16., 0.]]])
+ np.testing.assert_array_equal(sys3.den, [[[1., 0., -2., 2., 32., 0.]]])
+ np.testing.assert_array_equal(sys4.num, [[[-1., 7., -16., 16., 0.]]])
+ np.testing.assert_array_equal(sys4.den, [[[1., 0., 2., -8., 8., 0.]]])
+
if __name__ == "__main__":
unittest.main()
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:14:46 UTC (rev 62)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:14:51 UTC (rev 63)
@@ -120,11 +120,18 @@
# Compute poles and zeros
def poles(self):
+ """Compute the poles of a state space system."""
+
return sp.roots(sp.poly(self.A))
def zeros(self):
+ """Compute the zeros of a state space system."""
+
+ if self.inputs > 1 or self.outputs > 1:
+ raise NotImplementedError("StateSpace.zeros is currently \
+implemented only for SISO systems.")
+
den = sp.poly1d(sp.poly(self.A))
-
# Compute the numerator based on zeros
#! TODO: This is currently limited to SISO systems
num = sp.poly1d(\
@@ -273,12 +280,13 @@
# Already a state space system; just return it
return sys
elif isinstance(sys, xferfcn.TransferFunction):
- pass # TODO: convert SS to TF
+ raise NotImplementedError("Transfer function to state space conversion \
+is not implemented yet.")
elif (isinstance(sys, (int, long, float, complex))):
# Generate a simple state space system of the desired dimension
# The following Doesn't work due to inconsistencies in ltisys:
# return StateSpace([[]], [[]], [[]], sp.eye(outputs, inputs))
- return StateSpace(0, zeros((1, inputs)), zeros((outputs, 1)),
+ return StateSpace(0., zeros((1, inputs)), zeros((outputs, 1)),
sys * sp.eye(outputs, inputs))
else:
raise TypeError("Can't convert given type to StateSpace system.")
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:46 UTC (rev 62)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:51 UTC (rev 63)
@@ -108,8 +108,8 @@
raise ValueError("The numerator has %i output(s), but the \
denominator has %i\noutput(s)." % (outputs, len(den)))
- # Make sure that each row has the same number of columns.
for i in range(outputs):
+ # Make sure that each row has the same number of columns.
if len(num[i]) != inputs:
raise ValueError("Row 0 of the numerator matrix has %i \
elements, but row %i\nhas %i." % (inputs, i, len(num[i])))
@@ -117,17 +117,29 @@
raise ValueError("Row 0 of the denominator matrix has %i \
elements, but row %i\nhas %i." % (inputs, i, len(den[i])))
- # Check that we don't have any zero denominators.
+ # TODO: Right now these checks are only done during construction.
+ # It might be worthwhile to think of a way to perform checks if the
+ # user modifies the transfer function after construction.
for j in range(inputs):
- iszero = True
+ # Check that we don't have any zero denominators.
+ zeroden = True
for k in den[i][j]:
if k:
- iszero = False
+ zeroden = False
break
- if iszero:
+ if zeroden:
raise ValueError("Input %i, output %i has a zero \
denominator." % (j + 1, i + 1))
+ # If we have zero numerators, set the denominator to 1.
+ zeronum = True
+ for k in num[i][j]:
+ if k:
+ zeronum = False
+ break
+ if zeronum:
+ den[i][j] = sp.ones(1)
+
self.num = num
self.den = den
Lti2.__init__(self, inputs, outputs)
@@ -246,7 +258,7 @@
# Convert the second argument to a transfer function.
if not isinstance(other, xTransferFunction):
- other = convertToTransferFunction(other, self.outputs, self.outputs)
+ other = convertToTransferFunction(other, self.inputs, self.inputs)
# Check that the input-output sizes are consistent.
if self.inputs != other.outputs:
@@ -352,11 +364,31 @@
pass
- def feedback(sys1, sys2, sign=-1):
+ def feedback(self, other, sign=-1):
"""Feedback interconnection between two transfer functions."""
- pass
+ other = convertToTransferFunction(other)
+ if self.inputs > 1 or self.outputs > 1 or \
+ other.inputs > 1 or other.outputs > 1:
+ raise NotImplementedError("xTransferFunction.feedback is currently \
+only implemented for SISO functions.")
+
+ num1 = self.num[0][0]
+ den1 = self.den[0][0]
+ num2 = other.num[0][0]
+ den2 = other.den[0][0]
+
+ num = sp.polymul(num1, den2)
+ den = sp.polyadd(sp.polymul(den2, den1), -sign * sp.polymul(num2, num1))
+
+ return xTransferFunction(num, den)
+
+ # For MIMO or SISO systems, the analytic expression is
+ # self / (1 - sign * other * self)
+ # But this does not work correctly because the state size will be too
+ # large.
+
# This is the old TransferFunction class. It will be superceded by the
# xTransferFunction class (which will be renamed TransferFunction) when it is
# completed.
@@ -593,9 +625,15 @@
if isinstance(sys, xTransferFunction):
return sys
elif isinstance(sys, statesp.StateSpace):
- pass #TODO: convert TF to SS
+ raise NotImplementedError("State space to transfer function conversion \
+is not implemented yet.")
elif isinstance(sys, (int, long, float, complex)):
- coeff = sp.eye(outputs, inputs)
- return xTransferFunction(sys * coeff, coeff)
+ # Make an identity system.
+ num = [[[0] for j in range(inputs)] for i in range(outputs)]
+ den = [[[1] for j in range(inputs)] for i in range(outputs)]
+ for i in range(min(inputs, outputs)):
+ num[i][i] = [sys]
+
+ return xTransferFunction(num, den)
else:
- raise TypeError("Can't convert given type to StateSpace system.")
+ raise TypeError("Can't convert given type to xTransferFunction system.")
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:14:53
|
Revision: 62
http://python-control.svn.sourceforge.net/python-control/?rev=62&view=rev
Author: kkchen
Date: 2011-02-08 22:14:46 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Revised ss, tf, ss2tf, and tf2ss in matlab.py to suit our new classes.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/matlab.py
Modified: branches/control-0.4a/src/matlab.py
===================================================================
--- branches/control-0.4a/src/matlab.py 2011-02-08 22:14:42 UTC (rev 61)
+++ branches/control-0.4a/src/matlab.py 2011-02-08 22:14:46 UTC (rev 62)
@@ -60,8 +60,8 @@
# Control system library
import ctrlutil
import freqplot
-from statesp import StateSpace, rss_generate
-from xferfcn import TransferFunction
+from statesp import StateSpace, rss_generate, convertToStateSpace
+from xferfcn import TransferFunction, convertToTransferFunction
from exception import *
# Import MATLAB-like functions that can be used as-is
@@ -100,17 +100,17 @@
ss/getDelayModel - access internal delay model (state space only)
Conversions
- tf - conversion to transfer function
+* tf - conversion to transfer function
zpk - conversion to zero/pole/gain
- ss - conversion to state space
+* ss - conversion to state space
frd - conversion to frequency data
c2d - continuous to discrete conversion
d2c - discrete to continuous conversion
d2d - resample discrete-time model
upsample - upsample discrete-time LTI systems
- ss2tf - state space to transfer function
+* ss2tf - state space to transfer function
ss2zpk - transfer function to zero-pole-gain
- tf2ss - transfer function to state space
+* tf2ss - transfer function to state space
tf2zpk - transfer function to zero-pole-gain
zpk2ss - zero-pole-gain to state space
zpk2tf - zero-pole-gain to transfer function
@@ -257,54 +257,92 @@
* unwrap - unwrap a phase angle to give a continuous curve
"""
-# Create a state space system from appropriate matrices
-def ss(A, B, C, D):
- """Create a state space system from A, B, C, D"""
- return StateSpace(A, B, C, D)
+def ss(*args):
+ """Create a state space system.
-# Functions for creating a transfer function
-def tf(num, den):
- """Create a SISO transfer function given the numerator and denominator"""
- return TransferFunction(num, den)
+ Usage
+ =====
+ ss(A, B, C, D)
+ ss(sys)"""
+
+ if len(args) == 4:
+ return StateSpace(args[0], args[1], args[2], args[3])
+ elif len(args) == 1:
+ sys = args[0]
+ if isinstance(sys, StateSpace):
+ return sys
+ elif isinstance(sys, TransferFunction):
+ return tf2ss(sys)
+ else:
+ raise TypeError("ss(sys): sys must be a StateSpace or \
+TransferFunction object.")
+ else:
+ raise ValueError("Needs 1 or 4 arguments; received %i." % len(args))
-# Function for converting state space to transfer function
-def ss2tf(*args, **keywords):
- """Transform a state space system to a transfer function
+def tf(*args):
+ """Create a transfer function system.
+
+ Usage
+ =====
+ tf(num, den)
+ tf(sys)
+ num and den can be scalars or vectors for SISO systems, or lists of lists of
+ vectors for MIMO systems."""
+
+ if len(args) == 2:
+ return TransferFunction(args[0], args[1])
+ elif len(args) == 1:
+ sys = args[0]
+ if isinstance(sys, StateSpace):
+ return ss2tf(sys)
+ elif isinstance(sys, TransferFunction):
+ return sys
+ else:
+ raise TypeError("tf(sys): sys must be a StateSpace or \
+TransferFunction object.")
+ else:
+ raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
+
+def ss2tf(*args):
+ """Transform a state space system to a transfer function.
+
Usage
=====
ss2tf(A, B, C, D)
- ss2tf(sys) - sys should have attributes A, B, C, D
- """
- if (len(args) == 4):
+ ss2tf(sys) - sys should have attributes A, B, C, D"""
+
+ if len(args) == 4:
# Assume we were given the A, B, C, D matrix
- return TransferFunction(*args)
- elif (len(args) == 1):
- # Assume we were given a system object (lti or StateSpace)
+ return convertToTransferFunction(StateSpace(args[0], args[1], args[2],
+ args[3]))
+ elif len(args) == 1:
sys = args[0]
- return TransferFunction(sys.A, sys.B, sys.C, sys.D)
+ if not isinstance(sys, StateSpace):
+ raise TypeError("ss2tf(sys): sys must be a StateSpace object.")
+ return convertToTransferFunction(sys)
else:
- raise ValueError, "Needs 1 or 4 arguments."
+ raise ValueError("Needs 1 or 4 arguments; received %i." % len(args))
-# Function for converting transfer function to state space
-def tf2ss(*args, **keywords):
- """Transform a transfer function to a state space system
+def tf2ss(*args):
+ """Transform a transfer function to a state space system.
Usage
=====
tf2ss(num, den)
- ss2tf(sys) - sys should be a system object (lti or TransferFunction)
- """
- if (len(args) == 2):
+ tf2ss(sys) - sys should be a system object (lti or TransferFunction)"""
+
+ if len(args) == 2:
# Assume we were given the num, den
- return TransferFunction(*args)
- elif (len(args) == 1):
- # Assume we were given a system object (lti or TransferFunction)
+ return convertToStateSpace(TransferFunction(args[0], args[1]))
+ elif len(args) == 1:
sys = args[0]
- #! Should check to make sure object is a transfer function
- return StateSpace(sys.A, sys.B, sys.C, sys.D)
+ if not isinstance(sys, TransferFunction):
+ raise TypeError("tf2ss(sys): sys must be a TransferFunction \
+object.")
+ return convertToStateSpace(sys)
else:
- raise ValueError, "Needs 1 or 2 arguments."
+ raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
def rss(states=1, inputs=1, outputs=1):
"""Create a stable continuous random state space object."""
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:14:48
|
Revision: 61
http://python-control.svn.sourceforge.net/python-control/?rev=61&view=rev
Author: kkchen
Date: 2011-02-08 22:14:42 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Various edits:
- Put some functions in matlab.py and marked as complete.
- Minor fix to convertToStateSpace in statesp.py.
- Wrote convertToTransferFunction and put it to use in xferfcn.py.
- Changed a ValueError to TypeError; fixed in xferfcn.py and TestXferFcn.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestXferFcn.py
branches/control-0.4a/src/matlab.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:38 UTC (rev 60)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:42 UTC (rev 61)
@@ -14,7 +14,7 @@
def testBadInputType(self):
"""Give the constructor invalid input types."""
- self.assertRaises(ValueError, xTransferFunction, [[0., 1.], [2., 3.]],
+ self.assertRaises(TypeError, xTransferFunction, [[0., 1.], [2., 3.]],
[[5., 2.], [3., 0.]])
def testInconsistentDimension(self):
@@ -334,7 +334,7 @@
# Tests for xTransferFunction.freqresp.
- def testFRespSISO(self):
+ def testFreqRespSISO(self):
"""Evaluate the magnitude and phase of a SISO system at multiple
frequencies."""
@@ -351,7 +351,7 @@
np.testing.assert_array_almost_equal(phase, truephase)
np.testing.assert_array_almost_equal(omega, trueomega)
- def testFRespMIMO(self):
+ def testFreqRespMIMO(self):
"""Evaluate the magnitude and phase of a MIMO system at multiple
frequencies."""
Modified: branches/control-0.4a/src/matlab.py
===================================================================
--- branches/control-0.4a/src/matlab.py 2011-02-08 22:14:38 UTC (rev 60)
+++ branches/control-0.4a/src/matlab.py 2011-02-08 22:14:42 UTC (rev 61)
@@ -69,7 +69,7 @@
from freqplot import nyquist, nichols, gangof4
from bdalg import series, parallel, negate, feedback
from pzmap import pzmap
-from statefbk import ctrb, obsv, place, lqr
+from statefbk import ctrb, obsv, gram, place, lqr
from delay import pade
__doc__ = """
@@ -158,8 +158,8 @@
* 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
- lti/evalfr - evaluate frequency response at given frequency
+* lti/freqresp - frequency response over a frequency grid
+* lti/evalfr - evaluate frequency response at given frequency
Model simplification
minreal - minimal realization and pole/zero cancellation
@@ -193,7 +193,7 @@
canon - canonical forms of state-space models
* ctrb - controllability matrix
* obsv - observability matrix
- gram - controllability and observability gramians
+* gram - controllability and observability gramians
ss/prescale - optimal scaling of state-space models.
balreal - gramian-based input/output balancing
ss/xperm - reorder states.
@@ -317,13 +317,22 @@
return rss_generate(states, inputs, outputs, 'd')
def pole(sys):
+ """Return system poles."""
+
return sys.poles()
-# Frequency response is handled by the system object
-def freqresp(H, omega):
- """Return the frequency response for an object H at frequency omega"""
- return H.freqresp(omega)
+def evalfr(sys, omega):
+ """Evaluate the transfer function of an LTI system at a single frequency
+ omega."""
+ return sys.evalfr(omega)
+
+def freqresp(sys, omega):
+ """Return the frequency response for an LTI object at a list of frequencies
+ omega."""
+
+ return sys.freqresp(omega)
+
# Bode plots
def bode(*args, **keywords):
"""Bode plot of the frequency response
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:14:38 UTC (rev 60)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:14:42 UTC (rev 61)
@@ -267,7 +267,7 @@
# in the case of a scalar system
#
def convertToStateSpace(sys, inputs=1, outputs=1):
- """Convert a system to state space form (if needed)"""
+ """Convert a system to state space form (if needed)."""
if isinstance(sys, StateSpace):
# Already a state space system; just return it
@@ -278,11 +278,10 @@
# Generate a simple state space system of the desired dimension
# The following Doesn't work due to inconsistencies in ltisys:
# return StateSpace([[]], [[]], [[]], sp.eye(outputs, inputs))
- return StateSpace(-1, zeros((1, inputs)), zeros((outputs, 1)),
- sp.eye(outputs, inputs))
-
+ return StateSpace(0, zeros((1, inputs)), zeros((outputs, 1)),
+ sys * sp.eye(outputs, inputs))
else:
- raise TypeError("can't convert given type to StateSpace system")
+ raise TypeError("Can't convert given type to StateSpace system.")
def rss_generate(states, inputs, outputs, type):
"""This does the actual random state space generation expected from rss and
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:38 UTC (rev 60)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:42 UTC (rev 61)
@@ -91,7 +91,7 @@
else:
# If the user passed in anything else, then it's unclear what
# the meaning is.
- raise ValueError("The numerator and denominator inputs must be \
+ raise TypeError("The numerator and denominator inputs must be \
scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or \
MIMO).")
[num, den] = data
@@ -205,7 +205,7 @@
# Convert the second argument to a transfer function.
if not isinstance(other, xTransferFunction):
- other = ss2tf(other)
+ other = convertToTransferFunction(other, self.inputs, self.outputs)
# Check that the input-output sizes are consistent.
if self.inputs != other.inputs:
@@ -246,7 +246,7 @@
# Convert the second argument to a transfer function.
if not isinstance(other, xTransferFunction):
- other = ss2tf(other)
+ other = convertToTransferFunction(other, self.outputs, self.outputs)
# Check that the input-output sizes are consistent.
if self.inputs != other.outputs:
@@ -291,7 +291,7 @@
# Convert the second argument to a transfer function.
if not isinstance(other, xTransferFunction):
- other = ss2tf(other)
+ other = convertToTransferFunction(other, 1, 1)
num = sp.polymul(self.num[0][0], other.den[0][0])
den = sp.polymul(self.den[0][0], other.num[0][0])
@@ -586,8 +586,16 @@
den = sp.polymul(den1, den2)
return num, den
-
-def ss2tf(sys):
- """Convert a state space object to a transfer function object."""
- pass
+def convertToTransferFunction(sys, inputs=1, outputs=1):
+ """Convert a system to transfer function form (if needed.)"""
+
+ if isinstance(sys, xTransferFunction):
+ return sys
+ elif isinstance(sys, statesp.StateSpace):
+ pass #TODO: convert TF to SS
+ elif isinstance(sys, (int, long, float, complex)):
+ coeff = sp.eye(outputs, inputs)
+ return xTransferFunction(sys * coeff, coeff)
+ else:
+ raise TypeError("Can't convert given type to StateSpace system.")
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:14:44
|
Revision: 60
http://python-control.svn.sourceforge.net/python-control/?rev=60&view=rev
Author: kkchen
Date: 2011-02-08 22:14:38 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Implemented StateSpace.evalfr and StateSpace.freqresp in statesp.py; added
appropriate tests in TestStateSp.py.
Also made miscellaneous fixes in xferfcn.py and TestXferFcn.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestStateSp.py
branches/control-0.4a/src/TestXferFcn.py
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestStateSp.py
===================================================================
--- branches/control-0.4a/src/TestStateSp.py 2011-02-08 22:14:32 UTC (rev 59)
+++ branches/control-0.4a/src/TestStateSp.py 2011-02-08 22:14:38 UTC (rev 60)
@@ -2,8 +2,53 @@
import numpy as np
import matlab
+from statesp import StateSpace
import unittest
+class TestStateSpace(unittest.TestCase):
+ """Tests for the StateSpace class."""
+
+ 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."""
@@ -73,4 +118,4 @@
self.assertTrue(abs(z) < 1)
if __name__ == "__main__":
- unittest.main()
\ No newline at end of file
+ unittest.main()
Modified: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:32 UTC (rev 59)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:38 UTC (rev 60)
@@ -377,7 +377,7 @@
np.testing.assert_array_almost_equal(mag, truemag)
np.testing.assert_array_almost_equal(phase, truephase)
- np.testing.assert_array_almost_equal(omega, trueomega)
+ np.testing.assert_array_equal(omega, trueomega)
if __name__ == "__main__":
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:14:32 UTC (rev 59)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:14:38 UTC (rev 60)
@@ -43,6 +43,7 @@
import scipy as sp
from scipy import concatenate, zeros
+from numpy.linalg import solve
import xferfcn
from lti2 import Lti2
@@ -92,34 +93,31 @@
str += "D = " + self.D.__str__() + "\n"
return str
+ def evalfr(self, freq):
+ """Method for evaluating a system at one frequency."""
+
+ fresp = self.C * solve(freq * 1.j * sp.eye(self.states) - self.A,
+ self.B) + self.D
+ return fresp
+
# Method for generating the frequency response of the system
def freqresp(self, omega=None):
"""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
+ # Preallocate outputs.
+ numfreq = len(omega)
+ mag = sp.empty((self.outputs, self.inputs, numfreq))
+ phase = sp.empty((self.outputs, self.inputs, numfreq))
+ fresp = sp.empty((self.outputs, self.inputs, numfreq), dtype=complex)
- # Compute the denominator from the A matrix
- den = sp.poly1d(sp.poly(self.A))
+ for k in range(numfreq):
+ fresp[:, :, k] = self.evalfr(omega[k])
- # Compute the numerator based on zeros
- #! TODO: This is currently limited to SISO systems
- num = sp.poly1d(\
- sp.poly(self.A - sp.dot(self.B, self.C)) + (self.D[0] - 1) * den)
-
- # Generate the frequency response at each frequency
- fresp = map(lambda w: num(w*1j) / den(w*1j), omega)
- mag = sp.sqrt(sp.multiply(fresp, sp.conjugate(fresp)))
+ mag = abs(fresp)
phase = sp.angle(fresp)
return mag, phase, omega
- # Method for evaluating a system at one frequency
- def evalfr(self, freq):
- #! TODO: Not implemented
- return None
-
# Compute poles and zeros
def poles(self):
return sp.roots(sp.poly(self.A))
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:32 UTC (rev 59)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:38 UTC (rev 60)
@@ -145,8 +145,6 @@
if mimo:
outstr += "\nInput %i to output %i:" % (i + 1, j + 1)
- lablen = 0
-
# Convert the numerator and denominator polynomials to strings.
numstr = _tfpolyToString(self.num[j][i]);
denstr = _tfpolyToString(self.den[j][i]);
@@ -158,11 +156,11 @@
# Center the numerator or denominator
if (len(numstr) < dashcount):
numstr = ' ' * \
- int(round((dashcount - len(numstr))/2) + lablen) + \
+ int(round((dashcount - len(numstr))/2)) + \
numstr
if (len(denstr) < dashcount):
denstr = ' ' * \
- int(round((dashcount - len(denstr))/2) + lablen) + \
+ int(round((dashcount - len(denstr))/2)) + \
denstr
outstr += "\n" + numstr + "\n" + dashes + "\n" + denstr + "\n"
@@ -229,7 +227,7 @@
return xTransferFunction(num, den)
def __radd__(self, other):
- """Add two transfer functions (parallel connection)"""
+ """Add two transfer functions (parallel connection)."""
return self + other;
@@ -244,7 +242,7 @@
return other + (-self)
def __mul__(self, other):
- """Multiply two transfer functions (serial connection)"""
+ """Multiply two transfer functions (serial connection)."""
# Convert the second argument to a transfer function.
if not isinstance(other, xTransferFunction):
@@ -278,13 +276,13 @@
return xTransferFunction(num, den)
def __rmul__(self, other):
- """Multiply two transfer functions (serial connection)"""
+ """Multiply two transfer functions (serial connection)."""
return self * other
# TODO: Division of MIMO transfer function objects is quite difficult.
def __div__(self, other):
- """Divide two transfer functions"""
+ """Divide two transfer functions."""
if self.inputs > 1 or self.outputs > 1 or \
other.inputs > 1 or other.outputs > 1:
@@ -302,7 +300,7 @@
# TODO: Division of MIMO transfer function objects is quite difficult.
def __rdiv__(self, other):
- """Reverse divide two transfer functions"""
+ """Reverse divide two transfer functions."""
if self.inputs > 1 or self.outputs > 1 or \
other.inputs > 1 or other.outputs > 1:
@@ -312,7 +310,7 @@
return other / self
def evalfr(self, freq):
- """Evaluate a transfer function at a single frequency"""
+ """Evaluate a transfer function at a single frequency."""
# Preallocate the output.
out = sp.empty((self.outputs, self.inputs), dtype=complex)
@@ -325,23 +323,22 @@
return out
# Method for generating the frequency response of the system
- def freqresp(self, omega):
- """Evaluate a transfer function at a list of frequencies"""
+ def freqresp(self, omega=None):
+ """Evaluate a transfer function at a list of frequencies."""
+ # Preallocate outputs.
numfreq = len(omega)
+ mag = sp.empty((self.outputs, self.inputs, numfreq))
+ phase = sp.empty((self.outputs, self.inputs, numfreq))
- # Preallocate outputs.
- mag = sp.empty((self.outputs, self.inputs, numfreq), dtype=complex)
- phase = sp.empty((self.outputs, self.inputs, numfreq), dtype=complex)
-
for i in range(self.outputs):
for j in range(self.inputs):
fresp = map(lambda w: sp.polyval(self.num[i][j], w * 1.j) / \
sp.polyval(self.den[i][j], w * 1.j), omega)
fresp = sp.array(fresp)
- mag[i][j] = abs(fresp)
- phase[i][j] = sp.angle(fresp)
+ mag[i, j] = abs(fresp)
+ phase[i, j] = sp.angle(fresp)
return mag, phase, omega
@@ -356,7 +353,7 @@
pass
def feedback(sys1, sys2, sign=-1):
- """Feedback interconnection between two transfer functions"""
+ """Feedback interconnection between two transfer functions."""
pass
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:14:38
|
Revision: 59
http://python-control.svn.sourceforge.net/python-control/?rev=59&view=rev
Author: kkchen
Date: 2011-02-08 22:14:32 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added xTransferFunction.evalfr and xTransferFunction.freqresp in xferfcn.py.
Also added unit tests for these in TextXferFcn.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestXferFcn.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:27 UTC (rev 58)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:32 UTC (rev 59)
@@ -306,5 +306,79 @@
np.testing.assert_array_equal(sys4.num, sys3.den)
np.testing.assert_array_equal(sys4.den, sys3.num)
+ # Tests for xTransferFunction.evalfr.
+
+ def testEvalFrSISO(self):
+ """Evaluate the frequency response of a SISO system at one frequency."""
+
+ sys = xTransferFunction([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 = xTransferFunction(num, den)
+ resp = [[0.147058823529412 + 0.0882352941176471j, -0.75, 1.],
+ [-0.083333333333333, -0.188235294117647 - 0.847058823529412j,
+ -1. - 8.j]]
+
+ np.testing.assert_array_almost_equal(sys.evalfr(2.), resp)
+
+ # Tests for xTransferFunction.freqresp.
+
+ def testFRespSISO(self):
+ """Evaluate the magnitude and phase of a SISO system at multiple
+ frequencies."""
+
+ sys = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
+
+ truemag = [[[4.63507337473906, 0.707106781186548, 0.0866592803995351]]]
+ truephase = [[[-2.89596891081488, -2.35619449019234,
+ -1.32655885133871]]]
+ trueomega = [0.1, 1., 10.]
+
+ mag, phase, omega = sys.freqresp(trueomega)
+
+ np.testing.assert_array_almost_equal(mag, truemag)
+ np.testing.assert_array_almost_equal(phase, truephase)
+ np.testing.assert_array_almost_equal(omega, trueomega)
+
+ def testFRespMIMO(self):
+ """Evaluate the magnitude and phase of a MIMO system at multiple
+ frequencies."""
+
+ num = [[[1., 2.], [0., 3.], [2., -1.]],
+ [[1.], [4., 0.], [1., -4., 3.]]]
+ den = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]],
+ [[3., 0., .0], [2., -1., -1.], [1.]]]
+ sys = xTransferFunction(num, den)
+
+ trueomega = [0.1, 1., 10.]
+ truemag = [[[0.496287094505259, 0.307147558416976, 0.0334738176210382],
+ [300., 3., 0.03], [1., 1., 1.]],
+ [[33.3333333333333, 0.333333333333333, 0.00333333333333333],
+ [0.390285696125482, 1.26491106406735, 0.198759144198533],
+ [3.01663720059274, 4.47213595499958, 104.92378186093]]]
+ truephase = [[[3.7128711165168e-4, 0.185347949995695, 1.30770596539255],
+ [-np.pi, -np.pi, -np.pi], [0., 0., 0.]],
+ [[-np.pi, -np.pi, -np.pi],
+ [-1.66852323415362, -1.89254688119154, -1.62050658356412],
+ [-0.132989648369409, -1.1071487177940, -2.7504672066207]]]
+
+ mag, phase, omega = sys.freqresp(trueomega)
+
+ np.testing.assert_array_almost_equal(mag, truemag)
+ np.testing.assert_array_almost_equal(phase, truephase)
+ np.testing.assert_array_almost_equal(omega, trueomega)
+
+
if __name__ == "__main__":
- unittest.main()
\ No newline at end of file
+ unittest.main()
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:27 UTC (rev 58)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:32 UTC (rev 59)
@@ -314,15 +314,37 @@
def evalfr(self, freq):
"""Evaluate a transfer function at a single frequency"""
- # return sp.polyval(self.num, freq*1j) / sp.polyval(self.den, freq*1j)
- pass
+ # Preallocate the output.
+ out = sp.empty((self.outputs, self.inputs), dtype=complex)
+ for i in range(self.outputs):
+ for j in range(self.inputs):
+ out[i][j] = sp.polyval(self.num[i][j], freq * 1.j) / \
+ sp.polyval(self.den[i][j], freq * 1.j)
+
+ return out
+
# Method for generating the frequency response of the system
def freqresp(self, omega):
"""Evaluate a transfer function at a list of frequencies"""
- pass
+ numfreq = len(omega)
+ # Preallocate outputs.
+ mag = sp.empty((self.outputs, self.inputs, numfreq), dtype=complex)
+ phase = sp.empty((self.outputs, self.inputs, numfreq), dtype=complex)
+
+ for i in range(self.outputs):
+ for j in range(self.inputs):
+ fresp = map(lambda w: sp.polyval(self.num[i][j], w * 1.j) / \
+ sp.polyval(self.den[i][j], w * 1.j), omega)
+ fresp = sp.array(fresp)
+
+ mag[i][j] = abs(fresp)
+ phase[i][j] = sp.angle(fresp)
+
+ return mag, phase, omega
+
def poles(self):
"""Compute poles of a transfer function."""
@@ -571,4 +593,4 @@
def ss2tf(sys):
"""Convert a state space object to a transfer function object."""
- pass
\ No newline at end of file
+ pass
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:14:33
|
Revision: 58
http://python-control.svn.sourceforge.net/python-control/?rev=58&view=rev
Author: kkchen
Date: 2011-02-08 22:14:27 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Bug fix in xTransferFunction._truncatecoeff; miscellaneous other fixes in xferfcn.py.
Also added tests for exceptions in TestXferFcn.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestXferFcn.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:22 UTC (rev 57)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:27 UTC (rev 58)
@@ -9,14 +9,85 @@
function class. Throughout these tests, we will give different input
formats to the xTranferFunction constructor, to try to break it."""
- def testTruncateCoeff(self):
+ # Tests for raising exceptions.
+
+ def testBadInputType(self):
+ """Give the constructor invalid input types."""
+
+ self.assertRaises(ValueError, xTransferFunction, [[0., 1.], [2., 3.]],
+ [[5., 2.], [3., 0.]])
+
+ def testInconsistentDimension(self):
+ """Give the constructor a numerator and denominator of different
+ sizes."""
+
+ self.assertRaises(ValueError, xTransferFunction, [[[1.]]],
+ [[[1.], [2., 3.]]])
+ self.assertRaises(ValueError, xTransferFunction, [[[1.]]],
+ [[[1.]], [[2., 3.]]])
+ self.assertRaises(ValueError, xTransferFunction, [[[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, xTransferFunction, 1.,
+ [[[1.]], [[2.], [3.]]])
+ self.assertRaises(ValueError, xTransferFunction, [[[1.]], [[2.], [3.]]],
+ 1.)
+
+ def testZeroDenominator(self):
+ """Give the constructor a transfer function with a zero denominator."""
+
+ self.assertRaises(ValueError, xTransferFunction, 1., 0.)
+ self.assertRaises(ValueError, xTransferFunction,
+ [[[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 = xTransferFunction([[[1., 2.]]], [[[4., 5.]]])
+ sys2 = xTransferFunction([[[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 = xTransferFunction([[[1., 2.], [4., 5.]], [[2., 5.], [4., 3.]]],
+ [[[6., 2.], [4., 1.]], [[6., 7.], [2., 4.]]])
+ sys2 = xTransferFunction([[[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 xTransferFunction._truncatecoeff
+
+ def testTruncateCoeff1(self):
"""Remove extraneous zeros in polynomial representations."""
sys1 = xTransferFunction([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 = xTransferFunction([0., 0., 0.], 1.)
+
+ np.testing.assert_array_equal(sys1.num, [[[0.]]])
+ np.testing.assert_array_equal(sys1.den, [[[1.]]])
+ # Tests for xTransferFunction.__neg__
+
def testNegScalar(self):
"""Negate a direct feedthrough system."""
@@ -53,7 +124,9 @@
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 xTransferFunction.__add__
+
def testAddScalar(self):
"""Add two direct feedthrough systems."""
@@ -100,6 +173,8 @@
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 xTransferFunction.__sub__
+
def testSubScalar(self):
"""Add two direct feedthrough systems."""
@@ -134,20 +209,22 @@
[[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.]]]
+ 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.], [-2., -1., 1.]],
+ [[-12., -9., 6., 0., 0.], [2., -1., -1], [1., 0.]]]
sys1 = xTransferFunction(num1, den1)
sys2 = xTransferFunction(num2, den2)
- sys3 = sys1 + sys2
+ 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 xTransferFunction.__mul__
+
def testMulScalar(self):
"""Multiply two direct feedthrough systems."""
@@ -204,6 +281,8 @@
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 xTransferFunction.__div__
+
def testDivScalar(self):
"""Divide two direct feedthrough systems."""
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:22 UTC (rev 57)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:27 UTC (rev 58)
@@ -172,18 +172,24 @@
"""Remove extraneous zero coefficients from polynomials in numerator and
denominator matrices."""
- # Beware: this is a shallow copy.
+ # Beware: this is a shallow copy. This should be okay.
data = [self.num, self.den]
for p in range(len(data)):
for i in range(self.outputs):
for j in range(self.inputs):
# Find the first nontrivial coefficient.
- k = 0
- while not data[p][i][j][k]:
- k += 1
-
- # Now truncate the trivial coefficients.
- data[p][i][j] = data[p][i][j][k:]
+ nonzero = None
+ for k in range(data[p][i][j].size):
+ if data[p][i][j][k]:
+ nonzero = k
+ break
+
+ if nonzero is None:
+ # The array is all zeros.
+ data[p][i][j] = sp.zeros(1)
+ else:
+ # Truncate the trivial coefficients.
+ data[p][i][j] = data[p][i][j][nonzero:]
[self.num, self.den] = data
def __neg__(self):
@@ -205,11 +211,11 @@
# Check that the input-output sizes are consistent.
if self.inputs != other.inputs:
- raise ValueError("The first summand has %i input(s), but the second \
-has %i." % (self.inputs, other.inputs))
+ raise ValueError("The first summand has %i input(s), but the \
+second has %i." % (self.inputs, other.inputs))
if self.outputs != other.outputs:
- raise ValueError("The first summand has %i output(s), but the second \
-has %i." % (self.outputs, other.outputs))
+ raise ValueError("The first summand has %i output(s), but the \
+second has %i." % (self.outputs, other.outputs))
# Preallocate the numerator and denominator of the sum.
num = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
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:14:28
|
Revision: 57
http://python-control.svn.sourceforge.net/python-control/?rev=57&view=rev
Author: kkchen
Date: 2011-02-08 22:14:22 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added tests to TestXferFcn.py. Changes in xferfcn.py:
- Shallow copy bug fix in xTransferFunction.__neg__
- Coefficient vectors are now converted to numpy arrays.
- Simplification of xTrasferFunction.__rdiv__
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestXferFcn.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:17 UTC (rev 56)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:22 UTC (rev 57)
@@ -6,7 +6,8 @@
class TestXferFcn(unittest.TestCase):
"""These are tests for functionality and correct reporting of the transfer
- function class."""
+ function class. Throughout these tests, we will give different input
+ formats to the xTranferFunction constructor, to try to break it."""
def testTruncateCoeff(self):
"""Remove extraneous zeros in polynomial representations."""
@@ -16,23 +17,58 @@
np.testing.assert_array_equal(sys1.num, [[[1., 2.]]])
np.testing.assert_array_equal(sys1.den, [[[3., 2., 1.]]])
- def testAddSISO1(self):
+ def testNegScalar(self):
+ """Negate a direct feedthrough system."""
+
+ sys1 = xTransferFunction(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 = xTransferFunction([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 = xTransferFunction(num1, den1)
+ sys2 = - sys1
+ sys3 = xTransferFunction(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])
+
+ def testAddScalar(self):
"""Add two direct feedthrough systems."""
- # Try different input formats, too.
sys1 = xTransferFunction(1., [[[1.]]])
- sys2 = xTransferFunction([2.], [1.])
+ sys2 = xTransferFunction(np.array([2.]), [1.])
sys3 = sys1 + sys2
np.testing.assert_array_equal(sys3.num, 3.)
np.testing.assert_array_equal(sys3.den, 1.)
- def testAddSISO2(self):
+ def testAddSISO(self):
"""Add two SISO systems."""
- # Try different input formats, too.
sys1 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
- sys2 = xTransferFunction([[[-1., 3.]]], [[[1., 0., -1.]]])
+ sys2 = xTransferFunction([[np.array([-1., 3.])]], [[[1., 0., -1.]]])
sys3 = sys1 + sys2
# If sys3.num is [[[0., 20., 4., -8.]]], then this is wrong!
@@ -63,16 +99,133 @@
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])
+
+ def testSubScalar(self):
+ """Add two direct feedthrough systems."""
+
+ sys1 = xTransferFunction(1., [[[1.]]])
+ sys2 = xTransferFunction(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 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys2 = xTransferFunction([[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., -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.]]]
- def testMulSISO1(self):
+ sys1 = xTransferFunction(num1, den1)
+ sys2 = xTransferFunction(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])
+
+ def testMulScalar(self):
"""Multiply two direct feedthrough systems."""
sys1 = xTransferFunction(2., [1.])
sys2 = xTransferFunction(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 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys2 = xTransferFunction([[[-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 = xTransferFunction(num1, den1)
+ sys2 = xTransferFunction(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])
+ def testDivScalar(self):
+ """Divide two direct feedthrough systems."""
+
+ sys1 = xTransferFunction(np.array([3.]), -4.)
+ sys2 = xTransferFunction(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 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys2 = xTransferFunction([[[-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)
+
if __name__ == "__main__":
unittest.main()
\ No newline at end of file
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:17 UTC (rev 56)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:22 UTC (rev 57)
@@ -49,6 +49,7 @@
# External function declarations
import scipy as sp
import scipy.signal as signal
+import copy
import bdalg as bd
import statesp
from lti2 import Lti2
@@ -58,7 +59,7 @@
main data members are 'num' and 'den', which are 2-D lists of arrays
containing MIMO numerator and denominator coefficients. For example,
- >>> num[2][5] = [1., 4., 8.]
+ >>> num[2][5] = numpy.array([1., 4., 8.])
means that the numerator of the transfer function from the 6th input to the
3rd output is set to s^2 + 4s + 8."""
@@ -67,27 +68,32 @@
"""This is the constructor. The default transfer function is 1 (unit
gain direct feedthrough)."""
- # Make num and den into lists of lists of arrays, if necessary.
+ # Make num and den into lists of lists of arrays, if necessary. Beware:
+ # this is a shallow copy! This should be okay, but be careful.
data = [num, den]
for i in range(len(data)):
if isinstance(data[i], (int, float, long, complex)):
# Convert scalar to list of list of array.
- data[i] = [[[data[i]]]]
+ data[i] = [[sp.array([data[i]])]]
elif isinstance(data[i], (list, tuple, sp.ndarray)) and \
isinstance(data[i][0], (int, float, long, complex)):
# Convert array to list of list of array.
- data[i] = [[data[i]]]
- elif isinstance(data[i], (list, tuple, sp.ndarray)) and \
- isinstance(data[i][0], (list, tuple, sp.ndarray)) and \
+ data[i] = [[sp.array(data[i])]]
+ elif isinstance(data[i], list) and \
+ isinstance(data[i][0], list) and \
isinstance(data[i][0][0], (list, tuple, sp.ndarray)) and \
isinstance(data[i][0][0][0], (int, float, long, complex)):
- # We already have the right format.
- pass
+ # We might already have the right format. Convert the
+ # coefficient vectors to arrays, if necessary.
+ for j in range(len(data[i])):
+ for k in range(len(data[i][j])):
+ data[i][j][k] = sp.array(data[i][j][k])
else:
# If the user passed in anything else, then it's unclear what
# the meaning is.
raise ValueError("The numerator and denominator inputs must be \
-scalars or arrays (for\nSISO), or lists of lists of arrays (for MIMO).")
+scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or \
+MIMO).")
[num, den] = data
inputs = len(num[0])
@@ -166,8 +172,8 @@
"""Remove extraneous zero coefficients from polynomials in numerator and
denominator matrices."""
+ # Beware: this is a shallow copy.
data = [self.num, self.den]
-
for p in range(len(data)):
for i in range(self.outputs):
for j in range(self.inputs):
@@ -177,15 +183,19 @@
k += 1
# Now truncate the trivial coefficients.
- data[p][i][j] = data[p][i][j][k:]
-
+ data[p][i][j] = data[p][i][j][k:]
[self.num, self.den] = data
def __neg__(self):
"""Negate a transfer function."""
- return xTransferFunction(-self.num, self.den)
+ num = copy.deepcopy(self.num)
+ for i in range(self.outputs):
+ for j in range(self.inputs):
+ num[i][j] *= -1
+ return xTransferFunction(num, self.den)
+
def __add__(self, other):
"""Add two transfer functions (parallel connection)."""
@@ -285,7 +295,7 @@
return xTransferFunction(num, den)
# TODO: Division of MIMO transfer function objects is quite difficult.
- def __rdiv__(self, sys):
+ def __rdiv__(self, other):
"""Reverse divide two transfer functions"""
if self.inputs > 1 or self.outputs > 1 or \
@@ -293,15 +303,8 @@
raise NotImplementedError("xTransferFunction.__rdiv__ is currently \
implemented only for SISO systems.")
- # Convert the second argument to a transfer function.
- if not isinstance(other, xTransferFunction):
- other = ss2tf(other)
-
- num = sp.polymul(self.den[0][0], other.num[0][0])
- den = sp.polymul(self.num[0][0], other.den[0][0])
+ return other / self
- return xTransferFunction(num, den)
-
def evalfr(self, freq):
"""Evaluate a transfer function at a single frequency"""
@@ -329,6 +332,9 @@
pass
+# This is the old TransferFunction class. It will be superceded by the
+# xTransferFunction class (which will be renamed TransferFunction) when it is
+# completed.
class TransferFunction(signal.lti):
"""The TransferFunction class is used to represent linear
input/output systems via its transfer function.
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:14:23
|
Revision: 56
http://python-control.svn.sourceforge.net/python-control/?rev=56&view=rev
Author: kkchen
Date: 2011-02-08 22:14:17 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Changes to xferfcn.py:
- Check that we don't have any zero denominators.
- Bug fix in _addSISO
- Exceptions more specific
Also added a simple SISO multiplication test to TestXferFcn.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestXferFcn.py
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:13 UTC (rev 55)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:17 UTC (rev 56)
@@ -57,13 +57,22 @@
sys1 = xTransferFunction(num1, den1)
sys2 = xTransferFunction(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])
+
+ def testMulSISO1(self):
+ """Multiply two direct feedthrough systems."""
+
+ sys1 = xTransferFunction(2., [1.])
+ sys2 = xTransferFunction(1., 4.)
+ sys3 = sys1 * sys2
+
+ np.testing.assert_array_equal(sys3.num, [[[2.]]])
+ np.testing.assert_array_equal(sys3.den, [[[4.]]])
if __name__ == "__main__":
unittest.main()
\ No newline at end of file
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:13 UTC (rev 55)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:17 UTC (rev 56)
@@ -67,7 +67,7 @@
"""This is the constructor. The default transfer function is 1 (unit
gain direct feedthrough)."""
- # Make num and den into 3-d numpy arrays, if necessary.
+ # Make num and den into lists of lists of arrays, if necessary.
data = [num, den]
for i in range(len(data)):
if isinstance(data[i], (int, float, long, complex)):
@@ -93,19 +93,34 @@
inputs = len(num[0])
outputs = len(num)
+ # Make sure the numerator and denominator matrices have consistent
+ # sizes.
if inputs != len(den[0]):
- raise ValueError("The numerator and denominator matrices must have \
-the same column\n(input) size.")
+ raise ValueError("The numerator has %i input(s), but the \
+denominator has %i\ninput(s)." % (inputs, len(den[0])))
if outputs != len(den):
- raise ValueError("The numerator and denominator matrices must have \
-the same row\n(output) size.")
- for i in range(1, outputs):
+ raise ValueError("The numerator has %i output(s), but the \
+denominator has %i\noutput(s)." % (outputs, len(den)))
+
+ # Make sure that each row has the same number of columns.
+ for i in range(outputs):
if len(num[i]) != inputs:
- raise ValueError("Each row of the numerator matrix must have \
-the same number of\nelements.")
+ raise ValueError("Row 0 of the numerator matrix has %i \
+elements, but row %i\nhas %i." % (inputs, i, len(num[i])))
if len(den[i]) != inputs:
- raise ValueError("Each row of the denominator matrix must have \
-the same number of\nelements.")
+ raise ValueError("Row 0 of the denominator matrix has %i \
+elements, but row %i\nhas %i." % (inputs, i, len(den[i])))
+
+ # Check that we don't have any zero denominators.
+ for j in range(inputs):
+ iszero = True
+ for k in den[i][j]:
+ if k:
+ iszero = False
+ break
+ if iszero:
+ raise ValueError("Input %i, output %i has a zero \
+denominator." % (j + 1, i + 1))
self.num = num
self.den = den
@@ -180,11 +195,11 @@
# Check that the input-output sizes are consistent.
if self.inputs != other.inputs:
- raise ValueError("The two systems to be added must have the same \
-input size.")
+ raise ValueError("The first summand has %i input(s), but the second \
+has %i." % (self.inputs, other.inputs))
if self.outputs != other.outputs:
- raise ValueError("The two systems to be added must have the same \
-output size.")
+ raise ValueError("The first summand has %i output(s), but the second \
+has %i." % (self.outputs, other.outputs))
# Preallocate the numerator and denominator of the sum.
num = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
@@ -192,7 +207,7 @@
for i in range(self.outputs):
for j in range(self.inputs):
- num[i][j], den[i][j] = addSISO(self.num[i][j], self.den[i][j],
+ num[i][j], den[i][j] = _addSISO(self.num[i][j], self.den[i][j],
other.num[i][j], other.den[i][j])
return xTransferFunction(num, den)
@@ -221,15 +236,15 @@
# Check that the input-output sizes are consistent.
if self.inputs != other.outputs:
- raise ValueError("C = A * B: A must have the same number of \
-columns (inputs) as B has\nrows (outputs).")
+ raise ValueError("C = A * B: A has %i column(s) (input(s)), but B \
+has %i row(s)\n(output(s))." % (self.inputs, other.outputs))
inputs = other.inputs
outputs = self.outputs
# Preallocate the numerator and denominator of the sum.
num = [[[0] for j in range(inputs)] for i in range(outputs)]
- den = [[[0] for j in range(inputs)] for i in range(outputs)]
+ den = [[[1] for j in range(inputs)] for i in range(outputs)]
# Temporary storage for the summands needed to find the (i, j)th element
# of the product.
@@ -241,7 +256,7 @@
for k in range(self.inputs): # Multiply & add.
num_summand[k] = sp.polymul(self.num[i][k], other.num[k][j])
den_summand[k] = sp.polymul(self.den[i][k], other.den[k][j])
- num[i][j], den[i][j] = addSISO(num[i][j], den[i][j],
+ num[i][j], den[i][j] = _addSISO(num[i][j], den[i][j],
num_summand[k], den_summand[k])
return xTransferFunction(num, den)
@@ -251,16 +266,42 @@
return self * other
- def __div__(self, sys):
+ # TODO: Division of MIMO transfer function objects is quite difficult.
+ def __div__(self, other):
"""Divide two transfer functions"""
- pass
+ if self.inputs > 1 or self.outputs > 1 or \
+ other.inputs > 1 or other.outputs > 1:
+ raise NotImplementedError("xTransferFunction.__div__ is currently \
+implemented only for SISO systems.")
+
+ # Convert the second argument to a transfer function.
+ if not isinstance(other, xTransferFunction):
+ other = ss2tf(other)
+
+ num = sp.polymul(self.num[0][0], other.den[0][0])
+ den = sp.polymul(self.den[0][0], other.num[0][0])
+ return xTransferFunction(num, den)
+
+ # TODO: Division of MIMO transfer function objects is quite difficult.
def __rdiv__(self, sys):
"""Reverse divide two transfer functions"""
- pass
+ if self.inputs > 1 or self.outputs > 1 or \
+ other.inputs > 1 or other.outputs > 1:
+ raise NotImplementedError("xTransferFunction.__rdiv__ is currently \
+implemented only for SISO systems.")
+
+ # Convert the second argument to a transfer function.
+ if not isinstance(other, xTransferFunction):
+ other = ss2tf(other)
+
+ num = sp.polymul(self.den[0][0], other.num[0][0])
+ den = sp.polymul(self.num[0][0], other.den[0][0])
+ return xTransferFunction(num, den)
+
def evalfr(self, freq):
"""Evaluate a transfer function at a single frequency"""
@@ -354,6 +395,7 @@
def __sub__(self, other):
"""Subtract two transfer functions"""
return self + (-other)
+
def __rsub__(self, other):
"""Subtract two transfer functions"""
return other + (-self)
@@ -505,7 +547,7 @@
thestr = newstr
return thestr
-def addSISO(num1, den1, num2, den2):
+def _addSISO(num1, den1, num2, den2):
"""Return num/den = num1/den1 + num2/den2, where each numerator and
denominator is a list of polynomial coefficients."""
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:14:19
|
Revision: 55
http://python-control.svn.sourceforge.net/python-control/?rev=55&view=rev
Author: kkchen
Date: 2011-02-08 22:14:13 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added xTransferFunction.__str__ in xferfcn.py. Appears to work well.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:08 UTC (rev 54)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:13 UTC (rev 55)
@@ -111,14 +111,43 @@
self.den = den
Lti2.__init__(self, inputs, outputs)
- self.truncatecoeff()
+ self._truncatecoeff()
def __str__(self):
"""String representation of the transfer function."""
+
+ mimo = self.inputs > 1 or self.outputs > 1
+ outstr = ""
+
+ for i in range(self.inputs):
+ for j in range(self.outputs):
+ if mimo:
+ outstr += "\nInput %i to output %i:" % (i + 1, j + 1)
+
+ lablen = 0
+
+ # Convert the numerator and denominator polynomials to strings.
+ numstr = _tfpolyToString(self.num[j][i]);
+ denstr = _tfpolyToString(self.den[j][i]);
+
+ # Figure out the length of the separating line
+ dashcount = max(len(numstr), len(denstr))
+ dashes = '-' * dashcount
+
+ # Center the numerator or denominator
+ if (len(numstr) < dashcount):
+ numstr = ' ' * \
+ int(round((dashcount - len(numstr))/2) + lablen) + \
+ numstr
+ if (len(denstr) < dashcount):
+ denstr = ' ' * \
+ int(round((dashcount - len(denstr))/2) + lablen) + \
+ denstr
+
+ outstr += "\n" + numstr + "\n" + dashes + "\n" + denstr + "\n"
+ return outstr
- pass
-
- def truncatecoeff(self):
+ def _truncatecoeff(self):
"""Remove extraneous zero coefficients from polynomials in numerator and
denominator matrices."""
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:14:14
|
Revision: 54
http://python-control.svn.sourceforge.net/python-control/?rev=54&view=rev
Author: kkchen
Date: 2011-02-08 22:14:08 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Worked on xTransferFunction.__add__ and __mul__ in xferfcn.py.
Also added tests in TestXferFcn.py.
Minor edit in statesp.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Added Paths:
-----------
branches/control-0.4a/src/TestXferFcn.py
Added: branches/control-0.4a/src/TestXferFcn.py
===================================================================
--- branches/control-0.4a/src/TestXferFcn.py (rev 0)
+++ branches/control-0.4a/src/TestXferFcn.py 2011-02-08 22:14:08 UTC (rev 54)
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+
+import numpy as np
+from xferfcn import xTransferFunction
+import unittest
+
+class TestXferFcn(unittest.TestCase):
+ """These are tests for functionality and correct reporting of the transfer
+ function class."""
+
+ def testTruncateCoeff(self):
+ """Remove extraneous zeros in polynomial representations."""
+
+ sys1 = xTransferFunction([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 testAddSISO1(self):
+ """Add two direct feedthrough systems."""
+
+ # Try different input formats, too.
+ sys1 = xTransferFunction(1., [[[1.]]])
+ sys2 = xTransferFunction([2.], [1.])
+ sys3 = sys1 + sys2
+
+ np.testing.assert_array_equal(sys3.num, 3.)
+ np.testing.assert_array_equal(sys3.den, 1.)
+
+ def testAddSISO2(self):
+ """Add two SISO systems."""
+
+ # Try different input formats, too.
+ sys1 = xTransferFunction([1., 3., 5], [1., 6., 2., -1])
+ sys2 = xTransferFunction([[[-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 = xTransferFunction(num1, den1)
+ sys2 = xTransferFunction(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])
+
+if __name__ == "__main__":
+ unittest.main()
\ No newline at end of file
Property changes on: branches/control-0.4a/src/TestXferFcn.py
___________________________________________________________________
Added: svn:executable
+ *
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:14:03 UTC (rev 53)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:14:08 UTC (rev 54)
@@ -195,7 +195,7 @@
# Check to make sure the dimensions are OK
if (self.outputs != other.inputs):
raise ValueError, "Number of first's outputs must match number \
- of second's inputs."
+of second's inputs."
# Concatenate the various arrays
A = concatenate((
@@ -236,7 +236,7 @@
# Check to make sure the dimensions are OK
if ((self.inputs != other.outputs) or (self.outputs != other.inputs)):
raise ValueError, "State space systems don't have compatible \
- inputs/outputs for feedback."
+inputs/outputs for feedback."
# note that if there is an algebraic loop then this matrix inversion
# won't work
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:03 UTC (rev 53)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:08 UTC (rev 54)
@@ -111,11 +111,32 @@
self.den = den
Lti2.__init__(self, inputs, outputs)
+ self.truncatecoeff()
+
def __str__(self):
"""String representation of the transfer function."""
pass
+ def truncatecoeff(self):
+ """Remove extraneous zero coefficients from polynomials in numerator and
+ denominator matrices."""
+
+ data = [self.num, self.den]
+
+ for p in range(len(data)):
+ for i in range(self.outputs):
+ for j in range(self.inputs):
+ # Find the first nontrivial coefficient.
+ k = 0
+ while not data[p][i][j][k]:
+ k += 1
+
+ # Now truncate the trivial coefficients.
+ data[p][i][j] = data[p][i][j][k:]
+
+ [self.num, self.den] = data
+
def __neg__(self):
"""Negate a transfer function."""
@@ -128,16 +149,22 @@
if not isinstance(other, xTransferFunction):
other = ss2tf(other)
+ # Check that the input-output sizes are consistent.
+ if self.inputs != other.inputs:
+ raise ValueError("The two systems to be added must have the same \
+input size.")
+ if self.outputs != other.outputs:
+ raise ValueError("The two systems to be added must have the same \
+output size.")
+
# Preallocate the numerator and denominator of the sum.
num = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
den = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
for i in range(self.outputs):
for j in range(self.inputs):
- num[i][j] = sp.polyadd(
- sp.polymul(self.num[i][j], other.den[i][j]),
- sp.polymul(other.num[i][j], self.den[i][j]))
- den[i][j] = sp.polymul(self.den[i][j], other.den[i][j])
+ num[i][j], den[i][j] = addSISO(self.num[i][j], self.den[i][j],
+ other.num[i][j], other.den[i][j])
return xTransferFunction(num, den)
@@ -156,11 +183,40 @@
return other + (-self)
- def __mul__(self, sys):
+ def __mul__(self, other):
"""Multiply two transfer functions (serial connection)"""
- pass
+ # Convert the second argument to a transfer function.
+ if not isinstance(other, xTransferFunction):
+ other = ss2tf(other)
+
+ # Check that the input-output sizes are consistent.
+ if self.inputs != other.outputs:
+ raise ValueError("C = A * B: A must have the same number of \
+columns (inputs) as B has\nrows (outputs).")
+ inputs = other.inputs
+ outputs = self.outputs
+
+ # Preallocate the numerator and denominator of the sum.
+ num = [[[0] for j in range(inputs)] for i in range(outputs)]
+ den = [[[0] for j in range(inputs)] for i in range(outputs)]
+
+ # Temporary storage for the summands needed to find the (i, j)th element
+ # of the product.
+ num_summand = [[] for k in range(self.inputs)]
+ den_summand = [[] for k in range(self.inputs)]
+
+ for i in range(outputs): # Iterate through rows of product.
+ for j in range(inputs): # Iterate through columns of product.
+ for k in range(self.inputs): # Multiply & add.
+ num_summand[k] = sp.polymul(self.num[i][k], other.num[k][j])
+ den_summand[k] = sp.polymul(self.den[i][k], other.den[k][j])
+ num[i][j], den[i][j] = addSISO(num[i][j], den[i][j],
+ num_summand[k], den_summand[k])
+
+ return xTransferFunction(num, den)
+
def __rmul__(self, other):
"""Multiply two transfer functions (serial connection)"""
@@ -371,11 +427,6 @@
# Return the result as a transfer function
return TransferFunction(num, den)
-def ss2tf(sys):
- """Convert a state space object to a transfer function object."""
-
- pass
-
# Utility function to convert a transfer function polynomial to a string
# Borrowed from poly1d library
def _tfpolyToString(coeffs, var='s'):
@@ -424,3 +475,17 @@
else:
thestr = newstr
return thestr
+
+def addSISO(num1, den1, num2, den2):
+ """Return num/den = num1/den1 + num2/den2, where each numerator and
+ denominator is a list of polynomial coefficients."""
+
+ num = sp.polyadd(sp.polymul(num1, den2), sp.polymul(num2, den1))
+ den = sp.polymul(den1, den2)
+
+ return num, den
+
+def ss2tf(sys):
+ """Convert a state space object to a transfer function object."""
+
+ pass
\ No newline at end of file
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:14:12
|
Revision: 53
http://python-control.svn.sourceforge.net/python-control/?rev=53&view=rev
Author: kkchen
Date: 2011-02-08 22:14:03 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
xTransferFunction changed so that num and den are lists of lists of arrays,
instead of 3-d numpy arrays. Also implemented xTransferFunction.__add__; seems
to work so far.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:13:58 UTC (rev 52)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:14:03 UTC (rev 53)
@@ -55,69 +55,92 @@
class xTransferFunction(Lti2):
"""The TransferFunction class is derived from the Lti2 parent class. The
- main data members are 'num' and 'den', which are 3-D numpy arrays of
- MIMO numerator and denominator coefficients. For instance,
+ main data members are 'num' and 'den', which are 2-D lists of arrays
+ containing MIMO numerator and denominator coefficients. For example,
+
+ >>> num[2][5] = [1., 4., 8.]
- >>> num[2, 5] = numpy.array([1, 4, 8])
-
means that the numerator of the transfer function from the 6th input to the
- 3rd output is set to s^2 + 4s + 8.
+ 3rd output is set to s^2 + 4s + 8."""
- Note that numpy requires all polynomials in within an array to have the same
- order. Be sure to pad leading coefficients with 0 as necessary."""
-
- def __init__(self, num=sp.array([[[1.]]]), den=sp.array([[[1.]]])):
+ def __init__(self, num=1, den=1):
"""This is the constructor. The default transfer function is 1 (unit
gain direct feedthrough)."""
-
+
# Make num and den into 3-d numpy arrays, if necessary.
data = [num, den]
for i in range(len(data)):
- # Convert non-array types to array first.
- data[i] = sp.array(data[i])
-
- if data[i].ndim == 0:
- # Convert scalar to 3-d array.
- data[i] = sp.array([[[data[i]]]])
- elif data[i].ndim == 1:
- # Convert SISO transfer function polynomial to 3-d array.
- data[i] = sp.array([[data[i]]])
- elif data[i].ndim == 3:
- # We already have a MIMO system.
+ if isinstance(data[i], (int, float, long, complex)):
+ # Convert scalar to list of list of array.
+ data[i] = [[[data[i]]]]
+ elif isinstance(data[i], (list, tuple, sp.ndarray)) and \
+ isinstance(data[i][0], (int, float, long, complex)):
+ # Convert array to list of list of array.
+ data[i] = [[data[i]]]
+ elif isinstance(data[i], (list, tuple, sp.ndarray)) and \
+ isinstance(data[i][0], (list, tuple, sp.ndarray)) and \
+ isinstance(data[i][0][0], (list, tuple, sp.ndarray)) and \
+ isinstance(data[i][0][0][0], (int, float, long, complex)):
+ # We already have the right format.
pass
else:
- # If the user passed in a anything else, then it's unclear
- # what the meaning is.
+ # If the user passed in anything else, then it's unclear what
+ # the meaning is.
raise ValueError("The numerator and denominator inputs must be \
- scalars or vectors (for SISO), or 3-d arrays (for MIMO).")
+scalars or arrays (for\nSISO), or lists of lists of arrays (for MIMO).")
[num, den] = data
- inputs = num.shape[1]
- outputs = num.shape[0]
+ inputs = len(num[0])
+ outputs = len(num)
- if inputs != den.shape[1]:
+ if inputs != len(den[0]):
raise ValueError("The numerator and denominator matrices must have \
- the same column (input) size.")
- if outputs != den.shape[0]:
+the same column\n(input) size.")
+ if outputs != len(den):
raise ValueError("The numerator and denominator matrices must have \
- the same row (output) size.")
-
+the same row\n(output) size.")
+ for i in range(1, outputs):
+ if len(num[i]) != inputs:
+ raise ValueError("Each row of the numerator matrix must have \
+the same number of\nelements.")
+ if len(den[i]) != inputs:
+ raise ValueError("Each row of the denominator matrix must have \
+the same number of\nelements.")
+
self.num = num
self.den = den
Lti2.__init__(self, inputs, outputs)
- def __str__(self): pass
+ def __str__(self):
+ """String representation of the transfer function."""
+ pass
+
def __neg__(self):
"""Negate a transfer function."""
return xTransferFunction(-self.num, self.den)
- def __add__(self, sys):
+ def __add__(self, other):
"""Add two transfer functions (parallel connection)."""
-
- pass
-
+
+ # Convert the second argument to a transfer function.
+ if not isinstance(other, xTransferFunction):
+ other = ss2tf(other)
+
+ # Preallocate the numerator and denominator of the sum.
+ num = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
+ den = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
+
+ for i in range(self.outputs):
+ for j in range(self.inputs):
+ num[i][j] = sp.polyadd(
+ sp.polymul(self.num[i][j], other.den[i][j]),
+ sp.polymul(other.num[i][j], self.den[i][j]))
+ den[i][j] = sp.polymul(self.den[i][j], other.den[i][j])
+
+ return xTransferFunction(num, den)
+
def __radd__(self, other):
"""Add two transfer functions (parallel connection)"""
@@ -348,24 +371,11 @@
# Return the result as a transfer function
return TransferFunction(num, den)
-# Function to create a transfer function from another type
-def convertToTransferFunction(sys):
- """Convert a system to a transfer fuction (if needed)"""
- if (isinstance(sys, TransferFunction)):
- # Already a transfer function; just return it
- return sys
+def ss2tf(sys):
+ """Convert a state space object to a transfer function object."""
- elif (isinstance(sys, statesp.StateSpace)):
- # State space system, convert using signal.lti
- return TransferFunction(sys.A, sys.B, sys.C, sys.D)
+ pass
- elif (isinstance(sys, (int, long, float, complex))):
- # Convert a number into a transfer function
- return TransferFunction(sys, 1)
-
- else:
- raise TypeError("can't convert given type to TransferFunction")
-
# Utility function to convert a transfer function polynomial to a string
# Borrowed from poly1d library
def _tfpolyToString(coeffs, var='s'):
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:14:04
|
Revision: 52
http://python-control.svn.sourceforge.net/python-control/?rev=52&view=rev
Author: kkchen
Date: 2011-02-08 22:13:58 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Changed input processing in xferfcn.py; wrote skeleton for xTransferFunction.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/xferfcn.py
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:13:53 UTC (rev 51)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:13:58 UTC (rev 52)
@@ -73,30 +73,25 @@
# Make num and den into 3-d numpy arrays, if necessary.
data = [num, den]
for i in range(len(data)):
- if isinstance(data[i], (int, long, float, complex)):
+ # Convert non-array types to array first.
+ data[i] = sp.array(data[i])
+
+ if data[i].ndim == 0:
# Convert scalar to 3-d array.
data[i] = sp.array([[[data[i]]]])
- elif isinstance(data[i], sp.ndarray):
- if data[i].ndim == 0:
- # Convert scalar to 3-d array.
- data[i] = sp.array([[[data[i]]]])
- elif data[i].ndim == 1:
- # Convert SISO transfer function polynomial to 3-d array.
- data[i] = sp.array([[data[i]]])
- elif data[i].ndim == 3:
- # We already have a MIMO system.
- pass
- else:
- # If the user passed in a anything else, then it's unclear
- # what the meaning is.
- raise ValueError("The numerator and denominator inputs \
- must be scalars or vectors (for SISO), or 3-d arrays (for \
- MIMO).")
+ elif data[i].ndim == 1:
+ # Convert SISO transfer function polynomial to 3-d array.
+ data[i] = sp.array([[data[i]]])
+ elif data[i].ndim == 3:
+ # We already have a MIMO system.
+ pass
+ else:
+ # If the user passed in a anything else, then it's unclear
+ # what the meaning is.
+ raise ValueError("The numerator and denominator inputs must be \
+ scalars or vectors (for SISO), or 3-d arrays (for MIMO).")
[num, den] = data
- print num
- print den
-
inputs = num.shape[1]
outputs = num.shape[0]
@@ -111,9 +106,80 @@
self.den = den
Lti2.__init__(self, inputs, outputs)
- print self.inputs
- print self.outputs
+ def __str__(self): pass
+
+ def __neg__(self):
+ """Negate a transfer function."""
+
+ return xTransferFunction(-self.num, self.den)
+
+ def __add__(self, sys):
+ """Add two transfer functions (parallel connection)."""
+
+ pass
+
+ def __radd__(self, other):
+ """Add two transfer functions (parallel connection)"""
+
+ return self + other;
+
+ def __sub__(self, other):
+ """Subtract two transfer functions."""
+
+ return self + (-other)
+
+ def __rsub__(self, other):
+ """Subtract two transfer functions."""
+
+ return other + (-self)
+ def __mul__(self, sys):
+ """Multiply two transfer functions (serial connection)"""
+
+ pass
+
+ def __rmul__(self, other):
+ """Multiply two transfer functions (serial connection)"""
+
+ return self * other
+
+ def __div__(self, sys):
+ """Divide two transfer functions"""
+
+ pass
+
+ def __rdiv__(self, sys):
+ """Reverse divide two transfer functions"""
+
+ pass
+
+ def evalfr(self, freq):
+ """Evaluate a transfer function at a single frequency"""
+
+ # return sp.polyval(self.num, freq*1j) / sp.polyval(self.den, freq*1j)
+ pass
+
+ # Method for generating the frequency response of the system
+ def freqresp(self, omega):
+ """Evaluate a transfer function at a list of frequencies"""
+
+ pass
+
+ def poles(self):
+ """Compute poles of a transfer function."""
+
+ pass
+
+ def zeros(self):
+ """Compute zeros of a transfer function."""
+
+ pass
+
+ def feedback(sys1, sys2, sign=-1):
+ """Feedback interconnection between two transfer functions"""
+
+ pass
+
class TransferFunction(signal.lti):
"""The TransferFunction class is used to represent linear
input/output systems via its transfer function.
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:14:00
|
Revision: 51
http://python-control.svn.sourceforge.net/python-control/?rev=51&view=rev
Author: kkchen
Date: 2011-02-08 22:13:53 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
*Class structure change*
Added Lti2 superclass in lti2.py.
Started the transfer function subclass in xferfcn.py, currently implemented as xTrasferFunction.
Edits to the StateSpace class in statesp.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/statesp.py
branches/control-0.4a/src/xferfcn.py
Added Paths:
-----------
branches/control-0.4a/src/lti2.py
Added: branches/control-0.4a/src/lti2.py
===================================================================
--- branches/control-0.4a/src/lti2.py (rev 0)
+++ branches/control-0.4a/src/lti2.py 2011-02-08 22:13:53 UTC (rev 51)
@@ -0,0 +1,9 @@
+class Lti2:
+ """The Lti2 is a parent class to the StateSpace and TransferFunction child
+ classes. It only contains the number of inputs and outputs, but this can be
+ expanded in the future."""
+
+ def __init__(self, inputs=1, outputs=1):
+ # Data members common to StateSpace and TransferFunction.
+ self.inputs = inputs
+ self.outputs = outputs
\ No newline at end of file
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:13:47 UTC (rev 50)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:13:53 UTC (rev 51)
@@ -42,46 +42,25 @@
# $Id: statesp.py 21 2010-06-06 17:29:42Z murrayrm $
import scipy as sp
-import scipy.signal as signal
-import xferfcn
from scipy import concatenate, zeros
+import xferfcn
+from lti2 import Lti2
-class Lti2:
- """The Lti2 is a parent class to the StateSpace and TransferFunction child
- classes. It only contains the number of inputs and outputs, but this can be
- expanded in the future."""
-
- def __init__(self, inputs=1, outputs=1):
- # Data members common to StateSpace and TransferFunction.
- self.inputs = inputs
- self.outputs = outputs
-
class StateSpace(Lti2):
"""The StateSpace class is used throughout the python-control library to
represent systems in state space form. This class is derived from the Lti2
base class."""
- def __init__(self, A, B, C, D):
- # Here we're going to convert inputs to matrices, if the user gave a non
- # 2-D array or matrix type.
- matrices = [A, B, C, D]
+ def __init__(self, A=0, B=0, C=0, D=1):
+ """StateSpace constructor. The default constructor is the unit gain
+ direct feedthrough system."""
+
+ # Here we're going to convert inputs to matrices, if the user gave a
+ # non-matrix type.
+ matrices = [A, B, C, D]
for i in range(len(matrices)):
- if (isinstance(matrices[i], (int, long, float, complex))):
- # Convert scalars to matrices, if necessary.
- matrices[i] = sp.matrix(matrices[i])
- elif isinstance(matrices[i], sp.ndarray):
- # Convert 0- or 1-D arrays to matrices, if necessary.
- if len(matrices[i].shape) < 2:
- matrices[i] = sp.matrix(matrices[i])
- elif len(matrices[i].shape) == 2:
- # If we're already a 2-D array or a matrix, then perfect!
- pass
- else:
- raise ValueError("A, B, C, and D cannot have > 2 \
- dimensions.")
- else:
- # If the user gave us a non-numeric type.
- raise ValueError("A, B, C, and D must be arrays or matrices.")
+ # Convert to matrix first, if necessary.
+ matrices[i] = sp.matrix(matrices[i])
[A, B, C, D] = matrices
self.A = A
Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py 2011-02-08 22:13:47 UTC (rev 50)
+++ branches/control-0.4a/src/xferfcn.py 2011-02-08 22:13:53 UTC (rev 51)
@@ -51,7 +51,69 @@
import scipy.signal as signal
import bdalg as bd
import statesp
+from lti2 import Lti2
+class xTransferFunction(Lti2):
+ """The TransferFunction class is derived from the Lti2 parent class. The
+ main data members are 'num' and 'den', which are 3-D numpy arrays of
+ MIMO numerator and denominator coefficients. For instance,
+
+ >>> num[2, 5] = numpy.array([1, 4, 8])
+
+ means that the numerator of the transfer function from the 6th input to the
+ 3rd output is set to s^2 + 4s + 8.
+
+ Note that numpy requires all polynomials in within an array to have the same
+ order. Be sure to pad leading coefficients with 0 as necessary."""
+
+ def __init__(self, num=sp.array([[[1.]]]), den=sp.array([[[1.]]])):
+ """This is the constructor. The default transfer function is 1 (unit
+ gain direct feedthrough)."""
+
+ # Make num and den into 3-d numpy arrays, if necessary.
+ data = [num, den]
+ for i in range(len(data)):
+ if isinstance(data[i], (int, long, float, complex)):
+ # Convert scalar to 3-d array.
+ data[i] = sp.array([[[data[i]]]])
+ elif isinstance(data[i], sp.ndarray):
+ if data[i].ndim == 0:
+ # Convert scalar to 3-d array.
+ data[i] = sp.array([[[data[i]]]])
+ elif data[i].ndim == 1:
+ # Convert SISO transfer function polynomial to 3-d array.
+ data[i] = sp.array([[data[i]]])
+ elif data[i].ndim == 3:
+ # We already have a MIMO system.
+ pass
+ else:
+ # If the user passed in a anything else, then it's unclear
+ # what the meaning is.
+ raise ValueError("The numerator and denominator inputs \
+ must be scalars or vectors (for SISO), or 3-d arrays (for \
+ MIMO).")
+ [num, den] = data
+
+ print num
+ print den
+
+ inputs = num.shape[1]
+ outputs = num.shape[0]
+
+ if inputs != den.shape[1]:
+ raise ValueError("The numerator and denominator matrices must have \
+ the same column (input) size.")
+ if outputs != den.shape[0]:
+ raise ValueError("The numerator and denominator matrices must have \
+ the same row (output) size.")
+
+ self.num = num
+ self.den = den
+ Lti2.__init__(self, inputs, outputs)
+
+ print self.inputs
+ print self.outputs
+
class TransferFunction(signal.lti):
"""The TransferFunction class is used to represent linear
input/output systems via its transfer function.
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:54
|
Revision: 50
http://python-control.svn.sourceforge.net/python-control/?rev=50&view=rev
Author: kkchen
Date: 2011-02-08 22:13:47 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Created Lti2 parent class and made StateSpace its child, in statesp.py.
Added .DS_Store to .hgignore.
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:13:44 UTC (rev 49)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:13:47 UTC (rev 50)
@@ -46,19 +46,21 @@
import xferfcn
from scipy import concatenate, zeros
-#
-# StateSpace class
-#
-# The StateSpace class is used throughout the control systems library to
-# represent systems in state space form. This class is derived from
-# the ltisys class defined in the scipy.signal package, allowing many
-# of the functions that already existing in that package to be used
-# directly.
-#
-class StateSpace:
- """The StateSpace class is used to represent linear input/output systems.
- """
+class Lti2:
+ """The Lti2 is a parent class to the StateSpace and TransferFunction child
+ classes. It only contains the number of inputs and outputs, but this can be
+ expanded in the future."""
+
+ def __init__(self, inputs=1, outputs=1):
+ # Data members common to StateSpace and TransferFunction.
+ self.inputs = inputs
+ self.outputs = outputs
+class StateSpace(Lti2):
+ """The StateSpace class is used throughout the python-control library to
+ represent systems in state space form. This class is derived from the Lti2
+ base class."""
+
def __init__(self, A, B, C, D):
# Here we're going to convert inputs to matrices, if the user gave a non
# 2-D array or matrix type.
@@ -88,8 +90,7 @@
self.D = D
self.states = A.shape[0]
- self.inputs = B.shape[1]
- self.outputs = C.shape[0]
+ Lti2.__init__(self, B.shape[1], C.shape[0])
# Check that the matrix sizes are consistent.
if self.states != A.shape[1]:
@@ -103,8 +104,9 @@
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):
+ """Style to use for printing."""
+
str = "A = " + self.A.__str__() + "\n\n"
str += "B = " + self.B.__str__() + "\n\n"
str += "C = " + self.C.__str__() + "\n\n"
@@ -113,7 +115,8 @@
# Method for generating the frequency response of the system
def freqresp(self, omega=None):
- """Compute the response of a system to a list of frequencies"""
+ """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
@@ -154,12 +157,14 @@
# Negation of a system
def __neg__(self):
- """Negate a state space system"""
+ """Negate a state space system."""
+
return StateSpace(self.A, self.B, -self.C, -self.D)
# Addition of two transfer functions (parallel interconnection)
def __add__(self, other):
- """Add two state space systems"""
+ """Add two state space systems."""
+
# Check for a couple of special cases
if (isinstance(other, (int, long, float, complex))):
# Just adding a scalar; put it in the D matrix
@@ -186,17 +191,20 @@
# Reverse addition - just switch the arguments
def __radd__(self, other):
- """Add two state space systems"""
+ """Add two state space systems."""
+
return self.__add__(other)
# Subtraction of two transfer functions (parallel interconnection)
def __sub__(self, other):
- """Subtract two state space systems"""
+ """Subtract two state space systems."""
+
return __add__(self, other.__neg__())
# Multiplication of two transfer functions (series interconnection)
def __mul__(self, other):
- """Serial interconnection between two state space systems"""
+ """Serial interconnection between two state space systems."""
+
# Check for a couple of special cases
if (isinstance(other, (int, long, float, complex))):
# Just multiplying by a scalar; change the output
@@ -225,6 +233,7 @@
# Just need to convert LH argument to a state space object
def __rmul__(self, other):
"""Serial interconnection between two state space systems"""
+
# Check for a couple of special cases
if (isinstance(other, (int, long, float, complex))):
# Just multiplying by a scalar; change the input
@@ -238,7 +247,8 @@
# Feedback around a state space system
def feedback(self, other, sign=-1):
- """Feedback interconnection between two state space systems"""
+ """Feedback interconnection between two state space systems."""
+
# Check for special cases
if (isinstance(other, (int, long, float, complex))):
# Scalar feedback, create state space system that is this case
@@ -281,6 +291,7 @@
#
def convertToStateSpace(sys, inputs=1, outputs=1):
"""Convert a system to state space form (if needed)"""
+
if isinstance(sys, StateSpace):
# Already a state space system; just return it
return sys
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:50
|
Revision: 49
http://python-control.svn.sourceforge.net/python-control/?rev=49&view=rev
Author: kkchen
Date: 2011-02-08 22:13:44 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added interface for balred and modred in modelsimp.py. Also changed * to dot in gram function.
Steven Brunton <sbr...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/modelsimp.py
branches/control-0.4a/src/statefbk.py
Modified: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py 2011-02-08 22:13:40 UTC (rev 48)
+++ branches/control-0.4a/src/modelsimp.py 2011-02-08 22:13:44 UTC (rev 49)
@@ -77,6 +77,112 @@
# Return the Hankel singular values
return hsv
+def modred(sys,ELIM,method):
+ """Model reduction of sys by eliminating the states in ELIM using a given method
+
+ Usage
+ =====
+ rsys = modred(sys,ELIM,method)
+
+ Inputs
+ ======
+ sys : original system to reduce
+ ELIM : vector of states to eliminate
+ method : method of removing states in ELIM (truncate or matchdc)
+
+ Outputs
+ =======
+ rsys : a reduced order model
+
+ """
+
+ #Check for ss system object, need a utility for this?
+
+ #TODO: Check for continous or discrete, only continuous supported right now
+ # if isCont():
+ # dico = 'C'
+ # elif isDisc():
+ # dico = 'D'
+ # else:
+ dico = 'C'
+
+ #TODO: Check system is stable, perhaps a utility in ctrlutil.py
+ # or a method of the StateSpace class?
+ D,V = np.linalg.eig(sys.A)
+ for e in D:
+ if e.real >= 0:
+ raise ValueError, "Oops, the system is unstable!"
+
+ if method=='matchdc':
+ print "matchdc"
+ elif method=='truncate':
+ print "truncate"
+ else:
+ raise ValueError, "Oops, method is not supported!"
+
+ #Compute rhs using the Slycot routine XXXXXX
+ #make sure Slycot is installed
+ #try:
+ # from slycot import XXXXXX
+ #except ImportError:
+ # raise ControlSlycot("can't find slycot module 'XXXXXX'")
+ rsys = 0.
+ return rsys
+
+def balred(sys,orders,elimination,method):
+ """Balanced reduced order model of sys of a given order. States are eliminated based on Hankel singular value.
+
+ Usage
+ =====
+ rsys = balred(sys,order,elimination,method)
+
+ Inputs
+ ======
+ sys : original system to reduce
+ orders : desired order of reduced order model (if a vector, returns a vector of systems)
+ elimination : if elimination is specified, use 'method'
+ method : method of removing states (truncate or matchdc)
+
+ Outputs
+ =======
+ rsys : a reduced order model
+
+ """
+
+ #Check for ss system object, need a utility for this?
+
+ #TODO: Check for continous or discrete, only continuous supported right now
+ # if isCont():
+ # dico = 'C'
+ # elif isDisc():
+ # dico = 'D'
+ # else:
+ dico = 'C'
+
+ #TODO: Check system is stable, perhaps a utility in ctrlutil.py
+ # or a method of the StateSpace class?
+ D,V = np.linalg.eig(sys.A)
+ for e in D:
+ if e.real >= 0:
+ raise ValueError, "Oops, the system is unstable!"
+
+ if method=='matchdc':
+ print "matchdc"
+ elif method=='truncate':
+ print "truncate"
+ else:
+ raise ValueError, "Oops, method is not supported!"
+
+ #Compute rhs using the Slycot routine XXXXXX
+ #make sure Slycot is installed
+ #try:
+ # from slycot import XXXXXX
+ #except ImportError:
+ # raise ControlSlycot("can't find slycot module 'XXXXXX'")
+ rsys = 0.
+ return rsys
+
+
def era(YY,m,n,nin,nout,r):
"""Calculate an ERA model of order r based on the impulse-response data YY
Modified: branches/control-0.4a/src/statefbk.py
===================================================================
--- branches/control-0.4a/src/statefbk.py 2011-02-08 22:13:40 UTC (rev 48)
+++ branches/control-0.4a/src/statefbk.py 2011-02-08 22:13:44 UTC (rev 49)
@@ -267,11 +267,11 @@
if type=='c':
print "controllable"
trana = 'T'
- C = -sys.B*sys.B.transpose()
+ C = -np.dot(sys.B,sys.B.transpose())
elif type=='o':
print "observable"
trana = 'N'
- C = -sys.C.transpose()*sys.C
+ C = -np.dot(sys.C.transpose(),sys.C)
else:
raise ValueError, "Oops, neither observable, nor controllable!"
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:46
|
Revision: 48
http://python-control.svn.sourceforge.net/python-control/?rev=48&view=rev
Author: kkchen
Date: 2011-02-08 22:13:40 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added pole() to matlab.py.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/matlab.py
Modified: branches/control-0.4a/src/matlab.py
===================================================================
--- branches/control-0.4a/src/matlab.py 2011-02-08 22:13:36 UTC (rev 47)
+++ branches/control-0.4a/src/matlab.py 2011-02-08 22:13:40 UTC (rev 48)
@@ -129,7 +129,7 @@
dcgain - steady-state (D.C.) gain
lti/bandwidth - system bandwidth
lti/norm - h2 and Hinfinity norms of LTI models
- lti/pole - system poles
+* lti/pole - system poles
lti/zero - system (transmission) zeros
lti/order - model order (number of states)
* pzmap - pole-zero map
@@ -316,6 +316,9 @@
return rss_generate(states, inputs, outputs, 'd')
+def pole(sys):
+ return sys.poles()
+
# Frequency response is handled by the system object
def freqresp(H, omega):
"""Return the frequency response for an object H at frequency omega"""
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:42
|
Revision: 47
http://python-control.svn.sourceforge.net/python-control/?rev=47&view=rev
Author: kkchen
Date: 2011-02-08 22:13:36 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added initial unittests for modred and balred to TestModelsimp.py and debugged hsvd.
Lauren Padilla <lpa...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestModelsimp.py
branches/control-0.4a/src/modelsimp.py
Modified: branches/control-0.4a/src/TestModelsimp.py
===================================================================
--- branches/control-0.4a/src/TestModelsimp.py 2011-02-08 22:13:30 UTC (rev 46)
+++ branches/control-0.4a/src/TestModelsimp.py 2011-02-08 22:13:36 UTC (rev 47)
@@ -24,5 +24,93 @@
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,np.matrix('3, 4'),'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(sys.A, Artrue)
+ np.testing.assert_array_almost_equal(sys.B, Brtrue)
+ np.testing.assert_array_almost_equal(sys.C, Crtrue)
+ np.testing.assert_array_almost_equal(sys.D, Drtrue)
+
+ 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,np.matrix('3, 4'),'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(sys.A, Artrue)
+ np.testing.assert_array_almost_equal(sys.B, Brtrue)
+ np.testing.assert_array_almost_equal(sys.C, Crtrue)
+ np.testing.assert_array_almost_equal(sys.D, Drtrue)
+
+ def testBalredMatchDC(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 = 1
+ rsys = balred(sys,orders,'elimination','matchdc')
+ Artrue = np.matrix('-0.566')
+ Brtrue = np.matrix('-0.414')
+ Crtrue = np.matrix('-0.5728')
+ Drtrue = np.matrix('0.1145')
+ np.testing.assert_array_almost_equal(sys.A, Artrue)
+ np.testing.assert_array_almost_equal(sys.B, Brtrue)
+ np.testing.assert_array_almost_equal(sys.C, Crtrue)
+ np.testing.assert_array_almost_equal(sys.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,'elimination','truncate')
+ Artrue = np.matrix('-1.95, 0.6203; 2.314, -0.8432')
+ Brtrue = np.matrix('0.7221; -0.6306')
+ Crtrue = np.matrix('1.132, -0.2667')
+ Drtrue = np.matrix('0.')
+ np.testing.assert_array_almost_equal(sys.A, Artrue)
+ np.testing.assert_array_almost_equal(sys.B, Brtrue)
+ np.testing.assert_array_almost_equal(sys.C, Crtrue)
+ np.testing.assert_array_almost_equal(sys.D, Drtrue)
+
+
+
if __name__ == '__main__':
unittest.main()
Modified: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py 2011-02-08 22:13:30 UTC (rev 46)
+++ branches/control-0.4a/src/modelsimp.py 2011-02-08 22:13:36 UTC (rev 47)
@@ -43,6 +43,7 @@
import numpy as np
import ctrlutil
from control.exception import *
+from statefbk 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
@@ -69,7 +70,7 @@
Wo = gram(sys,'o')
WoWc = np.dot(Wo, Wc)
- w, v = LA.eig(WoWc)
+ w, v = np.linalg.eig(WoWc)
hsv = np.sqrt(w)
hsv = np.matrix(hsv)
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:36
|
Revision: 46
http://python-control.svn.sourceforge.net/python-control/?rev=46&view=rev
Author: kkchen
Date: 2011-02-08 22:13:30 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Moved rss, drss to matlab.py.
Changes made in {matlab,statesp,TestStateSp}.py as necessary.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestStateSp.py
branches/control-0.4a/src/matlab.py
branches/control-0.4a/src/statesp.py
Modified: branches/control-0.4a/src/TestStateSp.py
===================================================================
--- branches/control-0.4a/src/TestStateSp.py 2011-02-08 22:13:26 UTC (rev 45)
+++ branches/control-0.4a/src/TestStateSp.py 2011-02-08 22:13:30 UTC (rev 46)
@@ -1,7 +1,7 @@
#!/usr/bin/env python
import numpy as np
-import statesp as ss
+import matlab
import unittest
class TestRss(unittest.TestCase):
@@ -22,7 +22,7 @@
for states in range(1, self.maxStates):
for inputs in range(1, self.maxIO):
for outputs in range(1, self.maxIO):
- sys = ss.rss(states, inputs, outputs)
+ sys = matlab.rss(states, inputs, outputs)
self.assertEqual(sys.states, states)
self.assertEqual(sys.inputs, inputs)
self.assertEqual(sys.outputs, outputs)
@@ -33,7 +33,7 @@
for states in range(1, self.maxStates):
for inputs in range(1, self.maxIO):
for outputs in range(1, self.maxIO):
- sys = ss.rss(states, inputs, outputs)
+ sys = matlab.rss(states, inputs, outputs)
p = sys.poles()
for z in p:
self.assertTrue(z.real < 0)
@@ -56,7 +56,7 @@
for states in range(1, self.maxStates):
for inputs in range(1, self.maxIO):
for outputs in range(1, self.maxIO):
- sys = ss.drss(states, inputs, outputs)
+ sys = matlab.drss(states, inputs, outputs)
self.assertEqual(sys.states, states)
self.assertEqual(sys.inputs, inputs)
self.assertEqual(sys.outputs, outputs)
@@ -67,7 +67,7 @@
for states in range(1, self.maxStates):
for inputs in range(1, self.maxIO):
for outputs in range(1, self.maxIO):
- sys = ss.drss(states, inputs, outputs)
+ sys = matlab.drss(states, inputs, outputs)
p = sys.poles()
for z in p:
self.assertTrue(abs(z) < 1)
Modified: branches/control-0.4a/src/matlab.py
===================================================================
--- branches/control-0.4a/src/matlab.py 2011-02-08 22:13:26 UTC (rev 45)
+++ branches/control-0.4a/src/matlab.py 2011-02-08 22:13:30 UTC (rev 46)
@@ -60,7 +60,7 @@
# Control system library
import ctrlutil
import freqplot
-from statesp import StateSpace
+from statesp import StateSpace, rss_generate
from xferfcn import TransferFunction
from exception import *
@@ -187,8 +187,8 @@
augstate - augment output by appending states
State-space (SS) models
- rss - random stable continuous-time state-space models
- drss - random stable discrete-time state-space models
+* rss - random stable continuous-time state-space models
+* drss - random stable discrete-time state-space models
ss2ss - state coordinate transformation
canon - canonical forms of state-space models
* ctrb - controllability matrix
@@ -306,6 +306,16 @@
else:
raise ValueError, "Needs 1 or 2 arguments."
+def rss(states=1, inputs=1, outputs=1):
+ """Create a stable continuous random state space object."""
+
+ return rss_generate(states, inputs, outputs, 'c')
+
+def drss(states=1, inputs=1, outputs=1):
+ """Create a stable discrete random state space object."""
+
+ return rss_generate(states, inputs, outputs, 'd')
+
# Frequency response is handled by the system object
def freqresp(H, omega):
"""Return the frequency response for an object H at frequency omega"""
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:13:26 UTC (rev 45)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:13:30 UTC (rev 46)
@@ -295,17 +295,7 @@
else:
raise TypeError("can't convert given type to StateSpace system")
-
-def rss(states=1, inputs=1, outputs=1):
- """Create a stable continuous random state space object."""
- return rss_generate(states, inputs, outputs, 'c')
-
-def drss(states=1, inputs=1, outputs=1):
- """Create a stable discrete random state space object."""
-
- return rss_generate(states, inputs, outputs, 'd')
-
def rss_generate(states, inputs, outputs, type):
"""This does the actual random state space generation expected from rss and
drss. type is 'c' for continuous systems and 'd' for discrete systems."""
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:34
|
Revision: 45
http://python-control.svn.sourceforge.net/python-control/?rev=45&view=rev
Author: kkchen
Date: 2011-02-08 22:13:26 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Minor changes to statesp.py related to scipy.signal.lti overhaul.
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:13:22 UTC (rev 44)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:13:26 UTC (rev 45)
@@ -281,16 +281,17 @@
#
def convertToStateSpace(sys, inputs=1, outputs=1):
"""Convert a system to state space form (if needed)"""
- if (isinstance(sys, StateSpace) or
- isinstance(sys, xferfcn.TransferFunction)):
+ if isinstance(sys, StateSpace):
# Already a state space system; just return it
return sys
-
+ elif isinstance(sys, xferfcn.TransferFunction):
+ pass # TODO: convert SS to TF
elif (isinstance(sys, (int, long, float, complex))):
# Generate a simple state space system of the desired dimension
# The following Doesn't work due to inconsistencies in ltisys:
# return StateSpace([[]], [[]], [[]], sp.eye(outputs, inputs))
- return StateSpace(-1, 0, 0, sp.eye(outputs, inputs))
+ return StateSpace(-1, zeros((1, inputs)), zeros((outputs, 1)),
+ sp.eye(outputs, inputs))
else:
raise TypeError("can't convert given type to StateSpace system")
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:28
|
Revision: 44
http://python-control.svn.sourceforge.net/python-control/?rev=44&view=rev
Author: kkchen
Date: 2011-02-08 22:13:22 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added tests for valid inputs to StateSpace in statesp.py.
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:13:17 UTC (rev 43)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:13:22 UTC (rev 44)
@@ -60,6 +60,28 @@
"""
def __init__(self, A, B, C, D):
+ # Here we're going to convert inputs to matrices, if the user gave a non
+ # 2-D array or matrix type.
+ matrices = [A, B, C, D]
+ for i in range(len(matrices)):
+ if (isinstance(matrices[i], (int, long, float, complex))):
+ # Convert scalars to matrices, if necessary.
+ matrices[i] = sp.matrix(matrices[i])
+ elif isinstance(matrices[i], sp.ndarray):
+ # Convert 0- or 1-D arrays to matrices, if necessary.
+ if len(matrices[i].shape) < 2:
+ matrices[i] = sp.matrix(matrices[i])
+ elif len(matrices[i].shape) == 2:
+ # If we're already a 2-D array or a matrix, then perfect!
+ pass
+ else:
+ raise ValueError("A, B, C, and D cannot have > 2 \
+ dimensions.")
+ else:
+ # If the user gave us a non-numeric type.
+ raise ValueError("A, B, C, and D must be arrays or matrices.")
+ [A, B, C, D] = matrices
+
self.A = A
self.B = B
self.C = C
@@ -69,11 +91,6 @@
self.inputs = B.shape[1]
self.outputs = C.shape[0]
- # Check that the inputs are arrays or matrices.
- if not (isinstance(A, sp.ndarray) and isinstance(B, sp.ndarray) and
- isinstance(C, sp.ndarray) and isinstance(D, sp.ndarray)):
- raise ValueError("A, B, C, and D must be arrays or matrices.")
-
# Check that the matrix sizes are consistent.
if self.states != A.shape[1]:
raise ValueError("A must be square.")
@@ -389,4 +406,4 @@
C = C * Cmask
D = D * Dmask
- return StateSpace(A, B, C, D)
\ No newline at end of file
+ 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:13:23
|
Revision: 43
http://python-control.svn.sourceforge.net/python-control/?rev=43&view=rev
Author: kkchen
Date: 2011-02-08 22:13:17 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Added drss to statesp.py; added drss unit tests to TestStateSp.py.
Also changed the permissions of Test*.py to 755.
Kevin K. Chen <kk...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/TestStateSp.py
branches/control-0.4a/src/statesp.py
Modified: branches/control-0.4a/src/TestStateSp.py
===================================================================
--- branches/control-0.4a/src/TestStateSp.py 2011-02-08 22:13:13 UTC (rev 42)
+++ branches/control-0.4a/src/TestStateSp.py 2011-02-08 22:13:17 UTC (rev 43)
@@ -5,34 +5,72 @@
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)
-
+ """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 = 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, self.maxStates):
+ for inputs in range(1, self.maxIO):
+ for outputs in range(1, self.maxIO):
+ sys = ss.rss(states, inputs, outputs)
+ p = sys.poles()
+ 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 = ss.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 = ss.drss(states, inputs, outputs)
+ p = sys.poles()
+ for z in p:
+ self.assertTrue(abs(z) < 1)
+
if __name__ == "__main__":
- unittest.main()
\ No newline at end of file
+ unittest.main()
\ No newline at end of file
Modified: branches/control-0.4a/src/statesp.py
===================================================================
--- branches/control-0.4a/src/statesp.py 2011-02-08 22:13:13 UTC (rev 42)
+++ branches/control-0.4a/src/statesp.py 2011-02-08 22:13:17 UTC (rev 43)
@@ -58,7 +58,7 @@
class StateSpace:
"""The StateSpace class is used to represent linear input/output systems.
"""
- # Initialization
+
def __init__(self, A, B, C, D):
self.A = A
self.B = B
@@ -68,7 +68,12 @@
self.states = A.shape[0]
self.inputs = B.shape[1]
self.outputs = C.shape[0]
-
+
+ # Check that the inputs are arrays or matrices.
+ if not (isinstance(A, sp.ndarray) and isinstance(B, sp.ndarray) and
+ isinstance(C, sp.ndarray) and isinstance(D, sp.ndarray)):
+ raise ValueError("A, B, C, and D must be arrays or matrices.")
+
# Check that the matrix sizes are consistent.
if self.states != A.shape[1]:
raise ValueError("A must be square.")
@@ -224,9 +229,11 @@
# Check to make sure the dimensions are OK
if ((self.inputs != other.outputs) or (self.outputs != other.inputs)):
- raise ValueError, "State space systems don't have compatible inputs/outputs for feedback."
+ raise ValueError, "State space systems don't have compatible \
+ inputs/outputs for feedback."
- # note that if there is an algebraic loop then this matrix inversion won't work
+ # note that if there is an algebraic loop then this matrix inversion
+ # won't work
# (I-D1 D2) or (I-D2 D1) will be singular
# the easiest way to get this is to have D1 = I, D2 = I
#! TODO: trap this error and report algebraic loop
@@ -237,9 +244,11 @@
E12 = inv(eye(self.inputs)+sign*other.D*self.D)
A = concatenate((
- concatenate(( self.A-sign*self.B*E12*other.D*self.C, -sign*self.B*E12*other.C ),axis=1),
- concatenate(( other.B*E21*self.C, other.A-sign*other.B*E21*self.D*other.C ),axis=1),
- ),axis=0)
+ concatenate(( self.A-sign*self.B*E12*other.D*self.C,
+ -sign*self.B*E12*other.C ),axis=1),
+ concatenate(( other.B*E21*self.C,
+ other.A-sign*other.B*E21*self.D*other.C ),axis=1),),
+ axis=0)
B = concatenate( (self.B*E12, other.B*E21*self.D), axis=0 )
C = concatenate( (E21*self.C, -sign*E21*self.D*other.C), axis=1 )
D = E21*self.D
@@ -270,17 +279,43 @@
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."""
+ """Create a stable continuous random state space object."""
+ return rss_generate(states, inputs, outputs, 'c')
+
+def drss(states=1, inputs=1, outputs=1):
+ """Create a stable discrete random state space object."""
+
+ return rss_generate(states, inputs, outputs, 'd')
+
+def rss_generate(states, inputs, outputs, type):
+ """This does the actual random state space generation expected from rss and
+ drss. type is 'c' for continuous systems and 'd' for discrete systems."""
+
import numpy
from numpy.random import rand, randn
+
+ # Probability of repeating a previous root.
+ pRepeat = 0.05
+ # Probability of choosing a real root. Note that when choosing a complex
+ # root, the conjugate gets chosen as well. So the expected proportion of
+ # real roots is pReal / (pReal + 2 * (1 - pReal)).
+ pReal = 0.6
+ # Probability that an element in B or C will not be masked out.
+ pBCmask = 0.8
+ # Probability that an element in D will not be masked out.
+ pDmask = 0.3
+ # Probability that D = 0.
+ pDzero = 0.5
# 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.
+
+ while i < states:
+ if rand() < pRepeat and i != 0 and i != states - 1:
+ # Small chance of copying poles, if we're not at the first or last
+ # element.
if poles[i-1].imag == 0:
# Copy previous real pole.
poles[i] = poles[i-1]
@@ -289,19 +324,24 @@
# 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
+ elif rand() < pReal or i == states - 1:
+ # No-oscillation pole.
+ if type == 'c':
+ poles[i] = -sp.exp(randn()) + 0.j
+ elif type == 'd':
+ poles[i] = 2. * rand() - 1.
i += 1
else:
- # Complex conjugate pair of poles.
- poles[i] = complex(-sp.exp(randn()), sp.exp(randn()))
+ # Complex conjugate pair of oscillating poles.
+ if type == 'c':
+ poles[i] = complex(-sp.exp(randn()), 3. * sp.exp(randn()))
+ elif type == 'd':
+ mag = rand()
+ phase = 2. * numpy.pi * rand()
+ poles[i] = complex(mag * numpy.cos(phase),
+ mag * numpy.sin(phase))
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))
@@ -315,7 +355,6 @@
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)
@@ -333,18 +372,21 @@
# 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.
+ Bmask = rand(states, inputs) < pBCmask
+ if sp.any(Bmask): # 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.
+ Cmask = rand(outputs, states) < pBCmask
+ if sp.any(Cmask): # Retry if we get all zeros.
break
- D_mask = rand(outputs, inputs) < 0.3
+ if rand() < pDzero:
+ Dmask = numpy.zeros((outputs, inputs))
+ else:
+ Dmask = rand(outputs, inputs) < pDmask
# Apply masks.
- B = B * B_mask
- C = C * C_mask
- D = D * D_mask
+ B = B * Bmask
+ C = C * Cmask
+ D = D * Dmask
- return StateSpace(A, B, C, D)
+ return StateSpace(A, B, C, D)
\ No newline at end of file
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:19
|
Revision: 42
http://python-control.svn.sourceforge.net/python-control/?rev=42&view=rev
Author: kkchen
Date: 2011-02-08 22:13:13 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
Initial gram algorithm in statefbk.py. Passes initial tests.
Steven Brunton <sbr...@pr...>
Modified Paths:
--------------
branches/control-0.4a/src/statefbk.py
Modified: branches/control-0.4a/src/statefbk.py
===================================================================
--- branches/control-0.4a/src/statefbk.py 2011-02-08 22:13:09 UTC (rev 41)
+++ branches/control-0.4a/src/statefbk.py 2011-02-08 22:13:13 UTC (rev 42)
@@ -246,13 +246,44 @@
Wc = gram(sys,'c')
Wo = gram(sys,'o')
"""
-
+
+ #Check for ss system object, need a utility for this?
+
+ #TODO: Check for continous or discrete, only continuous supported right now
+ # if isCont():
+ # dico = 'C'
+ # elif isDisc():
+ # dico = 'D'
+ # else:
+ dico = 'C'
+
+ #TODO: Check system is stable, perhaps a utility in ctrlutil.py
+ # or a method of the StateSpace class?
+ D,V = np.linalg.eig(sys.A)
+ for e in D:
+ if e.real >= 0:
+ raise ValueError, "Oops, the system is unstable!"
+
if type=='c':
print "controllable"
+ trana = 'T'
+ C = -sys.B*sys.B.transpose()
elif type=='o':
print "observable"
- else:
+ trana = 'N'
+ C = -sys.C.transpose()*sys.C
+ else:
raise ValueError, "Oops, neither observable, nor controllable!"
- gram = 0.
+ #Compute Gramian by the Slycot routine sb03md
+ #make sure Slycot is installed
+ try:
+ from slycot import sb03md
+ except ImportError:
+ raise ControlSlycot("can't find slycot module 'sb03md'")
+ n = sys.states
+ U = np.zeros((n,n))
+ out = sb03md(n, C, sys.A, U, dico, 'X', 'N', trana)
+ gram = out[0]
return gram
+
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.
