13107996232

Committed 03 Feb 2025 06:53AM UTC coverage: 94.731% (+0.02%) from 94.709%

Build # 13107996232

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1094 from murrayrm/userguide-22Dec2024

Updated user documentation (User Guide, Reference Manual)

Run Details

9673 of 10211 relevant lines covered (94.73%)

8.28 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

72.41

control/sysnorm.py

# sysnorm.py - functions for computing system norms
#
# Initial author: Henrik Sandberg
# Creation date: 21 Dec 2023

"""Functions for computing system norms."""

import warnings

import numpy as np
import numpy.linalg as la

import control as ct

__all__ = ['system_norm', 'norm']

#------------------------------------------------------------------------------

def _h2norm_slycot(sys, print_warning=True):
    """H2 norm of a linear system. For internal use. Requires Slycot.

    See Also
    --------
    slycot.ab13bd

    """
    # See: https://github.com/python-control/Slycot/issues/199
    try:
        from slycot import ab13bd
    except ImportError:
        ct.ControlSlycot("Can't find slycot module ab13bd")

    try:
        from slycot.exceptions import SlycotArithmeticError
    except ImportError:
        raise ct.ControlSlycot(
            "Can't find slycot class SlycotArithmeticError")

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

    n = A.shape[0]
    m = B.shape[1]
    p = C.shape[0]

    dico = 'C' if sys.isctime() else 'D'  # Continuous or discrete time
    jobn = 'H'  # H2 (and not L2 norm)

    if n == 0:
    # ab13bd does not accept empty A, B, C
        if dico == 'C':
            if any(D.flat != 0):
                if print_warning:
                    warnings.warn(
                        "System has a direct feedthrough term!", UserWarning)
                return float("inf")
            else:
                return 0.0
        elif dico == 'D':
            return np.sqrt(D@D.T)

    try:
        norm = ab13bd(dico, jobn, n, m, p, A, B, C, D)
    except SlycotArithmeticError as e:
        if e.info == 3:
            if print_warning:
                warnings.warn(
                    "System has pole(s) on the stability boundary!",
                    UserWarning)
            return float("inf")
        elif e.info == 5:
            if print_warning:
                warnings.warn(
                    "System has a direct feedthrough term!", UserWarning)
            return float("inf")
        elif e.info == 6:
            if print_warning:
                warnings.warn("System is unstable!", UserWarning)
            return float("inf")
        else:
            raise e
    return norm

#------------------------------------------------------------------------------

