remove np.matrix dependencies from module code

murrayrm · murrayrm · commit f010542d86d3 · 2019-04-15T23:26:53.000-07:00
diff --git a/control/bdalg.py b/control/bdalg.py
@@ -334,7 +334,7 @@ def connect(sys, Q, inputv, outputv):
                 K[inp,outp-1] = 1.
             elif outp < 0 and -outp >= -sys.outputs:
                 K[inp,-outp-1] = -1.
-    sys = sys.feedback(sp.matrix(K), sign=1)
+    sys = sys.feedback(sp.array(K), sign=1)
 
     # now trim
     Ytrim = sp.zeros( (len(outputv), sys.outputs) )
@@ -343,4 +343,4 @@ def connect(sys, Q, inputv, outputv):
         Utrim[u-1,i] = 1.
     for i,y in enumerate(outputv):
         Ytrim[i,y-1] = 1.
-    return sp.matrix(Ytrim)*sys*sp.matrix(Utrim)
+    return sp.array(Ytrim) * sys * sp.array(Utrim)
diff --git a/control/frdata.py b/control/frdata.py
@@ -51,7 +51,7 @@
 # External function declarations
 import numpy as np
 from numpy import angle, array, empty, ones, \
-    real, imag, matrix, absolute, eye, linalg, where, dot
+    real, imag, absolute, eye, linalg, where, dot
 from scipy.interpolate import splprep, splev
 from .lti import LTI
 
@@ -80,6 +80,10 @@ class FRD(LTI):
 
     """
 
+    # Allow NDarray * StateSpace to give StateSpace._rmul_() priority
+    # https://docs.scipy.org/doc/numpy/reference/arrays.classes.html#numpy.class.__array_priority__
+    __array_priority__ = 11     # override ndarray and matrix types
+
     epsw = 1e-8
 
     def __init__(self, *args, **kwargs):
@@ -433,12 +437,13 @@ def feedback(self, other=1, sign=-1):
         # TODO: vectorize this
         # TODO: handle omega re-mapping
         for k, w in enumerate(other.omega):
-            fresp[:, :, k] = self.fresp[:, :, k].view(type=matrix)* \
+            fresp[:, :, k] = np.dot(
+                self.fresp[:, :, k],
                 linalg.solve(
-                eye(self.inputs) +
-                other.fresp[:, :, k].view(type=matrix) *
-                self.fresp[:, :, k].view(type=matrix),
-                eye(self.inputs))
+                    eye(self.inputs) + np.dot(other.fresp[:, :, k], 
+                                              self.fresp[:, :, k]),
+                    eye(self.inputs))
+            )
 
         return FRD(fresp, other.omega, smooth=(self.ifunc is not None))
 
diff --git a/control/mateqn.py b/control/mateqn.py
@@ -41,7 +41,7 @@
 Author: Bjorn Olofsson
 """
 
-from scipy import shape, size, asarray, asmatrix, copy, zeros, eye, dot
+from scipy import shape, size, array, asarray, copy, zeros, eye, dot
 from scipy.linalg import eigvals, solve_discrete_are, solve
 from .exception import ControlSlycot, ControlArgument
 
@@ -703,8 +703,8 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True):
     if S is not None or E is not None or not stabilizing:
         return dare_old(A, B, Q, R, S, E, stabilizing)
     else:
-        Rmat = asmatrix(R)
-        Qmat = asmatrix(Q)
+        Rmat = array(R, ndmin=2)
+        Qmat = array(Q, ndmin=2)
         X = solve_discrete_are(A, B, Qmat, Rmat)
         G = solve(B.T.dot(X).dot(B) + Rmat, B.T.dot(X).dot(A))
         L = eigvals(A - B.dot(G))
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -47,7 +47,7 @@
 import numpy as np
 from .exception import ControlSlycot
 from .lti import isdtime, isctime
-from .statesp import StateSpace
+from .statesp import StateSpace, ssmatrix
 from .statefbk import gram
 
 __all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
@@ -96,7 +96,7 @@ def hsvd(sys):
     w, v = np.linalg.eig(WoWc)
 
     hsv = np.sqrt(w)
-    hsv = np.matrix(hsv)
+    hsv = np.array(hsv, ndmin=2)                # was np.matrix(hsv)
     hsv = np.sort(hsv)
     hsv = np.fliplr(hsv)
     # Return the Hankel singular values
