Skip to content

Commit a19c501

Browse files
Mark-YeatmanMark-Yeatman
authored
Passivity indices and support for discrete time systems. (#750)
This adds support for input and output passivity indices. This implementation ignores frequency bands, relative passivity index, combined I/O passivity, and directional passivity index. This also adds support for computing indices and "is passive or not" for discrete time systems. Authored-by: Mark <myeatman@oxefit.com>
1 parent 6f124ec commit a19c501

File tree

4 files changed

+339
-56
lines changed

4 files changed

+339
-56
lines changed

‎control/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,7 @@
7272
from .config import *
7373
from .sisotool import *
7474
from .iosys import *
75+
from .passivity import *
7576

7677
# Exceptions
7778
from .exception import *

‎control/passivity.py‎

Lines changed: 256 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -1,42 +1,72 @@
1-
'''
2-
Author: Mark Yeatman
3-
Date: May 15, 2022
4-
'''
1+
"""
2+
Functions for passive control.
3+
4+
Author: Mark Yeatman
5+
Date: July 17, 2022
6+
"""
57

68
import numpy as np
7-
from control import statesp as ss
9+
from control import statesp
10+
from control.exception import ControlArgument, ControlDimension
811

912
try:
1013
import cvxopt as cvx
11-
except ImportError as e:
14+
except ImportError:
1215
cvx = None
1316

17+
eps = np.nextafter(0, 1)
18+
19+
__all__ = ["get_output_fb_index", "get_input_ff_index", "ispassive",
20+
"solve_passivity_LMI"]
21+
22+
23+
def solve_passivity_LMI(sys, rho=None, nu=None):
24+
"""Compute passivity indices and/or solves feasiblity via a LMI.
1425
15-
def ispassive(sys):
16-
'''
17-
Indicates if a linear time invariant (LTI) system is passive
26+
Constructs a linear matrix inequality (LMI) such that if a solution exists
27+
and the last element of the solution is positive, the system `sys` is
28+
passive. Inputs of None for `rho` or `nu` indicate that the function should
29+
solve for that index (they are mutually exclusive, they can't both be
30+
None, otherwise you're trying to solve a nonconvex bilinear matrix
31+
inequality.) The last element of `solution` is either the output or input
32+
passivity index, for `rho` = None and `nu` = None respectively.
1833
19-
Constructs a linear matrix inequality and a feasibility optimization
20-
such that if a solution exists, the system is passive.
34+
The sources for the algorithm are:
2135
22-
The source for the algorithm is:
23-
McCourt, Michael J., and Panos J. Antsaklis.
24-
"Demonstrating passivity and dissipativity using computational methods." ISIS 8 (2013).
36+
McCourt, Michael J., and Panos J. Antsaklis
37+
"Demonstrating passivity and dissipativity using computational
38+
methods."
39+
40+
Nicholas Kottenstette and Panos J. Antsaklis
41+
"Relationships Between Positive Real, Passive Dissipative, & Positive
42+
Systems", equation 36.
2543
2644
Parameters
2745
----------
28-
sys: A continuous LTI system
29-
System to be checked.
46+
sys : LTI
47+
System to be checked
48+
rho : float or None
49+
Output feedback passivity index
50+
nu : float or None
51+
Input feedforward passivity index
3052
3153
Returns
3254
-------
33-
bool:
34-
The input system passive.
35-
'''
55+
solution : ndarray
56+
The LMI solution
57+
"""
3658
if cvx is None:
3759
raise ModuleNotFoundError("cvxopt required for passivity module")
3860

39-
sys = ss._convert_to_statespace(sys)
61+
if sys.ninputs != sys.noutputs:
62+
raise ControlDimension(
63+
"The number of system inputs must be the same as the number of "
64+
"system outputs.")
65+
66+
if rho is None and nu is None:
67+
raise ControlArgument("rho or nu must be given a numerical value.")
68+
69+
sys = statesp._convert_to_statespace(sys)
4070

4171
A = sys.A
4272
B = sys.B
@@ -45,43 +75,218 @@ def ispassive(sys):
4575

4676
# account for strictly proper systems
4777
[n, m] = D.shape
48-
D = D + np.nextafter(0, 1)*np.eye(n, m)
78+
D = D + eps * np.eye(n, m)
4979

5080
[n, _] = A.shape
51-
A = A - np.nextafter(0, 1)*np.eye(n)
52-
53-
def make_LMI_matrix(P):
54-
V = np.vstack((
55-
np.hstack((A.T @ P + P@A, P@B)),
56-
np.hstack((B.T@P, np.zeros_like(D))))
57-
)
58-
return V
59-
60-
matrix_list = []
61-
state_space_size = sys.nstates
62-
for i in range(0, state_space_size):
63-
for j in range(0, state_space_size):
64-
if j <= i:
65-
P = np.zeros_like(A)
66-
P[i, j] = 1.0
67-
P[j, i] = 1.0
68-
matrix_list.append(make_LMI_matrix(P).flatten())
69-
70-
coefficents = np.vstack(matrix_list).T
71-
72-
constants = -np.vstack((
73-
np.hstack((np.zeros_like(A), - C.T)),
74-
np.hstack((- C, -D - D.T)))
75-
)
81+
A = A - eps*np.eye(n)
82+
83+
def make_LMI_matrix(P, rho, nu, one):
84+
q = sys.noutputs
85+
Q = -rho*np.eye(q, q)
86+
S = 1.0/2.0*(one+rho*nu)*np.eye(q)
87+
R = -nu*np.eye(m)
88+
if sys.isctime():
89+
off_diag = P@B - (C.T@S + C.T@Q@D)
90+
return np.vstack((
91+
np.hstack((A.T @ P + P@A - C.T@Q@C, off_diag)),
92+
np.hstack((off_diag.T, -(D.T@Q@D + D.T@S + S.T@D + R)))
93+
))
94+
else:
95+
off_diag = A.T@P@B - (C.T@S + C.T@Q@D)
96+
return np.vstack((
97+
np.hstack((A.T @ P @ A - P - C.T@Q@C, off_diag)),
98+
np.hstack((off_diag.T, B.T@P@B-(D.T@Q@D + D.T@S + S.T@D + R)))
99+
))
100+
101+
def make_P_basis_matrices(n, rho, nu):
102+
"""Make list of matrix constraints for passivity LMI.
103+
104+
Utility function to make basis matrices for a LMI from a
105+
symmetric matrix P of size n by n representing a parametrized symbolic
106+
matrix
107+
"""
108+
matrix_list = []
109+
for i in range(0, n):
110+
for j in range(0, n):
111+
if j <= i:
112+
P = np.zeros((n, n))
113+
P[i, j] = 1
114+
P[j, i] = 1
115+
matrix_list.append(make_LMI_matrix(P, 0, 0, 0).flatten())
116+
zeros = eps*np.eye(n)
117+
if rho is None:
118+
matrix_list.append(make_LMI_matrix(zeros, 1, 0, 0).flatten())
119+
elif nu is None:
120+
matrix_list.append(make_LMI_matrix(zeros, 0, 1, 0).flatten())
121+
return matrix_list
122+
123+
124+
def P_pos_def_constraint(n):
125+
"""Make a list of matrix constraints for P >= 0.
126+
127+
Utility function to make basis matrices for a LMI that ensures
128+
parametrized symbolic matrix of size n by n is positive definite
129+
"""
130+
matrix_list = []
131+
for i in range(0, n):
132+
for j in range(0, n):
133+
if j <= i:
134+
P = np.zeros((n, n))
135+
P[i, j] = -1
136+
P[j, i] = -1
137+
matrix_list.append(P.flatten())
138+
if rho is None or nu is None:
139+
matrix_list.append(np.zeros((n, n)).flatten())
140+
return matrix_list
141+
142+
n = sys.nstates
143+
144+
# coefficents for passivity indices and feasibility matrix
145+
sys_matrix_list = make_P_basis_matrices(n, rho, nu)
76146

147+
# get constants for numerical values of rho and nu
148+
sys_constants = list()
149+
if rho is not None and nu is not None:
150+
sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, nu, 1.0)
151+
elif rho is not None:
152+
sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, eps, 1.0)
153+
elif nu is not None:
154+
sys_constants = -make_LMI_matrix(np.zeros_like(A), eps, nu, 1.0)
155+
156+
sys_coefficents = np.vstack(sys_matrix_list).T
157+
158+
# LMI to ensure P is positive definite
159+
P_matrix_list = P_pos_def_constraint(n)
160+
P_coefficents = np.vstack(P_matrix_list).T
161+
P_constants = np.zeros((n, n))
162+
163+
# cost function
77164
number_of_opt_vars = int(
78-
(state_space_size**2-state_space_size)/2 + state_space_size)
165+
(n**2-n)/2 + n)
79166
c = cvx.matrix(0.0, (number_of_opt_vars, 1))
80167

168+
#we're maximizing a passivity index, include it in the cost function
169+
if rho is None or nu is None:
170+
c = cvx.matrix(np.append(np.array(c), -1.0))
171+
172+
Gs = [cvx.matrix(sys_coefficents)] + [cvx.matrix(P_coefficents)]
173+
hs = [cvx.matrix(sys_constants)] + [cvx.matrix(P_constants)]
174+
81175
# crunch feasibility solution
82176
cvx.solvers.options['show_progress'] = False
83-
sol = cvx.solvers.sdp(c,
84-
Gs=[cvx.matrix(coefficents)],
85-
hs=[cvx.matrix(constants)])
177+
sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
178+
return sol["x"]
179+
180+
181+
def get_output_fb_index(sys):
182+
"""Return the output feedback passivity (OFP) index for the system.
183+
184+
The OFP is the largest gain that can be placed in positive feedback
185+
with a system such that the new interconnected system is passive.
186+
187+
Parameters
188+
----------
189+
sys : LTI
190+
System to be checked
191+
192+
Returns
193+
-------
194+
float
195+
The OFP index
196+
"""
197+
sol = solve_passivity_LMI(sys, nu=eps)
198+
if sol is None:
199+
raise RuntimeError("LMI passivity problem is infeasible")
200+
else:
201+
return sol[-1]
202+
203+
204+
def get_input_ff_index(sys):
205+
"""Return the input feedforward passivity (IFP) index for the system.
86206
87-
return (sol["x"] is not None)
207+
The IFP is the largest gain that can be placed in negative parallel
208+
interconnection with a system such that the new interconnected system is
209+
passive.
210+
211+
Parameters
212+
----------
213+
sys : LTI
214+
System to be checked.
215+
216+
Returns
217+
-------
218+
float
219+
The IFP index
220+
"""
221+
sol = solve_passivity_LMI(sys, rho=eps)
222+
if sol is None:
223+
raise RuntimeError("LMI passivity problem is infeasible")
224+
else:
225+
return sol[-1]
226+
227+
228+
def get_relative_index(sys):
229+
"""Return the relative passivity index for the system.
230+
231+
(not implemented yet)
232+
"""
233+
raise NotImplementedError("Relative passivity index not implemented")
234+
235+
236+
def get_combined_io_index(sys):
237+
"""Return the combined I/O passivity index for the system.
238+
239+
(not implemented yet)
240+
"""
241+
raise NotImplementedError("Combined I/O passivity index not implemented")
242+
243+
244+
def get_directional_index(sys):
245+
"""Return the directional passivity index for the system.
246+
247+
(not implemented yet)
248+
"""
249+
raise NotImplementedError("Directional passivity index not implemented")
250+
251+
252+
def ispassive(sys, ofp_index=0, ifp_index=0):
253+
r"""Indicate if a linear time invariant (LTI) system is passive.
254+
255+
Checks if system is passive with the given output feedback (OFP) and input
256+
feedforward (IFP) passivity indices.
257+
258+
Parameters
259+
----------
260+
sys : LTI
261+
System to be checked
262+
ofp_index : float
263+
Output feedback passivity index
264+
ifp_index : float
265+
Input feedforward passivity index
266+
267+
Returns
268+
-------
269+
bool
270+
The system is passive.
271+
272+
Notes
273+
-----
274+
Querying if the system is passive in the sense of
275+
276+
.. math:: V(x) >= 0 \land \dot{V}(x) <= y^T u
277+
278+
is equivalent to the default case of `ofp_index` = 0 and `ifp_index` = 0.
279+
Note that computing the `ofp_index` and `ifp_index` for a system, then
280+
using both values simultaneously as inputs to this function is not
281+
guaranteed to have an output of True (the system might not be passive with
282+
both indices at the same time).
283+
284+
For more details, see [1].
285+
286+
References
287+
----------
288+
.. [1] McCourt, Michael J., and Panos J. Antsaklis
289+
"Demonstrating passivity and dissipativity using computational
290+
methods."
291+
"""
292+
return solve_passivity_LMI(sys, rho=ofp_index, nu=ifp_index) is not None

0 commit comments

Comments
 (0)