Merge pull request #814 from Mark-Yeatman/main

murrayrm · web-flow · commit 2d0d513b02e5 · 2022-12-16T21:31:18.000-08:00
Removed epsilon perturbation value in solve_passivity_LMI. Fix associated unit test.
diff --git a/control/passivity.py b/control/passivity.py
@@ -14,7 +14,6 @@
 except ImportError:
     cvx = None
 
-eps = np.nextafter(0, 1)
 
 __all__ = ["get_output_fb_index", "get_input_ff_index",  "ispassive",
            "solve_passivity_LMI"]
@@ -28,8 +27,8 @@ def solve_passivity_LMI(sys, rho=None, nu=None):
     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.
+    inequality.) The last element of the output `solution` is either the output or input
+    passivity index, for `rho` = None and `nu` = None respectively.
 
     The sources for the algorithm are:
 
@@ -74,11 +73,8 @@ def solve_passivity_LMI(sys, rho=None, nu=None):
     D = sys.D
 
     # account for strictly proper systems
-    [n, m] = D.shape
-    D = D + eps * np.eye(n, m)
-
+    [_, m] = D.shape
     [n, _] = A.shape
-    A = A - eps*np.eye(n)
 
     def make_LMI_matrix(P, rho, nu, one):
         q = sys.noutputs
@@ -113,7 +109,7 @@ def make_P_basis_matrices(n, rho, nu):
                     P[i, j] = 1
                     P[j, i] = 1
                     matrix_list.append(make_LMI_matrix(P, 0, 0, 0).flatten())
-        zeros = eps*np.eye(n)
+        zeros = 0.0*np.eye(n)
         if rho is None:
             matrix_list.append(make_LMI_matrix(zeros, 1, 0, 0).flatten())
         elif nu is None:
@@ -149,9 +145,9 @@ def P_pos_def_constraint(n):
     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)
+        sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, 0.0, 1.0)
     elif nu is not None:
-        sys_constants = -make_LMI_matrix(np.zeros_like(A), eps, nu, 1.0)
+        sys_constants = -make_LMI_matrix(np.zeros_like(A), 0.0, nu, 1.0)
 
     sys_coefficents = np.vstack(sys_matrix_list).T
 
@@ -174,8 +170,15 @@ def P_pos_def_constraint(n):
 
     # crunch feasibility solution
     cvx.solvers.options['show_progress'] = False
-    sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
-    return sol["x"]
+    try:
+        sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
+        return sol["x"]
+
+    except ZeroDivisionError as e:
+        raise ValueError("The system is probably ill conditioned. "
+                         "Consider perturbing the system matrices by a small amount."
+        ) from e
+
 
 
 def get_output_fb_index(sys):
@@ -194,7 +197,7 @@ def get_output_fb_index(sys):
     float
         The OFP index
     """
-    sol = solve_passivity_LMI(sys, nu=eps)
+    sol = solve_passivity_LMI(sys, nu=0.0)
     if sol is None:
         raise RuntimeError("LMI passivity problem is infeasible")
     else:
@@ -218,7 +221,7 @@ def get_input_ff_index(sys):
     float
         The IFP index
     """
-    sol = solve_passivity_LMI(sys, rho=eps)
+    sol = solve_passivity_LMI(sys, rho=0.0)
     if sol is None:
         raise RuntimeError("LMI passivity problem is infeasible")
     else:
diff --git a/control/tests/passivity_test.py b/control/tests/passivity_test.py
@@ -99,16 +99,14 @@ def test_system_dimension():
 
 
 @pytest.mark.parametrize(
-    "systemmatrices, expected",
-    [((A, B, C, D*0.0), True),
+    "system_matrices, expected",
+    [((A, B, C, D*0), True),
      ((A_d, B, C, D), True),
-     pytest.param((A*1e12, B, C, D*0), True,
-                  marks=pytest.mark.xfail(reason="gh-761")),
      ((A, B*0, C*0, D), True),
      ((A*0, B, C, D), True),
      ((A*0, B*0, C*0, D*0), True)])
-def test_ispassive_edge_cases(systemmatrices, expected):
-    sys = ss(*systemmatrices)
+def test_ispassive_edge_cases(system_matrices, expected):
+    sys = ss(*system_matrices)
     assert passivity.ispassive(sys) == expected
 
 
