Merge pull request #203 from KybernetikJo/add_sb10yd

Add sb10yd

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)
 

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

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])

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)