Skip to content

Commit c71222a

Browse files
slivingstonslivingston
committed
Merge pull request #101 "Refactored matrix inversions (see #85 and #91)"
#101 Changes are from branch `master` of https://github.com/mp4096/python-control.git There was merge conflict in how a for-loop was refactored into `map` (here) vs. list comprehension (from PR #110). I compared the two alternatives using %timeit of Jupyter for matrices that would correspond to LTI systems with 10 state dimensions, 2 inputs, 2 outputs (so, the A matrix has shape (10, 10), B matrix has shape (10,2), etc.), and with 100 state dimensions, 20 inputs, 20 outputs, all using matrices from numpy.random.random((r,c)). The difference in timing performance does not appear significant. However, the case of `map` was slightly faster (approximately 500 to 900 ns less in duration), so I decided to use that one to resolve the merge conflict.
2 parents d8b07cb + 71fef4b commit c71222a

14 files changed

+152
-76
lines changed

‎ChangeLog

+3-3
Original file line numberDiff line numberDiff line change
@@ -148,7 +148,7 @@
148148
2012-11-03 Richard Murray <murray@altura.local>
149149

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

153153
* tests/modelsimp_test.py, tests/slycot_convert_test.py,
154154
tests/mateqn_test.py, tests/statefbk_test.py: updated test suites to
@@ -604,7 +604,7 @@
604604
initial_response, impulse_response and step_response.
605605

606606
* src/rlocus.py: changed RootLocus to root_locus for better
607-
compatability with PEP 8. Also updated unit tests and examples.
607+
compatibility with PEP 8. Also updated unit tests and examples.
608608

609609
2011-07-25 Richard Murray <murray@malabar.local>
610610

@@ -994,7 +994,7 @@
994994
2010-09-02 Richard Murray <murray@sumatra.local>
995995

996996
* src/statefbk.py (place): Use np.size() instead of len() for
997-
finding length of placed_eigs for better compatability with
997+
finding length of placed_eigs for better compatibility with
998998
different python versions [courtesy of Roberto Bucher]
999999

10001000
* src/delay.py (pade): New file for delay-based computations +

‎control/canonical.py

+11-5
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
from .statefbk import ctrb, obsv
88

99
from numpy import zeros, shape, poly
10-
from numpy.linalg import inv
10+
from numpy.linalg import solve, matrix_rank
1111

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

@@ -68,23 +68,29 @@ def reachable_form(xsys):
6868

6969
# Generate the system matrices for the desired canonical form
7070
zsys.B = zeros(shape(xsys.B))
71-
zsys.B[0, 0] = 1
71+
zsys.B[0, 0] = 1.0
7272
zsys.A = zeros(shape(xsys.A))
7373
Apoly = poly(xsys.A) # characteristic polynomial
7474
for i in range(0, xsys.states):
7575
zsys.A[0, i] = -Apoly[i+1] / Apoly[0]
7676
if (i+1 < xsys.states):
77-
zsys.A[i+1, i] = 1
77+
zsys.A[i+1, i] = 1.0
7878

7979
# Compute the reachability matrices for each set of states
8080
Wrx = ctrb(xsys.A, xsys.B)
8181
Wrz = ctrb(zsys.A, zsys.B)
8282

83+
if matrix_rank(Wrx) != xsys.states:
84+
raise ValueError("System not controllable to working precision.")
85+
8386
# Transformation from one form to another
84-
Tzx = Wrz * inv(Wrx)
87+
Tzx = solve(Wrx.T, Wrz.T).T # matrix right division, Tzx = Wrz * inv(Wrx)
88+
89+
if matrix_rank(Tzx) != xsys.states:
90+
raise ValueError("Transformation matrix singular to working precision.")
8591

8692
# Finally, compute the output matrix
87-
zsys.C = xsys.C * inv(Tzx)
93+
zsys.C = solve(Tzx.T, xsys.C.T).T # matrix right division, zsys.C = xsys.C * inv(Tzx)
8894

8995
return zsys, Tzx
9096

‎control/margins.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
margin.phase_crossover_frequencies
99
"""
1010

11-
# Python 3 compatability (needs to go here)
11+
# Python 3 compatibility (needs to go here)
1212
from __future__ import print_function
1313

1414
"""Copyright (c) 2011 by California Institute of Technology

