Skip to content

Commit f60c83a

Browse files
author
Henrik Sandberg
committed
Added:
* Use of warnings package. * Use routine statesp.linfnorm when slycot installed. * New routine internal _h2norm_slycot when slycot is installed.
1 parent e8cca34 commit f60c83a

File tree

1 file changed

+195
-90
lines changed

1 file changed

+195
-90
lines changed

‎control/sysnorm.py

Lines changed: 195 additions & 90 deletions
Original file line numberDiff line numberDiff line change
@@ -12,28 +12,93 @@
1212
"""
1313

1414
import numpy as np
15+
import scipy as sp
1516
import numpy.linalg as la
17+
import warnings
1618

1719
import control as ct
1820

1921
__all__ = ['norm']
2022

2123
#------------------------------------------------------------------------------
2224

23-
def norm(system, p=2, tol=1e-6, print_warning=True):
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!")
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!")
71+
return float("inf")
72+
elif e.info == 5:
73+
if print_warning:
74+
warnings.warn("System has a direct feedthrough term!")
75+
return float("inf")
76+
elif e.info == 6:
77+
if print_warning:
78+
warnings.warn("System is unstable!")
79+
return float("inf")
80+
else:
81+
raise e
82+
return norm
83+
84+
#------------------------------------------------------------------------------
85+
86+
def norm(system, p=2, tol=1e-10, print_warning=True, use_slycot=True):
2487
"""Computes norm of system.
2588
2689
Parameters
2790
----------
2891
system : LTI (:class:`StateSpace` or :class:`TransferFunction`)
2992
System in continuous or discrete time for which the norm should be computed.
3093
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.
94+
Type of norm to be computed. p=2 gives the H2 norm, and p='inf' gives the L-infinity norm.
3295
tol : float
33-
Relative tolerance for accuracy of L_infinity norm computation. Ignored
96+
Relative tolerance for accuracy of L-infinity norm computation. Ignored
3497
unless p='inf'.
3598
print_warning : bool
3699
Print warning message in case norm value may be uncertain.
100+
use_slycot : bool
101+
Use Slycot routines if available.
37102
38103
Returns
39104
-------
@@ -42,7 +107,7 @@ def norm(system, p=2, tol=1e-6, print_warning=True):
42107
43108
Notes
44109
-----
45-
Does not yet compute the L_infinity norm for discrete time systems with pole(s) in z=0.
110+
Does not yet compute the L-infinity norm for discrete time systems with pole(s) in z=0 unless Slycot is used.
46111
47112
Examples
48113
--------
@@ -58,123 +123,163 @@ def norm(system, p=2, tol=1e-6, print_warning=True):
58123
C = G.C
59124
D = G.D
60125

61-
#
62-
# H_2-norm computation
63-
#
126+
# -------------------
127+
# H2 norm computation
128+
# -------------------
64129
if p == 2:
130+
# --------------------
65131
# Continuous time case
132+
# --------------------
66133
if G.isctime():
134+
135+
# Check for cases with infinite norm
67136
poles_real_part = G.poles().real
68-
if any(np.isclose(poles_real_part, 0.0)):
137+
if any(np.isclose(poles_real_part, 0.0)): # Poles on imaginary axis
69138
if print_warning:
70-
print("Warning: Poles close to, or on, the imaginary axis. Norm value may be uncertain.")
139+
warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.")
71140
return float('inf')
72-
elif (D != 0).any() or any(poles_real_part > 0.0): # System unstable or has direct feedthrough?
141+
elif any(poles_real_part > 0.0): # System unstable
142+
if print_warning:
143+
warnings.warn("System is unstable!")
73144
return float('inf')
74-
else:
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
145+
elif any(D.flat != 0): # System has direct feedthrough
146+
if print_warning:
147+
warnings.warn("System has a direct feedthrough term!")
148+
return float('inf')
149+
150+
else:
151+
# Use slycot, if available, to compute (finite) norm
152+
if ct.slycot_check() and use_slycot:
153+
return _h2norm_slycot(G, print_warning)
80154

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')
155+
# Else use scipy
87156
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
157+
P = ct.lyap(A, B@B.T) # Solve for controllability Gramian
158+
159+
# System is stable to reach this point, and P should be positive semi-definite.
160+
# Test next is a precaution in case the Lyapunov equation is ill conditioned.
161+
if any(la.eigvals(P).real < 0.0):
162+
if print_warning:
163+
warnings.warn("There appears to be poles close to the imaginary axis. Norm value may be uncertain.")
164+
return float('inf')
92165
else:
93-
return norm_value
166+
norm_value = np.sqrt(np.trace(C@P@C.T)) # Argument in sqrt should be non-negative
167+
if np.isnan(norm_value):
168+
raise ct.ControlArgument("Norm computation resulted in NaN.")
169+
else:
170+
return norm_value
94171

172+
# ------------------
95173
# Discrete time case
174+
# ------------------
96175
elif G.isdtime():
176+
177+
# Check for cases with infinite norm
97178
poles_abs = abs(G.poles())
98-
if any(np.isclose(poles_abs, 1.0)):
179+
if any(np.isclose(poles_abs, 1.0)): # Poles on imaginary axis
99180
if print_warning:
100-
print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
181+
warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
101182
return float('inf')
102-
elif any(poles_abs > 1.0): # System unstable?
183+
elif any(poles_abs > 1.0): # System unstable
184+
if print_warning:
185+
warnings.warn("System is unstable!")
103186
return float('inf')
187+
104188
else:
105-
try:
189+
# Use slycot, if available, to compute (finite) norm
190+
if ct.slycot_check() and use_slycot:
191+
return _h2norm_slycot(G, print_warning)
192+
193+
# Else use scipy
194+
else:
106195
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
110196

