From 8e40e0f7c2c5f261e854485c2d276e12946d75f5 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Sat, 18 Jun 2022 00:50:19 -0400
Subject: [PATCH 01/36] Add support for discrete time system to ispassive.

---
 control/passivity.py            | 29 ++++++++++++++++++++---------
 control/tests/passivity_test.py | 20 ++++++++++++++++++--
 2 files changed, 38 insertions(+), 11 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index e833fbf96..a27b67be7 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -19,13 +19,18 @@ def ispassive(sys):
     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).
+    The sources for the algorithm are: 
+    
+    McCourt, Michael J., and Panos J. Antsaklis
+        "Demonstrating passivity and dissipativity using computational methods." 
+
+    Nicholas Kottenstette and Panos J. Antsaklis
+        "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
+        equation 36.
 
     Parameters
     ----------
-    sys: A continuous LTI system
+    sys: An LTI system
         System to be checked.
 
     Returns
@@ -51,11 +56,17 @@ def ispassive(sys):
     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
+        if sys.isctime():
+            return np.vstack((
+                np.hstack((A.T @ P + P@A, P@B)),
+                np.hstack((B.T@P, np.zeros_like(D))))
+            )
+        else:
+            return np.vstack((
+                np.hstack((2*A.T @ P @A - P , 2*A.T @ P@B)),
+                np.hstack(((2*A.T @ P@B).T, 2*B.T@P@B)))
+            )
+        
 
     matrix_list = []
     state_space_size = sys.nstates
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index 791d70b6c..43f6d134d 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -4,14 +4,14 @@
 '''
 import pytest
 import numpy
-from control import ss, passivity, tf
+from control import ss, passivity, tf, sample_system
 from control.tests.conftest import cvxoptonly
 
 
 pytestmark = cvxoptonly
 
 
-def test_ispassive():
+def test_ispassive_ctime():
     A = numpy.array([[0, 1], [-2, -2]])
     B = numpy.array([[0], [1]])
     C = numpy.array([[-1, 2]])
@@ -27,6 +27,22 @@ def test_ispassive():
 
     assert(not passivity.ispassive(sys))
 
+def test_ispassive_dtime():
+    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)
+    sys = sample_system(sys, 1, method='bilinear')
+    # 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]])

From 7d4b5f242e952b2508865812fcc8aa171e6455c7 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Sun, 10 Jul 2022 22:05:35 -0400
Subject: [PATCH 02/36] Save some work on implementing passivity indices, my
 source for the algorithm only considers square, minimal systems.

---
 control/passivity.py | 214 +++++++++++++++++++++++++++++++++++--------
 1 file changed, 174 insertions(+), 40 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index a27b67be7..874e04ddd 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -3,6 +3,7 @@
 Date: May 15, 2022
 '''
 
+from msilib.schema import Error
 import numpy as np
 from control import statesp as ss
 
@@ -11,8 +12,126 @@
 except ImportError as e:
     cvx = None
 
+def __parse_lti__(sys):
+    '''
+    Utility function to parse LTI input for passivity module functions
+    '''
+    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)
+
+    return (A, B, C, D)
+
+def __make_P_basis_matrices__(n, make_LMI_matrix_func):
+    '''
+    Utility function to make basis matrices for a LMI from a 
+    functional make_LMI_matrix_func and a symmetric matrix P of size n by n
+    representing a parametrized 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.0
+                P[j, i] = 1.0
+                matrix_list.append(make_LMI_matrix_func(P).flatten())
+    return matrix_list
+
+def __P_pos_def_constraint__(n):
+    '''
+    Utility function to make basis matrices for a LMI that ensures parametrized 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.0
+                P[j, i] = -1.0
+                matrix_list.append(P.flatten())
+    return matrix_list
 
-def ispassive(sys):
+def __solve_LMI_problem__(A, C, D, nu, rho, n, make_LMI_matrix):
+    '''
+    Construct LMI for constants in passivitiy problem, format matrices and input into SDP solver. 
+    '''
+    #LMI for passivity from A,B,C,D   
+    sys_matrix_list =__make_P_basis_matrices__(n, make_LMI_matrix)
+
+    sys_coefficents = np.vstack(sys_matrix_list).T
+
+    sys_constants = -np.vstack((
+        np.hstack((np.zeros_like(A),  - C.T)),
+        np.hstack((- C, -D - D.T)))
+    )
+
+    # if nu is not None:
+    #     sys_constants+= -np.vstack((
+    #         np.hstack((np.zeros_like(A),  np.zeros_like(C.T))),
+    #         np.hstack((np.zeros_like(C),  nu*np.eye(n) )))
+    #     )
+
+    # if rho is not None:
+    #     sys_constants+= -np.vstack((
+    #         np.hstack((rho*C.T@C,  rho*C.T@D )),
+    #         np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
+    #     )
+
+    # if rho is not None and nu is not None:
+    #     sys_constants+= -np.vstack((
+    #         np.hstack((np.zeros_like(A),  -0.5*nu*rho*C.T)),
+    #         np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
+    #     )
+
+    #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))
+
+    # #LMI for passivity indices
+    # rho_coefficents = []
+    # nu_coefficents = []
+    # if nu is not None and rho is None:
+    #     #pick out coefficents for rho
+    #     rho_coefficents_matrix = [np.vstack((
+    #             np.hstack((C.T@C, 0.5*nu*C.T + C.T@D)),
+    #             np.hstack(( (0.5*nu*C.T + C.T@D).T, D.T@D-nu*(D+D.T))))
+    #         ).flatten()]
+    #     rho_coefficents =np.vstack(rho_coefficents_matrix).T
+    # elif rho is not None and nu is None:
+    #     #pick out coefficents for nu
+    #     nu_coefficents_matrix = [np.vstack((
+    #             np.hstack((np.zeros_like(A), 0.5*rho*C.T)),
+    #             np.hstack(( (0.5*rho*C.T + rho*C.T@D).T, rho*D.T@D-nu*(D+D.T))))
+    #         ).flatten()]
+
+    number_of_opt_vars = int(
+        (n**2-n)/2 + n)
+    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(sys_coefficents)]+[cvx.matrix(P_coefficents)],
+                          hs=[cvx.matrix(sys_constants)]+[cvx.matrix(P_constants)])
+
+    return (sol["x"] is not None)
+
+
+def ispassive(sys, nu = None, rho = None):
     '''
     Indicates if a linear time invariant (LTI) system is passive
 
@@ -41,19 +160,10 @@ def ispassive(sys):
     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)
+    if not sys.isctime() and rho is not None and nu is not None:
+        raise Error("Passivity indices for discrete time systems not supported yet.")
 
-    [n, _] = A.shape
-    A = A - np.nextafter(0, 1)*np.eye(n)
+    (A,B,C,D) = __parse_lti__(sys)
 
     def make_LMI_matrix(P):
         if sys.isctime():
@@ -62,37 +172,61 @@ def make_LMI_matrix(P):
                 np.hstack((B.T@P, np.zeros_like(D))))
             )
         else:
-            return np.vstack((
-                np.hstack((2*A.T @ P @A - P , 2*A.T @ P@B)),
-                np.hstack(((2*A.T @ P@B).T, 2*B.T@P@B)))
+            return 2*np.vstack((
+                np.hstack((A.T @ P @A - P , A.T @ P@B)),
+                np.hstack(((A.T @ P@B).T, B.T@P@B)))
             )
-        
 
-    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())
+    return __solve_LMI_problem__(A, C, D, None, None, sys.nstates, make_LMI_matrix)
 
-    coefficents = np.vstack(matrix_list).T
+# def isQSRdissapative(sys):
+#     '''
+#     Indicates if a linear time invariant (LTI) system is dissapative
 
-    constants = -np.vstack((
-        np.hstack((np.zeros_like(A),  - C.T)),
-        np.hstack((- C, -D - D.T)))
-    )
+#     Constructs a linear matrix inequality and a feasibility optimization
+#     such that if a solution exists, the system is dissapative.
 
-    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))
+#     The sources for the algorithm are: 
+    
+#     McCourt, Michael J., and Panos J. Antsaklis
+#         "Demonstrating passivity and dissipativity using computational methods." 
+
+#     Nicholas Kottenstette and Panos J. Antsaklis
+#         "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
+#         equation 36.
+
+#     Parameters
+#     ----------
+#     sys: An LTI system
+#         System to be checked.
+
+#     Returns
+#     -------
+#     bool: 
+#         The input system dissapative.
+#     '''
+#     if cvx is None:
+#         raise ModuleNotFoundError("cvxopt required for passivity module")
+
+#     (A,B,C,D) = _parse_lti(sys)
+
+#     def make_LMI_matrix(P):
+#         None
+        
+#     matrix_list =_make_basis_matrices(sys.nstates, make_LMI_matrix)
 
-    # crunch feasibility solution
-    cvx.solvers.options['show_progress'] = False
-    sol = cvx.solvers.sdp(c,
-                          Gs=[cvx.matrix(coefficents)],
-                          hs=[cvx.matrix(constants)])
+#     coefficents = np.vstack(matrix_list).T
 
-    return (sol["x"] is not None)
+#     constants = None
+
+#     number_of_opt_vars = int(
+#         (sys.nstates**2-sys.nstates)/2 + sys.nstates)
+#     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)
\ No newline at end of file

From 941b67367846d3e1e52f8f39129952976d8eef26 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 11 Jul 2022 14:59:58 -0400
Subject: [PATCH 03/36] Catch cases with dimension mismatch between input and
 outputs.

---
 control/passivity.py            | 48 ++++++++++++++++++++-------------
 control/tests/passivity_test.py | 12 +++++++++
 2 files changed, 41 insertions(+), 19 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 874e04ddd..974fdb8d3 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -12,6 +12,7 @@
 except ImportError as e:
     cvx = None
 
+
 def __parse_lti__(sys):
     '''
     Utility function to parse LTI input for passivity module functions
@@ -32,6 +33,7 @@ def __parse_lti__(sys):
 
     return (A, B, C, D)
 
+
 def __make_P_basis_matrices__(n, make_LMI_matrix_func):
     '''
     Utility function to make basis matrices for a LMI from a 
@@ -42,12 +44,13 @@ def __make_P_basis_matrices__(n, make_LMI_matrix_func):
     for i in range(0, n):
         for j in range(0, n):
             if j <= i:
-                P = np.zeros((n,n))
+                P = np.zeros((n, n))
                 P[i, j] = 1.0
                 P[j, i] = 1.0
                 matrix_list.append(make_LMI_matrix_func(P).flatten())
     return matrix_list
 
+
 def __P_pos_def_constraint__(n):
     '''
     Utility function to make basis matrices for a LMI that ensures parametrized symbolic matrix 
@@ -57,18 +60,19 @@ def __P_pos_def_constraint__(n):
     for i in range(0, n):
         for j in range(0, n):
             if j <= i:
-                P = np.zeros((n,n))
+                P = np.zeros((n, n))
                 P[i, j] = -1.0
                 P[j, i] = -1.0
                 matrix_list.append(P.flatten())
     return matrix_list
 
+
 def __solve_LMI_problem__(A, C, D, nu, rho, n, make_LMI_matrix):
     '''
     Construct LMI for constants in passivitiy problem, format matrices and input into SDP solver. 
     '''
-    #LMI for passivity from A,B,C,D   
-    sys_matrix_list =__make_P_basis_matrices__(n, make_LMI_matrix)
+    # LMI for passivity from A,B,C,D
+    sys_matrix_list = __make_P_basis_matrices__(n, make_LMI_matrix)
 
     sys_coefficents = np.vstack(sys_matrix_list).T
 
@@ -95,8 +99,8 @@ def __solve_LMI_problem__(A, C, D, nu, rho, n, make_LMI_matrix):
     #         np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
     #     )
 
-    #LMI to ensure P is positive definite
-    P_matrix_list =__P_pos_def_constraint__(n)
+    # 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))
@@ -125,13 +129,14 @@ def __solve_LMI_problem__(A, C, D, nu, rho, n, make_LMI_matrix):
     # crunch feasibility solution
     cvx.solvers.options['show_progress'] = False
     sol = cvx.solvers.sdp(c,
-                          Gs=[cvx.matrix(sys_coefficents)]+[cvx.matrix(P_coefficents)],
+                          Gs=[cvx.matrix(sys_coefficents)] +
+                          [cvx.matrix(P_coefficents)],
                           hs=[cvx.matrix(sys_constants)]+[cvx.matrix(P_constants)])
 
     return (sol["x"] is not None)
 
 
-def ispassive(sys, nu = None, rho = None):
+def ispassive(sys, nu=None, rho=None):
     '''
     Indicates if a linear time invariant (LTI) system is passive
 
@@ -139,7 +144,7 @@ def ispassive(sys, nu = None, rho = None):
     such that if a solution exists, the system is passive.
 
     The sources for the algorithm are: 
-    
+
     McCourt, Michael J., and Panos J. Antsaklis
         "Demonstrating passivity and dissipativity using computational methods." 
 
@@ -161,9 +166,14 @@ def ispassive(sys, nu = None, rho = None):
         raise ModuleNotFoundError("cvxopt required for passivity module")
 
     if not sys.isctime() and rho is not None and nu is not None:
-        raise Error("Passivity indices for discrete time systems not supported yet.")
+        raise Exception(
+            "Passivity indices for discrete time systems not supported yet.")
+
+    if sys.ninputs != sys.noutputs:
+        raise Exception(
+            "The number of system inputs must be the same as the number of system outputs.")
 
-    (A,B,C,D) = __parse_lti__(sys)
+    (A, B, C, D) = __parse_lti__(sys)
 
     def make_LMI_matrix(P):
         if sys.isctime():
@@ -173,7 +183,7 @@ def make_LMI_matrix(P):
             )
         else:
             return 2*np.vstack((
-                np.hstack((A.T @ P @A - P , A.T @ P@B)),
+                np.hstack((A.T @ P @ A - P, A.T @ P@B)),
                 np.hstack(((A.T @ P@B).T, B.T@P@B)))
             )
 
@@ -186,13 +196,13 @@ def make_LMI_matrix(P):
 #     Constructs a linear matrix inequality and a feasibility optimization
 #     such that if a solution exists, the system is dissapative.
 
-#     The sources for the algorithm are: 
-    
+#     The sources for the algorithm are:
+
 #     McCourt, Michael J., and Panos J. Antsaklis
-#         "Demonstrating passivity and dissipativity using computational methods." 
+#         "Demonstrating passivity and dissipativity using computational methods."
 
 #     Nicholas Kottenstette and Panos J. Antsaklis
-#         "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
+#         "Relationships Between Positive Real, Passive Dissipative, & Positive Systems"
 #         equation 36.
 
 #     Parameters
@@ -202,7 +212,7 @@ def make_LMI_matrix(P):
 
 #     Returns
 #     -------
-#     bool: 
+#     bool:
 #         The input system dissapative.
 #     '''
 #     if cvx is None:
@@ -212,7 +222,7 @@ def make_LMI_matrix(P):
 
 #     def make_LMI_matrix(P):
 #         None
-        
+
 #     matrix_list =_make_basis_matrices(sys.nstates, make_LMI_matrix)
 
 #     coefficents = np.vstack(matrix_list).T
@@ -229,4 +239,4 @@ def make_LMI_matrix(P):
 #                           Gs=[cvx.matrix(coefficents)],
 #                           hs=[cvx.matrix(constants)])
 
-#     return (sol["x"] is not None)
\ No newline at end of file
+#     return (sol["x"] is not None)
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index 43f6d134d..36ed44aeb 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -27,6 +27,7 @@ def test_ispassive_ctime():
 
     assert(not passivity.ispassive(sys))
 
+
 def test_ispassive_dtime():
     A = numpy.array([[0, 1], [-2, -2]])
     B = numpy.array([[0], [1]])
@@ -44,6 +45,17 @@ def test_ispassive_dtime():
     assert(not passivity.ispassive(sys))
 
 
+def test_system_dimension():
+    A = numpy.array([[0, 1], [-2, -2]])
+    B = numpy.array([[0], [1]])
+    C = numpy.array([[-1, 2], [0, 1]])
+    D = numpy.array([[1.5], [1]])
+    sys = ss(A, B, C, D)
+
+    with pytest.raises(Exception):
+        passivity.ispassive(sys)
+
+
 A_d = numpy.array([[-2, 0], [0, 0]])
 A = numpy.array([[-3, 0], [0, -2]])
 B = numpy.array([[0], [1]])

From 1ef69f600d230c2a57bff97e6b7ef5e6b43bb824 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Thu, 14 Jul 2022 22:41:03 -0400
Subject: [PATCH 04/36] Add support for input/output passivity indices.

---
 control/passivity.py            | 98 +++++++++++++++++----------------
 control/tests/passivity_test.py | 12 +++-
 2 files changed, 63 insertions(+), 47 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 974fdb8d3..1e155b244 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -74,66 +74,72 @@ def __solve_LMI_problem__(A, C, D, nu, rho, n, make_LMI_matrix):
     # LMI for passivity from A,B,C,D
     sys_matrix_list = __make_P_basis_matrices__(n, make_LMI_matrix)
 
-    sys_coefficents = np.vstack(sys_matrix_list).T
 
     sys_constants = -np.vstack((
         np.hstack((np.zeros_like(A),  - C.T)),
         np.hstack((- C, -D - D.T)))
     )
 
-    # if nu is not None:
-    #     sys_constants+= -np.vstack((
-    #         np.hstack((np.zeros_like(A),  np.zeros_like(C.T))),
-    #         np.hstack((np.zeros_like(C),  nu*np.eye(n) )))
-    #     )
-
-    # if rho is not None:
-    #     sys_constants+= -np.vstack((
-    #         np.hstack((rho*C.T@C,  rho*C.T@D )),
-    #         np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
-    #     )
-
-    # if rho is not None and nu is not None:
-    #     sys_constants+= -np.vstack((
-    #         np.hstack((np.zeros_like(A),  -0.5*nu*rho*C.T)),
-    #         np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
-    #     )
+    if nu is not None:
+        m = D.shape[1]
+        sys_constants+= -np.vstack((
+            np.hstack((np.zeros_like(A),  np.zeros_like(C.T))),
+            np.hstack((np.zeros_like(C),  nu*np.eye(m) )))
+        )
+
+    if rho is not None:
+        sys_constants+= -np.vstack((
+            np.hstack((rho*C.T@C,  rho*C.T@D )),
+            np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
+        )
+
+    if rho is not None and nu is not None:
+        sys_constants+= -np.vstack((
+            np.hstack((np.zeros_like(A),  -0.5*nu*rho*C.T)),
+            np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
+        )
 
     # 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))
-
-    # #LMI for passivity indices
-    # rho_coefficents = []
-    # nu_coefficents = []
-    # if nu is not None and rho is None:
-    #     #pick out coefficents for rho
-    #     rho_coefficents_matrix = [np.vstack((
-    #             np.hstack((C.T@C, 0.5*nu*C.T + C.T@D)),
-    #             np.hstack(( (0.5*nu*C.T + C.T@D).T, D.T@D-nu*(D+D.T))))
-    #         ).flatten()]
-    #     rho_coefficents =np.vstack(rho_coefficents_matrix).T
-    # elif rho is not None and nu is None:
-    #     #pick out coefficents for nu
-    #     nu_coefficents_matrix = [np.vstack((
-    #             np.hstack((np.zeros_like(A), 0.5*rho*C.T)),
-    #             np.hstack(( (0.5*rho*C.T + rho*C.T@D).T, rho*D.T@D-nu*(D+D.T))))
-    #         ).flatten()]
-
+    
     number_of_opt_vars = int(
         (n**2-n)/2 + n)
     c = cvx.matrix(0.0, (number_of_opt_vars, 1))
 
+    #LMI for passivity indices
+    if nu is not None and rho is None:
+        #pick out coefficents for rho
+        rho_coefficents_matrix = np.vstack((
+                np.hstack((C.T@C, 0.5*nu*C.T + C.T@D)),
+                np.hstack(( (0.5*nu*C.T + C.T@D).T, D.T@D-nu*(D+D.T))))
+            )
+        sys_matrix_list.append(rho_coefficents_matrix.flatten())
+        c = cvx.matrix(np.append(np.array(c),-1.0))
+        P_matrix_list.append(np.zeros_like(A).flatten())
+    elif rho is not None and nu is None:
+        #pick out coefficents for nu
+        nu_coefficents_matrix = np.vstack((
+                np.hstack((np.zeros_like(A), 0.5*rho*C.T)),
+                np.hstack(( (0.5*rho*C.T + rho*C.T@D).T, rho*D.T@D)))
+            )
+        sys_matrix_list.append(nu_coefficents_matrix.flatten())
+        c = cvx.matrix(np.append(np.array(c),-1.0))
+        P_matrix_list.append(np.zeros_like(A).flatten())
+
+    sys_coefficents = np.vstack(sys_matrix_list).T
+    P_coefficents = np.vstack(P_matrix_list).T
+    P_constants = np.zeros((n, n))
+
+    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
-    sol = cvx.solvers.sdp(c,
-                          Gs=[cvx.matrix(sys_coefficents)] +
-                          [cvx.matrix(P_coefficents)],
-                          hs=[cvx.matrix(sys_constants)]+[cvx.matrix(P_constants)])
-
-    return (sol["x"] is not None)
+    sol = cvx.solvers.sdp(c, Gs = Gs, hs= hs)
+    if nu is None and rho is None:
+        return sol["x"] is not None
+    else:
+        return np.ravel(sol["x"])[-1]
 
 
 def ispassive(sys, nu=None, rho=None):
@@ -187,7 +193,7 @@ def make_LMI_matrix(P):
                 np.hstack(((A.T @ P@B).T, B.T@P@B)))
             )
 
-    return __solve_LMI_problem__(A, C, D, None, None, sys.nstates, make_LMI_matrix)
+    return __solve_LMI_problem__(A, C, D, nu, rho, sys.nstates, make_LMI_matrix)
 
 # def isQSRdissapative(sys):
 #     '''
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index 36ed44aeb..e785c1ae2 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -2,6 +2,7 @@
 Author: Mark Yeatman  
 Date: May 30, 2022
 '''
+from tokenize import Double
 import pytest
 import numpy
 from control import ss, passivity, tf, sample_system
@@ -11,7 +12,7 @@
 pytestmark = cvxoptonly
 
 
-def test_ispassive_ctime():
+def test_ispassive_ctime(capsys):
     A = numpy.array([[0, 1], [-2, -2]])
     B = numpy.array([[0], [1]])
     C = numpy.array([[-1, 2]])
@@ -21,6 +22,15 @@ def test_ispassive_ctime():
     # happy path is passive
     assert(passivity.ispassive(sys))
 
+    # happy path is passive
+    ret = passivity.ispassive(sys, rho = 0.00001)
+    assert(isinstance(ret, float))
+    assert(ret>0.0)
+
+    ret = passivity.ispassive(sys, nu = 0.00001)
+    assert(isinstance(ret, float))
+    assert(ret>0.0)
+
     # happy path not passive
     D = -D
     sys = ss(A, B, C, D)

From 30f859ab8e6ec2aec2344b39b4f5bbbade8c0611 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Thu, 14 Jul 2022 22:46:54 -0400
Subject: [PATCH 05/36] Clean up.

---
 control/passivity.py | 161 ++++++++++++++-----------------------------
 1 file changed, 51 insertions(+), 110 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 1e155b244..fec265260 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -66,15 +66,62 @@ def __P_pos_def_constraint__(n):
                 matrix_list.append(P.flatten())
     return matrix_list
 
-
-def __solve_LMI_problem__(A, C, D, nu, rho, n, make_LMI_matrix):
+def ispassive(sys, nu=None, rho=None):
     '''
-    Construct LMI for constants in passivitiy problem, format matrices and input into SDP solver. 
+    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 sources for the algorithm are: 
+
+    McCourt, Michael J., and Panos J. Antsaklis
+        "Demonstrating passivity and dissipativity using computational methods." 
+
+    Nicholas Kottenstette and Panos J. Antsaklis
+        "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
+        equation 36.
+
+    Parameters
+    ----------
+    sys: An LTI system
+        System to be checked.
+
+    Returns
+    -------
+    bool: 
+        The input system passive.
     '''
+    if cvx is None:
+        raise ModuleNotFoundError("cvxopt required for passivity module")
+
+    if not sys.isctime() and rho is not None and nu is not None:
+        raise Exception(
+            "Passivity indices for discrete time systems not supported yet.")
+
+    if sys.ninputs != sys.noutputs:
+        raise Exception(
+            "The number of system inputs must be the same as the number of system outputs.")
+
+    (A, B, C, D) = __parse_lti__(sys)
+
+    def make_LMI_matrix(P):
+        if sys.isctime():
+            return np.vstack((
+                np.hstack((A.T @ P + P@A, P@B)),
+                np.hstack((B.T@P, np.zeros_like(D))))
+            )
+        else:
+            return 2*np.vstack((
+                np.hstack((A.T @ P @ A - P, A.T @ P@B)),
+                np.hstack(((A.T @ P@B).T, B.T@P@B)))
+            )
+
+    n = sys.nstates
+
     # LMI for passivity from A,B,C,D
     sys_matrix_list = __make_P_basis_matrices__(n, make_LMI_matrix)
 
-
     sys_constants = -np.vstack((
         np.hstack((np.zeros_like(A),  - C.T)),
         np.hstack((- C, -D - D.T)))
@@ -140,109 +187,3 @@ def __solve_LMI_problem__(A, C, D, nu, rho, n, make_LMI_matrix):
         return sol["x"] is not None
     else:
         return np.ravel(sol["x"])[-1]
-
-
-def ispassive(sys, nu=None, rho=None):
-    '''
-    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 sources for the algorithm are: 
-
-    McCourt, Michael J., and Panos J. Antsaklis
-        "Demonstrating passivity and dissipativity using computational methods." 
-
-    Nicholas Kottenstette and Panos J. Antsaklis
-        "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
-        equation 36.
-
-    Parameters
-    ----------
-    sys: An LTI system
-        System to be checked.
-
-    Returns
-    -------
-    bool: 
-        The input system passive.
-    '''
-    if cvx is None:
-        raise ModuleNotFoundError("cvxopt required for passivity module")
-
-    if not sys.isctime() and rho is not None and nu is not None:
-        raise Exception(
-            "Passivity indices for discrete time systems not supported yet.")
-
-    if sys.ninputs != sys.noutputs:
-        raise Exception(
-            "The number of system inputs must be the same as the number of system outputs.")
-
-    (A, B, C, D) = __parse_lti__(sys)
-
-    def make_LMI_matrix(P):
-        if sys.isctime():
-            return np.vstack((
-                np.hstack((A.T @ P + P@A, P@B)),
-                np.hstack((B.T@P, np.zeros_like(D))))
-            )
-        else:
-            return 2*np.vstack((
-                np.hstack((A.T @ P @ A - P, A.T @ P@B)),
-                np.hstack(((A.T @ P@B).T, B.T@P@B)))
-            )
-
-    return __solve_LMI_problem__(A, C, D, nu, rho, sys.nstates, make_LMI_matrix)
-
-# def isQSRdissapative(sys):
-#     '''
-#     Indicates if a linear time invariant (LTI) system is dissapative
-
-#     Constructs a linear matrix inequality and a feasibility optimization
-#     such that if a solution exists, the system is dissapative.
-
-#     The sources for the algorithm are:
-
-#     McCourt, Michael J., and Panos J. Antsaklis
-#         "Demonstrating passivity and dissipativity using computational methods."
-
-#     Nicholas Kottenstette and Panos J. Antsaklis
-#         "Relationships Between Positive Real, Passive Dissipative, & Positive Systems"
-#         equation 36.
-
-#     Parameters
-#     ----------
-#     sys: An LTI system
-#         System to be checked.
-
-#     Returns
-#     -------
-#     bool:
-#         The input system dissapative.
-#     '''
-#     if cvx is None:
-#         raise ModuleNotFoundError("cvxopt required for passivity module")
-
-#     (A,B,C,D) = _parse_lti(sys)
-
-#     def make_LMI_matrix(P):
-#         None
-
-#     matrix_list =_make_basis_matrices(sys.nstates, make_LMI_matrix)
-
-#     coefficents = np.vstack(matrix_list).T
-
-#     constants = None
-
-#     number_of_opt_vars = int(
-#         (sys.nstates**2-sys.nstates)/2 + sys.nstates)
-#     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)

From caa2614b268059922781ce801e87b0e4d589f8ca Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Thu, 14 Jul 2022 22:48:58 -0400
Subject: [PATCH 06/36] Run autopep8 formatter.

---
 control/passivity.py            | 41 +++++++++++++++++----------------
 control/tests/passivity_test.py |  8 +++----
 2 files changed, 25 insertions(+), 24 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index fec265260..b24090f75 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -66,6 +66,7 @@ def __P_pos_def_constraint__(n):
                 matrix_list.append(P.flatten())
     return matrix_list
 
+
 def ispassive(sys, nu=None, rho=None):
     '''
     Indicates if a linear time invariant (LTI) system is passive
@@ -129,48 +130,48 @@ def make_LMI_matrix(P):
 
     if nu is not None:
         m = D.shape[1]
-        sys_constants+= -np.vstack((
+        sys_constants += -np.vstack((
             np.hstack((np.zeros_like(A),  np.zeros_like(C.T))),
-            np.hstack((np.zeros_like(C),  nu*np.eye(m) )))
+            np.hstack((np.zeros_like(C),  nu*np.eye(m))))
         )
 
     if rho is not None:
-        sys_constants+= -np.vstack((
-            np.hstack((rho*C.T@C,  rho*C.T@D )),
-            np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
+        sys_constants += -np.vstack((
+            np.hstack((rho*C.T@C,  rho*C.T@D)),
+            np.hstack(((rho*C.T@D).T, rho*D.T@D)))
         )
 
     if rho is not None and nu is not None:
-        sys_constants+= -np.vstack((
+        sys_constants += -np.vstack((
             np.hstack((np.zeros_like(A),  -0.5*nu*rho*C.T)),
-            np.hstack(( (rho*C.T@D).T, rho*D.T@D)))
+            np.hstack(((rho*C.T@D).T, rho*D.T@D)))
         )
 
     # LMI to ensure P is positive definite
     P_matrix_list = __P_pos_def_constraint__(n)
-    
+
     number_of_opt_vars = int(
         (n**2-n)/2 + n)
     c = cvx.matrix(0.0, (number_of_opt_vars, 1))
 
-    #LMI for passivity indices
+    # LMI for passivity indices
     if nu is not None and rho is None:
-        #pick out coefficents for rho
+        # pick out coefficents for rho
         rho_coefficents_matrix = np.vstack((
-                np.hstack((C.T@C, 0.5*nu*C.T + C.T@D)),
-                np.hstack(( (0.5*nu*C.T + C.T@D).T, D.T@D-nu*(D+D.T))))
-            )
+            np.hstack((C.T@C, 0.5*nu*C.T + C.T@D)),
+            np.hstack(((0.5*nu*C.T + C.T@D).T, D.T@D-nu*(D+D.T))))
+        )
         sys_matrix_list.append(rho_coefficents_matrix.flatten())
-        c = cvx.matrix(np.append(np.array(c),-1.0))
+        c = cvx.matrix(np.append(np.array(c), -1.0))
         P_matrix_list.append(np.zeros_like(A).flatten())
     elif rho is not None and nu is None:
-        #pick out coefficents for nu
+        # pick out coefficents for nu
         nu_coefficents_matrix = np.vstack((
-                np.hstack((np.zeros_like(A), 0.5*rho*C.T)),
-                np.hstack(( (0.5*rho*C.T + rho*C.T@D).T, rho*D.T@D)))
-            )
+            np.hstack((np.zeros_like(A), 0.5*rho*C.T)),
+            np.hstack(((0.5*rho*C.T + rho*C.T@D).T, rho*D.T@D)))
+        )
         sys_matrix_list.append(nu_coefficents_matrix.flatten())
-        c = cvx.matrix(np.append(np.array(c),-1.0))
+        c = cvx.matrix(np.append(np.array(c), -1.0))
         P_matrix_list.append(np.zeros_like(A).flatten())
 
     sys_coefficents = np.vstack(sys_matrix_list).T
@@ -182,7 +183,7 @@ def make_LMI_matrix(P):
 
     # crunch feasibility solution
     cvx.solvers.options['show_progress'] = False
-    sol = cvx.solvers.sdp(c, Gs = Gs, hs= hs)
+    sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
     if nu is None and rho is None:
         return sol["x"] is not None
     else:
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index e785c1ae2..b9bb6017e 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -23,13 +23,13 @@ def test_ispassive_ctime(capsys):
     assert(passivity.ispassive(sys))
 
     # happy path is passive
-    ret = passivity.ispassive(sys, rho = 0.00001)
+    ret = passivity.ispassive(sys, rho=0.00001)
     assert(isinstance(ret, float))
-    assert(ret>0.0)
+    assert(ret > 0.0)
 
-    ret = passivity.ispassive(sys, nu = 0.00001)
+    ret = passivity.ispassive(sys, nu=0.00001)
     assert(isinstance(ret, float))
-    assert(ret>0.0)
+    assert(ret > 0.0)
 
     # happy path not passive
     D = -D

From b249ac49a54b61de757c703996e81f031e08dacd Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Thu, 14 Jul 2022 23:09:25 -0400
Subject: [PATCH 07/36] Add getPassiveIndex function.

---
 control/passivity.py | 30 ++++++++++++++++++++++++++++--
 1 file changed, 28 insertions(+), 2 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index b24090f75..209463d8e 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -87,11 +87,15 @@ def ispassive(sys, nu=None, rho=None):
     ----------
     sys: An LTI system
         System to be checked.
+    nu: float
+        Concrete value for input passivity index. 
+    rho: float
+        Concrete value for output passivity index. 
 
     Returns
     -------
-    bool: 
-        The input system passive.
+    bool or float: 
+        The input system passive, or the passivity index "opposite" the input. 
     '''
     if cvx is None:
         raise ModuleNotFoundError("cvxopt required for passivity module")
@@ -188,3 +192,25 @@ def make_LMI_matrix(P):
         return sol["x"] is not None
     else:
         return np.ravel(sol["x"])[-1]
+
+def getPassiveIndex(sys, index_type = None):
+    '''
+    Returns the passivity index associated with the input string. 
+    Parameters
+    ----------
+    sys: An LTI system
+        System to be checked.
+    index_type: String
+        Must be 'input' or 'output'. Indicates which passivity index will be returned. 
+
+    Returns
+    -------
+    float: 
+        The passivity index 
+    '''
+    if index_type is None:
+        raise Exception("Must provide index_type of 'input' or 'output'.")
+    if index_type == "input":
+        return ispassive(sys, rho = 0.000001)
+    if index_type == "output":
+        return ispassive(sys, nu = 0.000001)
\ No newline at end of file

From 439134f1b95563e477ca580cddc282bc2a8ad302 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Thu, 14 Jul 2022 23:33:14 -0400
Subject: [PATCH 08/36] Test commit.

---
 control/passivity.py            |  8 ++++++--
 control/tests/passivity_test.py | 21 +++++++++++----------
 control/tests/rlocus_test.py    |  2 +-
 3 files changed, 18 insertions(+), 13 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 209463d8e..a299e44a0 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -211,6 +211,10 @@ def getPassiveIndex(sys, index_type = None):
     if index_type is None:
         raise Exception("Must provide index_type of 'input' or 'output'.")
     if index_type == "input":
-        return ispassive(sys, rho = 0.000001)
+        nu = ispassive(sys, rho = 0.000001)
+        print(nu)
+        return nu
     if index_type == "output":
-        return ispassive(sys, nu = 0.000001)
\ No newline at end of file
+        nu = ispassive(sys, nu = 0.000001)
+        print(nu)
+        return nu
\ No newline at end of file
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index b9bb6017e..a3b9bb5d3 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -11,8 +11,18 @@
 
 pytestmark = cvxoptonly
 
+def test_passivity_indices():
+    sys = tf([1,1,5,0.1],[1,2,3,4])
 
-def test_ispassive_ctime(capsys):
+    nu = passivity.getPassiveIndex(sys, 'input')
+    assert(isinstance(nu, float))
+    assert(nu > 0.0250)
+
+    rho = passivity.getPassiveIndex(sys, 'output')
+    assert(isinstance(rho, float))
+    assert(rho < 0.2583)
+
+def test_ispassive_ctime():
     A = numpy.array([[0, 1], [-2, -2]])
     B = numpy.array([[0], [1]])
     C = numpy.array([[-1, 2]])
@@ -22,15 +32,6 @@ def test_ispassive_ctime(capsys):
     # happy path is passive
     assert(passivity.ispassive(sys))
 
-    # happy path is passive
-    ret = passivity.ispassive(sys, rho=0.00001)
-    assert(isinstance(ret, float))
-    assert(ret > 0.0)
-
-    ret = passivity.ispassive(sys, nu=0.00001)
-    assert(isinstance(ret, float))
-    assert(ret > 0.0)
-
     # happy path not passive
     D = -D
     sys = ss(A, B, C, D)
diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py
index a0ecebb15..01dfa8a92 100644
--- a/control/tests/rlocus_test.py
+++ b/control/tests/rlocus_test.py
@@ -112,7 +112,7 @@ def test_root_locus_zoom(self):
         assert_array_almost_equal(zoom_x, zoom_x_valid)
         assert_array_almost_equal(zoom_y, zoom_y_valid)
 
-    @pytest.mark.timeout(2)
+    # @pytest.mark.timeout(2)
     def test_rlocus_default_wn(self):
         """Check that default wn calculation works properly"""
         #

From b7adf478d9101322a3f32fedf35e1a56717846ba Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Sat, 16 Jul 2022 21:44:52 -0400
Subject: [PATCH 09/36] Rework passivity indices so that unit tests pass.

---
 control/passivity.py            | 223 +++++++++++++++-----------------
 control/tests/passivity_test.py |  30 +++--
 2 files changed, 119 insertions(+), 134 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index a299e44a0..51f368365 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -6,35 +6,17 @@
 from msilib.schema import Error
 import numpy as np
 from control import statesp as ss
+from control.modelsimp import minreal
 
 try:
     import cvxopt as cvx
 except ImportError as e:
     cvx = None
 
+eps = np.nextafter(0, 1)
 
-def __parse_lti__(sys):
-    '''
-    Utility function to parse LTI input for passivity module functions
-    '''
-    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)
-
-    return (A, B, C, D)
 
-
-def __make_P_basis_matrices__(n, make_LMI_matrix_func):
+def __make_P_basis_matrices__(n, rho, nu, make_LMI_matrix_func):
     '''
     Utility function to make basis matrices for a LMI from a 
     functional make_LMI_matrix_func and a symmetric matrix P of size n by n
@@ -47,7 +29,9 @@ def __make_P_basis_matrices__(n, make_LMI_matrix_func):
                 P = np.zeros((n, n))
                 P[i, j] = 1.0
                 P[j, i] = 1.0
-                matrix_list.append(make_LMI_matrix_func(P).flatten())
+                matrix_list.append(make_LMI_matrix_func(P, 0, 0, 0).flatten())
+    P = eps*np.eye(n)
+    matrix_list.append(make_LMI_matrix_func(P, rho, nu, 0).flatten())
     return matrix_list
 
 
@@ -64,136 +48,135 @@ def __P_pos_def_constraint__(n):
                 P[i, j] = -1.0
                 P[j, i] = -1.0
                 matrix_list.append(P.flatten())
+    matrix_list.append(np.zeros((n, n)).flatten())
     return matrix_list
 
 
-def ispassive(sys, nu=None, rho=None):
-    '''
-    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 sources for the algorithm are: 
-
-    McCourt, Michael J., and Panos J. Antsaklis
-        "Demonstrating passivity and dissipativity using computational methods." 
-
-    Nicholas Kottenstette and Panos J. Antsaklis
-        "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
-        equation 36.
-
-    Parameters
-    ----------
-    sys: An LTI system
-        System to be checked.
-    nu: float
-        Concrete value for input passivity index. 
-    rho: float
-        Concrete value for output passivity index. 
-
-    Returns
-    -------
-    bool or float: 
-        The input system passive, or the passivity index "opposite" the input. 
-    '''
+def __ispassive__(sys, rho=None, nu=None):
     if cvx is None:
         raise ModuleNotFoundError("cvxopt required for passivity module")
 
-    if not sys.isctime() and rho is not None and nu is not None:
+    if sys.isdtime():
         raise Exception(
-            "Passivity indices for discrete time systems not supported yet.")
+            "Passivity for discrete time systems not supported yet.")
 
     if sys.ninputs != sys.noutputs:
         raise Exception(
             "The number of system inputs must be the same as the number of system outputs.")
 
-    (A, B, C, D) = __parse_lti__(sys)
+    if rho is None and nu is None:
+        raise Exception("rho or nu must be given a float value.")
 
-    def make_LMI_matrix(P):
-        if sys.isctime():
-            return np.vstack((
-                np.hstack((A.T @ P + P@A, P@B)),
-                np.hstack((B.T@P, np.zeros_like(D))))
-            )
-        else:
-            return 2*np.vstack((
-                np.hstack((A.T @ P @ A - P, A.T @ P@B)),
-                np.hstack(((A.T @ P@B).T, B.T@P@B)))
-            )
+    sys = ss._convert_to_statespace(sys)
 
-    n = sys.nstates
+    A = sys.A
+    B = sys.B
+    C = sys.C
+    D = sys.D
 
-    # LMI for passivity from A,B,C,D
-    sys_matrix_list = __make_P_basis_matrices__(n, make_LMI_matrix)
+    # account for strictly proper systems
+    [n, m] = D.shape
+    D = D + eps * np.eye(n, m)
 
+    [n, _] = A.shape
+    A = A - eps*np.eye(n)
+
+    def make_LMI_matrix(P, rho, nu, one):
+        off_diag = P@B - 1.0/2.0*(one+rho*nu)*C.T + rho*C.T*D
+        return np.vstack((
+            np.hstack((A.T @ P + P@A + rho*C.T@C,  off_diag)),
+            np.hstack((off_diag.T, rho*D.T@D -
+                       (one+rho*nu)*(D+D.T)+nu*np.eye(m)))
+        ))
+
+    n = sys.nstates
+
+    # LMI for passivity indices from A,B,C,D
+    sys_matrix_list = list()
     sys_constants = -np.vstack((
-        np.hstack((np.zeros_like(A),  - C.T)),
-        np.hstack((- C, -D - D.T)))
+        np.hstack((np.zeros_like(A),  np.zeros_like(C.T))),
+        np.hstack((np.zeros_like(C), np.zeros_like(D))))
     )
 
-    if nu is not None:
-        m = D.shape[1]
-        sys_constants += -np.vstack((
-            np.hstack((np.zeros_like(A),  np.zeros_like(C.T))),
-            np.hstack((np.zeros_like(C),  nu*np.eye(m))))
-        )
-
     if rho is not None:
-        sys_constants += -np.vstack((
-            np.hstack((rho*C.T@C,  rho*C.T@D)),
-            np.hstack(((rho*C.T@D).T, rho*D.T@D)))
-        )
+        sys_matrix_list = __make_P_basis_matrices__(
+            n, eps, 1.0, make_LMI_matrix)
+        sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, eps, 1.0)
+    else:
+        sys_matrix_list = __make_P_basis_matrices__(
+            n, 1.0, eps, make_LMI_matrix)
+        sys_constants = -make_LMI_matrix(np.zeros_like(A), eps, nu, 1.0)
 
-    if rho is not None and nu is not None:
-        sys_constants += -np.vstack((
-            np.hstack((np.zeros_like(A),  -0.5*nu*rho*C.T)),
-            np.hstack(((rho*C.T@D).T, rho*D.T@D)))
-        )
+    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))
-
-    # LMI for passivity indices
-    if nu is not None and rho is None:
-        # pick out coefficents for rho
-        rho_coefficents_matrix = np.vstack((
-            np.hstack((C.T@C, 0.5*nu*C.T + C.T@D)),
-            np.hstack(((0.5*nu*C.T + C.T@D).T, D.T@D-nu*(D+D.T))))
-        )
-        sys_matrix_list.append(rho_coefficents_matrix.flatten())
-        c = cvx.matrix(np.append(np.array(c), -1.0))
-        P_matrix_list.append(np.zeros_like(A).flatten())
-    elif rho is not None and nu is None:
-        # pick out coefficents for nu
-        nu_coefficents_matrix = np.vstack((
-            np.hstack((np.zeros_like(A), 0.5*rho*C.T)),
-            np.hstack(((0.5*rho*C.T + rho*C.T@D).T, rho*D.T@D)))
-        )
-        sys_matrix_list.append(nu_coefficents_matrix.flatten())
-        c = cvx.matrix(np.append(np.array(c), -1.0))
-        P_matrix_list.append(np.zeros_like(A).flatten())
-
-    sys_coefficents = np.vstack(sys_matrix_list).T
-    P_coefficents = np.vstack(P_matrix_list).T
-    P_constants = np.zeros((n, n))
+    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)]
+    hs = [cvx.matrix(sys_constants)] + [cvx.matrix(P_constants)]
 
     # crunch feasibility solution
     cvx.solvers.options['show_progress'] = False
     sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
-    if nu is None and rho is None:
-        return sol["x"] is not None
+    return sol["x"]
+
+
+def ispassive(sys, rho=None, nu=None):
+    '''
+    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 sources for the algorithm are: 
+
+    McCourt, Michael J., and Panos J. Antsaklis
+        "Demonstrating passivity and dissipativity using computational methods." 
+
+    Nicholas Kottenstette and Panos J. Antsaklis
+        "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
+        equation 36.
+
+    Parameters
+    ----------
+    sys: An LTI system
+        System to be checked.
+    nu: float
+        Concrete value for input passivity index. 
+    rho: float
+        Concrete value for output passivity index. 
+
+    Returns
+    -------
+    bool or float: 
+        The input system passive, or the passivity index "opposite" the input. 
+    '''
+    if rho is not None and nu is not None:
+        return __ispassive__(sys, rho, nu) is not None
+    elif rho is None and nu is not None:
+        rho = __ispassive__(sys, nu=nu)[-1]
+        print(rho)
+        return rho
+    elif nu is None and rho is not None:
+        nu = __ispassive__(sys, rho=rho)[-1]
+        print(nu)
+        return nu
     else:
-        return np.ravel(sol["x"])[-1]
+        rho = __ispassive__(sys, nu=eps)[-1]
+        nu = __ispassive__(sys, rho=eps)[-1]
+        print((rho, nu))
+        return rho > 0 or nu > 0
 
-def getPassiveIndex(sys, index_type = None):
+
+def getPassiveIndex(sys, index_type=None):
     '''
     Returns the passivity index associated with the input string. 
     Parameters
@@ -211,10 +194,8 @@ def getPassiveIndex(sys, index_type = None):
     if index_type is None:
         raise Exception("Must provide index_type of 'input' or 'output'.")
     if index_type == "input":
-        nu = ispassive(sys, rho = 0.000001)
-        print(nu)
+        nu = ispassive(sys, rho=eps)
         return nu
     if index_type == "output":
-        nu = ispassive(sys, nu = 0.000001)
-        print(nu)
-        return nu
\ No newline at end of file
+        rho = ispassive(sys, nu=eps)
+        return rho
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index a3b9bb5d3..fefe14921 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -11,8 +11,9 @@
 
 pytestmark = cvxoptonly
 
+
 def test_passivity_indices():
-    sys = tf([1,1,5,0.1],[1,2,3,4])
+    sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
 
     nu = passivity.getPassiveIndex(sys, 'input')
     assert(isinstance(nu, float))
@@ -22,6 +23,7 @@ def test_passivity_indices():
     assert(isinstance(rho, float))
     assert(rho < 0.2583)
 
+
 def test_ispassive_ctime():
     A = numpy.array([[0, 1], [-2, -2]])
     B = numpy.array([[0], [1]])
@@ -46,14 +48,8 @@ def test_ispassive_dtime():
     D = numpy.array([[1.5]])
     sys = ss(A, B, C, D)
     sys = sample_system(sys, 1, method='bilinear')
-    # 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))
+    with pytest.raises(Exception):
+        passivity.ispassive(sys)
 
 
 def test_system_dimension():
@@ -80,11 +76,8 @@ def test_system_dimension():
      ((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)])
+     ((A*0, B, C, D), True)])
 def test_ispassive_edge_cases(test_input, expected):
-
-    # strictly proper
     A = test_input[0]
     B = test_input[1]
     C = test_input[2]
@@ -93,6 +86,17 @@ def test_ispassive_edge_cases(test_input, expected):
     assert(passivity.ispassive(sys) == expected)
 
 
+def test_ispassive_all_zeros():
+    A = numpy.array([[0]])
+    B = numpy.array([[0]])
+    C = numpy.array([[0]])
+    D = numpy.array([[0]])
+    sys = ss(A, B, C, D)
+
+    with pytest.raises(Exception):
+        passivity.ispassive(sys)
+
+
 def test_transfer_function():
     sys = tf([1], [1, 2])
     assert(passivity.ispassive(sys))

From c5dc526651f6d4026179471f8a596ee137c52f0f Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Sat, 16 Jul 2022 21:47:08 -0400
Subject: [PATCH 10/36] Remove development print statement.

---
 control/passivity.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 51f368365..2c6df96a6 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -163,16 +163,13 @@ def ispassive(sys, rho=None, nu=None):
         return __ispassive__(sys, rho, nu) is not None
     elif rho is None and nu is not None:
         rho = __ispassive__(sys, nu=nu)[-1]
-        print(rho)
         return rho
     elif nu is None and rho is not None:
         nu = __ispassive__(sys, rho=rho)[-1]
-        print(nu)
         return nu
     else:
         rho = __ispassive__(sys, nu=eps)[-1]
         nu = __ispassive__(sys, rho=eps)[-1]
-        print((rho, nu))
         return rho > 0 or nu > 0
 
 

From 6db437d8af0c3907667d6f025bce29ac31e266f5 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Sat, 16 Jul 2022 23:15:49 -0400
Subject: [PATCH 11/36] Finally getting correct values for passivity indices,
 source paper had a typo.

---
 control/passivity.py | 25 +++++++++++++++++--------
 1 file changed, 17 insertions(+), 8 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 2c6df96a6..ba75e04fd 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -86,7 +86,7 @@ def make_LMI_matrix(P, rho, nu, one):
         return np.vstack((
             np.hstack((A.T @ P + P@A + rho*C.T@C,  off_diag)),
             np.hstack((off_diag.T, rho*D.T@D -
-                       (one+rho*nu)*(D+D.T)+nu*np.eye(m)))
+                       1.0/2.0*(one+rho*nu)*(D+D.T)+nu*np.eye(m)))
         ))
 
     n = sys.nstates
@@ -162,15 +162,24 @@ def ispassive(sys, rho=None, nu=None):
     if rho is not None and nu is not None:
         return __ispassive__(sys, rho, nu) is not None
     elif rho is None and nu is not None:
