Dec 31, 2016 · Jun 13, 2016 · Jun 13, 2016 · Jun 13, 2016 · Jul 3, 2016 · Jul 3, 2016
diff --git a/ChangeLog b/ChangeLog
 2012-11-03  Richard Murray  <murray@altura.local>

 * src/rlocus.py (_RLSortRoots): convert output of range() to
 explicit list for python 3compatability
 explicit list for python 3compatibility

 * tests/modelsimp_test.py, tests/slycot_convert_test.py,
 tests/mateqn_test.py, tests/statefbk_test.py: updated test suites to
 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>

 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 bettercompatability with
 finding length of placed_eigs for bettercompatibility 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
 from .statefbk import ctrb, obsv

 from numpy import zeros, shape, poly
 from numpy.linalg importinv
 from numpy.linalg importsolve, matrix_rank

 __all__ = ['canonical_form', 'reachable_form', 'observable_form']


    # 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
 margin.phase_crossover_frequencies
 """

 # Python 3compatability (needs to go here)
 # Python 3compatibility (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
 Implementation of the functions lyap, dlyap, care and dare
 for solution of Lyapunov and Riccati equations. """

 # Python 3compatability (needs to go here)
 # Python 3compatibility (needs to go here)
 from __future__ import print_function

 """Copyright (c) 2011, All rights reserved.
 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']

        # 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

        # 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
        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


        # 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

        # 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
 #
 # $Id$

 # Python 3compatability
 # Python 3compatibility
 from __future__ import print_function

 # External packages and modules


    #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:]:
    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
    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!")
    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
 # IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.

 # Python 3compatability
 # Python 3compatibility
 from __future__ import print_function

 import numpy as np
diff --git a/control/statefbk.py b/control/statefbk.py

    # Make sure the system is controllable
    ct = ctrb(A, B)
    ifsp.linalg.det(ct)== 0:
    ifnp.linalg.matrix_rank(ct)!= a.shape[0]:
        raise ValueError("System not reachable; pole placement invalid")

    # Compute the desired characteristic polynomial
    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
    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];


    #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())
    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

 """

 # Python 3compatability (needs to go here)
 # Python 3compatibility (needs to go here)
 from __future__ import print_function

 """Copyright (c) 2010 by California Institute of Technology

        # 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

        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
Original file line number	Diff line number	Diff line change
Expand Up		@@ -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 3compatability
		explicit list for python 3compatibility

		* tests/modelsimp_test.py, tests/slycot_convert_test.py,
		tests/mateqn_test.py, tests/statefbk_test.py: updated test suites to
Expand DownExpand Up		@@ -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>

Expand DownExpand Up		@@ -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 bettercompatability with
		finding length of placed_eigs for bettercompatibility with
		different python versions [courtesy of Roberto Bucher]

		* src/delay.py (pade): New file for delay-based computations +
Expand Down
Original file line number	Diff line number	Diff line change
Expand Up		@@ -7,7 +7,7 @@
		from .statefbk import ctrb, obsv

		from numpy import zeros, shape, poly
		from numpy.linalg importinv
		from numpy.linalg importsolve, matrix_rank

		__all__ = ['canonical_form', 'reachable_form', 'observable_form']

Expand DownExpand Up		@@ -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

Expand Down
Original file line number	Diff line number	Diff line change
Expand Up		@@ -8,7 +8,7 @@
		margin.phase_crossover_frequencies
		"""

		# Python 3compatability (needs to go here)
		# Python 3compatibility (needs to go here)
		from __future__ import print_function

		"""Copyright (c) 2011 by California Institute of Technology
Expand Down
Original file line number	Diff line number	Diff line change
Expand Up		@@ -5,7 +5,7 @@
		Implementation of the functions lyap, dlyap, care and dare
		for solution of Lyapunov and Riccati equations. """

		# Python 3compatability (needs to go here)
		# Python 3compatibility (needs to go here)
		from __future__ import print_function

		"""Copyright (c) 2011, All rights reserved.
Expand DownExpand Up		@@ -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']
Expand DownExpand Up		@@ -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
Expand DownExpand Up		@@ -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
Expand DownExpand Up		@@ -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

Expand DownExpand Up		@@ -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
Expand DownExpand Up		@@ -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
Expand Down
Original file line number	Diff line number	Diff line change
Expand Up		@@ -40,7 +40,7 @@
		#
		# $Id$

		# Python 3compatability
		# Python 3compatibility
		from __future__ import print_function

		# External packages and modules
Expand DownExpand Up		@@ -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:]:
Expand All		@@ -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 - A12A22.IA21
		Br = B1 - A12A22.IB2
		Cr = C1 - C2A22.IA21
		Dr = sys.D - C2A22.IB2

		# 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
Expand DownExpand Up		@@ -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!")
Expand DownExpand Up		@@ -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
Original file line number	Diff line number	Diff line change
Expand Up		@@ -34,7 +34,7 @@
		# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
		# POSSIBILITY OF SUCH DAMAGE.

		# Python 3compatability
		# Python 3compatibility
		from __future__ import print_function

		import numpy as np
Expand Down
Original file line number	Diff line number	Diff line change
Expand Up		@@ -127,7 +127,7 @@ def acker(A, B, poles):

		# Make sure the system is controllable
		ct = ctrb(A, B)
		ifsp.linalg.det(ct)== 0:
		ifnp.linalg.matrix_rank(ct)!= a.shape[0]:
		raise ValueError("System not reachable; pole placement invalid")

		# Compute the desired characteristic polynomial
Expand All		@@ -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
Expand DownExpand Up		@@ -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];

Expand DownExpand Up		@@ -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())
Expand All		@@ -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