python-control · slivingston · Dec 31, 2016 · Jun 13, 2016 · Jun 13, 2016 · Jun 13, 2016
diff --git a/ChangeLog b/ChangeLog
@@ -148,7 +148,7 @@
 2012-11-03  Richard Murray  <murray@altura.local>
 
 	* src/rlocus.py (_RLSortRoots): convert output of range() to
-	explicit list for python 3 compatability
+	explicit list for python 3 compatibility
 
 	* tests/modelsimp_test.py, tests/slycot_convert_test.py,
 	tests/mateqn_test.py, tests/statefbk_test.py: updated test suites to
@@ -604,7 +604,7 @@
 	initial_response, impulse_response and step_response.
 
 	* src/rlocus.py: changed RootLocus to root_locus for better
-	compatability with PEP 8.  Also updated unit tests and examples.
+	compatibility with PEP 8.  Also updated unit tests and examples.
 
 2011-07-25  Richard Murray  <murray@malabar.local>
 
@@ -994,7 +994,7 @@
 2010-09-02  Richard Murray  <murray@sumatra.local>
 
 	* src/statefbk.py (place): Use np.size() instead of len() for
-	finding length of placed_eigs for better compatability with
+	finding length of placed_eigs for better compatibility with
 	different python versions [courtesy of Roberto Bucher]
 
 	* src/delay.py (pade): New file for delay-based computations +

diff --git a/control/canonical.py b/control/canonical.py
@@ -7,7 +7,7 @@
 from .statefbk import ctrb, obsv
 
 from numpy import zeros, shape, poly
-from numpy.linalg import inv
+from numpy.linalg import solve, matrix_rank
 
 __all__ = ['canonical_form', 'reachable_form', 'observable_form']
 
@@ -68,23 +68,29 @@ def reachable_form(xsys):
 
     # Generate the system matrices for the desired canonical form
     zsys.B = zeros(shape(xsys.B))
-    zsys.B[0, 0] = 1
+    zsys.B[0, 0] = 1.0
     zsys.A = zeros(shape(xsys.A))
     Apoly = poly(xsys.A)                # characteristic polynomial
     for i in range(0, xsys.states):
         zsys.A[0, i] = -Apoly[i+1] / Apoly[0]
         if (i+1 < xsys.states):
-            zsys.A[i+1, i] = 1
+            zsys.A[i+1, i] = 1.0
 
     # Compute the reachability matrices for each set of states
     Wrx = ctrb(xsys.A, xsys.B)
     Wrz = ctrb(zsys.A, zsys.B)
 
+    if matrix_rank(Wrx) != xsys.states:
+        raise ValueError("System not controllable to working precision.")
+
     # Transformation from one form to another
-    Tzx = Wrz * inv(Wrx)
+    Tzx = solve(Wrx.T, Wrz.T).T  # matrix right division, Tzx = Wrz * inv(Wrx)
+
+    if matrix_rank(Tzx) != xsys.states:
+        raise ValueError("Transformation matrix singular to working precision.")
 
     # Finally, compute the output matrix
-    zsys.C = xsys.C * inv(Tzx)
+    zsys.C = solve(Tzx.T, xsys.C.T).T  # matrix right division, zsys.C = xsys.C * inv(Tzx)
 
     return zsys, Tzx
 

diff --git a/control/margins.py b/control/margins.py
@@ -8,7 +8,7 @@
 margin.phase_crossover_frequencies
 """
 
-# Python 3 compatability (needs to go here)
+# Python 3 compatibility (needs to go here)
 from __future__ import print_function
 
 """Copyright (c) 2011 by California Institute of Technology

diff --git a/control/mateqn.py b/control/mateqn.py
@@ -5,7 +5,7 @@
 Implementation of the functions lyap, dlyap, care and dare
 for solution of Lyapunov and Riccati equations. """
 
-# Python 3 compatability (needs to go here)
+# Python 3 compatibility (needs to go here)
 from __future__ import print_function
 
 """Copyright (c) 2011, All rights reserved.
