8000 Merge pull request #971 from henriks76/system-norms · python-control/python-control@e1e33e4 · GitHub
[go: up one dir, main page]
Skip to content

Commit e1e33e4

Browse files
bnavigatorbnavigator
authored
Merge pull request #971 from henriks76/system-norms
System norms 2
2 parents af4f8a7 + bad229a commit e1e33e4

File tree

4 files changed

+372
-0
lines changed

4 files changed

+372
-0
lines changed

‎.gitignore

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,9 @@ TAGS
3131
# Files created by Spyder
3232
.spyproject/
3333

34+
# Files created by or for VS Code (HS, 13 Jan, 2024)
35+
.vscode/
36+
3437
# Environments
3538
.env
3639
.venv

‎control/__init__.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,7 @@
101101
from .config import *
102102
from .sisotool import *
103103
from .passivity import *
104+
from .sysnorm import *
104105

105106
# Exceptions
106107
from .exception import *

‎control/sysnorm.py

Lines changed: 294 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,294 @@
1+
# -*- coding: utf-8 -*-
2+
"""sysnorm.py
3+
4+
Functions for computing system norms.
5+
6+
Routine in this module:
7+
8+
norm
9+
10+
Created on Thu Dec 21 08:06:12 2023
11+
Author: Henrik Sandberg
12+
"""
13+
14+
import numpy as np
15+
import scipy as sp
16+
import numpy.linalg as la
17+
import warnings
18+
19+
import control as ct
20+
21+
__all__ = ['norm']
22+
23+
#------------------------------------------------------------------------------
24+
25+
def _h2norm_slycot(sys, print_warning=True):
26+
"""H2 norm of a linear system. For internal use. Requires Slycot.
27+
28+
See also
29+
--------
30+
``slycot.ab13bd`` : the Slycot routine that does the calculation
31+
https://github.com/python-control/Slycot/issues/199 : Post on issue with ``ab13bf``
32+
"""
33+
34+
try:
35+
from slycot import ab13bd
36+
except ImportError:
37+
ct.ControlSlycot("Can't find slycot module ``ab13bd``!")
38+
39+
try:
40+
from slycot.exceptions import SlycotArithmeticError
41+
except ImportError:
42+
raise ct.ControlSlycot("Can't find slycot class ``SlycotArithmeticError``!")
43+
44+
A, B, C, D = ct.ssdata(ct.ss(sys))
45+
46+
n = A.shape[0]
47+
m = B.shape[1]
48+
p = C.shape[0]
49+
50+
dico = 'C' if sys.isctime() else 'D' # Continuous or discrete time
51+
jobn = 'H' # H2 (and not L2 norm)
52+
53+
if n == 0:
54+
# ab13bd does not accept empty A, B, C
55+
if dico == 'C':
56+
if any(D.flat != 0):
57+
if print_warning:
58+
warnings.warn("System has a direct feedthrough term!", UserWarning)
59+
return float("inf")
60+
else:
61+
return 0.0
62+
elif dico == 'D':
63+
return np.sqrt(D@D.T)
64+
65+
try:
66+
norm = ab13bd(dico, jobn, n, m, p, A, B, C, D)
67+
except SlycotArithmeticError as e:
68+
if e.info == 3:
69+
if print_warning:
70+
warnings.warn("System has pole(s) on the stability boundary!", UserWarning)
71+
return float("inf")
72+
elif e.info == 5:
73+
if print_warning:
74+
warnings.warn("System has a direct feedthrough term!", UserWarning)
75+
return float("inf")
76+
elif e.info == 6:
77+
if print_warning:
78+
warnings.warn("System is unstable!", UserWarning)
79+
return float("inf")
80+
else:
81+
raise e
82+
return norm
83+
84+
#------------------------------------------------------------------------------
85+
86+
def norm(system, p=2, tol=1e-6, print_warning=True, method=None):
87+
"""Computes norm of system.
88+
89+
Parameters
90+
----------
91+
system : LTI (:class:`StateSpace` or :class:`TransferFunction`)
92+
System in continuous or discrete time for which the norm should be computed.
93+
p : int or str
94+
Type of norm to be computed. ``p=2`` gives the H2 norm, and ``p='inf'`` gives the L-infinity norm.
95+
tol : float
96+
Relative tolerance for accuracy of L-infinity norm computation. Ignored
97+
unless p='inf'.
98+
print_warning : bool
99+
Print warning message in case norm value may be uncertain.
100+
method : str, optional
101+
Set the method used for computing the result. Current methods are
102+
'slycot' and 'scipy'. If set to None (default), try 'slycot' first
103+
and then 'scipy'.
104+
105+
Returns
106+
-------
107+
norm_value : float
108+
Norm value of system.
109+
110+
Notes
111+
-----
112+
Does not yet compute the L-infinity norm for discrete time systems with pole(s) in z=0 unless Slycot is used.
113+
114+
Examples
115+
--------
116+
>>> Gc = ct.tf([1], [1, 2, 1])
117+
>>> ct.norm(Gc, 2)
118+
0.5000000000000001
119+
>>> ct.norm(Gc, 'inf', tol=1e-11, method='scipy')
120+
1.000000000007276
121+
"""
122+
123+
if not isinstance(system, (ct.StateSpace, ct.TransferFunction)):
124+
raise TypeError('Parameter ``system``: must be a ``StateSpace`` or ``TransferFunction``')
125+
126+
G = ct.ss(system)
127+
A = G.A
128+
B = G.B
129+
C = G.C
130+
D = G.D
131+
132+
# Decide what method to use
133+
method = ct.mateqn._slycot_or_scipy(method)
134+
135+
# -------------------
136+
# H2 norm computation
137+
# -------------------
138+
if p == 2:
139+
# --------------------
140+
# Continuous time case
141+
# --------------------
142+
if G.isctime():
143+
144+
# Check for cases with infinite norm
145+
poles_real_part = G.poles().real
146+
if any(np.isclose(poles_real_part, 0.0)): # Poles on imaginary axis
147+
if print_warning:
148+
warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.", UserWarning)
149+
return float('inf')
150+
elif any(poles_real_part > 0.0): # System unstable
151+
if print_warning:
152+
warnings.warn("System is unstable!", UserWarning)
153+
return float('inf')
154+
elif any(D.flat != 0): # System has direct feedthrough
155+
if print_warning:
156+
warnings.warn("System has a direct feedthrough term!", UserWarning)
157+
return float('inf')
158+
159+
else:
160+
# Use slycot, if available, to compute (finite) norm
161+
if method == 'slycot':
162+
return _h2norm_slycot(G, print_warning)
163+
164+
# Else use scipy
165+
else:
166+
P = ct.lyap(A, B@B.T, method=method) # Solve for controllability Gramian
167+
168+
# System is stable to reach this point, and P should be positive semi-definite.
169+
# Test next is a precaution in case the Lyapunov equation is ill conditioned.
170+
if any(la.eigvals(P).real < 0.0):
171+
if print_warning:
172+
warnings.warn("There appears to be poles close to the imaginary axis. Norm value may be uncertain.", UserWarning)
173+
return float('inf')
174+
else:
175+
norm_value = np.sqrt(np.trace(C@P@C.T)) # Argument in sqrt should be non-negative
176+
if np.isnan(norm_value):
177+
raise ct.ControlArgument("Norm computation resulted in NaN.")
178+
else:
179+
return norm_value
180+
181+
# ------------------
182+
# Discrete time case
183+
# ------------------
184+
elif G.isdtime():
185+
186+
# Check for cases with infinite norm
187+
poles_abs = abs(G.poles())
188+
if any(np.isclose(poles_abs, 1.0)): # Poles on imaginary axis
189+
if print_warning:
190+
warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.", UserWarning)
191+
return float('inf')
192+
elif any(poles_abs > 1.0): # System unstable
193+
if print_warning:
194+
warnings.warn("System is unstable!", UserWarning)
195+
return float('inf')
196+
197+
else:
198+
# Use slycot, if available, to compute (finite) norm
199+
if method == 'slycot':
200+
return _h2norm_slycot(G, print_warning)
201+
202+
# Else use scipy
203+
else:
204+
P = ct.dlyap(A, B@B.T, method=method)
205+
206+
# System is stable to reach this point, and P should be positive semi-definite.
207+
# Test next is a precaution in case the Lyapunov equation is ill conditioned.
208+
if any(la.eigvals(P).real < 0.0):
209+
if print_warning:
210+
warnings.warn("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.", UserWarning)
211+
return float('inf')
212+
else:
213+
norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T)) # Argument in sqrt should be non-negative
214+
if np.isnan(norm_value):
215+
raise ct.ControlArgument("Norm computation resulted in NaN.")
216+
else:
217+
return norm_value
218+
219+
# ---------------------------
220+
# L-infinity norm computation
221+
# ---------------------------
222+
elif p == "inf":
223+
224+
# Check for cases with infinite norm
225+
poles = G.poles()
226+
if G.isdtime(): # Discrete time
227+
if any(np.isclose(abs(poles), 1.0)): # Poles on unit circle
228+
if print_warning:
229+
warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.", UserWarning)
230+
return float('inf')
231+
else: # Continuous time
232+
if any(np.isclose(poles.real, 0.0)): # Poles on imaginary axis
233+
if print_warning:
234+
warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.", UserWarning)
235+
return float('inf')
236+
237+
# Use slycot, if available, to compute (finite) norm
238+
if method == 'slycot':
239+
return ct.linfnorm(G, tol)[0]
240+
241+
# Else use scipy
242+
else:
243+
244+
# ------------------
245+
# Discrete time case
246+
# ------------------
247+
# Use inverse bilinear transformation of discrete time system to s-plane if no poles on |z|=1 or z=0.
248+
# Allows us to use test for continuous time systems next.
249+
if G.isdtime():
250+
Ad = A
251+
Bd = B
252+
Cd = C
253+
Dd = D
254+
if any(np.isclose(la.eigvals(Ad), 0.0)):
255+
raise ct.ControlArgument("L-infinity norm computation for discrete time system with pole(s) in z=0 currently not supported unless Slycot installed.")
256+
257+
# Inverse bilinear transformation
258+
In = np.eye(len(Ad))
259+
Adinv = la.inv(Ad+In)
260+
A = 2*(Ad-In)@Adinv
261+
B = 2*Adinv@Bd
262+
C = 2*Cd@Adinv
263+
D = Dd - Cd@Adinv@Bd
264+
265+
# --------------------
266+
# Continuous time case
267+
# --------------------
268+
def _Hamilton_matrix(gamma):
269+
10000 """Constructs Hamiltonian matrix. For internal use."""
270+
R = Ip*gamma**2 - D.T@D
271+
invR = la.inv(R)
272+
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]])
273+
274+
gaml = la.norm(D,ord=2) # Lower bound
275+
gamu = max(1.0, 2.0*gaml) # Candidate upper bound
276+
Ip = np.eye(len(D))
277+
278+
while any(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)): # Find actual upper bound
279+
gamu *= 2.0
280+
281+
while (gamu-gaml)/gamu > tol:
282+
gam = (gamu+gaml)/2.0
283+
if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):
284+
gaml = gam
285+
else:
286+
gamu = gam
287+
return gam
288+
289+
# ----------------------
290+
# Other norm computation
291+
# ----------------------
292+
else:
293+
raise ct.ControlArgument(f"Norm computation for p={p} currently not supported.")
294+

‎control/tests/sysnorm_test.py

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
# -*- coding: utf-8 -*-
2+
"""
3+
Tests for sysnorm module.
4+
5+
Created on Mon Jan 8 11:31:46 2024
6+
Author: Henrik Sandberg
7+
"""
8+
9+
import control as ct
10+
import numpy as np
11+
import pytest
12+
13+
14+
def test_norm_1st_order_stable_system():
15+
"""First-order stable continuous-time system"""
16+
s = ct.tf('s')
17+
18+
G1 = 1/(s+1)
19+
assert np.allclose(ct.norm(G1, p='inf'), 1.0) # Comparison to norm computed in MATLAB
20+
assert np.allclose(ct.norm(G1, p=2), 0.707106781186547) # Comparison to norm computed in MATLAB
21+
22+
Gd1 = ct.sample_system(G1, 0.1)
23+
assert np.allclose(ct.norm(Gd1, p='inf'), 1.0) # Comparison to norm computed in MATLAB
24+
assert np.allclose(ct.norm(Gd1, p=2), 0.223513699524858) # Comparison to norm computed in MATLAB
25+
26+
27+
def test_norm_1st_order_unstable_system():
28+
"""First-order unstable continuous-time system"""
29+
s = ct.tf('s')
30+
31+
G2 = 1/(1-s)
32+
assert np.allclose(ct.norm(G2, p='inf'), 1.0) # Comparison to norm computed in MATLAB
33+
with pytest.warns(UserWarning, match="System is unstable!"):
34+
assert ct.norm(G2, p=2) == float('inf') # Comparison to norm computed in MATLAB
35+
36+
Gd2 = ct.sample_system(G2, 0.1)
37+
assert np.allclose(ct.norm(Gd2, p='inf'), 1.0) # Comparison to norm computed in MATLAB
38+
with pytest.warns(UserWarning, match="System is unstable!"):
39+
assert ct.norm(Gd2, p=2) == float('inf') # Comparison to norm computed in MATLAB
40+
41+
def test_norm_2nd_order_system_imag_poles():
42+
"""Second-order continuous-time system with poles on imaginary axis"""
43+
s = ct.tf('s')
44+
45+
G3 = 1/(s**2+1)
46+
with pytest.warns(UserWarning, match="Poles close to, or on, the imaginary axis."):
47+
assert ct.norm(G3, p='inf') == float('inf') # Comparison to norm computed in MATLAB
48+
with pytest.warns(UserWarning, match="Poles close to, or on, the imaginary axis."):
49+
assert ct.norm(G3, p=2) == float('inf') # Comparison to norm computed in MATLAB
50+
51+
Gd3 = ct.sample_system(G3, 0.1)
52+
with pytest.warns(UserWarning, match="Poles close to, or on, the complex unit circle."):
53+
assert ct.norm(Gd3, p='inf') == float('inf') # Comparison to norm computed in MATLAB
54+
with pytest.warns(UserWarning, match="Poles close to, or on, the complex unit circle."):
55+
assert ct.norm(Gd3, p=2) == float('inf') # Comparison to norm computed in MATLAB
56+
57+
def test_norm_3rd_order_mimo_system():
58+
"""Third-order stable MIMO continuous-time system"""
59+
A = np.array([[-1.017041847539126, -0.224182952826418, 0.042538079149249],
60+
[-0.310374015319095, -0.516461581407780, -0.119195790221750],
61+
[-1.452723568727942, 1.7995860837102088, -1.491935830615152]])
62+
B = np.array([[0.312858596637428, -0.164879019209038],
63+
[-0.864879917324456, 0.627707287528727],
64+
[-0.030051296196269, 1.093265669039484]])
65+
C = np.array([[1.109273297614398, 0.077359091130425, -1.113500741486764],
66+
[-0.863652821988714, -1.214117043615409, -0.006849328103348]])
67+
D = np.zeros((2,2))
68+
G4 = ct.ss(A,B,C,D) # Random system generated in MATLAB
69+
assert np.allclose(ct.norm(G4, p='inf'), 4.276759162964244) # Comparison to norm computed in MATLAB
70+
assert np.allclose(ct.norm(G4, p=2), 2.237461821810309) # Comparison to norm computed in MATLAB
71+
72+
Gd4 = ct.sample_system(G4, 0.1)
73+
assert np.allclose(ct.norm(Gd4, p='inf'), 4.276759162964228) # Comparison to norm computed in MATLAB
74+
assert np.allclose(ct.norm(Gd4, p=2), 0.707434962289554) # Comparison to norm computed in MATLAB

0 commit comments

Comments
 (0)
0