.57 . 74 B473 SLCO13373 DRC-5 MATRIX AUGMENTATION AND STRUCTURE PRESERVATION K IN LINEARLY CONSTRAINED CONTROL PROBLEMS by Johannes Bisschop and Alexander Meeraus Technical Note No. 5 Research Project 671-58 January 1978 nevelopment Research Center World Bank 1818 H Street, N.W. Washington, D.C. 20433 ABSTRACT Matrix augmentation is used for the inversion of bases associat with large linearly constrained control problems. It is shown how an efficient data structure can be maintained by keeping all state variables in the basis, and then nullifying some of them explicitly by using additional constraints. The proposed methodology, together with a basis updating scheme based on augmentation, forms the skeleton for an in-core algorithm using either the revised simplex method or the generalized reduced gradient method. Key Words: Large-scale programming Linear control problems Matrix modifications Partitioning methods Basis inversion I. Introduction The study of matrix modifications or matrix tearing dates back to the works of Kron [4] and Sherman and Morrison [7] . Consider the specific modification of a mx m matrix B where some of the columns of B are replaced by new columns previously not part of B . Let the m x p matrix U contain the new columns to be inserted in B , and the th p x m matrix V contain unit row vectors where the i row of V contains a one in the position of the column of B which is to be replaced by the ith column of U . Then B = B + (U - BV')V is the modified matrix, and the solution of Bz = b can be related directly to the solution of an augmented system of the form -B U] x] b [: f [:1 = En(1) V 0 y 0 Note that the equations Vx = 0 in (1) nullify explicitly all variables corresponding to those columns of B that will be replaced. One can verify that the solution vector z of Bz = b can be found by solving system (1) using the partitioned form of the inverse. Specifically, the identities B = B + (U - BV')V , Vx = 0 an'd VV' = I can be used to show that z = x + V'y where y = (VB-U)-1 VB-lb and x = B-1 (b - Uy) It should be noted that the above modifications are exactly the ones used in updating the basis inverse between major reinversions in the revised simplex method of linear programming. The matrix VB'U can be viewed as a giant pivot, and may be inverted by any stable method. This observation has resulted in a compact updating procedure where the growth of new nonzero elements is only a function of the number of iterations between reinversions and totally independent of the size of the original problem [2] This updating procedure is especially of interest in the case of structured problems where the basis inverse representation may be different from the traditional LU factorization. An example is any problem containing generalized upper bounds or an angular structure. This paper will extend the use of matrix augmentations and show how they can be employed for the reinversion of the basis matrix whenever the original problem contains a dontrol structure. When used in conjunction with the updating procedure described in (2] , a compact representation of the basis can-be maintained throughout the solution process. 2. A Simplified Control Problem In order to demonstrate the usage of matrix augmentation for basis inversion we have chosen to begin with a simplified version of a more general control problem. We will not consider any objective function in this paper, as the ideas developed here will apply to both the revised simplex method for linear programming and the generalized reduced gradient . method that one may employ in the case of a nonlinear objective function. Consider the first-order control equations of the form -3- Gx + Ax + Ku = z (2) t+1 t tt where- t=0,l,2,...,T-l , G is n x n and nonsingular, A is n x n and K is n x m . In many applications G is either the identity matrix, or a block diagonal matrix involving many 1 x 1 blocks. In each of these cases a compact representation for G-1 can be found. The tableau for system (2) , spelled out for all t-values, assumes the special structure displayed in Figure 1 , where ^0 = z - Ax0 and x 0 is specified as the initial condition of (2) Column Readings: 'l z2 z3 u ul E RS G K I z0 A G K A G 2 A . Figure 1 Tableau for Simplified Control.Problem The repeated occurrence of the matrices G , A and K suggest a generalization of the concept of supersparseness presently employed in large linear programming systems (3] . The tableau in Figure 1 can be stored with pointers to unique matrices G , A and K rather than to unique numbers. Such a storage scheme for the overall tableau is extremely -4- compact and especially useful when T , the number of t-values, becomes very large. We will modify this storage scheme slightly by having pointers to the columns of B rather than to the entire matrix B . This relatively minor change turns out to be necessary when trying to find an efficient storage configuration that can be maintained while inverting any Tn x Tn nonsingular submatrix of the tableau. Usually a scrambing of columns and rows is needed to build a compact representation of the inverse of a large sparse matrix [3} . To preserve our suggested storage scheme, only column permutations involving columns of B should be considered. By employing the problem structure and the idea of matrix augmentation, it is still possible to build a compact representation of any basis and its inverse using column permutations involving columns of B only. Consider a basis B and its augmented, partitioned form displayed in Figure 2 The matrix is partitioned into The upper-left block B has column headings for all elements of the vectors xt+1 and for selected basic elements of the vectors 'ut , t=0,l... ,T-1 . Since all x-variables are frozen into the basis, some of them may be nullified explicitly via constraints of the form Va xt = 0 and/or V xt+ 0 . The upper-right block B contains column headings for all basic elements of the vectors ut that could not be added to block Ba . Exactly how many basic components of ut are added to the upper-left block is determined at the time of inversion. -5- Colmi a a u Headings: 1 0 xt+1 0ut ut 0 OL- G K 0K z0 A0 I t t I V a0 A % VY' 0 ___ 0 figure 2-: Augmented Basis for Simplified Control Problem The idea is to add as many basic components of ut to the upper-left block as there are xt+1 components to be nullified, subject to the condition that Va G-1 K is invertible. If there are enough xt1 t t + components to be nullified and Va G K is nonsingular, the columns t t with heading ut will be vacuous. Any left-over nullifying constraints pertaining to elements of xt+1 are put into the lower-left block B. The nonsingularity of the matrices B and Bea guarantee that the number of rows of BT equals the number of columns of B and that the matrix B (B c)l B is nonsingular. The size of this matrix is not known a priori, but is equal to the number of unsuccessful matches between -6 - incoming ut components and nullified xt+1 components presently in the basis. Our effort to add as many nullifying constraints as possible to every t-block in Ba stems from the observation that computing a row of V' G-1 Kt requires considerable fewer computations that computing a row of By (B")-1 B . The size of Va G1 K a is p x p , where .t t Pt Xt 0 < p < min(m,n) for t=0,l,...,T-l . The partitioned form of the inverse of the augmented basis matrix B may now rely upon either full matrix technology or sparse matrix techniques depending on such parameters as size, and density. It should be noted that under the suggested augmen- tation scheme the matrix G-1 is determined only once, and then used over and over again for the inversion of any basis. One may make the following changes to the simplified control problem without affecting the proposed augmentation scheme at the time of inversion. The basic control equations may contain higher-order t-lags in both the state and the control variables. The proposed storage scheme for this problem is in principle the same as the one used for the first- order control problem. The matrices/columns of the lagged state/control variables will be carried along without affecting the size and the organization of the augmented basis. Another change to the basic control equations is the introduction of polynomial t-dependencies in the matrices A and K . If, for example, A = A(t) = A + tA + t A + ... is a simple polynomial of t , one may decide to store the matrices A0 , Al , A2 , ... together with the algebraic relationship that connects them. Again, the size and the organization of the augmented basis is unaffected, and the coefficient matrix A(t) in the inverse representation is computed for each t as it is needed. -7- 3. An Extended Control Problem The following problem statement distinguishes between four groups of equations: (i) first-order control equations; (ii) time integral equations; (iii) equations involving state variables only; and (iv) all other equations not belbnging to any of the first three categories. The problem statement does not contain simple upper bounds since they are handled trivially in both the revised simplex method and the generalized reduced gradient method. Consider i) Gxt+1 + Axt + Kut = t T T-1 2 2 ii) Cxt + Dut + s = w t=1 t=0 3 t-l (3) iii) Ex + s 2 v t t T T-1 3 3 iv) F x + H u + s = w t=l t=O - 1 2 3 The symbols s , s and s serve as possible slack (surplus) variables. The tableau for system (3) , spelled out for all t-values, is displayed in Figure 3 , where for notational convenience it is assumed that each of the equations ii) , iii) and iv) contains a slack (surplus) variable. As before, the repeated occurrence of several matrices suggest a special storage scheme. We propose pointers to unique matrices G , A C , and Ft , pointers to unique columns of K , D , and H , and t t I pointers to unique rows of E . Again, we use the problem structure and the concept of matrix augmentation to build a compact representation of -8- any basis and its inverse. Figure 4 displays a basis B in its augmented and partitioned form. The matrix is partitioned into B B[ Coltmn 1 2 2 3 'Headings: xl z2 u u u sIs a s R S G K A G K z A K z2 C C D D ID I v1 E I al E I 2 S22 71 2 0 H1 "2 Figure 3 Tableau for Extended Control Problem The notation in Figure 4 is complex as it distinguishes betwoen different types of partitionings. A closer inspection reveals the fact that the augmentation scheme is essentially the same as that in -9- Column a 2h a 82a U U * Headings:. xl 0 1 x2 1 21 u1 u G Kn a K0 00 0 0 I A G K C&K 21 1 11 VTa 1V A / Y 0 V.1 EC 1 C O C 1 I D F Hy H w 1 0 2 1 1 0 Figure 4 Augmented Basis for Extended Control Problem Figure 2 for the simplified control problem. The upper-left block B has been extended to include column headings for all basic s2 variables, while each t-block in Ba contains as many nullifying constraints and constraints of type iii) in (3) as possible. All other column headings are put in block B . It is interesting to note that the number of rows of V does not necessarily match the number of t - 10 - basic ua components as is the case in the simplified control problem. t The number of rows of Va and Ea added together, however, must equal t t a 2a the total number fo columns associated with u t ands The partitioned form of the inverse of B requires knowledge of B itself, plus the inverse of the matrices 0 01 [0- [Vt] G-1 [ K 0 ) and (-BS - By (Ba)-1B ) 0 I E a t Again, these inversions may rely upon either full matrix technology or sparse matrix techniques depending on size and density. Previous comments regarding the introductions of t-lags and polynomial t-dependencies in the matrices A , K , C , D , and .E apply here as well. There are several variations of the proposed storage and partitioning schemes that may be of interest. One alternation is to use a copy of K and G-1K . This will result in fewer computations, but may increase the required core storage whenever the number of nonzero elements determining G-1 is much less than the number of nonzero elements created by G'K . Our observations are by no means exhaustive in this matter, but merely point in a direction for those interested in software development. The main purpose of this note has been to show that matrix augmentations can play an important role in the design of a compact representation of a basis matrix and its inverse for the class of large linearly constrained control problems. REFERENCES [1] J. Abadie, "Application of the GRG-Algorithm to optimal control problems", J. Abadie, ed., Integer and Nonlinear Pro- gramming, (North Holland 1970). [2] J. Bisschop and A. Meeraus, "Matrix augmentation and partitioning in the updating of the basis inverse", Mathematical Programming 13 (1977) 241-254. [3] J.E. Kalan, "Aspects of large-scale in-core linear programming", - Proceedings of the 1971 Annual Conference of the Association for Computing Machinery (1971). [4] G. Kron, Diakoptics , (McDonald: London, 1956). [5] A.I. Propoi and V.E. Krivonozhko, "The dynamic simplex-method", IIASA Research Memorandum EM-77-24, Laxenburg, Austria, (1977). (6] D.J. Rose and J.R. Bunch, "The role of partitioning in the numerical solution of sparse systems", D,J. Rose and R.A. Willoughby, eds., Sparse Matrices and Their Applications, (Plenum Press: New York, 1972). [7] J. Sherman and W.J. Morrison, "Adjustment of an inverse matrix corresponding to changes in the elements of a given column of a.given row of the original matrix", Annals of Mathematical Statistics 20 (1949) 621.