def system_norm(system, p=2, tol=1e-6, print_warning=True, method=None):
    """Computes the input/output norm of system.

    Parameters
    ----------
    system : LTI (`StateSpace` or `TransferFunction`)
        System in continuous or discrete time for which the norm should
        be computed.
    p : int or str
        Type of norm to be computed. `p` = 2 gives the H2 norm, and
        `p` = 'inf' gives the L-infinity norm.
    tol : float
        Relative tolerance for accuracy of L-infinity norm
        computation. Ignored unless `p` = 'inf'.
    print_warning : bool
        Print warning message in case norm value may be uncertain.
    method : str, optional
        Set the method used for computing the result.  Current methods are
        'slycot' and 'scipy'. If set to None (default), try 'slycot' first
        and then 'scipy'.

    Returns
    -------
    norm_value : float
        Norm value of system.

    Notes
    -----
    Does not yet compute the L-infinity norm for discrete-time systems
    with pole(s) at the origin unless Slycot is used.

    Examples
    --------
    >>> Gc = ct.tf([1], [1, 2, 1])
    >>> round(ct.norm(Gc, 2), 3)
    0.5
    >>> round(ct.norm(Gc, 'inf', tol=1e-5, method='scipy'), 3)
    np.float64(1.0)

    """
    if not isinstance(system, (ct.StateSpace, ct.TransferFunction)):
        raise TypeError(
            "Parameter `system`: must be a `StateSpace` or `TransferFunction`")

    G = ct.ss(system)
    A = G.A
    B = G.B
    C = G.C
    D = G.D

    # Decide what method to use
    method = ct.mateqn._slycot_or_scipy(method)

    # -------------------
    # H2 norm computation
    # -------------------
    if p == 2:
        # --------------------
        # Continuous time case
        # --------------------
        if G.isctime():

            # Check for cases with infinite norm
            poles_real_part = G.poles().real
            if any(np.isclose(poles_real_part, 0.0)):  # Poles on imaginary axis
                if print_warning:
                    warnings.warn(
                        "Poles close to, or on, the imaginary axis. "
                        "Norm value may be uncertain.", UserWarning)
                return float('inf')
            elif any(poles_real_part > 0.0):  # System unstable
                if print_warning:
                    warnings.warn("System is unstable!", UserWarning)
                return float('inf')
            elif any(D.flat != 0):  # System has direct feedthrough
                if print_warning:
                    warnings.warn(
                        "System has a direct feedthrough term!", UserWarning)
                return float('inf')

            else:
                # Use slycot, if available, to compute (finite) norm
                if method == 'slycot':
                    return _h2norm_slycot(G, print_warning)

                # Else use scipy
                else:
                    # Solve for controllability Gramian
                    P = ct.lyap(A, B@B.T, method=method)

                    # System is stable to reach this point, and P should be
                    # positive semi-definite.  Test next is a precaution in
                    # case the Lyapunov equation is ill conditioned.
                    if any(la.eigvals(P).real < 0.0):
                        if print_warning:
                            warnings.warn(
                                "There appears to be poles close to the "
                                "imaginary axis. Norm value may be uncertain.",
                                UserWarning)
                        return float('inf')
                    else:
                        # Argument in sqrt should be non-negative
                        norm_value = np.sqrt(np.trace(C@P@C.T))
                        if np.isnan(norm_value):
                            raise ct.ControlArgument(
                                "Norm computation resulted in NaN.")
                        else:
                            return norm_value

        # ------------------
        # Discrete time case
        # ------------------
        elif G.isdtime():

            # Check for cases with infinite norm
            poles_abs = abs(G.poles())
            if any(np.isclose(poles_abs, 1.0)):  # Poles on imaginary axis
                if print_warning:
                    warnings.warn(
                        "Poles close to, or on, the complex unit circle. "
                        "Norm value may be uncertain.", UserWarning)
                return float('inf')
            elif any(poles_abs > 1.0):  # System unstable
                if print_warning:
                    warnings.warn("System is unstable!", UserWarning)
                return float('inf')
            else:
                # Use slycot, if available, to compute (finite) norm
                if method == 'slycot':
                    return _h2norm_slycot(G, print_warning)

                # Else use scipy
                else:
                    P = ct.dlyap(A, B@B.T, method=method)

                # System is stable to reach this point, and P should be
                # positive semi-definite.  Test next is a precaution in
                # case the Lyapunov equation is ill conditioned.
                if any(la.eigvals(P).real < 0.0):
                    if print_warning:
                        warnings.warn(
                            "There appears to be poles close to the complex "
                            "unit circle. Norm value may be uncertain.",
                            UserWarning)
                    return float('inf')
                else:
                    # Argument in sqrt should be non-negative
                    norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T))
                    if np.isnan(norm_value):
                        raise ct.ControlArgument(
                            "Norm computation resulted in NaN.")
                    else:
                        return norm_value

    # ---------------------------
    # L-infinity norm computation
    # ---------------------------
    elif p == "inf":

        # Check for cases with infinite norm
        poles = G.poles()
        if G.isdtime():  # Discrete time
            if any(np.isclose(abs(poles), 1.0)):  # Poles on unit circle
                if print_warning:
                    warnings.warn(
                        "Poles close to, or on, the complex unit circle. "
                        "Norm value may be uncertain.", UserWarning)
                return float('inf')
        else:  # Continuous time
            if any(np.isclose(poles.real, 0.0)):  # Poles on imaginary axis
                if print_warning:
                    warnings.warn(
                        "Poles close to, or on, the imaginary axis. "
                        "Norm value may be uncertain.", UserWarning)
                return float('inf')

        # Use slycot, if available, to compute (finite) norm
        if method == 'slycot':
            return ct.linfnorm(G, tol)[0]

        # Else use scipy
        else:

            # ------------------
            # Discrete time case
            # ------------------
            # Use inverse bilinear transformation of discrete-time system
            # to s-plane if no poles on |z|=1 or z=0.  Allows us to use
            # test for continuous-time systems next.
            if G.isdtime():
                Ad = A
                Bd = B
                Cd = C
                Dd = D
                if any(np.isclose(la.eigvals(Ad), 0.0)):
                    raise ct.ControlArgument(
                        "L-infinity norm computation for discrete-time "
                        "system with pole(s) in z=0 currently not supported "
                        "unless Slycot installed.")

                # Inverse bilinear transformation
                In = np.eye(len(Ad))
                Adinv = la.inv(Ad+In)
                A = 2*(Ad-In)@Adinv
                B = 2*Adinv@Bd
                C = 2*Cd@Adinv
                D = Dd - Cd@Adinv@Bd

            # --------------------
            # Continuous time case
            # --------------------
            def _Hamilton_matrix(gamma):
                """Constructs Hamiltonian matrix. For internal use."""
                R = Ip*gamma**2 - D.T@D
                invR = la.inv(R)
                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]])

            gaml = la.norm(D,ord=2)    # Lower bound
            gamu = max(1.0, 2.0*gaml)  # Candidate upper bound
            Ip = np.eye(len(D))

            while any(np.isclose(
                    la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)):
                # Find actual upper bound
                gamu *= 2.0

            while (gamu-gaml)/gamu > tol:
                gam = (gamu+gaml)/2.0
                if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):
                    gaml = gam
                else:
                    gamu = gam
            return gam

    # ----------------------
    # Other norm computation
    # ----------------------
    else:
        raise ct.ControlArgument(
            f"Norm computation for p={p} currently not supported.")


