Skip to content

Commit 26d5310

Browse files
murrayrmmurrayrm
committed
minor code reformatting of sysnorm.py + unit tests
1 parent 18f7a4c commit 26d5310

File tree

2 files changed

+317
-86
lines changed

2 files changed

+317
-86
lines changed

‎control/sysnorm.py‎

Lines changed: 126 additions & 86 deletions
Original file line numberDiff line numberDiff line change
@@ -20,58 +20,64 @@
2020

2121
__all__ = ['norm']
2222

23-
#------------------------------------------------------------------------------
2423

2524
def _h2norm_slycot(sys, print_warning=True):
2625
"""H2 norm of a linear system. For internal use. Requires Slycot.
2726
2827
See also
2928
--------
3029
``slycot.ab13bd`` : the Slycot routine that does the calculation
31-
https://github.com/python-control/Slycot/issues/199 : Post on issue with ``ab13bf``
30+
https://github.com/python-control/Slycot/issues/199 : Post on issue
31+
with ``ab13bf``
32+
3233
"""
33-
34+
3435
try:
3536
from slycot import ab13bd
3637
except ImportError:
3738
ct.ControlSlycot("Can't find slycot module ``ab13bd``!")
3839

3940
try:
4041
from slycot.exceptions import SlycotArithmeticError
41-
except ImportError:
42-
raise ct.ControlSlycot("Can't find slycot class ``SlycotArithmeticError``!")
42+
except ImportError:
43+
raise ct.ControlSlycot(
44+
"Can't find slycot class ``SlycotArithmeticError``!")
4345

4446
A, B, C, D = ct.ssdata(ct.ss(sys))
4547

4648
n = A.shape[0]
4749
m = B.shape[1]
4850
p = C.shape[0]
49-
51+
5052
dico = 'C' if sys.isctime() else 'D' # Continuous or discrete time
51-
jobn = 'H' # H2 (and not L2 norm)
53+
jobn = 'H' # H2 (and not L2 norm)
5254

5355
if n == 0:
5456
# ab13bd does not accept empty A, B, C
5557
if dico == 'C':
5658
if any(D.flat != 0):
5759
if print_warning:
58-
warnings.warn("System has a direct feedthrough term!", UserWarning)
60+
warnings.warn(
61+
"System has a direct feedthrough term!", UserWarning)
5962
return float("inf")
6063
else:
6164
return 0.0
6265
elif dico == 'D':
6366
return np.sqrt(D@D.T)
64-
67+
6568
try:
6669
norm = ab13bd(dico, jobn, n, m, p, A, B, C, D)
6770
except SlycotArithmeticError as e:
6871
if e.info == 3:
6972
if print_warning:
70-
warnings.warn("System has pole(s) on the stability boundary!", UserWarning)
73+
warnings.warn(
74+
"System has pole(s) on the stability boundary!",
75+
UserWarning)
7176
return float("inf")
7277
elif e.info == 5:
7378
if print_warning:
74-
warnings.warn("System has a direct feedthrough term!", UserWarning)
79+
warnings.warn(
80+
"System has a direct feedthrough term!", UserWarning)
7581
return float("inf")
7682
elif e.info == 6:
7783
if print_warning:
@@ -81,214 +87,248 @@ def _h2norm_slycot(sys, print_warning=True):
8187
raise e
8288
return norm
8389

84-
#------------------------------------------------------------------------------
8590

