Skip to content

Commit 2f29a4c

Browse files
murrayrmmurrayrm
authored
Merge pull request #739 from Mark-Yeatman/passivity-tools
Add passivity module, ispassive function, and passivity_test. Introduces optional dependency cvxopt.
2 parents d98ed29 + c2255b0 commit 2f29a4c

File tree

7 files changed

+195
-7
lines changed

7 files changed

+195
-7
lines changed

‎.github/workflows/python-package-conda.yml

+5-1
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ on: [push, pull_request]
44

55
jobs:
66
test-linux:
7-
name: Python ${{ matrix.python-version }}${{ matrix.slycot && format(' with Slycot from {0}', matrix.slycot) || ' without Slycot' }}${{ matrix.pandas && ', with pandas' || '' }}${{ matrix.array-and-matrix == 1 && ', array and matrix' || '' }}
7+
name: Python ${{ matrix.python-version }}${{ matrix.slycot && format(' with Slycot from {0}', matrix.slycot) || ' without Slycot' }}${{ matrix.pandas && ', with pandas' || '' }}${{ matrix.array-and-matrix == 1 && ', array and matrix' || '' }}${{ matrix.cvxopt && format(' with cvxopt from {0}', matrix.cvxopt) || ' without cvxopt' }}
88
runs-on: ubuntu-latest
99

1010
strategy:
@@ -13,6 +13,7 @@ jobs:
1313
python-version: [3.7, 3.9]
1414
slycot: ["", "conda"]
1515
pandas: [""]
16+
cvxopt: ["", "conda"]
1617
array-and-matrix: [0]
1718
include:
1819
- python-version: 3.9
@@ -45,6 +46,9 @@ jobs:
4546
if [[ '${{matrix.pandas}}' == 'conda' ]]; then
4647
conda install pandas
4748
fi
49+
if [[ '${{matrix.cvxopt}}' == 'conda' ]]; then
50+
conda install -c conda-forge cvxopt
51+
fi
4852
4953
- name: Test with pytest
5054
env:

‎control/exception.py

+12
Original file line numberDiff line numberDiff line change
@@ -84,3 +84,15 @@ def pandas_check():
8484
except:
8585
pandas_installed = False
8686
return pandas_installed
87+
88+
# Utility function to see if cvxopt is installed
89+
cvxopt_installed = None
90+
def cvxopt_check():
91+
global cvxopt_installed
92+
if cvxopt_installed is None:
93+
try:
94+
import cvxopt
95+
cvxopt_installed = True
96+
except:
97+
cvxopt_installed = False
98+
return cvxopt_installed

‎control/lti.py

+13-5
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
"""
1414

1515
import numpy as np
16+
1617
from numpy import absolute, real, angle, abs
1718
from warnings import warn
1819
from . import config
@@ -21,6 +22,7 @@
2122
__all__ = ['poles', 'zeros', 'damp', 'evalfr', 'frequency_response',
2223
'freqresp', 'dcgain', 'pole', 'zero']
2324

25+
2426
class LTI(NamedIOSystem):
2527
"""LTI is a parent class to linear time-invariant (LTI) system objects.
2628
@@ -44,6 +46,7 @@ class LTI(NamedIOSystem):
4446
Note: dt processing has been moved to the NamedIOSystem class.
4547
4648
"""
49+
4750
def __init__(self, inputs=1, outputs=1, states=None, name=None, **kwargs):
4851
"""Assign the LTI object's numbers of inputs and ouputs."""
4952
super().__init__(
@@ -71,8 +74,7 @@ def _set_inputs(self, value):
7174

7275
#: Deprecated
7376
inputs = property(
74-
_get_inputs, _set_inputs, doc=
75-
"""
77+
_get_inputs, _set_inputs, doc="""
7678
Deprecated attribute; use :attr:`ninputs` instead.
7779
7880
The ``inputs`` attribute was used to store the number of system
@@ -94,8 +96,7 @@ def _set_outputs(self, value):
9496

9597
#: Deprecated
9698
outputs = property(
97-
_get_outputs, _set_outputs, doc=
98-
"""
99+
_get_outputs, _set_outputs, doc="""
99100
Deprecated attribute; use :attr:`noutputs` instead.
100101
101102
The ``outputs`` attribute was used to store the number of system
@@ -201,6 +202,11 @@ def _dcgain(self, warn_infinite):
201202
else:
202203
return zeroresp
203204

205+
def ispassive(self):
206+
# importing here prevents circular dependancy
207+
from control.passivity import ispassive
208+
return ispassive(self)
209+
204210
#
205211
# Deprecated functions
206212
#
@@ -321,7 +327,7 @@ def damp(sys, doprint=True):
321327
wn, damping, poles = sys.damp()
322328
if doprint:
323329
print('_____Eigenvalue______ Damping___ Frequency_')
324-
for p, d, w in zip(poles, damping, wn) :
330+
for p, d, w in zip(poles, damping, wn):
325331
if abs(p.imag) < 1e-12:
326332
print("%10.4g %10.4g %10.4g" %
327333
(p.real, 1.0, -p.real))
@@ -330,6 +336,7 @@ def damp(sys, doprint=True):
330336
(p.real, p.imag, d, w))
331337
return wn, damping, poles
332338

339+
333340
def evalfr(sys, x, squeeze=None):
334341
"""Evaluate the transfer function of an LTI system for complex frequency x.
335342
@@ -388,6 +395,7 @@ def evalfr(sys, x, squeeze=None):
388395
"""
389396
return sys.__call__(x, squeeze=squeeze)
390397

398+
391399
def frequency_response(sys, omega, squeeze=None):
392400
"""Frequency response of an LTI system at multiple angular frequencies.
393401

‎control/passivity.py

+87
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
'''
2+
Author: Mark Yeatman
3+
Date: May 15, 2022
4+
'''
5+
6+
import numpy as np
7+
from control import statesp as ss
8+
9+
try:
10+
import cvxopt as cvx
11+
except ImportError as e:
12+
cvx = None
13+
14+
15+
def ispassive(sys):
16+
'''
17+
Indicates if a linear time invariant (LTI) system is passive
18+
19+
Constructs a linear matrix inequality and a feasibility optimization
20+
such that if a solution exists, the system is passive.
21+
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).
25+
26+
Parameters
27+
----------
28+
sys: A continuous LTI system
29+
System to be checked.
30+
31+
Returns
32+
-------
33+
bool:
34+
The input system passive.
35+
'''
36+
if cvx is None:
37+
raise ModuleNotFoundError("cvxopt required for passivity module")
38+
39+
sys = ss._convert_to_statespace(sys)
40+
41+
A = sys.A
42+
B = sys.B
43+
C = sys.C
44+
D = sys.D
45+
46+
# account for strictly proper systems
47+
[n, m] = D.shape
48+
D = D + np.nextafter(0, 1)*np.eye(n, m)
49+
50+
[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+
)
76+
77+
number_of_opt_vars = int(
78+
(state_space_size**2-state_space_size)/2 + state_space_size)
79+
c = cvx.matrix(0.0, (number_of_opt_vars, 1))
80+
81+
# crunch feasibility solution
82+
cvx.solvers.options['show_progress'] = False
83+
sol = cvx.solvers.sdp(c,
84+
Gs=[cvx.matrix(coefficents)],
85+
hs=[cvx.matrix(constants)])
86+
87+
return (sol["x"] is not None)

‎control/tests/conftest.py

+2
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,8 @@
1717
# pytest.param(marks=)
1818
slycotonly = pytest.mark.skipif(not control.exception.slycot_check(),
1919
reason="slycot not installed")
20+
cvxoptonly = pytest.mark.skipif(not control.exception.cvxopt_check(),
21+
reason="cvxopt not installed")
2022
matrixfilter = pytest.mark.filterwarnings("ignore:.*matrix subclass:"
2123
"PendingDeprecationWarning")
2224
matrixerrorfilter = pytest.mark.filterwarnings("error:.*matrix subclass:"

‎control/tests/passivity_test.py

+74
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
'''
2+
Author: Mark Yeatman
3+
Date: May 30, 2022
4+
'''
5+
import pytest
6+
import numpy
7+
from control import ss, passivity, tf
8+
from control.tests.conftest import cvxoptonly
9+
10+
11+
pytestmark = cvxoptonly
12+
13+
14+
def test_ispassive():
15+
A = numpy.array([[0, 1], [-2, -2]])
16+
B = numpy.array([[0], [1]])
17+
C = numpy.array([[-1, 2]])
18+
D = numpy.array([[1.5]])
19+
sys = ss(A, B, C, D)
20+
21+
# happy path is passive
22+
assert(passivity.ispassive(sys))
23+
24+
# happy path not passive
25+
D = -D
26+
sys = ss(A, B, C, D)
27+
28+
assert(not passivity.ispassive(sys))
29+
30+
31+
A_d = numpy.array([[-2, 0], [0, 0]])
32+
A = numpy.array([[-3, 0], [0, -2]])
33+
B = numpy.array([[0], [1]])
34+
C = numpy.array([[-1, 2]])
35+
D = numpy.array([[1.5]])
36+
37+
38+
@pytest.mark.parametrize(
39+
"test_input,expected",
40+
[((A, B, C, D*0.0), True),
41+
((A_d, B, C, D), True),
42+
((A*1e12, B, C, D*0), True),
43+
((A, B*0, C*0, D), True),
44+
((A*0, B, C, D), True),
45+
((A*0, B*0, C*0, D*0), True)])
46+
def test_ispassive_edge_cases(test_input, expected):
47+
48+
# strictly proper
49+
A = test_input[0]
50+
B = test_input[1]
51+
C = test_input[2]
52+
D = test_input[3]
53+
sys = ss(A, B, C, D)
54+
assert(passivity.ispassive(sys) == expected)
55+
56+
57+
def test_transfer_function():
58+
sys = tf([1], [1, 2])
59+
assert(passivity.ispassive(sys))
60+
61+
sys = tf([1], [1, -2])
62+
assert(not passivity.ispassive(sys))
63+
64+
65+
def test_oo_style():
66+
A = numpy.array([[0, 1], [-2, -2]])
67+
B = numpy.array([[0], [1]])
68+
C = numpy.array([[-1, 2]])
69+
D = numpy.array([[1.5]])
70+
sys = ss(A, B, C, D)
71+
assert(sys.ispassive())
72+
73+
sys = tf([1], [1, 2])
74+
assert(sys.ispassive())

‎setup.py

+2-1
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,7 @@
4646
'matplotlib'],
4747
extras_require={
4848
'test': ['pytest', 'pytest-timeout'],
49-
'slycot': [ 'slycot>=0.4.0' ]
49+
'slycot': [ 'slycot>=0.4.0' ],
50+
'cvxopt': [ 'cvxopt>=1.2.0' ]
5051
}
5152
)

0 commit comments

Comments
 (0)