- Notifications
You must be signed in to change notification settings - Fork441
System norms#960
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to ourterms of service andprivacy statement. We’ll occasionally send you account related emails.
Already on GitHub?Sign in to your account
Closed
Uh oh!
There was an error while loading.Please reload this page.
Closed
System norms#960
Changes fromall commits
Commits
Show all changes
27 commits Select commitHold shift + click to select a range
6a0e136 New function for LTI system norm computation
a0fbc80 * Updated documentation of function norm
457c623 Update sysnorm.py
adcb486 Merge branch 'python-control:main' into system-norms
henriks76fd076c5 New function for LTI system norm computation
0fe2c57 * Updated documentation of function norm
cdf4bab Update sysnorm.py
34f9537 Do not track changes in VS Code setup.
b419f12 Lowered tolerances in tests.
9ecd594 Added:
14fce67 Merge branch 'system-norms' of https://github.com/henriks76/python-co…
5c38d4f escape latex labels
bnavigator99e56f8 Use Numpy API for pi instead of undocumented scipy.pi
bnavigator2c32913 Update notebooks (no rerun, no output change)
bnavigator85328d9 Merge pull request #965 from bnavigator/fix-examples
bnavigator7ace4bc New function for LTI system norm computation
5103448 * Updated documentation of function norm
4fb252e Update sysnorm.py
108817c Do not track changes in VS Code setup.
6683eb3 Lowered tolerances in tests.
32d38bf Added:
6f810ba Update sysnorm.py
793f0d6 Added:
12ba677 Merge branch 'system-norms' of https://github.com/henriks76/python-co…
7a5af50 Fixed merge error with __all__.
49f7e5f Update sysnorm.py
247cff5 Merge branch 'system-norms' of https://github.com/henriks76/python-co…
File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Uh oh!
There was an error while loading.Please reload this page.
Jump to
Jump to file
Failed to load files.
Loading
Uh oh!
There was an error while loading.Please reload this page.
Diff view
Diff view
There are no files selected for viewing
3 changes: 3 additions & 0 deletions.gitignore
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.Learn more about bidirectional Unicode characters
1 change: 1 addition & 0 deletionscontrol/__init__.py
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -101,6 +101,7 @@ | ||
| from .config import * | ||
| from .sisotool import * | ||
| from .passivity import * | ||
| from .sysnorm import * | ||
bnavigator marked this conversation as resolved. Show resolvedHide resolvedUh oh!There was an error while loading.Please reload this page. | ||
| # Exceptions | ||
| from .exception import * | ||
306 changes: 306 additions & 0 deletionscontrol/sysnorm.py
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others.Learn more. You might want to add a type check similar to ifnotisinstance(sys, (StateSpace,TransferFunction)):raiseTypeError('Parameter ``sys``: must be a ``StateSpace`` or ``TransferFunction``)') befor calling G=ct.ss(sys) | ||
| 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.") | ||
63 changes: 63 additions & 0 deletionscontrol/tests/sysnorm_test.py
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
Oops, something went wrong.
Uh oh!
There was an error while loading.Please reload this page.
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.