Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commit457c623

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 parenta0fbc80 commit457c623

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
importcontrolasct
1818

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

21-
defnorm(system,p=2,tol=1e-6):
23+
defnorm(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+
Normvalueof 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-
ifp==2:# H_2-norm
61+
#
62+
# H_2-norm computation
63+
#
64+
ifp==2:
65+
# Continuous time case
5866
ifG.isctime():
59-
if (D!=0).any()orany(G.poles().real>=0):
67+
poles_real_part=G.poles().real
68+
ifany(np.isclose(poles_real_part,0.0)):
69+
ifprint_warning:
70+
print("Warning: Poles close to, or on, the imaginary axis. Norm value may be uncertain.")
71+
returnfloat('inf')
72+
elif (D!=0).any()orany(poles_real_part>0.0):# System unstable or has direct feedthrough?
6073
returnfloat('inf')
6174
else:
62-
P=ct.lyap(A,B@B.T)
63-
returnnp.sqrt(np.trace(C@P@C.T))
75+
try:
76+
P=ct.lyap(A,B@B.T)
77+
exceptExceptionase:
78+
print(f"An error occurred solving the continuous time Lyapunov equation:{e}")
79+
returnNone
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+
ifany(la.eigvals(P)<0.0):
84+
ifprint_warning:
85+
print("Warning: There appears to be poles close to the imaginary axis. Norm value may be uncertain.")
86+
returnfloat('inf')
87+
else:
88+
norm_value=np.sqrt(np.trace(C@P@C.T))# Argument in sqrt should be non-negative
89+
ifnp.isnan(norm_value):
90+
print("Unknown error. Norm computation resulted in NaN.")
91+
returnNone
92+
else:
93+
returnnorm_value
94+
95+
# Discrete time case
6496
elifG.isdtime():
65-
ifany(abs(G.poles())>=1):
97+
poles_abs=abs(G.poles())
98+
ifany(np.isclose(poles_abs,1.0)):
99+
ifprint_warning:
100+
print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
101+
returnfloat('inf')
102+
elifany(poles_abs>1.0):# System unstable?
66103
returnfloat('inf')
67104
else:
68-
P=ct.dlyap(A,B@B.T)
69-
returnnp.sqrt(np.trace(C@P@C.T+D@D.T))
70-
71-
elifp=="inf":# L_infinity-norm
105+
try:
106+
P=ct.dlyap(A,B@B.T)
107+
exceptExceptionase:
108+
print(f"An error occurred solving the discrete time Lyapunov equation:{e}")
109+
returnNone
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+
ifany(la.eigvals(P)<0.0):
114+
ifprint_warning:
115+
print("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.")
116+
returnfloat('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+
ifnp.isnan(norm_value):
120+
print("Unknown error. Norm computation resulted in NaN.")
121+
returnNone
122+
else:
123+
returnnorm_value
124+
#
125+
# L_infinity-norm computation
126+
#
127+
elifp=="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
returnnp.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-
ifG.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+
ifG.isdtime():
79138
Ad=A
80139
Bd=B
81140
Cd=C
82141
Dd=D
83142
ifany(np.isclose(abs(la.eigvals(Ad)),1.0)):
143+
ifprint_warning:
144+
print("Warning: Poles close to, or on, the complex unit circle. Norm value may be uncertain.")
84145
returnfloat('inf')
85146
elifany(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
returnNone
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
ifany(np.isclose(la.eigvals(A).real,0.0)):
160+
ifprint_warning:
161+
print("Warning: Poles close to, or on, imaginary axis. Norm value may be uncertain.")
96162
returnfloat('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-
whileany(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real,0.0)):# Findan upper bound
168+
whileany(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real,0.0)):# Findactual 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
returngam
112178
else:
113-
print("Norm computation for p =",p,"currently not supported.")
179+
print(f"Norm computation for p={p}currently not supported.")
114180
returnNone

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp