Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UF00084181/00001
## Material Information- Title:
- Automatic generation of solution procedures for indexed equation sets using GENIE
- Creator:
- Allen, Gary Louis, 1948-
- Publication Date:
- 1974
- Language:
- English
- Physical Description:
- vii, 160 leaves. : illus. ; 28 cm.
## Subjects- Subjects / Keywords:
- Algorithms ( jstor )
Data models ( jstor ) Index numbers ( jstor ) Index sets ( jstor ) Mathematical procedures ( jstor ) Mathematical variables ( jstor ) Matrices ( jstor ) Perceptron convergence procedure ( jstor ) Subroutines ( jstor ) Tears ( jstor ) Algorithms ( lcsh ) Chemical Engineering thesis Ph. D Dissertations, Academic -- Chemical Engineering -- UF Electronic data processing -- Equations ( lcsh ) GENIE (Computer program) ( lcsh ) - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt )
## Notes- Thesis:
- Thesis -- University of Florida.
- Bibliography:
- Bibliography: leaves 158-159.
- General Note:
- Typescript.
- General Note:
- Vita.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 000580528 ( ALEPH )
14046231 ( OCLC ) ADA8633 ( NOTIS )
## UFDC Membership |

Full Text |

AUTOMATIC GENERATION OF SOLUTION PROCEDURES FOR INDEXED EQUATION SETS USING GENIE By GARY LOUIS ALLEN A DISSERTATION PRES:NTKD TO THE GRA JA'\E CO^i(CII. O THE UiIIVERSily OiF FLORIDA S ,ARTiA! FULFIL!E' :" OF THE RP-,UIR ,M F ?:TS FO DEGREE OF DOCTOR CF PiILOSOPHY UNIVERSITY OF FLORIDA 1974 To Jan I ACKNOWLEDGEMENTS The author wishes to thank the chairman of his supervisory committee, Dr. A. W. Westerberg, Associate Professor of Chemical Engineering, for suggesting the research topic and for his guidance through the course of the research. Thanks are also extended to Dr. S. W. Director, Associate Professor of Electrical Engineering, and Dr. F. D. Vickers, Associate Professor of Computer and Information Sciences, for serving on the supervisory committee. The financial support of the Graduate School and of the College of Engineering is gratefully acknowledged. Thanks are also extended to the National Science Foundation for the financial support afforded by grants GK-18633 and GK-41606. In addition to assistantship support, these grants provided computer funds for the development of the programs in GENIE. The author accepted these funds with mixed feelings, as he believes that in a free society research should be supported by private concerns, rather than by the goi'ernet. The a.ithor extends special thanks to nis wife, Jan, iho not only provijdd encouragement when it was ,, os:t neded1, but also served as copy editor and typist for this dissertation. TABLE OF CONTENTS ACKNOWLEDGEMENTS.................................... ABSTRACT.............................................. CHAPTERS: 1 INTRODUCTION.................................. 2 INDEXED EQUATIONS............................. 2.1 Conventions Used in Describing Indexed Equations ............................... 2.2 Function-Variable Incidence Matrices.... 2.3 Index Display Matrices.................. 2.4 Index Imbedding......................... 2.5 Output Set Assignments ................. 2.6 Acceleration of Variables.............. 3 SOLUTION PROCEDURES FOR INDEX EQUATION SETS... 3.1 Definitions............................. 3.2 An Example Problem...................... 3.3 Output Sets and Decision Variables...... 3.4 Index Imbedding ........................ 3.5 Decoupling ............................. 3.6 Blocking Factors............................... 4 ALGORITHMS FOR DERIVING SOLUTION PROCEDURES FOR SETS OF INDEXED EQUATIONS ............................... 4.1 Decoupling.................................. 4.2 Index Output Set Assignments .................. iii vi 4.3 Function-Variable Output Set Assignment......... 4.4 Minimum Weighted Tearing..... ................. 5 COMPUTER IMPLEMENTATION ........................... 5.1 Data Structures.............................. . 5.2 A Versatile Memory System ..................... 5.3 Flow Diagram and Subroutine Descriptions........ 6 CONVERGENCE PROPERTIES .................................. 6.1 Review......................................... 6.2 Modification of Solution Procedures............ 6.3 Convergence Properties of Combined Gauss-Seidel and Newton-Raphson................ 7 EXAMPLE PROBLEMS .............................. ...... 7.1 Distillation Model ............................. 7.2 Examples ........................... ......... . 7.3 Discussion................................... . 8 CONCLUSIONS AND RECOMMENDATIONS...................... APPENDICES: A THE SIM DATA STRUCTURE.............. B SUBROUTINE DESCRIPTIONS............. BIBL IOG APHY.................................. BIOGRAPHICAL SKETCH ........................... Page 75 77 80 80 88 96 102 102 105 107 114 114 1 ! 123 131 135 139 158 160 . .. . w D t Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy AUTOMATIC GENERATION OF SOLUTION PROCEDURES FOR INDEXED EQUATION SETS USING GENIE By Gary Louis Allen June, 1974 Chairman: Arthur W. Westerberg Major Department: Chemical Engineering Chemical engineering design and simulation problems generally involve solving simultaneously a large number of equations. Many algorithms are available for structurally analyzing a set of equations with the goal of developing a solution procedure for the equations. Some of these algorithms have been implemented in a package of computer programs called "GENDER," an acronym for General Engineering Design Routines. Many chemical engineering processes are ma.iherr,atica'iiy mindeled with indexed equations. The index ranges are often largo, resulting in very large numbers of equations to be solved simultaneously. Conventional algorithms, such as those in GENDER, require that each equation to be analyzed by explicitly represented. This in itself is a difficult chore for indexed equations with large index ranges. A second shortcoming of existing algorithms is the time required to analyze such a large set of equations. The GENIE system, which is an acronym for GENDER with Indexed Equations, provides a more efficient means of analyzing sets of indexed equations. A compact means for representing the structure of indexed equation sets is presented. This representation allows both existing and newly developed algorithms to be applied in order to generate a solution procedure for the indexed equations. An evaluation of the convergence characteristics of the proposed solution procedure is made, using order of magnitude estimates of variable values. Should a proposed solution procedure be found to be non-convergent, means are provided for modifying it to enhance the likelihood of convergence. Formulas are presented which allow the estimation of the convergence properties of the new solution procedure. The solution procedures generated are independent of the ranges on the indices. Conventional methods derive solution procedures for a particular set of equations. The solution procedures are generated in the form of F0RTRAN "do-loops," so the index ranges become the do-loop parameters and can change from application to application. Solution procedures have been derived for a typical chemical engineering problem, a mathematical model of a disti llation column. Two different problems were solved by the solution procedure. GEN;I successfully predicted convergence difficulties for or= of the problems and modified the solution procedure. The modified solution procedure did solve the problem. CHAPTER 1 INTRODUCTION The design or simulation of a chemical process requires finding the solution to a large set of equations (many of them non-linear). Oftentimes these equations represent either a physical recycling of a process stream from one unit to another, or an interaction among variables within a unit; both conditions force a simultaneous solution of the equations. Since the equations are, in general, non-linear, an iterative solution is necessary. Frequently convergence of the itera- tive solution depends on the manner in which the solution is carried out. Unless the engineer has special familiarity with the particular set of equations he is trying to solve, he essentially must choose blindly from among the possible methods apparent to him. An alternative to this is to develop algorithms suitable for implementation on a digital computer capable of analyzing the structural and numerical properties of the set of equations with the aim oF finding a means of solving the equations which is not only efficient, but also likely to converge to the solution. The method chosen to solve a set of equations, whether chosen by the engineer or analytical algorithms, is called the "'solu-tion procedure." Two commonly employed types or solution procedures are tearing procedures (such as Gauss-Seidel) and "iNewton" type procedures (such as Newton-Raphson). Because of the cost (usually in computer time) of solving large sets of equations simultaneously, it is desirable to avoid Newton-Raphson type solution procedures whenever possible. The alterna- tive, tearing procedures, requires that certain variable values be guessed or torn (these variables become the iterative variables to be converged), all variable values are calculated algebraically including the torn variables, and the new values of the torn variables are used for the next iteration. The method continues until the criteria for convergence are satisfied. A tearing solution procedure requires an output set assignment for the algebraic solutions. An output set assignment is the assign- ment to each equation of one of the variables occurring in that equation is such a way that each variable is assigned to only one equation. Not all variables need be assigned to equations; the unassigned variables are "decision" variables. Early studies of the nature of output set assignments were concerned with structural proper- ties (Steward, 1962,1965). Later studies addressed the problem of assigning an output set with the aim of enhancing the likelihood of convergence (Edie, 1970) and developed effective algorithms for assign- ing outputs when some criterion for optimality has been defined (Gupta, 1972; Gupta et al., 1974). An iterative solution procedure is made iterative by the presence of tear variables which must be converged. iMany algorithms have been developed for choosing tear variables in order to satisfy any of several criteria, such as minimum number of torn variables or minimum sum of weights of torn variables, where each variable is given a weight (Sargent and Westerberg, 1964; Christensen and Rudd, 1969; Barkley and Motard, 1972; Ramirez and Vestal, 1972; Pho and Lapidus, 1973). The choice of decision variables has received comparatively less attention although the problems of simplifying the resulting solution procedure (Lee, Christensen, and Rudd, 1966) and avoiding singular systems (Edie, 1970) have been treated. The problem of deciding in what order to solve the equations has been solved (Sargent and Westerberg, 1964) to the extent that acyclic (non-recycle) systems can be discovered and equations involved in various recycle calculations can be identified. All of the algorithms developed thus far require an explicit representation of each equation and variable before the analysis can begin. For units described by a large number of equations (for example, a 50 plate distillation column with three components, which results in 350 equations) the total number of equations becomes quite large, and application of the above algorithms becomes time consuming and requires large amounts of computer storage space. In chemical engineering design and simulation the proliferation of equations and variables is quite often dlia to staged processes, which are modeled with indexed equations and variables, which are the same on every stage except for index value. Great savings of time and computer storage are realized by analyzing the set of indexed equations, each equation writtt.n enly once, rather than having to expand the set of equations and write separate equations for different index values. Unfortunately., existing algorithms are not directly applicable to these indexed eqaLtioris. Another shortcoming of existing algorithms is that the solution procedures derived for indexed equations (written in expanded form) would not be valid for an arbitrary range on index values. Algorithms for generating solution procedures for non-indexed equation sets have been programmed for use on digital computers by Cunningham (1972, 1973). The package of computer programs is called GENDER, for General Engineering Design Routines. The GENDER system is designed to provide an integrated computer program package for auto- matically generating solution procedures. The goal of this dissertation is to provide to the GENDER system a means of generating solution procedures for sets of indexed equations. The computer programs developed to do this are called the GENIE system, GENIE being an acronym for GENDER with Indexed Equations. Chapters 2 to 4 detail modifications both to existing algorithms and to the means of representing the set of equations which allow an automatic derivation of a solution procedure. The way that the data are represented within the computer storage, i.e., the data structure, has a profound effect on both the total storage required and on the time required to execute an algorithm. Chapter 5, along with Appendix A, discuss the various data structures used. Although the algorithms developed in Chapters 2 to 4 allow a solution procedure based on structural considerations to be derived, there is, unfortunately, no guarantee that the solution procedure will converge unless the system of equations fully precedence orders, i.e., is completely acyclic. It is desirable to base a solution procedure not only on structural considerations, but also on numerical consider- ations. Chapter 6 discusses some of the relevant numerical consider- ations and gives some insight into the prediction of convergence properties. 5 Chapter 7 gives some worked example problems, including distil- lation problems, which illustrate the importance of numerical consider- ations. Conclusions and recommendations are given in Chapter 8. CHAPTER 2 INDEXED EQUATIONS In order to discuss indexed equations a terminology must first be developed. This is especially true of the indices themselves since no existing terminology seems well suited to them. This is done in the next section of this chapter. Also, because it is desirable to have algorithms which efficiently derive solution procedures it was necessary to restrict both the types of equations allowed and the types of solution procedures which can result. Since it has been the aim to treat equations describing physical systems, those equations are among those allowed. The remainder of the chapter is devoted to further means of representing the equations and to a discussion of the restrictions imposed on the problems. 2.1 Conventions Used in Describing Indexed Equations The set of material balance equations described dhove might be written as: M 1 +v .-1 .-v .-f. = 0 2-1 3i I1 ij i-,] i+l 1,j ij where i = stage number i = component number M = name of function type, i.e., Material balance ., = liquid flow rate of component j leaving stage i v. = vapor flow rate of component j leaving stage i f.. = feed flow rate of component j entering stage i. 1i7 Thus, M.. represents the steady state material balances for all components on all stages. The name of an indexed function is used to distinguish functions of one function type from functions of another function type. Therefore, material balance functions, designated "M," and energy balance functions, designated "E," would belong to different function types. Similarly, variables occurring in indexed functions are classed by "variable type." For example, component liquid flow rate is of a different variable type from component vapor flow rate. It should be noted that in modeling a multi-unit process, material balance equations in a given unit contain different function types from material balance equations in another unit, even though the units are similar, e.g., distillation columns. It is convenient to make a distinction between the indices which subscript a function type and those which subscript a variable type. Consider equation 2-1. The indices "i" and "j" which subscript the function type M are said to be function indices. Similarly, the indices which subscript the variables 1, v, and f are said to be variable indices. In order to structurally analyze a set of equations it is not necessary to know the exact form of the equations. It is sufficient to kno;: wh.ih variables occur in which equations. All existing algorithms which perform structural analysis use only this information. This being the cas, ai further simplification of the representation of indexed euuations is made, with the ain of facilitating i2e structural analysis of the equation set. The convention is adopted that function indices will be designated by the letter "i." Thus, M. would become M,. ., the subscript on the SJ 1112 L subscript being used to distinguish between the two different function indices. Variable indices are designated by the letter "j." Since in the original representation of the function,the indices on the variables were functions of the indices on the functions, the variable indices must be expressed as functions of the function indices. First, some restrictions will be made on the possible values which the function indices can assume. Function indices will always have their possible values defined by a "range" which consists of a lower limit, an upper limit, and an increment. An index range is very similar to a FORTRAN do-loop range and consequently makes the eventual writing of the solution procedure in the form of do-loops convenient. Examples of function index ranges are: a) ii: L = 1 b) i2: L = 2 c) i3: L = 1 U = 10 U = 29 U = 13 A= 1 A = 1 A = 2 In case a) the function index would have possible index values of all integers between 1 and 10 inclusive, in case b) i2 would have possible index values of all integers between 2 and 29 inclusive, and in cazs c) 13 would have possible index values of all odd integer values between 1 and 13 inclusive. Although restricting the function indices to this form reduces the possible values that function indices in a function index set can have, it does not exclude any of the index sets an engineer is likely to create in modeling a chemical process. The stages in a staged process are numbered consecutively; a distillation column which has a third and fifth stage has a fourth stage. Certainly the index sets which arise from staged processes and from discretization of models are describable in terms of index ranges. More flexibility is desired in describing variable indices than is available in describing function indices. Variable index sets have their possible values defined by a range similar to the function index range or by a "list." The variable index range differs from the function index range in that its upper and lower limits will be permitted to be functions of either the function index value itself or the function index range limits. This restriction allows for many different arrangements of indexed variables. The equation f. = 0 = x. +x -x i=1,2,...10 2-2 1 i-i" +1 represented by a tri-diagonal incidence matrix, would have its index sets described as: i: L = 1 j: L = i-i U = 10 U = i+l A= 1 A= where both the lower and upper limits of the variable index range are defined as offsets from the function index value. Other examples of variable index ranges are: ji: L = L. j2: L ij3: L = L. U = U U. U U. i 1 A =1 A = 1 A 1 These; three examples illustrate the other three basic forms of the variable index range. The first, jil, defines a lower triangular incidence matrix, the second, j,, defines an upper triangular incidence matrix, and th third, J3, defines a full incidence mritrix. A second means of defining variable indices, the variable index list, allows even greater flexibility. Suppose we have an indexed equation of the following form: f. = x. +x.-x 2-3 S1i-2 1 1+1 The variable indices for this function cannot be expressed as an index range depending on the function index. In order to describe the index relationships in an equation of this type, we use a list. The variable index list is a list of offsets from the function index value; it has arbitrary length. The variable index list for equation 2-3 would have length 3 and would be: j: ii-2 il- There is some overlap between index lists and ranges. For example, either a list or a range can be used to describe the tri-diagonal relationship in equation 2-2. The index list and range provide for a structural description of a wide range of indexed functions, and in particular seem to describe any which are likely to occur in modeling any chemical process. 2.2 Function-Variable Incidence Matrices Since the incidence matrix is widely used to display structural information about a set of equations, it would be desirable to be able to represent a set of indexed equations in an incidence matrix. Unfortunately, for a typical problem such as a 20 stage, 4 component distillation column the number of rows would be 180 and the number of columns would be even greater. Incidence matrices of this size make analysis by hand time consuming and error prone. Even computer analysis becomes expensive, especially when much larger problems are encountered. A convenient, compact representation of the structure of a set of indexed equations is afforded by the "Function-Variable Incidence Matrix" (FVIM). The FVIM is not a true incidence matrix, and because of that has some unusual properties. A Function-Variable Incidence Matrix is a matrix whose rows each correspond to a function type and whose columns each correspond to a variable type. An element a.. of an FVIM exists if the variable type corresponding to j occurs in the function type corresponding to i, otherwise the element does not exist. An element which exists does not assume a value in the conventional sense. Instead, the element assumes the names of the variable indices which correspond to that variable type's incidence in the appropriate function type. As an example consider the following set of equations: f = 0 = x +x -y. +y 2-4 1112 i1-1,i2 1ii2 112 t1-2+1 g = 0 = x. +y -y gil1i 112 i li 2 +1,1i2 The Function-Variable Incidence Matrix would be: x y g(il,i2) 3?j4 1 324 it __-- I where I iL-1 J2: il h3" O iL i]+! J4: i2 Js: i2 i +2 In the position corresponding to function type f and variable type x is the element jij4. From the definitions of the variable indices we can deduce that in function type f we would find variables x. il- 2 or more compactly noted, variables x.. S13 L I It should be noted that variable index names cannot be assigned casually. An example will serve to point out a possible pitfall. The problem occurs when there is more than one variable index. Consider the following set of equations: f. = 0 = x +x. -y +y 2-5 1112 11-1,i2 1li2 ,1i2 11,i2+1 One might be tempted to assign four variable index names as follows: iJl i-1 j2:' il J3: i2 j4: i2 i- ii+1 12+1 This, however, does not accurately reflect the structure of the equations 2-5 for f. does not contain the variables y. Variables 11i2 j2j4 y j2j include yi. i+y. y ,i2. ,y.i+,i2"1, the second and third of which do not occur in equation f. This error can be 1112 avoided by several methods. First note that the problem arose from attempting to define a single set of variable indices for the variable y. Since this is not possible in this case, this equation could not be represented as it is in a Function-Variable Incidence Matrix. By rewriting the variables y as y and z, this shortcoming can be overcome. It should be noted that this type of problem is unlikely to occur, and does not occur in any of the several chemical process models examined in the course of this work. The Function-Variable Incidence Matrix then provides the means for compactly displaying the structure of a wide variety of sets of indexed equations and variables. Further, with slight modifications, virtually all sets of indexed equations and variables can be included. The FVIM will subsequently be shown to be valuable in effecting a speedy structural analysis of large sets of indexed equations. 2.3 Index Display Matrices The method used until now to describe variable indices, i.e., lists and ranges, is not well suited for a structural analysis of the relationship between function indices and variable indices. Since many existing algorithms are either designed to treat incidence matrices, or can be easily modified to do so, it is desirable to somehow convert the lists and ranges to incidence matrix representation. To do this the rows of the incidence matrix are made to correspond to the values of the function index and the columns are made to correspond to the values of the variable index. An element a.. is zero unless the variable index 137 can have the value corresponding to the column j when the function index is equal to the value corresponding to the row i. The resulting matrix is called an "Index Display Matrix" (IDM). The number of rows in an Index Display Matrix is determined by the range for the function index upon which it depends. Thus,.for the function index ii and the variable index j, which depends on it defined as follows: ii: L = 1 ji(ii): iu-1 U = 20 ii A = 1 il+i the Index Display Matrix (with zero elements blank and non-zero elements denoted by 'x') would be that of Fig. 2-1. This representation of the IDM occupies considerable space and is certainly larger than is neces- sary to convey the structural pattern either to a person studying it or to an algorithm analyzing it. In addition, as the range of the function index changes, so does the Index Display Matrix. In order to reduce the size of Index Display Matrices to something small enough to allow I I 1 I 1 I I i I 2 2 0 1 2345 6 7 890 12 345678 9 0 I X X X 4 XXX 2 XXX 3 XXX 5 XXX 5I ~XXX 6 XXX 7 XXX 8 XXX 9i XXX 10 X XX V \e V 3 XXX 14 X XX 15 XXX 16 X XX" 17 X X X 18 X, XX 2O1 ^ ^ FIGURE 2-1 INDEX DISPLAY MATRIX efficient structural analysis and to free the Index Display Matrix from its dependence upon the function index range the concept of the "blocking factor" is introduced. The blocking factor for a function index is the number of rows which will be included in the IDM for variable indices which depend upon that function index. If the function index ii above were assigned a blocking factor of 3 then the IDM for j1 would be: 0 1 2 3 4 1 x xx 2 X X X 2 x x x 3 x xx This Index Display Matrix conveys the notion of a tri-diagonal matrix while occupying considerably less space both on paper and in any data structure used in computer implementation. Since the blocking factor is arbitrary and is defined by the person defining a problem it should be pointed out that the blocking factor must be large enough to avoid any ambiguity in describing the IDM. In particular a blocking factor of 1 should never be used. 2.4 Tndex Imbedding The value of Function-Variable Incidence Matrices and Index Display Matrices is that while they are considerably smaller than the actual incidence matrix which they represent, they contain all of the inorrmation contained in the incidence matrix. This can best be illus- trated by the fact that the incidence matrix can be constructed from the Function Variable Incidence Matrix and its Index Display Matrices. First it is necessary to introduce the concept of index ordering. In naming the rows of an incidence matrix for a set of indexed equations it is desirable to have an orderly way to "step through" the various index values. The reasons for this are first, to insure that each function is represented and that no functions are included more than once,and second, to provide a logical means of stepping through the various functions which is easily adaptable to an iterative solution procedure. The method adopted is similar to the nesting of FORTRAN "do-loops." There is one important deviation, however, in that the function type is also treated as an index. The reasons for this modifi- cation become clear upon examination of the problem in Fig. 2-2. If we consider the three indices (function type, ii, and i2) as do-loop indices, obviously there are six ways in which these indices can be nested. This gives rise to six different structures for the indidence matrix. Here we have assumed that the ordering of the functions which would be dictated by the do-loop would be the order of the functions occurring in sequential rows. Also, the variables would be ordered similarly, i.e., each variable index would have the same nested position as the function index upon which it depends and the variable type would be assumed to depend upon function type. While it is obvious that the index I precedes the index 2 in a normal ordering, there is no such normal ordering when considering the function or variable type. It is not known whether f precedes g or whether x precedes y and at this point it does not really matter and the functions and variables will be put into an arbitrary order. Later the actual order will be dictated by precedence ordering considerations, similar to those applied to non- indexed equations and variables. The description of the index orderings Function-Variable Incidence Matrix x y a) f(il,i2) jlj2 J4j3 g(il,i2) jlj3 jlJ3 Index Definitions i,: L = 1 U=3 A = 1 Jl: L = 1 U = ii A= J4: L = 1 U = U. ii- i: L = l j2: i2-1 FIGURE 2-2 AN EXA:''LE PROBLEM U=2 A= i3: i2 i2+1 Index Display Matrices J2 J4 123 123 b) 1 x c) 1 x xx i 2 x x i1 2 x x x 3 x x 3 x x x J2 J2 012 123 0 1 2 1 2 3 d) 1 x x e) 1 x x i2 2 x i2 2 xx FIGURE 2-2 (cont.) (i.e., the order of nesting) in terms of do-loops makes the transition to computer implemented solution procedures in FORTRAN using do-loops straightforward. The advantages of using the do-loop in a solution procedure are primarily the compactness of code and the speed of execution. Indeed, most solution procedures implemented in order to solve indexed equations use do-loops or code similar to do-loops. In order to fully realize the consequences of the nesting of the indices and the effect of the Function-Variable Incidence Matrix and Index Display Matrices on the actual incidence matrix, examine Fig. 2-3. This is the expanded incidence matrix for the problem in Fig. 2-2. First look at Fig. 2-3 a), the entire incidence matrix. The matrix is partitioned into four submatrices, the divisions being chosen along the boundaries between function types and between variable types. These are the outermost nested "indices" and hence the ones examined first. The structure of the partitioned matrix is that of a 2x2, non-zero elements in all positions. That is the same as the structure of the FVIM, which has its rows and columns labeled by function and variable type. The next matrix, actually the (1,1) partition of Fig. 2-3 a), is itself partitioned along boundaries between the indices nested next innermost (corresponding to i2). Figure 2-2 d) exhibits the s;me structure as Fig. 2-3 b) (in a blockwise nature). This is what might be expected since Fig. 2-2 d) is the IDM for j2, which is a subscript of x in f. The pattern now becomes clear and as would be expected, Fig. 2-3 c) exhibits: the same structure as Fig. 2-2 b). Figures 2-4 a)-e) are representations of 2-3 a) achieved by varying the index orderings among the five other possible orderings. The type of behavior described above is evident in all of these incidence matrices also. Knowledge of this NM1rot0c K -iy NN0 r to N r0 -- 0 H0 Z 1 x x >. >, ,., ,, >.: K. x -,. b) fl1 f21 f12 f22 fs2 O000 011e0 0 N 8 0 N t 0oN N x t >s >c x: x x x x f e x X XX X X iX X I v X x x XX Xxx S -I d oo00o 21 IX X f3! -XXX- FIGURE 2-3 EXPANDED INCIDENCE MATRIX 00- f X :X N ck- c3^": fli f21 fS1 f12 f22 f52 eil g21 312 g22 g52 XX XX XX XXXXXX XX XX X XX XXX x %If r % XXXXXX XXX X S X X X X XX XX XX XX XX X XX X X XXX X ,_,__ X XX Xx___x x x x _x| _ _ii__ i_ i i i X . o C~ to oM r 0 m w ro cj to Nr- K) 0 e J rOr C0 V K o) r o -- \ cr N O ro to ro fi x x xx X Xxx f12 X X XX XX X X f21 XX Xx xX XX X f22 X X XX XX XX XX f. IX X X X X X X X X X X X f 2 XX X X X X X XX gfj XX xx 919 x x x X gl2 xx xx c21 xx xx xx xx g22 XX XX XX XX g1 XX X XX XX XX XX if, ix i2 0- N,. -> o0 o3--I rO 0 to 2 -- c Q cq cQ3 c" cj o j to r no ro r ') fi x xx x x 2 X x X X x X X gxs x x x xx f2 XX Y, Xx 91,/ A AX1 t2 XX XX XX XX XX f221 XX XX XX x x T2 x x xx xx xx x x fa xx xx xx x xX f32 XX XX XX XX XX XX '32. X X X X X X X X X, T?3 XX v/ v "^ B, / ? / '*- \ g ~ r! ?x I ____ A x^ A ,, A ^ __^. ^ AT.-^. il, if, 2 FIGURE 2-4 ALTERNATE INCIDENCE MATRICES o0- -- c~ MC r0 -- WM --w "z 0-- J rl SX X X XXX g11 xxx x Xa xxx x w x x x S xxx xxxx 9 x x x x z21 X XX X X X X X gz x x x X X X X X X f, xxx x Xxx x X X | fi xx x x x x x x x x x X^. xx Xx X 531 X X X X X X X X X X fx2 X XX X X XX X X XX X 9Q32 XXXX XXXX XXXX 'i l2,i'f fi f21 f31 g l 921 g31 f!2 f22 f32 g 12 g22 g 32 i ,if i FIGURE 2-4 (cont.) o_ 0 c\ _r- ro X XXX XXX x x xxI x xxx XX XX XXX XXX X X X X X X X X XX X X X X X_ XX XXXXX XXX x x X X X X xX XX XX XXX XXX x x x XX XX XX >X xx xx x xx X X X X X X X X X x d) e) i2viI lif FIGURE 2-4 (cont.) o 0 --- cJ 0 C C C N NN NN N tO"M tont - c 3 N M o CJ o to o- C< c tor X X >c >4,X > x >% x >% x >, x >% x A xx x x X x xx xx XXXX XXX x XXX X XXXX X X X X X XXXX X X X X I Xxxx X X X X XXXX XXXX X X X XXXXXX X X X .. x. x x XX x x x x x fi f21 g21 f3! t3i S}2 f?22 f32 932 behavior, called "Index Imbedding," allows the structure of the full incidence matrix to be predicted, and in fact visualized from knowledge of the Function-Variable Incidence Matrix and the Index Display Matrix. 2.5 Output Set Assignments A notion central to iterative solution of a set of equations is that of the output set assignment. Briefly an output set assignment is the assignment of one variable to each equation such that no variable is assigned to more than one equation. Variables which are not assigned are decision variables and have their values defined prior to the initiation of a solution procedure. There are, however, many output set assignments for a set of indexed equations which could be derived by existing algorithms operating on the full incidence matrix which would be extremely difficult, if not impossible, to solve using do-loops. As an example consider Fig. 2-5, where the circled elements are the out- puts. There is no reasonable or logical way to write a solution pro- cedure for that output set assignment using do-loops. If, however, we assign outputs to the Function-Variable Incidence Matrix and to the Index Display Matrices, a logical, concise solution procedure for a set of indexed equations can be implemented using do- loops. This means that, whenever a function of a given type is being solved, the output variable will always be of the type assigned to that function in the FVIM. Similarly, if a function being solved has index ii with a given value, the variable which is the output will have its index value corresponding to ii known. The IDM's which need to have outputs assigned (called "index outputs") are those which occur in output set assigned incidences of the FVIM. To illustrate the restric- cj C'] (\J (V (j 7'"* ri. II xxCi xx. .I'. s A / '^ 1 '*" i/; A 4> (I) A_ IR 2-5 > 0 a X A ;FI GRE I2 /.* ^ ..' >. A t ~ v \ *' ,^ * ^ ^~\ \'\" \ "f E 1 .,-*.} A A A ; *Vn- FIGUJRE 2-5 ARBITRARY OUTPUT SET ASSIGNMENT Fl 21 2 Gii Gi2 .i ;-' tions on output set assignments consider the system of equations from Fig. 2-5, defined structurally by f(ili2) g(i ,i2) x y JlJ4 J2J3 j274 J193 ji(il): L = L. U = i A= 1 j3(i2): i2 j2(il): L = L. U = U A 11 j4(i2): i2 U A with IDM shown index outputs the following fl1 f12 f21 f22 911 912 921 922 =2 =l i2+1 in Fig. 2-6. Assigning x to f indicated by the circles on the incidence matrix: X2 xll X12 X21 X22 Y11 Y12 Y13 Y91 and y IDM's ( x x x x ( x X X x x x ) x x x ( x x x x x x x ( x x x ( x x X X Y22 Y23 to g and making the in Fig. 2-6, we have with ii: L U A =1 = 2 =1 i2: L = 1 I 2 v J3 L 2 3 x o 12 I1 2 I x x X x J4 I 2 FIGURE 2-6 INDEX DISPLAY MATRICES A\ CX The variables marked with the "D" are decision variables. This is evident because they have index j3 equal to 1 and, by examining the index outputs, we see that for j3 the value 1 must be a decision. This output set assignment easily converts to a do-loop type iterative solution procedure such as: C READ DECISIONS AND TEARS READ (5,5000) Y 5000 FORMAT (10F8.0) 10 Do 20 IF=l,2 DO 20 11=1,2 DO 20 12=1,2 IF (IF .EQ. 1) X(I1,I2) = F(X,Y,I1,I2) IF (IF .EQ. 2) Y(I1,I2+1) = G(X,Y,I1,I2) 20 CONTINUE C CHECK CONVERGENCE If y not converged, go to 10 The convergence check would probably check the change in the torn variables y from iteration to iteration, convergence being defined as an arbitrarily small change in y. The functions "F" and "G" would be whatever f and g happen to be. Although these restrictions on the allowable output set assign- ments have been necessary to facilitate the writing of iterative solution procedures using do-loops, they do not weaken to any extent the solution procedure deriving power of the algorithms which operate under them. In fact these restrictions point the way, as it were, to the types of output set assignments engineers traditionally have employed. For example, in modeling a distillation column the material balance equations would either be solved for component liquid or vapor flow rates, not some liquid and some vapor flow rates. Similarly the flow rates solved for would always be those leaving the tray over which the material balance is being made. These are exactly the type of restric- tions which are being made on the FVIM and index outputs. In addition, because the outputs are made with small matrices (FVIM's and IDM's) existing algorithms can be used for the most part, which avoids the necessity of including in a computer program package, such as GENDER, separate routines for output set assigning indexed and non-indexed equations. Finally, by assigning the outputs to the small sub-matrices much effort is saved when compared to output set assigning the full incidence matrix. 2.6 Acceleration of Variables A possible consequence of defining variable indices with a range is a phenomenon called "acceleration of variables." Acceleration of variables is said to occur when the number of variable indices in a problem increases faster than the number of function indices as the range on a function index is increased. Acceleration of variables occurs at the lowest level of decomposition, i.e., on function indices and their associated variable indices. It is, of course, possible to have the number of variables in an entire problem increase faster than the number of equations and have no acceleration of variables. Accel- eration of variables is actually a possible result of the manner in which the indices are defined. To illustrate, consider the following function index and variable index: ii: L = 1 U = 2n+l A = 2 j1: L = L. U = U. A= For n=3 the IDM is 1 3 ii 5 7 1 2 3 4 5 6 7 x X X x x x x x x x x x x x x x with three more variables than equations. If n is increased to 4, i.e., if one equation is added, the IDM becomes j1 1 2 3 4 5 6 1 3 ii 5 7 9 7 8 9 which has four more variables than equations. Acceleration of variables produces the condition where a solution procedure must be derived for particular sets of function index ranges, and thus precludes any universal applicability of solution procedures. For this reason it is desirable to avoid acceleration of variables. The following theorem indicates how that can be insured. x X x X x x x x x x x x x x x x x x x x x x x x x X X X X XX X 31 Define the function index range as: L = M U = A.(N)-(A.-M) A = A. 1 where N and A. are positive constant integers, M a constant integer. Define the variable index range (it is only with the range that the problem can occur) as: L = aL.+(l-a)i+ki aE{O,l} U. = bU.+(1-b)i+k." bc{O,1} 3 2 3 A = A. 3 2 where k. and k .u are integers and A. is a positive integer, again all 3 3 3 constants. Theorem: If A. = A. there is no acceleration of variables. 3 1 Proof: Assume A. = A.. The number of functions is: j 1 A.(N)-(A .-M)-M Hf A. 1 = N The number of variables U. -L. n= max min +1 v A. U. U3max =bU,+(1-b)imax+k.u i 3 but imax = U. so, U = bU +(l-b) U+k . gmax 2 2 g = U .+k, . similarly, L = L .+k .' rmin J 7 U .+k .-(L .+k ) so n +1 A. U.-L. k.u-k . S2 + 3 +1 A. A. 1 2 k .-k . =n f+ 3 3 but since k.u, k and A. are all constants the difference between number of functions and the number of variables is always the same and hence there is no acceleration of variables. There can be no acceleration of variables for a variable index defined by a list since the difference between the number of functions (or function index values) and variables (or variable index values) is always the number of list elements minus one. CHAPTER 3 SOLUTION PROCEDURES FOR INDEXED EQUATION SETS The last chapter detailed some restrictions which are placed on the solution procedure developed. Other restrictions, along with their motivations, will be presented in this chapter. In order that efficient algorithms be written which can derive a solution procedure for indexed equations, it is desirable that theorems be developed which can be used to eliminate many of the possible choices open to the algorithms. This chapter presents many such useful theorems, along with comments on their application to problems. Most of the properties are related to a single example problem presented in the beginning of the chapter. 3.1 Definitions Before proceeding with the example it is helpful to define some of the terms used throughout the chapter. 3.1.1 FVIM Outputs As discussed in the last chapter, the FVIM must be output set assigned. The variable types chosen as outputs are called the "FVIM outputs." 3.1.2 FVIM Decisions In general, the number of variable types in an FVIM will be greater than the number of function types. When this is the case, some of the variable types cannot be assigned as outputs. These variable types must then be decision variables. The variable types which are decisions are called the "FVIM decisions." 3.1.3 Index Outputs The concept of index outputs was also introduced in Chapter 2. When the IDM is treated as an incidence matrix and output set assigned, the assignments made are called the "index outputs." 3.1.4 Index Decisions Just as the FVIM may have both outputs and decisions, so may the IDM. If any columns in an IDM remain unassigned after the index outputs are chosen those columns are the "index decisions." 3.1.5 Imbedded Loops The concept of index imbedding was introduced in the last chapter. The natural way to achieve index imbedding in a FORTRAN program is to define each index in a "do-loop" and nest the do-loops. Figure 3-1 illustrates a typical nesting of do-loops. The innermost loop, loop 3, is executed seven times for each pass through the next innermost loop, loop 2. Similarly, loop 2 is executed five times for each pass through loop 1. The index imbedding represented by this nesting of do-loops would be 13 imbedded in 12 which is in turn imbedded in Il. In terms of incidence matrices, each 12 incidence would actually represent an 13 IDM. The nested do-loops are termed "imbedded loops." The indices are said to be "nested inside" or "nested outside" the other indices. 3.1.6 Decomposed Problem and Expanded Problem As was seen in the last chapter, the FVIM and IDM's contain all of the structural information necessary to describe a set of indexed equations. The problem representation in terms of the FVIM and IDM's is called the "decomposed problem." When represented in terms of the 1 / 30 20 10 ---Q. I D D00( ---0 AT= 1210, 1 12= I,5, I3= I, 7,v -0 CONTINUE -20 CONTINUE -30 CNNTINUEU loop I loop 2 loop 3 end loop 5 end loop 2 end loop I FIGURE 3-1 NESTED DO-LOOPS full incidence matrix it is called the "expanded problem." 3.1.7 General Solution Procedure Existing algorithms generate a solution procedure for a particular set of equations. It is the aim of the algorithms here, not only to treat indexed equation sets, but to generate solution procedures which are independent of the index limits. A solution procedure which satisfies this requirement is called a "general solution procedure." 3.1.8 Function Ordering There is a natural ordering of function indices. They are either incremented or decremented in steps of equal size. There is no such natural ordering for the function types. In fact, the order of the functions is usually determined from precedence ordering considerations. The order in which the function types appear in the FVIM is called the "function ordering." The function ordering can change as the solution procedure is generated. 3.2 An Example Problem Suppose that a solution procedure is desired for this set of equations: f(il,i2) = 0 = f[x(j1,j3),y(jl,ji),z(jl)] 3-1 g(il,i2) = 0 = g[x(j2,j4),y(j2,j4),z(j2)] The decomposed representation of these equations is shown in Fig. 3-Z, which also defines the indices shown in the equations. To solve these equations, which may be assumed to be non-linear, a solution procedure must be specified. Suppose the following solution procedure is chosen. Assign x to f and y to g. Assign the index outputs as shown in the IDM's. Both FV M f ( i,i, ) Index 1, L= Definitions jr: j1 -l II A= I i2: L=I U=2 j- : L= Li U= Ui2 A=I SD s '2 4 i27, ," w C-, ! "i" ; v/ 1! FIGURE 3-2 DECOMPOSED PROBLEM J2: ii i,+i J4: i2 j3 2 |x x indices will be incremented and the function ordering will be as it is in the FVII. The index nesting will be il,if(function type),i2. The expanded incidence matrix resulting from these choices is shown in Fig. 3-3. It so happens that the equations fully precedence order for the solution procedure chosen. This example problem illustrates many things described in Chapter 2. Assigning x to f and y to g is an example of an FVIM output set assignment. Since z was not assigned it was an FVIM decision variable. Similarly the IDM's contain index outputs and index decisions. All variables with index ji were decisions for ji equal zero, and those with index j2 were decisions for j2 equal one. The problem decoupled in the variable type y for index i1. This example will be referred to by other sections in this chapter because it illustrates so many valuable concepts. 3.3 Output Sets and Decision Variables The choice of an output set and decision variable set almost entirely defines a solution procedure. This section discusses these choices as applied both to function and variable types and to indices. 3.3.1 Function-Variable Output Sets In the example problem a variable type was assigned to a function type. When this is the case, implementing the solution procedure is. simplified by the fact that each function type only has to be solved for one variable type. This is particularly helpful in computer implemented solution procedures, such as those to be generated by GENIE. For purposes of deriving solution procedures for sets of index equations, the requirement is made that outputs be assigned such that all outputs DDD D-D I I, I D DD - CIA" cQ CM Cki 0 0 0 C0 0 NC% C4 (C C\J C I fO rC O tO' -0 t~ C S x x x xlx) x x x xxxx xxx xx xxx x @ x Xx x x X X x __x x x _____ xxxx @xYx X X XX XXXX _x xx xx x X VX ," V ^7 A ^f\,.f) /A / ,, V XX v /^? v I. L~ Y 1, D = DECISION VARIABLE FIGURE 3-3 EXPANDED INCIDENCE MATRIX f l f12 gTI g2 f22 9 g22 T32 f32 gs32 of a particular function type are of the same variable type. This eliminates the need for more than one representation of a function type in a computer program. It also avoids having to determine during solution which output variable type would be required for a particular set of function index values. This restriction is a reasonable one to make, and actually is widely used in accepted solution procedures for typical modeling problems in chemical engineering. One shortcoming might appear to be the inability to assign different outputs to the equations describing different sections of a process, such as a distillation column. The rectifying and stripping sections are often treated separately and this restriction would seem to prevent that. This is not the case, however, since it is possible to describe the rectifying and stripping sections as separate problems with the appropriate connection equations. This has the added advantage that different solution procedures can be derived for the two sections. It should be noted that each output variable in the example problem had the same number of indices as the function for which it was the output. Clearly, z could not have been chosen as the output for either f or g since there are six equations of each type, but only five variables of type z. If the number of indices on a function type exceeds the number of indices on a variable type which occurs in that function, the total number of functions can be made to exceed the total number of variables by choosing appropriate index limits. For this reason, a variable type cannot be assigned as the output for a function type which has fewer indices than the function type. If the number of function indices were less than the number of variable indices for an output selection there could be many more variables than functions. While an output set assignment would be possible, determining the output variable index values for a particular set of function index values would be difficult or impossible. For this reason it will be required that the output variable type has the same number of indices as the function type for which it is the output. Suppose next that each of the variable indices for an output variable type does not depend upon a different function index. Then there must be a function index which does not appear in the output variable type. The range on this index could be increased until there are more functions than variables. This would mean that the variable type could not be the output. Therefore, each of the variable indices must depend upon a different function index. The following restrictions, then, are imposed on the FVIM output set assignment: 1. A separate variable type must be the output for each function type. 2. The number of indices for an output variable type must be the same as for the function type for which it is the output. 3. For an output variable type, each variable index must depend upon a different function index of the function type for which the variable type is the output. An interesting observation about the restrictions imposed on the output choice is that it is possible to have a set of equations which in expanded form has an output set assignment, and yet not be able to assign outputs to the FVIM. As an example, consider the constraint equation for the sum of the mole fractions in the liquid phase on a stage in a distillation column. NC LM. 0 = x..-l 3-2 j=-1 where: LM = Function Type i = Stage number j = Component number NC = Number of components x.. = Liquid mole fraction of component j on stage i The function type, LM, has one index. The only variable type is x, which has two indices. By the restrictions stated above, there is no possible output set assignment. There are two possible solutions to this problem. The first is to rearrange algebraically the equations in an attempt to eliminate the problem. The second is to assign, from the other variable types in the problem, a different variable type, with the correct number of indices, as an "implicit output." This variable would have to be connected to the function through the other function types (which means that through algebraic substitutions the implicit output could be made to appear in the function, although this is not done). A technique such as Newton-Raphson is then employed to converge the implicit output. 3.3.2 Variable Type Decisions In the example problem the variable type z was a decision variable type. The only restriction on the choice of decision variable types is that they be made so as not to reduce the number of variables below the number of equations. If the FVIM outputs are chosen before the decisions there can be no problem. If, however, any decision variable types are to be chosen before the FVIM outputs are assigned care must be taken to insure that the FVIM remains output set assignable. For the example problem, had x or y been chosen as the decision variable type, an output set assignment could not have been made. FVIM decisions may not be chosen so that the restrictions on FVIM outputs must be violated in assigning outputs. Consider the FVIM in Fig. 3-4 a). There are three function types and four variable types. Choosing x as the decision variable type results in the FVIM in Fig. 3-4 b) and choosing z results in Fig. 3-4 c). The FVIM in b) is not output set assignable according to the rules whereas the one in c) is. The decisions must be chosen so that there are the same number of variable type as function types with any given number of indices. 3.3.3 Index Outputs Index outputs are the outputs assigned in the Index Display Matrices. They indicate which of the variables of the output variable type will be the output for particular function index values. In extending index outputs from the IDM to a general solution procedure knowledge of the following properties is necessary. Theorem 3-1: All Index Display Matrices, when treated as incidence matrices are output set assignable provided they are not null. Proof: The proof is by induction and is divided into two parts, one for a list and one for a range. List: There must be at least one list element, which is an offset from the function index value. This is known from the definition of variable index lists. Therefore, for one function index value there is an output set assignment. Now suppose that there are output set assignments for an IDM with k rows. Row k+l introduces a new function index value, 1 a) f(il ,i2) J, J4 J2 J4 j2 i3 g (1,12) Ja J5 j I 4 J h(i ) j j5 ij j2 = . b) f(igi2) h(i ,) c) .'(i g(i, ,i2) ,'2) h(i, ) i J4 i2 I j3 j5 ij 3 1 2 w x y JI J4 J2 J4 J2 j3 j5 iJ 1 4 2 i5 j3 FIGURE 3-4 FVIM DECISION CHOICES larger than any of the other function index values. The largest offset in the list, when added to the new index value, must produce a new variable index value which can be assigned as an output. The first part is proved. Range: Define an index range as follows: i: L = L. j: L = a(i)+(l-a)L.+kz a{O0,1} U = U. U = b(i)+(1-b)U.+k bk{0,1} 2 1 U A = A. A = A. 2 1 where k and k are offsets. Note that L. cannot be greater than U. as this could result in a null IDM (for k=kU for example). Assume L. = U. i 2 then L. = k and U, = k In order that the IDM be non-null k must be greater than or equal to U k For one row the IDM is output set assignable since it is non-null. Now assume that the IDM is output set assignable for k rows. The number of rows can be expressed as: U.-L r i k A. and the number of columns as b(imax)+(l-b)U.+k -a(imin)-(l-a)L -k~ uc A --+ 1 3--3 but imadx = U, and imin = L. so eqn. 3-3 becomes Ui -L +k -ki k u + 1 a. k -k. k A. Now suppose that a new row is added. The number of columns becomes: U.+A.-L.+k -k Ck+1 + 1 3-4 1 = k + so there is a new column which can be assigned as the output for the new row. Theorem 3-1 is proved. Corollary: It is always possible to specify the index outputs as offsets from the index values in an IDM. Proof: For the list choose any list element as the index output offset. For the range note that for L.=U. the output can be expressed as an 7 1 offset from the function index value, say k Now assume that for L. U. is increased by A.. A new function index, the new U., 1 1 has been introduced. At the same time a new variable index value, U. +k not previously assigned as an output is introduced, which is offset from U. by the same factor as all previous outputs. The corollary is proved. The theorem and its corollary not only guarantee that index outputs are possible for arbitrary limits, but also that there is an index output set assignment which makes the determination of output variable index values easily obtainable from the function index values. 3.3.4 Index Decisions In an IDM the columns not assigned as index outputs are termed index decisions. Index decisions result from the lower limit offset for a variable index range being less than the upper limit offset, or from there being more than one list element for the variable index list. Theorem 3-2: For a general solution procedure, any index decisions must be declared as offsets from the function index upper and lower limits. Proof: Since for a general solution procedure the index range is arbi- trary, the only function index values which can be guaranteed to occur are L. and U.. Thus the only variable index values which can be 7 1 guaranteed to exist are offsets from L. and U.. (Remember that an 1 1 offset can be equal to zero.) The index decisions must therefore be chosen as offsets from L. and U.. Theorem 3-2 is proved. 1 2 To see the effect of index decisions on a problem recall the example problem at the beginning of this chapter. For j, the value zero (actually L. -1) was an index decision. All variables with index j, 11 equal to zero were decision variables. Similarly, all variables with index j2 equal to one were decision variables. 3.4 Index Imbedding As was seen in Chapter 2, the order of nesting chosen for the indices has a profound effect on the appearance of the expanded incidence matrix. Also, as in the case of the example problem at the beginning of this chapter, it will affect the efficiency of the solution procedure. First, index ordering (the order of nesting) can affect the number of function evaluations necessary to reach convergence. To see this, consider the incidence matrices c) and d) in Fig. 3-5, which are generated from the Index Display Matrices a) and b) with different index orderings. Suppose that these incidence matrices represent a set of equations to be solved for the outputs indicated by the circles. To FIGURE 3-5 DIFFERENT INDEX ORDERINGS FOR INDEX IMBEDDING 6ii Iv Ex xxB L-\ -{4 ('A ^-l* 'A*ut* ^^'w<^^^ t: S.? V 6 (1S/ A 5 bs) -, i~XXX ,. c--~caum0a 3. 0, V\` 'A, A- -X )r '< i A' A, '\ J S',: ". ,.. X X XX IkX XX 7 , Ki 'A. 'A I 'A 'A 'A .', 5.- ', s s s 'A A isi \f 'V v' '1 ' 1 ' X( j /* S 59 - - s. x'x x I. x i 6 f ' I i '^ ."t Ii t r S'V 'P v V . * -J \. a ' 49 Cm w. wVZX-I-llW YU- rm --r-_-^ --~- --Y t d) solve for the (1,1) variable in matrix c) the (1,5) and (1,9) variables must be torn and converged. To do this involves evaluating all of the functions, since the order of the rows in the incidence matrix is the order of solving the equations. To solve for the same variable in the second matrix involves evaluating only three functions until conver- gence. Since the equations in both incidence matrices are the same, the extra function evaluations in the first method would be wasted. The same number of iterations would be required to converge the first variable for both methods, but the number of function evaluations would not be the same. Theorem 3-3: Let il be a function index whose IDM does not fully precedence order and i2 be a function index whose IDM does fully precedence order. If the two indices are adjacent in index ordering, the number of function evaluations necessary for convergence of the set of functions they describe for the ordering i2il will be less than for the ordering ili2. Proof: Convergence is necessary for blocks of variables with i2 constant, iI ranging over its allowable values. Call the set of variables in a block with a particular value of i2 an is block. Let the number of iterations required to converge each i2 block be k, . '2 For the ordering i2il the number of function evaluations to converge the entire set of variables would be M2 M2 NI = (k i2. ) 1= M1 k. 3-5 i2=1 i= =1 Mi and M2 are the number of rows (and columns) in the IDM's for i1 and i2 respectively. Equation 3-5 states that the total number of function evaluations is the sum of the number of iterations for each i2 block times the number of functions in each i, block, over all such blocks. Now suppose that the order of the indices is ili2. The solution procedure would converge the first variable in the first i, block simultaneously with the first variable in each ii block. This is equivalent to converging all variables in the first i2 block. The number of function evaluations would be N2 = kiM1M2 3-6 since all functions must be evaluated. Equation 3-6 represents only the number of function evaluations required to converge the first variable in each iI block. The total number of function evaluations required would be M2 N = MM, I k. 3-7 T i 2 i2=1 N is greater than N1 by a factor of M2. Since the IDM for i2 does not fully precedence order, M2 must be greater than one. The theorem is proved. This theorem indicates the desirability of nesting indices whose IDM's fully precedence order outside those whose IDM's do not. A problem arises when a function index is to be nested outside the index for function type. Consider the incidence matrices in Fig. 3-6 a) and b). Any of the four partitions in the incidence matrix a) themselves have legal output set assignments. In the incidence matrix b), however, the (1,0) and (2,1) blocks do not have legal output set assignments. This condition is caused by ii being nested outside i The function g contains only variable index j2, and hence no variables with index value equal to zero. The result is that if the index output is to be assigned inde- pendent of function type it must be chosen from a special IDM. This IDM Nt? 0 *- C jO O cj fl) 1 / / SV v" | *"* ''' di ) } : t, n I si p 1 o c) -- C~ Cxl iro to A. s*** .^,; >' ^ >' t -- ) :,*.;, -ixiE. mea- r e-aev *mm "I |e.e f. '' 4 I I. I H fi i> *j :* 'f Fj I h F .. Li V i, t- 02~ tt~~' '"-~ FIGURE 3-6 INDEX DISPLAY MATRICES is the result of performing a logical "AND" operation on all IDM's which occur in incidences of the FVIM which are assigned as outputs. For example, refer to Fig. 3-6. If y is assigned to f and x to g, either ii or i1+i can be the index output. If, however, x is assigned to f and y to g, only ii can be the index output. The logical "AND" then provides a means for discovering which index outputs can be chosen inde- pendent of function type. When only structural considerations are being made, this can save time. This is because the logical "AND" can be performed faster than an output set can be assigned and because only one output set assignment is required in this case. An important result is that index ordering does not affect the number of tears in the problem. Consider the two different index orderings in Fig. 3-7. Both incidence matrices can be solved with four tears. Theorem 3-4: If two function indices are adjacent in ordering, the number of tears for the expanded incidence matrix is independent of the ordering of the indices. Proof: Let the IDM's for the two indices be called a and b. Let the number of tears for a be NT and for b be N ,. Let Nca and Ncb be the number of columns (and rows) in the IDM's a and b respectively. Assume b is nested inside a. Then the total number of tears is N = NTa x No. of columns for each occurence of a + NTb x No. of occurences of b not already torn -= Ta x Nb + NTb(ca NTa = NTa x N b + N Tb x Nca NTbNTa 3- which is symmetric in NTa, N Tb NTc and Ncb. Thus, the number of tears is independent of the index ordering. The theorem is proved. a) i T->* t f X Xi V S v B r, "-I '.^a J T T T T 4 o XX / " x x x x xxx x x xY X X X X NX X X X X v ! Si, GX XEX X3 I
x x":x x: j ^_^ x x' FIGURE 3-7 TWO DIFFERENT INDEX ORDERINGS b) tJ L [ /u ',7"'7:'!~ As shown earlier, a function index may have several different choices for the index output. It has also been shown that an index which fully precedence orders should be nested outside those which do not. If a function index has only one possible index output, it must fully precedence order. Since all IDM's are output set assignable, for a function index to fully precedence order, with only one index output, all of its IDM's must have the same index output and that index output must be the only possible one for each of the IDM's. This is equivalent to saying that the logical "OR" of all IDM's for a function index produces an IDM which fully precedence orders. This, then, provides an easy way to discover function indices which will give a full precedence ordering of blocks. It has been stated that function indices will either be incre- mented or decremented. Some variable indices are best treated by incrementing and some by decrementing. For example, a lower triangular IDM would fully precedence order for an incremented index and an upper triangular IDi would fully precedence order for a decremented index. Therefore, for some functions it may be desirable to increment a particular index and for others to decrement it. If the index in question is nested inside function type a simple test of function type can be made to determine which of the two methods is to be used. If, however, the function index is nested outside function type, no such test is possible and the index must be either incremented or decre- mented. 3.5 Decoupling Freedom in the choice of index decisions often affords the opportunity to choose a solution procedure much simpler than would be expected from an examination of the decomposed problem. Consider the decomposed problem in Fig. 3-8. At first glance it would appear that in order to solve for either x or y in f, the variable type not chosen as the output must be torn. This is not the case. By assigning y to f and x to g, and by making the index output assignments ii and i2+1, the incidence matrix in Fig. 3-9 results. For the outputs shown, the first nine columns represent decision variables. The first block of variables of output type y can be calculated without tearing any variables of type x. In fact, the entire set of equations can be solved without tearing any variables of type x. This phenomenon is called decouplingg." Decoupling is made possible by a judicious choice of index outputs. This problem is said to have decoupled in the variable x. Decoupling of a function type f in its output variable type x occurs when the function ordering would appear to require that the variable x be torn to solve functions not of type f, but the variable type x actually does not need to be torn for that reason. Theorem 3-5: A variable type x will decouple a function type f, for which it is the output, if and only if the following conditions exist: 1) A variable index for x must have index output offsets strictly greater than (or less than) the possible index output offsets for all other incidences of that variable in any other function type. 2) The index in 1) must be nested outside the index for function type. Proof: The proof will be presented for the "greater than" case. The proof for the "less than" case is analogous. If part The proof will be by induction. Assume 1) and 2) hold for some set of indexed equations. Let the index in question be called x y f(igi2) i j2 j4 ij g(il'i2 )Jl J3 J, J3 i i I x %x x x j2 X X 2 r x x j4 Tx x xB IXXX %x X x~ 13 x x 2x xX FIGURE 3-8 A DECOMPOSED PROBLEM fl 911 921 931 g, f12 f32 912 922 932 EXAMPLE OF DECOUPLING 0 0 0 J J c N NCJ CUN N re o tO M - cr cn -N N I x< X x X X X X X >< x x x x x xx x A xx x x xxx !xx Sx xx x x x x x xx Sx x R x XX X XXX Y x x x xx x x o xx xX X X X X x (X [ xx x xxxx y v w ,,,'" x x x x Ix x x x( x x xI F E 3if 1 FIGURE 3-9 ii. For il=L. all incidences of x in functions other than f which have 11 an index depending on ii, have index values less than the lowest valued index output (by condition 1), and thus are decisions. Since 1i is nested outside function type (by condition 2), the first pass for il can be completed without tearing x for functions other than f. Now assume that k passes through the i1 loop have been completed without tearing any x's for functions other than f. The index outputs for il in x, for iteration k, were greater than any of the other index offsets for ii in x. Thus these index outputs of x are greater than or equal to any of the index values for x in functions other than f for iteration k+l. The functions other than f, then, do not require any tear of x for iteration k+l. The if part is proved. Only if part The proof will be by contradiction. Assume 1) does not hold. The index output offset for x in f is less than or equal to some other index value for il in an incidence of x in a function, other than f, which must be solved before f. The variable x in this function must be torn for iteration k since there are index values of x present which will be outputs for iteration k. This contradicts the definition of decoupling, so condition 1) is required for decoupling. Assume 2) does not hold. Then a function not f, which contains the variable x must be solved for all values of ii before any function f is solved. This would require tearing x which would contradict the definition of decoupling, thus condition 2) is required. The theorem is proved. This theorem provides the basis for algorithms which discover conditions which allow decoupling. 3.6 Blocking Factors The blocking factor defines the size of the Index Display Matrices for the various function indices. It is desired to perform analyses on the IDM's and have the results of those analyses be applicable to an expanded incidence matrix with arbitrary index limits. It is possible, by employing different blocking factors, to analyze an IDM in the same manner twice and reach different results. For example, the blocking factor can affect the ratio between the number of tears and the number of rows in a problem. (See Fig. 3-10.) In this figure, a) and b) are Index Display Matrices representing a tri-diagonal matrix with first and last columns declared as decisions and omitted from the figure. For a blocking factor of 2, either of the two possible output set assignments results in one tear for the two rows. When the blocking factor is increased from 2 to 3, two of three possible output set assignments only require one tear for the three rows. The require- ment that the order of solution be dictated by the IDM (i.e., be top to bottom or bottom to top) is relaxed for the analysis within the blocks. When the problem is expanded, the IDM for an index is made up of many blocks. These blocks must be solved in either ascending or descending order. It is not always the case that the solution procedure derived for a block will hold for the expanded IDM. Consider the two expanded IDM's in Fig. 3-11, each representing one of the one-tear, blocking factor of 3, solution procedures from Fig. 3-10. The first IDM reflects a solution procedure which still allows for solution of the equations by solving the (1,1) block then the (2,2) block. The second, however, does not allow for solving either block before the other. When a solution b) X X X X T T vR l x LX@X xx x8 T = TEAR VARIABLE FIGURE 3-10 ALTERNATE BLOCKING FACTORS a) X X XX T xi IxI T x x xx X 6) A A l i ! xx xx xl x x x x FIGURE 3-11 EXPANDED INDEX DISPLAY MATRICES I i x procedure for a blocking factor remains effective as the IDM expands, the solution procedure is said to "extend." Theorem 3-6: Call a block with blocking factor equal to k a k-group. For the solution of k-groups in ascending order the following condition is necessary and sufficient for extension of a solution procedure. The decision variables offset from U. must all be k greater than tear variables within the last k-group. An analogous condition exists for solution of k-groups in descending order. Proof: The proof for ascending order only will be presented. Analogous arguments hold for the case of descending order. Sufficient: Suppose that all decision variables offset from U. are offset from the k-group tear variables by k. When the next k-group is added, the decisions offset from U. are no longer decision variables. 1 In order to solve the next to last k-group, values for these variables must be available. The values are available since these variables have become tear variables. Thus if a solution procedure holds for n k-groups it holds for n+l k-groups. It must hold for one k-group, otherwise it would not be a solution procedure. The condition is sufficient by induction. Necessary: Assume that some decision variable offset from U. is not k greater than a tear variable in the last k-group. When the next k-group is added, this variable's value must be available to solve the next to last k-group. It is no longer a decision and is not a tear, but is an output of a function in the last k-group and hence its value is not yet calculated. For its value to be known it must be a decision or tear variable, but this is a contradiction. Thus the solution procedure does not extend and the necessary condition holds by contradiction. The theorem is proved. This theorem provides an easy method for determining which solution procedures to IDM's are worth analyzing and which are not. If there are no decision variables (decision indices actually) involved, the analysis is even simpler. For a variable index defined by a list with no decisions, letting L.=U. proves there is only one list element. Thus for the index there is a full precedence ordering. For a variable index defined by a range, four cases must be considered. First consider L. and U. offset from L. and U. by the same factor (to result in no 3 3 1 1 decisions). In this case a solution procedure does not extend since, for any blocking factor k, k-1 tears are made in the first block and k tears must be made in the second block. Next consider L. offset from L. -7 and U. offset from i, again by the same factor. A full precedence 3 ordering exists as this results in a lower triangular IDM. Similarly L offset from i and U. offset from U. results in full precedence order. -7 .7 With both limits offset from i by the same amount the IDM is diagonal and fully precedence orders. Knowledge of the behavior of solution procedures from index definitions is very useful in choosing among the various solution procedures. Many possible solution procedures can be rejected without extensive analysis because they are known to be inefficient for the problem at hand. CHAPTER 4 ALGORITHMS FOR DERIVING SOLUTION PROCEDURES FOR SETS OF INDEXED EQUATIONS The algorithms implemented in the GENDER system for deriving solution procedures are capable of treating non-indexed equations only. In deriving a solution procedure GENDER performs the following analyses. First the acyclic algorithm is applied to determine which parts of the problem are acyclic and which contain recycle calculations. An output set assignment is then made. After the outputs are assigned the precedence ordering is found and the groups of functions and variables requiring simultaneous solution are identified. The last step is to choose the tear variables for each group. The treatment of indexed equations by GENIE is similar in approach to that of GENDER. First, the FVIM is analyzed by the acyclic algorithm to discover the acyclic nature of the set of equations. The decoupling algorithm is then applied to the set of function and variable types involved in recycle calculations in an attempt to identify further acyclic structure. Outputs are next assigned to both function types and function indices. A precedence ordering of function types is determined and groups are identified. Finally, the tear variables within each group are determined by a minimum weighted tearing algo- rithm. The algorithms presented in this chapter perform the particular analyses necessary for deriving a solution procedure for indexed equations. The manner in which these algorithms are used is explained in section 5.3. Some of the algorithms are new, while some are adap- tations of algorithms previously developed. Where existing algorithms are used, the alterations to the problems necessary to allow their use are detailed. 4.1 Decoupling Decoupling can result in significant simplification of the solu- tion procedure for a set of indexed equations. Since the presence of decoupling depends upon the function-variable and index output set assignments made, it is desirable to choose these in such a way that decoupling will result. In addition, the function ordering also affects the degree to which decoupling is possible. Finally, the direction in which the function indices are incremented (ascending or descending), influences any decoupling. For a particular direction of index incrementing, some IDM's may fully precedence order while others may not. All IDM's for a function index must fully precedence order for it to be a candidate for a decoupling index. Recall that decoupling results in a full precedence ordering of function types. This is equivalent to a full precedence ordering for an index, the index in this case being function type. It has been shown that nesting a fully precedence ordered index inside one which does not fully precedence order results in increased function evaluations. This would be the case if a decoupling index did not fully precedence order. What is required, then, is an algorithm which will analyze the FVIM and IDM's to choose the function ordering, function outputs, index outputs, and direction of index incrementing in such a way as to decouple the functions if possible. The following algorithm does that. (An example which follows illustrates the algorithm.) I. Produce the FVIM. Flag all function indices whose IDM's do not precedence order. II. A. Break all assignments made by this algorithm. B. Choose a set of directions for incrementing the function indices not previously analyzed. If there are none, go to III.C. C. Untag all function and variable types and function indices. III. A. If the set of unassigned function and variable types does not fully precedence order, go to IV.A., otherwise, continue. B. A full precedence ordering has been discovered. Stop. C. Decoupling will not produce a full precedence ordering. Stop. IV. A. Choose an untagged, unassigned variable type. If there are none, go to II.A. Tag the variable type. Untag all function indices and unassigned function types. B. Choose an untagged, unassigned function type. If there are none, go to IV.A. Tag the function type. C. If the variable type chosen is not a legal output of the function type chosen, go to IV.B. V. A. Choose the first untagged function index. If there is none, go to IV.A. Otherwise, go to VI.A. B. Choose the next untagged function index. If there is none, go to IV.B. VI. A. Does the variable index depending upon the chosen func- tion index decouple the function type from the unassigned functions? Tag the function indices which fail for all functions. If no decoupling, go to V.B. B. Assign the chosen variable type to the chosen function type. Give the function type a position in the function ordering immediately after all of the presently unassigned functions. C. Make the index output set assignments which produced the decoupling. Flag the function index as a decoupling index. Go to II.C. The decoupling algorithm performs an exhaustive search over all combinations of directions of index incrementing. Additionally, when- ever an output assignment is made, all remaining function and variable types are untagged and the algorithm restarted on the remaining unassigned function and variable types. Generally, an exhaustive search of this type would be extremely long. In the case of the decoupling algorithm, the search will be relatively short since the FVIM will usually have only a few function and variable types and there are usually only one or two function indices. For some of the function index directions it may be possible to reject certain variable types at the outset. Remember that the aim of decoupling is to manipulate the choice of solution procedure so that variables which would otherwise have to be torn do not have to be torn. Suppose that the function index ii is incremented in ascending order. Suppose, also, that the variable type x is being examined to determine if it can be used to decouple a function type from the others. If x has any FVIM incidence, say in function f, with a variable index defined by a range with upper limit offset from U. decoupling will not hold for that function index. The reason is that, for iz=Li the first pass through the "ii loop," the function f will require values of x for all values of ii. This being the case these values must be supplied by tearing x, hence decoupling could not hold. An analogous condition holds for the function index being decremented rather than incremented. Another "shortcut" is to note that, if for some variable type all variable indices depending upon a particular function index are the same, decoupling cannot occur for that function index. Figure 4-1 presents an example problem to which the decoupling algorithm will be applied. The first step is to produce the FVIM, which is at the top of the figure. No assignments need be broken as none have been made. Choose ascending order for both indices initially. The function and variable types are untagged. Now the algorithm begins in earnest. The FVIM does not fully precedence order, so the algorithm skips to step IV.A. The following steps are then carried out: IV. A. Choose x. Tag x. IV. B. Choose f. Tag f. IV. C. Legal output; continue. V. A. Choose ii. VI. A. No decoupling (all variable indices are j1). Tag i V. B. Choose i2. VI. A. No decoupling (U. =U. ). Tag i . 34 12 V. B. No more function indices. IV. B. Choose g. Tag g. IV. C. Legal output; continue. 1 f(iI ,i2) g(i, ,i2) ig L= U=3 A= U2 L=I U=2 A=I i ii j3 i2 J2 L= Li L =Li, A= I FIGURE 4-1 EXAMPLE PROBLEM FOR DECOUPLING J J4 jz J3 i J j3 j4 V. A. No untagged indices. IV. A. Choose y. Tag y. Untag f and g. Untag i1 and i2. IV. B. Choose f. Tag f. IV. C. Legal output; continue. V. A. Choose ii. VI. A. No decoupling. (No possible index output for j2 greater than all index values for ji.) V. B. Choose i2. VI. A. No decoupling (fails for all functions since U j=U ). Tag i2. IV. B. Choose g. Tag g. IV. C. Legal output; continue. V. A. Choose ii. VI. A. Decoupling. VI. B. Assign y to g. Solve g after f. VI. C. Assign iii+ as index output for i1 in g. II. C. Untag everything. III. A. There is a full precedence ordering of unassigned function and variable types (that is, after column For y and row for g deleted). III. B. Stop. The expanded incidence matrix for this example, arranged to reflect the decoupling, is shown in Fig. 4-2. The index output for ii in f is chosen to be i1+1 and for i2 the index output is chosen to be i2, all choices being arbitrary. There are only four tear variables in this problem, and the convergence loops (each 2x2) are the smallest possible. DECISIONS x x x xI I , ,x X jx x x x xx x r V- xx xxx x x x x I x x x I xxx x@q FIGURE 4-2 A DECOUPLED PROBLEM fi Ti2 ci g2 -C fa 2g, 931 g32 Although full decoupling may not be possible for a set of indexed equations, a partial decoupling may be. A partial decoupling results when some function types are decoupled from the set of equations, but not enough to render the FVIM acyclic. As a result a subset of the original function and variable types must be solved simultaneously. If full decoupling is not possible, there will, in general, be a partial decoupling for each set of index incrementing directions. Since a partial decoupling results in a smaller set of equations which must be solved simultaneously than does no decoupling, it would be desirable to choose one of the partial decoupling schemes for the solution procedure. It is possible to characterize the "degree of decoupling" by the expected number of functions decoupled (where the expected number for a function is calculated from its expected index ranges, not the blocking factor). The different partial decoupling schemes could then be compared quantitatively and the one exhibiting the greatest degree of decoupling chosen for the solution procedure. This is the strategy adopted by GENIE. 4.2 Index Output Set Assignments Index outputs are assigned by treating the IDM's for each variable index as an incidence matrix. The IDM's are output set assigned by conventional output set assignment methods. The resulting output set assignment is then available for defining the index outputs, should that be necessary. As indicated in Chapter 3, it is sometimes desirable to assign index outputs to the IDM resulting from a logical "OR" or a logical "AND" being performed over all IDM's for a function index. Both of these special IDM's are produced by the same method. An IDM for the function index in question is chosen. The logical "AND" or logical "OR" of each of the other IDM's (for that function index) and the chosen IDM is then stored in place of the chosen IDM. By successively operating on all of the IDM's required, the end result is the logical "AND" or "OR" of all desired IDM's. The resulting IDM is then available for analysis. In the case of the logical "AND," if it can be output set assigned, the result is an output set assignment which can be used for all IDM's depending on that function index. In the case of the logical "OR," if an output set can be found to fully precedence order, then the function index fully precedence orders. When this is the case this function index should be nested outside those which do not fully precedence order. The output set assignments are made by an algorithm which allows the incidence matrix elements to have weights and then minimizes (or maximizes) the sum of the weights of the output elements (Gupta et al., 1974). In order to fully utilize that capability, weights should be assigned to the IDM elements. If this is not done an arbitrary output set is chosen. The weights should be assigned to the elements in a manner which indicates to the algorithm which are the preferred outputs. One such possible weighting scheme could be to assign weights to reflect the non-linearity of an equation in a certain term. For example, consider the following set of indexed equations: f. = 0 = x +2x. 3(lnx -- )+3x 1-4 4-1 11 1- 11 1 x. il+1 11 For any value of ii, it would be considerably easier to solve eqn. 4-1 for x. or x.1 than for x This can be reflected in the IDM by 11-1 11+1 11 assigning the incidences corresponding to i3 a much higher weight than the incidences corresponding to either i-i or i1+i. The resulting IDM, for a blocking factor of 6, would then appear something like the one in Fig. 4-3 a) (for minimization of the sum of the output weights). The outputs could then be as in Fig. 4-3 b) or c). The assignment of weights is a subjective matter, but it is one area where engineering intuition can be employed to guide the selection of the solution procedure. 4.3 Function-Variable Output Set Assignment The assignment of variable type outputs is made by treating the FVIM as an incidence matrix and using the same algorithm as was used for the index outputs. The aim is to employ the maximum output product criterion developed by Edie (1970). The maximum product criterion is used in an attempt to enhance the convergence properties of the system of equations. It states that the outputs should be assigned in such a way as to maximize the product of the sensitivities of the functions to their chosen outputs. Then sensitivities become weights for the inci- dence matrix elements. The maximum product of the output weights can be achieved by minimizing the negative of the logarithms of the output weights. The real problem is to assign meaningful weights to the elements in the FVIM. What is actually desired is to assign outputs .to the expanded problem, each element having a weight proportional to its corresponding element in the Jacobian. If this were done the restric- tions on FVIM output set assignments would most likely be violated. The procedure adopted for calculating the FVIM weights is the following. The Jacobian is calculated for the expanded set of b) __ 9 9 C It 9G c) I 9 FIGURE 4-3 91 i S9 1 0 90 19 0 INDEX DISPLAY MATRICES WITH WEIGHTS 1 9 1 9 9 1 9 1 1 9 1 9 1 191 equations, with index limits defined by the blocking factors. For the calculation of the Jacobian, estimated variable values are used. The index ordering for the Jacobian has function type nested outermost, as indicated in Fig. 4-4. For each function-variable partition in the Jacobian [Fig. 4-4 a)], an output product is computed by multiplying the elements corresponding to the index outputs for that function in that variable. If a variable type is not a legal output for a function type, the corresponding product is set to a very large positive number. These products are then the weights used in the FVIM output set assignment, as is shown in Fig. 4-4 b). While this procedure does not necessarily generate the true maximum product output set assignment, it does generate the most nearly maximum product output set assignment possible under the restrictions for FVIM outputs. 4.4 Minimum Weighted Tearing Tear variables for a solution procedure with assigned output variables are chosen using the minimum weighted tearing algorithm of Pho and Lapidus (1973). If all variable weights are equal, the tearing scheme chosen will be one with the minimum number of torn variable types. The tear variables, however, are chosen by analysis of the FVIM. The columns in the FVIM do not necessarily all represent the same number of variables. The number of variables of a variable type is determined by the function index definition (not the variable index definitions, which may include decisions). The expected index ranges can be used to calculate the expected number of variables for each variable type. This number can then be used as the weight for the x V z af af af a xl a vT a zT a xT aT a zT a XT avT aZT ifi ,.- ,in x y z ,-4 - y( -I -5 -I ___ l - FIGURE 4-4 FUNCTION-VARIABLE OUTPUT SET ASSIGNMENT b) variable type. The weighted tearing algorithm can then be used to determine the minimum number of actual torn variables directly from the FVIM. The choice of tear variables, like the choice of output variables, is restricted to variable type. Another possible means of assigning weights to variable types is to have the user assign them. If the user has any familiarity with the problem, he may know that some variables are more difficult to guess accurately than others. He could then assign high weights to those difficult to guess. The GENIE system presently allows either of these two methods of assigning weights. An alternative method of assigning tear variables would be to perform a sensitivity analysis on the functions and variables. Then, the tears could be chosen so that the functions would be relatively insensitive to bad guesses for the tear variables. This appears to involve an expensive combinatorial analysis and is not at present implemented. CHAPTER 5 COMPUTER IMPLEMENTATION Since a detailed description of the means of computer implemen- tation of GENIE would require a detailed description of GENDER, the description of the computer implementation given here will be restricted to a discussion of the data structures developed for indexed equations, a discussion of the memory system developed for the GENDER system, and a general flow diagram and description of key subroutines. 5.1 Data Structures In designing data structures for use in the GENDER system the goal has been to use as little space as possible while producing a structure which is compatible with the algorithms which will use it. In order to do this there has been an attempt to minimize the number of list pro- cessing type data structures and to use the addressing capabilities of FORTRAN to minimize searches. In fact, there is only one list of a type which must be searched in the data structures which contain all of the index data for the functions and variables. 5.1.1 Index Data Structures The first major data structure necessary to the analysis of indexed equations is the data structure which contains the information about the function and variable indices. This data structure is the "INDEX" COMMON area. Two of the vectors in this data structure contain definitions of the function indices (FIDV) and of the variable indices 1 (VIDV), where both are vectors of integer variables. The first of these, FIDV for Function Index Definition Vector, is set up as illustrated in Fig. 5-1. The length indicates how many words of the vector are occupied by data. The pointer to VS points to a cell within the vector "VS" which contains the name of the index. The next three entries in an FIDV cell, L., Ui, A. are the lower and upper limits and increment for the index. The range limits are suggested values only since the solution procedures developed are for arbitrary index ranges. The suggested range indicates typical values for the indices in problems which will be solved by the generated solution procedure. The entries L.B and U define the blocking factor and are I 1 the range limits used to produce the Index Display Matrices. The "FPO pointer" is a circular linked list pointer which links together all of the function indices which fully precedence order. A scalar variable "IFPO" serves as a pointer into the circular list so that it can be searched. The last entry in an FIDV cell is a pointer to "IDL," the index decision list vector. If this variable is equal to zero it does not point to anything. The vector IDL will be explained later. The second index data vector is VIDV, the Variable Index Defini- tion Vector, which is illustrated in Fig. 5-2. As in most vectors used in the GENDER system, the first word of VIDV is a length used to indicate how much of VIDV is occupied by data. The size of a VIDV cell is not fixed; however, the first six words always serve the same function. Each of the VIDV cells represents a different variable index. The first word of a VIDV cell is a pointer to a cell within FIDV and thus indicates upon which function index the variable index in question depends. The second word of the cell serves as a pointer to "IDMPV," F I DV Length Pointer to VS Lij Lil Ui Ai1 FPO Pointer Pointer to IDL Pointer to VS Li2 Ui2 Ai2 L 1 FPO Pointer Pointer to IDL - -- 1 etc. FIGURE 5-1 FUNCTION INDEX DEFINITION VECTOR FIDV Cell for i1 FIDV Cell for i2 / VIDV Range L Flag L Offset U Flag U Offset Increment List No. of Offsets First Offset Weight Second Offset Weight etc. FIGURE 5-2 VARIABLE INDEX DEFINITION VECTOR Length Pointer to FIDV Pointer to IDMPV Output Flag Output Offset Range or List Flag Direction Flag Z hng the Index Display Matrix Pointer Vector; for, as mentioned earlier, analysis on variable indices will largely be performed on their Index Display Matrices. The third word serves as a flag which indicates whether or not the variable index has been output set assigned. If it has, the next word indicates what the output is if it is a simple offset from the function index. If the output is not a simple offset the Index Display Matrix must be examined to determine the index output. The fifth word of the VIDV cell is a flag indicating whether or not the variable index is defined by a range or a list. The sixth word indicates whether the index is to be incremented in ascending or descending order. The rest of the VIDV cell depends on whether the index is defined by a range or a list. If the index is defined by a range, the length of the VIDV cell is fixed, as there are five more entries for a range defined index. The first entry indicates whether the lower limit of the range is offset from the lower limit of the function index or from the function index itself. The next entry indicates what that offset is. The next two entries serve the same purpose for the variable index range upper limit. The last entry is the range increment. If the index is defined by a list, the length of the VIDV cell is variable. This is because the number of list elements is variable. Within the list definition section of the VIDV cell, the first entry indicates how many list elements there are. Subsequently each list element occupies two entries, the first of these being the offset from the function index value and the second being a weight which is used by some of the analysis routines. To facilitate analyses performed on function and variable indices and to prevent searches, there are two vectors of pointers, one associated with FIDV and the other with VIDV. The form of both vectors (VIPV and FIPV) is the same, and they will be described together. As usual the first word is the length of the data in the vector. The second word is a pointer to the first FIDV or VIDV cell, the third word is a pointer to the second cell, etc. The entire data structure as described so far fits together in the fashion illustrated in Fig. 5-3. Satellite data structures to the index representation data structures are the vectors "IDMPV," the Index Display Matrix Pointer Vector, and "IDL," the Index Decision List. The vector IDMPV is accessed through VIDV and merely provides the name of the Index Display Matrix for each of the variable indices in VIDV. The name is not stored in VIDV because that would complicate searching through the Index Display Matrix should that be necessary. The index decision list is the vector which contains the index decision declarations. Index decisions are offsets either from the lower or upper limit of the function index. Each function index has associated with it an IDL cell which is composed of two sub-cells, as shown in Fig. 5-4. If a function index has no index decisions, its IDL cell will have no entries. Each IDL cell has three length specifications, the first stating the length of the entire cell. The other specifies the length of each of the sub- cells. The first sub-cell, the L offset sub-cell, specifies the offsets of decisions from the function index lower limit, and the second specifies the offsets of decisions from the function index upper limit. The index decisions can be declared in two ways; either the user can declare them or algorithms can assign them. Declared decision variables can never be changed, while derived decisions can be changed. In order VIPV etc. FIGURE 5-3 INDEX DATA STRUCTURE VIDV ID L Length I Length Length L Offset 1 L Offset 2 etc. Length U Offset 1 U Offset 2 Setc. etc. L Offset Sub-cell U Offset Sub-cell FIGURE 5-4 INDEX DECISION LIST IDL Cell to distinguish between these two types of decisions, derived decision offsets are incremented by 10,000. Since decision offsets will never be that large, large offsets are easily recognized as derived and are adjusted to their true value when necessary. 5.1.2 Function-Variable Incidence Matrix Representation It is necessary that the Function-Variable Incidence Matrix be stored in a data structure which allows for easy access by the algo- rithms. Since the FVIM is a type of incidence matrix it is reasonable to store it in the data structure used to store incidence matrices. The incidence matrices are stored in a data structure called the SIM data structure (Cunningham, 1973). Appendix A contains a description of the SIM data structure. The principal difference between the normal incidence matrix and the FVIM is the fact that an element in the FVIM is actually a list of variable indices. To allow for that it is necessary to modify the element entry in the normal SIM structure to accommodate an extra word. This is easily done since the length of an entry is variable in the SIM data structure. This extra word serves as a pointer to the vector "VILIST." Thus, each FVIM entry has associated with it a pointer to VILIST. The structure of VILIST is illustrated in Fig. 5-5. Each FVIM element pointer points to a different VILIST cell. The VILIST cells are merely lists of pointers to variable index definition cells. This results in each element in an FVIM having associated with it a list of variable indices, which is what is required. 5.2 A Versatile Memory System A part of the design philosophy of the GENDER system has been to construct the different data structures necessary in separate labeled COMM0N areas. For example, a special data structure was developed for VILIST Length Length Pointer to VIDV Pointer to VIDV Pointer to VIDV Length Pointer to VIDV Pointer to VIDV etc. VILIST Cell FIGURE 5-5 VARIABLE INDEX LIST VECTOR the storage of sparse incidence matrices. This data structure is built into a set of labeled COMMON areas called "SIM COMMON." At any time only one sparse incidence matrix (SIM) can occupt this common area. Quite often during the analysis performed by the GENIE system it becomes necessary to temporarily stop working with one SIM and start working with another, without losing the data in the first. Since all routines operate on SIM's in the SIM COMMON area, it is necessary to move the first SIM out of common and at some later time restore it to the common area. This type of movement is not restricted to SIM's and is desirable for many different data structures. It became necessary to design a general purpose memory system, easily assessible to any of a variety of routines which need to store data in a relocatable fashion, so that the data can be sestored to the appropriate data structure when required. The memory system which accomplishes this is itself a labeled COMMON area called, appropriately enough, "MEMORY." The memory system consists of two major storage areas, or vectors, the larger being "MEMORY" and the smaller being "MEMDIR." MEMORY is the area of core into which data to be saved are moved. It is a long vector (in the first version 5,000 words) and could conceivably be much longer. MEMDIR is a directory to the vector MEMORY which contains all of the informa- tion necessary to relocate a data item stored in MEMORY. A third major storage area associated with the memory system is a random access mass memory data set onto which data from MEMORY are written should MEMORY ever not have enough space available to store a data item. It should be noted that many different data items will in general be stored in MEMORY at the same time, and that if enough are stored MEMORY will become full. It is perhaps easiest to consider ME: !,' to be a buffer for the mass memory device. The memory system then can be considered to be a mass memory system using relatively slow speed mass memory with a higher speed core resident beffer. Access to stored data is made through the buffer. The memory system handles all access to the mass memroy. Although to the user of the memory system it appears that all data are kept in MEMORY, an understanding of the interaction between MEMORY and mass memory helps in understanding the system. The random access data set has associated with it a logical record length (1,000 words in the case of the memory system). The data set is divided up into a number of these logical records. In this particular implementation, written for IBM 370/165, the random access data set has data written to or read from a particular record on a disk. This makes it possible to access a particular record. Since the records are numbered sequen- tially, the words in the first record can be thought of as being words 1 to 1,000 and the words in the second record as 1,001 to 2,000, etc. Figure 5-6 illustrates this addressing scheme. If the address of a word of data is known, then the record of the data set it is stored in can easily be determined. In the memory system the first record is reserved for a special function which will be described later and is never used for data storage. Data storage begins with the second record. Data which are moved to MEMORY have associated with them absolute addresses, which are the word addresses that the data would have were they stored onto the random access data set. This can most easily be visualized by imagining that the vector MEMORY is superimposed on the random access file in Fig. 5-6. The first word of MEMORY would have absolute address 1,001, the RANDOM ACCESS DATA SET Word 1 Word 2 + Word 1000 - Word 1001 + Word 1002 Word 2000 FIGURE 5-6 STRUCTURE OF RANDOM ASSECC DATA SET FOR MEMORY SYSTEM Record 1 Record 2 I e. etc. second word 1,002, etc. The vector ME i RY has two states, one indi- cating that the data in MEMORY are not actually stored in the random access locations corresponding to their absolute addresses, the other indicating that the data in MEMORY are also in random access memory. Until data in MEMORY have been written to the random access file, the first word of MEMORY will correspond to word 1,001 of the random access file. Suppose, however, MEMORY is 5,000 words long and that 4,500 words of MEMORY are being used to store data and a request to store a data item of length 1,000 words is made. Since there are only 500 words available for data, some of the data in MEMORY will be written to the random access file. Actually all of MEMORY is written to direct access. Since data are accessed by 1,000 word records, the first 4,000 words of MEMORY are considered as written and the next 500 words are moved to the first 500 locations in MEMORY. The first word of MEMORY then corresponds to absolute location 5,001. This would correspond to superimposing MEMORY over the random access data set beginning with the sixth record. In this manner considerably more words than can be stored in MEMORY can be stored and addressed. When a data item is stored into MEMORY it has a directory entry created for it so that it can be referenced through the directory. A directory entry consists of six words, allocated to the following purposes: WORD CONTENTS 1 Number of Characters in Name 2 First 3 Characters in Name 3 Second 3 Characters in Name |