remove _dare_slycot and make slycot/scipy structure moer uniform

murrayrm · murrayrm · commit f1f5375978e4 · 2021-12-26T16:56:45.000-08:00
diff --git a/control/mateqn.py b/control/mateqn.py
@@ -3,7 +3,7 @@
 # Implementation of the functions lyap, dlyap, care and dare
 # for solution of Lyapunov and Riccati equations.
 #
-# Author: Bjorn Olofsson
+# Original author: Bjorn Olofsson
 
 # Copyright (c) 2011, All rights reserved.
 
@@ -162,6 +162,7 @@ def lyap(A, Q, C=None, E=None, method=None):
         _check_shape("Q", Q, n, n, square=True, symmetric=True)
 
         if method == 'scipy':
+            # Solve the Lyapunov equation using SciPy
             return sp.linalg.solve_continuous_lyapunov(A, -Q)
 
         # Solve the Lyapunov equation by calling Slycot function sb03md
@@ -177,6 +178,7 @@ def lyap(A, Q, C=None, E=None, method=None):
         _check_shape("C", C, n, m)
 
         if method == 'scipy':
+            # Solve the Sylvester equation using SciPy
             return sp.linalg.solve_sylvester(A, Q, -C)
 
         # Solve the Sylvester equation by calling the Slycot function sb04md
@@ -293,6 +295,7 @@ def dlyap(A, Q, C=None, E=None, method=None):
         _check_shape("Q", Q, n, n, square=True, symmetric=True)
 
         if method == 'scipy':
+            # Solve the Lyapunov equation using SciPy
             return sp.linalg.solve_discrete_lyapunov(A, Q)
 
         # Solve the Lyapunov equation by calling the Slycot function sb03md
@@ -396,24 +399,6 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None,
     # Decide what method to use
     method = _slycot_or_scipy(method)
 
-    if method == 'slycot':
-        # Make sure we can import required slycot routines
-        try:
-            from slycot import sb02md
-        except ImportError:
-            raise ControlSlycot("Can't find slycot module 'sb02md'")
-
-        try:
-            from slycot import sb02mt
-        except ImportError:
-            raise ControlSlycot("Can't find slycot module 'sb02mt'")
-
-        # Make sure we can find the required slycot routine
-        try:
-            from slycot import sg02ad
-        except ImportError:
-            raise ControlSlycot("Can't find slycot module 'sg02ad'")
-
     # Reshape input arrays
     A = np.array(A, ndmin=2)
     B = np.array(B, ndmin=2)
@@ -447,9 +432,16 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None,
             E, _ = np.linalg.eig(A - B @ K)
             return _ssmatrix(X), E, _ssmatrix(K)
 
-        # Create back-up of arrays needed for later computations
-        R_ba = copy(R)
-        B_ba = copy(B)
+        # Make sure we can import required slycot routines
+        try:
+            from slycot import sb02md
+        except ImportError:
+            raise ControlSlycot("Can't find slycot module 'sb02md'")
+
+        try:
+            from slycot import sb02mt
+        except ImportError:
+            raise ControlSlycot("Can't find slycot module 'sb02mt'")
 
         # Solve the standard algebraic Riccati equation by calling Slycot
         # functions sb02mt and sb02md
@@ -459,7 +451,7 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None,
         X, rcond, w, S_o, U, A_inv = sb02md(n, A, G, Q, 'C', sort=sort)
 
         # Calculate the gain matrix G
-        G = solve(R_ba, B_ba.T) @ X
+        G = solve(R, B.T) @ X
 
         # Return the solution X, the closed-loop eigenvalues L and
         # the gain matrix G
@@ -486,11 +478,11 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None,
             eigs, _ = sp.linalg.eig(A - B @ K, E)
             return _ssmatrix(X), eigs, _ssmatrix(K)
 
-        # Create back-up of arrays needed for later computations
-        R_b = copy(R)
-        B_b = copy(B)
-        E_b = copy(E)
-        S_b = copy(S)
+        # Make sure we can find the required slycot routine
+        try:
+            from slycot import sg02ad
+        except ImportError:
+            raise ControlSlycot("Can't find slycot module 'sg02ad'")
 
         # Solve the generalized algebraic Riccati equation by calling the
         # Slycot function sg02ad
@@ -505,7 +497,7 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None,
         L = np.array([(alfar[i] + alfai[i]*1j) / beta[i] for i in range(n)])
 
         # Calculate the gain matrix G
-        G = solve(R_b, B_b.T @ X @ E_b + S_b.T)
+        G = solve(R, B.T @ X @ E + S.T)
 
         # Return the solution X, the closed-loop eigenvalues L and
         # the gain matrix G
@@ -589,14 +581,11 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True, method=None,
         _check_shape(S_s, S, n, m)
 
     # Figure out how to solve the problem
-    if method == 'scipy' and not stabilizing:
-        raise ControlArgument(
-            "method='scipy' not valid when stabilizing is not True")
-
-    elif method == 'slycot':
-        return _dare_slycot(A, B, Q, R, S, E, stabilizing)
+    if method == 'scipy':
+        if not stabilizing:
+            raise ControlArgument(
+                "method='scipy' not valid when stabilizing is not True")
 
-    else:
         X = sp.linalg.solve_discrete_are(A, B, Q, R, e=E, s=S)
         if S is None:
             G = solve(B.T @ X @ B + R, B.T @ X @ A)
@@ -609,28 +598,12 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True, method=None,
 
         return _ssmatrix(X), L, _ssmatrix(G)
 
-
-def _dare_slycot(A, B, Q, R, S=None, E=None, stabilizing=True):
-    # Make sure we can import required slycot routines
-    try:
-        from slycot import sb02md
-    except ImportError:
-        raise ControlSlycot("Can't find slycot module 'sb02md'")
-
-    try:
-        from slycot import sb02mt
-    except ImportError:
-        raise ControlSlycot("Can't find slycot module 'sb02mt'")
-
+    # Make sure we can import required slycot routine
     try:
         from slycot import sg02ad
     except ImportError:
         raise ControlSlycot("Can't find slycot module 'sg02ad'")
 
-    # Determine main dimensions
-    n = A.shape[0]
-    m = B.shape[1]
-
     # Initialize optional matrices
     S = np.zeros((n, m)) if S is None else np.array(S, ndmin=2)
     E = np.eye(A.shape[0]) if E is None else np.array(E, ndmin=2)
