8000 * Updated documentation of function norm · python-control/python-control@e8cca34 · GitHub
[go: up one dir, main page]
Skip to content {"props":{"docsUrl":"https://docs.github.com/get-started/accessibility/keyboard-shortcuts"}}

Commit e8cca34

Browse files
author
Henrik Sandberg
committed
* 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.
1 parent a1f27db commit e8cca34

File tree

3 files changed

+194
-23
lines changed

3 files changed

+194
-23
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/sysnorm.py

Lines changed: 128 additions & 23 deletions
+
Norm value of system (float) or None if computation could not be completed.
Original file line numberDiff line numberDiff line change
@@ -1,75 +1,180 @@
11
# -*- coding: utf-8 -*-
2-
"""
3-
Created on Thu Dec 21 08:06:12 2023
2+
"""sysnorm.py
3+
4+
Functions for computing system norms.
5+
6+
Routine in this module:
7+
8+
norm
49
5-
@author: hsan
10+
Created on Thu Dec 21 08:06:12 2023
11+
Author: Henrik Sandberg
612
"""
713

814
import numpy as np
915
import numpy.linalg as la
1016

1117
import control as ct
1218

19+
__all__ = ['norm']
20+
1321
#------------------------------------------------------------------------------
1422

15-
def norm(system, p=2, tol=1e-6):
16-
"""Computes H_2 (p=2) or L_infinity (p="inf", tolerance tol) norm of system."""
23+
def norm(system, p=2, tol=1e-6, print_warning=True):
24+
"""Computes norm of system.
25+
26+
Parameters
27+
----------
28+
system : LTI (:class:`StateSpace` or :class:`TransferFunction`)
29+
System in continuous or discrete time for which the norm should be computed.
30+
p : int or str
31+
Type of norm to be computed. p=2 gives the H_2 norm, and p='inf' gives the L_infinity norm.
32+
tol : float
33+
Relative tolerance for accuracy of L_infinity norm computation. Ignored
34+
unless p='inf'.
35+
print_warning : bool
36+
Print warning message in case norm value may be uncertain.
37+
38+
Returns
39+
-------
40+
norm_value : float or NoneType
41
42+
43+
Notes
44+
-----
45+
Does not yet compute the L_infinity norm for discrete time systems with pole(s) in z=0.
46+
47+
Examples
48+
--------
49+
>>> Gc = ct.tf([1], [1, 2, 1])
50+
>>> ct.norm(Gc,2)
51+
0.5000000000000001
52+
>>> ct.norm(Gc,'inf',tol=1e-10)
53+
1.0000000000582077
54+
"""
1755
G = ct.ss(system)
1856
A = G.A
1957
B = G.B
2058
C = G.C
2159
D = G.D
2260

23-
if p == 2: # H_2-norm
61+
#
62+
# H_2-norm computation
63+
#
64+
if p == 2:
65+
# Continuous time case
2466
if G.isctime():
25-
if (D != 0).any() or any(G.poles().real >= 0):
67+
poles_real_part = G.poles().real
68+
if any(np.isclose(poles_real_part, 0.0)):
69+
if print_warning:
70+
print("Warning: Poles close to, or on, the imaginary axis. Norm value may be uncertain.")
71+
return float('inf')
72+
elif (D != 0).any() or any(poles_real_part > 0.0): # System unstable or has direct feedthrough?
2673
return float('inf')
2774
else:
28-
P = ct.lyap(A, B@B.T)
29-
return np.sqrt(np.trace(C@P@C.T))
75+
try:
76+
P = ct.lyap(A, B@B.T)
77+
except Exception as e:
78+
print(f"An error occurred solving the continuous time Lyapunov equation: {e}")
79+
return None
80+
81+
# System is stable to reach this point, and P should be positive semi-definite.
82+
# Test next is a precaution in case the Lyapunov equation is ill conditioned.
83+
if any(la.eigvals(P) < 0.0):
84+
if print_warning:
85+
print("Warning: There appears to be poles close to the imaginary axis. Norm value may be uncertain.")
86+
return float('inf')
87+
else:
88+
norm_value = np.sqrt(np.trace(C@P@C.T)) # Argument in sqrt should be non-negative
89+
if np.isnan(norm_value):
90+
print("Unknown error. Norm computation resulted in NaN.")
91+
return None
92+
else:
93+
return norm_value
94+
95+
# Discrete time case
3096
elif G.isdtime():
31-
if any(abs(G.poles()) >= 1):
97+
poles_abs = abs(G.poles())
98+
if any(np.isclose(poles_abs, 1.0)):
99+
if print_warning:
100+
print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
101+
return float('inf')
102+
elif any(poles_abs > 1.0): # System unstable?
32103
return float('inf')
33104
else:
34-
P = ct.dlyap(A, B@B.T)
35-
return np.sqrt(np.trace(C@P@C.T + D@D.T))
36-
37-
elif p == "inf": # L_infinity-norm
38-
def Hamilton_matrix(gamma):
39-
"""Constructs Hamiltonian matrix."""
105+
try:
106+
P = ct.dlyap(A, B@B.T)
107+
except Exception as e:
108+
print(f"An error occurred solving the discrete time Lyapunov equation: {e}")
109+
return None
110+
111+
# System is stable to reach this point, and P should be positive semi-definite.
112+
# Test next is a precaution in case the Lyapunov equation is ill conditioned.
113+
if any(la.eigvals(P) < 0.0):
114+
if print_warning:
115+
print("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.")
116+
return float('inf')
117+
else:
118+
norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T)) # Argument in sqrt should be non-negative
119+
if np.isnan(norm_value):
120+
print("Unknown error. Norm computation resulted in NaN.")
121+
return None
122+
else:
123+
return norm_value
124+
#
125+
# L_infinity-norm computation
126+
#
127+
elif p == "inf":
128+
def _Hamilton_matrix(gamma):
129+
"""Constructs Hamiltonian matrix. For internal use."""
40130
R = Ip*gamma**2 - D.T@D
41131
invR = la.inv(R)
42132
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]])
43133

44-
if G.isdtime(): # Bilinear transformation to s-plane
134+
# Discrete time case
135+
# Use inverse bilinear transformation of discrete time system to s-plane if no poles on |z|=1 or z=0.
136+
# Allows us to use test for continuous time systems next.
137+
if G.isdtime():
45138
Ad = A
46139
Bd = B
47140
Cd = C
48141
Dd = D
142+
if any(np.isclose(abs(la.eigvals(Ad)), 1.0)):
143+
if print_warning:
144+
print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
145+
return float('inf')
146+
elif any(np.isclose(la.eigvals(Ad), 0.0)):
147+
print("L_infinity norm computation for discrete time system with pole(s) at z=0 currently not supported.")
148+
return None
149+
150+
# Inverse bilinear transformation
49151
In = np.eye(len(Ad))
50152
Adinv = la.inv(Ad+In)
51153
A = 2*(Ad-In)@Adinv
52154
B = 2*Adinv@Bd
53155
C = 2*Cd@Adinv
54156
D = Dd - Cd@Adinv@Bd
55-
157+
158+
# Continus time case
56159
if any(np.isclose(la.eigvals(A).real, 0.0)):
160+
if print_warning:
161+
print("Warning: Poles close to, or on, imaginary axis. Norm value may be uncertain.")
57162
return float('inf')
58163

59-
gaml = la.norm(D,ord=2) # Lower bound
60-
gamu = max(1.0, 2.0*gaml) # Candidate upper bound
164+
gaml = la.norm(D,ord=2) # Lower bound
165+
gamu = max(1.0, 2.0*gaml) # Candidate upper bound
61166
Ip = np.eye(len(D))
62167

63-
while any(np.isclose(la.eigvals(Hamilton_matrix(gamu)).real, 0.0)): # Find an upper bound
168+
while any(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)): # Find actual upper bound
64169
gamu *= 2.0
65170

66171
while (gamu-gaml)/gamu > tol:
67172
gam = (gamu+gaml)/2.0
68-
if any(np.isclose(la.eigvals(Hamilton_matrix(gam)).real, 0.0)):
173+
if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):
69174
gaml = gam
70175
else:
71176
gamu = gam
72177
return gam
73178
else:
74-
# Norm computation only supported for p=2 and p='inf'
179+
print(f"Norm computation for p={p} currently not supported.")
75180
return None

‎control/tests/sysnorm_test.py

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
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+
12+
def test_norm_1st_order_stable_system():
13+
"""First-order stable continuous-time system"""
14+
s = ct.tf('s')
15+
G1 = 1/(s+1)
16+
assert np.allclose(ct.norm(G1, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
17+
assert np.allclose(ct.norm(G1, p=2), 0.707106781186547) # Comparison to norm computed in MATLAB
18+
19+
Gd1 = ct.sample_system(G1, 0.1)
20+
assert np.allclose(ct.norm(Gd1, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
21+
assert np.allclose(ct.norm(Gd1, p=2), 0.223513699524858) # Comparison to norm computed in MATLAB
22+
23+
24+
def test_norm_1st_order_unstable_system():
25+
"""First-order unstable continuous-time system"""
26+
s = ct.tf('s')
27+
G2 = 1/(1-s)
28+
assert np.allclose(ct.norm(G2, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
29+
assert ct.norm(G2, p=2) == float('inf') # Comparison to norm computed in MATLAB
30+
31+
Gd2 = ct.sample_system(G2, 0.1)
32+
assert np.allclose(ct.norm(Gd2, p='inf', tol=1e-9), 1.0) # Comparison to norm computed in MATLAB
33+
assert ct.norm(Gd2, p=2) == float('inf') # Comparison to norm computed in MATLAB
34+
35+
def test_norm_2nd_order_system_imag_poles():
36+
"""Second-order continuous-time system with poles on imaginary axis"""
37+
s = ct.tf('s')
38+
G3 = 1/(s**2+1)
39+
assert ct.norm(G3, p='inf') == float('inf') # Comparison to norm computed in MATLAB
40+
assert ct.norm(G3, p=2) == float('inf') # Comparison to norm computed in MATLAB
41+
42+
Gd3 = ct.sample_system(G3, 0.1)
43+
assert ct.norm(Gd3, p='inf') == float('inf') # Comparison to norm computed in MATLAB
44+
assert ct.norm(Gd3, p=2) == float('inf') # Comparison to norm computed in MATLAB
45+
46+
def test_norm_3rd_order_mimo_system():
47+
"""Third-order stable MIMO continuous-time system"""
48+
A = np.array([[-1.017041847539126, -0.224182952826418, 0.042538079149249],
49+
[-0.310374015319095, -0.516461581407780, -0.119195790221750],
50+
[-1.452723568727942, 1.7995860837102088, -1.491935830615152]])
51+
B = np.array([[0.312858596637428, -0.164879019209038],
52+
[-0.864879917324456, 0.627707287528727],
53+
[-0.030051296196269, 1.093265669039484]])
54+
C = np.array([[1.109273297614398, 0.077359091130425, -1.113500741486764],
55+
[-0.863652821988714, -1.214117043615409, -0.006849328103348]])
56+
D = np.zeros((2,2))
57+
G4 = ct.ss(A,B,C,D) # Random system generated in MATLAB
58+
assert np.allclose(ct.norm(G4, p='inf', tol=1e-9), 4.276759162964244) # Comparison to norm computed in MATLAB
59+
assert np.allclose(ct.norm(G4, p=2), 2.237461821810309) # Comparison to norm computed in MATLAB
60+
61+
Gd4 = ct.sample_system(G4, 0.1)
62+
assert np.allclose(ct.norm(Gd4, p='inf', tol=1e-9), 4.276759162964228) # Comparison to norm computed in MATLAB
63+
assert np.allclose(ct.norm(Gd4, p=2), 0.707434962289554) # Comparison to norm computed in MATLAB

0 commit comments

Comments
 (0)
0