-        rho = __ispassive__(sys, nu=nu)[-1]
-        return rho
+        sol = __ispassive__(sys, nu=nu)
+        if sol is not None:
+            return sol[-1]
+        else:
+            return None
     elif nu is None and rho is not None:
-        nu = __ispassive__(sys, rho=rho)[-1]
-        return nu
+        sol = __ispassive__(sys, rho=rho)
+        if sol is not None:
+            return sol[-1]
+        else:
+            return None
     else:
-        rho = __ispassive__(sys, nu=eps)[-1]
-        nu = __ispassive__(sys, rho=eps)[-1]
-        return rho > 0 or nu > 0
+        sol_rho = __ispassive__(sys, nu=eps)
+        sol_nu = __ispassive__(sys, rho=eps)
+        if sol_rho is not None and sol_nu is not None:
+            return sol_rho[-1] > 0 or sol_nu[-1] > 0
+        else:
+            return False
 
 
 def getPassiveIndex(sys, index_type=None):

From b331e750df0accf988780e178b59e74ea712c6f8 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Sun, 17 Jul 2022 23:05:52 -0400
Subject: [PATCH 12/36] Add passivity indices for discrete system and unit
 tests.

---
 control/passivity.py            | 44 +++++++++++++++++++-----
 control/tests/passivity_test.py | 59 ++++++++++++++++++++++++---------
 2 files changed, 79 insertions(+), 24 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index ba75e04fd..c3a44e8b5 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -56,9 +56,9 @@ def __ispassive__(sys, rho=None, nu=None):
     if cvx is None:
         raise ModuleNotFoundError("cvxopt required for passivity module")
 
-    if sys.isdtime():
-        raise Exception(
-            "Passivity for discrete time systems not supported yet.")
+    # if sys.isdtime():
+    #     raise Exception(
+    #         "Passivity for discrete time systems not supported yet.")
 
     if sys.ninputs != sys.noutputs:
         raise Exception(
@@ -82,12 +82,38 @@ def __ispassive__(sys, rho=None, nu=None):
     A = A - eps*np.eye(n)
 
     def make_LMI_matrix(P, rho, nu, one):
-        off_diag = P@B - 1.0/2.0*(one+rho*nu)*C.T + rho*C.T*D
-        return np.vstack((
-            np.hstack((A.T @ P + P@A + rho*C.T@C,  off_diag)),
-            np.hstack((off_diag.T, rho*D.T@D -
-                       1.0/2.0*(one+rho*nu)*(D+D.T)+nu*np.eye(m)))
-        ))
+        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, -(D.T@Q@D + D.T@S + S.T@D + R)))
+            ))
+
+    # def make_LMI_matrix(P, rho, nu, one):
+    #     if sys.isctime():
+    #         off_diag = P@B - 1.0/2.0*(one+rho*nu)*C.T + rho*C.T*D
+    #         return np.vstack((
+    #             np.hstack((A.T @ P + P@A + rho*C.T@C,  off_diag)),
+    #             np.hstack((off_diag.T, rho*D.T@D -
+    #                     1.0/2.0*(one+rho*nu)*(D+D.T)+nu*np.eye(m)))
+    #         ))
+    #     else:
+    #         off_diag = A.T@P@B - 1.0/2.0*(one+rho*nu)*C.T + rho*C.T*D
+    #         return np.vstack((
+    #             np.hstack((A.T @ P  @ A - P + rho*C.T@C,  off_diag)),
+    #             np.hstack((off_diag.T, rho*D.T@D -
+    #                     1.0/2.0*(one+rho*nu)*(D+D.T)+nu*np.eye(m)))
+    #         ))
 
     n = sys.nstates
 
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index fefe14921..b82184316 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -5,25 +5,13 @@
 from tokenize import Double
 import pytest
 import numpy
-from control import ss, passivity, tf, sample_system
+from control import ss, passivity, tf, sample_system, parallel, feedback
 from control.tests.conftest import cvxoptonly
 
 
 pytestmark = cvxoptonly
 
 
-def test_passivity_indices():
-    sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
-
-    nu = passivity.getPassiveIndex(sys, 'input')
-    assert(isinstance(nu, float))
-    assert(nu > 0.0250)
-
-    rho = passivity.getPassiveIndex(sys, 'output')
-    assert(isinstance(rho, float))
-    assert(rho < 0.2583)
-
-
 def test_ispassive_ctime():
     A = numpy.array([[0, 1], [-2, -2]])
     B = numpy.array([[0], [1]])
@@ -48,8 +36,49 @@ def test_ispassive_dtime():
     D = numpy.array([[1.5]])
     sys = ss(A, B, C, D)
     sys = sample_system(sys, 1, method='bilinear')
-    with pytest.raises(Exception):
-        passivity.ispassive(sys)
+    assert(passivity.ispassive(sys))
+
+
+def test_passivity_indices_ctime():
+    sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
+
+    nu = passivity.getPassiveIndex(sys, 'input')
+    rho = passivity.getPassiveIndex(sys, 'output')
+
+    assert(isinstance(nu, float))
+
+    sys_ff_nu = parallel(-nu, sys)
+    sys_fb_rho = feedback(rho, sys, sign=1)
+
+    assert(sys_ff_nu.ispassive())
+    assert(sys_fb_rho.ispassive())
+
+    sys_ff_nu = parallel(-nu-1e-6, sys)
+    sys_fb_rho = feedback(rho+1e-6, sys, sign=1)
+
+    assert(not sys_ff_nu.ispassive())
+    assert(not sys_fb_rho.ispassive())
+
+
+def test_passivity_indices_dtime():
+    sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
+    sys = sample_system(sys, Ts=0.01, alpha=1, method="bilinear")
+    nu = passivity.getPassiveIndex(sys, 'input')
+    rho = passivity.getPassiveIndex(sys, 'output')
+
+    assert(isinstance(nu, float))
+
+    sys_ff_nu = parallel(-nu, sys)
+    sys_fb_rho = feedback(rho, sys, sign=1)
+
+    assert(sys_ff_nu.ispassive())
+    assert(sys_fb_rho.ispassive())
+
+    sys_ff_nu = parallel(-nu-1e-6, sys)
+    sys_fb_rho = feedback(rho+1e-6, sys, sign=1)
+
+    assert(not sys_ff_nu.ispassive())
+    assert(not sys_fb_rho.ispassive())
 
 
 def test_system_dimension():

From 8d7c75b2e54ad7dfa97b34cf6d4df9eab41ff702 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Sun, 17 Jul 2022 23:27:41 -0400
Subject: [PATCH 13/36] Remove commented code.

---
 control/passivity.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/control/passivity.py b/control/passivity.py
index c3a44e8b5..7cba4840b 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -1,6 +1,6 @@
 '''
 Author: Mark Yeatman  
-Date: May 15, 2022
+Date: July 17, 2022
 '''
 
 from msilib.schema import Error

From f837542c52ff6feac1978a00ec33b55dbe0f3aea Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Sun, 17 Jul 2022 23:34:25 -0400
Subject: [PATCH 14/36] Actually remove commented code.

---
 control/passivity.py | 20 --------------------
 1 file changed, 20 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 7cba4840b..abafbb332 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -56,10 +56,6 @@ def __ispassive__(sys, rho=None, nu=None):
     if cvx is None:
         raise ModuleNotFoundError("cvxopt required for passivity module")
 
-    # if sys.isdtime():
-    #     raise Exception(
-    #         "Passivity for discrete time systems not supported yet.")
-
     if sys.ninputs != sys.noutputs:
         raise Exception(
             "The number of system inputs must be the same as the number of system outputs.")
@@ -99,22 +95,6 @@ def make_LMI_matrix(P, rho, nu, one):
                 np.hstack((off_diag.T, -(D.T@Q@D + D.T@S + S.T@D + R)))
             ))
 
-    # def make_LMI_matrix(P, rho, nu, one):
-    #     if sys.isctime():
-    #         off_diag = P@B - 1.0/2.0*(one+rho*nu)*C.T + rho*C.T*D
-    #         return np.vstack((
-    #             np.hstack((A.T @ P + P@A + rho*C.T@C,  off_diag)),
-    #             np.hstack((off_diag.T, rho*D.T@D -
-    #                     1.0/2.0*(one+rho*nu)*(D+D.T)+nu*np.eye(m)))
-    #         ))
-    #     else:
-    #         off_diag = A.T@P@B - 1.0/2.0*(one+rho*nu)*C.T + rho*C.T*D
-    #         return np.vstack((
-    #             np.hstack((A.T @ P  @ A - P + rho*C.T@C,  off_diag)),
-    #             np.hstack((off_diag.T, rho*D.T@D -
-    #                     1.0/2.0*(one+rho*nu)*(D+D.T)+nu*np.eye(m)))
-    #         ))
-
     n = sys.nstates
 
     # LMI for passivity indices from A,B,C,D

From b7a79c17132419c1d771d3c05df9ffc58d7ca318 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 18 Jul 2022 10:03:00 -0400
Subject: [PATCH 15/36] Code cleanup.

---
 control/passivity.py            | 1 -
 control/tests/passivity_test.py | 1 -
 control/tests/rlocus_test.py    | 2 +-
 3 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index abafbb332..237a70046 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -3,7 +3,6 @@
 Date: July 17, 2022
 '''
 
-from msilib.schema import Error
 import numpy as np
 from control import statesp as ss
 from control.modelsimp import minreal
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index b82184316..db185a5a7 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -2,7 +2,6 @@
 Author: Mark Yeatman  
 Date: May 30, 2022
 '''
-from tokenize import Double
 import pytest
 import numpy
 from control import ss, passivity, tf, sample_system, parallel, feedback
diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py
index 01dfa8a92..a0ecebb15 100644
--- a/control/tests/rlocus_test.py
+++ b/control/tests/rlocus_test.py
@@ -112,7 +112,7 @@ def test_root_locus_zoom(self):
         assert_array_almost_equal(zoom_x, zoom_x_valid)
         assert_array_almost_equal(zoom_y, zoom_y_valid)
 
-    # @pytest.mark.timeout(2)
+    @pytest.mark.timeout(2)
     def test_rlocus_default_wn(self):
         """Check that default wn calculation works properly"""
         #

From 0bf1003b94f94b5b68858b07915441cfe00b7480 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 18 Jul 2022 10:17:03 -0400
Subject: [PATCH 16/36] Use more specific exceptions.

---
 control/passivity.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 237a70046..36866a792 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -5,7 +5,7 @@
 
 import numpy as np
 from control import statesp as ss
-from control.modelsimp import minreal
+from control.exception import ControlArgument, ControlDimension
 
 try:
     import cvxopt as cvx
@@ -56,11 +56,11 @@ def __ispassive__(sys, rho=None, nu=None):
         raise ModuleNotFoundError("cvxopt required for passivity module")
 
     if sys.ninputs != sys.noutputs:
