# sysnorm.py - functions for computing system norms # # Initial author: Henrik Sandberg # Creation date: 21 Dec 2023 """Functions for computing system norms.""" import warnings import numpy as np import numpy.linalg as la import control as ct __all__ = ['system_norm', 'norm'] #------------------------------------------------------------------------------ def _h2norm_slycot(sys, print_warning=True): """H2 norm of a linear system. For internal use. Requires Slycot. See Also -------- slycot.ab13bd """ # See: https://github.com/python-control/Slycot/issues/199 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!", UserWarning) 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!", UserWarning) return float("inf") elif e.info == 5: if print_warning: warnings.warn( "System has a direct feedthrough term!", UserWarning) return float("inf") elif e.info == 6: if print_warning: warnings.warn("System is unstable!", UserWarning) return float("inf") else: raise e return norm #------------------------------------------------------------------------------ def system_norm(system, p=2, tol=1e-6, print_warning=True, method=None): """Computes the input/output norm of system. Parameters ---------- system : LTI (`StateSpace` or `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) at the origin unless Slycot is used. Examples -------- >>> Gc = ct.tf([1], [1, 2, 1]) >>> round(ct.norm(Gc, 2), 3) 0.5 >>> round(ct.norm(Gc, 'inf', tol=1e-5, method='scipy'), 3) np.float64(1.0) """ 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 = ct.mateqn._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.", UserWarning) return float('inf') elif any(poles_real_part > 0.0): # System unstable if print_warning: warnings.warn("System is unstable!", UserWarning) return float('inf') elif any(D.flat != 0): # System has direct feedthrough if print_warning: warnings.warn( "System has a direct feedthrough term!", UserWarning) 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: # Solve for controllability Gramian P = ct.lyap(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( "There appears to be poles close to the " "imaginary axis. Norm value may be uncertain.", UserWarning) return float('inf') else: # Argument in sqrt should be non-negative norm_value = np.sqrt(np.trace(C@P@C.T)) 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.", UserWarning) return float('inf') elif any(poles_abs > 1.0): # System unstable if print_warning: warnings.warn("System is unstable!", UserWarning) 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( "There appears to be poles close to the complex " "unit circle. Norm value may be uncertain.", UserWarning) return float('inf') else: # Argument in sqrt should be non-negative norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T)) 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.", UserWarning) 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.", UserWarning) 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.") norm = system_norm