python-control · murrayrm · Apr 16, 2019 · Apr 16, 2019 · Apr 16, 2019 · Apr 16, 2019
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/config.py b/control/config.py
@@ -7,13 +7,20 @@
 # files.  For now, you can just choose between MATLAB and FBS default
 # values.
 
+import warnings
+import numpy as np
+from .statesp import StateSpaceMatrix
+
 # Bode plot defaults
 bode_dB = False                 # Bode plot magnitude units
 bode_deg = True                 # Bode Plot phase units
 bode_Hz = False                 # Bode plot frequency units
 bode_number_of_samples = None   # Bode plot number of samples
 bode_feature_periphery_decade = 1.0  # Bode plot feature periphery in decades
 
+# State space return type (change to StateSpaceMatrix at next major revision)
+ss_return_type = np.matrix
+
 
 def reset_defaults():
     """Reset configuration values to their default values."""
@@ -22,6 +29,7 @@ def reset_defaults():
     global bode_Hz; bode_Hz = False
     global bode_number_of_samples; bode_number_of_samples = None
     global bode_feature_periphery_decade; bode_feature_periphery_decade = 1.0
+    global ss_return_type; ss_return_type = np.matrix   # TODO: update
 
 
 # Set defaults to match MATLAB
@@ -52,3 +60,49 @@ def use_fbs_defaults():
     global bode_deg; bode_deg = True
     global bode_Hz; bode_Hz = False
 
+#
+# State space function return type
+#
+# These functions are used to set the return type for state space functions
+# that return a matrix.  In the original version of python-control, these
+# functions returned a numpy.matrix object.  To handle the eventual
+# deprecation of the numpy.matrix type, the StateSpaceMatrix object type,
+# which implements matrix multiplications and related numpy.matrix operations
+# was created.  To avoid breaking existing code, the return type from
+# functions that used to return an np.matrix object now call the
+# get_ss_return_type function to obtain the return type.  The return type can
+# also be set using the `return_type` keyword, overriding the default state
+# space matrix return type.
+#
+
+
+# Set state space return type
+def set_ss_return_type(subtype, warn=True):
+    global ss_return_type
+    # If return_type is np.matrix, issue a pending deprecation warning
+    if (subtype is np.matrix and warn):
+        warnings.warn("Return type numpy.matrix is soon to be deprecated.",
+                      stacklevel=2)
+    ss_return_type = subtype
+
+
+# Get the state space return type
+def get_ss_return_type(subtype=None, warn=True):
+    global ss_return_type
+    return_type = ss_return_type if subtype is None else subtype
+    # If return_type is np.matrix, issue a pending deprecation warning
+    if (return_type is np.matrix and warn):
+        warnings.warn("Returning numpy.matrix, soon to be deprecated; "
+                      "make sure calling code can handle nparray.",
+                      stacklevel=2)
+    return return_type
+
+
+# Function to turn on/off use of np.matrix type
+def use_numpy_matrix(flag=True, warn=True):
+    if flag and warn:
+        warnings.warn("Return type numpy.matrix is soon to be deprecated.",
+                      stacklevel=2)
+        set_ss_return_type(np.matrix, warn=False)
+    else:
+        set_ss_return_type(StateSpaceMatrix)
diff --git a/control/frdata.py b/control/frdata.py
@@ -52,7 +52,7 @@
 from warnings import warn
 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
 
@@ -81,6 +81,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):
@@ -436,12 +440,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
@@ -49,14 +49,15 @@
 from .lti import isdtime, isctime
 from .statesp import StateSpace
 from .statefbk import gram
+from .config import get_ss_return_type
 
 __all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
 
 # Hankel Singular Value Decomposition
 #   The following returns the Hankel singular values, which are singular values
 #of the matrix formed by multiplying the controllability and observability
 #grammians
-def hsvd(sys):
+def hsvd(sys, return_type=None):
     """Calculate the Hankel singular values.
 
     Parameters
@@ -66,9 +67,12 @@ def hsvd(sys):
 
     Returns
     -------
-    H : Matrix
+    H : matrix
         A list of Hankel singular values
 
+    return_type: nparray subtype, optional (default = numpy.matrix)
+        Set the ndarray subtype for the return value
+
     See Also
     --------
     gram
@@ -96,11 +100,12 @@ 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
-    return hsv
+
+    # Return the Hankel singular values (casting type, if needed)
+    return hsv.view(type=get_ss_return_type(return_type))
 
 def modred(sys, ELIM, method='matchdc'):
     """
@@ -409,16 +414,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))