Skip to content

Commit 5925723

Browse files
murrayrmmurrayrm
committed
* 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
1 parent 62cc545 commit 5925723

8 files changed

Lines changed: 160 additions & 37 deletions

File tree

‎ChangeLog‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,36 @@
1+
2012-08-29 Richard Murray <murray@altura.local>
2+
3+
* src/xferfcn.py (TransferFunction._common_den): fixed up code for
4+
case where poles have multiplicity > 1. Set end limits of check to
5+
avoid overruning the list + combine complex conjugate pole pairs
6+
from the outside in to prevent small errors from growing.
7+
8+
2012-08-26 Richard Murray <murray@altura.local>
9+
10+
* tests/convert_test.py (TestConvert.testConvert): replaced previous
11+
test of transfer function coefficients with a frequency response
12+
test. This is necessary because various degenerate conditions will
13+
generate a situation where the transfer functions are of different
14+
forms but still equal. Eg 0/1 = 1e-17 s + 1e-16 / (s^2 + s + 1).
15+
Still getting errors, but looks like actual problem in conversion.
16+
17+
* src/statesp.py (_mimo2siso): set default value for warn_conversion
18+
to False
19+
20+
* tests/convert_test.py (TestConvert.testConvert): Added check to
21+
make sure that we don't create problems with uncontrollable or
22+
unobservable systems.
23+
124
2012-08-25 Richard Murray <murray@altura.local>
225

26+
* src/xferfcn.py (TransferFunction._common_den): identified bug in
27+
the way that common denominators are computed; see comments in code
28+
regarding complex conjugate pairs.
29+
30+
* src/statefbk.py (place): repalced nd.size(placed_eigs) with
31+
len(placed_eigs) to fix intermittent problems pointed out by Roberto
32+
Bucher in control-0.3c.
33+
334
* tests/rlocus_test.py (TestRootLocus.testRootLocus): added sort()
435
to test to get rid of problems in which ordering was generating an
536
error.

‎src/__init__.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -74,9 +74,9 @@
7474
from xferfcn import TransferFunction
7575

7676
# Import some of the more common (and benign) MATLAB shortcuts
77+
# By default, don't import conflicting commands here
7778
from matlab import ss, tf, ss2tf, tf2ss, drss
7879
from matlab import pole, zero, evalfr, freqresp, dcgain
7980
from matlab import nichols, rlocus, margin
8081
# bode and nyquist come directly from freqplot.py
8182
from matlab import step, impulse, initial, lsim
82-

‎src/matlab.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -984,6 +984,7 @@ def bode(*args, **keywords):
984984

985985
# Warn about unimplemented plotstyles
986986
#! TODO: remove this when plot styles are implemented in bode()
987+
#! TODO: uncomment unit test code that tests this out
987988
if (len(plotstyle) != 0):
988989
print("Warning (matabl.bode): plot styles not implemented");
989990

‎src/statefbk.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -93,7 +93,7 @@ def place(A, B, p):
9393

9494
# Call SLICOT routine to place the eigenvalues
9595
A_z,w,nfp,nap,nup,F,Z = \
96-
sb01bd(B_mat.shape[0], B_mat.shape[1], np.size(placed_eigs), alpha,
96+
sb01bd(B_mat.shape[0], B_mat.shape[1], len(placed_eigs), alpha,
9797
A_mat, B_mat, placed_eigs, 'C');
9898

9999
# Return the gain matrix, with MATLAB gain convention

‎src/statesp.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -472,7 +472,7 @@ def _convertToStateSpace(sys, **kw):
472472
index = [len(den) - 1 for i in range(sys.outputs)]
473473
# Repeat the common denominator along the rows.
474474
den = array([den for i in range(sys.outputs)])
475-
# TODO: transfer function to state space conversion is still buggy!
475+
#! TODO: transfer function to state space conversion is still buggy!
476476
#print num
477477
#print shape(num)
478478
ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=0.0)
@@ -624,7 +624,7 @@ def _rss_generate(states, inputs, outputs, type):
624624
return StateSpace(A, B, C, D)
625625

