python-control · roryyorke · May 2, 2022 · Apr 30, 2022 · May 1, 2022 · May 1, 2022
diff --git a/slycot/__init__.py b/slycot/__init__.py
@@ -15,7 +15,7 @@
     # Analysis routines (15/40 wrapped)
     from .analysis import ab01nd, ab05md, ab05nd, ab07nd, ab08nd, ab08nz
     from .analysis import ab09ad, ab09ax, ab09bd, ab09md, ab09nd
-    from .analysis import ab13bd, ab13dd, ab13ed, ab13fd
+    from .analysis import ab13bd, ab13dd, ab13ed, ab13fd, ab13md
 
 
     # Data analysis routines (0/7 wrapped)

diff --git a/slycot/analysis.py b/slycot/analysis.py
@@ -17,6 +17,8 @@
 #       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 #       MA 02110-1301, USA.
 
+import numpy as np
+
 from . import _wrapper
 from .exceptions import raise_if_slycot_error, SlycotParameterError
 
@@ -1628,7 +1630,118 @@ def ab13fd(n, A, tol = 0.0):
     raise_if_slycot_error(out[-1], arg_list, ab13fd.__doc__)
     return out[0], out[1]
 
+
+def ab13md(n, Z, nblock, itype, x=None):
+    """mubound, d, g, xout = ab13md(n, Z, nblock, itype, [x])
+
+    Find an upper bound for the structured singular value of complex
+    matrix Z and given block diagonal structure.
+
+    Parameters
+    ----------
+    n : integer
+      Order of Z; n=Z.shape[0].
+
+    Z : (n,n) complex array
+      Matrix to find structured singular value upper bound of
+
+    nblock : (m,) integer array
+      The size of the block diagonals of the uncertainty structure;
+      i.e., nblock(i)=p means that the ith block is pxp.
+
+    itype : (m,) integer array
+      The type of each block diagonal uncertainty defined in nblock.
+      itype(i)==1 means that the ith block is real, while itype(i)==2
+      means the the ith block is complex.  Real blocks must be 1x1,
+      i.e., if itype(i)==1, ntype(i) must be 1.
+
+    x : (q,) real array or None
+      If not None, must be the output of a previous call to ab13md.
+      The previous call must have been with the same values of n,
+      nblock, and itype; and the previous call's Z should be "close"
+      to the current call's Z.
+
+      q is determined by the block structure; see SLICOT AB13MD for
+      details.
+
+    Returns
+    -------
+    mubound : non-negative real scalar
+      Upper bound on structure singular value for given arguments
+
+    d, g : (n,) real arrays
+      Real arrays such that if D=np.diag(g), G=np.diag(G), and ZH = Z.T.conj(), then
+        ZH @ D**2 @ Z + 1j * (G@Z - ZH@G) - mu**2 * D**2
+      will be negative semi-definite.
+
+    xout : (q,) real array
+      For use as ``x`` argument in subsequent call to ``ab13md``.
+
+    For scalar Z and real uncertainty (ntype=1, itype=1), returns 0
+    instead of abs(Z).
+
+    Raises
+    ------
+    SlycotArithmeticError
+        :info = 1: Block sizes must be positive
+        :info = 2: Block sizes must sum to n
+        :info = 3: Real blocks must be of size 1
+        :info = 4: Block types must be 1 or 2
+        :info = 5: Error in linear equation solution
+        :info = 6: Error in eigenvalue or singular value computation
+
+    Notes
+    -----
+    This wraps SLICOT routine AB13MD, which implements the upper bound
+    of [1].
+
+    References
+    ----------
+    .. [1] Fan, M.K.H., Tits, A.L., and Doyle, J.C., "Robustness in
+       the presence of mixed parametric uncertainty and unmodeled
+       dynamics," IEEE Trans. Automatic Control, vol. AC-36, 1991,
+       pp. 25-38.
+
+    """
+    hidden = ' (hidden by the wrapper)'
+    arg_list = ['fact' + hidden, 'n' + hidden, 'z', 'ldz' + hidden,
+                'm' + hidden, 'nblock', 'itype', 'x', 'bound', 'd', 'g',
+                'iwork' + hidden, 'dwork' + hidden, 'ldwork' + hidden,
+                'zwork' + hidden, 'lzwork' + hidden, 'info' + hidden]
+
+    # prepare the "x" input and output
+
+    # x, in SLICOT, needs to be length m+mr-1. m is the length of
+    # nblock (and itype), and mr is the number of real blocks.
+
+    # In analysis.pyf x is specified as length 2*m-1, since I couldn't
+    # figure out how to express the m+mr-1 constraint there.
+
+    # The code here is to arrange for the user-visible part of x,
+    # which is length m+mr-1, to be packed into the 2*m-1-length array
+    # to pass the SLICOT routine.
+
+    m = len(nblock)
+    mr = np.count_nonzero(1==itype)
+
+    if x is None:
+        fact='N'
+        x = np.empty(2*m-1)
+    else:
+        fact='F'
+        if len(x) != m+mr-1:
+            raise ValueError(f'Require len(x)==m+mr-1, but len(x)={len(x)}, m={m}, mr={mr}')
+        x = np.concatenate([x,np.zeros(2*m-1-len(x))])
+
+    x, bound, d, g, info = _wrapper.ab13md(fact, n, Z, nblock, itype, x)
+
+    raise_if_slycot_error(info, arg_list, ab13md.__doc__)
+
+    return bound, d, g, x[:m+mr-1]
+
+
 def ag08bd(l,n,m,p,A,E,B,C,D,equil='N',tol=0.0,ldwork=None):
+
     """ Af,Ef,nrank,niz,infz,kronr,infe,kronl = ag08bd(l,n,m,p,A,E,B,C,D,[equil,tol,ldwork])
 
     To extract from the system pencil

diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
@@ -351,6 +351,25 @@ subroutine ab13dd(dico,jobe,equil,jobd,n,m,p,fpeak,a,lda,e,lde,b,ldb,c,ldc,d,ldd
     integer intent(hide),depend(n,m,p) :: lcwork = max(1,(n+m)*(n+p)+2*min(p,m)+max(p,m))
     integer intent(out) :: info
 end subroutine ab13dd
+subroutine ab13md(fact, n, z, ldz, m, nblock, itype, x, bound, d, g, iwork, dwork, ldwork, zwork, lzwork, info ) ! in AB13MD.f
+  character intent(in) :: fact
+  integer check(n>=0) :: n
+  complex*16 intent(in),dimension(n,n),depend(n) :: z
+  integer intent(hide),depend(z) :: ldz = shape(z,0)
+  integer intent(required, in) :: m
+  integer intent(in),dimension(m),depend(m) :: nblock
+  integer intent(in),dimension(m),depend(m) :: itype
+  double precision intent(in,out),dimension(2*m-1) :: x ! dim m+mr-1; mr<=m
+  double precision intent(out) :: bound
+  double precision intent(out),dimension(n),depend(n) :: d
+  double precision intent(out),dimension(n),depend(n) :: g
+  integer intent(hide,cache),dimension(max(4*m-2,n)),depend(m,n) :: iwork
+  double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
+  integer intent(hide,in) :: ldwork=2*n*n*m - n*n + 9*m*m + n*m + 11*n + 33*m - 11
+  complex*16 intent(hide,cache),dimension(lzwork) :: zwork
+  integer intent(hide,in),depend(n,m) :: lzwork=6*n*n*m + 12*n*n + 6*m + 6*n - 3
+  integer intent(out) :: info
+end subroutine ab13dd
 subroutine ab13ed(n,a,lda,low,high,tol,dwork,ldwork,info) ! in AB13ED.f
     integer check(n>=0) :: n
     double precision dimension(n,n),depend(n) :: a

diff --git a/slycot/tests/test_ab01.py b/slycot/tests/test_ab01.py
diff --git a/slycot/tests/test_ab13md.py b/slycot/tests/test_ab13md.py
@@ -0,0 +1,139 @@
+import numpy as np
+from numpy.testing import assert_allclose, assert_array_less
+
+import pytest
+
+from slycot import ab13md
+
+# References:
+# [1] Skogestand & Postlethwaite, Multivariable Feedback Control, 1996
+# [2] slycot/src/SLICOT-Reference/examples
+
+def slicot_example():
+    # [2] AB13MD.dat & TAB13MD.f
+    n = 6
+    nblock = np.array([1, 1, 2, 1, 1])
+    itype = np.array([1, 1, 2, 2, 2])
+    # this unpleasant looking array is the result of text editor
+    # search-and-replace on the Z array in AB13MD.dat
+    Z = np.array([
+        complex(-1.0e0,6.0e0),  complex(2.0e0,-3.0e0),  complex(3.0e0,8.0e0),
+        complex(3.0e0,8.0e0),   complex(-5.0e0,-9.0e0), complex(-6.0e0,2.0e0),
+        complex(4.0e0,2.0e0),   complex(-2.0e0,5.0e0),  complex(-6.0e0,-7.0e0),
+        complex(-4.0e0,11.0e0), complex(8.0e0,-7.0e0),  complex(12.0e0,-1.0e0),
+        complex(5.0e0,-4.0e0),  complex(-4.0e0,-8.0e0), complex(1.0e0,-3.0e0),
+        complex(-6.0e0,14.0e0), complex(2.0e0,-5.0e0),  complex(4.0e0,16.0e0),
+        complex(-1.0e0,6.0e0),  complex(2.0e0,-3.0e0),  complex(3.0e0,8.0e0),
+        complex(3.0e0,8.0e0),   complex(-5.0e0,-9.0e0), complex(-6.0e0,2.0e0),
+        complex(4.0e0,2.0e0),   complex(-2.0e0,5.0e0),  complex(-6.0e0,-7.0e0),
+        complex(-4.0e0,11.0e0), complex(8.0e0,-7.0e0),  complex(12.0e0,-1.0e0),
+        complex(5.0e0,-4.0e0),  complex(-4.0e0,-8.0e0), complex(1.0e0,-3.0e0),
+        complex(-6.0e0,14.0e0), complex(2.0e0,-5.0e0),  complex(4.0e0,16.0e0),
+    ])
+    Z = np.reshape(Z, (n, n))
+    return n, Z, nblock, itype
+
+
+def test_cached_inputoutput():
+    # check x, cached working area, input and output, and error
+    n, Z, nblock, itype = slicot_example()
+
+    m = len(nblock)
+    mr = np.count_nonzero(1==itype)
+
+    mu0, d0, g0, x0 = ab13md(n, Z, nblock, itype)
+    assert m+mr-1 == len(x0)
+
+    mu1, d1, g1, x1 = ab13md(n, Z, nblock, itype, x0)
+
+    assert_allclose(mu1, mu0)
+
+    with pytest.raises(ValueError):
+        mu0, d, g, x = ab13md(n, Z, nblock, itype, np.ones(mr+mr))
+
+
+class TestReference:
+    # check a few reference cases
+    def test_complex_scalar(self):
+        # [1] (8.74)
+        n = 1
+        nblock = np.array([1])
+        itype = np.array([2]) # complex
+
+        z = np.array([[1+2j]])
+        mu = ab13md(n,z,nblock,itype)[0]
+        assert_allclose(mu, abs(z))
+
+
+    @pytest.mark.xfail(reason="https://github.com/SLICOT/SLICOT-Reference/issues/4")
+    def test_real_scalar_real_uncertainty(self):
+        # [1] (8.75)
+        n=1
+        nblock=np.array([1])
+        itype=np.array([1]) # real
+        z = np.array([[5.34]])
+        mu = ab13md(n,z,nblock,itype)[0]
+        assert_allclose(mu, abs(z))
+
+
+    def test_complex_scalar_real_uncertainty(self):
+        # [1] (8.75)
+        n=1
+        nblock=np.array([1])
+        itype=np.array([1]) # real
+
+        z = np.array([[6.78j]])
+        mu = ab13md(n,z,nblock,itype)[0]
+        assert_allclose(mu, 0)
+
+
+    def test_sp85_part1(self):
+        # [1] Example 8.5 part 1, unstructured uncertainty
+        M = np.array([[2, 2], [-1, -1]])
+        muref = 3.162
+
+        n = M.shape[0]
+        nblock=np.array([2])
+        itype=np.array([2])
+
+        mu = ab13md(n, M, nblock, itype)[0]
+        assert_allclose(mu, muref, rtol=5e-4)
+
+
+    def test_sp85_part2(self):
+        # [1] Example 8.5 part 2, structured uncertainty
+        M = np.array([[2, 2], [-1, -1]])
+        muref = 3.000
+
+        n = M.shape[0]
+        nblock=np.array([1, 1])
+        itype=np.array([2, 2])
+
+        mu = ab13md(n, M, nblock, itype)[0]
+        assert_allclose(mu, muref, rtol=5e-4)
+
+    def test_slicot(self):
+        # besides muref, check that we've output d and g correctly
+
+        muref =  0.4174753408e+02 # [2] AB13MD.res
+
+        n, Z, nblock, itype = slicot_example()
+
+        mu, d, g, x = ab13md(n, Z, nblock, itype)
+
+        assert_allclose(mu, muref)
+
+        ZH = Z.T.conj()
+        D = np.diag(d)
+        G = np.diag(g)
+
+        # this matrix should be negative semi-definite
+        negsemidef = (ZH @ D**2 @ Z
+                      + 1j * (G@Z - ZH@G)
+                      - mu**2 * D**2)
+
+        # a check on the algebra, not ab13md!
+        assert_allclose(negsemidef, negsemidef.T.conj())
+
+        evals = np.linalg.eigvalsh(negsemidef)
+        assert max(evals) < np.finfo(float).eps**0.5