GV

GeezerVison

The World Wide Web, restored to 800 pixels of glory.

Home|About|Gaming cabinet|Archives|How it works|Starter sites

On this site

  • Sign in
  • ActionsCI/CD and automation
  • CodespacesCloud development environments
  • All Features
  • Documentation
  • Training
  • Certifications
  • Build with Copilot
  • Customize Copilot
  • Integrations
  • Trending repos
  • GitHub Sponsors
  • Developer articles
  • Enterprises

This view

  • Original page
  • New address

GeezerVison is showing github.com — System norms 2 by henriks76 · Pull Request #971 · python-control/python-control

System norms 2 #971

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.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bnavigator merged 6 commits into python-control:main from henriks76:system-norms
Feb 18, 2024

Uh oh!

There was an error while loading. Please reload this page.
Merged

System norms 2 #971

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
a1f27db
New function for LTI system norm computation
Jan 7, 2024
e8cca34
* Updated documentation of function norm
Jan 8, 2024
f60c83a
Added:
Jan 21, 2024
99e49d9
Added in sysnorm.py:
Jan 28, 2024
b39710f
* In sysnorm: Warnings now of type UserWarning
Feb 11, 2024
bad229a
* sysnorm: Changed default tolerance in norm to 1e-6.
Feb 18, 2024
File filter
Filter by extension
Clear filters
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 3 .gitignore
Show comments
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions 1 control/__init__.py
Show comments
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,7 @@
from .config import *
from .sisotool import *
from .passivity import *
from .sysnorm import *

# Exceptions
from .exception import *
Expand Down
294 changes: 294 additions & 0 deletions 294 control/sysnorm.py
Show comments
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,294 @@
# -*- 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 _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!", UserWarning)
return float("inf")
else:
return 0.0
elif dico == 'D':
return np.sqrt([email protected])

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 norm(system, p=2, tol=1e-6, 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 = 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:
P = ct.lyap(A, [email protected], 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.", UserWarning)
return float('inf')
else:
norm_value = np.sqrt(np.trace(C@[email protected])) # 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.", 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, [email protected], 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.", UserWarning)
return float('inf')
else:
norm_value = np.sqrt(np.trace(C@[email protected] + [email protected])) # 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.", 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@[email protected]@C, B@[email protected]], [-C.T@(Ip+D@[email protected])@C, -(A+B@[email protected]@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.")

74 changes: 74 additions & 0 deletions 74 control/tests/sysnorm_test.py
Show comments
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
# -*- 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
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'), 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'), 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'), 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'), 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)
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)
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"""
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'), 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'), 4.276759162964228) # Comparison to norm computed in MATLAB
assert np.allclose(ct.norm(Gd4, p=2), 0.707434962289554) # Comparison to norm computed in MATLAB
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.
You are visitor 01402503 — GeezerVison does not store the pages you look up.
Best viewed with Netscape Navigator 4.7 or Internet Explorer 6 at 800×600. No cookies required.