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__importprint_function

 # External packages and modules


 #Check system is stable
 D,V=np.linalg.eig(sys.A)
 foreinD:
 ife.real>=0:
 raiseValueError("Oops, the system is unstable!")
 ifnp.any(np.linalg.eigvals(sys.A).real>=0.0):
 raiseValueError("Oops, the system is unstable!")

 ELIM=np.sort(ELIM)
 NELIM= []
 # Create list of elements not to eliminate (NELIM)
 foriinrange(0,len(sys.A)):
 ifinotinELIM:
 NELIM.append(i)
 NELIM= [iforiinrange(len(sys.A))ifinotinELIM]
 # A1 is a matrix of all columns of sys.A not to eliminate
 A1=sys.A[:,NELIM[0]]
 foriinNELIM[1:]:
 B1=sys.B[NELIM,:]
 B2=sys.B[ELIM,:]

 A22I=np.linalg.inv(A22)

 ifmethod=='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
 ifnp.linalg.matrix_rank(A22)!=len(ELIM):
 raiseValueError("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
 elifmethod=='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
 foreinD:
 ife.real>=0:
 raiseValueError("Oops, the system is unstable!")
 ifnp.any(np.linalg.eigvals(sys.A).real>=0.0):
 raiseValueError("Oops, the system is unstable!")

 ifmethod=='matchdc':
 raiseValueError ("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]

 returnH

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]:
 raiseValueError("System not reachable; pole placement invalid")

 # Compute the desired characteristic polynomial
 pmat=p[n-1]*a**0
 foriinnp.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
 returnK
 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)
 foreinD:
 ife.real>=0:
 raiseValueError("Oops, the system is unstable!")
 ifnp.any(np.linalg.eigvals(sys.A).real>=0.0):
 raiseValueError("Oops, the system is unstable!")

 iftype=='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
 returngram

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__importprint_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)
		foreinD:
		ife.real>=0:
		raiseValueError("Oops, the system is unstable!")
		ifnp.any(np.linalg.eigvals(sys.A).real>=0.0):
		raiseValueError("Oops, the system is unstable!")

		ELIM=np.sort(ELIM)
		NELIM= []
		# Create list of elements not to eliminate (NELIM)
		foriinrange(0,len(sys.A)):
		ifinotinELIM:
		NELIM.append(i)
		NELIM= [iforiinrange(len(sys.A))ifinotinELIM]
		# A1 is a matrix of all columns of sys.A not to eliminate
		A1=sys.A[:,NELIM[0]]
		foriinNELIM[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)

		ifmethod=='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
		ifnp.linalg.matrix_rank(A22)!=len(ELIM):
		raiseValueError("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
		elifmethod=='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
		foreinD:
		ife.real>=0:
		raiseValueError("Oops, the system is unstable!")
		ifnp.any(np.linalg.eigvals(sys.A).real>=0.0):
		raiseValueError("Oops, the system is unstable!")

		ifmethod=='matchdc':
		raiseValueError ("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]

		returnH
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]:
		raiseValueError("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
		foriinnp.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
		returnK
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)
		foreinD:
		ife.real>=0:
		raiseValueError("Oops, the system is unstable!")
		ifnp.any(np.linalg.eigvals(sys.A).real>=0.0):
		raiseValueError("Oops, the system is unstable!")

		iftype=='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
		returngram