‎control/mateqn.py

+15-16
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
Implementation of the functions lyap, dlyap, care and dare
66
for solution of Lyapunov and Riccati equations. """
77

8-
# Python 3 compatability (needs to go here)
8+
# Python 3 compatibility (needs to go here)
99
from __future__ import print_function
1010

1111
"""Copyright (c) 2011, All rights reserved.
@@ -41,9 +41,8 @@
4141
Author: Bjorn Olofsson
4242
"""
4343

44-
from numpy.linalg import inv
4544
from scipy import shape, size, asarray, asmatrix, copy, zeros, eye, dot
46-
from scipy.linalg import eigvals, solve_discrete_are
45+
from scipy.linalg import eigvals, solve_discrete_are, solve
4746
from .exception import ControlSlycot, ControlArgument
4847

4948
__all__ = ['lyap', 'dlyap', 'dare', 'care']
@@ -557,9 +556,9 @@ def care(A,B,Q,R=None,S=None,E=None):
557556

558557
# Calculate the gain matrix G
559558
if size(R_b) == 1:
560-
G = dot(dot(1/(R_ba),asarray(B_ba).T) , X)
559+
G = dot(dot(1/(R_ba), asarray(B_ba).T), X)
561560
else:
562-
G = dot(dot(inv(R_ba),asarray(B_ba).T) , X)
561+
G = dot(solve(R_ba, asarray(B_ba).T), X)
563562

564563
# Return the solution X, the closed-loop eigenvalues L and
565564
# the gain matrix G
@@ -660,9 +659,9 @@ def care(A,B,Q,R=None,S=None,E=None):
660659

661660
# Calculate the gain matrix G
662661
if size(R_b) == 1:
663-
G = dot(1/(R_b),dot(asarray(B_b).T,dot(X,E_b))+asarray(S_b).T)
662+
G = dot(1/(R_b), dot(asarray(B_b).T, dot(X,E_b)) + asarray(S_b).T)
664663
else:
665-
G = dot(inv(R_b),dot(asarray(B_b).T,dot(X,E_b))+asarray(S_b).T)
664+
G = solve(R_b, dot(asarray(B_b).T, dot(X, E_b)) + asarray(S_b).T)
666665

667666
# Return the solution X, the closed-loop eigenvalues L and
668667
# the gain matrix G
@@ -699,7 +698,7 @@ def dare(A,B,Q,R,S=None,E=None):
699698
Rmat = asmatrix(R)
700699
Qmat = asmatrix(Q)
701700
X = solve_discrete_are(A, B, Qmat, Rmat)
702-
G = inv(B.T.dot(X).dot(B) + Rmat) * B.T.dot(X).dot(A)
701+
G = solve(B.T.dot(X).dot(B) + Rmat, B.T.dot(X).dot(A))
703702
L = eigvals(A - B.dot(G))
704703
return X, L, G
705704

@@ -825,11 +824,11 @@ def dare_old(A,B,Q,R,S=None,E=None):
825824

826825
# Calculate the gain matrix G
827826
if size(R_b) == 1:
828-
G = dot( 1/(dot(asarray(B_ba).T,dot(X,B_ba))+R_ba) , \
829-
dot(asarray(B_ba).T,dot(X,A_ba)) )
827+
G = dot(1/(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba), \
828+
dot(asarray(B_ba).T, dot(X, A_ba)))
830829
else:
831-
G = dot( inv(dot(asarray(B_ba).T,dot(X,B_ba))+R_ba) , \
832-
dot(asarray(B_ba).T,dot(X,A_ba)) )
830+
G = solve(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba, \
831+
dot(asarray(B_ba).T, dot(X, A_ba)))
833832

834833
# Return the solution X, the closed-loop eigenvalues L and
835834
# the gain matrix G
@@ -930,11 +929,11 @@ def dare_old(A,B,Q,R,S=None,E=None):
930929

931930
# Calculate the gain matrix G
932931
if size(R_b) == 1:
933-
G = dot( 1/(dot(asarray(B_b).T,dot(X,B_b))+R_b) , \
934-
dot(asarray(B_b).T,dot(X,A_b)) + asarray(S_b).T)
932+
G = dot(1/(dot(asarray(B_b).T, dot(X,B_b)) + R_b), \
933+
dot(asarray(B_b).T, dot(X,A_b)) + asarray(S_b).T)
935934
else:
936-
G = dot( inv(dot(asarray(B_b).T,dot(X,B_b))+R_b) , \
937-
dot(asarray(B_b).T,dot(X,A_b)) + asarray(S_b).T)
935+
G = solve(dot(asarray(B_b).T, dot(X,B_b)) + R_b, \
936+
dot(asarray(B_b).T, dot(X,A_b)) + asarray(S_b).T)
938937

939938
# Return the solution X, the closed-loop eigenvalues L and
940939
# the gain matrix G

‎control/modelsimp.py

+26-24
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@
4040
#
4141
# $Id$
4242

43-
# Python 3 compatability
43+
# Python 3 compatibility
4444
from __future__ import print_function
4545

4646
# External packages and modules
@@ -149,16 +149,12 @@ def modred(sys, ELIM, method='matchdc'):
149149

150150

151151
#Check system is stable
152-
D,V = np.linalg.eig(sys.A)
153-
for e in D:
154-
if e.real >= 0:
155-
raise ValueError("Oops, the system is unstable!")
152+
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
153+
raise ValueError("Oops, the system is unstable!")
154+
156155
ELIM = np.sort(ELIM)
157-
NELIM = []
158156
# Create list of elements not to eliminate (NELIM)
159-
for i in range(0,len(sys.A)):
160-
if i not in ELIM:
161-
NELIM.append(i)
157+
NELIM = [i for i in range(len(sys.A)) if i not in ELIM]
162158
# A1 is a matrix of all columns of sys.A not to eliminate
163159
A1 = sys.A[:,NELIM[0]]
164160
for i in NELIM[1:]:
@@ -177,14 +173,26 @@ def modred(sys, ELIM, method='matchdc'):
177173
B1 = sys.B[NELIM,:]
178174
B2 = sys.B[ELIM,:]
179175

180-
A22I = np.linalg.inv(A22)
181-
182176
if method=='matchdc':
183177
# if matchdc, residualize
184-
Ar = A11 - A12*A22.I*A21
185-
Br = B1 - A12*A22.I*B2
186-
Cr = C1 - C2*A22.I*A21
187-
Dr = sys.D - C2*A22.I*B2
178+
179+
# Check if the matrix A22 is invertible
180+
if np.linalg.matrix_rank(A22) != len(ELIM):
181+
raise ValueError("Matrix A22 is singular to working precision.")
182+
183+
# Now precompute A22\A21 and A22\B2 (A22I = inv(A22))
184+
# We can solve two linear systems in one pass, since the
185+
# coefficients matrix A22 is the same. Thus, we perform the LU
186+
# decomposition (cubic runtime complexity) of A22 only once!
187+
# The remaining back substitutions are only quadratic in runtime.
188+
A22I_A21_B2 = np.linalg.solve(A22, np.concatenate((A21, B2), axis=1))
189+
A22I_A21 = A22I_A21_B2[:, :A21.shape[1]]
190+
A22I_B2 = A22I_A21_B2[:, A21.shape[1]:]
191+
192+
Ar = A11 - A12*A22I_A21
193+
Br = B1 - A12*A22I_B2
194+
Cr = C1 - C2*A22I_A21
195+
Dr = sys.D - C2*A22I_B2
188196
elif method=='truncate':
189197
# if truncate, simply discard state x2
190198
Ar = A11
@@ -243,12 +251,8 @@ def balred(sys, orders, method='truncate'):
243251
dico = 'C'
244252

245253
#Check system is stable
246-
D,V = np.linalg.eig(sys.A)
247-
# print D.shape
248-
# print D
249-
for e in D:
250-
if e.real >= 0:
251-
raise ValueError("Oops, the system is unstable!")
254+
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
255+
raise ValueError("Oops, the system is unstable!")
252256

253257
if method=='matchdc':
254258
raise ValueError ("MatchDC not yet supported!")
@@ -373,8 +377,6 @@ def markov(Y, U, M):
373377
UU = np.hstack((UU, Ulast))
374378

375379
# Invert and solve for Markov parameters
376-
H = UU.I
377-
H = np.dot(H, Y)
380+
H = np.linalg.lstsq(UU, Y)[0]
378381

379382
return H
380-

‎control/phaseplot.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@
3434
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
3535
# POSSIBILITY OF SUCH DAMAGE.
3636

37-
# Python 3 compatability
37+
# Python 3 compatibility
3838
from __future__ import print_function
3939

4040
import numpy as np

‎control/statefbk.py

+7-8
Original file line numberDiff line numberDiff line change
@@ -127,7 +127,7 @@ def acker(A, B, poles):
127127

128128
# Make sure the system is controllable
129129
ct = ctrb(A, B)
130-
if sp.linalg.det(ct) == 0:
130+
if np.linalg.matrix_rank(ct) != a.shape[0]:
131131
raise ValueError("System not reachable; pole placement invalid")
132132

133133
# Compute the desired characteristic polynomial
@@ -138,7 +138,7 @@ def acker(A, B, poles):
138138
pmat = p[n-1]*a**0
139139
for i in np.arange(1,n):
140140
pmat = pmat + p[n-i-1]*a**i
141-
K = sp.linalg.inv(ct) * pmat
141+
K = np.linalg.solve(ct, pmat)
142142

143143
K = K[-1][:] # Extract the last row
144144
return K
@@ -242,7 +242,8 @@ def lqr(*args, **keywords):
242242
X,rcond,w,S,U,A_inv = sb02md(nstates, A_b, G, Q_b, 'C')
243243

244244
# Now compute the return value
245-
K = np.dot(np.linalg.inv(R), (np.dot(B.T, X) + N.T));
245+
# We assume that R is positive definite and, hence, invertible
246+
K = np.linalg.solve(R, np.dot(B.T, X) + N.T);
246247
S = X;
247248
E = w[0:nstates];
248249

@@ -354,10 +355,9 @@ def gram(sys,type):
354355

355356
#TODO: Check system is stable, perhaps a utility in ctrlutil.py
356357
# or a method of the StateSpace class?
357-
D,V = np.linalg.eig(sys.A)
358-
for e in D:
359-
if e.real >= 0:
360-
raise ValueError("Oops, the system is unstable!")
358+
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
359+
raise ValueError("Oops, the system is unstable!")
360+
361361
if type=='c':
362362
tra = 'T'
363363
C = -np.dot(sys.B,sys.B.transpose())
@@ -379,4 +379,3 @@ def gram(sys,type):
379379
X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
380380
gram = X
381381
return gram
382-

‎control/statesp.py

+10-5
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
99
"""
1010