111197
# System is stable to reach this point, and P should be positive semi-definite.
112198
# Test next is a precaution in case the Lyapunov equation is ill conditioned.
113-
if any(la.eigvals(P) < 0.0):
199+
if any(la.eigvals(P).real < 0.0):
114200
if print_warning:
115-
print("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.")
201+
warnings.warn("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.")
116202
return float('inf')
117203
else:
118204
norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T)) # Argument in sqrt should be non-negative
119205
if np.isnan(norm_value):
120-
print("Unknown error. Norm computation resulted in NaN.")
121-
return None
206+
raise ct.ControlArgument("Norm computation resulted in NaN.")
122207
else:
123-
return norm_value
124-
#
125-
# L_infinity-norm computation
126-
#
208+
return norm_value
209+
210+
# ---------------------------
211+
# L-infinity norm computation
212+
# ---------------------------
127213
elif p == "inf":
128-
def _Hamilton_matrix(gamma):
129-
"""Constructs Hamiltonian matrix. For internal use."""
130-
R = Ip*gamma**2 - D.T@D
131-
invR = la.inv(R)
132-
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]])
133-
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():
138-
Ad = A
139-
Bd = B
140-
Cd = C
141-
Dd = D
142-
if any(np.isclose(abs(la.eigvals(Ad)), 1.0)):
214+
215+
# Check for cases with infinite norm
216+
poles = G.poles()
217+
if G.isdtime(): # Discrete time
218+
if any(np.isclose(abs(poles), 1.0)): # Poles on unit circle
143219
if print_warning:
144-
print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
220+
warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
221+
return float('inf')
222+
else: # Continuous time
223+
if any(np.isclose(poles.real, 0.0)): # Poles on imaginary axis
224+
if print_warning:
225+
warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.")
145226
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
151-
In = np.eye(len(Ad))
152-
Adinv = la.inv(Ad+In)
153-
A = 2*(Ad-In)@Adinv
154-
B = 2*Adinv@Bd
155-
C = 2*Cd@Adinv
156-
D = Dd - Cd@Adinv@Bd
157-
158-
# Continus time case
159-
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.")
162-
return float('inf')
163227

164-
gaml = la.norm(D,ord=2) # Lower bound
165-
gamu = max(1.0, 2.0*gaml) # Candidate upper bound
166-
Ip = np.eye(len(D))
167-
168-
while any(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)): # Find actual upper bound
169-
gamu *= 2.0
228+
# Use slycot, if available, to compute (finite) norm
229+
if ct.slycot_check() and use_slycot:
230+
return ct.linfnorm(G, tol)[0]
170231

171-
while (gamu-gaml)/gamu > tol:
172-
gam = (gamu+gaml)/2.0
173-
if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):
174-
gaml = gam
175-
else:
176-
gamu = gam
177-
return gam
232+
# Else use scipy
233+
else:
234+
235+
# ------------------
236+
# Discrete time case
237+
# ------------------
238+
# Use inverse bilinear transformation of discrete time system to s-plane if no poles on |z|=1 or z=0.
239+
# Allows us to use test for continuous time systems next.
240+
if G.isdtime():
241+
Ad = A
242+
Bd = B
243+
Cd = C
244+
Dd = D
245+
if any(np.isclose(la.eigvals(Ad), 0.0)):
246+
raise ct.ControlArgument("L-infinity norm computation for discrete time system with pole(s) in z=0 currently not supported unless Slycot installed.")
247+
248+
# Inverse bilinear transformation
249+
In = np.eye(len(Ad))
250+
Adinv = la.inv(Ad+In)
251+
A = 2*(Ad-In)@Adinv
252+
B = 2*Adinv@Bd
253+
C = 2*Cd@Adinv
254+
D = Dd - Cd@Adinv@Bd
255+
256+
# --------------------
257+
# Continuous time case
258+
# --------------------
259+
def _Hamilton_matrix(gamma):
260+
"""Constructs Hamiltonian matrix. For internal use."""
261+
R = Ip*gamma**2 - D.T@D
262+
invR = la.inv(R)
263+
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]])
264+
265+
gaml = la.norm(D,ord=2) # Lower bound
266+
gamu = max(1.0, 2.0*gaml) # Candidate upper bound
267+
Ip = np.eye(len(D))
268+
269+
while any(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)): # Find actual upper bound
270+
gamu *= 2.0
271+
272+
while (gamu-gaml)/gamu > tol:
273+
gam = (gamu+gaml)/2.0
274+
if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):
275+
gaml = gam
276+
else:
277+
gamu = gam
278+
return gam
279+
280+
# ----------------------
281+
# Other norm computation
282+
# ----------------------
178283
else:
179-
print(f"Norm computation for p={p} currently not supported.")
180-
return None
284+
raise ct.ControlArgument(f"Norm computation for p={p} currently not supported.")
285+

0 commit comments

Comments
 (0)