From a1f27dbdbccde9ac898bc662f647a488984da5b0 Mon Sep 17 00:00:00 2001 From: Henrik Sandberg Date: Sun, 7 Jan 2024 21:53:32 +0100 Subject: [PATCH 1/6] New function for LTI system norm computation --- control/__init__.py | 1 + control/sysnorm.py | 75 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+) create mode 100644 control/sysnorm.py diff --git a/control/__init__.py b/control/__init__.py index 120d16325..5a9e05e95 100644 --- 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 new file mode 100644 index 000000000..5caae0918 --- /dev/null +++ b/control/sysnorm.py @@ -0,0 +1,75 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Dec 21 08:06:12 2023 + +@author: hsan +""" + +import numpy as np +import numpy.linalg as la + +import control as ct + +#------------------------------------------------------------------------------ + +def norm(system, p=2, tol=1e-6): + """Computes H_2 (p=2) or L_infinity (p="inf", tolerance tol) norm of system.""" + 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 to s-plane + Ad = A + Bd = B + Cd = C + Dd = D + 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: + # Norm computation only supported for p=2 and p='inf' + return None From e8cca3431b33181c33c643738caabe0be98c5115 Mon Sep 17 00:00:00 2001 From: Henrik Sandberg Date: Mon, 8 Jan 2024 13:51:33 +0100 Subject: [PATCH 2/6] * Updated documentation of function norm * Added control/tests/sysnorm_test.py * Added additional tests and warning messages for systems with poles close to stability boundary * Added __all__ * Do not track changes in VS Code setup. --- .gitignore | 3 + control/sysnorm.py | 151 ++++++++++++++++++++++++++++------ control/tests/sysnorm_test.py | 63 ++++++++++++++ 3 files changed, 194 insertions(+), 23 deletions(-) create mode 100644 control/tests/sysnorm_test.py diff --git a/.gitignore b/.gitignore index 1b10a3585..4a6aa3cc0 100644 --- 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/sysnorm.py b/control/sysnorm.py index 5caae0918..2065f8721 100644 --- a/control/sysnorm.py +++ b/control/sysnorm.py @@ -1,8 +1,14 @@ # -*- coding: utf-8 -*- -""" -Created on Thu Dec 21 08:06:12 2023 +"""sysnorm.py + +Functions for computing system norms. + +Routine in this module: + +norm -@author: hsan +Created on Thu Dec 21 08:06:12 2023 +Author: Henrik Sandberg """ import numpy as np @@ -10,66 +16,165 @@ import control as ct +__all__ = ['norm'] + #------------------------------------------------------------------------------ -def norm(system, p=2, tol=1e-6): - """Computes H_2 (p=2) or L_infinity (p="inf", tolerance tol) norm of system.""" +def norm(system, p=2, tol=1e-6, print_warning=True): + """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'. + print_warning : bool + Print warning message in case norm value may be uncertain. + + Returns + ------- + norm_value : float or NoneType + Norm value of system (float) or None if computation could not be completed. + + 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 + # + # H_2-norm computation + # + if p == 2: + # Continuous time case if G.isctime(): - if (D != 0).any() or any(G.poles().real >= 0): + poles_real_part = G.poles().real + if any(np.isclose(poles_real_part, 0.0)): + if print_warning: + print("Warning: Poles close to, or on, the imaginary axis. Norm value may be uncertain.") + return float('inf') + elif (D != 0).any() or any(poles_real_part > 0.0): # System unstable or has direct feedthrough? return float('inf') else: - P = ct.lyap(A, B@B.T) - return np.sqrt(np.trace(C@P@C.T)) + try: + P = ct.lyap(A, B@B.T) + except Exception as e: + print(f"An error occurred solving the continuous time Lyapunov equation: {e}") + return None + + # 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) < 0.0): + if print_warning: + print("Warning: 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): + print("Unknown error. Norm computation resulted in NaN.") + return None + else: + return norm_value + + # Discrete time case elif G.isdtime(): - if any(abs(G.poles()) >= 1): + poles_abs = abs(G.poles()) + if any(np.isclose(poles_abs, 1.0)): + if print_warning: + print("Warning: 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? 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.""" + try: + P = ct.dlyap(A, B@B.T) + except Exception as e: + print(f"An error occurred solving the discrete time Lyapunov equation: {e}") + return None + + # 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) < 0.0): + if print_warning: + print("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): + print("Unknown error. Norm computation resulted in NaN.") + return None + else: + return norm_value + # + # L_infinity-norm computation + # + elif p == "inf": + 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]]) - if G.isdtime(): # Bilinear transformation to s-plane + # 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(abs(la.eigvals(Ad)), 1.0)): + if print_warning: + print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.") + 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 + + # 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 - + + # Continus time case if any(np.isclose(la.eigvals(A).real, 0.0)): + if print_warning: + print("Warning: Poles close to, or on, imaginary axis. Norm value may be uncertain.") return float('inf') - gaml = la.norm(D,ord=2) # Lower bound - gamu = max(1.0, 2.0*gaml) # Candidate upper bound + 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 + 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)): + if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)): gaml = gam else: gamu = gam return gam else: - # Norm computation only supported for p=2 and p='inf' + print(f"Norm computation for p={p} currently not supported.") return None diff --git a/control/tests/sysnorm_test.py b/control/tests/sysnorm_test.py new file mode 100644 index 000000000..917e98d04 --- /dev/null +++ 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 From f60c83ad11de1cfab20b8bfbf37ff7c2a3a68d5f Mon Sep 17 00:00:00 2001 From: Henrik Sandberg Date: Sun, 21 Jan 2024 19:29:32 +0100 Subject: [PATCH 3/6] Added: * Use of warnings package. * Use routine statesp.linfnorm when slycot installed. * New routine internal _h2norm_slycot when slycot is installed. --- control/sysnorm.py | 285 +++++++++++++++++++++++++++++++-------------- 1 file changed, 195 insertions(+), 90 deletions(-) diff --git a/control/sysnorm.py b/control/sysnorm.py index 2065f8721..a25ef305f 100644 --- a/control/sysnorm.py +++ b/control/sysnorm.py @@ -12,7 +12,9 @@ """ import numpy as np +import scipy as sp import numpy.linalg as la +import warnings import control as ct @@ -20,7 +22,68 @@ #------------------------------------------------------------------------------ -def norm(system, p=2, tol=1e-6, print_warning=True): +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, use_slycot=True): """Computes norm of system. Parameters @@ -28,12 +91,14 @@ def norm(system, p=2, tol=1e-6, print_warning=True): 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. + 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 + 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. + use_slycot : bool + Use Slycot routines if available. Returns ------- @@ -42,7 +107,7 @@ def norm(system, p=2, tol=1e-6, print_warning=True): Notes ----- - Does not yet compute the L_infinity norm for discrete time systems with pole(s) in z=0. + Does not yet compute the L-infinity norm for discrete time systems with pole(s) in z=0 unless Slycot is used. Examples -------- @@ -58,123 +123,163 @@ def norm(system, p=2, tol=1e-6, print_warning=True): C = G.C D = G.D - # - # H_2-norm computation - # + # ------------------- + # 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)): + if any(np.isclose(poles_real_part, 0.0)): # Poles on imaginary axis if print_warning: - print("Warning: Poles close to, or on, the imaginary axis. Norm value may be uncertain.") + warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.") return float('inf') - elif (D != 0).any() or any(poles_real_part > 0.0): # System unstable or has direct feedthrough? + elif any(poles_real_part > 0.0): # System unstable + if print_warning: + warnings.warn("System is unstable!") return float('inf') - else: - try: - P = ct.lyap(A, B@B.T) - except Exception as e: - print(f"An error occurred solving the continuous time Lyapunov equation: {e}") - return None + 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 ct.slycot_check() and use_slycot: + return _h2norm_slycot(G, print_warning) - # 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) < 0.0): - if print_warning: - print("Warning: There appears to be poles close to the imaginary axis. Norm value may be uncertain.") - return float('inf') + # Else use scipy else: - norm_value = np.sqrt(np.trace(C@P@C.T)) # Argument in sqrt should be non-negative - if np.isnan(norm_value): - print("Unknown error. Norm computation resulted in NaN.") - return None + P = ct.lyap(A, B@B.T) # 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: - return norm_value + 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)): + if any(np.isclose(poles_abs, 1.0)): # Poles on imaginary axis if print_warning: - print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.") + 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? + elif any(poles_abs > 1.0): # System unstable + if print_warning: + warnings.warn("System is unstable!") return float('inf') + else: - try: + # Use slycot, if available, to compute (finite) norm + if ct.slycot_check() and use_slycot: + return _h2norm_slycot(G, print_warning) + + # Else use scipy + else: P = ct.dlyap(A, B@B.T) - except Exception as e: - print(f"An error occurred solving the discrete time Lyapunov equation: {e}") - return None # 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) < 0.0): + if any(la.eigvals(P).real < 0.0): if print_warning: - print("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.") + 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): - print("Unknown error. Norm computation resulted in NaN.") - return None + raise ct.ControlArgument("Norm computation resulted in NaN.") else: - return norm_value - # - # L_infinity-norm computation - # + return norm_value + + # --------------------------- + # L-infinity norm computation + # --------------------------- elif p == "inf": - 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]]) - - # 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(abs(la.eigvals(Ad)), 1.0)): + + # 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: - print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.") + 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') - 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 - - # 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 - - # Continus time case - if any(np.isclose(la.eigvals(A).real, 0.0)): - if print_warning: - print("Warning: Poles close to, or on, imaginary axis. Norm value may be uncertain.") - 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 actual upper bound - gamu *= 2.0 + # Use slycot, if available, to compute (finite) norm + if ct.slycot_check() and use_slycot: + return ct.linfnorm(G, tol)[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 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: - print(f"Norm computation for p={p} currently not supported.") - return None + raise ct.ControlArgument(f"Norm computation for p={p} currently not supported.") + From 99e49d966083b836a4237cdee4ea04de14ad5ff5 Mon Sep 17 00:00:00 2001 From: Henrik Sandberg Date: Sun, 28 Jan 2024 15:57:32 +0100 Subject: [PATCH 4/6] Added in sysnorm.py: * type check when calling ct.norm * metod argument in ct.norm (slycot or scipy) using function _slycot_or_scipy from ct.mateqn Change in sysnorm_test.py * set print_warning=False --- control/sysnorm.py | 43 +++++++++++++++++++++-------------- control/tests/sysnorm_test.py | 32 +++++++++++++------------- 2 files changed, 42 insertions(+), 33 deletions(-) diff --git a/control/sysnorm.py b/control/sysnorm.py index a25ef305f..f3e117dee 100644 --- a/control/sysnorm.py +++ b/control/sysnorm.py @@ -27,19 +27,19 @@ def _h2norm_slycot(sys, print_warning=True): See also -------- - slycot.ab13bd : the Slycot routine that does the calculation - https://github.com/python-control/Slycot/issues/199 : Post on issue with ab13bf + ``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'!") + 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'!") + raise ct.ControlSlycot("Can't find slycot class ``SlycotArithmeticError``!") A, B, C, D = ct.ssdata(ct.ss(sys)) @@ -83,7 +83,7 @@ def _h2norm_slycot(sys, print_warning=True): #------------------------------------------------------------------------------ -def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True): +def norm(system, p=2, tol=1e-10, print_warning=True, method=None): """Computes norm of system. Parameters @@ -97,13 +97,15 @@ def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True): unless p='inf'. print_warning : bool Print warning message in case norm value may be uncertain. - use_slycot : bool - Use Slycot routines if available. + 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 or NoneType - Norm value of system (float) or None if computation could not be completed. + norm_value : float + Norm value of system. Notes ----- @@ -112,17 +114,24 @@ def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True): Examples -------- >>> Gc = ct.tf([1], [1, 2, 1]) - >>> ct.norm(Gc,2) + >>> ct.norm(Gc, 2) 0.5000000000000001 - >>> ct.norm(Gc,'inf',tol=1e-10) - 1.0000000000582077 + >>> 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 = ct.mateqn._slycot_or_scipy(method) + # ------------------- # H2 norm computation # ------------------- @@ -149,12 +158,12 @@ def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True): else: # Use slycot, if available, to compute (finite) norm - if ct.slycot_check() and use_slycot: + if method == 'slycot': return _h2norm_slycot(G, print_warning) # Else use scipy else: - P = ct.lyap(A, B@B.T) # Solve for controllability Gramian + 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. @@ -187,12 +196,12 @@ def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True): else: # Use slycot, if available, to compute (finite) norm - if ct.slycot_check() and use_slycot: + if method == 'slycot': return _h2norm_slycot(G, print_warning) # Else use scipy else: - P = ct.dlyap(A, B@B.T) + 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. @@ -226,7 +235,7 @@ def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True): return float('inf') # Use slycot, if available, to compute (finite) norm - if ct.slycot_check() and use_slycot: + if method == 'slycot': return ct.linfnorm(G, tol)[0] # Else use scipy diff --git a/control/tests/sysnorm_test.py b/control/tests/sysnorm_test.py index 917e98d04..74e8de09d 100644 --- a/control/tests/sysnorm_test.py +++ b/control/tests/sysnorm_test.py @@ -13,35 +13,35 @@ 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 + assert np.allclose(ct.norm(G1, p='inf', tol=1e-9, print_warning=False), 1.0) # Comparison to norm computed in MATLAB + assert np.allclose(ct.norm(G1, p=2, print_warning=False), 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 + assert np.allclose(ct.norm(Gd1, p='inf', tol=1e-9, print_warning=False), 1.0) # Comparison to norm computed in MATLAB + assert np.allclose(ct.norm(Gd1, p=2, print_warning=False), 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 + assert np.allclose(ct.norm(G2, p='inf', tol=1e-9, print_warning=False), 1.0) # Comparison to norm computed in MATLAB + assert ct.norm(G2, p=2, print_warning=False) == 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 + assert np.allclose(ct.norm(Gd2, p='inf', tol=1e-9, print_warning=False), 1.0) # Comparison to norm computed in MATLAB + assert ct.norm(Gd2, p=2, print_warning=False) == 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 + assert ct.norm(G3, p='inf', print_warning=False) == float('inf') # Comparison to norm computed in MATLAB + assert ct.norm(G3, p=2, print_warning=False) == 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 + assert ct.norm(Gd3, p='inf', print_warning=False) == float('inf') # Comparison to norm computed in MATLAB + assert ct.norm(Gd3, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB def test_norm_3rd_order_mimo_system(): """Third-order stable MIMO continuous-time system""" @@ -55,9 +55,9 @@ def test_norm_3rd_order_mimo_system(): [-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 + assert np.allclose(ct.norm(G4, p='inf', tol=1e-9, print_warning=False), 4.276759162964244) # Comparison to norm computed in MATLAB + assert np.allclose(ct.norm(G4, p=2, print_warning=False), 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 + assert np.allclose(ct.norm(Gd4, p='inf', tol=1e-9, print_warning=False), 4.276759162964228) # Comparison to norm computed in MATLAB + assert np.allclose(ct.norm(Gd4, p=2, print_warning=False), 0.707434962289554) # Comparison to norm computed in MATLAB From b39710fcfdc20f1a665c1e9669bfc4bb92763e50 Mon Sep 17 00:00:00 2001 From: Henrik Sandberg Date: Sun, 11 Feb 2024 13:30:45 +0100 Subject: [PATCH 5/6] * In sysnorm: Warnings now of type UserWarning * In sysnorm_test: Use pytest.warns context manager --- control/sysnorm.py | 28 +++++++++++------------ control/tests/sysnorm_test.py | 43 ++++++++++++++++++++++------------- 2 files changed, 41 insertions(+), 30 deletions(-) diff --git a/control/sysnorm.py b/control/sysnorm.py index f3e117dee..54e0055d4 100644 --- a/control/sysnorm.py +++ b/control/sysnorm.py @@ -55,7 +55,7 @@ def _h2norm_slycot(sys, print_warning=True): if dico == 'C': if any(D.flat != 0): if print_warning: - warnings.warn("System has a direct feedthrough term!") + warnings.warn("System has a direct feedthrough term!", UserWarning) return float("inf") else: return 0.0 @@ -67,15 +67,15 @@ def _h2norm_slycot(sys, print_warning=True): except SlycotArithmeticError as e: if e.info == 3: if print_warning: - warnings.warn("System has pole(s) on the stability boundary!") + 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!") + warnings.warn("System has a direct feedthrough term!", UserWarning) return float("inf") elif e.info == 6: if print_warning: - warnings.warn("System is unstable!") + warnings.warn("System is unstable!", UserWarning) return float("inf") else: raise e @@ -91,7 +91,7 @@ def norm(system, p=2, tol=1e-10, print_warning=True, method=None): 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. + 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'. @@ -145,15 +145,15 @@ def norm(system, p=2, tol=1e-10, print_warning=True, method=None): 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.") + 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!") + 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!") + warnings.warn("System has a direct feedthrough term!", UserWarning) return float('inf') else: @@ -169,7 +169,7 @@ def norm(system, p=2, tol=1e-10, print_warning=True, method=None): # 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.") + warnings.warn("There appears to be poles close to the imaginary axis. Norm value may be uncertain.", UserWarning) return float('inf') else: norm_value = np.sqrt(np.trace(C@P@C.T)) # Argument in sqrt should be non-negative @@ -187,11 +187,11 @@ def norm(system, p=2, tol=1e-10, print_warning=True, method=None): 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.") + 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!") + warnings.warn("System is unstable!", UserWarning) return float('inf') else: @@ -207,7 +207,7 @@ def norm(system, p=2, tol=1e-10, print_warning=True, method=None): # 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.") + warnings.warn("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.", UserWarning) return float('inf') else: norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T)) # Argument in sqrt should be non-negative @@ -226,12 +226,12 @@ def norm(system, p=2, tol=1e-10, print_warning=True, method=None): 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.") + 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.") + 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 diff --git a/control/tests/sysnorm_test.py b/control/tests/sysnorm_test.py index 74e8de09d..e7cfc88a9 100644 --- a/control/tests/sysnorm_test.py +++ b/control/tests/sysnorm_test.py @@ -8,40 +8,51 @@ import control as ct import numpy as np +import pytest + 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, print_warning=False), 1.0) # Comparison to norm computed in MATLAB - assert np.allclose(ct.norm(G1, p=2, print_warning=False), 0.707106781186547) # Comparison to norm computed in MATLAB + 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, print_warning=False), 1.0) # Comparison to norm computed in MATLAB - assert np.allclose(ct.norm(Gd1, p=2, print_warning=False), 0.223513699524858) # Comparison to norm computed in MATLAB + 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, print_warning=False), 1.0) # Comparison to norm computed in MATLAB - assert ct.norm(G2, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB + assert np.allclose(ct.norm(G2, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB + with pytest.warns(UserWarning, match="System is unstable!"): + 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, print_warning=False), 1.0) # Comparison to norm computed in MATLAB - assert ct.norm(Gd2, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB + assert np.allclose(ct.norm(Gd2, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB + with pytest.warns(UserWarning, match="System is unstable!"): + 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', print_warning=False) == float('inf') # Comparison to norm computed in MATLAB - assert ct.norm(G3, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB + with pytest.warns(UserWarning, match="Poles close to, or on, the imaginary axis."): + assert ct.norm(G3, p='inf') == float('inf') # Comparison to norm computed in MATLAB + with pytest.warns(UserWarning, match="Poles close to, or on, the imaginary axis."): + 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', print_warning=False) == float('inf') # Comparison to norm computed in MATLAB - assert ct.norm(Gd3, p=2, print_warning=False) == float('inf') # Comparison to norm computed in MATLAB + with pytest.warns(UserWarning, match="Poles close to, or on, the complex unit circle."): + assert ct.norm(Gd3, p='inf') == float('inf') # Comparison to norm computed in MATLAB + with pytest.warns(UserWarning, match="Poles close to, or on, the complex unit circle."): + 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""" @@ -55,9 +66,9 @@ def test_norm_3rd_order_mimo_system(): [-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, print_warning=False), 4.276759162964244) # Comparison to norm computed in MATLAB - assert np.allclose(ct.norm(G4, p=2, print_warning=False), 2.237461821810309) # Comparison to norm computed 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, print_warning=False), 4.276759162964228) # Comparison to norm computed in MATLAB - assert np.allclose(ct.norm(Gd4, p=2, print_warning=False), 0.707434962289554) # Comparison to norm computed in MATLAB + 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 From bad229a6006ae93b57a76246dd03318509ffe7a6 Mon Sep 17 00:00:00 2001 From: Henrik Sandberg Date: Sun, 18 Feb 2024 17:14:32 +0100 Subject: [PATCH 6/6] * sysnorm: Changed default tolerance in norm to 1e-6. * sysnorm_test: Use default tolerance in the tests. --- control/sysnorm.py | 2 +- control/tests/sysnorm_test.py | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/control/sysnorm.py b/control/sysnorm.py index 54e0055d4..49873dbaf 100644 --- a/control/sysnorm.py +++ b/control/sysnorm.py @@ -83,7 +83,7 @@ def _h2norm_slycot(sys, print_warning=True): #------------------------------------------------------------------------------ -def norm(system, p=2, tol=1e-10, print_warning=True, method=None): +def norm(system, p=2, tol=1e-6, print_warning=True, method=None): """Computes norm of system. Parameters diff --git a/control/tests/sysnorm_test.py b/control/tests/sysnorm_test.py index e7cfc88a9..68edad230 100644 --- a/control/tests/sysnorm_test.py +++ b/control/tests/sysnorm_test.py @@ -16,11 +16,11 @@ def test_norm_1st_order_stable_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='inf'), 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='inf'), 1.0) # Comparison to norm computed in MATLAB assert np.allclose(ct.norm(Gd1, p=2), 0.223513699524858) # Comparison to norm computed in MATLAB @@ -29,12 +29,12 @@ def test_norm_1st_order_unstable_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 np.allclose(ct.norm(G2, p='inf'), 1.0) # Comparison to norm computed in MATLAB with pytest.warns(UserWarning, match="System is unstable!"): 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 np.allclose(ct.norm(Gd2, p='inf'), 1.0) # Comparison to norm computed in MATLAB with pytest.warns(UserWarning, match="System is unstable!"): assert ct.norm(Gd2, p=2) == float('inf') # Comparison to norm computed in MATLAB @@ -66,9 +66,9 @@ def test_norm_3rd_order_mimo_system(): [-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='inf'), 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='inf'), 4.276759162964228) # Comparison to norm computed in MATLAB assert np.allclose(ct.norm(Gd4, p=2), 0.707434962289554) # Comparison to norm computed in MATLAB