-        raise Exception(
+        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 Exception("rho or nu must be given a float value.")
+        raise ControlArgument("rho or nu must be given a float value.")
 
     sys = ss._convert_to_statespace(sys)
 
@@ -203,7 +203,7 @@ def getPassiveIndex(sys, index_type=None):
         The passivity index 
     '''
     if index_type is None:
-        raise Exception("Must provide index_type of 'input' or 'output'.")
+        raise ControlArgument("Must provide index_type of 'input' or 'output'.")
     if index_type == "input":
         nu = ispassive(sys, rho=eps)
         return nu

From 7132945514b954bd49b55ca51211dca3c2ab6364 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 18 Jul 2022 10:51:30 -0400
Subject: [PATCH 17/36] Better doc string, use more specific exceptions in
 pytest.

---
 control/passivity.py            | 5 +++--
 control/tests/passivity_test.py | 6 +++---
 2 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 36866a792..a7a85af15 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -161,8 +161,9 @@ def ispassive(sys, rho=None, nu=None):
 
     Returns
     -------
-    bool or float: 
-        The input system passive, or the passivity index "opposite" the input. 
+    bool, float, or None: 
+        The input system is passive, or the passivity index "opposite" the input. 
+        If the problem is unfeasiable when requesting the passivity index, returns None.
     '''
     if rho is not None and nu is not None:
         return __ispassive__(sys, rho, nu) is not None
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index db185a5a7..a6b7ee8c0 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -6,7 +6,7 @@
 import numpy
 from control import ss, passivity, tf, sample_system, parallel, feedback
 from control.tests.conftest import cvxoptonly
-
+from control.exception import ControlDimension
 
 pytestmark = cvxoptonly
 
@@ -87,7 +87,7 @@ def test_system_dimension():
     D = numpy.array([[1.5], [1]])
     sys = ss(A, B, C, D)
 
-    with pytest.raises(Exception):
+    with pytest.raises(ControlDimension):
         passivity.ispassive(sys)
 
 
@@ -121,7 +121,7 @@ def test_ispassive_all_zeros():
     D = numpy.array([[0]])
     sys = ss(A, B, C, D)
 
-    with pytest.raises(Exception):
+    with pytest.raises(ValueError):
         passivity.ispassive(sys)
 
 

From 0d322a3a32c7731cbba54c4fa6d8f11570ef4b37 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 18 Jul 2022 12:55:44 -0400
Subject: [PATCH 18/36] Rename getPassiveIndex to get_passivity_index.

---
 control/passivity.py            | 2 +-
 control/tests/passivity_test.py | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index a7a85af15..d6fcbfc1b 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -188,7 +188,7 @@ def ispassive(sys, rho=None, nu=None):
             return False
 
 
-def getPassiveIndex(sys, index_type=None):
+def get_passivity_index(sys, index_type=None):
     '''
     Returns the passivity index associated with the input string. 
     Parameters
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index a6b7ee8c0..0889d9746 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -41,8 +41,8 @@ def test_ispassive_dtime():
 def test_passivity_indices_ctime():
     sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
 
-    nu = passivity.getPassiveIndex(sys, 'input')
-    rho = passivity.getPassiveIndex(sys, 'output')
+    nu = passivity.get_passivity_index(sys, 'input')
+    rho = passivity.get_passivity_index(sys, 'output')
 
     assert(isinstance(nu, float))
 
@@ -62,8 +62,8 @@ def test_passivity_indices_ctime():
 def test_passivity_indices_dtime():
     sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
     sys = sample_system(sys, Ts=0.01, alpha=1, method="bilinear")
-    nu = passivity.getPassiveIndex(sys, 'input')
-    rho = passivity.getPassiveIndex(sys, 'output')
+    nu = passivity.get_passivity_index(sys, 'input')
+    rho = passivity.get_passivity_index(sys, 'output')
 
     assert(isinstance(nu, float))
 

From 0609a171a4a82a0f00f252b9c42b95c4011ef918 Mon Sep 17 00:00:00 2001
From: Mark <myeatman@oxefit.com>
Date: Thu, 21 Jul 2022 23:22:20 -0500
Subject: [PATCH 19/36] Fix unit discrete time unit tests/ implementation.

---
 control/passivity.py            |  2 +-
 control/tests/passivity_test.py | 18 +++++++++---------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index d6fcbfc1b..bba474a5a 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -91,7 +91,7 @@ def make_LMI_matrix(P, rho, nu, one):
             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, -(D.T@Q@D + D.T@S + S.T@D + R)))
+                np.hstack((off_diag.T, B.T@P@B-(D.T@Q@D + D.T@S + S.T@D + R)))
             ))
 
     n = sys.nstates
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index 0889d9746..343a286eb 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -46,14 +46,14 @@ def test_passivity_indices_ctime():
 
     assert(isinstance(nu, float))
 
-    sys_ff_nu = parallel(-nu, sys)
-    sys_fb_rho = feedback(rho, sys, sign=1)
+    sys_ff_nu = parallel(sys, -nu)
+    sys_fb_rho = feedback(sys, rho, sign=1)
 
     assert(sys_ff_nu.ispassive())
     assert(sys_fb_rho.ispassive())
 
-    sys_ff_nu = parallel(-nu-1e-6, sys)
-    sys_fb_rho = feedback(rho+1e-6, sys, sign=1)
+    sys_ff_nu = parallel(sys, -nu-1e-6)
+    sys_fb_rho = feedback(sys, rho+1e-6,  sign=1)
 
     assert(not sys_ff_nu.ispassive())
     assert(not sys_fb_rho.ispassive())
@@ -61,20 +61,20 @@ def test_passivity_indices_ctime():
 
 def test_passivity_indices_dtime():
     sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
-    sys = sample_system(sys, Ts=0.01, alpha=1, method="bilinear")
+    sys = sample_system(sys, Ts=0.1)
     nu = passivity.get_passivity_index(sys, 'input')
     rho = passivity.get_passivity_index(sys, 'output')
 
     assert(isinstance(nu, float))
 
-    sys_ff_nu = parallel(-nu, sys)
-    sys_fb_rho = feedback(rho, sys, sign=1)
+    sys_ff_nu = parallel(sys, -nu)
+    sys_fb_rho = feedback(sys, rho-1e-6,  sign=1)
 
     assert(sys_ff_nu.ispassive())
     assert(sys_fb_rho.ispassive())
 
-    sys_ff_nu = parallel(-nu-1e-6, sys)
-    sys_fb_rho = feedback(rho+1e-6, sys, sign=1)
+    sys_ff_nu = parallel(sys, -nu-1e-6)
+    sys_fb_rho = feedback(sys, rho, sign=1)
 
     assert(not sys_ff_nu.ispassive())
     assert(not sys_fb_rho.ispassive())

From c1fc4942522b8767df95750f9ca5c7f1a3185b17 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 25 Jul 2022 21:33:42 -0400
Subject: [PATCH 20/36] Change function signature on is_passive, split
 get_passivity_index into get_output_fb_index and get_input_ff_index.

---
 control/passivity.py            | 105 ++++++++++++++++----------------
 control/tests/passivity_test.py |  36 +++++------
 2 files changed, 69 insertions(+), 72 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index bba474a5a..2177839b3 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -51,7 +51,20 @@ def __P_pos_def_constraint__(n):
     return matrix_list
 
 
-def __ispassive__(sys, rho=None, nu=None):
+def __solve_passivity_LMI__(sys, rho=None, nu=None):
+    '''
+    Constructs a linear matrix inequality such that if a solution exists 
+    and the last element of the solution is positive, the system is passive.
+
+    The sources for the algorithm are: 
+
+    McCourt, Michael J., and Panos J. Antsaklis
+        "Demonstrating passivity and dissipativity using computational methods." 
+
+    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")
 
@@ -134,80 +147,64 @@ def make_LMI_matrix(P, rho, nu, one):
     return sol["x"]
 
 
-def ispassive(sys, rho=None, nu=None):
+def get_output_fb_index(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 sources for the algorithm are: 
+    Returns the output feedback passivity index for the input system. This is largest gain that can be placed in
+    positive feedback with a system such that the new interconnected system is passive.
+    Parameters
+    ----------
+    sys: An LTI system
+        System to be checked.
 
-    McCourt, Michael J., and Panos J. Antsaklis
-        "Demonstrating passivity and dissipativity using computational methods." 
+    Returns
+    -------
+    float: 
+        The OFP index 
+    '''
+    sol = __solve_passivity_LMI__(sys, nu=eps)
+    if sol is not None:
+        return sol[-1]
+    else:
+        return -np.inf
 
-    Nicholas Kottenstette and Panos J. Antsaklis
-        "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
-        equation 36.
 
+def get_input_ff_index(sys, index_type=None):
+    '''
+    Returns the input feedforward passivity (IFP) index for the input system. This 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: An LTI system
         System to be checked.
-    nu: float
-        Concrete value for input passivity index. 
-    rho: float
-        Concrete value for output passivity index. 
 
     Returns
     -------
-    bool, float, or None: 
-        The input system is passive, or the passivity index "opposite" the input. 
-        If the problem is unfeasiable when requesting the passivity index, returns None.
+    float: 
+        The IFP index 
     '''
-    if rho is not None and nu is not None:
-        return __ispassive__(sys, rho, nu) is not None
-    elif rho is None and nu is not None:
-        sol = __ispassive__(sys, nu=nu)
-        if sol is not None:
-            return sol[-1]
-        else:
-            return None
-    elif nu is None and rho is not None:
-        sol = __ispassive__(sys, rho=rho)
-        if sol is not None:
-            return sol[-1]
-        else:
-            return None
+
+    sol = __solve_passivity_LMI__(sys, rho=eps)
+    if sol is not None:
+        return sol[-1]
     else:
-        sol_rho = __ispassive__(sys, nu=eps)
-        sol_nu = __ispassive__(sys, rho=eps)
-        if sol_rho is not None and sol_nu is not None:
-            return sol_rho[-1] > 0 or sol_nu[-1] > 0
-        else:
-            return False
+        return -np.inf
 
 
-def get_passivity_index(sys, index_type=None):
+def ispassive(sys, rho=None, nu=None):
     '''
-    Returns the passivity index associated with the input string. 
+    Indicates if a linear time invariant (LTI) system is passive
+
     Parameters
     ----------
     sys: An LTI system
         System to be checked.
-    index_type: String
-        Must be 'input' or 'output'. Indicates which passivity index will be returned. 
 
     Returns
     -------
-    float: 
-        The passivity index 
+    bool, float, or None: 
+        The input system is passive, or the passivity index "opposite" the input. 
+        If the problem is unfeasiable when requesting the passivity index, returns None.
     '''
-    if index_type is None:
-        raise ControlArgument("Must provide index_type of 'input' or 'output'.")
-    if index_type == "input":
-        nu = ispassive(sys, rho=eps)
-        return nu
-    if index_type == "output":
-        rho = ispassive(sys, nu=eps)
-        return rho
+    output_fb_index = get_output_fb_index(sys)
+    input_ff_index = get_input_ff_index(sys)
+    return output_fb_index >= 0 or input_ff_index >= 0
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index 343a286eb..eaf7d417a 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -41,40 +41,40 @@ def test_ispassive_dtime():
 def test_passivity_indices_ctime():
     sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
 
-    nu = passivity.get_passivity_index(sys, 'input')
-    rho = passivity.get_passivity_index(sys, 'output')
+    iff_index = passivity.get_input_ff_index(sys)
+    ofb_index = passivity.get_output_fb_index(sys)
 
-    assert(isinstance(nu, float))
+    assert(isinstance(ofb_index, float))
 
-    sys_ff_nu = parallel(sys, -nu)
-    sys_fb_rho = feedback(sys, rho, sign=1)
+    sys_ff = parallel(sys, -iff_index)
+    sys_fb = feedback(sys, ofb_index, sign=1)
 
-    assert(sys_ff_nu.ispassive())
-    assert(sys_fb_rho.ispassive())
+    assert(sys_ff.ispassive())
+    assert(sys_fb.ispassive())
 
-    sys_ff_nu = parallel(sys, -nu-1e-6)
-    sys_fb_rho = feedback(sys, rho+1e-6,  sign=1)
+    sys_ff = parallel(sys, -iff_index-1e-6)
+    sys_fb = feedback(sys, ofb_index+1e-6,  sign=1)
 
-    assert(not sys_ff_nu.ispassive())
-    assert(not sys_fb_rho.ispassive())
+    assert(not sys_ff.ispassive())
+    assert(not sys_fb.ispassive())
 
 
 def test_passivity_indices_dtime():
     sys = tf([1, 1, 5, 0.1], [1, 2, 3, 4])
     sys = sample_system(sys, Ts=0.1)
-    nu = passivity.get_passivity_index(sys, 'input')
-    rho = passivity.get_passivity_index(sys, 'output')
+    iff_index = passivity.get_input_ff_index(sys)
+    ofb_index = passivity.get_output_fb_index(sys)
 
-    assert(isinstance(nu, float))
+    assert(isinstance(iff_index, float))
 
-    sys_ff_nu = parallel(sys, -nu)
-    sys_fb_rho = feedback(sys, rho-1e-6,  sign=1)
+    sys_ff_nu = parallel(sys, -iff_index)
+    sys_fb_rho = feedback(sys, ofb_index-1e-6,  sign=1)
 
     assert(sys_ff_nu.ispassive())
     assert(sys_fb_rho.ispassive())
 
-    sys_ff_nu = parallel(sys, -nu-1e-6)
-    sys_fb_rho = feedback(sys, rho, sign=1)
+    sys_ff_nu = parallel(sys, -iff_index-1e-6)
+    sys_fb_rho = feedback(sys, ofb_index, sign=1)
 
     assert(not sys_ff_nu.ispassive())
     assert(not sys_fb_rho.ispassive())

From ab8795098897f85f576d6e8579e2b00e1e3dc477 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 25 Jul 2022 21:41:04 -0400
Subject: [PATCH 21/36] Change weird double negative on conditional statements.

---
 control/passivity.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 2177839b3..946fa57d7 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -162,10 +162,10 @@ def get_output_fb_index(sys):
         The OFP index 
     '''
     sol = __solve_passivity_LMI__(sys, nu=eps)
-    if sol is not None:
-        return sol[-1]
-    else:
+    if sol is None:
         return -np.inf
+    else:
+        return sol[-1]
 
 
 def get_input_ff_index(sys, index_type=None):
@@ -184,10 +184,10 @@ def get_input_ff_index(sys, index_type=None):
     '''
 
     sol = __solve_passivity_LMI__(sys, rho=eps)
-    if sol is not None:
-        return sol[-1]
-    else:
+    if sol is None:
         return -np.inf
+    else:
+        return sol[-1]
 
 
 def ispassive(sys, rho=None, nu=None):

From e92f13ab94ff36919cb7fea1935ff6850031d6c3 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 25 Jul 2022 21:44:03 -0400
Subject: [PATCH 22/36] Make relevant passivity module function names conform
 to convention for weak "internal use" indication.

---
 control/passivity.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 946fa57d7..13b64355a 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -15,7 +15,7 @@
 eps = np.nextafter(0, 1)
 
 
-def __make_P_basis_matrices__(n, rho, nu, make_LMI_matrix_func):
+def _make_P_basis_matrices(n, rho, nu, make_LMI_matrix_func):
     '''
     Utility function to make basis matrices for a LMI from a 
     functional make_LMI_matrix_func and a symmetric matrix P of size n by n
@@ -34,7 +34,7 @@ def __make_P_basis_matrices__(n, rho, nu, make_LMI_matrix_func):
     return matrix_list
 
 
-def __P_pos_def_constraint__(n):
+def _P_pos_def_constraint(n):
     '''
     Utility function to make basis matrices for a LMI that ensures parametrized symbolic matrix 
     of size n by n is positive definite.
@@ -51,7 +51,7 @@ def __P_pos_def_constraint__(n):
     return matrix_list
 
 
-def __solve_passivity_LMI__(sys, rho=None, nu=None):
+def _solve_passivity_LMI(sys, rho=None, nu=None):
     '''
     Constructs a linear matrix inequality such that if a solution exists 
     and the last element of the solution is positive, the system is passive.
@@ -117,18 +117,18 @@ def make_LMI_matrix(P, rho, nu, one):
     )
 
     if rho is not None:
-        sys_matrix_list = __make_P_basis_matrices__(
+        sys_matrix_list = _make_P_basis_matrices(
             n, eps, 1.0, make_LMI_matrix)
         sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, eps, 1.0)
     else:
-        sys_matrix_list = __make_P_basis_matrices__(
+        sys_matrix_list = _make_P_basis_matrices(
             n, 1.0, eps, make_LMI_matrix)
         sys_constants = -make_LMI_matrix(np.zeros_like(A), eps, 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_matrix_list = _P_pos_def_constraint(n)
     P_coefficents = np.vstack(P_matrix_list).T
     P_constants = np.zeros((n, n))
 
@@ -161,7 +161,7 @@ def get_output_fb_index(sys):
     float: 
         The OFP index 
     '''
-    sol = __solve_passivity_LMI__(sys, nu=eps)
+    sol = _solve_passivity_LMI(sys, nu=eps)
     if sol is None:
         return -np.inf
     else:
@@ -183,7 +183,7 @@ def get_input_ff_index(sys, index_type=None):
         The IFP index 
     '''
 
-    sol = __solve_passivity_LMI__(sys, rho=eps)
+    sol = _solve_passivity_LMI(sys, rho=eps)
     if sol is None:
         return -np.inf
     else:

From e5c0f1b203a4188d58f50712db4ede028058ce32 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Mon, 25 Jul 2022 21:57:06 -0400
Subject: [PATCH 23/36] Add that last bit of test coverage.

---
 control/tests/passivity_test.py | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index eaf7d417a..5f6b40b38 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -6,7 +6,7 @@
 import numpy
 from control import ss, passivity, tf, sample_system, parallel, feedback
 from control.tests.conftest import cvxoptonly
-from control.exception import ControlDimension
+from control.exception import ControlArgument, ControlDimension
 
 pytestmark = cvxoptonly
 
@@ -124,6 +124,15 @@ def test_ispassive_all_zeros():
     with pytest.raises(ValueError):
         passivity.ispassive(sys)
 
+def test_rho_and_nu_are_none():
+    A = numpy.array([[0]])
+    B = numpy.array([[0]])
+    C = numpy.array([[0]])
+    D = numpy.array([[0]])
+    sys = ss(A, B, C, D)
+
+    with pytest.raises(ControlArgument):
+        passivity._solve_passivity_LMI(sys)
 
 def test_transfer_function():
     sys = tf([1], [1, 2])

From 8bd658f680f2c6d213d098cddaf6d6d3096941d3 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Tue, 26 Jul 2022 21:13:51 -0400
Subject: [PATCH 24/36] Clean up some unused function inputs, update doc
 string.

---
 control/passivity.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 13b64355a..74a4e40ac 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -168,7 +168,7 @@ def get_output_fb_index(sys):
         return sol[-1]
 
 
-def get_input_ff_index(sys, index_type=None):
+def get_input_ff_index(sys):
     '''
     Returns the input feedforward passivity (IFP) index for the input system. This is the largest gain that can be 
     placed in negative parallel interconnection with a system such that the new interconnected system is passive.
@@ -190,7 +190,7 @@ def get_input_ff_index(sys, index_type=None):
         return sol[-1]
 
 
-def ispassive(sys, rho=None, nu=None):
+def ispassive(sys):
     '''
     Indicates if a linear time invariant (LTI) system is passive
 
@@ -201,9 +201,8 @@ def ispassive(sys, rho=None, nu=None):
 
     Returns
     -------
-    bool, float, or None: 
-        The input system is passive, or the passivity index "opposite" the input. 
-        If the problem is unfeasiable when requesting the passivity index, returns None.
+    bool: 
+        The input system is passive.
     '''
     output_fb_index = get_output_fb_index(sys)
     input_ff_index = get_input_ff_index(sys)

From ffd0837c54da3bafc4fc1046f88386a784001bc8 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Tue, 26 Jul 2022 22:17:29 -0400
Subject: [PATCH 25/36] Add passivity module to __init__.py and Sphinx
 documentation.

---
 control/__init__.py | 1 +
 doc/control.rst     | 1 +
 2 files changed, 2 insertions(+)

diff --git a/control/__init__.py b/control/__init__.py
index ad2685273..e0edc96a2 100644
--- a/control/__init__.py
+++ b/control/__init__.py
@@ -72,6 +72,7 @@
 from .config import *
 from .sisotool import *
 from .iosys import *
+from .passivity import *
 
 # Exceptions
 from .exception import *
diff --git a/doc/control.rst b/doc/control.rst
index fc6618d24..d067ab93c 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -88,6 +88,7 @@ Control system analysis
     pzmap
     root_locus
     sisotool
+    passivity
 
 Matrix computations
 ===================

From 9234188ef96f4f02c5ad92b58891d271eae55c2f Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Tue, 26 Jul 2022 22:21:17 -0400
Subject: [PATCH 26/36] Fix import name collision.

---
 control/passivity.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index 74a4e40ac..f5b07b550 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -4,7 +4,7 @@
 '''
 
 import numpy as np
-from control import statesp as ss
+from control import statesp
 from control.exception import ControlArgument, ControlDimension
 
 try:
@@ -75,7 +75,7 @@ def _solve_passivity_LMI(sys, rho=None, nu=None):
     if rho is None and nu is None:
         raise ControlArgument("rho or nu must be given a float value.")
 
-    sys = ss._convert_to_statespace(sys)
+    sys = statesp._convert_to_statespace(sys)
 
     A = sys.A
     B = sys.B

From 5a3245f6217a0d8e1bb115416e0201907338c559 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Wed, 27 Jul 2022 12:03:19 -0400
Subject: [PATCH 27/36] Improve sphinx docs and function doc strings.

---
 control/passivity.py | 36 ++++++++++++++++++++++--------------
 doc/control.rst      |  6 +++++-
 2 files changed, 27 insertions(+), 15 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index f5b07b550..c92a7c9ac 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -14,9 +14,11 @@
 
 eps = np.nextafter(0, 1)
 
+__all__ = ["get_output_fb_index", "get_input_ff_index", "ispassive"]
 
 def _make_P_basis_matrices(n, rho, nu, make_LMI_matrix_func):
-    '''
+    '''Make list of matrix constraints for passivity LMI
+
     Utility function to make basis matrices for a LMI from a 
     functional make_LMI_matrix_func and a symmetric matrix P of size n by n
     representing a parametrized symbolic matrix
@@ -35,7 +37,8 @@ def _make_P_basis_matrices(n, rho, nu, make_LMI_matrix_func):
 
 
 def _P_pos_def_constraint(n):
-    '''
+    '''Make list of matrix constraints for P >= 0
+
     Utility function to make basis matrices for a LMI that ensures parametrized symbolic matrix 
     of size n by n is positive definite.
     '''
@@ -52,9 +55,11 @@ def _P_pos_def_constraint(n):
 
 
 def _solve_passivity_LMI(sys, rho=None, nu=None):
-    '''
-    Constructs a linear matrix inequality such that if a solution exists 
-    and the last element of the solution is positive, the system is passive.
+    '''Compute passivity indices via a linear matrix inequality (LMI)
+
+    Constructs an LMI such that if a solution exists and the last element of the 
+    solution is positive, the system is passive. The last element is either the 
+    input or output passivity index, for nu=None and rho=None respectively.
 
     The sources for the algorithm are: 
 
@@ -124,7 +129,7 @@ def make_LMI_matrix(P, rho, nu, one):
         sys_matrix_list = _make_P_basis_matrices(
             n, 1.0, eps, make_LMI_matrix)
         sys_constants = -make_LMI_matrix(np.zeros_like(A), eps, nu, 1.0)
-
+    
     sys_coefficents = np.vstack(sys_matrix_list).T
 
     # LMI to ensure P is positive definite
@@ -148,9 +153,11 @@ def make_LMI_matrix(P, rho, nu, one):
 
 
 def get_output_fb_index(sys):
-    '''
-    Returns the output feedback passivity index for the input system. This is largest gain that can be placed in
-    positive feedback with a system such that the new interconnected system is passive.
+    '''Returns the output feedback passivity (OFP) index for the input system. 
+    
+    The OFP is largest gain that can be placed in positive feedback 
+    with a system such that the new interconnected system is passive.
+
     Parameters
     ----------
     sys: An LTI system
@@ -169,9 +176,11 @@ def get_output_fb_index(sys):
 
 
 def get_input_ff_index(sys):
-    '''
-    Returns the input feedforward passivity (IFP) index for the input system. This is the largest gain that can be 
-    placed in negative parallel interconnection with a system such that the new interconnected system is passive.
+    '''Returns the input feedforward passivity (IFP) index for the input system. 
+    
+    The 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: An LTI system
@@ -191,8 +200,7 @@ def get_input_ff_index(sys):
 
 
 def ispassive(sys):
-    '''
-    Indicates if a linear time invariant (LTI) system is passive
+    '''Indicates if a linear time invariant (LTI) system is passive
 
     Parameters
     ----------
diff --git a/doc/control.rst b/doc/control.rst
index d067ab93c..172790f83 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -80,6 +80,9 @@ Control system analysis
     describing_function
     evalfr
     freqresp
+    get_input_ff_index
+    get_output_fb_index
+    ispassive
     margin
     stability_margins
     phase_crossover_frequencies
@@ -88,7 +91,8 @@ Control system analysis
     pzmap
     root_locus
     sisotool
-    passivity
+
+
 
 Matrix computations
 ===================

From 01b55913eba04a478012a4613413fd0a3cf4bf31 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Wed, 27 Jul 2022 23:03:07 -0400
Subject: [PATCH 28/36] Pull some helper functions into _solve_passivity_LMI.
 Add more detailed and clear comments and doc strings. Fix correctness of
 _solve_passivity_LMI for inputing values for both input and output passivity
 indices. Allow checking of passivity using ispassive given values for input
 and ouput passivity indices.

---
 control/passivity.py            | 124 ++++++++++++++++----------------
 control/tests/passivity_test.py |  30 +++-----
 2 files changed, 74 insertions(+), 80 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index c92a7c9ac..4081ded2b 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -16,50 +16,15 @@
 
 __all__ = ["get_output_fb_index", "get_input_ff_index", "ispassive"]
 
-def _make_P_basis_matrices(n, rho, nu, make_LMI_matrix_func):
-    '''Make list of matrix constraints for passivity LMI
-
-    Utility function to make basis matrices for a LMI from a 
-    functional make_LMI_matrix_func and a symmetric matrix P of size n by n
-    representing a parametrized 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.0
-                P[j, i] = 1.0
-                matrix_list.append(make_LMI_matrix_func(P, 0, 0, 0).flatten())
-    P = eps*np.eye(n)
-    matrix_list.append(make_LMI_matrix_func(P, rho, nu, 0).flatten())
-    return matrix_list
-
-
-def _P_pos_def_constraint(n):
-    '''Make list of matrix constraints for P >= 0
-
-    Utility function to make basis matrices for a LMI that ensures parametrized 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.0
-                P[j, i] = -1.0
-                matrix_list.append(P.flatten())
-    matrix_list.append(np.zeros((n, n)).flatten())
-    return matrix_list
-
 
 def _solve_passivity_LMI(sys, rho=None, nu=None):
-    '''Compute passivity indices via a linear matrix inequality (LMI)
+    '''Computes passivity indices and/or solves feasiblity via a linear matrix inequality (LMI).
 
     Constructs an LMI such that if a solution exists and the last element of the 
-    solution is positive, the system is passive. The last element is either the 
-    input or output passivity index, for nu=None and rho=None respectively.
+    solution is positive, the system is passive. Inputs of None for rho or nu indicates 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 bilinear matrix inequality.) The last element is either the 
+    output or input passivity index, for rho=None and nu=None respectively.
 
     The sources for the algorithm are: 
 
@@ -78,7 +43,7 @@ def _solve_passivity_LMI(sys, rho=None, nu=None):
             "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 float value.")
+        raise ControlArgument("rho or nu must be given a numerical value.")
 
     sys = statesp._convert_to_statespace(sys)
 
@@ -112,28 +77,64 @@ def make_LMI_matrix(P, rho, nu, one):
                 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):
+        '''Makes 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 parametrized 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 = eps*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):
+        '''Makes a list of matrix constraints for P >= 0.
+
+        Utility function to make basis matrices for a LMI that ensures parametrized 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
 
-    # LMI for passivity indices from A,B,C,D
-    sys_matrix_list = list()
-    sys_constants = -np.vstack((
-        np.hstack((np.zeros_like(A),  np.zeros_like(C.T))),
-        np.hstack((np.zeros_like(C), np.zeros_like(D))))
-    )
+    # coefficents for passivity indices and feasibility matrix 
+    sys_matrix_list = make_P_basis_matrices(n, rho, nu)
 
-    if rho is not None:
-        sys_matrix_list = _make_P_basis_matrices(
-            n, eps, 1.0, make_LMI_matrix)
+    # 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, eps, 1.0)
-    else:
-        sys_matrix_list = _make_P_basis_matrices(
-            n, 1.0, eps, make_LMI_matrix)
+    elif nu is not None:
         sys_constants = -make_LMI_matrix(np.zeros_like(A), eps, 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_matrix_list = P_pos_def_constraint(n)
     P_coefficents = np.vstack(P_matrix_list).T
     P_constants = np.zeros((n, n))
 
@@ -141,7 +142,10 @@ def make_LMI_matrix(P, rho, nu, one):
     number_of_opt_vars = int(
         (n**2-n)/2 + n)
     c = cvx.matrix(0.0, (number_of_opt_vars, 1))
-    c = cvx.matrix(np.append(np.array(c), -1.0))
+
+    #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)]
@@ -155,7 +159,7 @@ def make_LMI_matrix(P, rho, nu, one):
 def get_output_fb_index(sys):
     '''Returns the output feedback passivity (OFP) index for the input system. 
     
-    The OFP is largest gain that can be placed in positive feedback 
+    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
@@ -199,8 +203,8 @@ def get_input_ff_index(sys):
         return sol[-1]
 
 
-def ispassive(sys):
-    '''Indicates if a linear time invariant (LTI) system is passive
+def ispassive(sys, ofp_index = 0, ifp_index = 0):
+    '''Indicates if a linear time invariant (LTI) system is passive.
 
     Parameters
     ----------
@@ -212,6 +216,4 @@ def ispassive(sys):
     bool: 
         The input system is passive.
     '''
-    output_fb_index = get_output_fb_index(sys)
-    input_ff_index = get_input_ff_index(sys)
-    return output_fb_index >= 0 or input_ff_index >= 0
+    return _solve_passivity_LMI(sys, rho = ofp_index, nu = ifp_index) is not None
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index 5f6b40b38..d7fdb29cf 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -67,17 +67,17 @@ def test_passivity_indices_dtime():
 
     assert(isinstance(iff_index, float))
 
-    sys_ff_nu = parallel(sys, -iff_index)
-    sys_fb_rho = feedback(sys, ofb_index-1e-6,  sign=1)
+    sys_ff = parallel(sys, -iff_index)
+    sys_fb = feedback(sys, ofb_index,  sign=1)
 
-    assert(sys_ff_nu.ispassive())
-    assert(sys_fb_rho.ispassive())
+    assert(sys_ff.ispassive())
+    assert(sys_fb.ispassive())
 
-    sys_ff_nu = parallel(sys, -iff_index-1e-6)
-    sys_fb_rho = feedback(sys, ofb_index, sign=1)
+    sys_ff = parallel(sys, -iff_index-1e-2)
+    sys_fb = feedback(sys, ofb_index+1e-2, sign=1)
 
-    assert(not sys_ff_nu.ispassive())
-    assert(not sys_fb_rho.ispassive())
+    assert(not sys_ff.ispassive())
+    assert(not sys_fb.ispassive())
 
 
 def test_system_dimension():
@@ -104,7 +104,8 @@ def test_system_dimension():
      ((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, C, D), True),
+     ((A*0, B*0, C*0, D*0), True)])
 def test_ispassive_edge_cases(test_input, expected):
     A = test_input[0]
     B = test_input[1]
@@ -114,16 +115,6 @@ def test_ispassive_edge_cases(test_input, expected):
     assert(passivity.ispassive(sys) == expected)
 
 
-def test_ispassive_all_zeros():
-    A = numpy.array([[0]])
-    B = numpy.array([[0]])
-    C = numpy.array([[0]])
-    D = numpy.array([[0]])
-    sys = ss(A, B, C, D)
-
-    with pytest.raises(ValueError):
-        passivity.ispassive(sys)
-
 def test_rho_and_nu_are_none():
     A = numpy.array([[0]])
     B = numpy.array([[0]])
@@ -134,6 +125,7 @@ def test_rho_and_nu_are_none():
     with pytest.raises(ControlArgument):
         passivity._solve_passivity_LMI(sys)
 
+
 def test_transfer_function():
     sys = tf([1], [1, 2])
     assert(passivity.ispassive(sys))

From d9cb82860e19ad0be4624bf70019d91583375a1e Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Wed, 27 Jul 2022 23:09:51 -0400
Subject: [PATCH 29/36] Add function I/O description to doc string for
 _solve_passivity_LMI.

---
 control/passivity.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/control/passivity.py b/control/passivity.py
index 4081ded2b..d9caa3db2 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -34,6 +34,20 @@ def _solve_passivity_LMI(sys, rho=None, nu=None):
     Nicholas Kottenstette and Panos J. Antsaklis
         "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
         equation 36.
+
+    Parameters
+    ----------
+    sys: An LTI system
+        System to be checked.
+    rho: Float or None
+        Output feedback passivity index
+    nu: Float or None
+        Input feedforward passivity index
+        
+    Returns
+    -------
+    numpy array: 
+        The LMI solution
     '''
     if cvx is None:
         raise ModuleNotFoundError("cvxopt required for passivity module")

From 32e9b24319f524a56c09dcb5e2c458fa7c67969e Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Wed, 27 Jul 2022 23:36:06 -0400
Subject: [PATCH 30/36] Change _solve_passivity_LMI to solve_passivity_LMI,
 include in __all__. Add even more detail to doc strings.

---
 control/passivity.py            | 29 ++++++++++++++++++++++-------
 control/tests/passivity_test.py |  2 +-
 2 files changed, 23 insertions(+), 8 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index d9caa3db2..ca094901d 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -14,16 +14,17 @@
 
 eps = np.nextafter(0, 1)
 
-__all__ = ["get_output_fb_index", "get_input_ff_index", "ispassive"]
+__all__ = ["get_output_fb_index", "get_input_ff_index", 
+    "ispassive", "solve_passivity_LMI"]
 
 
-def _solve_passivity_LMI(sys, rho=None, nu=None):
+def solve_passivity_LMI(sys, rho=None, nu=None):
     '''Computes passivity indices and/or solves feasiblity via a linear matrix inequality (LMI).
 
     Constructs an LMI such that if a solution exists and the last element of the 
     solution is positive, the system is passive. Inputs of None for rho or nu indicates 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 bilinear matrix inequality.) The last element is either the 
+    otherwise you're trying to solve a nonconvex bilinear matrix inequality.) The last element is either the 
     output or input passivity index, for rho=None and nu=None respectively.
 
     The sources for the algorithm are: 
@@ -186,7 +187,7 @@ def get_output_fb_index(sys):
     float: 
         The OFP index 
     '''
-    sol = _solve_passivity_LMI(sys, nu=eps)
+    sol = solve_passivity_LMI(sys, nu=eps)
     if sol is None:
         return -np.inf
     else:
@@ -210,7 +211,7 @@ def get_input_ff_index(sys):
         The IFP index 
     '''
 
-    sol = _solve_passivity_LMI(sys, rho=eps)
+    sol = solve_passivity_LMI(sys, rho=eps)
     if sol is None:
         return -np.inf
     else:
@@ -220,14 +221,28 @@ def get_input_ff_index(sys):
 def ispassive(sys, ofp_index = 0, ifp_index = 0):
     '''Indicates 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. Querying if the system is passive in the sense of V(x)>=0 and \\dot{V}(x) <= y.T*u, 
+    is equiavlent to the default case of ofp_index = 0, ifp_index = 0.
+
+    Note that computing the ofp_index and ifp_index for a system, then using both values simultaneously 
+    to 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:
+        McCourt, Michael J., and Panos J. Antsaklis
+            "Demonstrating passivity and dissipativity using computational methods." 
+
     Parameters
     ----------
     sys: An LTI system
         System to be checked.
+    ofp_index: float
+        Output feedback passivity index.
+    ifp_index: float
+        Input feedforward passivity index.
 
     Returns
     -------
     bool: 
-        The input system is passive.
+        The system is passive.
     '''
-    return _solve_passivity_LMI(sys, rho = ofp_index, nu = ifp_index) is not None
+    return solve_passivity_LMI(sys, rho = ofp_index, nu = ifp_index) is not None
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
index d7fdb29cf..4c95c96b9 100644
--- a/control/tests/passivity_test.py
+++ b/control/tests/passivity_test.py
@@ -123,7 +123,7 @@ def test_rho_and_nu_are_none():
     sys = ss(A, B, C, D)
 
     with pytest.raises(ControlArgument):
-        passivity._solve_passivity_LMI(sys)
+        passivity.solve_passivity_LMI(sys)
 
 
 def test_transfer_function():

From f4809d982927d45ad3d1c52bb552ceaa72b26162 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Wed, 27 Jul 2022 23:54:04 -0400
Subject: [PATCH 31/36] Change functions to compute passivity indices to return
 RuntimeErrors instead of -inf when infeasible. (I'm not sure how to trigger
 this code for unit testing, mostly a just-in-case.

---
 control/passivity.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index ca094901d..5c9f6e74f 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -189,7 +189,7 @@ def get_output_fb_index(sys):
     '''
     sol = solve_passivity_LMI(sys, nu=eps)
     if sol is None:
-        return -np.inf
+        raise RuntimeError("LMI passivity problem is infeasible")
     else:
         return sol[-1]
 
@@ -213,7 +213,7 @@ def get_input_ff_index(sys):
 
     sol = solve_passivity_LMI(sys, rho=eps)
     if sol is None:
-        return -np.inf
+        raise RuntimeError("LMI passivity problem is infeasible")
     else:
         return sol[-1]
 

From 4a95fa33daf1007f9895b6de401b175553756175 Mon Sep 17 00:00:00 2001
From: Mark <mark.x.yeatman@gmail.com>
Date: Thu, 28 Jul 2022 14:59:14 -0400
Subject: [PATCH 32/36] Add reminders for future implementation.

---
 control/passivity.py | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/control/passivity.py b/control/passivity.py
index 5c9f6e74f..b5d6b29fb 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -218,6 +218,30 @@ def get_input_ff_index(sys):
         return sol[-1]
 
 
+def get_relative_index(sys):
+    """Returns the relative passivity index for the system.
+
+    (not implemented yet)
+    """
+    raise NotImplemented("Relative passivity index not implemented")
+
+
+def get_combined_io_index(sys):
+    """Returns the combined I/O passivity index for the system.
+
+    (not implemented yet)
+    """
+    raise NotImplemented("Combined I/O passivity index not implemented")
+
+
+def get_directional_index(sys):
+    """Returns the directional passivity index for the system.
+
+    (not implemented yet)
+    """
+    raise NotImplemented("Directional passivity index not implemented")
+
+    
 def ispassive(sys, ofp_index = 0, ifp_index = 0):
     '''Indicates if a linear time invariant (LTI) system is passive.
 

From 2c2d25af7f49ac22b6fdf539cb7d09920305fb6d Mon Sep 17 00:00:00 2001
From: Mark <mark.x.yeatman@gmail.com>
Date: Thu, 28 Jul 2022 15:26:58 -0400
Subject: [PATCH 33/36] Update documentation so that math renders and
 organization is better.

---
 control/passivity.py | 37 ++++++++++++++++++++++++-------------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index b5d6b29fb..6995608b8 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -243,18 +243,11 @@ def get_directional_index(sys):
 
     
 def ispassive(sys, ofp_index = 0, ifp_index = 0):
-    '''Indicates 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. Querying if the system is passive in the sense of V(x)>=0 and \\dot{V}(x) <= y.T*u, 
-    is equiavlent to the default case of ofp_index = 0, ifp_index = 0.
-
-    Note that computing the ofp_index and ifp_index for a system, then using both values simultaneously 
-    to 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:
-        McCourt, Michael J., and Panos J. Antsaklis
-            "Demonstrating passivity and dissipativity using computational methods." 
-
+    '''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: An LTI system
@@ -263,10 +256,28 @@ def ispassive(sys, ofp_index = 0, ifp_index = 0):
         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`, `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." 
     '''
     return solve_passivity_LMI(sys, rho = ofp_index, nu = ifp_index) is not None

From 0b0d6921c682bd80ab89b3c9b312900d01c91c73 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Fri, 29 Jul 2022 20:38:22 -0400
Subject: [PATCH 34/36] Fix doc string.

---
 control/passivity.py |   4 +-
 passivity_test.ipynb | 105 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 107 insertions(+), 2 deletions(-)
 create mode 100644 passivity_test.ipynb

diff --git a/control/passivity.py b/control/passivity.py
index 6995608b8..ee757c4f5 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -243,7 +243,7 @@ def get_directional_index(sys):
 
     
 def ispassive(sys, ofp_index = 0, ifp_index = 0):
-    '''Indicate if a linear time invariant (LTI) system is passive.
+    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.
@@ -279,5 +279,5 @@ def ispassive(sys, ofp_index = 0, ifp_index = 0):
     ----------
     .. [1] McCourt, Michael J., and Panos J. Antsaklis
           "Demonstrating passivity and dissipativity using computational methods." 
-    '''
+    """
     return solve_passivity_LMI(sys, rho = ofp_index, nu = ifp_index) is not None
diff --git a/passivity_test.ipynb b/passivity_test.ipynb
new file mode 100644
index 000000000..48a538566
--- /dev/null
+++ b/passivity_test.ipynb
@@ -0,0 +1,105 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-inf\n",
+      "-inf\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "c:\\Users\\Mark\\anaconda3\\lib\\site-packages\\scipy\\signal\\lti_conversion.py:109: RuntimeWarning: invalid value encountered in subtract\n",
+      "  C = num[:, 1:] - outer(num[:, 0], den[1:])\n",
+      "c:\\Users\\Mark\\python-control\\control\\passivity.py:77: RuntimeWarning: invalid value encountered in matmul\n",
+      "  np.hstack((off_diag.T, B.T@P@B-(D.T@Q@D + D.T@S + S.T@D + R)))\n",
+      "c:\\Users\\Mark\\python-control\\control\\passivity.py:77: RuntimeWarning: invalid value encountered in add\n",
+      "  np.hstack((off_diag.T, B.T@P@B-(D.T@Q@D + D.T@S + S.T@D + R)))\n"
+     ]
+    },
+    {
+     "ename": "ArithmeticError",
+     "evalue": "1",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[1;31mArithmeticError\u001b[0m                           Traceback (most recent call last)",
+      "\u001b[1;32mc:\\Users\\Mark\\python-control\\passivity_test.ipynb Cell 1\u001b[0m in \u001b[0;36m<cell line: 19>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=15'>16</a>\u001b[0m sys_ff_nu \u001b[39m=\u001b[39m parallel(sys, \u001b[39m-\u001b[39miff_index)\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=16'>17</a>\u001b[0m sys_fb_rho \u001b[39m=\u001b[39m feedback(sys, ofb_index,  sign\u001b[39m=\u001b[39m\u001b[39m1\u001b[39m)\n\u001b[1;32m---> <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=18'>19</a>\u001b[0m \u001b[39massert\u001b[39;00m(sys_ff_nu\u001b[39m.\u001b[39;49mispassive())\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=19'>20</a>\u001b[0m \u001b[39massert\u001b[39;00m(sys_fb_rho\u001b[39m.\u001b[39mispassive())\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=21'>22</a>\u001b[0m sys_ff_nu \u001b[39m=\u001b[39m parallel(sys, \u001b[39m-\u001b[39miff_index\u001b[39m-\u001b[39m\u001b[39m1e-2\u001b[39m)\n",
+      "File \u001b[1;32mc:\\Users\\Mark\\python-control\\control\\lti.py:208\u001b[0m, in \u001b[0;36mLTI.ispassive\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    205\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mispassive\u001b[39m(\u001b[39mself\u001b[39m):\n\u001b[0;32m    206\u001b[0m     \u001b[39m# importing here prevents circular dependancy\u001b[39;00m\n\u001b[0;32m    207\u001b[0m     \u001b[39mfrom\u001b[39;00m \u001b[39mcontrol\u001b[39;00m\u001b[39m.\u001b[39;00m\u001b[39mpassivity\u001b[39;00m \u001b[39mimport\u001b[39;00m ispassive\n\u001b[1;32m--> 208\u001b[0m     \u001b[39mreturn\u001b[39;00m ispassive(\u001b[39mself\u001b[39;49m)\n",
+      "File \u001b[1;32mc:\\Users\\Mark\\python-control\\control\\passivity.py:219\u001b[0m, in \u001b[0;36mispassive\u001b[1;34m(sys, ofp_index, ifp_index)\u001b[0m\n\u001b[0;32m    206\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mispassive\u001b[39m(sys, ofp_index \u001b[39m=\u001b[39m \u001b[39m0\u001b[39m, ifp_index \u001b[39m=\u001b[39m \u001b[39m0\u001b[39m):\n\u001b[0;32m    207\u001b[0m     \u001b[39m'''Indicates if a linear time invariant (LTI) system is passive.\u001b[39;00m\n\u001b[0;32m    208\u001b[0m \n\u001b[0;32m    209\u001b[0m \u001b[39m    Parameters\u001b[39;00m\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    217\u001b[0m \u001b[39m        The input system is passive.\u001b[39;00m\n\u001b[0;32m    218\u001b[0m \u001b[39m    '''\u001b[39;00m\n\u001b[1;32m--> 219\u001b[0m     \u001b[39mreturn\u001b[39;00m _solve_passivity_LMI(sys, rho \u001b[39m=\u001b[39;49m ofp_index, nu \u001b[39m=\u001b[39;49m ifp_index) \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m\n",
+      "File \u001b[1;32mc:\\Users\\Mark\\python-control\\control\\passivity.py:155\u001b[0m, in \u001b[0;36m_solve_passivity_LMI\u001b[1;34m(sys, rho, nu)\u001b[0m\n\u001b[0;32m    153\u001b[0m \u001b[39m# crunch feasibility solution\u001b[39;00m\n\u001b[0;32m    154\u001b[0m cvx\u001b[39m.\u001b[39msolvers\u001b[39m.\u001b[39moptions[\u001b[39m'\u001b[39m\u001b[39mshow_progress\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m \u001b[39mFalse\u001b[39;00m\n\u001b[1;32m--> 155\u001b[0m sol \u001b[39m=\u001b[39m cvx\u001b[39m.\u001b[39;49msolvers\u001b[39m.\u001b[39;49msdp(c, Gs\u001b[39m=\u001b[39;49mGs, hs\u001b[39m=\u001b[39;49mhs)\n\u001b[0;32m    156\u001b[0m \u001b[39mreturn\u001b[39;00m sol[\u001b[39m\"\u001b[39m\u001b[39mx\u001b[39m\u001b[39m\"\u001b[39m]\n",
+      "File \u001b[1;32mc:\\Users\\Mark\\anaconda3\\lib\\site-packages\\cvxopt\\coneprog.py:4126\u001b[0m, in \u001b[0;36msdp\u001b[1;34m(c, Gl, hl, Gs, hs, A, b, kktsolver, solver, primalstart, dualstart, **kwargs)\u001b[0m\n\u001b[0;32m   4123\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m   4124\u001b[0m     ds \u001b[39m=\u001b[39m \u001b[39mNone\u001b[39;00m\n\u001b[1;32m-> 4126\u001b[0m sol \u001b[39m=\u001b[39m conelp(c, G, h, dims, A \u001b[39m=\u001b[39;49m A, b \u001b[39m=\u001b[39;49m b, primalstart \u001b[39m=\u001b[39;49m ps, dualstart \u001b[39m=\u001b[39;49m ds, kktsolver \u001b[39m=\u001b[39;49m kktsolver, options \u001b[39m=\u001b[39;49m options)\n\u001b[0;32m   4127\u001b[0m \u001b[39mif\u001b[39;00m sol[\u001b[39m'\u001b[39m\u001b[39ms\u001b[39m\u001b[39m'\u001b[39m] \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m   4128\u001b[0m     sol[\u001b[39m'\u001b[39m\u001b[39msl\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m \u001b[39mNone\u001b[39;00m\n",
+      "File \u001b[1;32mc:\\Users\\Mark\\anaconda3\\lib\\site-packages\\cvxopt\\coneprog.py:1033\u001b[0m, in \u001b[0;36mconelp\u001b[1;34m(c, G, h, dims, A, b, primalstart, dualstart, kktsolver, xnewcopy, xdot, xaxpy, xscal, ynewcopy, ydot, yaxpy, yscal, **kwargs)\u001b[0m\n\u001b[0;32m   1026\u001b[0m \u001b[39m# Compute initial scaling W:\u001b[39;00m\n\u001b[0;32m   1027\u001b[0m \u001b[39m#\u001b[39;00m\n\u001b[0;32m   1028\u001b[0m \u001b[39m#     W * z = W^{-T} * s = lambda\u001b[39;00m\n\u001b[0;32m   1029\u001b[0m \u001b[39m#     dg * tau = 1/dg * kappa = lambdag.\u001b[39;00m\n\u001b[0;32m   1031\u001b[0m \u001b[39mif\u001b[39;00m iters \u001b[39m==\u001b[39m \u001b[39m0\u001b[39m:\n\u001b[1;32m-> 1033\u001b[0m     W \u001b[39m=\u001b[39m misc\u001b[39m.\u001b[39;49mcompute_scaling(s, z, lmbda, dims, mnl \u001b[39m=\u001b[39;49m \u001b[39m0\u001b[39;49m)\n\u001b[0;32m   1035\u001b[0m     \u001b[39m#     dg = sqrt( kappa / tau )\u001b[39;00m\n\u001b[0;32m   1036\u001b[0m     \u001b[39m#     dgi = sqrt( tau / kappa )\u001b[39;00m\n\u001b[0;32m   1037\u001b[0m     \u001b[39m#     lambda_g = sqrt( tau * kappa )\u001b[39;00m\n\u001b[0;32m   1038\u001b[0m     \u001b[39m#\u001b[39;00m\n\u001b[0;32m   1039\u001b[0m     \u001b[39m# lambda_g is stored in the last position of lmbda.\u001b[39;00m\n\u001b[0;32m   1041\u001b[0m     dg \u001b[39m=\u001b[39m math\u001b[39m.\u001b[39msqrt( kappa \u001b[39m/\u001b[39m tau )\n",
+      "File \u001b[1;32mc:\\Users\\Mark\\anaconda3\\lib\\site-packages\\cvxopt\\misc.py:387\u001b[0m, in \u001b[0;36mcompute_scaling\u001b[1;34m(s, z, lmbda, dims, mnl)\u001b[0m\n\u001b[0;32m    385\u001b[0m \u001b[39m# Factor sk = Ls*Ls'; store Ls in ds[inds[k]:inds[k+1]].\u001b[39;00m\n\u001b[0;32m    386\u001b[0m blas\u001b[39m.\u001b[39mcopy(s, Ls, offsetx \u001b[39m=\u001b[39m ind2, n \u001b[39m=\u001b[39m m\u001b[39m*\u001b[39m\u001b[39m*\u001b[39m\u001b[39m2\u001b[39m) \n\u001b[1;32m--> 387\u001b[0m lapack\u001b[39m.\u001b[39;49mpotrf(Ls, n \u001b[39m=\u001b[39;49m m, ldA \u001b[39m=\u001b[39;49m m)\n\u001b[0;32m    389\u001b[0m \u001b[39m# Factor zs[k] = Lz*Lz'; store Lz in dz[inds[k]:inds[k+1]].\u001b[39;00m\n\u001b[0;32m    390\u001b[0m blas\u001b[39m.\u001b[39mcopy(z, Lz, offsetx \u001b[39m=\u001b[39m ind2, n \u001b[39m=\u001b[39m m\u001b[39m*\u001b[39m\u001b[39m*\u001b[39m\u001b[39m2\u001b[39m) \n",
+      "\u001b[1;31mArithmeticError\u001b[0m: 1"
+     ]
+    }
+   ],
+   "source": [
+    "import matplotlib.pyplot \n",
+    "import numpy\n",
+    "import math\n",
+    "import sys as system\n",
+    "from control import passivity, ss, tf, parallel, feedback, sample_system\n",
+    "from scipy import optimize\n",
+    "\n",
+    "sys = tf([1], [1, -2])\n",
+    "print(sys)\n",
+    "sys = sample_system(sys, Ts=1)\n",
+    "iff_index = passivity.get_input_ff_index(sys)\n",
+    "ofb_index = passivity.get_output_fb_index(sys)\n",
+    "\n",
+    "print(iff_index)\n",
+    "print(ofb_index)\n",
+    "\n",
+    "sys_ff_nu = parallel(sys, -iff_index)\n",
+    "sys_fb_rho = feedback(sys, ofb_index,  sign=1)\n",
+    "\n",
+    "assert(sys_ff_nu.ispassive())\n",
+    "assert(sys_fb_rho.ispassive())\n",
+    "\n",
+    "sys_ff_nu = parallel(sys, -iff_index-1e-2)\n",
+    "sys_fb_rho = feedback(sys, ofb_index+1e-2, sign=1)\n",
+    "\n",
+    "assert(not sys_ff_nu.ispassive())\n",
+    "assert(not sys_fb_rho.ispassive())\n",
+    "\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.9.12 ('base')",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.12"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "ccecdc02aa416300a64487fe6971a8398ebaf053d2b2768f769ac21dbbf13f86"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

From fdaab428aedfe4c08d9e37537f2082aa7a295ef2 Mon Sep 17 00:00:00 2001
From: mark-yeatman <mark.x.yeatman@gmail.com>
Date: Fri, 29 Jul 2022 20:40:30 -0400
Subject: [PATCH 35/36] Remove test code.

---
 passivity_test.ipynb | 105 -------------------------------------------
 1 file changed, 105 deletions(-)
 delete mode 100644 passivity_test.ipynb

diff --git a/passivity_test.ipynb b/passivity_test.ipynb
deleted file mode 100644
index 48a538566..000000000
--- a/passivity_test.ipynb
+++ /dev/null
@@ -1,105 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 23,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "-inf\n",
-      "-inf\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "c:\\Users\\Mark\\anaconda3\\lib\\site-packages\\scipy\\signal\\lti_conversion.py:109: RuntimeWarning: invalid value encountered in subtract\n",
-      "  C = num[:, 1:] - outer(num[:, 0], den[1:])\n",
-      "c:\\Users\\Mark\\python-control\\control\\passivity.py:77: RuntimeWarning: invalid value encountered in matmul\n",
-      "  np.hstack((off_diag.T, B.T@P@B-(D.T@Q@D + D.T@S + S.T@D + R)))\n",
-      "c:\\Users\\Mark\\python-control\\control\\passivity.py:77: RuntimeWarning: invalid value encountered in add\n",
-      "  np.hstack((off_diag.T, B.T@P@B-(D.T@Q@D + D.T@S + S.T@D + R)))\n"
-     ]
-    },
-    {
-     "ename": "ArithmeticError",
-     "evalue": "1",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[1;31mArithmeticError\u001b[0m                           Traceback (most recent call last)",
-      "\u001b[1;32mc:\\Users\\Mark\\python-control\\passivity_test.ipynb Cell 1\u001b[0m in \u001b[0;36m<cell line: 19>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=15'>16</a>\u001b[0m sys_ff_nu \u001b[39m=\u001b[39m parallel(sys, \u001b[39m-\u001b[39miff_index)\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=16'>17</a>\u001b[0m sys_fb_rho \u001b[39m=\u001b[39m feedback(sys, ofb_index,  sign\u001b[39m=\u001b[39m\u001b[39m1\u001b[39m)\n\u001b[1;32m---> <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=18'>19</a>\u001b[0m \u001b[39massert\u001b[39;00m(sys_ff_nu\u001b[39m.\u001b[39;49mispassive())\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=19'>20</a>\u001b[0m \u001b[39massert\u001b[39;00m(sys_fb_rho\u001b[39m.\u001b[39mispassive())\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/Mark/python-control/passivity_test.ipynb#ch0000000?line=21'>22</a>\u001b[0m sys_ff_nu \u001b[39m=\u001b[39m parallel(sys, \u001b[39m-\u001b[39miff_index\u001b[39m-\u001b[39m\u001b[39m1e-2\u001b[39m)\n",
-      "File \u001b[1;32mc:\\Users\\Mark\\python-control\\control\\lti.py:208\u001b[0m, in \u001b[0;36mLTI.ispassive\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    205\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mispassive\u001b[39m(\u001b[39mself\u001b[39m):\n\u001b[0;32m    206\u001b[0m     \u001b[39m# importing here prevents circular dependancy\u001b[39;00m\n\u001b[0;32m    207\u001b[0m     \u001b[39mfrom\u001b[39;00m \u001b[39mcontrol\u001b[39;00m\u001b[39m.\u001b[39;00m\u001b[39mpassivity\u001b[39;00m \u001b[39mimport\u001b[39;00m ispassive\n\u001b[1;32m--> 208\u001b[0m     \u001b[39mreturn\u001b[39;00m ispassive(\u001b[39mself\u001b[39;49m)\n",
-      "File \u001b[1;32mc:\\Users\\Mark\\python-control\\control\\passivity.py:219\u001b[0m, in \u001b[0;36mispassive\u001b[1;34m(sys, ofp_index, ifp_index)\u001b[0m\n\u001b[0;32m    206\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mispassive\u001b[39m(sys, ofp_index \u001b[39m=\u001b[39m \u001b[39m0\u001b[39m, ifp_index \u001b[39m=\u001b[39m \u001b[39m0\u001b[39m):\n\u001b[0;32m    207\u001b[0m     \u001b[39m'''Indicates if a linear time invariant (LTI) system is passive.\u001b[39;00m\n\u001b[0;32m    208\u001b[0m \n\u001b[0;32m    209\u001b[0m \u001b[39m    Parameters\u001b[39;00m\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    217\u001b[0m \u001b[39m        The input system is passive.\u001b[39;00m\n\u001b[0;32m    218\u001b[0m \u001b[39m    '''\u001b[39;00m\n\u001b[1;32m--> 219\u001b[0m     \u001b[39mreturn\u001b[39;00m _solve_passivity_LMI(sys, rho \u001b[39m=\u001b[39;49m ofp_index, nu \u001b[39m=\u001b[39;49m ifp_index) \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m\n",
-      "File \u001b[1;32mc:\\Users\\Mark\\python-control\\control\\passivity.py:155\u001b[0m, in \u001b[0;36m_solve_passivity_LMI\u001b[1;34m(sys, rho, nu)\u001b[0m\n\u001b[0;32m    153\u001b[0m \u001b[39m# crunch feasibility solution\u001b[39;00m\n\u001b[0;32m    154\u001b[0m cvx\u001b[39m.\u001b[39msolvers\u001b[39m.\u001b[39moptions[\u001b[39m'\u001b[39m\u001b[39mshow_progress\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m \u001b[39mFalse\u001b[39;00m\n\u001b[1;32m--> 155\u001b[0m sol \u001b[39m=\u001b[39m cvx\u001b[39m.\u001b[39;49msolvers\u001b[39m.\u001b[39;49msdp(c, Gs\u001b[39m=\u001b[39;49mGs, hs\u001b[39m=\u001b[39;49mhs)\n\u001b[0;32m    156\u001b[0m \u001b[39mreturn\u001b[39;00m sol[\u001b[39m\"\u001b[39m\u001b[39mx\u001b[39m\u001b[39m\"\u001b[39m]\n",
-      "File \u001b[1;32mc:\\Users\\Mark\\anaconda3\\lib\\site-packages\\cvxopt\\coneprog.py:4126\u001b[0m, in \u001b[0;36msdp\u001b[1;34m(c, Gl, hl, Gs, hs, A, b, kktsolver, solver, primalstart, dualstart, **kwargs)\u001b[0m\n\u001b[0;32m   4123\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m   4124\u001b[0m     ds \u001b[39m=\u001b[39m \u001b[39mNone\u001b[39;00m\n\u001b[1;32m-> 4126\u001b[0m sol \u001b[39m=\u001b[39m conelp(c, G, h, dims, A \u001b[39m=\u001b[39;49m A, b \u001b[39m=\u001b[39;49m b, primalstart \u001b[39m=\u001b[39;49m ps, dualstart \u001b[39m=\u001b[39;49m ds, kktsolver \u001b[39m=\u001b[39;49m kktsolver, options \u001b[39m=\u001b[39;49m options)\n\u001b[0;32m   4127\u001b[0m \u001b[39mif\u001b[39;00m sol[\u001b[39m'\u001b[39m\u001b[39ms\u001b[39m\u001b[39m'\u001b[39m] \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m   4128\u001b[0m     sol[\u001b[39m'\u001b[39m\u001b[39msl\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m \u001b[39mNone\u001b[39;00m\n",
-      "File \u001b[1;32mc:\\Users\\Mark\\anaconda3\\lib\\site-packages\\cvxopt\\coneprog.py:1033\u001b[0m, in \u001b[0;36mconelp\u001b[1;34m(c, G, h, dims, A, b, primalstart, dualstart, kktsolver, xnewcopy, xdot, xaxpy, xscal, ynewcopy, ydot, yaxpy, yscal, **kwargs)\u001b[0m\n\u001b[0;32m   1026\u001b[0m \u001b[39m# Compute initial scaling W:\u001b[39;00m\n\u001b[0;32m   1027\u001b[0m \u001b[39m#\u001b[39;00m\n\u001b[0;32m   1028\u001b[0m \u001b[39m#     W * z = W^{-T} * s = lambda\u001b[39;00m\n\u001b[0;32m   1029\u001b[0m \u001b[39m#     dg * tau = 1/dg * kappa = lambdag.\u001b[39;00m\n\u001b[0;32m   1031\u001b[0m \u001b[39mif\u001b[39;00m iters \u001b[39m==\u001b[39m \u001b[39m0\u001b[39m:\n\u001b[1;32m-> 1033\u001b[0m     W \u001b[39m=\u001b[39m misc\u001b[39m.\u001b[39;49mcompute_scaling(s, z, lmbda, dims, mnl \u001b[39m=\u001b[39;49m \u001b[39m0\u001b[39;49m)\n\u001b[0;32m   1035\u001b[0m     \u001b[39m#     dg = sqrt( kappa / tau )\u001b[39;00m\n\u001b[0;32m   1036\u001b[0m     \u001b[39m#     dgi = sqrt( tau / kappa )\u001b[39;00m\n\u001b[0;32m   1037\u001b[0m     \u001b[39m#     lambda_g = sqrt( tau * kappa )\u001b[39;00m\n\u001b[0;32m   1038\u001b[0m     \u001b[39m#\u001b[39;00m\n\u001b[0;32m   1039\u001b[0m     \u001b[39m# lambda_g is stored in the last position of lmbda.\u001b[39;00m\n\u001b[0;32m   1041\u001b[0m     dg \u001b[39m=\u001b[39m math\u001b[39m.\u001b[39msqrt( kappa \u001b[39m/\u001b[39m tau )\n",
-      "File \u001b[1;32mc:\\Users\\Mark\\anaconda3\\lib\\site-packages\\cvxopt\\misc.py:387\u001b[0m, in \u001b[0;36mcompute_scaling\u001b[1;34m(s, z, lmbda, dims, mnl)\u001b[0m\n\u001b[0;32m    385\u001b[0m \u001b[39m# Factor sk = Ls*Ls'; store Ls in ds[inds[k]:inds[k+1]].\u001b[39;00m\n\u001b[0;32m    386\u001b[0m blas\u001b[39m.\u001b[39mcopy(s, Ls, offsetx \u001b[39m=\u001b[39m ind2, n \u001b[39m=\u001b[39m m\u001b[39m*\u001b[39m\u001b[39m*\u001b[39m\u001b[39m2\u001b[39m) \n\u001b[1;32m--> 387\u001b[0m lapack\u001b[39m.\u001b[39;49mpotrf(Ls, n \u001b[39m=\u001b[39;49m m, ldA \u001b[39m=\u001b[39;49m m)\n\u001b[0;32m    389\u001b[0m \u001b[39m# Factor zs[k] = Lz*Lz'; store Lz in dz[inds[k]:inds[k+1]].\u001b[39;00m\n\u001b[0;32m    390\u001b[0m blas\u001b[39m.\u001b[39mcopy(z, Lz, offsetx \u001b[39m=\u001b[39m ind2, n \u001b[39m=\u001b[39m m\u001b[39m*\u001b[39m\u001b[39m*\u001b[39m\u001b[39m2\u001b[39m) \n",
-      "\u001b[1;31mArithmeticError\u001b[0m: 1"
-     ]
-    }
-   ],
-   "source": [
-    "import matplotlib.pyplot \n",
-    "import numpy\n",
-    "import math\n",
-    "import sys as system\n",
-    "from control import passivity, ss, tf, parallel, feedback, sample_system\n",
-    "from scipy import optimize\n",
-    "\n",
-    "sys = tf([1], [1, -2])\n",
-    "print(sys)\n",
-    "sys = sample_system(sys, Ts=1)\n",
-    "iff_index = passivity.get_input_ff_index(sys)\n",
-    "ofb_index = passivity.get_output_fb_index(sys)\n",
-    "\n",
-    "print(iff_index)\n",
-    "print(ofb_index)\n",
-    "\n",
-    "sys_ff_nu = parallel(sys, -iff_index)\n",
-    "sys_fb_rho = feedback(sys, ofb_index,  sign=1)\n",
-    "\n",
-    "assert(sys_ff_nu.ispassive())\n",
-    "assert(sys_fb_rho.ispassive())\n",
-    "\n",
-    "sys_ff_nu = parallel(sys, -iff_index-1e-2)\n",
-    "sys_fb_rho = feedback(sys, ofb_index+1e-2, sign=1)\n",
-    "\n",
-    "assert(not sys_ff_nu.ispassive())\n",
-    "assert(not sys_fb_rho.ispassive())\n",
-    "\n"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3.9.12 ('base')",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.9.12"
-  },
-  "orig_nbformat": 4,
-  "vscode": {
-   "interpreter": {
-    "hash": "ccecdc02aa416300a64487fe6971a8398ebaf053d2b2768f769ac21dbbf13f86"
-   }
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}

From 2f887786d529e55a1bda74a9759bc6702cadec9e Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Sat, 30 Jul 2022 22:56:50 +0200
Subject: [PATCH 36/36] remove some lint, enforce numpydoc

---
 control/passivity.py | 169 +++++++++++++++++++++++--------------------
 1 file changed, 89 insertions(+), 80 deletions(-)

diff --git a/control/passivity.py b/control/passivity.py
index ee757c4f5..2ec1a7683 100644
--- a/control/passivity.py
+++ b/control/passivity.py
@@ -1,7 +1,9 @@
-'''
-Author: Mark Yeatman  
+"""
+Functions for passive control.
+
+Author: Mark Yeatman
 Date: July 17, 2022
-'''
+"""
 
 import numpy as np
 from control import statesp
@@ -9,53 +11,57 @@
 
 try:
     import cvxopt as cvx
-except ImportError as e:
+except ImportError:
     cvx = None
 
 eps = np.nextafter(0, 1)
 
-__all__ = ["get_output_fb_index", "get_input_ff_index", 
-    "ispassive", "solve_passivity_LMI"]
+__all__ = ["get_output_fb_index", "get_input_ff_index",  "ispassive",
+           "solve_passivity_LMI"]
 
 
 def solve_passivity_LMI(sys, rho=None, nu=None):
-    '''Computes passivity indices and/or solves feasiblity via a linear matrix inequality (LMI).
+    """Compute passivity indices and/or solves feasiblity via a LMI.
 
-    Constructs an LMI such that if a solution exists and the last element of the 
-    solution is positive, the system is passive. Inputs of None for rho or nu indicates 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 is either the 
-    output or input passivity index, for rho=None and nu=None respectively.
+    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 `solution` is either the output or input
+    passivity index, for `rho` = None and  `nu` = None respectively.
 
-    The sources for the algorithm are: 
+    The sources for the algorithm are:
 
     McCourt, Michael J., and Panos J. Antsaklis
-        "Demonstrating passivity and dissipativity using computational methods." 
+        "Demonstrating passivity and dissipativity using computational
+        methods."
 
     Nicholas Kottenstette and Panos J. Antsaklis
-        "Relationships Between Positive Real, Passive Dissipative, & Positive Systems" 
-        equation 36.
+        "Relationships Between Positive Real, Passive Dissipative, & Positive
+        Systems", equation 36.
 
     Parameters
     ----------
-    sys: An LTI system
-        System to be checked.
-    rho: Float or None
+    sys : LTI
+        System to be checked
+    rho : float or None
         Output feedback passivity index
-    nu: Float or None
+    nu : float or None
         Input feedforward passivity index
-        
+
     Returns
     -------
-    numpy array: 
+    solution : ndarray
         The LMI solution
-    '''
+    """
     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.")
+            "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.")
@@ -93,11 +99,12 @@ def make_LMI_matrix(P, rho, nu, one):
             ))
 
     def make_P_basis_matrices(n, rho, nu):
-        '''Makes list of matrix constraints for passivity LMI.
+        """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 parametrized symbolic matrix
-        '''
+        Utility function to make basis matrices for a LMI from a
+        symmetric matrix P of size n by n representing a parametrized symbolic
+        matrix
+        """
         matrix_list = []
         for i in range(0, n):
             for j in range(0, n):
@@ -115,11 +122,11 @@ def make_P_basis_matrices(n, rho, nu):
 
 
     def P_pos_def_constraint(n):
-        '''Makes a list of matrix constraints for P >= 0.
+        """Make a list of matrix constraints for P >= 0.
 
-        Utility function to make basis matrices for a LMI that ensures parametrized symbolic matrix 
-        of size n by n is positive definite.
-        '''
+        Utility function to make basis matrices for a LMI that ensures
+        parametrized symbolic matrix of size n by n is positive definite
+        """
         matrix_list = []
         for i in range(0, n):
             for j in range(0, n):
@@ -134,7 +141,7 @@ def P_pos_def_constraint(n):
 
     n = sys.nstates
 
-    # coefficents for passivity indices and feasibility matrix 
+    # coefficents 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
@@ -145,7 +152,7 @@ def P_pos_def_constraint(n):
         sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, eps, 1.0)
     elif nu is not None:
         sys_constants = -make_LMI_matrix(np.zeros_like(A), eps, nu, 1.0)
-    
+
     sys_coefficents = np.vstack(sys_matrix_list).T
 
     # LMI to ensure P is positive definite
@@ -172,21 +179,21 @@ def P_pos_def_constraint(n):
 
 
 def get_output_fb_index(sys):
-    '''Returns the output feedback passivity (OFP) index for the input system. 
-    
-    The OFP is the largest gain that can be placed in positive feedback 
+    """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: An LTI system
-        System to be checked.
+    sys : LTI
+        System to be checked
 
     Returns
     -------
-    float: 
-        The OFP index 
-    '''
+    float
+        The OFP index
+    """
     sol = solve_passivity_LMI(sys, nu=eps)
     if sol is None:
         raise RuntimeError("LMI passivity problem is infeasible")
@@ -195,22 +202,22 @@ def get_output_fb_index(sys):
 
 
 def get_input_ff_index(sys):
-    '''Returns the input feedforward passivity (IFP) index for the input system. 
-    
-    The IFP is the largest gain that can be placed in negative parallel interconnection 
-    with a system such that the new interconnected system is passive.
+    """Return the input feedforward passivity (IFP) index for the system.
+
+    The 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: An LTI system
+    sys : LTI
         System to be checked.
 
     Returns
     -------
-    float: 
-        The IFP index 
-    '''
-
+    float
+        The IFP index
+    """
     sol = solve_passivity_LMI(sys, rho=eps)
     if sol is None:
         raise RuntimeError("LMI passivity problem is infeasible")
@@ -219,65 +226,67 @@ def get_input_ff_index(sys):
 
 
 def get_relative_index(sys):
-    """Returns the relative passivity index for the system.
+    """Return the relative passivity index for the system.
 
     (not implemented yet)
     """
-    raise NotImplemented("Relative passivity index not implemented")
+    raise NotImplementedError("Relative passivity index not implemented")
 
 
 def get_combined_io_index(sys):
-    """Returns the combined I/O passivity index for the system.
+    """Return the combined I/O passivity index for the system.
 
     (not implemented yet)
     """
-    raise NotImplemented("Combined I/O passivity index not implemented")
+    raise NotImplementedError("Combined I/O passivity index not implemented")
 
 
 def get_directional_index(sys):
-    """Returns the directional passivity index for the system.
+    """Return the directional passivity index for the system.
 
     (not implemented yet)
     """
-    raise NotImplemented("Directional passivity index not implemented")
+    raise NotImplementedError("Directional passivity index not implemented")
 
-    
-def ispassive(sys, ofp_index = 0, ifp_index = 0):
+
+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: An LTI system
-        System to be checked.
-    ofp_index: float
-        Output feedback passivity index.
-    ifp_index: float
-        Input feedforward passivity index.
-        
+    sys : LTI
+        System to be checked
+    ofp_index : float
+        Output feedback passivity index
+    ifp_index : float
+        Input feedforward passivity index
+
     Returns
     -------
-    bool: 
+    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`, `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).
-    
+
+    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." 
+          "Demonstrating passivity and dissipativity using computational
+          methods."
     """
-    return solve_passivity_LMI(sys, rho = ofp_index, nu = ifp_index) is not None
+    return solve_passivity_LMI(sys, rho=ofp_index, nu=ifp_index) is not None