626626
# Convert a MIMO system to a SISO system
627-
def _mimo2siso(sys, input, output, warn_conversion):
627+
def _mimo2siso(sys, input, output, warn_conversion=False):
628628
#pylint: disable=W0622
629629
"""
630630
Convert a MIMO system to a SISO system. (Convert a system with multiple

‎src/xferfcn.py‎

Lines changed: 33 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -638,14 +638,34 @@ def _common_den(self, imag_tol=None):
638638
if abs(poles[n].imag) > 10 * eps:
639639
# To prevent buildup of imaginary part error, handle complex
640640
# pole pairs together.
641-
quad = polymul([1., -poles[n]], [1., -poles[n+1]])
642-
assert all(quad.imag < 10 * eps), \
643-
"The quadratic has a nontrivial imaginary part: %g" \
644-
% quad.imag.max()
645-
quad = quad.real
646-
647-
den = polymul(den, quad)
648-
n += 2
641+
#
642+
# Because we might have repeated real parts of poles
643+
# and the fact that we are using lexigraphical
644+
# ordering, we can't just combine adjacent poles.
645+
# Instead, we have to figure out the multiplicity
646+
# first, then multiple the pairs from the outside in.
647+
648+
# Figure out the multiplicity
649+
m = 1; # multiplicity count
650+
while (n+m < len(poles) and
651+
poles[n].real == poles[n+m].real and
652+
poles[n].imag * poles[n+m].imag > 0):
653+
m += 1
654+
655+
if (m > 1):
656+
print "Found pole with multiplicity %d" % m
657+
# print "Poles = ", poles
658+
659+
# Multiple pairs from the outside in
660+
for i in range(m):
661+
quad = polymul([1., -poles[n]], [1., -poles[n+2*(m-i)-1]])
662+
assert all(quad.imag < 10 * eps), \
663+
"Quadratic has a nontrivial imaginary part: %g" \
664+
% quad.imag.max()
665+
666+
den = polymul(den, quad.real)
667+
n += 1 # move to next pair
668+
n += m # skip past conjugate pairs
649669
else:
650670
den = polymul(den, [1., -poles[n].real])
651671
n += 1
@@ -654,16 +674,19 @@ def _common_den(self, imag_tol=None):
654674
num = deepcopy(self.num)
655675
if isinstance(den,float):
656676
den = array([den])
677+
657678
for i in range(self.outputs):
658679
for j in range(self.inputs):
659680
# The common denominator has leading coefficient 1. Scale out
660681
# the existing denominator's leading coefficient.
661682
assert self.den[i][j][0], "The i = %i, j = %i denominator has \
662683
a zero leading coefficient." % (i, j)
663684
num[i][j] = num[i][j] / self.den[i][j][0]
685+
664686
# Multiply in the missing poles.
665687
for p in missingpoles[i][j]:
666688
num[i][j] = polymul(num[i][j], [1., -p])
689+
667690
# Pad all numerator polynomials with zeros so that the numerator arrays
668691
# are the same size as the denominator.
669692
for i in range(self.outputs):
@@ -672,11 +695,12 @@ def _common_den(self, imag_tol=None):
672695
if(pad>0):
673696
num[i][j] = insert(num[i][j], zeros(pad),
674697
zeros(pad))
698+
675699
# Finally, convert the numerator to a 3-D array.
676700
num = array(num)
677701
# Remove trivial imaginary parts. Check for nontrivial imaginary parts.
678702
if any(abs(num.imag) > sqrt(eps)):
679-
print ("Warning: The numerator has a nontrivial nontrivial part: %g"
703+
print ("Warning: The numerator has a nontrivial imaginary part: %g"
680704
% abs(num.imag).max())
681705
num = num.real
682706

‎tests/convert_test.py‎

Lines changed: 89 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@
1616

1717
import unittest
1818
import numpy as np
19+
import control
1920
import control.matlab as matlab
2021

2122
class TestConvert(unittest.TestCase):
@@ -27,9 +28,9 @@ def setUp(self):
2728
# Number of times to run each of the randomized tests.
2829
self.numTests = 1 #almost guarantees failure
2930
# Maximum number of states to test + 1
30-
self.maxStates = 20
31+
self.maxStates = 4
3132
# Maximum number of inputs and outputs to test + 1
32-
self.maxIO = 10
33+
self.maxIO = 5
3334
# Set to True to print systems to the output.
3435
self.debug = False
3536

@@ -42,20 +43,35 @@ def printSys(self, sys, ind):
4243

4344
def testConvert(self):
4445
"""Test state space to transfer function conversion."""
45-
#Currently it only tests that a TF->SS->TF generates an unchanged TF
4646
verbose = self.debug
47+
from control.statesp import _mimo2siso
4748

4849
#print __doc__
4950

51+
# Machine precision for floats.
52+
eps = np.finfo(float).eps
53+
5054
for states in range(1, self.maxStates):
5155
for inputs in range(1, self.maxIO):
5256
for outputs in range(1, self.maxIO):
53-
#start with a random SS system and transform to TF
54-
#then back to SS, check that the matrices are the same.
57+
# start with a random SS system and transform to TF then
58+
# back to SS, check that the matrices are the same.
5559
ssOriginal = matlab.rss(states, inputs, outputs)
5660
if (verbose):
5761
self.printSys(ssOriginal, 1)
5862

63+
# Make sure the system is not degenerate
64+
Cmat = control.ctrb(ssOriginal.A, ssOriginal.B)
65+
if (np.linalg.matrix_rank(Cmat) != states):
66+
if (verbose):
67+
print " skipping (not reachable)"
68+
continue
69+
Omat = control.obsv(ssOriginal.A, ssOriginal.C)
70+
if (np.linalg.matrix_rank(Omat) != states):
71+
if (verbose):
72+
print " skipping (not observable)"
73+
continue
74+
5975
tfOriginal = matlab.tf(ssOriginal)
6076
if (verbose):
6177
self.printSys(tfOriginal, 2)
@@ -67,27 +83,77 @@ def testConvert(self):
6783
tfTransformed = matlab.tf(ssTransformed)
6884
if (verbose):
6985
self.printSys(tfTransformed, 4)
70-
86+
87+
# Check to see if the state space systems have same dim
88+
if (ssOriginal.states != ssTransformed.states):
89+
print "WARNING: state space dimension mismatch: " + \
90+
"%d versus %d" % \
91+
(ssOriginal.states, ssTransformed.states)
92+
93+
# Now make sure the frequency responses match
94+
# Since bode() only handles SISO, go through each I/O pair
95+
# For phase, take sine and cosine to avoid +/- 360 offset
7196
for inputNum in range(inputs):
7297
for outputNum in range(outputs):
73-
np.testing.assert_array_almost_equal(\
74-
tfOriginal.num[outputNum][inputNum], \
75-
tfTransformed.num[outputNum][inputNum], \
76-
err_msg='numerator mismatch')
98+
if (verbose):
99+
print "Checking input %d, output %d" \
100+
% (inputNum, outputNum)
101+
ssorig_mag, ssorig_phase, ssorig_omega = \
102+
control.bode(_mimo2siso(ssOriginal, \
103+
inputNum, outputNum), \
104+
deg=False, Plot=False)
105+
ssorig_real = ssorig_mag * np.cos(ssorig_phase)
106+
ssorig_imag = ssorig_mag * np.sin(ssorig_phase)
107+
108+
#
109+
# Make sure TF has same frequency response
110+
#
111+
num = tfOriginal.num[outputNum][inputNum]
112+
den = tfOriginal.den[outputNum][inputNum]
113+
tforig = control.tf(num, den)
114+
115+
tforig_mag, tforig_phase, tforig_omega = \
116+
control.bode(tforig, ssorig_omega, \
117+
deg=False, Plot=False)
118+
119+
tforig_real = tforig_mag * np.cos(tforig_phase)
120+
tforig_imag = tforig_mag * np.sin(tforig_phase)
121+
np.testing.assert_array_almost_equal( \
122+
ssorig_real, tforig_real)
123+
np.testing.assert_array_almost_equal( \
124+
ssorig_imag, tforig_imag)
125+
126+
#
127+
# Make sure xform'd SS has same frequency response
128+
#
129+
ssxfrm_mag, ssxfrm_phase, ssxfrm_omega = \
130+
control.bode(_mimo2siso(ssTransformed, \
131+
inputNum, outputNum), \
132+
ssorig_omega, \
133+
deg=False, Plot=False)
134+
ssxfrm_real = ssxfrm_mag * np.cos(ssxfrm_phase)
135+
ssxfrm_imag = ssxfrm_mag * np.sin(ssxfrm_phase)
136+
np.testing.assert_array_almost_equal( \
137+
ssorig_real, ssxfrm_real)
138+
np.testing.assert_array_almost_equal( \
139+
ssorig_imag, ssxfrm_imag)
140+
141+
#
142+
# Make sure xform'd TF has same frequency response
143+
#
144+
num = tfTransformed.num[outputNum][inputNum]
145+
den = tfTransformed.den[outputNum][inputNum]
146+
tfxfrm = control.tf(num, den)
147+
tfxfrm_mag, tfxfrm_phase, tfxfrm_omega = \
148+
control.bode(tfxfrm, ssorig_omega, \
149+
deg=False, Plot=False)
77150

78-
np.testing.assert_array_almost_equal(\
79-
tfOriginal.den[outputNum][inputNum], \
80-
tfTransformed.den[outputNum][inputNum],
81-
err_msg='denominator mismatch')
82-
83-
#To test the ss systems is harder because they aren't the same
84-
#realization. This could be done with checking that they have the
85-
#same freq response with bode, but apparently it doesn't work
86-
#the way it should right now:
87-
## Bode should work like this:
88-
#[mag,phase,freq]=bode(sys)
89-
#it doesn't seem to......
90-
#This should be added.
151+
tfxfrm_real = tfxfrm_mag * np.cos(tfxfrm_phase)
152+
tfxfrm_imag = tfxfrm_mag * np.sin(tfxfrm_phase)
153+
np.testing.assert_array_almost_equal( \
154+
ssorig_real, tfxfrm_real)
155+
np.testing.assert_array_almost_equal( \
156+
ssorig_imag, tfxfrm_imag)
91157

92158
def suite():
93159
return unittest.TestLoader().loadTestsFromTestCase(TestConvert)

‎tests/matlab_test.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -238,7 +238,8 @@ def testBode(self):
238238
w = logspace(-3, 3);
239239
bode(self.siso_ss1, w)
240240
bode(self.siso_ss1, self.siso_tf2, w)
241-
bode(self.siso_ss1, '-', self.siso_tf1, 'b--', self.siso_tf2, 'k.')
241+
# Not yet implemented
242+
# bode(self.siso_ss1, '-', self.siso_tf1, 'b--', self.siso_tf2, 'k.')
242243

243244
def testRlocus(self):
244245
rlocus(self.siso_ss1)

0 commit comments

Comments
 (0)