Skip to content

Commit 061e1e6

Browse files
author
Lauren Padilla
committed
A few changes to unittests and usage of slycot routines.
1 parent 396ac41 commit 061e1e6

File tree

5 files changed

+56
-54
lines changed

5 files changed

+56
-54
lines changed

‎src/TestSlycot.py

Lines changed: 48 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -8,19 +8,30 @@
88

99

1010
class TestSlycot(unittest.TestCase):
11+
"""TestSlycot compares transfer function and state space conversions for
12+
various numbers of inputs,outputs and states.
13+
1. Usually passes for SISO systems of any state dim, occasonally, there will be a dimension mismatch if the original randomly generated ss system is not minimal because td04ad returns a minimal system.
14+
15+
2. For small systems with many inputs, n<<m, the tests fail because td04ad returns a minimal ss system which has fewer states than the original system. It is typical for systems with many more inputs than states to have extraneous states.
16+
17+
3. For systems with larger dimensions, n~>5 and with 2 or more outputs the conversion to statespace (td04ad) intermittently results in an equivalent realization of higher order than the original tf order. We think this has to do with minimum realization tolerances in the Fortran. The algorithm doesn't recognize that two denominators are identical and so it creates a system with nearly duplicate eigenvalues and double the state dimension. This should not be a problem in the python-control usage because the common_den() method finds repeated roots within a tolerance that we specify.
18+
19+
Matlab: Matlab seems to force its statespace system output to have order less than or equal to the order of denominators provided, avoiding the problem of very large state dimension we describe in 3. It does however, still have similar problems with pole/zero cancellation such as we encounter in 2, where a statespace system may have fewer states than the original order of transfer function.
20+
"""
1121
def setUp(self):
1222
"""Define some test parameters."""
13-
self.numTests = 1
14-
self.maxStates = 10 #seems to fail rarely with 4, sometimes with 5, frequently with 6. Could it be a problem with the subsystems?
15-
self.maxIO = 10
16-
23+
self.numTests = 5
24+
self.maxStates = 10
25+
self.maxI = 1
26+
self.maxO = 1
27+
1728
def testTF(self):
1829
""" Directly tests the functions tb04ad and td04ad through direct comparison of transfer function coefficients.
1930
Similar to TestConvert, but tests at a lower level.
2031
"""
2132
for states in range(1, self.maxStates):
22-
for inputs in range(1, self.maxIO+1):
23-
for outputs in range(1, self.maxIO+1):
33+
for inputs in range(1, self.maxI+1):
34+
for outputs in range(1, self.maxO+1):
2435
for testNum in range(self.numTests):
2536

2637
ssOriginal = matlab.rss(states, inputs, outputs)
@@ -34,23 +45,25 @@ def testTF(self):
3445

3546
tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
3647
tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad('R',states,inputs,outputs,\
37-
ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=1e-10)
48+
ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0)
3849

3950
ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\
40-
= td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=1e-8)
51+
= td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=0.0)
4152

4253
tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\
4354
tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad('R',\
4455
ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\
45-
ssTransformed_D,tol1=1e-10)
46-
print 'size(Trans_A)=',ssTransformed_A.shape
47-
print 'Trans_nr=',ssTransformed_nr
48-
print 'tfOrig_index=',tfOriginal_index
49-
print 'tfOrig_ucoeff=',tfOriginal_ucoeff
50-
print 'tfOrig_dcoeff=',tfOriginal_dcoeff
51-
print 'tfTrans_index=',tfTransformed_index
52-
print 'tfTrans_ucoeff=',tfTransformed_ucoeff
53-
print 'tfTrans_dcoeff=',tfTransformed_dcoeff
56+
ssTransformed_D,tol1=0.0)
57+
#print 'size(Trans_A)=',ssTransformed_A.shape
58+
print '===== Transformed SS =========='
59+
print matlab.ss(ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D)
60+
#print 'Trans_nr=',ssTransformed_nr
61+
#print 'tfOrig_index=',tfOriginal_index
62+
#print 'tfOrig_ucoeff=',tfOriginal_ucoeff
63+
#print 'tfOrig_dcoeff=',tfOriginal_dcoeff
64+
#print 'tfTrans_index=',tfTransformed_index
65+
#print 'tfTrans_ucoeff=',tfTransformed_ucoeff
66+
#print 'tfTrans_dcoeff=',tfTransformed_dcoeff
5467
#Compare the TF directly, must match
5568
#numerators
5669
np.testing.assert_array_almost_equal(tfOriginal_ucoeff,tfTransformed_ucoeff,decimal=3)
@@ -61,24 +74,24 @@ def testFreqResp(self):
6174
"""Compare the bode reponses of the SS systems and TF systems to the original SS
6275
They generally are different realizations but have same freq resp.
6376
Currently this test may only be applied to SISO systems.
64-
77+
"""
6578
for states in range(1,self.maxStates):
6679
for testNum in range(self.numTests):
67-
for inputs in range(1,self.maxIO+1):
68-
for outputs in range(1,self.maxIO+1):
80+
for inputs in range(1,1):
81+
for outputs in range(1,1):
6982
ssOriginal = matlab.rss(states, inputs, outputs)
7083

