diff --git a/slycot/__init__.py b/slycot/__init__.py index 32b4a5db..9d8704a8 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -53,12 +53,12 @@ # Nonlinear Systems (0/16 wrapped) - # Synthesis routines ((15+1)/131 wrapped), sb03md57 is not part of slicot + # Synthesis routines ((16+1)/131 wrapped), sb03md57 is not part of slicot from .synthesis import (sb01bd, sb02md, sb02mt, sb02od, sb03md, sb03md57, sb03od, sb04md, sb04qd, - sb10ad, sb10dd, sb10fd, sb10hd, + sb10ad, sb10dd, sb10fd, sb10hd, sb10yd, sg02ad, sg03ad, sg03bd) diff --git a/slycot/src/synthesis.pyf b/slycot/src/synthesis.pyf index a84cd03f..1f797763 100644 --- a/slycot/src/synthesis.pyf +++ b/slycot/src/synthesis.pyf @@ -575,6 +575,27 @@ subroutine sb10jd(n,m,np,a,lda,b,ldb,c,ldc,d,ldd,e,lde,nsys,dwork,ldwork,info) ! integer required intent(in) :: ldwork integer intent(out) :: info end subroutine sb10jd +subroutine sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,a,lda,b,c,d,tol,iwork,dwork,ldwork,zwork,lzwork,info) ! in SB10YD.f + integer intent(in) :: discfl + integer intent(in) :: flag + integer intent(in) :: lendat + double precision intent(in), dimension(lendat), depend(lendat) :: rfrdat + double precision intent(in), dimension(lendat), depend(lendat) :: ifrdat + double precision intent(in), dimension(lendat), depend(lendat) :: omega + integer intent(in,out,copy), :: n + double precision intent(out), dimension(lda,n) :: a + integer, intent(hide) :: lda=max(n,1) + double precision intent(out), dimension(n,1) :: b + double precision intent(out), dimension(1,n) :: c + double precision intent(out), dimension(1,1) :: d + double precision :: tol + integer intent(hide, cache), dimension(max(2,2*n+1)), depend(n) :: iwork + double precision intent(hide, cache), dimension(ldwork), depend(ldwork) :: dwork + integer intent(in) :: ldwork + complex*16 intent(hide, cache), dimension(lzwork), depend(lzwork) :: zwork + integer intent(in):: lzwork + integer intent(out):: info +end subroutine sb10yd subroutine sg03ad(dico,job,fact,trans,uplo,n,a,lda,e,lde,q,ldq,z,ldz,x,ldx,scale,sep,ferr,alphar,alphai,beta,iwork,dwork,ldwork,info) ! in SG03AD.f character :: dico character :: job diff --git a/slycot/synthesis.py b/slycot/synthesis.py index 0ffbbd93..5f4e06e2 100644 --- a/slycot/synthesis.py +++ b/slycot/synthesis.py @@ -25,7 +25,6 @@ from . import _wrapper from .exceptions import raise_if_slycot_error, SlycotParameterError - def sb01bd(n,m,np,alpha,A,B,w,dico,tol=0.0,ldwork=None): """ A_z,w,nfp,nap,nup,F,Z = sb01bd(n,m,np,alpha,A,B,w,dico,[tol,ldwork]) @@ -1700,6 +1699,119 @@ def sb10jd(n,m,np,A,B,C,D,E,ldwork=None): return A[:nsys,:nsys],B[:nsys,:m],C[:np, :nsys],D[:np, :m] +def sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,ldwork=None): + """ A,B,C,D = sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,[ldwork]) + + To fit a supplied frequency response data with a stable, minimum + phase SISO (single-input single-output) system represented by its + matrices A, B, C, D. It handles both discrete- and continuous-time + cases. + + :: + + dx/dt = A*x + B*u + y = C*x + D*u + + :: + + x[n+1] = A*x[n] + B*u[n] + y[n] = C*x[n] + D*u[n] . + + Parameters + ---------- + discfl : int + Indicatres the type of the system, as follows: + = 0: continuous-time system; + = 1: discrete-time system. + flag : int + If flag = 0, then the system zeros and poles are not constrained. + If flag = 1, then the system zeros and poles will have negative + real parts in the continuous-time case, or moduli less than 1 in + the discrete-time case. Consequently, flag must be equal to 1 in + mu-synthesis routines. + lendat : int + The length of the vectors rfrdat, ifrdat and omega. + length >= 2. + rfrdat : dimension (lendat), array_like + The real part of the frequency data to be fitted. + ifrdat : dimension (lendat), array_like + The imaginary part of the frequency data to be fitted. + omega : dimension (lendat), array_like + The frequencies corresponding to rfrdat and ifrdat. + These value must be nonnegative and monotonically increasing. + Additionally, for discrete-time systems they must be between 0 and PI. + n : int + On entry, the desired order of the system to be fitted. + n <= lendat-1. + tol : int, optional + The length of the cache array. + ldwork : int + With None it will be automatically calculated. + For details see SLICOT help. + + Returns + ------- + n : int + The order of the obtained system. The value of n + could only be modified if n > 0 and flag = 1. + A : (n, n) ndarray + The computed matrix A. + matrix A. + B : (n, 1) ndarray + The computed column vector B. + C : (1, n) ndarray + The computed row vector C. + D : (1, 1) ndarray + The computed scalar D. + + Raises + ------ + SlycotArithmeticError + :info == 0: successful exit; + :info < 0: if info = -i, the i-th argument had an illegal value + :info = 1: if the discret --> continous transformation cannot be made; + :info = 2: if the system poles cannot be found; + :info = 3: if the inverse system cannot be found, i.e., D is (close to) zero; + :info = 4: if the system zeros cannot be found; + :info = 5: if the state-space representation of the new transfer function T(s) cannot found; + :info = 6: if the continous --> discrete transformation cannot be made. + The iteration for computing singular value + decomposition did not converge. + """ + + hidden = ' (hidden by the wrapper)' + arg_list = ['discfl', 'flag', 'lendat', 'rfrdat', 'ifrdat', 'omega', + 'n', 'A', 'lda' + hidden, 'B', 'C', 'D', + 'TOL', 'IWORK' + hidden, 'DWORK' + hidden, 'LDWORK', + 'ZWORK' + hidden, 'LZWORK', 'INFO' + hidden] + + if ldwork is None: + lw1 = 2*lendat + 4*2048 + lw2 = lendat + 6*2048 + mn = min(2*lendat,2*n+1) + if n > 0: + lw3 = 2*lendat*(2*n+1) + max(2*lendat,2*n+1) + max(mn+6*n+4,2*mn+1) + elif n == 0: + lw3 = 4*lendat + 5 + if flag == 1: + lw4 = max(n*n+5*n,6*n+1+min(1,n)) + elif flag == 0: + lw4 = 0 + ldwork = max(2, lw1, lw2, lw3, lw4) + if n > 0: + lzwork = lendat*(2*n+3) + elif n == 0: + lzwork = lendat + + out = _wrapper.sb10yd( + discfl, flag, lendat, rfrdat, ifrdat, omega, + n, + tol,ldwork,lzwork) + + raise_if_slycot_error(out[-1], arg_list, sb10yd.__doc__) + + return out[:-1] + def sg03ad(dico,job,fact,trans,uplo,N,A,E,Q,Z,X,ldwork=None): """ A,E,Q,Z,X,scale,sep,ferr,alphar,alphai,beta = sg03ad(dico,job,fact,trans,uplo,N,A,E,Q,Z,X,[ldwork]) diff --git a/slycot/tests/test_sb10yd.py b/slycot/tests/test_sb10yd.py new file mode 100644 index 00000000..8574bd51 --- /dev/null +++ b/slycot/tests/test_sb10yd.py @@ -0,0 +1,173 @@ +import numpy as np +from scipy import signal + +from slycot import synthesis + + +class Test_sb10yd(): + + # TODO: There might be better systems/filters to do these tests. + + def test_sb10yd_cont_exec_case_n0(self): + """A simple execution test. Case n=0. + """ + sys_tf = signal.TransferFunction(1, 1) + num, den = sys_tf.num, sys_tf.den + + omega, H = signal.freqs(num, den) + + real_H_resp = np.real(H) + imag_H_resp = np.imag(H) + + n = 0 + dico = 0 # 0 for continuous time + flag = 0 # 0 for no constraints on the poles + n_id, *_ = synthesis.sb10yd( + dico, flag, len(omega), + real_H_resp, imag_H_resp, omega, n, tol=0) + + # Because flag = 0, we expect n == n_id + np.testing.assert_equal(n, n_id) + + def test_sb10yd_cont_exec(self): + """A simple execution test. Case n=2. + """ + A = np.array([[0.0, 1.0], [-0.5, -0.1]]) + B = np.array([[0.0], [1.0]]) + C = np.array([[1.0, 0.0]]) + D = np.zeros((1, 1)) + + sys_tf = signal.ss2tf(A, B, C, D) + num, den = sys_tf + + omega, H = signal.freqs(num.squeeze(), den) + + real_H_resp = np.real(H) + imag_H_resp = np.imag(H) + + n = 2 + dico = 0 # 0 for continuous time + flag = 0 # 0 for no constraints on the poles + n_id, *_ = synthesis.sb10yd( + dico, flag, len(omega), + real_H_resp, imag_H_resp, omega, n, tol=0) + + # Because flag = 0, we expect n == n_id + np.testing.assert_equal(n, n_id) + + def test_sb10yd_cont_allclose(self): + """Compare given and identified frequency response. Case n=2. + """ + + A = np.array([[0.0, 1.0], [-0.5, -0.1]]) + B = np.array([[0.0], [1.0]]) + C = np.array([[1.0, 0.0]]) + D = np.zeros((1, 1)) + + sys_tf = signal.ss2tf(A, B, C, D) + num, den = sys_tf + + omega, H = signal.freqs(num.squeeze(), den) + + real_H_resp = np.real(H) + imag_H_resp = np.imag(H) + + n = 2 + dico = 0 # 0 for continuous time + flag = 0 # 0 for no constraints on the poles + n_id, A_id, B_id, C_id, D_id = synthesis.sb10yd( + dico, flag, len(omega), + real_H_resp, imag_H_resp, omega, n, tol=0) + + sys_tf_id = signal.ss2tf(A_id, B_id, C_id, D_id) + num_id, den_id = sys_tf_id + w_id, H_id = signal.freqs(num_id.squeeze(), den_id, worN=omega) + + np.testing.assert_allclose(abs(H_id), abs(H), rtol=0, atol=0.1) + + def test_sb10yd_disc_exec_case_n0(self): + """A simple execution test. Case n=0. + """ + + sys_tf = signal.TransferFunction(1, 1, dt=0.1) + num, den = sys_tf.num, sys_tf.den + + omega, H = signal.freqz(num.squeeze(), den) + + real_H_resp = np.real(H) + imag_H_resp = np.imag(H) + + n = 0 + dico = 1 # 0 for discrete time + flag = 0 # 0 for no constraints on the poles + n_id, *_ = synthesis.sb10yd( + dico, flag, len(omega), + real_H_resp, imag_H_resp, omega, n, tol=0) + + # Because flag = 0, we expect n == n_id + np.testing.assert_equal(n, n_id) + + def test_sb10yd_disc_exec(self): + """A simple execution test. Case n=2. + """ + + A = np.array([[0.0, 1.0], [-0.5, -0.1]]) + B = np.array([[0.0], [1.0]]) + C = np.array([[1.0, 0.0]]) + D = np.zeros((1, 1)) + + sys_tf = signal.ss2tf(A, B, C, D) + num, den = sys_tf + + dt = 0.1 + num, den, dt = signal.cont2discrete((num, den), dt, method="zoh") + print(den) + + omega, H = signal.freqz(num.squeeze(), den) + + real_H_resp = np.real(H) + imag_H_resp = np.imag(H) + + n = 2 + dico = 1 # 0 for discrete time + flag = 0 # 0 for no constraints on the poles + n_id, *_ = synthesis.sb10yd( + dico, flag, len(omega), + real_H_resp, imag_H_resp, omega, n, tol=0) + + # Because flag = 0, we expect n == n_id + np.testing.assert_equal(n, n_id) + + def test_sb10yd_disc_allclose(self): + """Compare given and identified frequency response. Case n=2. + """ + + A = np.array([[0.0, 1.0], [-0.5, -0.1]]) + B = np.array([[0.0], [1.0]]) + C = np.array([[1.0, 0.0]]) + D = np.zeros((1, 1)) + + sys_tf = signal.ss2tf(A, B, C, D) + num, den = sys_tf + + dt = 0.01 + num, den, dt = signal.cont2discrete((num, den), dt, method="zoh") + + omega, H = signal.freqz(num.squeeze(), den) + + real_H_resp = np.real(H) + imag_H_resp = np.imag(H) + + n = 2 + dico = 1 # 0 for discrete time + flag = 0 # 0 for no constraints on the poles + n_id, A_id, B_id, C_id, D_id = synthesis.sb10yd( + dico, flag, len(omega), + real_H_resp, imag_H_resp, omega, n, tol=0) + + sys_id = signal.dlti(A_id, B_id, C_id, D_id, dt=dt) + sys_tf_id = signal.TransferFunction(sys_id) + num_id, den_id = sys_tf_id.num, sys_tf_id.den + w_id, H_id = signal.freqz(num_id.squeeze(), den_id, worN=omega) + + np.testing.assert_allclose(abs(H_id), abs(H), rtol=0, atol=1.0)