13096124094

Committed 02 Feb 2025 05:58AM UTC coverage: 94.661% (+0.002%) from 94.659%

Build # 13096124094

Build Type

Pull #1116

github

Committed by

web-flow

Commit Message

Merge 4fc70f7b0 into 2cb0520b6

Pull Request Pull Request #1116: Update _ssmatrix and _check_shape for consistent usage

Run Details

9663 of 10208 relevant lines covered (94.66%)

8.27 hits per line

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

72.41

control/sysnorm.py

# -*- coding: utf-8 -*-
"""sysnorm.py

Functions for computing system norms.

Routine in this module:

norm

Created on Thu Dec 21 08:06:12 2023
Author: Henrik Sandberg
"""

import numpy as np
import scipy as sp
import numpy.linalg as la
import warnings

import control as ct

__all__ = ['norm']

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

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

    See also
    --------
    ``slycot.ab13bd`` : the Slycot routine that does the calculation
    https://github.com/python-control/Slycot/issues/199 : Post on issue with ``ab13bf``
    """
    
    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 norm(system, p=2, tol=1e-6, print_warning=True, method=None):
    """Computes norm of system.
    
    Parameters
    ----------
    system : LTI (:class:`StateSpace` or :class:`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) in z=0 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:
                    P = ct.lyap(A, B@B.T, method=method)  # Solve for controllability Gramian
                                        
                    # 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:
                        norm_value = np.sqrt(np.trace(C@P@C.T))  # Argument in sqrt should be non-negative
                        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("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.", UserWarning)
                    return float('inf')
                else:
                    norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T))  # Argument in sqrt should be non-negative               
                    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.")           


1	# -- coding: utf-8 --
2	"""sysnorm.py
3
4	Functions for computing system norms.
5
6	Routine in this module:
7
8	norm
9
10	Created on Thu Dec 21 08:06:12 2023
11	Author: Henrik Sandberg
12	"""
13
14	import numpy as np	9✔
15	import scipy as sp	9✔
16	import numpy.linalg as la	9✔
17	import warnings	9✔
18
19	import control as ct	9✔
20
21	__all__ = ['norm']	9✔
22
23	#------------------------------------------------------------------------------
24
25	def _h2norm_slycot(sys, print_warning=True):	9✔
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:	5✔
35	from slycot import ab13bd	5✔
36	except ImportError:	×
37	ct.ControlSlycot("Can't find slycot module ``ab13bd``!")	×
38
39	try:	5✔
40	from slycot.exceptions import SlycotArithmeticError	5✔
41	except ImportError:	×
42	raise ct.ControlSlycot("Can't find slycot class ``SlycotArithmeticError``!")	×
43
44	A, B, C, D = ct.ssdata(ct.ss(sys))	5✔
45
46	n = A.shape[0]	5✔
47	m = B.shape[1]	5✔
48	p = C.shape[0]	5✔
49
50	dico = 'C' if sys.isctime() else 'D' # Continuous or discrete time	5✔
51	jobn = 'H' # H2 (and not L2 norm)	5✔
52
53	if n == 0:	5✔
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!", UserWarning)	×
59	return float("inf")	×
60	else:
61	return 0.0	×
62	elif dico == 'D':	×
63	return np.sqrt(D@D.T)	×
64
65	try:	5✔
66	norm = ab13bd(dico, jobn, n, m, p, A, B, C, D)	5✔
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!", UserWarning)	×
71	return float("inf")	×
72	elif e.info == 5:	×
73	if print_warning:	×
74	warnings.warn("System has a direct feedthrough term!", UserWarning)	×
75	return float("inf")	×
76	elif e.info == 6:	×
77	if print_warning:	×
78	warnings.warn("System is unstable!", UserWarning)	×
79	return float("inf")	×
80	else:
81	raise e	×
82	return norm	5✔
83
84	#------------------------------------------------------------------------------
85
86	def norm(system, p=2, tol=1e-6, print_warning=True, method=None):	9✔
87	"""Computes norm of system.
88
89	Parameters
90	----------
91	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
99	Print warning message in case norm value may be uncertain.
100	method : str, optional
101	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
105	Returns
106	-------
107	norm_value : float
108	Norm value of system.
109
110	Notes
111	-----
112	Does not yet compute the L-infinity norm for discrete time systems with pole(s) in z=0 unless Slycot is used.
113
114	Examples
115	--------
116	>>> Gc = ct.tf([1], [1, 2, 1])
117	>>> round(ct.norm(Gc, 2), 3)
118	0.5
119	>>> round(ct.norm(Gc, 'inf', tol=1e-5, method='scipy'), 3)
120	np.float64(1.0)
121	"""
122
123	if not isinstance(system, (ct.StateSpace, ct.TransferFunction)):	9✔
124	raise TypeError('Parameter ``system``: must be a ``StateSpace`` or ``TransferFunction``')	×
125
126	G = ct.ss(system)	9✔
127	A = G.A	9✔
128	B = G.B	9✔
129	C = G.C	9✔
130	D = G.D	9✔
131
132	# Decide what method to use
133	method = ct.mateqn._slycot_or_scipy(method)	9✔
134
135	# -------------------
136	# H2 norm computation
137	# -------------------
138	if p == 2:	9✔
139	# --------------------
140	# Continuous time case
141	# --------------------
142	if G.isctime():	9✔
143
144	# Check for cases with infinite norm
145	poles_real_part = G.poles().real	9✔
146	if any(np.isclose(poles_real_part, 0.0)): # Poles on imaginary axis	9✔
147	if print_warning:	9✔
148	warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.", UserWarning)	9✔
149	return float('inf')	9✔
150	elif any(poles_real_part > 0.0): # System unstable	9✔
151	if print_warning:	9✔
152	warnings.warn("System is unstable!", UserWarning)	9✔
153	return float('inf')	9✔
154	elif any(D.flat != 0): # System has direct feedthrough	9✔
155	if print_warning:	×
156	warnings.warn("System has a direct feedthrough term!", UserWarning)	×
157	return float('inf')	×
158
159	else:
160	# Use slycot, if available, to compute (finite) norm
161	if method == 'slycot':	9✔
162	return _h2norm_slycot(G, print_warning)	5✔
163
164	# Else use scipy
165	else:
166	P = ct.lyap(A, B@B.T, method=method) # Solve for controllability Gramian	4✔
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):	4✔
171	if print_warning:	×
172	warnings.warn("There appears to be poles close to the imaginary axis. Norm value may be uncertain.", UserWarning)	×
173	return float('inf')	×
174	else:
175	norm_value = np.sqrt(np.trace(C@P@C.T)) # Argument in sqrt should be non-negative	4✔
176	if np.isnan(norm_value):	4✔
177	raise ct.ControlArgument("Norm computation resulted in NaN.")	×
178	else:
179	return norm_value	4✔
180
181	# ------------------
182	# Discrete time case
183	# ------------------
184	elif G.isdtime():	9✔
185
186	# Check for cases with infinite norm
187	poles_abs = abs(G.poles())	9✔
188	if any(np.isclose(poles_abs, 1.0)): # Poles on imaginary axis	9✔
189	if print_warning:	9✔
190	warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.", UserWarning)	9✔
191	return float('inf')	9✔
192	elif any(poles_abs > 1.0): # System unstable	9✔
193	if print_warning:	9✔
194	warnings.warn("System is unstable!", UserWarning)	9✔
195	return float('inf')	9✔
196
197	else:
198	# Use slycot, if available, to compute (finite) norm
199	if method == 'slycot':	9✔
200	return _h2norm_slycot(G, print_warning)	5✔
201
202	# Else use scipy
203	else:
204	P = ct.dlyap(A, B@B.T, method=method)	4✔
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.
208	if any(la.eigvals(P).real < 0.0):	4✔
209	if print_warning:	×
210	warnings.warn("Warning: There appears to be poles close to the complex unit circle. Norm value may be uncertain.", UserWarning)	×
211	return float('inf')	×
212	else:
213	norm_value = np.sqrt(np.trace(C@P@C.T + D@D.T)) # Argument in sqrt should be non-negative	4✔
214	if np.isnan(norm_value):	4✔
215	raise ct.ControlArgument("Norm computation resulted in NaN.")	×
216	else:
217	return norm_value	4✔
218
219	# ---------------------------
220	# L-infinity norm computation
221	# ---------------------------
222	elif p == "inf":	9✔
223
224	# Check for cases with infinite norm
225	poles = G.poles()	9✔
226	if G.isdtime(): # Discrete time	9✔
227	if any(np.isclose(abs(poles), 1.0)): # Poles on unit circle	9✔
228	if print_warning:	9✔
229	warnings.warn("Poles close to, or on, the complex unit circle. Norm value may be uncertain.", UserWarning)	9✔
230	return float('inf')	9✔
231	else: # Continuous time
232	if any(np.isclose(poles.real, 0.0)): # Poles on imaginary axis	9✔
233	if print_warning:	9✔
234	warnings.warn("Poles close to, or on, the imaginary axis. Norm value may be uncertain.", UserWarning)	9✔
235	return float('inf')	9✔
236
237	# Use slycot, if available, to compute (finite) norm
238	if method == 'slycot':	9✔
239	return ct.linfnorm(G, tol)[0]	5✔
240
241	# Else use scipy
242	else:
243
244	# ------------------
245	# Discrete time case
246	# ------------------
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.
249	if G.isdtime():	4✔
250	Ad = A	4✔
251	Bd = B	4✔
252	Cd = C	4✔
253	Dd = D	4✔
254	if any(np.isclose(la.eigvals(Ad), 0.0)):	4✔
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
257	# Inverse bilinear transformation
258	In = np.eye(len(Ad))	4✔
259	Adinv = la.inv(Ad+In)	4✔
260	A = 2*(Ad-In)@Adinv	4✔
261	B = 2*Adinv@Bd	4✔
262	C = 2*Cd@Adinv	4✔
263	D = Dd - Cd@Adinv@Bd	4✔
264
265	# --------------------
266	# Continuous time case
267	# --------------------
268	def _Hamilton_matrix(gamma):	4✔
269	"""Constructs Hamiltonian matrix. For internal use."""
270	R = Ipgamma*2 - D.T@D	4✔
271	invR = la.inv(R)	4✔
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]])	4✔
273
274	gaml = la.norm(D,ord=2) # Lower bound	4✔
275	gamu = max(1.0, 2.0*gaml) # Candidate upper bound	4✔
276	Ip = np.eye(len(D))	4✔
277
278	while any(np.isclose(la.eigvals(_Hamilton_matrix(gamu)).real, 0.0)): # Find actual upper bound	4✔
279	gamu *= 2.0	4✔
280
281	while (gamu-gaml)/gamu > tol:	4✔
282	gam = (gamu+gaml)/2.0	4✔
283	if any(np.isclose(la.eigvals(_Hamilton_matrix(gam)).real, 0.0)):	4✔
284	gaml = gam	4✔
285	else:
286	gamu = gam	4✔
287	return gam	4✔
288
289	# ----------------------
290	# Other norm computation
291	# ----------------------
292	else:
293	raise ct.ControlArgument(f"Norm computation for p={p} currently not supported.")	×
294

python-control / python-control / 13096124094

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