norm = system_norm

1	# sysnorm.py - functions for computing system norms
2	#
3	# Initial author: Henrik Sandberg
4	# Creation date: 21 Dec 2023
5
6	"""Functions for computing system norms."""
7
8	import warnings	9✔
9
10	import numpy as np	9✔
11	import numpy.linalg as la	9✔
12
13	import control as ct	9✔
14
15	__all__ = ['system_norm', 'norm']	9✔
16
17	#------------------------------------------------------------------------------
18
19	def _h2norm_slycot(sys, print_warning=True):	9✔
20	"""H2 norm of a linear system. For internal use. Requires Slycot.
21
22	See Also
23	--------
24	slycot.ab13bd
25
26	"""
27	# See: https://github.com/python-control/Slycot/issues/199
28	try:	5✔
29	from slycot import ab13bd	5✔
30	except ImportError:	×
31	ct.ControlSlycot("Can't find slycot module ab13bd")	×
32
33	try:	5✔
34	from slycot.exceptions import SlycotArithmeticError	5✔
35	except ImportError:	×
36	raise ct.ControlSlycot(	×
37	"Can't find slycot class SlycotArithmeticError")
38
39	A, B, C, D = ct.ssdata(ct.ss(sys))	5✔
40
41	n = A.shape[0]	5✔
42	m = B.shape[1]	5✔
43	p = C.shape[0]	5✔
44
45	dico = 'C' if sys.isctime() else 'D' # Continuous or discrete time	5✔
46	jobn = 'H' # H2 (and not L2 norm)	5✔
47
48	if n == 0:	5✔
49	# ab13bd does not accept empty A, B, C
50	if dico == 'C':	×
51	if any(D.flat != 0):	×
52	if print_warning:	×
53	warnings.warn(	×
54	"System has a direct feedthrough term!", UserWarning)
55	return float("inf")	×
56	else:
57	return 0.0	×
58	elif dico == 'D':	×
59	return np.sqrt(D@D.T)	×
60
61	try:	5✔
62	norm = ab13bd(dico, jobn, n, m, p, A, B, C, D)	5✔
63	except SlycotArithmeticError as e:	×
64	if e.info == 3:	×
65	if print_warning:	×
66	warnings.warn(	×
67	"System has pole(s) on the stability boundary!",
68	UserWarning)
69	return float("inf")	×
70	elif e.info == 5:	×
71	if print_warning:	×
72	warnings.warn(	×
73	"System has a direct feedthrough term!", UserWarning)
74	return float("inf")	×
75	elif e.info == 6:	×
76	if print_warning:	×
77	warnings.warn("System is unstable!", UserWarning)	×
78	return float("inf")	×
79	else:
80	raise e	×
81	return norm	5✔
82
83	#------------------------------------------------------------------------------
84
85	def system_norm(system, p=2, tol=1e-6, print_warning=True, method=None):	9✔
86	"""Computes the input/output norm of system.
87
88	Parameters
89	----------
90	system : LTI (`StateSpace` or `TransferFunction`)
91	System in continuous or discrete time for which the norm should
92	be computed.
93	p : int or str
94	Type of norm to be computed. `p` = 2 gives the H2 norm, and
95	`p` = 'inf' gives the L-infinity norm.
96	tol : float
97	Relative tolerance for accuracy of L-infinity norm
98	computation. Ignored unless `p` = 'inf'.
99	print_warning : bool
100	Print warning message in case norm value may be uncertain.
101	method : str, optional
102	Set the method used for computing the result. Current methods are
103	'slycot' and 'scipy'. If set to None (default), try 'slycot' first
104	and then 'scipy'.
105
106	Returns
107	-------
108	norm_value : float
109	Norm value of system.
110
111	Notes
112	-----
113	Does not yet compute the L-infinity norm for discrete-time systems
114	with pole(s) at the origin unless Slycot is used.
115
116	Examples
117	--------
118	>>> Gc = ct.tf([1], [1, 2, 1])
119	>>> round(ct.norm(Gc, 2), 3)
120	0.5
121	>>> round(ct.norm(Gc, 'inf', tol=1e-5, method='scipy'), 3)
122	np.float64(1.0)
123
124	"""
125	if not isinstance(system, (ct.StateSpace, ct.TransferFunction)):	9✔
126	raise TypeError(	×
127	"Parameter `system`: must be a `StateSpace` or `TransferFunction`")
128
129	G = ct.ss(system)	9✔
130	A = G.A	9✔
131	B = G.B	9✔
132	C = G.C	9✔
133	D = G.D	9✔
134
135	# Decide what method to use
136	method = ct.mateqn._slycot_or_scipy(method)	9✔
137
138	# -------------------
139	# H2 norm computation
140	# -------------------
141	if p == 2:	9✔
142	# --------------------
143	# Continuous time case
144	# --------------------
145	if G.isctime():	9✔
146
147	# Check for cases with infinite norm
148	poles_real_part = G.poles().real	9✔
149	if any(np.isclose(poles_real_part, 0.0)): # Poles on imaginary axis	9✔
150	if print_warning:	9✔
151	warnings.warn(	9✔
152	"Poles close to, or on, the imaginary axis. "
153	"Norm value may be uncertain.", UserWarning)
154	return float('inf')	9✔
155	elif any(poles_real_part > 0.0): # System unstable	9✔
156	if print_warning:	9✔
157	warnings.warn("System is unstable!", UserWarning)	9✔
158	return float('inf')	9✔
159	elif any(D.flat != 0): # System has direct feedthrough	9✔
160	if print_warning:	×
161	warnings.warn(	×
162	"System has a direct feedthrough term!", UserWarning)
163	return float('inf')	×
164
165	else:
166	# Use slycot, if available, to compute (finite) norm
167	if method == 'slycot':	9✔
168	return _h2norm_slycot(G, print_warning)	5✔
169
170	# Else use scipy
171	else:
172	# Solve for controllability Gramian
173	P = ct.lyap(A, B@B.T, method=method)	4✔
174
175	# System is stable to reach this point, and P should be
176	# positive semi-definite. Test next is a precaution in
177	# case the Lyapunov equation is ill conditioned.
178	if any(la.eigvals(P).real < 0.0):	4✔
179	if print_warning:	×
180	warnings.warn(	×
181	"There appears to be poles close to the "
182	"imaginary axis. Norm value may be uncertain.",
183	UserWarning)
184	return float('inf')	×
185	else:
186	# Argument in sqrt should be non-negative
187	norm_value = np.sqrt(np.trace(C@P@C.T))	4✔
188	if np.isnan(norm_value):	4✔
189	raise ct.ControlArgument(	×
190	"Norm computation resulted in NaN.")
191	else:
192	return norm_value	4✔
193
194	# ------------------
195	# Discrete time case
196	# ------------------
197	elif G.isdtime():	9✔
198
199	# Check for cases with infinite norm
200	poles_abs = abs(G.poles())	9✔
201	if any(np.isclose(poles_abs, 1.0)): # Poles on imaginary axis	9✔
202	if print_warning:	9✔
203	warnings.warn(	9✔
204	"Poles close to, or on, the complex unit circle. "
205	"Norm value may be uncertain.", UserWarning)
206	return float('inf')	9✔
207	elif any(poles_abs > 1.0): # System unstable	9✔
208	if print_warning:	9✔
209	warnings.warn("System is unstable!", UserWarning)	9✔
210	return float('inf')	9✔
211	else:
212	# Use slycot, if available, to compute (finite) norm
213	if method == 'slycot':	9✔
214	return _h2norm_slycot(G, print_warning)	5✔
215
216	# Else use scipy
217	else:
218	P = ct.dlyap(A, B@B.T, method=method)	4✔
219
220	# System is stable to reach this point, and P should be
221	# positive semi-definite. Test next is a precaution in
222	# case the Lyapunov equation is ill conditioned.
223	if any(la.eigvals(P).real < 0.0):	4✔
224	if print_warning:	×
225	warnings.warn(	×
226	"There appears to be poles close to the complex "
227	"unit circle. Norm value may be uncertain.",
228	UserWarning)
229	return float('inf')	×
230	else:
231	# Argument in sqrt should be non-negative
232	norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T))	4✔
233	if np.isnan(norm_value):	4✔
234	raise ct.ControlArgument(	×
235	"Norm computation resulted in NaN.")
236	else:
237	return norm_value	4✔
238
239	# ---------------------------
240	# L-infinity norm computation
241	# ---------------------------
242	elif p == "inf":	9✔
243
244	# Check for cases with infinite norm
245	poles = G.poles()	9✔
246	if G.isdtime(): # Discrete time	9✔
247	if any(np.isclose(abs(poles), 1.0)): # Poles on unit circle	9✔
248	if print_warning:	9✔
249	warnings.warn(	9✔
250	"Poles close to, or on, the complex unit circle. "
251	"Norm value may be uncertain.", UserWarning)
252	return float('inf')	9✔
253	else: # Continuous time
254	if any(np.isclose(poles.real, 0.0)): # Poles on imaginary axis	9✔
255	if print_warning:	9✔
256	warnings.warn(	9✔
257	"Poles close to, or on, the imaginary axis. "
258	"Norm value may be uncertain.", UserWarning)
259	return float('inf')	9✔
260
261	# Use slycot, if available, to compute (finite) norm
262	if method == 'slycot':	9✔
263	return ct.linfnorm(G, tol)[0]	5✔
264
265	# Else use scipy
266	else:
267
268	# ------------------
269	# Discrete time case
270	# ------------------
271	# Use inverse bilinear transformation of discrete-time system
272	# to s-plane if no poles on \|z\|=1 or z=0. Allows us to use
273	# test for continuous-time systems next.
274	if G.isdtime():	4✔
275	Ad = A	4✔
276	Bd = B	4✔
277	Cd = C	4✔
278	Dd = D	4✔
279	if any(np.isclose(la.eigvals(Ad), 0.0)):	4✔
280	raise ct.ControlArgument(	×
281	"L-infinity norm computation for discrete-time "
282	"system with pole(s) in z=0 currently not supported "
283	"unless Slycot installed.")
284
285	# Inverse bilinear transformation
286	In = np.eye(len(Ad))	4✔
287	Adinv = la.inv(Ad+In)	4✔
288	A = 2*(Ad-In)@Adinv	4✔
289	B = 2*Adinv@Bd	4✔
290	C = 2*Cd@Adinv	4✔
291	D = Dd - Cd@Adinv@Bd	4✔
292
293	# --------------------
294	# Continuous time case
295	# --------------------
296	def _Hamilton_matrix(gamma):	4✔
297	"""Constructs Hamiltonian matrix. For internal use."""
298	R = Ipgamma*2 - D.T@D	4✔
299	invR = la.inv(R)	4✔
300	return np.block([	4✔
301	[A+B@invR@D.T@C, B@invR@B.T],
302	[-C.T@(Ip+D@invR@D.T)@C, -(A+B@invR@D.T@C).T]])
303
304	gaml = la.norm(D,ord=2) # Lower bound	4✔
305	gamu = max(1.0, 2.0*gaml) # Candidate upper bound	4✔
306	Ip = np.eye(len(D))	4✔
307
308	while any(np.isclose(	4✔
309	la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)):
310	# Find actual upper bound
311	gamu *= 2.0	4✔
312
313	while (gamu-gaml)/gamu > tol:	4✔
314	gam = (gamu+gaml)/2.0	4✔
315	if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):	4✔
316	gaml = gam	4✔
317	else:
318	gamu = gam	4✔
319	return gam	4✔
320
321	# ----------------------
322	# Other norm computation
323	# ----------------------
324	else:
325	raise ct.ControlArgument(	×
326	f"Norm computation for p={p} currently not supported.")
327
328
329	norm = system_norm	9✔

python-control / python-control / 13107996232

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous