-
Notifications
You must be signed in to change notification settings - Fork 453
Expand file tree
/
Copy pathpassivity.py
More file actions
303 lines (234 loc) · 8.92 KB
/
passivity.py
File metadata and controls
303 lines (234 loc) · 8.92 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
# passivity.py - functions for passive control
#
# Initial author: Mark Yeatman
# Creation date: July 17, 2022
"""Functions for passive control."""
import numpy as np
from control import statesp
from control.exception import ControlArgument, ControlDimension
try:
import cvxopt as cvx
except ImportError:
cvx = None
__all__ = ["get_output_fb_index", "get_input_ff_index", "ispassive",
"solve_passivity_LMI"]
def solve_passivity_LMI(sys, rho=None, nu=None):
"""Compute passivity indices and/or solves feasibility via a LMI.
Constructs a linear matrix inequality (LMI) such that if a solution
exists and the last element of the solution is positive, the system
`sys` is passive. Inputs of None for `rho` or `nu` indicate that the
function should solve for that index (they are mutually exclusive, they
can't both be None, otherwise you're trying to solve a nonconvex
bilinear matrix inequality.) The last element of the output `solution`
is either the output or input passivity index, for `rho` = None and
`nu` = None, respectively.
Parameters
----------
sys : LTI
System to be checked.
rho : float or None
Output feedback passivity index.
nu : float or None
Input feedforward passivity index.
Returns
-------
solution : ndarray
The LMI solution.
References
----------
.. [1] McCourt, Michael J., and Panos J. Antsaklis, "Demonstrating
passivity and dissipativity using computational methods."
.. [2] Nicholas Kottenstette and Panos J. Antsaklis,
"Relationships Between Positive Real, Passive Dissipative, &
Positive Systems", equation 36.
"""
if cvx is None:
raise ModuleNotFoundError("cvxopt required for passivity module")
if sys.ninputs != sys.noutputs:
raise ControlDimension(
"The number of system inputs must be the same as the number of "
"system outputs.")
if rho is None and nu is None:
raise ControlArgument("rho or nu must be given a numerical value.")
sys = statesp._convert_to_statespace(sys)
A = sys.A
B = sys.B
C = sys.C
D = sys.D
# account for strictly proper systems
[_, m] = D.shape
[n, _] = A.shape
def make_LMI_matrix(P, rho, nu, one):
q = sys.noutputs
Q = -rho*np.eye(q, q)
S = 1.0/2.0*(one+rho*nu)*np.eye(q)
R = -nu*np.eye(m)
if sys.isctime():
off_diag = P@B - (C.T@S + C.T@Q@D)
return np.vstack((
np.hstack((A.T @ P + P@A - C.T@Q@C, off_diag)),
np.hstack((off_diag.T, -(D.T@Q@D + D.T@S + S.T@D + R)))
))
else:
off_diag = A.T@P@B - (C.T@S + C.T@Q@D)
return np.vstack((
np.hstack((A.T @ P @ A - P - C.T@Q@C, off_diag)),
np.hstack((off_diag.T, B.T@P@B-(D.T@Q@D + D.T@S + S.T@D + R)))
))
def make_P_basis_matrices(n, rho, nu):
"""Make list of matrix constraints for passivity LMI.
Utility function to make basis matrices for a LMI from a
symmetric matrix P of size n by n representing a parameterized
symbolic matrix.
"""
matrix_list = []
for i in range(0, n):
for j in range(0, n):
if j <= i:
P = np.zeros((n, n))
P[i, j] = 1
P[j, i] = 1
matrix_list.append(make_LMI_matrix(P, 0, 0, 0).flatten())
zeros = 0.0*np.eye(n)
if rho is None:
matrix_list.append(make_LMI_matrix(zeros, 1, 0, 0).flatten())
elif nu is None:
matrix_list.append(make_LMI_matrix(zeros, 0, 1, 0).flatten())
return matrix_list
def P_pos_def_constraint(n):
"""Make a list of matrix constraints for P >= 0.
Utility function to make basis matrices for a LMI that ensures
parameterized symbolic matrix of size n by n is positive definite
"""
matrix_list = []
for i in range(0, n):
for j in range(0, n):
if j <= i:
P = np.zeros((n, n))
P[i, j] = -1
P[j, i] = -1
matrix_list.append(P.flatten())
if rho is None or nu is None:
matrix_list.append(np.zeros((n, n)).flatten())
return matrix_list
n = sys.nstates
# coefficients for passivity indices and feasibility matrix
sys_matrix_list = make_P_basis_matrices(n, rho, nu)
# get constants for numerical values of rho and nu
sys_constants = list()
if rho is not None and nu is not None:
sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, nu, 1.0)
elif rho is not None:
sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, 0.0, 1.0)
elif nu is not None:
sys_constants = -make_LMI_matrix(np.zeros_like(A), 0.0, nu, 1.0)
sys_coefficents = np.vstack(sys_matrix_list).T
# LMI to ensure P is positive definite
P_matrix_list = P_pos_def_constraint(n)
P_coefficents = np.vstack(P_matrix_list).T
P_constants = np.zeros((n, n))
# cost function
number_of_opt_vars = int(
(n**2-n)/2 + n)
c = cvx.matrix(0.0, (number_of_opt_vars, 1))
#we're maximizing a passivity index, include it in the cost function
if rho is None or nu is None:
c = cvx.matrix(np.append(np.array(c), -1.0))
Gs = [cvx.matrix(sys_coefficents)] + [cvx.matrix(P_coefficents)]
hs = [cvx.matrix(sys_constants)] + [cvx.matrix(P_constants)]
# crunch feasibility solution
cvx.solvers.options['show_progress'] = False
try:
sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
return sol["x"]
except ZeroDivisionError as e:
raise ValueError(
"The system is probably ill conditioned. Consider perturbing "
"the system matrices by a small amount."
) from e
def get_output_fb_index(sys):
"""Return the output feedback passivity (OFP) index for the system.
The OFP is the largest gain that can be placed in positive feedback
with a system such that the new interconnected system is passive.
Parameters
----------
sys : LTI
System to be checked.
Returns
-------
float
The OFP index.
"""
sol = solve_passivity_LMI(sys, nu=0.0)
if sol is None:
raise RuntimeError("LMI passivity problem is infeasible")
else:
return sol[-1]
def get_input_ff_index(sys):
"""Input feedforward passivity (IFP) index for a system.
The input feedforward passivity (IFP) is the largest gain that can be
placed in negative parallel interconnection with a system such that the
new interconnected system is passive.
Parameters
----------
sys : LTI
System to be checked.
Returns
-------
float
The IFP index.
"""
sol = solve_passivity_LMI(sys, rho=0.0)
if sol is None:
raise RuntimeError("LMI passivity problem is infeasible")
else:
return sol[-1]
def get_relative_index(sys):
"""Return the relative passivity index for the system.
(not implemented yet)
"""
raise NotImplementedError("Relative passivity index not implemented")
def get_combined_io_index(sys):
"""Return the combined I/O passivity index for the system.
(not implemented yet)
"""
raise NotImplementedError("Combined I/O passivity index not implemented")
def get_directional_index(sys):
"""Return the directional passivity index for the system.
(not implemented yet)
"""
raise NotImplementedError("Directional passivity index not implemented")
def ispassive(sys, ofp_index=0, ifp_index=0):
r"""Indicate if a linear time invariant (LTI) system is passive.
Checks if system is passive with the given output feedback (OFP)
and input feedforward (IFP) passivity indices.
Parameters
----------
sys : LTI
System to be checked.
ofp_index : float
Output feedback passivity index.
ifp_index : float
Input feedforward passivity index.
Returns
-------
bool
The system is passive.
Notes
-----
Querying if the system is passive in the sense of
.. math:: V(x) >= 0 \land \dot{V}(x) <= y^T u
is equivalent to the default case of `ofp_index` = 0 and `ifp_index` =
0. Note that computing the `ofp_index` and `ifp_index` for a system,
then using both values simultaneously as inputs to this function is not
guaranteed to have an output of True (the system might not be passive
with both indices at the same time).
For more details, see [1]_.
References
----------
.. [1] McCourt, Michael J., and Panos J. Antsaklis "Demonstrating
passivity and dissipativity using computational methods."
Technical Report of the ISIS Group at the University of Notre
Dame. ISIS-2013-008, Aug. 2013.
"""
return solve_passivity_LMI(sys, rho=ofp_index, nu=ifp_index) is not None