Skip to content

Commit a0fbc80

Browse files
author
Henrik Sandberg
committed
* Updated documentation of function norm
* Added control/tests/sysnorm_test.py
1 parent 6a0e136 commit a0fbc80

File tree

2 files changed

+111
-9
lines changed

2 files changed

+111
-9
lines changed

‎control/sysnorm.py‎

Lines changed: 48 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,14 @@
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.
45
5-
@author: hsan
6+
Routines in this module:
7+
8+
norm()
9+
10+
Created on Thu Dec 21 08:06:12 2023
11+
Author: Henrik Sandberg
612
"""
713

814
import numpy as np
@@ -13,7 +19,35 @@
1319
#------------------------------------------------------------------------------
1420

1521
def norm(system, p=2, tol=1e-6):
16-
"""Computes H_2 (p=2) or L_infinity (p="inf", tolerance tol) norm of system."""
22+
"""Computes norm of system.
23+
24+
Parameters
25+
----------
26+
system : LTI (:class:`StateSpace` or :class:`TransferFunction`)
27+
System in continuous or discrete time for which the norm should be computed.
28+
p : int or str
29+
Type of norm to be computed. p=2 gives the H_2 norm, and p='inf' gives the L_infinity norm.
30+
tol : float
31+
Relative tolerance for accuracy of L_infinity norm computation. Ignored
32+
unless p='inf'.
33+
34+
Returns
35+
-------
36+
norm : float
37+
Norm of system
38+
39+
Notes
40+
-----
41+
Does not yet compute the L_infinity norm for discrete time systems with pole(s) in z=0.
42+
43+
Examples
44+
--------
45+
>>> Gc = ct.tf([1], [1, 2, 1])
46+
>>> ct.norm(Gc,2)
47+
0.5000000000000001
48+
>>> ct.norm(Gc,'inf',tol=1e-10)
49+
1.0000000000582077
50+
"""
1751
G = ct.ss(system)
1852
A = G.A
1953
B = G.B
@@ -35,17 +69,22 @@ def norm(system, p=2, tol=1e-6):
3569
return np.sqrt(np.trace(C@P@C.T + D@D.T))
3670

3771
elif p == "inf": # L_infinity-norm
38-
def Hamilton_matrix(gamma):
72+
def _Hamilton_matrix(gamma):
3973
"""Constructs Hamiltonian matrix."""
4074
R = Ip*gamma**2 - D.T@D
4175
invR = la.inv(R)
4276
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]])
4377

44-
if G.isdtime(): # Bilinear transformation to s-plane
78+
if G.isdtime(): # Bilinear transformation of discrete time system to s-plane if no poles at |z|=1 or z=0
4579
Ad = A
4680
Bd = B
4781
Cd = C
4882
Dd = D
83+
if any(np.isclose(abs(la.eigvals(Ad)), 1.0)):
84+
return float('inf')
85+
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.")
87+
return None
4988
In = np.eye(len(Ad))
5089
Adinv = la.inv(Ad+In)
5190
A = 2*(Ad-In)@Adinv
@@ -60,16 +99,16 @@ def Hamilton_matrix(gamma):
6099
gamu = max(1.0, 2.0*gaml) # Candidate upper bound
61100
Ip = np.eye(len(D))
62101

63-
while any(np.isclose(la.eigvals(Hamilton_matrix(gamu)).real, 0.0)): # Find an upper bound
102+
while any(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)): # Find an upper bound
64103
gamu *= 2.0
65104

66105
while (gamu-gaml)/gamu > tol:
67106
gam = (gamu+gaml)/2.0
68-
if any(np.isclose(la.eigvals(Hamilton_matrix(gam)).real, 0.0)):
107+
if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):
69108
gaml = gam
70109
else:
71110
gamu = gam
72111
return gam
73112
else:
74-
# Norm computation only supported for p=2 and p='inf'
113+
print("Norm computation for p =", p, "currently not supported.")
75114
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, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB
17+
assert np.allclose(ct.norm(G1, p=2), 0.707106781186547, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB
21+
assert np.allclose(ct.norm(Gd1, p=2), 0.223513699524858, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB
59+
assert np.allclose(ct.norm(G4, p=2), 2.237461821810309, rtol=1e-09, atol=1e-08) # 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, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB
63+
assert np.allclose(ct.norm(Gd4, p=2), 0.707434962289554, rtol=1e-09, atol=1e-08) # Comparison to norm computed in MATLAB

0 commit comments

Comments
 (0)