8691
def norm(system, p=2, tol=1e-6, print_warning=True, method=None):
8792
"""Computes norm of system.
88-
93+
8994
Parameters
9095
----------
9196
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
97+
System in continuous or discrete time for which the norm should be
98+
computed.
99+
p : int or str, optional
100+
Type of norm to be computed. ``p=2`` (default) gives the H2 norm, and
101+
``p='inf'`` gives the L-infinity norm.
102+
tol : float, optional
103+
Relative tolerance for accuracy of L-infinity norm
104+
computation. Ignored unless ``p='inf'``. Default value, 1e-6.
105+
print_warning : bool, optional
99106
Print warning message in case norm value may be uncertain.
100107
method : str, optional
101108
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-
109+
``'slycot'`` and ``'scipy'``. If set to ``None`` (default), try
110+
``'slycot'`` first and then ``'scipy'``.
111+
105112
Returns
106113
-------
107114
norm_value : float
108115
Norm value of system.
109-
116+
110117
Notes
111118
-----
112-
Does not yet compute the L-infinity norm for discrete time systems with pole(s) in z=0 unless Slycot is used.
113-
119+
Does not yet compute the L-infinity norm for discrete time systems with
120+
pole(s) in z=0 unless Slycot is used.
121+
114122
Examples
115123
--------
116124
>>> Gc = ct.tf([1], [1, 2, 1])
117125
>>> round(ct.norm(Gc, 2), 3)
118126
0.5
119127
>>> round(ct.norm(Gc, 'inf', tol=1e-5, method='scipy'), 3)
120128
np.float64(1.0)
129+
121130
"""
122-
131+
123132
if not isinstance(system, (ct.StateSpace, ct.TransferFunction)):
124-
raise TypeError('Parameter ``system``: must be a ``StateSpace`` or ``TransferFunction``')
125-
133+
raise TypeError(
134+
"Parameter ``system``: must be a ``StateSpace`` or "
135+
"``TransferFunction``")
136+
126137
G = ct.ss(system)
127138
A = G.A
128139
B = G.B
129140
C = G.C
130141
D = G.D
131-
142+
132143
# Decide what method to use
133144
method = ct.mateqn._slycot_or_scipy(method)
134-
145+
135146
# -------------------
136147
# H2 norm computation
137148
# -------------------
138-
if p == 2:
149+
if p == 2:
139150
# --------------------
140151
# Continuous time case
141152
# --------------------
142153
if G.isctime():
143-
154+
144155
# Check for cases with infinite norm
145-
poles_real_part = G.poles().real
156+
poles_real_part = G.poles().real
146157
if any(np.isclose(poles_real_part, 0.0)): # Poles on imaginary axis
147158
if print_warning:
148-
warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.", UserWarning)
159+
warnings.warn(
160+
"Poles close to, or on, the imaginary axis. "
161+
"Norm value may be uncertain.", UserWarning)
149162
return float('inf')
150163
elif any(poles_real_part > 0.0): # System unstable
151164
if print_warning:
152165
warnings.warn("System is unstable!", UserWarning)
153166
return float('inf')
154167
elif any(D.flat != 0): # System has direct feedthrough
155168
if print_warning:
156-
warnings.warn("System has a direct feedthrough term!", UserWarning)
157-
return float('inf')
158-
159-
else:
169+
warnings.warn(
170+
"System has a direct feedthrough term!", UserWarning)
171+
return float('inf')
172+
173+
else:
160174
# Use slycot, if available, to compute (finite) norm
161175
if method == 'slycot':
162-
return _h2norm_slycot(G, print_warning)
163-
164-
# Else use scipy
176+
return _h2norm_slycot(G, print_warning)
177+
178+
# Else use scipy
165179
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):
180+
# Solve for controllability Gramian
181+
P = ct.lyap(A, B@B.T, method=method)
182+
183+
# System is stable to reach this point, and P should
184+
# be positive semi-definite. Test next is a precaution
185+
# in case the Lyapunov equation is ill conditioned.
186+
if any(la.eigvals(P).real < 0.0):
171187
if print_warning:
172-
warnings.warn("There appears to be poles close to the imaginary axis. Norm value may be uncertain.", UserWarning)
188+
warnings.warn(
189+
"There appears to be poles close to the "
190+
"imaginary axis. Norm value may be uncertain.",
191+
UserWarning)
173192
return float('inf')
174193
else:
175-
norm_value = np.sqrt(np.trace(C@P@C.T)) # Argument in sqrt should be non-negative
194+
# Argument in sqrt should be non-negative
195+
norm_value = np.sqrt(np.trace(C@P@C.T))
176196
if np.isnan(norm_value):
177-
raise ct.ControlArgument("Norm computation resulted in NaN.")
197+
raise ct.ControlArgument(
198+
"Norm computation resulted in NaN.")
178199
else:
179200
return norm_value
180-
201+
181202
# ------------------
182203
# Discrete time case
183204
# ------------------
184205
elif G.isdtime():
185-
206+
186207
# Check for cases with infinite norm
187208
poles_abs = abs(G.poles())
188209
if any(np.isclose(poles_abs, 1.0)): # Poles on imaginary axis
189210
if print_warning:
190-
warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.", UserWarning)
211+
warnings.warn(
212+
"Poles close to, or on, the complex unit circle. "
213+
"Norm value may be uncertain.", UserWarning)
191214
return float('inf')
192215
elif any(poles_abs > 1.0): # System unstable
193216
if print_warning:
194217
warnings.warn("System is unstable!", UserWarning)
195218
return float('inf')
196-
219+
197220
else:
198221
# Use slycot, if available, to compute (finite) norm
199222
if method == 'slycot':
200-
return _h2norm_slycot(G, print_warning)
201-
202-
# Else use scipy
223+
return _h2norm_slycot(G, print_warning)
224+
225+
# Else use scipy
203226
else:
204227
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.
228+
229+
# System is stable to reach this point, and P should be
230+
# positive semi-definite. Test next is a precaution in
231+
# case the Lyapunov equation is ill conditioned.
208232
if any(la.eigvals(P).real < 0.0):
209233
if print_warning:
210-
warnings.warn("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.", UserWarning)
234+
warnings.warn(
235+
"Poles close to the complex unit circle. "
236+
"Norm value may be uncertain.", UserWarning)
211237
return float('inf')
212238
else:
213-
norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T)) # Argument in sqrt should be non-negative
239+
# Argument in sqrt should be non-negative
240+
norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T))
214241
if np.isnan(norm_value):
215-
raise ct.ControlArgument("Norm computation resulted in NaN.")
242+
raise ct.ControlArgument(
243+
"Norm computation resulted in NaN.")
216244
else:
217-
return norm_value
218-
245+
return norm_value
246+
219247
# ---------------------------
220248
# L-infinity norm computation
221249
# ---------------------------
222-
elif p == "inf":
223-
250+
elif p == "inf":
251+
224252
# Check for cases with infinite norm
225253
poles = G.poles()
226254
if G.isdtime(): # Discrete time
227255
if any(np.isclose(abs(poles), 1.0)): # Poles on unit circle
228256
if print_warning:
229-
warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.", UserWarning)
257+
warnings.warn(
258+
"Poles close to, or on, the complex unit circle. "
259+
"Norm value may be uncertain.", UserWarning)
230260
return float('inf')
231261
else: # Continuous time
232262
if any(np.isclose(poles.real, 0.0)): # Poles on imaginary axis
233263
if print_warning:
234-
warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.", UserWarning)
264+
warnings.warn(
265+
"Poles close to, or on, the imaginary axis. "
266+
"Norm value may be uncertain.", UserWarning)
235267
return float('inf')
236-
268+
237269
# Use slycot, if available, to compute (finite) norm
238270
if method == 'slycot':
239271
return ct.linfnorm(G, tol)[0]
240-
272+
241273
# Else use scipy
242274
else:
243-
244-
# ------------------
275+
276+
# ------------------
245277
# Discrete time case
246278
# ------------------
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.
279+
# Use inverse bilinear transformation of discrete time system
280+
# to s-plane if no poles on |z|=1 or z=0. Allows us to use test
281+
# for continuous time systems next.
249282
if G.isdtime():
250283
Ad = A
251284
Bd = B
252285
Cd = C
253286
Dd = D
254287
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-
288+
raise ct.ControlArgument(
289+
"L-infinity norm computation for discrete time "
290+
"system with pole(s) in z=0 currently not "
291+
"supported unless Slycot installed.")
292+
257293
# Inverse bilinear transformation
258294
In = np.eye(len(Ad))
259295
Adinv = la.inv(Ad+In)
260296
A = 2*(Ad-In)@Adinv
261297
B = 2*Adinv@Bd
262298
C = 2*Cd@Adinv
263299
D = Dd - Cd@Adinv@Bd
264-
300+
265301
# --------------------
266302
# Continuous time case
267303
# --------------------
268304
def _Hamilton_matrix(gamma):
269305
"""Constructs Hamiltonian matrix. For internal use."""
270306
R = Ip*gamma**2 - D.T@D
271307
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]])
308+
return np.block(
309+
[[A+B@invR@D.T@C, B@invR@B.T],
310+
[-C.T@(Ip+D@invR@D.T)@C, -(A+B@invR@D.T@C).T]])
273311

274312
gaml = la.norm(D,ord=2) # Lower bound
275313
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
314+
Ip = np.eye(len(D))
315+
316+
# Find actual upper bound
317+
while any(np.isclose(la.eigvals(
318+
_Hamilton_matrix(gamu)).real, 0.0)):
279319
gamu *= 2.0
280-
320+
281321
while (gamu-gaml)/gamu > tol:
282322
gam = (gamu+gaml)/2.0
283323
if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):
284324
gaml = gam
285325
else:
286326
gamu = gam
287327
return gam
288-
328+
289329
# ----------------------
290330
# Other norm computation
291331
# ----------------------
292332
else:
293-
raise ct.ControlArgument(f"Norm computation for p={p} currently not supported.")
294-
333+
raise ct.ControlArgument(
334+
f"Norm computation for p={p} currently not supported.")

0 commit comments

Comments
 (0)