Skip to content

Add sb10fd #118

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.

Sign up for GitHub

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

Merged
repagh merged 13 commits into from
May 11, 2020
Merged

Add sb10fd #118

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
f528c17
add sb10fd, an alternative routine for H-infinity control
repagh May 5, 2020
ec69a7a
add failing test case from python-control issue #367
repagh May 5, 2020
8bdb769
adapt testing and doc for sb10fd to new structure
repagh May 11, 2020
3955b01
Update slycot/synthesis.py
repagh May 11, 2020
7c087dc
Update slycot/synthesis.py
repagh May 11, 2020
88be4c8
Update slycot/synthesis.py
repagh May 11, 2020
825a84a
Update slycot/synthesis.py
repagh May 11, 2020
1c9d2b4
Update slycot/synthesis.py
repagh May 11, 2020
4443cc9
Update slycot/tests/test_sb.py
repagh May 11, 2020
79cb1c0
Update slycot/synthesis.py
repagh May 11, 2020
055c463
some final tweaks to docstring
repagh May 11, 2020
a9d9f9e
Apply suggestions from code review
repagh May 11, 2020
d546be0
Update slycot/__init__.py [skip ci]
repagh May 11, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
add sb10fd, an alternative routine for H-infinity control
  • Loading branch information
@repagh
repagh committed May 10, 2020
commit f528c17eed1e10839445ab2a95e2fed30abb1d17
2 changes: 1 addition & 1 deletion slycot/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@
# Synthesis routines (14/50 wrapped)
from .synthesis import sb01bd,sb02md,sb02mt,sb02od,sb03md,sb03od
from .synthesis import sb04md,sb04qd,sb10ad,sb10dd,sb10hd,sg03ad
from .synthesis import sg02ad, sg03bd
from .synthesis import sg02ad, sg03bd, sb10fd

# Transformation routines (9/40 wrapped)
from .transform import tb01id, tb03ad, tb04ad
Expand Down
31 changes: 31 additions & 0 deletions slycot/src/synthesis.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -490,6 +490,37 @@ subroutine sb10dd(n,m,np,ncon,nmeas,gamma,a,lda,b,ldb,c,ldc,d,ldd,ak,ldak,bk,ldb
logical intent(hide,cache), dimension(2*n), depend(n) :: bwork
integer intent(out) :: info
end subroutine sb10dd
subroutine sb10fd(n,m,np,ncon,nmeas,gamma,a,lda,b,ldb,c,ldc,d,ldd,ak,ldak,bk,ldbk,ck,ldck,dk,lddk,rcond,tol,iwork,dwork,ldwork,bwork,info) ! in SB10FD.f
integer intent(in),check(n>=0) :: n
integer intent(in),check(m>=0) :: m
integer intent(in),check(np>=0) :: np
integer intent(in),check(m>=ncon && ncon>=0 && np-nmeas>=ncon),depend(m,ncon):: ncon
integer intent(in),check(np>=nmeas && nmeas>=0 && m-ncon>=nmeas),depend(np,ncon):: nmeas
double precision intent(in),check(gamma>=0.0) :: gamma
double precision intent(in),dimension(n,n),depend(n) :: a
integer intent(hide),depend(a) :: lda=shape(a,0)
double precision intent(in),dimension(n,m),depend(n,m) :: b
integer intent(hide),depend(b) :: ldb=shape(b,0)
double precision intent(in),dimension(np,n),depend(np,n) :: c
integer intent(hide),depend(c) :: ldc=shape(c,0)
double precision intent(in),dimension(np,m),depend(np,m) :: d
integer intent(hide),depend(d) :: ldd=shape(d,0)
double precision intent(out),dimension(n,n) :: ak
integer intent(hide),depend(ak) :: ldak=shape(ak,0)
double precision intent(out),dimension(n,nmeas) :: bk
integer intent(hide),depend(bk) :: ldbk=shape(bk,0)
double precision intent(out),dimension(ncon,n) :: ck
integer intent(hide),depend(ck) :: ldck=shape(ck,0)
double precision intent(out),dimension(ncon,nmeas) :: dk
integer intent(hide),depend(dk) :: lddk=shape(dk,0)
double precision intent(out),dimension(4) :: rcond
double precision intent(in) :: tol
integer intent(hide,cache),dimension(max(2*max(n,m-ncon),2*max(np-nmeas,ncon)),n*n),depend(n,m,ncon,np,nmeas) :: iwork
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
integer intent(in) :: ldwork
logical intent(hide,cache),dimension(2*n),depend(n) :: bwork
integer intent(out) :: info
end subroutine sb10fd
subroutine sb10hd(n,m,np,ncon,nmeas,a,lda,b,ldb,c,ldc,d,ldd,ak,ldak,bk,ldbk,ck,ldck,dk,lddk,rcond,tol,iwork,dwork,ldwork,bwork,info) ! in :python-control:SB10HD.f
integer check(n>0) :: n
integer check(m>0) :: m
Expand Down
209 changes: 209 additions & 0 deletions slycot/synthesis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2448,3 +2448,212 @@ def sg03bd(n,m,A,E,Q,Z,B,dico,fact='N',trans='N',ldwork=None):
alpha.real = alphar[0:n]
alpha.imag = alphai[0:n]
return U,scale,alpha/beta

def sb10fd(n,m,np,ncon,nmeas,gamma,A,B,C,D,tol=0.0,ldwork=None):
""" AK,BK,CK,DK,rcond = \
sb10fd(n,m,np,ncon,nmeas,gamma,A,B,C,D,[tol,ldwork])

To compute the matrices of an H-infinity (sub)optimal n-state
controller

| AK | BK |
K = |----|----|,
| CK | DK |

using modified Glover's and Doyle's 1988 formulas, for the system

| A | B1 B2 | | A | B |
P = |----|---------| = |---|---|
| C1 | D11 D12 | | C | D |
| C2 | D21 D22 |

and for a given value of gamma, where B2 has as column size the
number of control inputs (NCON) and C2 has as row size the number
of measurements (NMEAS) being provided to the controller.

It is assumed that

(A1) (A,B2) is stabilizable and (C2,A) is detectable,

(A2) D12 is full column rank and D21 is full row rank,

(A3) | A-j*omega*I B2 | has full column rank for all omega,
| C1 D12 |

(A4) | A-j*omega*I B1 | has full row rank for all omega.
| C2 D21 |

Required arguments
------------------

n : int
The order of the system. (size of matrix A).
m : int
The column size of the matrix B
np : int
The row size of the matrix C
ncon : int
The number of control inputs. m >= ncon >= 0, np-nmeas >= ncon.
nmeas : int
The number of measurements. np >= nmeas >= 0, m-ncon >= nmeas.
gamma : float
The value of gamma. It is assumed that gamma is
sufficiently large so that the controller is admissible.
gamma >= 0.
A : rank-2 array('d'), shape (n,n)
B : rank-2 array('d'), shape (n,m)
C : rank-2 array('d'), shape (np,n)
D : rank-2 array('d'), shape (np,m)

Optional arguments
------------------

tol : float
Tolerance used for controlling the accuracy of the applied
transformations for computing the normalized form in
SLICOT Library routine SB10PD. Transformation matrices
whose reciprocal condition numbers are less than tol are
not allowed. If tol <= 0, then a default value equal to
sqrt(eps) is used, where eps is the relative machine
precision.

ldwork : int
The dimension of the cache array.
LDWORK >= N*M + NP*(N+M) + M2*M2 + NP2*NP2 +
max(1,LW1,LW2,LW3,LW4,LW5,LW6), where
LW1 = (N+NP1+1)*(N+M2) + max(3*(N+M2)+N+NP1,5*(N+M2)),
LW2 = (N+NP2)*(N+M1+1) + max(3*(N+NP2)+N+M1,5*(N+NP2)),
LW3 = M2 + NP1*NP1 + max(NP1*max(N,M1),3*M2+NP1,5*M2),
LW4 = NP2 + M1*M1 + max(max(N,NP1)*M1,3*NP2+M1,5*NP2),
LW5 = 2*N*N + N*(M+NP) +
max(1,M*M + max(2*M1,3*N*N+max(N*M,10*N*N+12*N+5)),
NP*NP + max(2*NP1,3*N*N +
max(N*NP,10*N*N+12*N+5))),
LW6 = 2*N*N + N*(M+NP) +
max(1, M2*NP2 + NP2*NP2 + M2*M2 +
max(D1*D1 + max(2*D1, (D1+D2)*NP2),
D2*D2 + max(2*D2, D2*M2), 3*N,
N*(2*NP2 + M2) +
max(2*N*M2, M2*NP2 +
max(M2*M2+3*M2, NP2*(2*NP2+
M2+max(NP2,N)))))),
with D1 = NP1 - M2, D2 = M1 - NP2,
NP1 = NP - NP2, M1 = M - M2.
For good performance, LDWORK must generally be larger.
Denoting Q = max(M1,M2,NP1,NP2), an upper bound is
2*Q*(3*Q+2*N)+max(1,(N+Q)*(N+Q+6),Q*(Q+max(N,Q,5)+1),
2*N*(N+2*Q)+max(1,4*Q*Q+
max(2*Q,3*N*N+max(2*N*Q,10*N*N+12*N+5)),
Q*(3*N+3*Q+max(2*N,4*Q+max(N,Q))))).
if the default (None) value is used, the size for good performance
is automatically used, when LDWORK is set to zero, the minimum
cache size will be used.

Return objects
--------------

Ak : rank-2 array('d'), shape (n,n)
The controller state matrix Ak.
Bk : rank-2 array('d') with bound s(n,nmeas)
The controller input matrix Bk.
Ck : rank-2 array('d'), shape (ncon,n)
The controller output matrix Ck.
Dk : rank-2 array('d'), shape (ncon,nmeas)
The controller input/output matrix Dk.
rcond : rank-1 array('d'), shape(4,)
rcond[1] contains the reciprocal condition number of the
control transformation matrix
rcond[2] contains the reciprocal condition number of the
measurement transformation matrix
rcond[3] contains an estimate of the reciprocal condition
number of the X-Riccati equation
rcond[4] contains an estimate of the reciprocal condition
number of the Y-Riccati equation

Raises
------

SlycotParameterError : e
:e.info = -i: the i-th argument had an illegal value;
SlycotArithmeticError : e
:e.info = 1:
The matrix | A-j*omega*I B2 | had no full
. | C1 D12 |
column rank in respect to the tolerance eps
:e.info = 2:
The matrix | A-j*omega*I B1 | had not full row
. | C2 D21 |
rank in respect to the tolerance EPS
:e.info = 3:
The matrix D12 has no full column rank in
respect to the tolerance tol
:e.info = 4:
The matrix D21 had no full row rank in respect
to the tolerance tol
:e.info = 5:
The singular value decomposition (SVD) algorithm
did not converge (when computing the SVD of one of
the matrices |A B2 |, |A B1 |, D12 or D21).
. |C1 D12| |C2 D21|
:e.info = 6:
The controller is not admissible (too small value
of gamma)
:e.info = 7:
The X-Riccati equation was not solved
successfully (the controller is not admissible or
there are numerical difficulties)
:e.info = 8:
The Y-Riccati equation was not solved
successfully (the controller is not admissible or
there are numerical difficulties)
:e.info = 9:
The determinant of Im2 + Tu*D11HAT*Ty*D22 is zero
"""
hidden = ' (hidden by the wrapper)'
arg_list = ('n', 'm', 'np', 'ncon', 'nmeas', 'gamma',
'A', 'LDA'+hidden, 'B', 'LDB'+hidden,
'C', 'LDC'+hidden, 'D', 'LDD'+hidden,
'AK'+hidden, 'LDAK'+hidden, 'BK'+hidden, 'LDBK'+hidden,
'CK'+hidden, 'LDCK'+hidden, 'DK'+hidden, 'LDDK'+hidden,
'RCOND'+hidden, 'tol', 'IWORK'+hidden, 'DWORK'+hidden,
'ldwork', 'BWORK'+hidden, 'INFO'+hidden)

if ldwork is None:
# M2 = NCON NP2=NMEAS M1 = M - M2
q = max(m-ncon, ncon, np-nmeas,nmeas)
ldwork = 2*q*(3*q+2*n) + max(
1,(n+q)*(n+q+6),q*(1+max(n,q,5)+1),
2*n*(n+2*q)+max(1,4*q*q+
max(2*q,3*n*n+max(2*n*q,10*n*n+12*n+5)),
q*(3*n+3*q+max(2*n,4*q+max(n,q)))))
elif ldwork == 0:
np1 = np - nmeas
np2 = nmeas
m1 = ncon
m2 = m - ncon
d1 = np1 - m2
d2 = m1 - np2
lw1 = (n+np1+1)*(n+m2) + max(3*(n+m2)+n+np1,5*(n+m2))
lw2 = (n+np2)*(n+m1+1) + max(3*(n+np2)+n+m1,5*(n+np2))
lw3 = m2 + np1*np1 + max(np1*max(n,m1),3*m2+np1,5*m2)
lw4 = np2 + m1*m1 + max(max(n,np1)*m1,3*np2+m1,5*np2)
lw5 = 2*n*n + n*(m+np) + \
max(1,m*m + max(2*m1,3*n*n+max(n*m,10*n*n+12*n+5)),
np*np + max(2*np1,3*n*n +
max(n*np,10*n*n+12*n+5)))
lw6 = 2*n*n + n*(m+np) + \
max(1, m2*np2 + np2*np2 + m2*m2 +
max(d1*d1 + max(2*d1, (d1+d2)*np2),
d2*d2 + max(2*d2, d2*m2), 3*n,
n*(2*np2 + m2) +
max(2*n*m2, m2*np2 +
max(m2*m2+3*m2, np2*(2*np2+
m2+max(np2,n))))))
ldwork = n*m + np*(n+m) + ncon*ncon + nmeas*nmeas + \
max(1,lw1,lw2,lw3,lw4,lw5,lw6)

out = _wrapper.sb10fd(n,m,np,ncon,nmeas,gamma,A,B,C,D,tol,ldwork)

raise_if_slycot_error(arg_list, out[-1], sb10fd.__doc__)

return out[:-1]
65 changes: 65 additions & 0 deletions slycot/tests/test_sb.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
import pytest
from .test_exceptions import assert_docstring_parse


def test_sb02mt():
"""Test if sb02mt is callable

Expand Down Expand Up @@ -104,6 +105,70 @@ def test_sb10jd():
assert_allclose(D_r, Dexp, atol=1e-5)


def test_sb10fd():
A = array(((-1.0, 0.0, 4.0, 5.0, -3.0, -2.0),
(-2.0, 4.0, -7.0, -2.0, 0.0, 3.0),
(-6.0, 9.0, -5.0, 0.0, 2.0, -1.0),
(-8.0, 4.0, 7.0, -1.0, -3.0, 0.0),
( 2.0, 5.0, 8.0, -9.0, 1.0, -4.0),
( 3.0, -5.0, 8.0, 0.0, 2.0, -6.0)))
B = array(((-3.0, -4.0, -2.0, 1.0, 0.0),
( 2.0, 0.0, 1.0, -5.0, 2.0),
(-5.0, -7.0, 0.0, 7.0, -2.0),
( 4.0, -6.0, 1.0, 1.0, -2.0),
(-3.0, 9.0, -8.0, 0.0, 5.0),
( 1.0, -2.0, 3.0, -6.0, -2.0)))
C = array(((1.0, -1.0, 2.0, -4.0, 0.0, -3.0),
(-3.0, 0.0, 5.0, -1.0, 1.0, 1.0),
(-7.0, 5.0, 0.0, -8.0, 2.0, -2.0),
( 9.0, -3.0, 4.0, 0.0, 3.0, 7.0),
( 0.0, 1.0, -2.0, 1.0, -6.0, -2.0)))
D = array((( 1.0, -2.0, -3.0, 0.0, 0.0),
( 0.0, 4.0, 0.0, 1.0, 0.0),
( 5.0, -3.0, -4.0, 0.0, 1.0),
( 0.0, 1.0, 0.0, 1.0, -3.0),
( 0.0, 0.0, 1.0, 7.0, 1.0)))

gamma, tol = 15.0, 0.00000001
n, m, np, ncon, nmeas = 6, 5, 5, 2, 2

assert_raises(ValueError, synthesis.sb10fd,
n, m, np, ncon, nmeas, gamma, A, B, C, D, tol, 1)
Ak, Bk, Ck, Dk, rcond = synthesis.sb10fd(
n, m, np, ncon, nmeas, gamma, A, B, C, D, tol, 900)
Ak, Bk, Ck, Dk, rcond = synthesis.sb10fd(
n, m, np, ncon, nmeas, gamma, A, B, C, D, tol, 0)
Ak, Bk, Ck, Dk, rcond = synthesis.sb10fd(
n, m, np, ncon, nmeas, gamma, A, B, C, D, tol)

Ak_ref = array((
( -2.8043, 14.7367, 4.6658, 8.1596, 0.0848, 2.5290),
( 4.6609, 3.2756, -3.5754, -2.8941, 0.2393, 8.2920),
(-15.3127, 23.5592, -7.1229, 2.7599, 5.9775, -2.0285),
(-22.0691, 16.4758, 12.5523, -16.3602, 4.4300, -3.3168),
( 30.6789, -3.9026, -1.3868, 26.2357, -8.8267, 10.4860),
( -5.7429, 0.0577, 10.8216, -11.2275, 1.5074, -10.7244)))
Bk_ref = array((
( -0.1581, -0.0793),
( -0.9237, -0.5718),
( 0.7984, 0.6627),
( 0.1145, 0.1496),
( -0.6743, -0.2376),
( 0.0196, -0.7598)))
Ck_ref = array((
( -0.2480, -0.1713, -0.0880, 0.1534, 0.5016, -0.0730),
( 2.8810, -0.3658, 1.3007, 0.3945, 1.2244, 2.5690)))
Dk_ref = array((
( 0.0554, 0.1334),
( -0.3195, 0.0333)))

assert_allclose(Ak, Ak_ref, rtol=1e-3)
assert_allclose(Bk, Bk_ref, rtol=1e-3)
assert_allclose(Ck, Ck_ref, rtol=1e-3)
assert_allclose(Dk, Dk_ref, rtol=1e-3, atol=1e-4)
assert_allclose(rcond, (1.0, 1.0, 0.011241, 0.80492e-3), rtol=1e-4)


@pytest.mark.parametrize(
'fun, exception_class, erange, checkvars',
((synthesis.sb01bd, SlycotArithmeticError, 2, {}),
Expand Down