@@ -406,16 +406,17 @@ def markov(Y, U, M):
     """
 
     # Convert input parameters to matrices (if they aren't already)
-    Ymat = np.mat(Y)
-    Umat = np.mat(U)
+    Ymat = np.array(Y)
+    Umat = np.array(U)
     n = np.size(U)
 
     # Construct a matrix of control inputs to invert
     UU = Umat
     for i in range(1, M-1):
-        newCol = np.vstack((0, UU[0:n-1,i-2]))
+        #! TODO: second index on UU doesn't seem right; could be neg or pos??
+        newCol = np.vstack((0, np.reshape(UU[0:n-1,i-2], (-1,1))))
         UU = np.hstack((UU, newCol))
-    Ulast = np.vstack((0, UU[0:n-1,M-2]))
+    Ulast = np.vstack((0, np.reshape(UU[0:n-1,M-2], (-1,1))))
     for i in range(n-1,0,-1):
         Ulast[i] = np.sum(Ulast[0:i-1])
     UU = np.hstack((UU, Ulast))
diff --git a/control/statefbk.py b/control/statefbk.py
@@ -238,9 +238,9 @@ def acker(A, B, poles):
         Gains such that A - B K has given eigenvalues
 
     """
-    # Convert the inputs to matrices
-    a = np.mat(A)
-    b = np.mat(B)
+    # Convert the inputs to matrices (arrays)
+    a = np.array(A)
+    b = np.array(B)
 
     # Make sure the system is controllable
     ct = ctrb(A, B)
@@ -252,9 +252,9 @@ def acker(A, B, poles):
 
     # Place the poles using Ackermann's method
     n = np.size(p)
-    pmat = p[n-1]*a**0
+    pmat = p[n-1] * np.linalg.matrix_power(a, 0)
     for i in np.arange(1,n):
-        pmat = pmat + p[n-i-1]*a**i
+        pmat = pmat + np.dot(p[n-i-1], np.linalg.matrix_power(a, i))
     K = np.linalg.solve(ct, pmat)
 
     K = K[-1][:]                # Extract the last row
@@ -388,13 +388,13 @@ def ctrb(A,B):
     """
 
     # Convert input parameters to matrices (if they aren't already)
-    amat = np.mat(A)
-    bmat = np.mat(B)
+    amat = np.array(A)
+    bmat = np.array(B)
     n = np.shape(amat)[0]
     # Construct the controllability matrix
     ctrb = bmat
     for i in range(1, n):
-        ctrb = np.hstack((ctrb, amat**i*bmat))
+        ctrb = np.hstack((ctrb, np.dot(np.linalg.matrix_power(amat, i), bmat)))
     return ctrb
 
 def obsv(A, C):
@@ -417,14 +417,14 @@ def obsv(A, C):
    """
 
     # Convert input parameters to matrices (if they aren't already)
-    amat = np.mat(A)
-    cmat = np.mat(C)
+    amat = np.array(A)
+    cmat = np.array(C)
     n = np.shape(amat)[0]
 
     # Construct the controllability matrix
     obsv = cmat
     for i in range(1, n):
-        obsv = np.vstack((obsv, cmat*amat**i))
+        obsv = np.vstack((obsv, np.dot(cmat, np.linalg.matrix_power(amat, i))))
     return obsv
 
 def gram(sys,type):
diff --git a/control/timeresp.py b/control/timeresp.py
@@ -310,7 +310,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
             #   [ u(dt) ] = exp [  0     0    I ] [  u0   ]
             #   [u1 - u0]       [  0     0    0 ] [u1 - u0]
 
-            M = np.bmat([[A * dt, B * dt, np.zeros((n_states, n_inputs))],
+            M = np.block([[A * dt, B * dt, np.zeros((n_states, n_inputs))],
                          [np.zeros((n_inputs, n_states + n_inputs)),
                           np.identity(n_inputs)],
                          [np.zeros((n_inputs, n_states + 2 * n_inputs))]])
diff --git a/control/xferfcn.py b/control/xferfcn.py
@@ -91,6 +91,10 @@ class TransferFunction(LTI):
     discrete time system with unspecified sampling time.
     """
 
+    # Allow NDarray * StateSpace to give StateSpace._rmul_() priority
+    # https://docs.scipy.org/doc/numpy/reference/arrays.classes.html#numpy.class.__array_priority__
+    __array_priority__ = 11     # override ndarray and matrix types
+
     def __init__(self, *args):
         """TransferFunction(num, den[, dt])
 
