Skip to content

Commit a618aeb

Browse files
roryyorkeroryyorke
authored
Merge pull request #10 from rabraker/master
Added wrapper for TB05AD, which finds the frequency response of a state-space system.
2 parents 5d21c69 + 05242f1 commit a618aeb

File tree

5 files changed

+535
-4
lines changed

5 files changed

+535
-4
lines changed

‎slycot/__init__.py‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -31,10 +31,11 @@
3131
from .synthesis import sg02ad, sg03bd
3232

3333
# Transformation routines (9/40 wrapped)
34-
from .transform import tb01id,tb03ad,tb04ad
35-
from .transform import tc04ad,tc01od
36-
from .transform import tf01md,tf01rd
37-
from .transform import td04ad,tb01pd
34+
from .transform import tb01id, tb03ad, tb04ad
35+
from .transform import tb05ad
36+
from .transform import tc04ad, tc01od
37+
from .transform import tf01md, tf01rd
38+
from .transform import td04ad, tb01pd
3839

3940
from numpy.testing import Tester
4041
test = Tester().test

‎slycot/examples.py‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -240,3 +240,31 @@ def tb01pd_example():
240240
print('reduced order', out[-2])
241241
print(out)
242242

243+
244+
def tb05ad_example():
245+
"""
246+
Example of calculating the frequency response using tb05ad
247+
on a second-order system with a natural frequency of 10 rad/s
248+
and damping ratio of 1.05.
249+
"""
250+
import numpy as np
251+
A = np.array([[0.0, 1.0],
252+
[-100.0, -20.1]])
253+
B = np.array([[0.],[100]])
254+
C = np.array([[1., 0.]])
255+
n = np.shape(A)[0]
256+
m = np.shape(B)[1]
257+
p = np.shape(C)[0]
258+
259+
jw_s = [1j*11, 1j*15]
260+
at, bt, ct, g_1, hinvb,info = slycot.tb05ad(n, m, p, jw_s[0],
261+
A, B, C, job='NG')
262+
g_2, hinv2, info = slycot.tb05ad(n, m, p, jw_s[1], at, bt, ct, job='NH')
263+
print('--- Example for tb05ad...')
264+
print('Frequency response for (A, B, C)')
265+
print('-------------------------')
266+
print('Frequency | Response')
267+
print('%s | %s '%(jw_s[0], g_1[0, 0]))
268+
print('%s | %s '%(jw_s[1], g_2[0, 0]))
269+
270+

