Analysis Techniques for Large-Scale Electrical Systems - LU Decomposition & Sparse Matrices
This content covers the essentials of LU decomposition, Gaussian elimination, upper triangular matrices, LU decomposition theorem, corollaries, applications, and composite matrix representation in solving large-scale electrical systems. It discusses the process of triangularization in solving linear equations and various operations on matrices to facilitate efficient solutions.
Download Presentation
Please find below an Image/Link to download the presentation.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author.If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.
You are allowed to download the files provided on this website for personal or commercial use, subject to the condition that they are used lawfully. All files are the property of their respective owners.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author.
E N D
Presentation Transcript
ECE 530 Analysis Techniques for Large-Scale Electrical Systems Lecture 6: LU Decomposition; Sparse Matrices Prof. Hao Zhu Dept. of Electrical and Computer Engineering University of Illinois at Urbana-Champaign haozhu@illinois.edu 9/14/2015 1
Announcements Homework 2 posted, due Sep 23(Wed) 2
Gaussian Elimination A triangularization procedure to solving Ax = b Step k = 1, ,n a) eliminate rows 1, , k-1 from row k: kj kj km a = a a k k km b = b a (m) (m ) (m ) (m) mj 1 1 a j = m+ 1 ... , ,n (m) (m ) (m ) (m) m 1 1 b m = 1 2 ... , , , k - 1 b) normalize row k: a a a (k kj ) 1 (k) kj k + , ,n = = 1 ... j (k kk ) 1 b a (k k (k kk ) 1 (k) k b = ) 1 3
Upper Triangular Matrix The triangular decomposition scheme results in an upper triangular matrix U (1) 12 (1) 13 (1) 14 (1) 1n (1) 1 a a a a b 1 (2) 2 (2) 21 (2) 24 (2) 2 0 a a a b 1 n (3) 3 (3) 34 (3) 3 0 0 a a b 1 ?, ? = n (4) 4 (4) 4 0 0 0 a b 1 n (n) n 0 0 0 0 b 1 4
LU Decomposition Theorem We can prove that any nonsingular matrix A can be factorized into with the lower triangular matrix L A = LU ?11 ?21 (0) 0 ?22 ?32 ??2 0 0 0 0 (0) (1) (1) (2) ? = 0 0 ?33 ?31 ??1 (1) (2) ??3 Equivalently, ? = ? 1?; L contains the record of the operations preformed on A to obtain U 5
One Corollary The vector obtained after n steps of Gaussian elimination T ^ b (n) n (1) 1 (2) 2 b ,b , ,b = ... is related to the original b by the transformation ? = ? 1? 6
LU Decomposition Application As a result of this theorem we can rewrite Ax = LUx = b y = Ux Ly = b Define Then There exist multiple solution pairs since U could have non-unity diagonals Once A has been factored, we can solve for x by first solving for y, a process known as forward substitution, then solving for x using back substitution Note that y is exactly ? 7
LU Composite Matrix We can store the L and U triangular factors into a composite n by n matrix denoted by [L\U] u u un 11 12 13 1 u un 21 22 23 2 L\U = n n n nn 1 2 3 Hence i j i,j L\U = i,j u i< j i,j 8
LU Composite Matrix, contd The diagonal terms of U are not stored since they are known to be 1 s by the steps of the Gaussian elimination Using the composite matrix [L\U], there is no need to include the vector b during the process of Gaussian elimination since it has all the required information to compute A-1 b for any b 9
LDU Decomposition In the previous case we required that the diagonals of U be unity, while there was no such restriction on the diagonals of L An alternative decomposition is A = LDU L = LD with where D is the diagonal part of L, and the lower triangular matrix is modified to require unity for the diagonals 10
Symmetric Matrix Factorization The LDU formulation is quite useful for the case of a symmetric matrix = = = = U L A = U DU T A A = LDU A T T T U DL A T T Hence only the upper triangular elements and the diagonal elements need to be stored, reducing storage by almost a factor of 2 11
Pivoting An immediate problem that can occur with Gaussian elimination is the issue of zeros on the diagonal; for example 0 1 2 3 A = This problem can be solved by a process known as pivoting, which involves the interchange of either both rows and columns (full pivoting) or just the rows (partial pivoting) Partial pivoting is much easier to implement, and actually can be shown to work quite well 12
Pivoting, cont. In the previous example the (partial) pivot would just be to interchange the two rows 2 0 1 3 A = For LU decomposition, we need to keep track of the interchanged rows! Partial pivoting can be helpful in improving numerical stability even when the diagonals are not zero When factoring row k interchange rows so the new diagonal is the largest element in column k for rows j >= k 13
LU Algorithm Without Pivoting Processing by row We will use the more common approach of having ones on the diagonals of L. Also in the common, diagonally dominant power system problems pivoting is not needed Below algorithm is in row form For i := 2 to n Do Begin // This is the row being processed For j := 1 to i-1 Do Begin // Rows subtracted from row i A[i,j] = A[i,j]/A[j,j] // This is the scaling For k := j+1 to n Do Begin // Go through each column in i A[i,k] = A[i,k] - A[i,j]*A[j,k] End; End; End; 14
LU Example Starting matrix 20 12 12 3 5 6 A = 5 4 8 First row is unchanged; start with i=2 Result with i=2, j=1; done with row 2 20 12 0.25 9 4 3 5 A = 7.25 8 15
LU Example, cont. Result with i=3, j=1; 20 0.25 0.2 12 9 5.4 5 A = 7.25 7 Result with i=3, j=2; done with row 3; done! 20 0.25 0.2 12 9 0.6 5 A = 7.25 2.65 16
LU Example, cont. Original matrix is used to hold L and U 20 12 0.25 9 0.2 0.6 = = 5 = A = LU 7.25 2.65 1 0 1 0.6 0 0 1 With this approach the original A matrix has been replaced by the factored values! L 0.25 0.2 20 0 0 12 9 0 5 U 7.25 2.65 17
Forward Substitution b = Ly Forward substitution solves with values in b being over written (replaced by the y values) For i := 2 to n Do Begin // This is the row being processed For j := 1 to i-1 Do Begin b[i] = b[i] - A[i,j]*b[j] // This is just using the L matrix End; End; 18
Forward Substitution Example 10 Let = 20 b 30 1 0 1 0.6 0 0 1 = L From before 0.25 0.2 = = = [1] 10 [2] [3] y y y = 20 ( 0.25)*10 30 ( 0.2)*10 ( 0.6)*22.5 22.5 = 45.5 19
Backward Substitution Backward substitution solves (with values of y contained in the b vector as a result of the forward substitution) For i := n to 1 Do Begin // This is the row being processed For j := i+1 to n Do Begin b[i] = b[i] - A[i,j]*b[j] // This is just using the U matrix End; b[i] = b[i]/A[i,i] // The A[i,i] values are 0 if it is nonsingular End y = Ux 20
Backward Substitution Example 10 Let = 22.5 y 45.5 20 0 0 12 9 0 5 = U From before 7.25 2.65 = = = = [3] [2] x (1/ 2.65)*45.5 17.17 (1/9)* 22.5 ( 7.25)*17.17 x ( ) = 16.33 ( ) = [1] x (1/ 20)* 10 ( 5)*17.17 ( 12)*16.33 14.59 21
Computational Complexity Computational complexity indicates how the number of numerical operations scales with the size of the problem Computational complexity is expressed using the Big O notation; assume a problem of size n Adding the number of elements in a vector is O(n) Adding two n by n full matrices is O(n2) Multiplying two n by n full matrices is O(n3) Inverting an n by n full matrix, or doing Gaussian elimination is O(n3) Solving the traveling salesman problem by brute-force search is O(n!) 22
Computational Complexity Knowing the computational complexity of a problem can help to determine whether it can be solved (at least using a particular method) Scaling factors do not affect the computation complexity an algorithm that takes n3/2 operations has the same computational complexity of one the takes n3/10 operations (though obviously the second one is faster!) With O(n3) factoring a full matrix becomes computationally intractable quickly! A 100 by 100 matrix takes a million operations (give or take) A 1000 by 1000 matrix takes a billion operations A 10,000 by 10,000 matrix takes a trillion operations! 23
Sparse Systems This motivates us to exploit sparse matrix structure to reduce computational complexity For a sparse system, only nonzero elements need to be stored in the computer since no arithmetic operations are performed on the 0 s The triangularization scheme is adapted to solve sparse systems in order to preserve the sparsity as much as possible 24
Sparse Matrix History A nice overview of sparse matrix history is by Iain Duff at http://www.siam.org/meetings/la09/talks/duff.pdf Sparse matrices developed simultaneously in several different disciplines in the early 1960 s with power systems definitely one of the key players (Bill Tinney from BPA) In addition to the 1967 IEEE Proc. paper by Tinney see N. Sato, W.F. Tinney, Techniques for Exploiting the Sparsity of the Network Admittance Matrix, Power App. and Syst., pp 944-950, December 1963 25
Computational Complexity The computational order of factoring a sparse matrix, or doing a forward/backward substitution depends on the matrix structure Full matrix is O(n3) A diagonal matrix is O(n); that is, just invert each element For power system problems the classic paper is F. L. Alvarado, Computational complexity in power systems, IEEE Transactions on Power Apparatus and Systems, May/June 1976 O(n1.4) for factoring, O(n1.2) for forward/backward For a 100,000 by 100,000 matrix changes computation for factoring from 1 quadrillion to 10 million! 26
Inverse of A Sparse Matrix The inverse of a sparse matrix is NOT in general a sparse matrix We never (or at least very, very, very seldom) explicitly invert a sparse matrix Individual columns of the inverse of a sparse matrix can be obtained by solving x = A-1b with b set to all zeros except for a single nonzero in the position of the desired column If a few desired elements of A-1 are desired (such as the diagonal values) they can usually be computed quite efficiently using sparse vector methods (a topic we ll be considering soon) We can t invert a singular matrix (with sparse or not) 27
Scope of Sparse Matrices Sparse matrices arise in many areas, and can have domain specific structures Symmetric matrices Structurally symmetric matrices Tridiagnonal matrices Banded matrices ECE 530 is focused on problems that arise in the electric power; not a general sparse matrix course 28
Computer Architecture Aspects With modern computers the processor speed is many times faster than the time it takes to access data in main memory Some instructions can be processed in parallel Caches are used to provide quicker access to more commonly used data Caches are smaller than main memory Different cache levels are used with the quicker caches, like L1, have faster speeds but smaller sizes; L1 might be 64K, whereas the slower L2 might be 1M Data structures can have a significant impact on sparse matrix computation 29
Full Matrix versus Sparse Matrix Storage Full matrices are easily stored in arrays with just one variable needed to store each value since the value s row and column are implicitly available from its matrix position With sparse matrices two or three elements are needed to store each value The zero values are not explicitly stored The value itself, its row number and its column number Storage can be reduced by storing all the elements in a particular row or column together Because large matrices are often quite sparse, the total storage is still substantially reduced 30
Sparse Matrix Usage Can Determine the Optimal Storage How a sparse matrix is used can determine the best storage scheme to use Row versus column access; does structure change Is the matrix essentially used only once? That is, its structure and values are assumed new each time used Is the matrix structure constant, with its values changed This would be common in the N-R power flow, in which the structure doesn t change each iteration, but its values do Is the matrix structure and values constant, with just the b vector in Ax=b changing Quite common in transient stability solutions 31
Numerical Precision Required numerical precision determines type of variables used to represent numbers Specified as number of bytes, and whether signed or not For Integers One byte is either 0 to 255 or -128 to 127 Two bytes is either smallint (-32,768 to 32,767) or word (0 to 65,536) Four bytes is either Integer (-2,147,483,648 to 2,147,483,647) or Cardinal (0 to 4,294,967,295) This is usually sufficient for power system row/column numbers Eight bytes (Int64) if four bytes is not enough 32
Numerical Precision, cont. For floating point values using choice is between four bytes (single precision) or eight bytes (double precision); extended precision has ten bytes Single precision allows for 6 to 7 significant digits Double precision allows for 15 to 17 significant digits Extended allows for about 18 significant digits More bytes requires more storage Computational impacts depend on the underlying device; on PCs there isn t much impact; GPUs can be 3 to 8 times slower for double precision For most power problems double precision is best 33
General Sparse Matrix Storage Coordinate format: three vectors, each dimensioned to number of elements AA: Stores the values, usually in power system analysis as double precision values (8 bytes) JR: Stores the row number; for power problems usually as an integer (4 bytes) JC: Stores the column number, again as an integer New elements could easily be added, but costly to delete elements An unordered approach is not preferred since elements used next computationally aren t necessarily nearby Usually ordered by row preferred for our methods 34
Coordinate Format Example Assume 5 0 0 0 4 0 0 0 3 4 3 2 = A 4 3 2 10 1 1 2 = = = AA 5 4 4 3 3 2 4 3 2 10 Then JR 2 3 3 4 4 4 4 JC 1 4 2 4 3 4 1 2 3 4 Note, this example is a symmetric matrix, but the technique is general 35
Compressed Sparse Row Storage If elements are ordered (as was case for previous example) storage can be further reduced by noting we do not need to continually store each row number A common method for storing sparse matrices is known as the Compressed Sparse Row (CSR) format Values are stored row by row Has three vector arrays: AA: Stores the values as before JA: Stores the column index (done by JC in previous example) IA: Stores the pointer to the index of the beginning of each row 36
CSR Format Example Assume as before 5 0 0 0 4 0 0 0 3 4 3 2 = A 4 3 2 10 1 1 2 = = = 2 10 Then AA 5 4 4 3 3 2 4 3 JA 2 3 3 1 2 3 4 IA 1 3 5 7 37
CSR Comments The CSR format reduces the storage requirements by taking advantage of needing only one element per row The CSR format has good advantages for computation when using cache since (as we shall see) during matrix operations we are often sequentially going through the vectors An alternative approach is Compressed Sparse Column (CSC), which stores the values by column It is difficult to add values using CSR/CSC We ll mostly use the linked list approach here, which makes matrix manipulation simpler 38
