From 05242f17ef81485c70b73d1cd1d22b3409527698 Mon Sep 17 00:00:00 2001 From: arnold Date: Tue, 28 Mar 2017 18:41:06 -0600 Subject: [PATCH] Added a wrapper for slicot routine tb05ad.f. This adds the slycot python function tb05ad in transform.py. --- slycot/__init__.py | 9 +- slycot/examples.py | 28 +++++ slycot/src/transform.pyf | 94 ++++++++++++++ slycot/tests/test_tb05ad.py | 165 ++++++++++++++++++++++++ slycot/transform.py | 243 ++++++++++++++++++++++++++++++++++++ 5 files changed, 535 insertions(+), 4 deletions(-) create mode 100644 slycot/tests/test_tb05ad.py diff --git a/slycot/__init__.py b/slycot/__init__.py index 59e5ea04..d9d23741 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -31,10 +31,11 @@ from .synthesis import sg02ad, sg03bd # Transformation routines (9/40 wrapped) - from .transform import tb01id,tb03ad,tb04ad - from .transform import tc04ad,tc01od - from .transform import tf01md,tf01rd - from .transform import td04ad,tb01pd + from .transform import tb01id, tb03ad, tb04ad + from .transform import tb05ad + from .transform import tc04ad, tc01od + from .transform import tf01md, tf01rd + from .transform import td04ad, tb01pd from numpy.testing import Tester test = Tester().test diff --git a/slycot/examples.py b/slycot/examples.py index a9ddebed..90ee67e9 100644 --- a/slycot/examples.py +++ b/slycot/examples.py @@ -240,3 +240,31 @@ def tb01pd_example(): print('reduced order', out[-2]) print(out) + +def tb05ad_example(): + """ + Example of calculating the frequency response using tb05ad + on a second-order system with a natural frequency of 10 rad/s + and damping ratio of 1.05. + """ + import numpy as np + A = np.array([[0.0, 1.0], + [-100.0, -20.1]]) + B = np.array([[0.],[100]]) + C = np.array([[1., 0.]]) + n = np.shape(A)[0] + m = np.shape(B)[1] + p = np.shape(C)[0] + + jw_s = [1j*11, 1j*15] + at, bt, ct, g_1, hinvb,info = slycot.tb05ad(n, m, p, jw_s[0], + A, B, C, job='NG') + g_2, hinv2, info = slycot.tb05ad(n, m, p, jw_s[1], at, bt, ct, job='NH') + print('--- Example for tb05ad...') + print('Frequency response for (A, B, C)') + print('-------------------------') + print('Frequency | Response') + print('%s | %s '%(jw_s[0], g_1[0, 0])) + print('%s | %s '%(jw_s[1], g_2[0, 0])) + + diff --git a/slycot/src/transform.pyf b/slycot/src/transform.pyf index be1d49bf..cb35e310 100644 --- a/slycot/src/transform.pyf +++ b/slycot/src/transform.pyf @@ -136,6 +136,100 @@ subroutine tb04ad_c(rowcol,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,dcoeff,lddc integer optional :: ldwork = max(1,n*(n+1)+max(n*p+2*n+max(n,p),max(3*p,m))) integer intent(out) :: info end subroutine tb04ad_c +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) + ! The case where A is is a general matrix that should be balanced + ! and converted to upper Hessenberg form. + fortranname tb05ad + character intent(hide) :: baleig = 'A' + character intent(hide) :: inita = 'G' + integer check(n>0) :: n + integer check(m>0) :: m + integer check(p>0) :: p + complex*16 intent(in) :: freq + double precision intent(in,out,copy),dimension(n,n),depend(n) :: a + integer intent(hide),depend(a) :: lda=shape(a,0) + double precision intent(in,out,copy),dimension(n,m),depend(n,m) :: b + integer intent(hide),depend(b) :: ldb=shape(b,0) + double precision intent(in,out,copy),dimension(p,n),depend(n,p) :: c + integer intent(hide),depend(c) :: ldc=shape(c,0) + double precision intent(out) :: rcond + complex*16 intent(out),dimension(ldg,m),depend(ldg,m) :: g + integer intent(hide),depend(p) :: ldg = p + double precision intent(out),dimension(n),depend(n):: evre + double precision intent(out),dimension(n),depend(n):: evim + complex*16 intent(out),dimension(ldhinv,m),depend(ldhinv,m) :: hinvb + integer intent(hide),depend(n) :: ldhinv = n + ! cache variables + integer intent(hide,cache),dimension(n) :: iwork + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork + integer optional,depend(n) :: ldwork = 2*n + complex*16 intent(hide,cache),dimension(lzwork),depend(lzwork) :: zwork + integer optional,depend(n) :: lzwork = n*n+2*n + integer intent(out) :: info +end subroutine tb05ad_ag +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) + ! The case where A is is a general matrix that should not be + ! balanced but is only converted to upper Hessenberg form. + fortranname tb05ad + character intent(hide) :: baleig = 'N' + character intent(hide) :: inita = 'G' + integer check(n>0) :: n + integer check(m>0) :: m + integer check(p>0) :: p + complex*16 intent(in) :: freq + double precision intent(in,out,copy),dimension(n,n),depend(n) :: a + integer intent(hide),depend(a) :: lda=shape(a,0) + double precision intent(in,out,copy),dimension(n,m),depend(n,m) :: b + integer intent(hide),depend(b) :: ldb=shape(b,0) + double precision intent(in,out,copy),dimension(p,n),depend(n,p) :: c + integer intent(hide),depend(c) :: ldc=shape(c,0) + double precision intent(hide) :: rcond + complex*16 intent(out),dimension(ldg,m),depend(ldg,m) :: g + integer intent(hide),depend(p) :: ldg = p + double precision intent(hide),dimension(n),depend(n):: evre + double precision intent(hide),dimension(n),depend(n):: evim + complex*16 intent(out),dimension(ldhinv,m),depend(ldhinv,m) :: hinvb + integer intent(hide),depend(n) :: ldhinv = n + ! cache variables + integer intent(hide,cache),dimension(n) :: iwork + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork + integer optional,depend(n) :: ldwork = 2*n + complex*16 intent(hide,cache),dimension(lzwork),depend(lzwork) :: zwork + integer optional,depend(n) :: lzwork = n*n+2*n + integer intent(out) :: info +end subroutine tb05ad_ng + +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) + ! The case where A is already balanced and hessenberg + fortranname tb05ad + character intent(hide) :: baleig = 'N' + character intent(hide) :: inita = 'H' + integer check(n>0) :: n + integer check(m>0) :: m + integer check(p>0) :: p + complex*16 intent(in) :: freq + double precision intent(in),dimension(n,n),depend(n) :: a + integer intent(hide),depend(a) :: lda=shape(a,0) + double precision intent(in),dimension(n,m),depend(n,m) :: b + integer intent(hide),depend(b) :: ldb=shape(b,0) + double precision intent(in),dimension(p,n),depend(n,p) :: c + integer intent(hide),depend(c) :: ldc=shape(c,0) + double precision intent(hide) :: rcond + complex*16 intent(out),dimension(ldg,m),depend(ldg,m) :: g + integer intent(hide),depend(p) :: ldg = p + double precision intent(hide),dimension(n),depend(n):: evre + double precision intent(hide),dimension(n),depend(n):: evim + complex*16 intent(out),dimension(ldhinv,m),depend(ldhinv,m) :: hinvb + integer intent(hide),depend(n) :: ldhinv = n + ! cache variables + integer intent(hide,cache),dimension(n) :: iwork + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork + integer optional,depend(n) :: ldwork = 2*n + complex*16 intent(hide,cache),dimension(lzwork),depend(lzwork) :: zwork + integer optional,depend(n) :: lzwork = n*n+2*n + integer intent(out) :: info +end subroutine tb05ad_h + subroutine tc01od_l(leri,m,p,indlim,pcoeff,ldpco1,ldpco2,qcoeff,ldqco1,ldqco2,info) ! in TC01OD.f fortranname tc01od character intent(hide) :: leri = 'L' diff --git a/slycot/tests/test_tb05ad.py b/slycot/tests/test_tb05ad.py new file mode 100644 index 00000000..d06146cf --- /dev/null +++ b/slycot/tests/test_tb05ad.py @@ -0,0 +1,165 @@ +# =================================================== +# tb05ad tests +import unittest +from slycot import transform +import numpy as np + +from numpy.testing import assert_raises, assert_almost_equal + + +# set the random seed so we can get consistent results. +np.random.seed(40) +CASES = {} + +# This is a known failure for tb05ad when running job 'AG' +CASES['fail1'] = {'A': np.array([[-0.5, 0., 0., 0. ], + [ 0., -1., 0. , 0. ], + [ 1., 0., -0.5, 0. ], + [ 0., 1., 0., -1. ]]), + 'B': np.array([[ 1., 0.], + [ 0., 1.], + [ 0., 0.], + [ 0., 0.]]), + 'C': np.array([[ 0., 1., 1., 0.], + [ 0., 1., 0., 1.], + [ 0., 1., 1., 1.]])} + +n = 20 +p = 10 +m = 14 + +CASES['pass1'] = {'A': np.random.randn(n, n), + 'B': np.random.randn(n, m), + 'C': np.random.randn(p, n)} + + +class test_tb05ad(unittest.TestCase): + + def test_tb05ad_ng(self): + """ + Test that tb05ad with job 'NG' computes the correct + frequency response. + """ + for key in CASES: + sys = CASES[key] + self.check_tb05ad_AG_NG(sys, 10*1j, 'NG') + + @unittest.expectedFailure + def test_tb05ad_ag_failure(self): + """ Test tb05ad and job 'AG' (i.e., balancing enabled) fails + on certain A matrices. + """ + self.check_tb05ad_AG_NG(CASES['fail1'], 10*1j, 'AG') + + def test_tb05ad_nh(self): + """Test that tb05ad with job = 'NH' computes the correct + frequency response after conversion to Hessenberg form. + + First call tb05ad with job='NH' to transform to upper Hessenberg + form which outputs the transformed system. + Subsequently, call tb05ad with job='NH' using this transformed system. + """ + jomega = 10*1j + for key in CASES: + sys = CASES[key] + sys_transformed = self.check_tb05ad_AG_NG(sys, jomega, 'NG') + self.check_tb05ad_NH(sys_transformed, sys, jomega) + + + def test_tb05ad_errors(self): + """ + Test tb05ad error handling. We give wrong inputs and + and check that this raises an error. + """ + self.check_tb05ad_errors(CASES['pass1']) + + def check_tb05ad_AG_NG(self, sys, jomega, job): + """ + Check that tb05ad computes the correct frequency response when + running jobs 'AG' and/or 'NG'. + + Inputs + ------ + + sys: A a dict of system matrices with keys 'A', 'B', and 'C'. + jomega: A complex scalar, which is the frequency we are + evaluating the system at. + job: A string, either 'AG' or 'NH' + + Returns + ------- + sys_transformed: A dict of the system matrices which have been + transformed according the job. + """ + n, m = sys['B'].shape + p = sys['C'].shape[0] + result = transform.tb05ad(n, m, p, jomega, + sys['A'], sys['B'], sys['C'], job=job) + g_i = result[3] + hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) + g_i_solve = sys['C'].dot(hinvb) + assert_almost_equal(g_i_solve, g_i) + sys_transformed = {'A': result[0], 'B': result[1], 'C': result[2]} + return sys_transformed + + def check_tb05ad_NH(self, sys_transformed, sys, jomega): + """ + Check tb05ad, computes the correct frequency response when + job='NH' and we supply system matrices 'A', 'B', and 'C' + which have been transformed by a previous call to tb05ad. + We check we get the same result as computing C(sI - A)^-1B + with the original system. + + Inputs + ------ + + sys_transformed: A a dict of the transformed (A in upper + hessenberg form) system matrices with keys + 'A', 'B', and 'C'. + + sys: A dict of the original un-transformed system matrices. + + jomega: A complex scalar, which is the frequency to evaluate at. + + """ + + n, m = sys_transformed['B'].shape + p = sys_transformed['C'].shape[0] + result = transform.tb05ad(n, m, p, jomega, sys_transformed['A'], + sys_transformed['B'], sys_transformed['C'], + job='NH') + g_i = result[0] + hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) + g_i_solve = sys['C'].dot(hinvb) + assert_almost_equal(g_i_solve, g_i) + + def check_tb05ad_errors(self, sys): + """ + Check the error handling of tb05ad. We give wrong inputs and + and check that this raises an error. + """ + n, m = sys['B'].shape + p = sys['C'].shape[0] + jomega = 10*1j + # test error handling + # wrong size A + assert_raises(ValueError, transform.tb05ad, n+1, m, p, + jomega, sys['A'], sys['B'], sys['C'], job='NH') + # wrong size B + assert_raises(ValueError, transform.tb05ad, n, m+1, p, + jomega, sys['A'], sys['B'], sys['C'], job='NH') + # wrong size C + assert_raises(ValueError, transform.tb05ad, n, m, p+1, + jomega, sys['A'], sys['B'], sys['C'], job='NH') + # unrecognized job + assert_raises(ValueError, transform.tb05ad, n, m, p, jomega, + sys['A'], sys['B'], sys['C'], job='a') + + + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestConvert) + + +if __name__ == "__main__": + unittest.main() diff --git a/slycot/transform.py b/slycot/transform.py index 1185b694..dfa6cfda 100644 --- a/slycot/transform.py +++ b/slycot/transform.py @@ -333,6 +333,249 @@ def tb04ad(n,m,p,A,B,C,D,tol1=0.0,tol2=0.0,ldwork=None): kdcoef = max(index)+1 return A[:Nr,:Nr],B[:Nr,:m],C[:p,:Nr],Nr,index,dcoeff[:porm,:kdcoef],ucoeff[:porm,:porp,:kdcoef] + +def tb05ad(n, m, p, jomega, A, B, C, job='NG'): + """tb05ad(n, m, p, jomega, A, B, C, job='NG') + + To find the complex frequency response matrix (transfer matrix) + G(freq) of the state-space representation (A,B,C) given by + -1 + G(freq) = C * ((freq*I - A) ) * B + + where A, B and C are real N-by-N, N-by-M and P-by-N matrices + respectively and freq is a complex scalar. + + Required Arguments + ------------------ + + n : integer + The number of states, i.e. the order of the state + transition matrix A. + + m : integer + The number of inputs, i.e. the number of columns in the + matrix B. + + p : integer + The number of outputs, i.e. the number of rows in the + matrix C. + + freq complex + The frequency freq at which the frequency response matrix + (transfer matrix) is to be evaluated. For continuous time + systems, this is j*omega, where omega is the frequency to + be evaluated. For discrete time systems, + freq = exp(j*omega*Ts) + + A : double precision array, dimension (n,n). + On entry, this array must contain the state transition + matrix A. + + + B : double precision array, dimension (n,m). + On entry, this array must contain the input/state matrix B. + + + C : double precision array, dimension (p,n) + On entry, of this array must contain the state/output matrix C. + + + job : string, 'AG', 'NG', or 'NH' + If job = 'AG' (i.e., 'all', 'general matrix'), the A matrix is + first balanced. The balancing transformation + is then appropriately applied to matrices B and C. The A matrix + is (again) transformed to an upper Hessenberg representation and + the B and C matrices are also transformed. In addition, + the condition number of the problem is calculated as well as the + eigenvalues of A. + + If job='NG' (i.e., 'none', 'general matrix'), no balancing is done. + Neither the condition number nor the eigenvalues are calculated. + The routine still transforms A into upper Hessenberg form. The + matrices B and C are also appropriately transformed. + + If job = 'NH' (i.e., 'none', 'hessenberg matrix'), the function + assumes the matrices have already been transformed into Hessenberg + form, i.e., by a previous function call tb05ad. If this not the + case, the routine will return a wrong result without warning. + + Returns + ------- + if job = 'AG': + -------------- + At: The A matrix which has been both balanced and + transformed to upper Hessenberg form. The balancing + transforms A according to + A1 = P^-1 * A * P. + The transformation to upper Hessenberg form then yields + At = Q^T * (P^-1 * A * P ) * Q. + Note that the lower triangle of At is in general not zero. + Rather, it contains information on the orthogonal matrix Q + used to transform A1 to Hessenberg form. See docs for lappack + DGEHRD(): + http://www.netlib.org/lapack/explore-3.1.1-html/dgehrd.f.html + However, it does not apparently contain information on P, the + matrix used in the balancing procedure. + + Bt: The matrix B transformed according to + Bt = Q^T * P^-1 * B. + + Ct: The matrix C transformed according to + Ct = C * P * Q + + rcond: RCOND contains an estimate of the reciprocal of the + condition number of matrix H with respect to inversion, where + H = (j*freq * I - A) + + g_jw: complex p-by-m array, which contains the frequency response + matrix G(freq). + + ev: Eigenvalues of the matrix A. + + hinvb : complex n-by-m array, which contains the product + -1 + H B. + + if job = 'NG': + -------------- + At: The matrix A transformed to upper Hessenberg form according + to + At = Q^T * A * Q. + The lower triangle is not zero. It containts info on the + orthoganal transformation. See docs for linpack DGEHRD() + http://www.netlib.org/lapack/explore-3.1.1-html/dgehrd.f.html + + Bt: The matrix B transformed according to + Bt = Q^T * B. + + Ct: The matrix C transformed according to + Ct = C * Q + g_jw: complex array with dim p-by-m which contains the frequency + response matrix G(freq). + + hinvb : complex array with dimension p-by-m. + This array contains the + -1 + product H B. + + if job = 'NH' + -------------- + g_jw: complex p-by-m array which contains the frequency + response matrix G(freq). + + hinvb : complex p-by-m array which contains the + -1 + product H B. + + + Raises + ------ + ValueError : e + e.info contains information about the exact type of exception. + < 0 : if info = -i, the ith argument had an illegal value; + = 1 : More than 30 iterations were required to isolate the + eigenvalues of A. The computation is continued ?. + = 2 : Either FREQ is too near to an eigenvalue of A, or RCOND + is less than the machine precision EPS. + + Example + ------- + >>> A = np.array([[0.0, 1.0], + [-100.0, -20.1]]) + >>> B = np.array([[0.],[100]]) + >>> C = np.array([[1., 0.]]) + >>> n = np.shape(A)[0] + >>> m = np.shape(B)[1] + >>> p = np.shape(C)[0] + >>> jw_s = [1j*10, 1j*15] + >>> at, bt, ct, g_1, hinvb, info = slycot.tb05ad(n, m, p, jw_s[0], + A, B, C, job='NG') + >>> g_2, hinv2,info = slycot.tb05ad(n, m, p, jw_s[1], at, bt, ct, job='NH') + + """ + def error_handler(out, arg_list): + if out[-1] < 0: + # Conform fortran 1-based argument indexing to + # to python zero indexing. + error_text = ("The following argument had an illegal value: " + + arg_list[-out[-1]-1]) + e = ValueError(error_text) + e.info = out[-1] + raise e + if out[-1] == 1: + error_text = ("More than 30 iterations are required " + "to isolate the eigenvalue of A; the computations " + "are continued.") + e = ValueError(error_text) + e.info = out[-1] + raise e + if out[-1] == 2: + error_text = ("Either FREQ is too near to an eigenvalue of A, or " + "RCOND is less than the machine precision EPS.") + e = ValueError(error_text) + e.info = out[-1] + raise e + + hidden = ' (hidden by the wrapper)' + arg_list = ['baleig'+hidden, 'inita'+hidden, 'n', 'm', 'p', 'freq', 'a', + 'lda'+hidden, 'b', 'ldb'+hidden, 'c', 'ldc'+hidden, 'rcond', + 'g', 'ldg'+hidden, 'evre', 'evim', 'hinvb', 'ldhinv'+hidden, + 'iwork'+hidden, 'dwork'+hidden, 'ldwork'+hidden, + 'zwork'+hidden, 'lzwork'+hidden, 'info'+hidden] + # Fortran function prototype: + # TB05AD(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) + + # Sanity check on matrix dimensions + if A.shape != (n, n): + e = ValueError("The shape of A is (" + str(A.shape[0]) + "," + + str(A.shape[1]) + "), but expected (" + str(n) + + "," + str(n) + ")") + raise e + + if B.shape != (n, m): + e = ValueError("The shape of B is (" + str(B.shape[0]) + "," + + str(B.shape[1]) + "), but expected (" + str(n) + + "," + str(m) + ")") + raise e + if C.shape != (p, n): + e = ValueError("The shape of C is (" + str(C.shape[0]) + "," + + str(C.shape[1]) + "), but expected (" + str(p) + + "," + str(n) + ")") + raise e + + # ---------------------------------------------------- + # Checks done, do computation. + if job == 'AG': + out = _wrapper.tb05ad_ag(n, m, p, jomega, A, B, C) + error_handler(out, arg_list) + At, Bt, Ct, rcond, g_jw, evre, evim, hinvb = out[:-1] + ev = _np.zeros(n, 'complex64') + ev.real = evre + ev.imag = evim + info = out[-1] + return At, Bt, Ct, g_jw, rcond, ev, hinvb, info + elif job == 'NG': + # use tb05ad_ng, for 'NONE' , and 'General', because balancing + # (option 'A' for 'ALL') seems to have a bug. + out = _wrapper.tb05ad_ng(n, m, p, jomega, A, B, C) + error_handler(out, arg_list) + At, Bt, Ct, g_jw, hinvb = out[:-1] + info = out[-1] + return At, Bt, Ct, g_jw, hinvb, info + elif job == 'NH': + out = _wrapper.tb05ad_nh(n, m, p, jomega, A, B, C) + error_handler(out, arg_list) + g_i, hinvb = out[:-1] + info = out[-1] + return g_i, hinvb, info + else: + error_text = ("Unrecognized job. Expected job = 'AG' or " + "job='NG' or job = 'NH' but received job=%s"%job) + e = ValueError(error_text) + raise e + + def td04ad(rowcol,m,p,index,dcoeff,ucoeff,tol=0.0,ldwork=None): """ nr,A,B,C,D = td04ad(m,p,index,dcoeff,ucoeff,[tol,ldwork])