Hi, I followed the approach mentioned above, got the following error:

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (-1,1,)

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/pi/Coding/Stability/TF_Estimate.py", line 56, in <module>
    A, B, C, D = sb10yd(0,0,9,real_H1_resp,imag_H1_resp,frequency,1,-1)
  File "/home/pi/Coding/Stability/venv/lib/python3.7/site-packages/slycot-0.4.0-py3.7-linux-armv7l.egg/slycot/synthesis.py", line 1766, in sb10yd
    A,B,C,D,info = _wrapper.sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,ldwork,lzwork)
ValueError: failed in converting hidden `a' of _wrapper.sb10yd to C/Fortran array

My example is:

omega = 2*np.pi
frequency = omega*np.arange(1,10,1)
H1 = control.frd(H,frequency)
real_H1_resp = [4.74192334e-03, 1.18688787e-03, 5.27621687e-04, 2.96810037e-04, 1.89965190e-04, 1.31922823e-04, 9.69240212e-05, 7.42080156e-05, 5.86337982e-05]#H1.fresp.real
imag_H1_resp = [-0.11917753, -0.05965949, -0.03978174, -0.0298386,  -0.02387173, -0.01989349, -0.01705176, -0.01492041, -0.01326265]#H1.fresp.imag
print(H1.fresp)
print(real_H1_resp)
print(imag_H1_resp)
print(frequency)

A, B, C, D = sb10yd(0,0,9,real_H1_resp,imag_H1_resp,frequency,1,-1)

My synthesis.pyf:

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),check(discfl>=0 && discfl<=1) :: discfl
    integer intent(in),check(discfl>=0 && discfl<=1) :: flag
    integer intent(in),check(lendat>=2) :: 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(inout),check(n>0 && n<=lendat-1) :: n
    double precision intent(out),dimension(lda,n),depend(n) :: a
    integer intent(hide),check(shape(a,0)==lda),depend(a) :: lda=shape(a,0)
    double precision intent(out),dimension(n) :: b
    double precision intent(out),dimension(n) :: c
    double precision intent(out),dimension(1) :: d
    double precision :: tol
    integer intent(hide,cache),dimension(max(2,2*n)),depend(n) :: iwork
    double precision intent(hide,cache),dimension(ldwork) :: dwork
    integer optional,check(ldwork>=max(3,6*n)),depend(n) :: ldwork
    complex*16 intent(hide,cache),dimension(lzwork) :: zwork
    integer optional,:: lzwork
    integer intent(out):: info
end subroutine sb10yd

My synthesis.py

def sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,ldwork=None,lzwork=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']

    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)
    A,B,C,D,info = _wrapper.sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,ldwork,lzwork)

    raise_if_slycot_error(info, arg_list, sb10yd.__doc__)

    return A[:nsys,:nsys],B[:nsys],C[:nsys],D[:1]

Hope you can help, thanks.
Bo

SLICOT SB10YD wrapper #153

Description

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions