python-control · henriks76 · Jan 7, 2024 · Jan 8, 2024 · Jan 9, 2024 · Jan 21, 2024
diff --git a/.gitignore b/.gitignore
@@ -31,6 +31,9 @@ TAGS
 # Files created by Spyder
 .spyproject/
 
+# Files created by or for VS Code (HS, 13 Jan, 2024)
+.vscode/
+
 # Environments
 .env
 .venv

diff --git a/control/__init__.py b/control/__init__.py
@@ -101,6 +101,7 @@
 from .config import *
 from .sisotool import *
 from .passivity import *
+from .sysnorm import *
 
 # Exceptions
 from .exception import *

diff --git a/control/sysnorm.py b/control/sysnorm.py
@@ -0,0 +1,306 @@
+# -*- coding: utf-8 -*-
+"""sysnorm.py
+
+Functions for computing system norms.
+
+Routine in this module:
+
+norm
+
+Created on Thu Dec 21 08:06:12 2023
+Author: Henrik Sandberg
+"""
+
+import numpy as np
+import scipy as sp
+import numpy.linalg as la
+import warnings
+
+import control as ct
+
+__all__ = ['norm']
+
+#------------------------------------------------------------------------------
+
+def _slycot_or_scipy(method):
+    """ Copied from ct.mateqn. For internal use."""
+
+    if method == 'slycot' or (method is None and ct.slycot_check()):
+        return 'slycot'
+    elif method == 'scipy' or (method is None and not ct.slycot_check()):
+        return 'scipy'
+    else:
+        raise ct.ControlArgument(f"Unknown argument '{method}'.")
+
+#------------------------------------------------------------------------------
+
+def _h2norm_slycot(sys, print_warning=True):
+    """H2 norm of a linear system. For internal use. Requires Slycot.
+
+    See also
+    --------
+    ``slycot.ab13bd`` : the Slycot routine that does the calculation
+    https://github.com/python-control/Slycot/issues/199 : Post on issue with ``ab13bf``
+    """
+
+    try:
+        from slycot import ab13bd
+    except ImportError:
+        ct.ControlSlycot("Can't find slycot module ``ab13bd``!")
+
+    try:
+        from slycot.exceptions import SlycotArithmeticError
+    except ImportError: 
+        raise ct.ControlSlycot("Can't find slycot class ``SlycotArithmeticError``!")
+
+    A, B, C, D = ct.ssdata(ct.ss(sys))
+
+    n = A.shape[0]
+    m = B.shape[1]
+    p = C.shape[0]
+
+    dico = 'C' if sys.isctime() else 'D'  # Continuous or discrete time
+    jobn = 'H'  # H2 (and not L2 norm)  
+
+    if n == 0:
+    # ab13bd does not accept empty A, B, C
+        if dico == 'C':
+            if any(D.flat != 0):
+                if print_warning:
+                    warnings.warn("System has a direct feedthrough term!")
+                return float("inf")
+            else:
+                return 0.0
+        elif dico == 'D':
+            return np.sqrt(D@D.T)
+
+    try:
+        norm = ab13bd(dico, jobn, n, m, p, A, B, C, D)
+    except SlycotArithmeticError as e:
+        if e.info == 3:
+            if print_warning:
+                warnings.warn("System has pole(s) on the stability boundary!")
+            return float("inf")
+        elif e.info == 5:
+            if print_warning:
+                warnings.warn("System has a direct feedthrough term!")
+            return float("inf")
+        elif e.info == 6:
+            if print_warning:
+                warnings.warn("System is unstable!")
+            return float("inf")
+        else:
+            raise e
+    return norm
+
+#------------------------------------------------------------------------------
+
+def norm(system, p=2, tol=1e-10, print_warning=True, method=None):
+    """Computes norm of system.
+
+    Parameters
+    ----------
+    system : LTI (:class:`StateSpace` or :class:`TransferFunction`)
+        System in continuous or discrete time for which the norm should be computed.
+    p : int or str
+        Type of norm to be computed. p=2 gives the H2 norm, and p='inf' gives the L-infinity norm.
+    tol : float
+        Relative tolerance for accuracy of L-infinity norm computation. Ignored
+        unless p='inf'.
+    print_warning : bool
+        Print warning message in case norm value may be uncertain.
+    method : str, optional
+        Set the method used for computing the result.  Current methods are
+        'slycot' and 'scipy'. If set to None (default), try 'slycot' first
+        and then 'scipy'.
+
+    Returns
+    -------
+    norm_value : float
+        Norm value of system.
+
+    Notes
+    -----
+    Does not yet compute the L-infinity norm for discrete time systems with pole(s) in z=0 unless Slycot is used.
+
+    Examples
+    --------
+    >>> Gc = ct.tf([1], [1, 2, 1])
+    >>> ct.norm(Gc, 2)
+    0.5000000000000001
+    >>> ct.norm(Gc, 'inf', tol=1e-11, method='scipy')
+    1.000000000007276
+    """
+
+    if not isinstance(system, (ct.StateSpace, ct.TransferFunction)):
+        raise TypeError('Parameter ``system``: must be a ``StateSpace`` or ``TransferFunction``')
+
+    G = ct.ss(system)
+    A = G.A
+    B = G.B
+    C = G.C
+    D = G.D
+
+    # Decide what method to use
+    method = _slycot_or_scipy(method)
+
+    # -------------------
+    # H2 norm computation
+    # -------------------
+    if p == 2:  
+        # --------------------
+        # Continuous time case
+        # --------------------
+        if G.isctime():
+
+            # Check for cases with infinite norm
+            poles_real_part = G.poles().real    
+            if any(np.isclose(poles_real_part, 0.0)):  # Poles on imaginary axis
+                if print_warning:
+                    warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.")
+                return float('inf')
+            elif any(poles_real_part > 0.0):  # System unstable
+                if print_warning:
+                    warnings.warn("System is unstable!")
+                return float('inf')
+            elif any(D.flat != 0):  # System has direct feedthrough
+                if print_warning:
+                    warnings.warn("System has a direct feedthrough term!")
+                return float('inf')   
+
+            else:      
+                # Use slycot, if available, to compute (finite) norm
+                if method == 'slycot':
+                    return _h2norm_slycot(G, print_warning)  
+
+                # Else use scipy 
+                else:
+                    P = ct.lyap(A, B@B.T, method=method)  # Solve for controllability Gramian
+
+                    # System is stable to reach this point, and P should be positive semi-definite. 
+                    # Test next is a precaution in case the Lyapunov equation is ill conditioned.
+                    if any(la.eigvals(P).real < 0.0):  
+                        if print_warning:
+                            warnings.warn("There appears to be poles close to the imaginary axis. Norm value may be uncertain.")
+                        return float('inf')
+                    else:
+                        norm_value = np.sqrt(np.trace(C@P@C.T))  # Argument in sqrt should be non-negative
+                        if np.isnan(norm_value):
+                            raise ct.ControlArgument("Norm computation resulted in NaN.")           
+                        else:
+                            return norm_value
+
+        # ------------------
+        # Discrete time case
+        # ------------------
+        elif G.isdtime():
+
+            # Check for cases with infinite norm
+            poles_abs = abs(G.poles())
+            if any(np.isclose(poles_abs, 1.0)):  # Poles on imaginary axis
+                if print_warning:
+                    warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
+                return float('inf')
+            elif any(poles_abs > 1.0):  # System unstable
+                if print_warning:
+                    warnings.warn("System is unstable!")
+                return float('inf')
+
+            else:
+                # Use slycot, if available, to compute (finite) norm
+                if method == 'slycot':
+                    return _h2norm_slycot(G, print_warning)  
+
+                # Else use scipy 
+                else:
+                    P = ct.dlyap(A, B@B.T, method=method)
+
+                # System is stable to reach this point, and P should be positive semi-definite. 
+                # Test next is a precaution in case the Lyapunov equation is ill conditioned.
+                if any(la.eigvals(P).real < 0.0):
+                    if print_warning:
+                        warnings.warn("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.")
+                    return float('inf')
+                else:
+                    norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T))  # Argument in sqrt should be non-negative               
+                    if np.isnan(norm_value):
+                        raise ct.ControlArgument("Norm computation resulted in NaN.")    
+                    else:
+                        return norm_value   
+
+    # ---------------------------
+    # L-infinity norm computation
+    # ---------------------------
+    elif p == "inf":   
+
+        # Check for cases with infinite norm
+        poles = G.poles()
+        if G.isdtime():  # Discrete time
+            if any(np.isclose(abs(poles), 1.0)):  # Poles on unit circle
+                if print_warning:
+                    warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
+                return float('inf')
+        else:  # Continuous time
+            if any(np.isclose(poles.real, 0.0)):  # Poles on imaginary axis
+                if print_warning:
+                    warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.")
+                return float('inf')
+
+        # Use slycot, if available, to compute (finite) norm
+        if method == 'slycot':
+            return ct.linfnorm(G, tol)[0]
+
+        # Else use scipy
+        else:
+
+            # ------------------   
+            # Discrete time case
+            # ------------------
+            # Use inverse bilinear transformation of discrete time system to s-plane if no poles on |z|=1 or z=0.
+            # Allows us to use test for continuous time systems next.
+            if G.isdtime():
+                Ad = A
+                Bd = B
+                Cd = C
+                Dd = D
+                if any(np.isclose(la.eigvals(Ad), 0.0)):
+                    raise ct.ControlArgument("L-infinity norm computation for discrete time system with pole(s) in z=0 currently not supported unless Slycot installed.")            
+
+                # Inverse bilinear transformation
+                In = np.eye(len(Ad))
+                Adinv = la.inv(Ad+In)
+                A = 2*(Ad-In)@Adinv
+                B = 2*Adinv@Bd
+                C = 2*Cd@Adinv
+                D = Dd - Cd@Adinv@Bd
+
+            # --------------------
+            # Continuous time case
+            # --------------------
+            def _Hamilton_matrix(gamma):
+                """Constructs Hamiltonian matrix. For internal use."""
+                R = Ip*gamma**2 - D.T@D
+                invR = la.inv(R)
+                return np.block([[A+B@invR@D.T@C, B@invR@B.T], [-C.T@(Ip+D@invR@D.T)@C, -(A+B@invR@D.T@C).T]])        
+
+            gaml = la.norm(D,ord=2)    # Lower bound
+            gamu = max(1.0, 2.0*gaml)  # Candidate upper bound
+            Ip = np.eye(len(D))    
+
+            while any(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)):  # Find actual upper bound
+                gamu *= 2.0
+
+            while (gamu-gaml)/gamu > tol:
+                gam = (gamu+gaml)/2.0
+                if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):
+                    gaml = gam
+                else:
+                    gamu = gam
+            return gam
+
+    # ----------------------
+    # Other norm computation
+    # ----------------------
+    else:
+        raise ct.ControlArgument(f"Norm computation for p={p} currently not supported.")           
+
diff --git a/control/tests/sysnorm_test.py b/control/tests/sysnorm_test.py
@@ -0,0 +1,63 @@
+# -*- coding: utf-8 -*-
+"""
+Tests for sysnorm module.
+
+Created on Mon Jan  8 11:31:46 2024
+Author: Henrik Sandberg
+"""
+
+import control as ct
+import numpy as np
+
+def test_norm_1st_order_stable_system():
+    """First-order stable continuous-time system"""
+    s = ct.tf('s')
+    G1 = 1/(s+1)
+    assert np.allclose(ct.norm(G1, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
+    assert np.allclose(ct.norm(G1, p=2), 0.707106781186547) # Comparison to norm computed in MATLAB
+
+    Gd1 = ct.sample_system(G1, 0.1)
+    assert np.allclose(ct.norm(Gd1, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
+    assert np.allclose(ct.norm(Gd1, p=2), 0.223513699524858) # Comparison to norm computed in MATLAB
+
+
+def test_norm_1st_order_unstable_system():
+    """First-order unstable continuous-time system"""
+    s = ct.tf('s')
+    G2 = 1/(1-s)
+    assert np.allclose(ct.norm(G2, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
+    assert ct.norm(G2, p=2) == float('inf') # Comparison to norm computed in MATLAB
+
+    Gd2 = ct.sample_system(G2, 0.1)
+    assert np.allclose(ct.norm(Gd2, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
+    assert ct.norm(Gd2, p=2) == float('inf') # Comparison to norm computed in MATLAB
+
+def test_norm_2nd_order_system_imag_poles():
+    """Second-order continuous-time system with poles on imaginary axis"""
+    s = ct.tf('s')
+    G3 = 1/(s**2+1)
+    assert ct.norm(G3, p='inf') == float('inf') # Comparison to norm computed in MATLAB
+    assert ct.norm(G3, p=2) == float('inf') # Comparison to norm computed in MATLAB
+
+    Gd3 = ct.sample_system(G3, 0.1)
+    assert ct.norm(Gd3, p='inf') == float('inf') # Comparison to norm computed in MATLAB
+    assert ct.norm(Gd3, p=2) == float('inf') # Comparison to norm computed in MATLAB
+
+def test_norm_3rd_order_mimo_system():
+    """Third-order stable MIMO continuous-time system"""
+    A = np.array([[-1.017041847539126,  -0.224182952826418,   0.042538079149249],
+                  [-0.310374015319095,  -0.516461581407780,  -0.119195790221750],
+                  [-1.452723568727942,   1.7995860837102088,  -1.491935830615152]])
+    B = np.array([[0.312858596637428,  -0.164879019209038],
+                  [-0.864879917324456,   0.627707287528727],
+                  [-0.030051296196269,   1.093265669039484]])
+    C = np.array([[1.109273297614398,   0.077359091130425,  -1.113500741486764],
+                  [-0.863652821988714,  -1.214117043615409,  -0.006849328103348]])
+    D = np.zeros((2,2))
+    G4 = ct.ss(A,B,C,D) # Random system generated in MATLAB
+    assert np.allclose(ct.norm(G4, p='inf', tol=1e-9), 4.276759162964244) # Comparison to norm computed in MATLAB
+    assert np.allclose(ct.norm(G4, p=2), 2.237461821810309) # Comparison to norm computed in MATLAB
+
+    Gd4 = ct.sample_system(G4, 0.1)
+    assert np.allclose(ct.norm(Gd4, p='inf', tol=1e-9), 4.276759162964228) # Comparison to norm computed in MATLAB
+    assert np.allclose(ct.norm(Gd4, p=2), 0.707434962289554) # Comparison to norm computed in MATLAB