‎slycot/src/transform.pyf‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -136,6 +136,100 @@ subroutine tb04ad_c(rowcol,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,dcoeff,lddc
136136
integer optional :: ldwork = max(1,n*(n+1)+max(n*p+2*n+max(n,p),max(3*p,m)))
137137
integer intent(out) :: info
138138
end subroutine tb04ad_c
139+
subroutine tb05ad_ag(baleig,inita,n,m,p,freq,a,lda,b,ldb,c,ldc,rcond,g,ldg,evre,evim,hinvb,ldhinv,iwork,dwork,ldwork,zwork,lzwork,info)
140+
! The case where A is is a general matrix that should be balanced
141+
! and converted to upper Hessenberg form.
142+
fortranname tb05ad
143+
character intent(hide) :: baleig = 'A'
144+
character intent(hide) :: inita = 'G'
145+
integer check(n>0) :: n
146+
integer check(m>0) :: m
147+
integer check(p>0) :: p
148+
complex*16 intent(in) :: freq
149+
double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
150+
integer intent(hide),depend(a) :: lda=shape(a,0)
151+
double precision intent(in,out,copy),dimension(n,m),depend(n,m) :: b
152+
integer intent(hide),depend(b) :: ldb=shape(b,0)
153+
double precision intent(in,out,copy),dimension(p,n),depend(n,p) :: c
154+
integer intent(hide),depend(c) :: ldc=shape(c,0)
155+
double precision intent(out) :: rcond
156+
complex*16 intent(out),dimension(ldg,m),depend(ldg,m) :: g
157+
integer intent(hide),depend(p) :: ldg = p
158+
double precision intent(out),dimension(n),depend(n):: evre
159+
double precision intent(out),dimension(n),depend(n):: evim
160+
complex*16 intent(out),dimension(ldhinv,m),depend(ldhinv,m) :: hinvb
161+
integer intent(hide),depend(n) :: ldhinv = n
162+
! cache variables
163+
integer intent(hide,cache),dimension(n) :: iwork
164+
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
165+
integer optional,depend(n) :: ldwork = 2*n
166+
complex*16 intent(hide,cache),dimension(lzwork),depend(lzwork) :: zwork
167+
integer optional,depend(n) :: lzwork = n*n+2*n
168+
integer intent(out) :: info
169+
end subroutine tb05ad_ag
170+
subroutine tb05ad_ng(baleig,inita,n,m,p,freq,a,lda,b,ldb,c,ldc,rcond,g,ldg,evre,evim,hinvb,ldhinv,iwork,dwork,ldwork,zwork,lzwork,info)
171+
! The case where A is is a general matrix that should not be
172+
! balanced but is only converted to upper Hessenberg form.
173+
fortranname tb05ad
174+
character intent(hide) :: baleig = 'N'
175+
character intent(hide) :: inita = 'G'
176+
integer check(n>0) :: n
177+
integer check(m>0) :: m
178+
integer check(p>0) :: p
179+
complex*16 intent(in) :: freq
180+
double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
181+
integer intent(hide),depend(a) :: lda=shape(a,0)
182+
double precision intent(in,out,copy),dimension(n,m),depend(n,m) :: b
183+
integer intent(hide),depend(b) :: ldb=shape(b,0)
184+
double precision intent(in,out,copy),dimension(p,n),depend(n,p) :: c
185+
integer intent(hide),depend(c) :: ldc=shape(c,0)
186+
double precision intent(hide) :: rcond
187+
complex*16 intent(out),dimension(ldg,m),depend(ldg,m) :: g
188+
integer intent(hide),depend(p) :: ldg = p
189+
double precision intent(hide),dimension(n),depend(n):: evre
190+
double precision intent(hide),dimension(n),depend(n):: evim
191+
complex*16 intent(out),dimension(ldhinv,m),depend(ldhinv,m) :: hinvb
192+
integer intent(hide),depend(n) :: ldhinv = n
193+
! cache variables
194+
integer intent(hide,cache),dimension(n) :: iwork
195+
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
196+
integer optional,depend(n) :: ldwork = 2*n
197+
complex*16 intent(hide,cache),dimension(lzwork),depend(lzwork) :: zwork
198+
integer optional,depend(n) :: lzwork = n*n+2*n
199+
integer intent(out) :: info
200+
end subroutine tb05ad_ng
201+
202+
subroutine tb05ad_nh(baleig,inita,n,m,p,freq,a,lda,b,ldb,c,ldc,rcond,g,ldg,evre,evim,hinvb,ldhinv,iwork,dwork,ldwork,zwork,lzwork,info)
203+
! The case where A is already balanced and hessenberg
204+
fortranname tb05ad
205+
character intent(hide) :: baleig = 'N'
206+
character intent(hide) :: inita = 'H'
207+
integer check(n>0) :: n
208+
integer check(m>0) :: m
209+
integer check(p>0) :: p
210+
complex*16 intent(in) :: freq
211+
double precision intent(in),dimension(n,n),depend(n) :: a
212+
integer intent(hide),depend(a) :: lda=shape(a,0)
213+
double precision intent(in),dimension(n,m),depend(n,m) :: b
214+
integer intent(hide),depend(b) :: ldb=shape(b,0)
215+
double precision intent(in),dimension(p,n),depend(n,p) :: c
216+
integer intent(hide),depend(c) :: ldc=shape(c,0)
217+
double precision intent(hide) :: rcond
218+
complex*16 intent(out),dimension(ldg,m),depend(ldg,m) :: g
219+
integer intent(hide),depend(p) :: ldg = p
220+
double precision intent(hide),dimension(n),depend(n):: evre
221+
double precision intent(hide),dimension(n),depend(n):: evim
222+
complex*16 intent(out),dimension(ldhinv,m),depend(ldhinv,m) :: hinvb
223+
integer intent(hide),depend(n) :: ldhinv = n
224+
! cache variables
225+
integer intent(hide,cache),dimension(n) :: iwork
226+
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
227+
integer optional,depend(n) :: ldwork = 2*n
228+
complex*16 intent(hide,cache),dimension(lzwork),depend(lzwork) :: zwork
229+
integer optional,depend(n) :: lzwork = n*n+2*n
230+
integer intent(out) :: info
231+
end subroutine tb05ad_h
232+
139233
subroutine tc01od_l(leri,m,p,indlim,pcoeff,ldpco1,ldpco2,qcoeff,ldqco1,ldqco2,info) ! in TC01OD.f
140234
fortranname tc01od
141235
character intent(hide) :: leri = 'L'

‎slycot/tests/test_tb05ad.py‎

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
# ===================================================
2+
# tb05ad tests
3+
import unittest
4+
from slycot import transform
5+
import numpy as np
6+
7+
from numpy.testing import assert_raises, assert_almost_equal
8+
9+
10+
# set the random seed so we can get consistent results.
11+
np.random.seed(40)
12+
CASES = {}
13+
14+
# This is a known failure for tb05ad when running job 'AG'
15+
CASES['fail1'] = {'A': np.array([[-0.5, 0., 0., 0. ],
16+
[ 0., -1., 0. , 0. ],
17+
[ 1., 0., -0.5, 0. ],
18+
[ 0., 1., 0., -1. ]]),
19+
'B': np.array([[ 1., 0.],
20+
[ 0., 1.],
21+
[ 0., 0.],
22+
[ 0., 0.]]),
23+
'C': np.array([[ 0., 1., 1., 0.],
24+
[ 0., 1., 0., 1.],
25+
[ 0., 1., 1., 1.]])}
26+
27+
n = 20
28+
p = 10
29+
m = 14
30+
31+
CASES['pass1'] = {'A': np.random.randn(n, n),
32+
'B': np.random.randn(n, m),
33+
'C': np.random.randn(p, n)}
34+
35+
36+
class test_tb05ad(unittest.TestCase):
37+
38+
def test_tb05ad_ng(self):
39+
"""
40+
Test that tb05ad with job 'NG' computes the correct
41+
frequency response.
42+
"""
43+
for key in CASES:
44+
sys = CASES[key]
45+
self.check_tb05ad_AG_NG(sys, 10*1j, 'NG')
46+
47+
@unittest.expectedFailure
48+
def test_tb05ad_ag_failure(self):
49+
""" Test tb05ad and job 'AG' (i.e., balancing enabled) fails
50+
on certain A matrices.
51+
"""
52+
self.check_tb05ad_AG_NG(CASES['fail1'], 10*1j, 'AG')
53+
54+
def test_tb05ad_nh(self):
55+
"""Test that tb05ad with job = 'NH' computes the correct
56+
frequency response after conversion to Hessenberg form.
57+
58+
First call tb05ad with job='NH' to transform to upper Hessenberg
59+
form which outputs the transformed system.
60+
Subsequently, call tb05ad with job='NH' using this transformed system.
61+
"""
62+
jomega = 10*1j
63+
for key in CASES:
64+
sys = CASES[key]
65+
sys_transformed = self.check_tb05ad_AG_NG(sys, jomega, 'NG')
66+
self.check_tb05ad_NH(sys_transformed, sys, jomega)
67+
68+
69+
def test_tb05ad_errors(self):
70+
"""
71+
Test tb05ad error handling. We give wrong inputs and
72+
and check that this raises an error.
73+
"""
74+
self.check_tb05ad_errors(CASES['pass1'])
75+
76+
def check_tb05ad_AG_NG(self, sys, jomega, job):
77+
"""
78+
Check that tb05ad computes the correct frequency response when
79+
running jobs 'AG' and/or 'NG'.
80+
81+
Inputs
82+
------
83+
84+
sys: A a dict of system matrices with keys 'A', 'B', and 'C'.
85+
jomega: A complex scalar, which is the frequency we are
86+
evaluating the system at.
87+
job: A string, either 'AG' or 'NH'
88+
89+
Returns
90+
-------
91+
sys_transformed: A dict of the system matrices which have been
92+
transformed according the job.
93+
"""
94+
n, m = sys['B'].shape
95+
p = sys['C'].shape[0]
96+
result = transform.tb05ad(n, m, p, jomega,
97+
sys['A'], sys['B'], sys['C'], job=job)
98+
g_i = result[3]
99+
hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B'])
100+
g_i_solve = sys['C'].dot(hinvb)
101+
assert_almost_equal(g_i_solve, g_i)
102+
sys_transformed = {'A': result[0], 'B': result[1], 'C': result[2]}
103+
return sys_transformed
104+
105+
def check_tb05ad_NH(self, sys_transformed, sys, jomega):
106+
"""
107+
Check tb05ad, computes the correct frequency response when
108+
job='NH' and we supply system matrices 'A', 'B', and 'C'
109+
which have been transformed by a previous call to tb05ad.
110+
We check we get the same result as computing C(sI - A)^-1B
111+
with the original system.
112+
113+
Inputs
114+
------
115+
116+
sys_transformed: A a dict of the transformed (A in upper
117+
hessenberg form) system matrices with keys
118+
'A', 'B', and 'C'.
119+
120+
sys: A dict of the original un-transformed system matrices.
121+
122+
jomega: A complex scalar, which is the frequency to evaluate at.
123+
124+
"""
125+
126+
n, m = sys_transformed['B'].shape
127+
p = sys_transformed['C'].shape[0]
128+
result = transform.tb05ad(n, m, p, jomega, sys_transformed['A'],
129+
sys_transformed['B'], sys_transformed['C'],
130+
job='NH')
131+
g_i = result[0]
132+
hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B'])
133+
g_i_solve = sys['C'].dot(hinvb)
134+
assert_almost_equal(g_i_solve, g_i)
135+
136+
def check_tb05ad_errors(self, sys):
137+
"""
138+
Check the error handling of tb05ad. We give wrong inputs and
139+
and check that this raises an error.
140+
"""
141+
n, m = sys['B'].shape
142+
p = sys['C'].shape[0]
143+
jomega = 10*1j
144+
# test error handling
145+
# wrong size A
146+
assert_raises(ValueError, transform.tb05ad, n+1, m, p,
147+
jomega, sys['A'], sys['B'], sys['C'], job='NH')
148+
# wrong size B
149+
assert_raises(ValueError, transform.tb05ad, n, m+1, p,
150+
jomega, sys['A'], sys['B'], sys['C'], job='NH')
151+
# wrong size C
152+
assert_raises(ValueError, transform.tb05ad, n, m, p+1,
153+
jomega, sys['A'], sys['B'], sys['C'], job='NH')
154+
# unrecognized job
155+
assert_raises(ValueError, transform.tb05ad, n, m, p, jomega,
156+
sys['A'], sys['B'], sys['C'], job='a')
157+
158+
159+
160+
def suite():
161+
return unittest.TestLoader().loadTestsFromTestCase(TestConvert)
162+
163+
164+
if __name__ == "__main__":
165+
unittest.main()

0 commit comments

Comments
 (0)