Add sb01yd, first commit

python-control · bnavigator · Aug 26, 2023 · Aug 1, 2023 · Aug 3, 2023 · Aug 3, 2023
commit aa22020271b07ca75849b0faa598c43234828afc
diff --git 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/sb10yd.pyf b/slycot/src/sb10yd.pyf
@@ -0,0 +1,28 @@
+!    -*- f90 -*-
+! Note: the context of this file is case sensitive.
+
+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 :: discfl
+    integer :: flag
+    integer :: lendat
+    double precision dimension(*) :: rfrdat
+    double precision dimension(*) :: ifrdat
+    double precision dimension(*) :: omega
+    integer :: n
+    double precision dimension(lda,*) :: a
+    integer, optional,check(shape(a, 0) == lda),depend(a) :: lda=shape(a, 0)
+    double precision dimension(*) :: b
+    double precision dimension(*) :: c
+    double precision dimension(*) :: d
+    double precision :: tol
+    integer dimension(*) :: iwork
+    double precision dimension(*) :: dwork
+    integer :: ldwork
+    complex*16 dimension(*) :: zwork
+    integer :: lzwork
+    integer :: info
+end subroutine sb10yd
+
+! This file was auto-generated with f2py (version:1.25.1).
+! See:
+! https://web.archive.org/web/20140822061353/http://cens.ioc.ee/projects/f2py2e
diff --git a/slycot/src/synthesis.pyf b/slycot/src/synthesis.pyf
@@ -1,4 +1,25 @@
 !    -*- f90 -*-
+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 sb01bd(dico,n,m,np,alpha,a,lda,b,ldb,wr,wi,nfp,nap,nup,f,ldf,z,ldz,tol,dwork,ldwork,iwarn,info) ! in SB01BD.f
     character :: dico
     integer required,check(n>=0) :: n

diff --git a/slycot/synthesis.py b/slycot/synthesis.py
@@ -25,6 +25,101 @@
 from . import _wrapper
 from .exceptions import raise_if_slycot_error, SlycotParameterError
 
+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,lzwork])
+
+    To fit frequency response data with a stable, minimum phase SISO system
+
+    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       : double precision array, dimension (lendat)
+        The real part of the frequency data to be fitted.
+    ifrdat       : double precision array, dimension (lendat)
+        The imaginary part of the frequency data to be fitted.
+    omega        : double precision array, dimension (lendat)
+        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            : integer
+        On entry, the desired order of the system to be fitted.
+        n <= lendat-1.
+    tol          : int, optional
+        The length of the cache array.
+        ldwork >= max( 1, 2*n*n + 2*n + n*MAX( 5, n + m + np ) ).
+        For good performance, ldwork must generally be larger.
+
+    Returns
+    -------
+    A            : (n, n) double precision array
+        The leading n-by-n part of this array contains the
+        matrix A.
+    B            : (n) double precision array
+        The computed vector B.
+    C            : (n) double precision array
+        The computed vector C.
+    D            : (1) double precision array
+        The computed scalar D.
+
+    Raises
+    ------
+    SlycotArithmeticError
+        :info == 0:  successful exit;
+        :info < 0:   if infor = -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__,  locals())
+
+    return out[:-1]
 
 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])

diff --git a/slycot/tests/test_sb10yd.py b/slycot/tests/test_sb10yd.py
@@ -0,0 +1,63 @@
+import unittest
+from slycot import synthesis
+import numpy as np
+from scipy import signal
+
+from numpy.testing import assert_almost_equal, assert_equal
+
+class test_sb10yd(unittest.TestCase):
+
+    def test_sb10yd_exec(self):
+        """Test execution. 
+        """
+
+        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
+        n_id, *_ = synthesis.sb10yd(
+            0, 0, len(omega), 
+            real_H_resp, imag_H_resp, omega, n, tol=0)
+
+        np.testing.assert_equal(n, n_id)
+
+    def test_sb10yd_allclose(self):
+        """Compare given and identified frequency response.
+        """
+
+        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
+        n_id, A_id, B_id, C_id, D_id = synthesis.sb10yd(
+            0, 0, 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),abs(H_id),rtol=0.3,atol=0)
+
+if __name__ == "__main__":
+    unittest.main()