11-
# Python 3 compatability (needs to go here)
11+
# Python 3 compatibility (needs to go here)
1212
from __future__ import print_function
1313

1414
"""Copyright (c) 2010 by California Institute of Technology
@@ -122,7 +122,7 @@ def __init__(self, *args):
122122
else:
123123
raise ValueError("Needs 1 or 4 arguments; received %i." % len(args))
124124

125-
A, B, C, D = [matrix(M) for M in (A, B, C, D)]
125+
A, B, C, D = map(matrix, [A, B, C, D])
126126

127127
# TODO: use super here?
128128
LTI.__init__(self, inputs=D.shape[1], outputs=D.shape[0], dt=dt)
@@ -447,11 +447,16 @@ def feedback(self, other=1, sign=-1):
447447

448448
F = eye(self.inputs) - sign * D2 * D1
449449
if matrix_rank(F) != self.inputs:
450-
raise ValueError("I - sign * D2 * D1 is singular.")
450+
raise ValueError("I - sign * D2 * D1 is singular to working precision.")
451451

452452
# Precompute F\D2 and F\C2 (E = inv(F))
453-
E_D2 = solve(F, D2)
454-
E_C2 = solve(F, C2)
453+
# We can solve two linear systems in one pass, since the
454+
# coefficients matrix F is the same. Thus, we perform the LU
455+
# decomposition (cubic runtime complexity) of F only once!
456+
# The remaining back substitutions are only quadratic in runtime.
457+
E_D2_C2 = solve(F, concatenate((D2, C2), axis=1))
458+
E_D2 = E_D2_C2[:, :other.inputs]
459+
E_C2 = E_D2_C2[:, other.inputs:]
455460

456461
T1 = eye(self.outputs) + sign * D1 * E_D2
457462
T2 = eye(self.inputs) + sign * E_D2 * D1

0 commit comments

Comments
 (0)