7184
tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
7285
tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad('R',states,inputs,outputs,\
73-
ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=1e-10)
86+
ssOriginal.A,ssOriginal.B,ssOriginal.C,ssOriginal.D,tol1=0.0)
7487

7588
ssTransformed_nr, ssTransformed_A, ssTransformed_B, ssTransformed_C, ssTransformed_D\
76-
= td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=1e-8)
89+
= td04ad('R',inputs,outputs,tfOriginal_index,tfOriginal_dcoeff,tfOriginal_ucoeff,tol=0.0)
7790

7891
tfTransformed_Actrb, tfTransformed_Bctrb, tfTransformed_Cctrb, tfTransformed_nctrb,\
7992
tfTransformed_index, tfTransformed_dcoeff, tfTransformed_ucoeff = tb04ad('R',\
8093
ssTransformed_nr,inputs,outputs,ssTransformed_A, ssTransformed_B, ssTransformed_C,\
81-
ssTransformed_D,tol1=1e-10)
94+
ssTransformed_D,tol1=0.0)
8295

8396
numTransformed = np.array(tfTransformed_ucoeff)
8497
denTransformed = np.array(tfTransformed_dcoeff)
@@ -88,21 +101,21 @@ def testFreqResp(self):
88101
ssTransformed = matlab.ss(ssTransformed_A,ssTransformed_B,ssTransformed_C,ssTransformed_D)
89102
for inputNum in range(inputs):
90103
for outputNum in range(outputs):
91-
#[ssOriginalMag,ssOriginalPhase,freq] = matlab.bode(ssOriginal,Plot=False)
104+
[ssOriginalMag,ssOriginalPhase,freq] = matlab.bode(ssOriginal,Plot=False)
92105
[tfOriginalMag,tfOriginalPhase,freq] = matlab.bode(matlab.tf(numOriginal[outputNum][inputNum],denOriginal[outputNum]),Plot=False)
93-
#[ssTransformedMag,ssTransformedPhase,freq] = matlab.bode(ssTransformed,freq,Plot=False)
106+
[ssTransformedMag,ssTransformedPhase,freq] = matlab.bode(ssTransformed,freq,Plot=False)
94107
[tfTransformedMag,tfTransformedPhase,freq] = matlab.bode(matlab.tf(numTransformed[outputNum][inputNum],denTransformed[outputNum]),freq,Plot=False)
95-
print 'numOrig=',numOriginal[outputNum][inputNum]
96-
print 'denOrig=',denOriginal[outputNum]
97-
print 'numTrans=',numTransformed[outputNum][inputNum]
98-
print 'denTrans=',denTransformed[outputNum]
99-
#np.testing.assert_array_almost_equal(ssOriginalMag,tfOriginalMag,decimal=3)
100-
#np.testing.assert_array_almost_equal(ssOriginalPhase,tfOriginalPhase,decimal=3)
101-
#np.testing.assert_array_almost_equal(ssOriginalMag,ssTransformedMag,decimal=3)
102-
#np.testing.assert_array_almost_equal(ssOriginalPhase,ssTransformedPhase,decimal=3)
103-
#np.testing.assert_array_almost_equal(tfOriginalMag,tfTransformedMag,decimal=3)
108+
#print 'numOrig=',numOriginal[outputNum][inputNum]
109+
#print 'denOrig=',denOriginal[outputNum]
110+
#print 'numTrans=',numTransformed[outputNum][inputNum]
111+
#print 'denTrans=',denTransformed[outputNum]
112+
np.testing.assert_array_almost_equal(ssOriginalMag,tfOriginalMag,decimal=3)
113+
np.testing.assert_array_almost_equal(ssOriginalPhase,tfOriginalPhase,decimal=3)
114+
np.testing.assert_array_almost_equal(ssOriginalMag,ssTransformedMag,decimal=3)
115+
np.testing.assert_array_almost_equal(ssOriginalPhase,ssTransformedPhase,decimal=3)
116+
np.testing.assert_array_almost_equal(tfOriginalMag,tfTransformedMag,decimal=3)
104117
np.testing.assert_array_almost_equal(tfOriginalPhase,tfTransformedPhase,decimal=2)
105-
"""
118+
106119
#These are here for once the above is made into a unittest.
107120
def suite():
108121
return unittest.TestLoader().loadTestsFromTestCase(TestSlycot)

‎src/modelsimp.py

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -222,19 +222,10 @@ def balred(sys,orders,method='truncate'):
222222
raise ControlSlycot("can't find slycot subroutine ab09ad")
223223
job = 'B' # balanced (B) or not (N)
224224
equil = 'N' # scale (S) or not (N)
225-
ordsel = 'F' # fixed truncation level (F) or find the truncation level given tol (A)
226225
n = np.size(sys.A,0)
227226
m = np.size(sys.B,1)
228227
p = np.size(sys.C,0)
229-
nr = orders
230-
tol = 0.
231-
out = ab09ad(dico,job,equil,ordsel, n, m, p, nr, sys.A, sys.B, sys.C,tol)
232-
Ar = out[0][0:nr,0:nr]
233-
Br = out[1][0:nr,0:m]
234-
Cr = out[2][0:p,0:nr]
235-
hsv = out[3]
236-
iwarn = out[4]
237-
info = out[5]
228+
Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=orders,tol=0.0)
238229

239230
rsys = StateSpace(Ar, Br, Cr, sys.D)
240231
else:

‎src/statefbk.py

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -177,7 +177,7 @@ def lqr(*args, **keywords):
177177
sb02mt(nstates, ninputs, B, R, A, Q, N, jobl='N');
178178

179179
# Call the SLICOT function
180-
X,rcond,w,S,U = sb02md(nstates, A_b, G, Q_b, 'C')
180+
X,rcond,w,S,U,A_inv = sb02md(nstates, A_b, G, Q_b, 'C')
181181

182182
# Now compute the return value
183183
K = np.dot(np.linalg.inv(R), (np.dot(B.T, X) + N.T));
@@ -288,10 +288,10 @@ def gram(sys,type):
288288
if e.real >= 0:
289289
raise ValueError, "Oops, the system is unstable!"
290290
if type=='c':
291-
trana = 'T'
291+
tra = 'T'
292292
C = -np.dot(sys.B,sys.B.transpose())
293293
elif type=='o':
294-
trana = 'N'
294+
tra = 'N'
295295
C = -np.dot(sys.C.transpose(),sys.C)
296296
else:
297297
raise ValueError, "Oops, neither observable, nor controllable!"
@@ -305,7 +305,7 @@ def gram(sys,type):
305305
n = sys.states
306306
U = np.zeros((n,n))
307307
A = sys.A
308-
out = sb03md(n, C, A, U, dico, 'X', 'N', trana)
309-
gram = out[0]
308+
X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
309+
gram = X
310310
return gram
311311

‎src/statesp.py

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -163,7 +163,6 @@ def _remove_useless_states(self):
163163
useless = []
164164

165165
# Search for useless states.
166-
tol = 1e-16
167166
for i in range(self.states):
168167
if (all(self.A[i, :] == zeros((1, self.states))) and
169168
all(self.B[i, :] == zeros((1, self.inputs)))):
@@ -474,7 +473,7 @@ def _convertToStateSpace(sys, **kw):
474473
is still buggy! Advise converting state space sys back to tf to verify the transformation was correct."
475474
#print num
476475
#print shape(num)
477-
ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=1e-8)
476+
ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=0.0)
478477

479478
states = ssout[0]
480479
return StateSpace(ssout[1][:states, :states],

‎src/xferfcn.py

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -212,7 +212,6 @@ def _truncatecoeff(self):
212212
"""
213213

214214
# Beware: this is a shallow copy. This should be okay.
215-
tol = 1e-16
216215
data = [self.num, self.den]
217216
for p in range(len(data)):
218217
for i in range(self.outputs):
@@ -717,7 +716,7 @@ def _convertToTransferFunction(sys, **kw):
717716
print "Warning: state space to transfer function conversion by tb04ad \
718717
is still buggy!"
719718
tfout = tb04ad('R',sys.states, sys.inputs, sys.outputs, sys.A, sys.B, sys.C,
720-
sys.D,tol1=1e-10)
719+
sys.D,tol1=0.0)
721720

722721
# Preallocate outputs.
723722
num = [[[] for j in range(sys.inputs)] for i in range(sys.outputs)]

0 commit comments

Comments
 (0)