Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add sb10yd #203

Merged
merged 11 commits into from
Aug 26, 2023
4 changes: 2 additions & 2 deletions slycot/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
21 changes: 21 additions & 0 deletions slycot/src/synthesis.pyf
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
114 changes: 113 additions & 1 deletion slycot/synthesis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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])

Expand Down Expand Up @@ -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:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you please cover this case in the unit tests?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have to look into that.

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])
Expand Down
173 changes: 173 additions & 0 deletions slycot/tests/test_sb10yd.py
Original file line number Diff line number Diff line change
@@ -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)
Loading