@@ -41,9 +41,8 @@
 Author: Bjorn Olofsson
 """
 
-from numpy.linalg import inv
 from scipy import shape, size, asarray, asmatrix, copy, zeros, eye, dot
-from scipy.linalg import eigvals, solve_discrete_are
+from scipy.linalg import eigvals, solve_discrete_are, solve
 from .exception import ControlSlycot, ControlArgument
 
 __all__ = ['lyap', 'dlyap', 'dare', 'care']
@@ -557,9 +556,9 @@ def care(A,B,Q,R=None,S=None,E=None):
 
         # Calculate the gain matrix G
         if size(R_b) == 1:
-            G = dot(dot(1/(R_ba),asarray(B_ba).T) , X)
+            G = dot(dot(1/(R_ba), asarray(B_ba).T), X)
         else:
-            G = dot(dot(inv(R_ba),asarray(B_ba).T) , X)
+            G = dot(solve(R_ba, asarray(B_ba).T), X)
 
         # Return the solution X, the closed-loop eigenvalues L and
         # the gain matrix G
@@ -660,9 +659,9 @@ def care(A,B,Q,R=None,S=None,E=None):
 
         # Calculate the gain matrix G
         if size(R_b) == 1:
-            G = dot(1/(R_b),dot(asarray(B_b).T,dot(X,E_b))+asarray(S_b).T)
+            G = dot(1/(R_b), dot(asarray(B_b).T, dot(X,E_b)) + asarray(S_b).T)
         else:
-            G = dot(inv(R_b),dot(asarray(B_b).T,dot(X,E_b))+asarray(S_b).T)
+            G = solve(R_b, dot(asarray(B_b).T, dot(X, E_b)) + asarray(S_b).T)
 
         # Return the solution X, the closed-loop eigenvalues L and
         # the gain matrix G
@@ -699,7 +698,7 @@ def dare(A,B,Q,R,S=None,E=None):
         Rmat = asmatrix(R)
         Qmat = asmatrix(Q)
         X = solve_discrete_are(A, B, Qmat, Rmat)
-        G = inv(B.T.dot(X).dot(B) + Rmat) * B.T.dot(X).dot(A)
+        G = solve(B.T.dot(X).dot(B) + Rmat, B.T.dot(X).dot(A))
         L = eigvals(A - B.dot(G))
         return X, L, G
 
@@ -825,11 +824,11 @@ def dare_old(A,B,Q,R,S=None,E=None):
 
         # Calculate the gain matrix G
         if size(R_b) == 1:
-            G = dot( 1/(dot(asarray(B_ba).T,dot(X,B_ba))+R_ba) , \
-                dot(asarray(B_ba).T,dot(X,A_ba)) )
+            G = dot(1/(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba), \
+                dot(asarray(B_ba).T, dot(X, A_ba)))
         else:
-            G = dot( inv(dot(asarray(B_ba).T,dot(X,B_ba))+R_ba) , \
-                dot(asarray(B_ba).T,dot(X,A_ba)) )
+            G = solve(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba, \
+                dot(asarray(B_ba).T, dot(X, A_ba)))
 
         # Return the solution X, the closed-loop eigenvalues L and
         # the gain matrix G
@@ -930,11 +929,11 @@ def dare_old(A,B,Q,R,S=None,E=None):
 
         # Calculate the gain matrix G
         if size(R_b) == 1:
-            G = dot( 1/(dot(asarray(B_b).T,dot(X,B_b))+R_b) , \
-                dot(asarray(B_b).T,dot(X,A_b)) + asarray(S_b).T)
+            G = dot(1/(dot(asarray(B_b).T, dot(X,B_b)) + R_b), \
+                dot(asarray(B_b).T, dot(X,A_b)) + asarray(S_b).T)
         else:
-            G = dot( inv(dot(asarray(B_b).T,dot(X,B_b))+R_b) , \
-                dot(asarray(B_b).T,dot(X,A_b)) + asarray(S_b).T)
+            G = solve(dot(asarray(B_b).T, dot(X,B_b)) + R_b, \
+                dot(asarray(B_b).T, dot(X,A_b)) + asarray(S_b).T)
 
         # Return the solution X, the closed-loop eigenvalues L and
         # the gain matrix G

diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -40,7 +40,7 @@
 #
 # $Id$
 
-# Python 3 compatability
+# Python 3 compatibility
 from __future__ import print_function
 
 # External packages and modules
@@ -149,16 +149,12 @@ def modred(sys, ELIM, method='matchdc'):
 
 
     #Check system is stable
-    D,V = np.linalg.eig(sys.A)
-    for e in D:
-        if e.real >= 0:
-            raise ValueError("Oops, the system is unstable!")
+    if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
+        raise ValueError("Oops, the system is unstable!")
+
     ELIM = np.sort(ELIM)
-    NELIM = []
     # Create list of elements not to eliminate (NELIM)
-    for i in range(0,len(sys.A)):
-        if i not in ELIM:
-            NELIM.append(i)
+    NELIM = [i for i in range(len(sys.A)) if i not in ELIM]
     # A1 is a matrix of all columns of sys.A not to eliminate
     A1 = sys.A[:,NELIM[0]]
     for i in NELIM[1:]:
@@ -177,14 +173,26 @@ def modred(sys, ELIM, method='matchdc'):
     B1 = sys.B[NELIM,:]
     B2 = sys.B[ELIM,:]
 
-    A22I = np.linalg.inv(A22)
-
     if method=='matchdc':
         # if matchdc, residualize
-        Ar = A11 - A12*A22.I*A21
-        Br = B1 - A12*A22.I*B2
-        Cr = C1 - C2*A22.I*A21
-        Dr = sys.D - C2*A22.I*B2
+
+        # Check if the matrix A22 is invertible
+        if np.linalg.matrix_rank(A22) != len(ELIM):
+            raise ValueError("Matrix A22 is singular to working precision.")
+
+        # Now precompute A22\A21 and A22\B2 (A22I = inv(A22))
+        # We can solve two linear systems in one pass, since the
+        # coefficients matrix A22 is the same. Thus, we perform the LU
+        # decomposition (cubic runtime complexity) of A22 only once!
+        # The remaining back substitutions are only quadratic in runtime.
+        A22I_A21_B2 = np.linalg.solve(A22, np.concatenate((A21, B2), axis=1))
+        A22I_A21 = A22I_A21_B2[:, :A21.shape[1]]
+        A22I_B2 = A22I_A21_B2[:, A21.shape[1]:]
+
+        Ar = A11 - A12*A22I_A21
+        Br = B1 - A12*A22I_B2
+        Cr = C1 - C2*A22I_A21
+        Dr = sys.D - C2*A22I_B2
     elif method=='truncate':
         # if truncate, simply discard state x2
         Ar = A11
@@ -243,12 +251,8 @@ def balred(sys, orders, method='truncate'):
     dico = 'C'
 
     #Check system is stable
-    D,V = np.linalg.eig(sys.A)
-    # print D.shape
-    # print D
-    for e in D:
-        if e.real >= 0:
-            raise ValueError("Oops, the system is unstable!")
+    if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
+        raise ValueError("Oops, the system is unstable!")
 
     if method=='matchdc':
         raise ValueError ("MatchDC not yet supported!")
@@ -373,8 +377,6 @@ def markov(Y, U, M):
     UU = np.hstack((UU, Ulast))
 
     # Invert and solve for Markov parameters
-    H = UU.I
-    H = np.dot(H, Y)
+    H = np.linalg.lstsq(UU, Y)[0]
 
     return H
-
diff --git a/control/phaseplot.py b/control/phaseplot.py
@@ -34,7 +34,7 @@
 # IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
 
-# Python 3 compatability
+# Python 3 compatibility
 from __future__ import print_function
 
 import numpy as np

diff --git a/control/statefbk.py b/control/statefbk.py
@@ -127,7 +127,7 @@ def acker(A, B, poles):
 
     # Make sure the system is controllable
     ct = ctrb(A, B)
-    if sp.linalg.det(ct) == 0:
+    if np.linalg.matrix_rank(ct) != a.shape[0]:
         raise ValueError("System not reachable; pole placement invalid")
 
     # Compute the desired characteristic polynomial
@@ -138,7 +138,7 @@ def acker(A, B, poles):
     pmat = p[n-1]*a**0
     for i in np.arange(1,n):
         pmat = pmat + p[n-i-1]*a**i
-    K = sp.linalg.inv(ct) * pmat
+    K = np.linalg.solve(ct, pmat)
 
     K = K[-1][:]                # Extract the last row
     return K
@@ -242,7 +242,8 @@ def lqr(*args, **keywords):
     X,rcond,w,S,U,A_inv = sb02md(nstates, A_b, G, Q_b, 'C')
 
     # Now compute the return value
-    K = np.dot(np.linalg.inv(R), (np.dot(B.T, X) + N.T));
+    # We assume that R is positive definite and, hence, invertible
+    K = np.linalg.solve(R, np.dot(B.T, X) + N.T);
     S = X;
     E = w[0:nstates];
 
@@ -354,10 +355,9 @@ def gram(sys,type):
 
     #TODO: Check system is stable, perhaps a utility in ctrlutil.py
         # or a method of the StateSpace class?
-    D,V = np.linalg.eig(sys.A)
-    for e in D:
-        if e.real >= 0:
-            raise ValueError("Oops, the system is unstable!")
+    if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
+        raise ValueError("Oops, the system is unstable!")
+
     if type=='c':
         tra = 'T'
         C = -np.dot(sys.B,sys.B.transpose())
@@ -379,4 +379,3 @@ def gram(sys,type):
     X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
     gram = X
     return gram
-
diff --git a/control/statesp.py b/control/statesp.py
@@ -8,7 +8,7 @@
 
 """
 
-# Python 3 compatability (needs to go here)
+# Python 3 compatibility (needs to go here)
 from __future__ import print_function
 
 """Copyright (c) 2010 by California Institute of Technology
@@ -124,12 +124,7 @@ def __init__(self, *args):
 
         # Here we're going to convert inputs to matrices, if the user gave a
         # non-matrix type.
-        #! TODO: [A, B, C, D] = map(matrix, [A, B, C, D])?
-        matrices = [A, B, C, D]
-        for i in range(len(matrices)):
-            # Convert to matrix first, if necessary.
-            matrices[i] = matrix(matrices[i])
-        [A, B, C, D] = matrices
+        A, B, C, D = map(matrix, [A, B, C, D])
 
         LTI.__init__(self, B.shape[1], C.shape[0], dt)
         self.A = A
@@ -453,11 +448,16 @@ def feedback(self, other=1, sign=-1):
 
         F = eye(self.inputs) - sign * D2 * D1
         if matrix_rank(F) != self.inputs:
-            raise ValueError("I - sign * D2 * D1 is singular.")
+            raise ValueError("I - sign * D2 * D1 is singular to working precision.")
 
         # Precompute F\D2 and F\C2 (E = inv(F))
-        E_D2 = solve(F, D2)
-        E_C2 = solve(F, C2)
+        # We can solve two linear systems in one pass, since the
+        # coefficients matrix F is the same. Thus, we perform the LU
+        # decomposition (cubic runtime complexity) of F only once!
+        # The remaining back substitutions are only quadratic in runtime.
+        E_D2_C2 = solve(F, concatenate((D2, C2), axis=1))
+        E_D2 = E_D2_C2[:, :other.inputs]
+        E_C2 = E_D2_C2[:, other.inputs:]
 
         T1 = eye(self.outputs) + sign * D1 * E_D2
         T2 = eye(self.inputs) + sign * E_D2 * D1