python-control · murrayrm · Jun 16, 2022 · May 31, 2022 · Jun 6, 2022 · Jun 6, 2022
diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml
@@ -4,7 +4,7 @@ on: [push, pull_request]
 
 jobs:
   test-linux:
-    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' || '' }}
+    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' }}
     runs-on: ubuntu-latest
 
     strategy:
@@ -13,6 +13,7 @@ jobs:
         python-version: [3.7, 3.9]
         slycot: ["", "conda"]
         pandas: [""]
+        cvxopt: ["", "conda"]
         array-and-matrix: [0]
         include:
           - python-version: 3.9
@@ -45,6 +46,9 @@ jobs:
         if [[ '${{matrix.pandas}}' == 'conda' ]]; then
           conda install pandas
         fi
+        if [[ '${{matrix.cvxopt}}' == 'conda' ]]; then
+          conda install -c conda-forge cvxopt
+        fi
 
     - name: Test with pytest
       env:

diff --git a/control/exception.py b/control/exception.py
@@ -84,3 +84,15 @@ def pandas_check():
         except:
             pandas_installed = False
     return pandas_installed
+
+# Utility function to see if cvxopt is installed
+cvxopt_installed = None
+def cvxopt_check():
+    global cvxopt_installed
+    if cvxopt_installed is None:
+        try:
+            import cvxopt
+            cvxopt_installed = True
+        except:
+            cvxopt_installed = False
+    return cvxopt_installed
diff --git a/control/lti.py b/control/lti.py
@@ -13,6 +13,7 @@
 """
 
 import numpy as np
+
 from numpy import absolute, real, angle, abs
 from warnings import warn
 from . import config
@@ -21,6 +22,7 @@
 __all__ = ['poles', 'zeros', 'damp', 'evalfr', 'frequency_response',
            'freqresp', 'dcgain', 'pole', 'zero']
 
+
 class LTI(NamedIOSystem):
     """LTI is a parent class to linear time-invariant (LTI) system objects.
 
