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
@@ -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,61 @@
+'''
+Author: Mark Yeatman  
+Date: May 15, 2022
+'''
+
+from . import statesp as ss
+import numpy as np
+import cvxopt as cvx
+
+
+def is_passive(sys):
+    '''
+    Indicates if a linear time invarient system is passive
+
+    Constructs a linear matrix inequality and a feasibility optimization
+    such that is 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).
+    '''
+
+    A = sys.A
+    B = sys.B
+    C = sys.C
+    D = sys.D
+
+    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
+
+    P = np.zeros_like(A)
+    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 = P*0.0
+                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
+    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,25 @@
+'''
+Author: Mark Yeatman  
+Date: May 30, 2022
+'''
+
+import numpy
+from control import ss, passivity
+from sympy import Matrix
+from control.tests.conftest import cvxoptonly
+
+@cvxoptonly
+def test_is_passive():
+    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(passivity.is_passive(sys))
+
+    D = -D
+    sys = ss(A, B, C, D)
+
+    assert(not passivity.is_passive(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' ]
     }
 )