murrayrm · Jun 16, 2022 · May 31, 2022 · Jun 6, 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
@@ -46,6 +47,9 @@ jobs:
         if [[ '${{matrix.pandas}}' == 'conda' ]]; then
           conda install -c conda-forge 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/passivity.py b/control/passivity.py
@@ -0,0 +1,81 @@
+'''
+Author: Mark Yeatman  
+Date: May 15, 2022
+'''
+
+import numpy as np
+
+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")
+
+    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)
+
+    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
@@ -18,6 +18,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")
 noscipy0 = pytest.mark.skipif(StrictVersion(sp.__version__) < "1.0",
                               reason="requires SciPy 1.0 or greater")
 nopython2 = pytest.mark.skipif(sys.version_info < (3, 0),

diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
@@ -0,0 +1,63 @@
+'''
+Author: Mark Yeatman  
+Date: May 30, 2022
+'''
+
+import numpy
+from control import ss, passivity
+from control.tests.conftest import cvxoptonly
+
+
+@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))
+
+
+@cvxoptonly
+def test_ispassive_edge_cases():
+    A = numpy.array([[0, 1], [-2, -2]])
+    B = numpy.array([[0], [1]])
+    C = numpy.array([[-1, 2]])
+    D = numpy.array([[1.5]])
+
+    D *= 0
+
+    # strictly proper
+    sys = ss(A, B, C, D)
+    assert(passivity.ispassive(sys))
+
+    # ill conditioned
+    A = A*1e12
+    sys = ss(A, B, C, D)
+    assert(passivity.ispassive(sys))
+
+    # different combinations of zero A,B,C,D are 0
+    B *= 0
+    C *= 0
+    assert(passivity.ispassive(sys))
+
+    A *= 0
+    B = numpy.array([[0], [1]])
+    C = numpy.array([[-1, 2]])
+    D = numpy.array([[1.5]])
+    assert(passivity.ispassive(sys))
+
+    B *= 0
+    C *= 0
+    assert(passivity.ispassive(sys))
+
+    A *= 0
+    assert(passivity.ispassive(sys))
diff --git a/setup.py b/setup.py
@@ -43,6 +43,7 @@
                       'matplotlib'],
     extras_require={
        'test': ['pytest', 'pytest-timeout'],
-       'slycot': [ 'slycot>=0.4.0' ]
+       'slycot': [ 'slycot>=0.4.0' ],
+       'cvxopt': [ 'cvxopt>=1.2.0' ]
     }
 )