Menu
▾
▴
[python-control-discuss] SF.net SVN: python-control:[185] trunk
|
From: <mur...@us...> - 2012-08-30 05:44:39
|
Revision: 185
http://python-control.svn.sourceforge.net/python-control/?rev=185&view=rev
Author: murrayrm
Date: 2012-08-30 05:44:32 +0000 (Thu, 30 Aug 2012)
Log Message:
-----------
* Fixed bug in statespace -> transfer function conversation in the way that complex conjugate pairs of poles with multiplicity > 1 were handled
* Fixed bugs in tests/convert_test in the way that state space and transfer function conversions were tested; added tf to ss comparison
* Replace nd.size() with len() in place() to fix bug reported by R. Bucher
* See ChangeLog for more detailed descriptions
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/__init__.py
trunk/src/matlab.py
trunk/src/statefbk.py
trunk/src/statesp.py
trunk/src/xferfcn.py
trunk/tests/convert_test.py
trunk/tests/matlab_test.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2012-08-25 17:45:56 UTC (rev 184)
+++ trunk/ChangeLog 2012-08-30 05:44:32 UTC (rev 185)
@@ -1,5 +1,36 @@
+2012-08-29 Richard Murray <murray@altura.local>
+
+ * src/xferfcn.py (TransferFunction._common_den): fixed up code for
+ case where poles have multiplicity > 1. Set end limits of check to
+ avoid overruning the list + combine complex conjugate pole pairs
+ from the outside in to prevent small errors from growing.
+
+2012-08-26 Richard Murray <murray@altura.local>
+
+ * tests/convert_test.py (TestConvert.testConvert): replaced previous
+ test of transfer function coefficients with a frequency response
+ test. This is necessary because various degenerate conditions will
+ generate a situation where the transfer functions are of different
+ forms but still equal. Eg 0/1 = 1e-17 s + 1e-16 / (s^2 + s + 1).
+ Still getting errors, but looks like actual problem in conversion.
+
+ * src/statesp.py (_mimo2siso): set default value for warn_conversion
+ to False
+
+ * tests/convert_test.py (TestConvert.testConvert): Added check to
+ make sure that we don't create problems with uncontrollable or
+ unobservable systems.
+
2012-08-25 Richard Murray <murray@altura.local>
+ * src/xferfcn.py (TransferFunction._common_den): identified bug in
+ the way that common denominators are computed; see comments in code
+ regarding complex conjugate pairs.
+
+ * src/statefbk.py (place): repalced nd.size(placed_eigs) with
+ len(placed_eigs) to fix intermittent problems pointed out by Roberto
+ Bucher in control-0.3c.
+
* tests/rlocus_test.py (TestRootLocus.testRootLocus): added sort()
to test to get rid of problems in which ordering was generating an
error.
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2012-08-25 17:45:56 UTC (rev 184)
+++ trunk/src/__init__.py 2012-08-30 05:44:32 UTC (rev 185)
@@ -74,9 +74,9 @@
from xferfcn import TransferFunction
# Import some of the more common (and benign) MATLAB shortcuts
+# By default, don't import conflicting commands here
from matlab import ss, tf, ss2tf, tf2ss, drss
from matlab import pole, zero, evalfr, freqresp, dcgain
from matlab import nichols, rlocus, margin
# bode and nyquist come directly from freqplot.py
from matlab import step, impulse, initial, lsim
-
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2012-08-25 17:45:56 UTC (rev 184)
+++ trunk/src/matlab.py 2012-08-30 05:44:32 UTC (rev 185)
@@ -984,6 +984,7 @@
# Warn about unimplemented plotstyles
#! TODO: remove this when plot styles are implemented in bode()
+ #! TODO: uncomment unit test code that tests this out
if (len(plotstyle) != 0):
print("Warning (matabl.bode): plot styles not implemented");
Modified: trunk/src/statefbk.py
===================================================================
--- trunk/src/statefbk.py 2012-08-25 17:45:56 UTC (rev 184)
+++ trunk/src/statefbk.py 2012-08-30 05:44:32 UTC (rev 185)
@@ -93,7 +93,7 @@
# Call SLICOT routine to place the eigenvalues
A_z,w,nfp,nap,nup,F,Z = \
- sb01bd(B_mat.shape[0], B_mat.shape[1], np.size(placed_eigs), alpha,
+ sb01bd(B_mat.shape[0], B_mat.shape[1], len(placed_eigs), alpha,
A_mat, B_mat, placed_eigs, 'C');
# Return the gain matrix, with MATLAB gain convention
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2012-08-25 17:45:56 UTC (rev 184)
+++ trunk/src/statesp.py 2012-08-30 05:44:32 UTC (rev 185)
@@ -472,7 +472,7 @@
index = [len(den) - 1 for i in range(sys.outputs)]
# Repeat the common denominator along the rows.
den = array([den for i in range(sys.outputs)])
- # TODO: transfer function to state space conversion is still buggy!
+ #! TODO: transfer function to state space conversion is still buggy!
#print num
#print shape(num)
ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=0.0)
@@ -624,7 +624,7 @@
return StateSpace(A, B, C, D)
# Convert a MIMO system to a SISO system
-def _mimo2siso(sys, input, output, warn_conversion):
+def _mimo2siso(sys, input, output, warn_conversion=False):
#pylint: disable=W0622
"""
Convert a MIMO system to a SISO system. (Convert a system with multiple
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2012-08-25 17:45:56 UTC (rev 184)
+++ trunk/src/xferfcn.py 2012-08-30 05:44:32 UTC (rev 185)
@@ -638,14 +638,34 @@
if abs(poles[n].imag) > 10 * eps:
# To prevent buildup of imaginary part error, handle complex
# pole pairs together.
- quad = polymul([1., -poles[n]], [1., -poles[n+1]])
- assert all(quad.imag < 10 * eps), \
- "The quadratic has a nontrivial imaginary part: %g" \
- % quad.imag.max()
- quad = quad.real
+ #
+ # Because we might have repeated real parts of poles
+ # and the fact that we are using lexigraphical
+ # ordering, we can't just combine adjacent poles.
+ # Instead, we have to figure out the multiplicity
+ # first, then multiple the pairs from the outside in.
- den = polymul(den, quad)
- n += 2
+ # Figure out the multiplicity
+ m = 1; # multiplicity count
+ while (n+m < len(poles) and
+ poles[n].real == poles[n+m].real and
+ poles[n].imag * poles[n+m].imag > 0):
+ m += 1
+
+ if (m > 1):
+ print "Found pole with multiplicity %d" % m
+ # print "Poles = ", poles
+
+ # Multiple pairs from the outside in
+ for i in range(m):
+ quad = polymul([1., -poles[n]], [1., -poles[n+2*(m-i)-1]])
+ assert all(quad.imag < 10 * eps), \
+ "Quadratic has a nontrivial imaginary part: %g" \
+ % quad.imag.max()
+
+ den = polymul(den, quad.real)
+ n += 1 # move to next pair
+ n += m # skip past conjugate pairs
else:
den = polymul(den, [1., -poles[n].real])
n += 1
@@ -654,6 +674,7 @@
num = deepcopy(self.num)
if isinstance(den,float):
den = array([den])
+
for i in range(self.outputs):
for j in range(self.inputs):
# The common denominator has leading coefficient 1. Scale out
@@ -661,9 +682,11 @@
assert self.den[i][j][0], "The i = %i, j = %i denominator has \
a zero leading coefficient." % (i, j)
num[i][j] = num[i][j] / self.den[i][j][0]
+
# Multiply in the missing poles.
for p in missingpoles[i][j]:
num[i][j] = polymul(num[i][j], [1., -p])
+
# Pad all numerator polynomials with zeros so that the numerator arrays
# are the same size as the denominator.
for i in range(self.outputs):
@@ -672,11 +695,12 @@
if(pad>0):
num[i][j] = insert(num[i][j], zeros(pad),
zeros(pad))
+
# Finally, convert the numerator to a 3-D array.
num = array(num)
# Remove trivial imaginary parts. Check for nontrivial imaginary parts.
if any(abs(num.imag) > sqrt(eps)):
- print ("Warning: The numerator has a nontrivial nontrivial part: %g"
+ print ("Warning: The numerator has a nontrivial imaginary part: %g"
% abs(num.imag).max())
num = num.real
Modified: trunk/tests/convert_test.py
===================================================================
--- trunk/tests/convert_test.py 2012-08-25 17:45:56 UTC (rev 184)
+++ trunk/tests/convert_test.py 2012-08-30 05:44:32 UTC (rev 185)
@@ -16,6 +16,7 @@
import unittest
import numpy as np
+import control
import control.matlab as matlab
class TestConvert(unittest.TestCase):
@@ -27,9 +28,9 @@
# Number of times to run each of the randomized tests.
self.numTests = 1 #almost guarantees failure
# Maximum number of states to test + 1
- self.maxStates = 20
+ self.maxStates = 4
# Maximum number of inputs and outputs to test + 1
- self.maxIO = 10
+ self.maxIO = 5
# Set to True to print systems to the output.
self.debug = False
@@ -42,20 +43,35 @@
def testConvert(self):
"""Test state space to transfer function conversion."""
- #Currently it only tests that a TF->SS->TF generates an unchanged TF
verbose = self.debug
+ from control.statesp import _mimo2siso
#print __doc__
+ # Machine precision for floats.
+ eps = np.finfo(float).eps
+
for states in range(1, self.maxStates):
for inputs in range(1, self.maxIO):
for outputs in range(1, self.maxIO):
- #start with a random SS system and transform to TF
- #then back to SS, check that the matrices are the same.
+ # start with a random SS system and transform to TF then
+ # back to SS, check that the matrices are the same.
ssOriginal = matlab.rss(states, inputs, outputs)
if (verbose):
self.printSys(ssOriginal, 1)
+ # Make sure the system is not degenerate
+ Cmat = control.ctrb(ssOriginal.A, ssOriginal.B)
+ if (np.linalg.matrix_rank(Cmat) != states):
+ if (verbose):
+ print " skipping (not reachable)"
+ continue
+ Omat = control.obsv(ssOriginal.A, ssOriginal.C)
+ if (np.linalg.matrix_rank(Omat) != states):
+ if (verbose):
+ print " skipping (not observable)"
+ continue
+
tfOriginal = matlab.tf(ssOriginal)
if (verbose):
self.printSys(tfOriginal, 2)
@@ -67,27 +83,77 @@
tfTransformed = matlab.tf(ssTransformed)
if (verbose):
self.printSys(tfTransformed, 4)
-
+
+ # Check to see if the state space systems have same dim
+ if (ssOriginal.states != ssTransformed.states):
+ print "WARNING: state space dimension mismatch: " + \
+ "%d versus %d" % \
+ (ssOriginal.states, ssTransformed.states)
+
+ # Now make sure the frequency responses match
+ # Since bode() only handles SISO, go through each I/O pair
+ # For phase, take sine and cosine to avoid +/- 360 offset
for inputNum in range(inputs):
for outputNum in range(outputs):
- np.testing.assert_array_almost_equal(\
- tfOriginal.num[outputNum][inputNum], \
- tfTransformed.num[outputNum][inputNum], \
- err_msg='numerator mismatch')
+ if (verbose):
+ print "Checking input %d, output %d" \
+ % (inputNum, outputNum)
+ ssorig_mag, ssorig_phase, ssorig_omega = \
+ control.bode(_mimo2siso(ssOriginal, \
+ inputNum, outputNum), \
+ deg=False, Plot=False)
+ ssorig_real = ssorig_mag * np.cos(ssorig_phase)
+ ssorig_imag = ssorig_mag * np.sin(ssorig_phase)
+
+ #
+ # Make sure TF has same frequency response
+ #
+ num = tfOriginal.num[outputNum][inputNum]
+ den = tfOriginal.den[outputNum][inputNum]
+ tforig = control.tf(num, den)
+
+ tforig_mag, tforig_phase, tforig_omega = \
+ control.bode(tforig, ssorig_omega, \
+ deg=False, Plot=False)
+
+ tforig_real = tforig_mag * np.cos(tforig_phase)
+ tforig_imag = tforig_mag * np.sin(tforig_phase)
+ np.testing.assert_array_almost_equal( \
+ ssorig_real, tforig_real)
+ np.testing.assert_array_almost_equal( \
+ ssorig_imag, tforig_imag)
+
+ #
+ # Make sure xform'd SS has same frequency response
+ #
+ ssxfrm_mag, ssxfrm_phase, ssxfrm_omega = \
+ control.bode(_mimo2siso(ssTransformed, \
+ inputNum, outputNum), \
+ ssorig_omega, \
+ deg=False, Plot=False)
+ ssxfrm_real = ssxfrm_mag * np.cos(ssxfrm_phase)
+ ssxfrm_imag = ssxfrm_mag * np.sin(ssxfrm_phase)
+ np.testing.assert_array_almost_equal( \
+ ssorig_real, ssxfrm_real)
+ np.testing.assert_array_almost_equal( \
+ ssorig_imag, ssxfrm_imag)
+
+ #
+ # Make sure xform'd TF has same frequency response
+ #
+ num = tfTransformed.num[outputNum][inputNum]
+ den = tfTransformed.den[outputNum][inputNum]
+ tfxfrm = control.tf(num, den)
+ tfxfrm_mag, tfxfrm_phase, tfxfrm_omega = \
+ control.bode(tfxfrm, ssorig_omega, \
+ deg=False, Plot=False)
- np.testing.assert_array_almost_equal(\
- tfOriginal.den[outputNum][inputNum], \
- tfTransformed.den[outputNum][inputNum],
- err_msg='denominator mismatch')
-
- #To test the ss systems is harder because they aren't the same
- #realization. This could be done with checking that they have the
- #same freq response with bode, but apparently it doesn't work
- #the way it should right now:
- ## Bode should work like this:
- #[mag,phase,freq]=bode(sys)
- #it doesn't seem to......
- #This should be added.
+ tfxfrm_real = tfxfrm_mag * np.cos(tfxfrm_phase)
+ tfxfrm_imag = tfxfrm_mag * np.sin(tfxfrm_phase)
+ np.testing.assert_array_almost_equal( \
+ ssorig_real, tfxfrm_real)
+ np.testing.assert_array_almost_equal( \
+ ssorig_imag, tfxfrm_imag)
def suite():
return unittest.TestLoader().loadTestsFromTestCase(TestConvert)
Modified: trunk/tests/matlab_test.py
===================================================================
--- trunk/tests/matlab_test.py 2012-08-25 17:45:56 UTC (rev 184)
+++ trunk/tests/matlab_test.py 2012-08-30 05:44:32 UTC (rev 185)
@@ -238,7 +238,8 @@
w = logspace(-3, 3);
bode(self.siso_ss1, w)
bode(self.siso_ss1, self.siso_tf2, w)
- bode(self.siso_ss1, '-', self.siso_tf1, 'b--', self.siso_tf2, 'k.')
+# Not yet implemented
+# bode(self.siso_ss1, '-', self.siso_tf1, 'b--', self.siso_tf2, 'k.')
def testRlocus(self):
rlocus(self.siso_ss1)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
