henriks76 · Jan 7, 2024 · Jan 8, 2024 · Jan 9, 2024 · Jan 21, 2024 · Jan 7, 2024
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,114 @@
+# -*- coding: utf-8 -*-
+"""sysnorm.py
+
+Functions for computing system norms.
+
+Routines in this module:
+
+norm()
+
+Created on Thu Dec 21 08:06:12 2023
+Author: Henrik Sandberg
+"""
+
+import numpy as np
+import numpy.linalg as la
+
+import control as ct
+
+#------------------------------------------------------------------------------
+
+def norm(system, p=2, tol=1e-6):
+    """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 H_2 norm, and p='inf' gives the L_infinity norm.
+    tol : float
+        Relative tolerance for accuracy of L_infinity norm computation. Ignored
+        unless p='inf'.
+
+    Returns
+    -------
+    norm : float
+        Norm of system
+
+    Notes
+    -----
+    Does not yet compute the L_infinity norm for discrete time systems with pole(s) in z=0.
+
+    Examples
+    --------
+    >>> Gc = ct.tf([1], [1, 2, 1])
+    >>> ct.norm(Gc,2)
+    0.5000000000000001
+    >>> ct.norm(Gc,'inf',tol=1e-10)
+    1.0000000000582077
+    """
+    G = ct.ss(system)
+    A = G.A
+    B = G.B
+    C = G.C
+    D = G.D
+
+    if p == 2:  # H_2-norm
+        if G.isctime():
+            if (D != 0).any() or any(G.poles().real >= 0):
+                return float('inf')
+            else:
+                P = ct.lyap(A, B@B.T)
+                return np.sqrt(np.trace(C@P@C.T))
+        elif G.isdtime():
+            if any(abs(G.poles()) >= 1):
+                return float('inf')
+            else:
+                P = ct.dlyap(A, B@B.T)
+                return np.sqrt(np.trace(C@P@C.T + D@D.T))
+
+    elif p == "inf":    # L_infinity-norm
+        def _Hamilton_matrix(gamma):
+            """Constructs Hamiltonian matrix."""
+            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]])    
+
+        if G.isdtime(): # Bilinear transformation of discrete time system to s-plane if no poles at |z|=1 or z=0
+            Ad = A
+            Bd = B
+            Cd = C
+            Dd = D
+            if any(np.isclose(abs(la.eigvals(Ad)), 1.0)):
+                return float('inf')
+            elif any(np.isclose(la.eigvals(Ad), 0.0)):
+                print("L_infinity norm computation for discrete time system with pole(s) at z = 0 currently not supported.")            
+                return None
+            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
+
+        if any(np.isclose(la.eigvals(A).real, 0.0)):
+            return float('inf')
+
+        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 an 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
+    else:
+        print("Norm computation for p =", p, "currently not supported.")           
+        return None
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, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB
+    assert np.allclose(ct.norm(G1, p=2), 0.707106781186547, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB
+    assert np.allclose(ct.norm(Gd1, p=2), 0.223513699524858, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB
+    assert np.allclose(ct.norm(G4, p=2), 2.237461821810309, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB
+    assert np.allclose(ct.norm(Gd4, p=2), 0.707434962289554, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB