Skip to content

Commit 457c623

Browse files
author
Henrik Sandberg
committed
Update sysnorm.py
* Added additional tests and warning messages for systems with poles close to stability boundary * Added __all__ * Added more comments
1 parent a0fbc80 commit 457c623

File tree

1 file changed

+88
-22
lines changed

1 file changed

+88
-22
lines changed

‎control/sysnorm.py‎

Lines changed: 88 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,9 @@
33
44
Functions for computing system norms.
55
6-
Routines in this module:
6+
Routine in this module:
77
8-
norm()
8+
norm
99
1010
Created on Thu Dec 21 08:06:12 2023
1111
Author: Henrik Sandberg
@@ -16,9 +16,11 @@
1616

1717
import control as ct
1818

19+
__all__ = ['norm']
20+
1921
#------------------------------------------------------------------------------
2022

21-
def norm(system, p=2, tol=1e-6):
23+
def norm(system, p=2, tol=1e-6, print_warning=True):
2224
"""Computes norm of system.
2325
2426
Parameters
@@ -30,11 +32,13 @@ def norm(system, p=2, tol=1e-6):
3032
tol : float
3133
Relative tolerance for accuracy of L_infinity norm computation. Ignored
3234
unless p='inf'.
35+
print_warning : bool
36+
Print warning message in case norm value may be uncertain.
3337
3438
Returns
3539
-------
36-
norm : float
37-
Norm of system
40+
norm_value : float or NoneType
41+
Norm value of system (float) or None if computation could not be completed.
3842
3943
Notes
4044
-----
@@ -54,52 +58,114 @@ def norm(system, p=2, tol=1e-6):
5458
C = G.C
5559
D = G.D
5660

57-
if p == 2: # H_2-norm
61+
#
62+
# H_2-norm computation
63+
#
64+
if p == 2:
65+
# Continuous time case
5866
if G.isctime():
59-
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?
6073
return float('inf')
6174
else:
62-
P = ct.lyap(A, B@B.T)
63-
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
6496
elif G.isdtime():
65-
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?
66103
return float('inf')
67104
else:
68-
P = ct.dlyap(A, B@B.T)
69-
return np.sqrt(np.trace(C@P@C.T + D@D.T))
70-
71-
elif p == "inf": # L_infinity-norm
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":
72128
def _Hamilton_matrix(gamma):
73-
"""Constructs Hamiltonian matrix."""
129+
"""Constructs Hamiltonian matrix. For internal use."""
74130
R = Ip*gamma**2 - D.T@D
75131
invR = la.inv(R)
76132
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]])
77133

78-
if G.isdtime(): # Bilinear transformation of discrete time system to s-plane if no poles at |z|=1 or z=0
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():
79138
Ad = A
80139
Bd = B
81140
Cd = C
82141
Dd = D
83142
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.")
84145
return float('inf')
85146
elif any(np.isclose(la.eigvals(Ad), 0.0)):
86-
print("L_infinity norm computation for discrete time system with pole(s) at z = 0 currently not supported.")
147+
print("L_infinity norm computation for discrete time system with pole(s) at z=0 currently not supported.")
87148
return None
149+
150+
# Inverse bilinear transformation
88151
In = np.eye(len(Ad))
89152
Adinv = la.inv(Ad+In)
90153
A = 2*(Ad-In)@Adinv
91154
B = 2*Adinv@Bd
92155
C = 2*Cd@Adinv
93156
D = Dd - Cd@Adinv@Bd
94-
157+
158+
# Continus time case
95159
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.")
96162
return float('inf')
97163

98-
gaml = la.norm(D,ord=2) # Lower bound
99-
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
100166
Ip = np.eye(len(D))
101167

102-
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
103169
gamu *= 2.0
104170

105171
while (gamu-gaml)/gamu > tol:
@@ -110,5 +176,5 @@ def _Hamilton_matrix(gamma):
110176
gamu = gam
111177
return gam
112178
else:
113-
print("Norm computation for p =", p, "currently not supported.")
179+
print(f"Norm computation for p={p} currently not supported.")
114180
return None

0 commit comments

Comments
 (0)