@@ -44,6 +46,7 @@ class LTI(NamedIOSystem):
     Note: dt processing has been moved to the NamedIOSystem class.
 
     """
+
     def __init__(self, inputs=1, outputs=1, states=None, name=None, **kwargs):
         """Assign the LTI object's numbers of inputs and ouputs."""
         super().__init__(
@@ -71,8 +74,7 @@ def _set_inputs(self, value):
 
     #: Deprecated
     inputs = property(
-        _get_inputs, _set_inputs, doc=
-        """
+        _get_inputs, _set_inputs, doc="""
         Deprecated attribute; use :attr:`ninputs` instead.
 
         The ``inputs`` attribute was used to store the number of system
@@ -94,8 +96,7 @@ def _set_outputs(self, value):
 
     #: Deprecated
     outputs = property(
-        _get_outputs, _set_outputs, doc=
-        """
+        _get_outputs, _set_outputs, doc="""
         Deprecated attribute; use :attr:`noutputs` instead.
 
         The ``outputs`` attribute was used to store the number of system
@@ -201,6 +202,11 @@ def _dcgain(self, warn_infinite):
         else:
             return zeroresp
 
+    def ispassive(self):
+        # importing here prevents circular dependancy
+        from control.passivity import ispassive
+        return ispassive(self)
+
     #
     # Deprecated functions
     #
@@ -321,7 +327,7 @@ def damp(sys, doprint=True):
     wn, damping, poles = sys.damp()
     if doprint:
         print('_____Eigenvalue______ Damping___ Frequency_')
-        for p, d, w in zip(poles, damping, wn) :
+        for p, d, w in zip(poles, damping, wn):
             if abs(p.imag) < 1e-12:
                 print("%10.4g            %10.4g %10.4g" %
                       (p.real, 1.0, -p.real))
@@ -330,6 +336,7 @@ def damp(sys, doprint=True):
                       (p.real, p.imag, d, w))
     return wn, damping, poles
 
+
 def evalfr(sys, x, squeeze=None):
     """Evaluate the transfer function of an LTI system for complex frequency x.
 
@@ -388,6 +395,7 @@ def evalfr(sys, x, squeeze=None):
     """
     return sys.__call__(x, squeeze=squeeze)
 
+
 def frequency_response(sys, omega, squeeze=None):
     """Frequency response of an LTI system at multiple angular frequencies.
 

diff --git a/control/passivity.py b/control/passivity.py
@@ -0,0 +1,87 @@
+'''
+Author: Mark Yeatman  
+Date: May 15, 2022
+'''
+
+import numpy as np
+from control import statesp as ss
+
+try:
+    import cvxopt as cvx
+except ImportError as e:
+    cvx = None
+
+
+def ispassive(sys):
+    '''
+    Indicates if a linear time invariant (LTI) system is passive
+
+    Constructs a linear matrix inequality and a feasibility optimization
+    such that if a solution exists, the system is passive.
+
+    The source for the algorithm is: 
+    McCourt, Michael J., and Panos J. Antsaklis. 
+        "Demonstrating passivity and dissipativity using computational methods." ISIS 8 (2013).
+
+    Parameters
+    ----------
+    sys: A continuous LTI system
+        System to be checked.
+
+    Returns
+    -------
+    bool: 
+        The input system passive.
+    '''
+    if cvx is None:
+        raise ModuleNotFoundError("cvxopt required for passivity module")
+
+    sys = ss._convert_to_statespace(sys)
+
+    A = sys.A
+    B = sys.B
+    C = sys.C
+    D = sys.D
+
+    # account for strictly proper systems
+    [n, m] = D.shape
+    D = D + np.nextafter(0, 1)*np.eye(n, m)
+
+    [n, _] = A.shape
+    A = A - np.nextafter(0, 1)*np.eye(n)
+
+    def make_LMI_matrix(P):
+        V = np.vstack((
+            np.hstack((A.T @ P + P@A, P@B)),
+            np.hstack((B.T@P, np.zeros_like(D))))
+        )
+        return V
+
+    matrix_list = []
+    state_space_size = sys.nstates
+    for i in range(0, state_space_size):
+        for j in range(0, state_space_size):
+            if j <= i:
+                P = np.zeros_like(A)
+                P[i, j] = 1.0
+                P[j, i] = 1.0
+                matrix_list.append(make_LMI_matrix(P).flatten())
+
+    coefficents = np.vstack(matrix_list).T
+
+    constants = -np.vstack((
+        np.hstack((np.zeros_like(A),  - C.T)),
+        np.hstack((- C, -D - D.T)))
+    )
+
+    number_of_opt_vars = int(
+        (state_space_size**2-state_space_size)/2 + state_space_size)
+    c = cvx.matrix(0.0, (number_of_opt_vars, 1))
+
+    # crunch feasibility solution
+    cvx.solvers.options['show_progress'] = False
+    sol = cvx.solvers.sdp(c,
+                          Gs=[cvx.matrix(coefficents)],
+                          hs=[cvx.matrix(constants)])
+
+    return (sol["x"] is not None)
diff --git a/control/tests/conftest.py b/control/tests/conftest.py
@@ -17,6 +17,8 @@
 # pytest.param(marks=)
 slycotonly = pytest.mark.skipif(not control.exception.slycot_check(),
                                 reason="slycot not installed")
+cvxoptonly = pytest.mark.skipif(not control.exception.cvxopt_check(),
+                                reason="cvxopt not installed")
 matrixfilter = pytest.mark.filterwarnings("ignore:.*matrix subclass:"
                                           "PendingDeprecationWarning")
 matrixerrorfilter = pytest.mark.filterwarnings("error:.*matrix subclass:"

diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
@@ -0,0 +1,74 @@
+'''
+Author: Mark Yeatman  
+Date: May 30, 2022
+'''
+import pytest
+import numpy
+from control import ss, passivity, tf
+from control.tests.conftest import cvxoptonly
+
+
+pytestmark = cvxoptonly
+
+
+def test_ispassive():
+    A = numpy.array([[0, 1], [-2, -2]])
+    B = numpy.array([[0], [1]])
+    C = numpy.array([[-1, 2]])
+    D = numpy.array([[1.5]])
+    sys = ss(A, B, C, D)
+
+    # happy path is passive
+    assert(passivity.ispassive(sys))
+
+    # happy path not passive
+    D = -D
+    sys = ss(A, B, C, D)
+
+    assert(not passivity.ispassive(sys))
+
+
+A_d = numpy.array([[-2, 0], [0, 0]])
+A = numpy.array([[-3, 0], [0, -2]])
+B = numpy.array([[0], [1]])
+C = numpy.array([[-1, 2]])
+D = numpy.array([[1.5]])
+
+
+@pytest.mark.parametrize(
+    "test_input,expected",
+    [((A, B, C, D*0.0), True),
+     ((A_d, B, C, D), True),
+     ((A*1e12, B, C, D*0), True),
+     ((A, B*0, C*0, D), True),
+     ((A*0, B, C, D), True),
+     ((A*0, B*0, C*0, D*0), True)])
+def test_ispassive_edge_cases(test_input, expected):
+
+    # strictly proper
+    A = test_input[0]
+    B = test_input[1]
+    C = test_input[2]
+    D = test_input[3]
+    sys = ss(A, B, C, D)
+    assert(passivity.ispassive(sys) == expected)
+
+
+def test_transfer_function():
+    sys = tf([1], [1, 2])
+    assert(passivity.ispassive(sys))
+
+    sys = tf([1], [1, -2])
+    assert(not passivity.ispassive(sys))
+
+
+def test_oo_style():
+    A = numpy.array([[0, 1], [-2, -2]])
+    B = numpy.array([[0], [1]])
+    C = numpy.array([[-1, 2]])
+    D = numpy.array([[1.5]])
+    sys = ss(A, B, C, D)
+    assert(sys.ispassive())
+
+    sys = tf([1], [1, 2])
+    assert(sys.ispassive())
diff --git a/setup.py b/setup.py
@@ -46,6 +46,7 @@
                       'matplotlib'],
     extras_require={
        'test': ['pytest', 'pytest-timeout'],
-       'slycot': [ 'slycot>=0.4.0' ]
+       'slycot': [ 'slycot>=0.4.0' ],
+       'cvxopt': [ 'cvxopt>=1.2.0' ]
     }
 )