User Guide - Computing and Software
Contents
1. Figure 5 2 Structure of illPosed2 In this example the missing equations and variables are made obvious by the shaded entries red in color print However for a larger example we may not be able to easily distinguish missing equations and variables A call of getMissingEqnsVars yields gt gt meqn mvar getMissingEqnsVars sadata meqn mvar Oo i w I Finally we show how DAESA displays several diagnostic blocks for a SIP DAE We use illPosed3 from the examples directory see Figure function f illPosed3 t z x1 x4 f 1 2 f 3 f 4 f 5 f 6 end z 1 x2 z 2 x3 z 4 x5 z 5 x6 x1 x2 x1 3 x2 x1 x2 x3 x4 z 3 z 6 xb xb xb x573 x5 xb x6 5 x6 E x6 R x674 x673 x6 Figure 5 3 MATLAB function for evaluating illPosed3 example Calling showStruct with disptype set as blocks yields the diagnostic block structure seen in Figure 5 3 Chemical Akzo Nobel 35 Indices of Equations ILLPOSED3 Diagnostic BTF Structurally ill posed Shaded under and over determined Indices of Variables Figure 5 4 Diagnostic block triangularized structure of illPosed3 Here there is a structurally under determined set of variables x4 and x3 in equation 0 fs a structurally over determined set of equations 0 f3 fa fe in variables 15 and zg and a structurally well determined system in variables x and2. L c i where the state variables of the ith pendulum i gt 1 are z y and 4 G gt 0 is gravity L gt 0 is the length of the first pendulum and c is a constant For i gt 2 the length of the ith pendulum depends on A A system of p pendula has size n 3p and index 2p 1 2 5 The function for evaluating is in Figure Here the number of pendula p is passed as a parameter In the for loop we evaluate the equations for pendulum 2 p The daeSA function can be called as p 5 n 3 p sadata daeSA multiplependula n p In Figure 5 10 we show the original and fine block triangularizations of this problem the coarse and fine triangularizations are the same 5 4 Multiple pendula 38 function f multiplependula t x p G 9 8 L 1 0 1 first pendulum Cc f 1 Dif x 1 2 x 1 x 3 4 0 A121 f 2 Dif x 2 2 x 2 x 3 G 4 0 y A1y1 40 ue 2 f 3 x 1 72 x 2 72 L 2 FS kb 4 pendulum gt 1 for i 2 p xi 3 i 2 yi 3 i 1 4 x y indices li 3 i lil 3 i 3 4 Ai Aj 1 indices f xi Dif x xi 2 x xi x li tisa A f yi Dif x yi 2 x yi x li G 40 y A y G f li x xi 72 x yi 2 L cxx 1i1 72 y lo E E y L cX a end end Figure 5 9 DAESA function for evaluating the multiple pendula problem MULTIPLEPENDULA Size 15 structural Index 11 DOF 10 MULTIPLEPENDULA Fine BTF Shaded structural n
3. 6 6 4 2 3 0 cl 0 0 2 0 0 0 dl 2 2 0 0 3 0 which when compared to the results of Figure should give the reader a better understanding of these functions The local offsets are returned in the original equation order and not in the permuted block triangular equation order shown by calling showStruct with disptype set as fineblocks If we wanted to change the ordering the commands using getBTF detailed below should be used We may have a large solution scheme produced from the printInitData and printConstr func tions of 2 3 3 and as such may want a more compact representation of the data To do this we 31 5 1 Well posed DAE example 32 use getInitData and getConstr to produce gt gt iv getInitData sadata iv 2 2 0 1 4 0 gt gt constr getConstr sadata constr 4 4 6 0 1 3 so for example variable one and its derivative are needed initial conditions for the solution process as already shown with printInitData in 82 3 3 Finally we consider the result of applying getBTF and getQLdata to the modified2pendula example pe pv cb fb getBTF sadata pe 5 4 6 3 2 1 pv 5 6 4 1 2 3 cb 1 4 7 fb 1 2 3 4 7 gt gt cql fql getQLdata sadata cql 0 1 fql 0 1 0 1 Recalling Figure in particular the blocks in the coarse and fine cases should give the user a good understanding of the outputs of these functions Finally if we wanted to produce now the local offsets in the BTF or
4. The reason for the extra components in the non quasilinear case is as follows Consider a particular independent variable value t If a set of values X in is consistent with some solution of the DAE at t then in the linear case that solution is unique In the nonlinear case there may be several solutions consistent with these values however including the next level of derivatives in X restores local uniqueness Chapter 4 DAESA functions DAESA exploits MATLAB s operator overloading to process the DAE given by a user supplied function for evaluating the f in 1 1 In particular this is the method DAESA uses to extract the signature matrix and determines for each equation if it is quasilinear in the leading derivatives After the signature matrix is constructed DAESA finds out if the problem is structurally well posed and if so calculates the offsets of the problem and then determines structural index and DOF Since it knows the structure of the analyzed DAE DAESA constructs a block triangular form of it finds local offsets and determines block by block quasilinearity Based on the offsets and linearity information DAESA deduces which variables and derivatives of them need to be initialized and what the constraints are The SA is performed by the function daeSA It returns an object of the class SAdata which encapsulates all the data obtained from the SA Each of the remaining DAESA functions takes an object of this class as a p
5. fineblocks When called on sadata returned by daeSA above these functions produce the plots in Figure 2 2 MODIFIED2PENDULA Size 6 structural Index 7 DOF 5 MODIFIED2PENDULA Fine BTF Shaded structural nonzeros in system Jacobian J MODIFIED2PENDULA Coarse BTF Size 6 structural index 7 DOF 5 Boxed HVT Size 6 structural index 7 DOF 5 Shaded structural nonzeros in system Jacobian J Shaded structural nonzeros in system Jacobian J Boxed positions that contribute to det J Boxed positions that contribute to det J Indices of Variables Indices of Variables 1 2 3 4 5 6 c Indices of Variables 5 6 P 1 2 3 c 5 6 4 1 2 3 amp aids 1 o 4 s B I n 5 B B 0 0 o o i 2 an ea 0 B 0 Eu B WS n c B jj t E Ij E g3 p A 6 FED a al gs o pl Boe 5 5 5 g A 9 ME B E 6 g b p 2 6 9 E 9 25 e Bo FE B ols 22 a oz 6 E a o 2 1 E Dl 1 2 o 4 d 6 6 4 2 8 O d3 0 2 6 6 4 d 3 0 0 2 2 0 d 3 0 2 6 6 4 a original structure b coarse block triangularized structure c fine block triangularized structure Figure 2 2 Structure of 2 1 and its block triangularizations In Figure 2 2 the oo entries are not displayed The top of the figures indicate index and DOF The shaded entries green and yellow in color print are where the system Jacobian 83 1 J is structurally nonzero a The boxed entries denote a highest value transversal HVT e and c the boxed entries denote the ent
6. 2 3 Extract structural analysis data 7 That is the constraints are A 0 CI lt 3 fh E 0 C2 lt 3 f 0 C3 lt 0 Te 0 C6 lt 2 The default equation names f1 2 can be changed like the variable names in 82 3 3 using the key fcnnames instead of varnames Also as in 82 3 3 omitting the key pair outfile and gt constr txt will produce an output in the command window 2 3 5 Solution scheme A solution scheme is reported by printSolScheme It reports how to solve the DAE what variables to initialize what equations to solve and for which variables by stages using its fine block triangularized structure In particular it shows how to compute values for the derivatives of each variable Full solution scheme The call printSolScheme sadata varnames vars gt outfile solscheme txt detail full with vars set as above produces file solscheme txt with Solution scheme for modified2pendula problem Initialization summary x x y var u V v v y STAGE k 6 1 block Block 4 6 Solve nonlinear equation give trial values O f3 for x y STAGE k 5 1 block Block 4 6 Using x y Solve linear equation give trial values 0 f3 for x y STAGE k 4 1 block Block 4 6 Using x x y y Solve linear 3x3 system O fi f2 3 for x y lam 2 3 Extract structural analysis data STAGE k 3 2 bloc
7. alently Taylor coefficients the reason for the raised offsets of pendulum 1 is clear In the absence of coupling the natural scheme is at stage k 2 find x y and u v at k 1 find x y and u v at k 0 find x y A and vw v u and so on However u v must satisfy F 0 With the coupling F involves x which in the above scheme has not been found yet Similarly u v must satisfy F 0 which involves x which has not been found yet and so on This is cured by shifting pendulum 1 back one stage so that the scheme becomes k find then find 3 L Y 2 y u v 1 a y uy MW IN H ui 0O x y A 0 and so on Hence provided pendulum 1 s equations are solved before those of pendulum 2 at each k stage pendulums 1 and 2 can be solved as separate size 3 systems but each derivative of z is available just when needed by pendulum 2 Because of the shift pendulum 1 behaves as if the overall stages k 3 2 1 are its local stages k k 1 2 1 0 associated with local offsets amp 0 0 2 and d 2 2 0 which come from analyzing it as a stand alone system This shows that the relation between the minimal initial values problem and the sequencing of a solution scheme involves the DAE s block lower triangular structure Hence we can take advantage of a DAEs block triangularization by considering local offsets and local solution stages in each block to red
8. 3 N S NEDIALKOV AND J D PRYCE Solving differential algebraic equations by Taylor series 1 Computing Taylor coefficients BIT 45 2005 pp 561 591 4 Solving differential algebraic equations by Taylor series II Computing the System Jaco bian BIT 47 2007 pp 121 135 5 Solving differential algebraic equations by Taylor series 111 The DAETS code JNAIAM 3 2008 pp 61 80 ISSN 1790A R8140 6 DAETS user guide tech rep Department of Computing and Software McMaster Uni versity Hamilton Ontario Canada L8S 4K1 2008 2009 7 C C PANTELIDES The consistent initialization of differential algebraic systems SIAM J Sci Stat Comput 9 1988 pp 213 231 8 J PRYCE N S NEDIALKOV AND G TAN DAESA a Matlab tool for structural analysis of DAEs Theory 2013 Accepted for publication in ACM TOMS 9 J D PRYCE A simple structural analysis method for DAEs BIT 41 2001 pp 364 394 10 G REissIG W S MARTINSON AND P I BARTON Differential algebraic equations of index 1 may have an arbitrarily high structural index SIAM J Sci Comput 21 1999 pp 1987 1990 42
9. Before using DAESA one should set paths to the above directories This can be done by executing startup m in directory DAESA or by setting them manually through MATLAB s GUI Example programs The results in this user guide are produced by the following scripts in the examples directory file used in DAE function in examples DAEs sa modified2pendula m 5 1 modified2pendula m sa illPosed m illPosed1 m illPosed2 m illPosed3 m sa chemakzo m akzonobel m sa multiplependula m multiplependula m The remaining files in the examples and examples DAEs directories are used in 2 Once DAESA is on MATLAB s path it can be used within MATLAB from an arbitrary working directory If this is done whenever an outfile is specified in the print functions of 82 3 3 the file path if not an absolute path will be relative to the current working directory 39 Appendix A Supported standard functions The following two tables list all standard functions currently supported by DAESA DAESA name resulting function DAESA name resulting function n power where the power is sinh hyperbolic sine a real number cosh hyperbolic cosine exp exponential tanh hyperbolic tangent log natural logarithm sech hyperbolic secant log10 base 10 logarithm csch hyperbolic cosecant sqrt square root coth hyperbolic cotangent sin sine asinh hyperbolic arcsine cos cosine acosh hyperbolic arccosine tan tangent atanh hyperbolic arctangent
10. differentiation index on some problems see e g Al 9 We would have a certificate that the SA is successful if the system Jacobian 83 1 is non singular at a consistent point see also 4 The present tool does not compute consistent points and does not evaluate the system Jacobian it performs symbolic type analysis of DAEs We plan to incorporate the evaluation of this Jacobian in a next version of DAESA This user guide is organized as follows Chapter 2 provides a quick start to using DAESA we show basic SA with it Chapter 3 gives an overview of the theory behind DAESA Chapter 4 describes the functions of DAESA Chapter 5 presents several examples of analyzing DAEs How to obtain and install DAESA is described in Chapter 6 Chapter 2 Quick start In this chapter we illustrate how DAESA performs basic SA more examples are in g5 Using DAESA involves three steps 1 specifying a DAE 82 1 2 calling the main SA function 82 2 and 3 extracting structural data 82 3 2 1 Specify a DAE A DAE is given by a function with arguments and return value as follows function f daefcn t x options Here t and x are the independent and dependent variables respectively and options is an optional list of one or more arguments x is a column n vector where x j is variable x and f is a column n vector with f i containing the evaluation of fi We apply our tool on the artificially modified two pendulum pr
11. is X annotated with the offsets c dj the names of the functions and variables and the positions of a HVT The oo entries are left blank for readability since in larger systems typically almost all entries are oo Also these are forbidden positions since an HVT of a well posed DAE cannot have an entry in a oo position Also shown in the tableau is the number F of degrees of freedom DOF given by F Val E Y d Soa It is the number of independent initial values IVs it requires which is the same as the maximum number of IVs that can be specified fixed This is illustrated below for the simple pendulum example an index 3 DAE which has two HVTs marked and Example 1 Simple pendulum f x g y yrA G 3 2 heu E The signature tableau is T y ci 2 0 1 0 gt J o 0 0 2 f g h 2 2 0 DOF 2 d S For n up to about 8 setting up X and finding a HVT by eye is usually easy Otherwise let DAESA do the analysis as suggested in 4 The n x n System Jacobian matrix is formed as c o i y a s DUET id m if d Ci dij J or equivalently J On 3 3 a Cu d 1 E 0 otherwise Equivalence of the equation 3 3 is given in 9 By f we mean df dt treating the variables and their derivatives as unknown functions of t For instance if f x 1 13 then fi x 2 13 2 2 similarly for higher derivatives The x and derivatives thereof are tre
12. or t e Currently DAESA supports only scalar operations Array operations will be implemented in a future version e In function daefcn the indices of x and f are bounded from 1 to n Each variable x j should be used at least once and each equation f i should be evaluated where i j 1 2 n As an example of such a function see Figure 4 2 Structural analysis 21 4 2 Structural analysis The structural analysis is performed by the daeSA function function sadata daeSA daefcn n options daefcn Input a function handle for evaluating the DAE see p Input positive integer problem size options Input optional list of one or more arguments that are passed to the definition of the DAE sadata Output an SAdata class object This function performs the structural analysis and encapsulates all the data in a return SAdata object sadata These SA data can be accessed by the functions described in daeSA can fail in the following cases 1 Problem size n in daeSA is not consistent with the number of variables or and equations encoded in daefcn This can happen if a variable or equation of index not in 1 n is accessed b a variable is not used and or an equation is not evaluated 2 daefcn contains functions or operations that are currently not supported 3 The DAE system defined by daefcn is over or underdetermined As long as all the given variable and equation indices are in 1 n the structure of the D
13. posed DAE Example ss mias as E EU hope mcd a Se ed 5 2 Ill posed DAE examples LL RR 5 3 Chemical Akzo Nobel vica ok he ew LEO RO ee ede unc om es 5 4 Multiple pendula 12 12 14 15 16 17 19 20 21 21 27 29 CONTENTS iii 6 Installation 39 A Supported standard functions 40 41 Bibliography 42 In the electronic version of this document every cross reference is a hyperlink For instance you can click on the entry Quick start above to jump to that section This also applies to page numbers in the Index and in the body of the text to chapter and section references and to equation numbers To return to where you just were use your PDF readers Back command List of Figures 21 DAESA function for evaluating 2 1 ese A 4 a oro DEI 5 eee vans 33 6 pe Re ee A T E ees 34 MET 34 5 4 Diagnostic block triangularized structure of illPosedB3 35 5 5 DAESA function for evaluating the chemical Akzo Nobel DAE 35 Se ROSE EBENE mE E WIN UU NR NIE 36 5 7 Akzo Nobel solution scheme 36 He 36 Verve cM MM 38 5 10 Structure of 5 1 and its fine block triangularization so 38 lv List of Tables a cala ii 10 ee ee be ee bee bee eee 15 a pang pa Be 37 a a a ps ia 40 Chapter 1 DAESA Overview DAESA Differential Algebraic Equations Structural Analyzer is a MATLAB tool for structural anal ysis SA of differential algebraic equatio
14. report the solution scheme in compact form or full which results in a more detailed output See g2 3 5 for more information on the detail option In the three functions above the output file fname can be specified as a full path e g c My Documents My_DAESA_outputs my_files_name txt If no full path is specified then the path is relative to the current MATLAB working directory For example Reports filename txt will print to a Reports folder in the current directory For all the above print functions the format is plain text Data output samples using these functions can be found in Chapter 5 Examples We now work through examples of some of the functions of the previous chapter We start with the SWP modified2pendula example 85 1 and then we show how DAESA visualizes problems that are structurally ill posed 85 2 Then we present results from applying DAESA on the chemical Akzo Nobel problem 1 85 3 and a problem consisting of several coupled pendula 85 4 In all following examples it will be assumed that sadata refers to the result of calling daeSA on the daefcn being explained 5 1 Well posed DAE example First we use the SWP case example Recall the result of showStruct on the modified2pendula example displayed in Figure We now apply the functions getHVT and getOffsets to this example Returning HVT getHVT sadata HVT 3 2 1 6 5 4 gt gt c d cl dl getOffsets sadata 4 4 6 0 0 2 d
15. y2 f1 yt Figure 5 8 Fine block triangularized structure of the chemical Akzo Nobel problem Figure 5 7 Akzo Nobel solution scheme AKZONOBEL Fine BTF Size 6 structural index 1 DOF 5 Shaded structural nonzeros in system Jacobian J Boxed positions that contribute to de Indices of Equations A w N a HD o Indices of Variables 1 2 3 4 5 6 fl 0 0 o 0 o fi 0 B 0 0 fl o 0 0 o fl o o o o 0 H B 0 0 lo 1 1 1 1 4d 0 1 1 1 1 4d 0 J Gi Gi 0 0 5 4 Multiple pendula 37 This scheme suggests that at stage k 1 we should give initial values for y for all 1 5 For stages k gt 0 we proceed as shown in Table block solve for 66 SP ww 5 5 po ge dd jp we 3 3 n yo 9 9 jo yor 1d je Jen Table 5 1 Solution scheme for the Akzo Nobel problem for k gt 0 In I this problem is classified as nonlinear and initial values for all variables and their first derivatives are given Our solution scheme reveals six linear equations if solved in a block wise manner and reports that only yi ys need to be initialized 5 4 Multiple pendula In this example we illustrate using a for loop in the function for evaluating a DAE system where the number of iterations is passed as a parameter Consider the chain of pendula first pendulum ith pendulum 0 27 0 x Aiti a 7 5 1 0 yitAu G 0 y Ay G 0 224 y2 LI 0 z2 y
16. AE even if ill posed can be shown by showStruct see p 4 3 Obtaining structural data If the analyzed DAE is structurally well posed we write SWP otherwise if it is structurally ill posed we write SIP An SAdata class object returned by daeSA contains data from which one can extract the following information about the DAE 1 swp a logical value that indicates if the DAE is structurally well posed 2 meqn mvar indices of missing equations and variables respectively 3 index structural index 4 3 Obtaining structural data 22 4 dof number of degrees of freedom 5 S signature matrix 6 HVT position of a highest value transversal T iv set of variables to be initialized 8 constr set of constraints 9 c d c1 d1 global and local offsets respectively 10 pe pv permutation row vectors of 1 n for the triangularized signature matrix 11 cb fb boundaries of the diagonal blocks of the coarse fine block triangularization respectively 12 cql fq1 quasi linearity data of diagonal blocks in the coarse and fine block triangularization respectively 13 fhandle DAE definition function handle 14 n problem size Each of the remaining functions take as input the object obtained for daeSA which we name sadata function swp isSWP sadata swp Output 1 if SWP and 0 otherwise function megn mvar getMissingEgnsVars sadata meqn mvar Output the indices of equations in row vector meqn and indices of varia
17. DAE function name 4 5 Data output 29 structurally ill posed indices of variables and equations e When disptype is blocks DAE function name structurally ill posed indices of variables and equations missing equations variables a diagnostic BTF as explained in the description of getBTF on p 25 and shown in Figure 5 1 in e When disptype is fineblocks the user will receive a message printed to the command window telling them the DAE is SIP and fine blocks cannot be displayed daefcn is structurally ill posed Fine blocks cannot be displayed 4 5 Data output For the following functions we assume SWP If SIP an error message is printed to the command window daefcn is structurally ill posed function printInitData sadata options This function reports the variables and derivatives of them that need to be initialized for solving the DAE This function can have in options one or two key value pairs where the keys are varnames and outfile varnames vars By default the variable names are x1 x2 xn If each of the variable names needs to be replaced by the same variable name say v followed by index number one can pass varnames v Alternatively if one would like to specify each variable varnames cavarnames Different variable names can be given in a cell array of n strings say cavarnames and then passed as varnames cavarnames T
18. DAESA User Guide Version 1 0 Ross McKenzie John D Pryce Cardiff University United Kingdom Guangning Tan Nedialko S Nedialkov McMaster University Canada Copyright 2013 R McKenzie N S Nedialkov J D Pryce and G Tan Preface DAESA Differential Algebraic Equations Structural Analyzer is a MATLAB tool for structural analy sis SA of differential algebraic equations DAEs It allows convenient translation of a DAE system into MATLAB and provides a small set of easy to use functions DAESA can analyze systems that are fully nonlinear high index and of any order It determines structural index degrees of freedom constraints variables to be initialized and suggests a solution scheme The structure of a DAE can be readily visualized by this tool It can also construct a block triangular form of the DAE which can be exploited to solve it efficiently in a block wise manner This code was originated by JDP in 2001 2003 Ning Liu M Sc 2006 Computing and Soft ware McMaster added new functionality and in particular the computation of signature matrices through operator overloading GT did a major rewrite and added many new features to DAESA Ian Washington Ph D candidate Chemical Engineering McMaster was the first user of DAESA he applied it to study chemical engineering problems and suggested many improvements The figure on the title page is from DAESA applied to one of these problems
19. We are hoping that this tool will be helpful to researchers and practitioners when studying systems containing differential and algebraic equations It could also be applied to pure differential and pure algebraic systems For bug reporting questions and suggestions on improving DAESA please contact the authors at daesatool gmail com We acknowledge with thanks the support given to RM and JDP by the Leverhulme Trust and the Engineering and Physical Sciences Research Council both of the UK and to GT and NSN by the Canadian Natural Sciences and Engineering Research Council and the McMaster Centre for Software Certification July 9 2013 Contents 1 DAESA Overview 2 Quick start 24 pera DAE ssc ed ecek perata p Hea A web kd Ron ES Rom Ed de d us pe Se eee ie uud Pee ee See E m AE ee ee ee a 23 1 OSA zai zou Ae NN 23 2 Index ang DOF 3 2 626224 eo be we EUR RON Ea PAGS BEES 2 00 Imibializationi summary kox4 amp xx x 4OECYUE 2 x x ceded Ewes Ewe x 23 4 Constraints 2 3 5 Solution schemel oos 3 Theory Overview epa tds y a dig dea a Se Spee Hida pw do we au Oe ar e beso cross A A ala aa a A 3 4 Advantage of block triangularization iL 2 2 sala usos da Rs CA BEG ARA Sc A ES e Re Siue RES 4 DAESA functions 4 1 Specifying a DAE 4 2 Structural analysis 4 3 Obtaining structural data amp lt 4 234 4 4 2 Emm He She na dE 4 4 Visualization 4 5 Data output 5 Examples 5 1 Well
20. arameter and extracts from it the data it needs In this chapter we present all the functions in DAESA In we show how a DAE is specified in a function In we introduce the main function daeSA that performs SA In 4 3 we describe all the functions and how to obtain structural data from the analysis performed 19 4 1 Specifying a DAE 20 4 1 Specifying a DAE function f daefcn t x options t Input independent variable t Input column n vector dependent variables where x j denotes variable x for j 1 n options Input optional list of one or more arguments that are passed to the definition of the DAE useful for parameters Output column n vector equations where f i contains the evaluation of f in 1 1 The code of daefcn may contain e unary operations e binary arithmetic operations e standard functions sin cos tan sinh cosh tanh sqrt exp log where denotes the real raised to a real power function For a complete list of supported functions see Appendix A e Dif var k the differentiation operator d dt that returns the kth derivative of var For instance Dif x 3 2 represents 14 d 23 dt Notes e The code should not contain branches that is if or switch statements Functions that are defined using these are not covered by the current theory e The code may freely introduce other variables that depend directly or indirectly on the inputs x and
21. ated 3 2 The solution method 14 as unrelated independent variables within f For instance if fi is z125x1 then O0 fi dx Op OE z rh If J thus defined is not identically singular the SA based method almost certainly succeeds To be precise if there is a consistent point i e a point that satisfies all the equations f 0 where J is nonsingular a solution to the DAE exists through that point it is locally analytic DAESA can permute the system Jacobian J to upper triangularized block form and identify the structural non zeros in it For the pendulum 3 3 gives the system Jacobian Of Ox 0 Of OX L dd J 0 0g 0y 8g 0OA 0 1 y Oh Ox Oh Oy O0 2r 2y 0 This is not identically singular since det J 2 1 y Indeed at a consistent point from the third equation of 3 2 det J 2L 0 5 An upper bound for the standard differentiation index v4 is given by the structural index 1 if some d is zero Vg max cj E i 0 otherwise For the commonest kinds of DAE vs v4 but the difference can be arbitrarily large 10 3 2 The solution method After finding the offsets and HVT for the DAE as described in the SA method solves the system in stages numbered as k with fen for all i such that k c gt 0 3 4 to solve for variables k d ai ta for all j such that k d gt 0 3 5 We start at a negative valued stage equal to max dj and at stage 0 we are solving an n x n sy
22. bles in row vector mvar that are missing in the DAE definition If no equation is missing meqn and if no variable is missing mvar function index getIndex sadata index Output structural index if SWP and NaN otherwise 4 3 Obtaining structural data 23 function dof getDOF sadata DOF Output DOF if SWP and NaN otherwise function S getSigma sadata S Output signature matrix in dense form function HVT getHVT sadata HVT Output if SWP a row n vector such that i HVT i i 1 n forms a HVT in the unpermuted signature matrix For example for the MOD2PEND problem 82 1 this function returns the vector 3 2 1 6 5 4 cf the HVT in Figure 2 2 a If SIP HVT is a row n vector of NaN s function c d getOffsets sadata function c d cl dl getOffsets sadata c d cl dl Output If SWP c and d are row n vectors containing the global equation and variable offsets respectively cl and dl are row n vectors containing the local equation and variable offsets respectively If SIP all output row vectors are Note The number of output arguments should be either 2 or 4 if not an error message will be printed 4 3 Obtaining structural data 24 function iv getInitData sadata iv Output If SWP iv is a non negative integer row n vector where iv j is the number of derivatives of variable j that need initial values when solving the DAE As a result the sum of
23. dering given by calling showStruct with disptype set as fineblocks then we would use the following commands gt gt cl cl pe cl 0 0 0 2 0 0 gt gt dl dl pv dl 3 0 0 2 2 0 5 2 IlI posed DAE examples 33 5 2 Ill posed DAE examples We now turn our attention to the SIP problem resulting from removing equation three in the modified2pendula example mentioned above By applying showStruct we obtain the following figures ILLPOSED1 ILLPOSED 1 Diagnostic BTF Structurally ill posed Structurally ill posed Missing equation variable Indices of Variables Indices of Variables 1 2 3 4 5 6 6 1 2 3 4 5 1 2 0 1 B D o ca 2 0 52 BE T w NN G8 RN d Id go B HH 5 5 o 4 2 0 9 4 A E 2 2 B5 3 0 25 D E 6 2 0 0 3 a original structure b diagnostic block triangularized structure Figure 5 1 Structure of illPosed1 and its diagnostic block triangularization If we apply getBTF to this example we obtain gt gt pe pv cb fb getBTF sadata pe 1 2 6 4 5 3 pv 6 1 2 3 4 5 cb fb 1 6 6 6 7 1 2 7 7 7 see also p 25 Now to show how DAESA displays missing equations and variables we consider illPosed2 which is the modified2pendula example omitting both equation three and variable y A call of showStruct yields Figure 5 2 Ill posed DAE examples 34 ILLPOSED2 Structurally ill posed Indices of Variables Indices of Equations mm EH
24. e cql i false and this block is non quasilinear Similarly for the fine block triangularization if fql i true then the ith diago nal block comprising equations of indices pe fb i fb i 1 1 in variables of indices pv fb i fb i 1 1 is quasilinear in the leading derivatives of these variables otherwise fql i false and this block is non quasilinear Assume SIP All output row vectors are function h getDAEfhandle sadata h Output returns a handle of the function passed to daeSA function n getSize sadata n Assume SWP Output returns problem size the number of equations and variables Assume SIP Output if the problem is missing variables or equations the result is still the same the function returns the number of variables and equations specified in the initial daeSA function call 4 4 Visualization 27 4 4 Visualization function showStruct sadata options Assume SWP This function displays the structure of the signature matrix in the current figure or in a new figure if none is open showStruct can have in options one to three optional key value pair arguments explained as follows gt disptype typeValue typeValue is a string that indicates how the structure the DAE should be displayed typeValue displays original default original structure of the DAE as given showStruct sadata disptype original is the same as showStruct sadata gt blocks coarse block
25. his results in the th variable name being replaced by cavarnamesti outfile fname By default the text is printed on the screen If optional outfile fname is passed the output is stored in a file with name fname 4 5 Data output 30 function printConstr sadata options This function reports the functions and derivatives of them that form the set of the constraints This function can have in options one or two key value pairs where the keys are fcnnames and outfile fcnnames fcnname fcnnames cafcnnames By default the function names are f1 f2 fn If function names are different from the default ones they can be replaced in the output in the same way as in printInitData For instance fcnnames g replaces each of the function names with g followed by index number fcnnames cafcnnames replaces the default function names with different function names where cafcnnames is a cell array of n strings outfile fname The same as above in printInitData function printSolScheme sadata options This function reports a solution scheme for the DAE options is an optional list of one to four key value pair arguments from varnames varname or cavarnames fcnnames fcnname or cafcnnames outfile fname detail detaillevel The first three pairs are as described in printInitData and printConstr The detaillevel can be compact default which causes this function to
26. ify A with its n x n incidence matrix aij where a 1 if i j A 0 otherwise A natural sparsity pattern for the DAE is the set where the entries of X are finite S i 7 0g gt co the sparsity pattern of 2 3 6 However a more informative BTF comes from the sparsity pattern of the Jacobian J It depends as does J on the valid offset vectors c d used So Solc d i j dj ci 0g the sparsity pattern of J 3 8 Since d cj oi holds on each HVT by 3 1 and implies o gt oo we have Sa e d C S for any c d Experience suggests that in applications a BTF based on So is usually significantly finer than one based on S We refer to the former as a fine BTF and we refer to the latter as coarse BTF for more detail see 3 4 Advantage of block triangularization 16 3 4 Advantage of block triangularization We now discuss by means of an example why it is advantageous to find a BTF First we note that at each k stage in the solution scheme it may be necessary to prescribe initial values see 82 3 5 The equations 3 4 3 5 can do without initial values when k is large enough that they become square linear systems The latest this happens is at stage k 0 However it can happen for parts of the DAE before one reaches stage k 0 as illustrated by 3 9 comprising two simple pendula with a coupling term 0 A z zA 0 B y yA G 0 0 r y L 0 D u up 0 E v
27. ix indices of equations and variables c d global offsets HVT position boxed in BW print yellow in color print structural non zeros in the system Jacobian shaded in BW print colored green in color print e With disptype blocks or disptype fineblocks the figure displays B permuted signature matrix B S pe pv pe pv permuted indices of equations and variables c d cl dl global and local offsets union of HVTs boxed in BW print yellow in color print structural non zeros in the system Jacobian shaded in BW print colored green in color print e With argument submat matValue or blocksubmat blkValue the figure only displays DAE function name submatrix or blocks information submatrix or blocks of the signature matrix indices of equations and variables union of HVTs boxed in BW print yellow in color print structural non zeros in the system Jacobian shaded in BW print colored green in color print If either the width or height of the matrix displayed is larger than 40 the figure does not show the indices and the offsets Contents of figure if SIP If any variable or equation of index in 1 n is not accessed the DAE is SIP and the equations or variables that make the DAE SIP will be highlighted dark shading in BW print colored red in color print and figure will display e When no disptype is specified
28. ks Block 4 6 Using x x x y y y lam Solve linear 3x3 system O fil I27 13 for x y lam Block 1 1 Solve nothing give initial value for V STAGE k 2 3 blocks Block 4 6 Using x x x x y y y y lam lam Solve linear 3x3 system 0 f1 f2 f3 for X DUE yov lam Block 3 3 Using lam lam lam v Solve nonlinear equation give trial value O 6 for u Block 1 1 Solve nothing give initial value for v STAGE k 1 3 blocks Block 4 6 Using X x x KILL xo tu y y ys y 5 lam Solve linear 3x3 system O fl f2 37 5 for x 5 y 5 lam Block 3 3 Using lam lam lam lam u v v Solve linear equation O f6 for vw Block 1 1 Solve nothing give initial value for vos STAGE k 0 4 blocks Block 4 6 Using X x x gt IRA x 5 y vo y lam lam lam lam Solve linear 3x3 system y lam lam 22 R yr yin y 5 O 1 2 3 6 for x 6 y 6 lam Block 3 3 Using lam lam lam lam lam u wu Solve linear equation 0 2 6 for v Block 2 2 Using u wu u Solve linear equation v v vo 2 3 Extract structural analysis data 9 O f4 for mu Block 1 1 Using v v v mu Solve nonlinear equation give trial val
29. nalysis Structural analysis is performed by the daeSA function function sadata daeSA daefcn n options Here daefcn is a function for evaluating the DAE n is the size of the problem and options is an optional list of one or more arguments that are passed by daeSA to daefcn The return object sadata encapsulates data obtained from analyzing the DAE represented in daefcn It can be accessed by the functions described in and For our problem we set values for the size and the additional parameters and then call daeSA n 6 G 9 8 L 1 0 c 0 1 sadata daeSA modified2pendula n G L c We shall use this sadata in the examples that follow 2 3 Extract structural analysis data 5 Remark 3 The result of the SA is not affected by the values of G L and c However we provide this mechanism for passing them to daefcn as such constants matter when computing a consistent point evaluating the system Jacobian and integrating the DAE DAESA does not perform these three tasks but we plan to implement them in the future 2 3 Extract structural analysis data One can visualize the structure of a DAE and obtain SA data as follows 2 3 1 Visualization To display the DAE structure we call showStruct sadata to produce the coarse block triangularized DAE structure we call showStruct sadata disptype blocks and to produce the fine block triangularized DAE structure we call showStruct sadata disptype
30. ns DAEs It allows convenient translation of a DAE into MATLAB and provides a set of currently 18 easy to use functions for determining and visualizing key structural properties of the DAE The package is applicable to DAE systems of the general form fi t x and derivatives of them 0 i 1 n 1 1 where t is the independent variable and the x t are n state variables The formulation 1 1 includes high order systems and systems that are nonlinear in leading derivatives Furthermore 1 1 includes systems of ordinary differential equations ODEs and pure algebraic systems In the next paragraphs a concept is set in slanting type on first occurrence with a forward reference to further details about it DAESA performs analysis that is similar to the one the C solver DAETS does 5 6 However DAETS is not suitable for rapid investigation of DAEs as it requires C knowledge and compiling the user code DAESA is a light weight easy to use tool for SA that provides convenient facilities for analyzing a DAE to find its structural index 83 1 the number of degrees of freedom henceforth referred to as the DOF 83 1 the constraints and required initial values 83 5 and a solution scheme 83 1 and also to reduce the DAE to a coarse or a fine block triangular form B TF 83 3 which can be exploited for efficient solution in a block wise fashion DAESA applies MATLAB s operator overloading to process a DAE given by a user su
31. oblem MOD2PEND an index 7 non quasilinear DAE see also 8 0 f z zA 0 fa y y G 0 fs x0 y L 2 1 0 Jj u uu 0 fs v vu G 0 fg u v L cA 2 In the original two pendulum problem fs v up G and fg u v L cAY 3 2 2 Perform structural analysis 4 The state variables are x y A u v and u G is gravity L gt 0 is the length of the first pendulum and c gt 0 is a given constant Figure shows an encoding of 2 1 function f modified2pendula t z G L c x z 1 y z 2 la z 3 u z 4 v z 5 mu z 6 first pendulum oO 00 N o ar De w N f 1 Dif x 2 x la ha xA 0 f 2 Dif y 2 y la G hy yA G 0 3 x 2ty 2 L 2 4 z y L 90 modified second pendulum 4 Dif u 2 u mu u up 0 1 5 Dif v 3 72 v mu G 0 vu G 0 11 f 6 u72 v72 L c la 72 Dif la 2 12 u 0 L cA 1 0 13 end gt s Figure 2 1 DAESA function for evaluating 2 1 The Dif var k operator returns the kth derivative of variable var Constants G L and c are passed as additional parameters In lines 2 3 the input variables z i are renamed for better readability Equations 2 1 are translated in lines 5 11 Remark 2 Specifying a DAE in DAESA is similar to how ODEs are specified in MATLAB e g for ode45 except that instead of returning y in y f t y we return a vector of the f s in 1 1 2 2 Perform structural a
32. onzeros in system Jacobian J Size 15 structural index 11 DOF 10 Boxed HVT Shaded structural nonzeros in system Jacobian J Boxed positions that contribute to det J Indices of Variables Indices of Variables 123 4 5 67 8 9101112131415 ci 13 14151011127 8 9 4 5 6 1 23 cj 8 o 8 15 o o B 2 2 2 Bg 8 14 Bo 0 0 3 a B 10 B p 0 0 4 A O 6 12 la o B 2 4 E a a 6 5 11 2 0 0 2 2 S 10 B p 0 2 Ss 6 ooo 8 8 S 7 A 0 4 9 oo 9 2 6 amp 8 2 0 0 4 5 8 B B 4 2 7 a o 0 4 Q n g 9 8 pg 6 a WW E EN 2 E 10 A 0 2 5 El dl 0 6 11 E 8 2 4 H fo 0 6 12 boo 4 3 o bl 2 10 13 B o 2 Bg o e 14 Bg o 1 El bo s 15 a a 2 dj220220220220220 dj10 10 8 886664442220 dj 2 2 0 4 4 2 6 6 4 8 8 6 1010 8 a original structure b fine block triangularized structure Figure 5 10 Structure of 5 1 and its fine block triangularization Chapter 6 Installation This package is available at http www cas mcmaster ca nedialk daesa as the zipped file DAESA 1 0 zip The implementation is distributed in the form of MATLAB pcode files When DAESA 1 0 zip is unzipped it creates directory DAESA 1 0 in the current directory with subdirectories as follows subdirectory contains src implementation pcode files of the implementation src interface m files that execute the corresponding pcode files examples the code for the examples in this User Guide examples DAEs the functions for evaluating the DAEs used in this User Guide
33. pplied function for evaluating the in 1 1 First using operator overloading DAESA extracts the DAE s signature matrix 83 1 from whose sparsity pattern the coarse block triangular form can be found It then finds out if the problem is structurally well posed and if so solves a linear assignment problem to calculate the global offsets 83 1 of the DAE which give the structural index and DOF 93 1 Using both the signature matrix and the offsets DAESA constructs the fine block triangular form finds local offsets 83 5 for the blocks and for each block determines if it is quasilinear in the leading derivatives 3 5 Based on the local offsets and linearity information DAESA deduces a minimal set of variables and derivatives of them that need to be initialized and a minimal set of inherent constraints The package provides functions for displaying the original sparsity structure of the DAE as well as for displaying its coarse and fine block triangularizations and functions for reporting the constraints initialization summary and a solution scheme for the DAE Remark 1 The structural index computed by DAESA is an upper bound on the differentiation index and in our experience it usually equals it That is DAESA usually finds the smallest possible set of differentiations Although successful on many problems of interest the underlying SA theory and DAESA respectively may fail to determine the correct structural and therefore
34. quation in block 1 1 does not require solving but at this stage we give an initial value for v which is used in later stages Stage k 2 We have three blocks 4 6 3 3 and 1 1 For block 4 6 we use the computed SE uu y y A X to solve 0 I n je for 9 y A For block 3 3 we use computed A X A v and give a trial value for u to solve 0 fs Finally for block 1 1 we give an initial value for v which is used in following stages Stage k 1 Again we have three blocks 4 6 3 3 and 1 1 For block 4 6 we use the computed a x q 0 y y y y y A X A to solve 0 fl fl E for 20 yO A For block 3 3 we use the computed A X A A u v v to solve 0 f for u Lastly for block 1 1 we give an initial value for v to be used in the positive number stages Stage k gt 0 We solve linear systems as summarized in Table 2 3 Extract structural analysis data 10 block using solve for Jab eos poet 465 gu ca mh gH Gh Gm gt yh 000 AA cS ACTA 3 3 WW NE Ln T u 3 Vi 9 9 uu utt s l Bon oo i A V v Ves v2 k k 3 1 1 dd 10 T U Table 2 1 MOD2PEND solution process for stages k gt 0 Compact solution scheme A more compact representation of the solution scheme is produced by printSolScheme sadata varnames vars outfile solschemecompact txt detail compact which results in the file solschemecompact txt Compact
35. ries that contribute to det J which is the union of all HVTs 4 Global offsets are denoted by c dj and local offsets are denoted by Gj d 2 3 Extract structural analysis data 6 2 3 2 Index and DOF The index and the DOF are obtained by index getIndex sadata DOF getDOF sadata which return 7 and 5 respectively 2 3 3 Initialization summary The function printInitData reports which variables and derivatives of them need to be initialized for a numerical solution of the analyzed DAE The call vars x y lam u v mu printInitData sadata varnames vars outfile initdata txt produces file initdata txt with modified2pendula problem Initialization summary X x y u v v v y That is x z y y u v v v and v need to be initialized By default printInitData outputs as variable names x followed by a variable number from 1 to n Here these names are replaced using the key pair varnames vars If the key pair outfile and initdata txt are omitted the function will print the output to the command window 2 3 4 Constraints The DAE constraints are reported by the function printConstr For example the call printConstr sadata outfile constr txt produces file constr txt with modified2pendula problem Constraints fl as OS E EIA dl A ES f27 412 FOP 28s de Et 975 ESITO I f3 5 5 6 6 f6
36. sec secant asech hyperbolic arcsecant csc cosecant acsch hyperbolic arccosecant cot cotangent acoth hyperbolic arccotangent asin arcsine acos arccosine atan arctangent asec arcsecant acsc arccosecant acot arccotangent Table A 1 Standard functions supported by DAESA 40 Index Functions daeSA 4 21 B7 getBTF getConstr getDA EfhandlegetDAEfhandle getDOF 6 getHVT 23 31 get Index 6 getInitData getMissingEqnsVars getOffsets getQLdata 26 82 getSigma 23 getSize isSWP printConstr 6 printInitData 6 29 printSolSchene 7 30 36 Compact Full showStruct Coarse blocks Diagnostic blocks Fine blocks Theory Block triangularization Constraints Degrees of freedom Highest value transversal Incidence matrix Initial values Linear assignment problem Offsets Canonical Global Local Quasilinearity 17 Signature matrix 1 12 16 Signature tableau Solution stages Sparsity pattern Structural index Structurally ill posed Structurally well posed System Jacobian Val 2 E Bibliography 1 F MAZZIA AND F IAVERNARO Test set for initial value problem solvers Tech Rep 40 Department of Mathematics University of Bari Italy 2003 http pitagora dm uniba it 2 N S NEDIALKOV J PRYCE AND G TAN DAESA a Matlab tool for structural analysis of DAEs Software 2013 Accepted for publication in ACM TOMS
37. sis by hand For further details see 3 4 9 Our method gives essentially the same result as does that of Pantelides 7 but is easier to use The steps to perform SA are as follows 1 Form the n x n signature matrix X 0 of the DAE where order of the derivative to which the jth variable v occurs in g 4 the ith equation f or oo if x does not occur in fi 2 Find a highest value transversal HVT of X A transversal T is a set of n positions in the matrix with one entry in each row and each column We denote Val T the sum of the c in those positions and we aim at finding T such that Val T is as large as possible The value of the signature matrix X written Val 2 is defined as the value of any HVT The value of any transversal and of X is either an integer or oo The DAE is structurally well posed SWP if Val 2 is finite that is if there exists at least one transversal all of whose gij are finite Otherwise the DAE is structurally ill posed SIP which implies there is some error in the problem formulation 3 Find the global offsets integer vectors c c1 Cn and d d d with all c gt 0 that satisfy dj ci gt o4 forall i j 3 1 12 3 1 Structural analysis 13 with equality holding on the HVT They are not unique but there exist canonical smallest offsets in the sense of c lt c if c e for all i We like to show the results by a signature tableau which
38. solution scheme for modified2pendula problem Initialization summary x x y y u V v v y k 6 f3 x y k 5 37 x y k 4 f1 2 3 x y lam k 3 fi 2 3 x y lam Q v ko SD f1 122 f3 gt gt gt gt Eb EA ya lam f6 u i o ov k 1 1 gt f2 3 5 x 5 y 5 lam f6 v eye k 0 1 gt gt 27797 3 6 x 6 y 6 lam t6 7 aa f4 mu Ses r3 On the right of is the list of variables for which we solve The brackets mark the system being solved and means no equations are solved but value s for variable s must be given The 2 3 Extract structural analysis data 11 notation denotes that the block of equations inside these brackets is non quasilinear in its highest order derivatives If no detail level is specified then the compact level is given by default As in omitting the key pair outfile and filename txt from the printSolScheme function will result in outputs being printed to the command window Chapter 3 Theory Overview In this chapter we first introduce the basic theory and terminology of structural analysis SA and then explain the solution method the coarse and fine block triangular forms and the quasi linearity analysis QLA 3 1 Structural analysis Here we give enough theory that a user can on small problems do the analy
39. stem Thus at any stage k lt 0 we are using derivatives found either at stage k or at a previous stage due to the incremental nature of the formula Use the notation J to denote the system Jacobian at stage k Then Jx is a submatrix of J in the sense of being obtained by selecting certain rows and columns of J not necessarily contiguous However see 9 if we order the equations and variables in descending order of offsets J becomes the leading my x nj submatrix of J for some my and ny such that my lt ny The Jj are nested in that both my and ny are non decreasing with k and equal n when k gt 0 3 3 The coarse and fine block triangularizations 15 For the simple pendulum the system Jacobian takes the needed form when the equations are ordered h f g as shown below J_2 and J which are equal are shown boxed ha 54 0 J a 0 f A 0 Jy gr The solution scheme for the simple pendulum is given in Table k solve equations for variables 2 h 25 1 M x y 0 9 h Dv y 1 T g h g y N Table 3 1 Solution stages for the simple pendulum The structural analysis method uses 3 4 to find Taylor coefficients to solve the DAE via Taylor series 3 3 3 The coarse and fine block triangularizations A block triangular form BTF of 1 1 is obtained from a sparsity pattern which is some subset A of 1 n the n x n matrix positions i j for i and j from 1 to n When appropriate we ident
40. the elements of iv is the total number of IVs needed For instance iv 2 3 means IVs for x2 25 1 are necessary iv 4 0 implies x4 does not participate in initialization If SIP iv function constr getConstr sadata constr Output If SWP constr is a non negative integer row n vector such that 0 je for 0 r lt constr i are constraints As a result the sum of the elements of constr gives the total number of constraints If SIP constr For instance constr 3 3 implies 0 fs f4 f7 are constraints constr 1 0 implies fi is not in the set of constraints If SIP constr 4 3 Obtaining structural data 25 function pe pv getBTF sadata function pe pv cb getBTF sadata function pe pv cb fb getBTF sadata Let S getSigma sadata Assume SWP pe pv Output pe pv are two row n vectors such that S pe pv is block upper triangularized signa ture matrix where the entries below the diagonal blocks are all oo The th equation in the permuted system is pe i and the jth variable in the system is pv Cj fb Output cb fb are row integer vectors that specify the boundaries of the diagonal blocks of the coarse and fine block triangularization respectively If B S pe pv the i j th block of the coarse block triangularization is B cb i cb it 1 1 cb j cb j 1 1 Similarly the i j th block of the fine block trian gularization is B fb i fb i 1 1 fb j fb j 1 1 Ass
41. triangularized structure of the DAE gt fineblocks fine block triangularized structure of the DAE gt submat matValue matValue is a row vector of size 2 or 4 specifying what submatrix of the signature matrix S to be displayed size displays 2 diagonal submatrix comprising rows matValue 1 to matValue 2 and columns matValue 1 to matValue 2 submatrix comprising rows matValue 1 to matValue 2 and columns matValue 3 to matValue 4 blocksubmat blkValue blkValue is a row vector of size 2 or 4 specifying what blocks of the signature matrix S to be displayed size displays 2 diagonal blocks comprising block rows blkValue 1 to blkValue 2 and block columns blkValue 1 to blkValue 2 submatrix comprising block rows blkValue 1 to blkValue 2 and block columns blkValue 3 to blkValue 4 Note to use the key pair blocksubmat blkValue the disptype must be set to give a block triangularization Assume SIP A similar set of outputs to the SWP case are available except the disptype fineblocks is no longer available and the disptype blocks now produces a diagnostic block triangularization as detailed below and in 2 4 4 Visualization 28 Contents of figure if SWP Without using arguments submat or blocksubmat the figure displays the DAE function name size structural index and DOF e In the default setting the figure also displays S signature matr
42. uce the number of needed initial values 3 5 Quasi linearity analysis This section illustrates how to do linearity analysis for the blocks we obtain by doing the decompo sition After permuting the DAE to upper triangularized form of b blocks where possible we calculate the local offsets and d for each block Then we analyze the linearity of the highest order derivatives HODs of the variables in each block Denote ax k 1 b such that ax 1 if the HODs within block 4 are linear and ax 0 otherwise With the permuted variables we associate them with an n vector such that j is the block of variable z j 1 n Similarly we associate the permuted equations with a vector m such that m i is the block of equation f i 1 n For variables xj the required IVs are E dj ay 3 Xie ica 3 5 Quasi linearity analysis 18 where d denotes local d offsets here For equation f in block k the set of constraints are Bee oda ir my where c denotes global c offsets Therefore the required IVs for the system are X X4 X5 T Xn and the constraints are F FE E5 Fn 0 where some of the X s can be the empty vector and some of the F s can be missing 3 10 In the quasilinear case an z whose d is zero it must have c lt 0 for all i and thus is a purely algebraic variable does not appear in the vector X Similarly an f with c 0 does not appear in F
43. ue O f5 for v In this output Block p q denotes the sub block in the fine block triangularized DAE comprising rows p to q and columns p to q For example Block 4 6 implies equations fs fa fi and variables x y A cf Figure 2 2 c At the head of the output is an initialization summary giving all the variables and derivatives of them that need initial values When such a value is needed in the following solution process it is marked by give trial value when it is needed but no equations are solved it is marked by give initial value Trial values are guesses by the user which a solver may change to satisfy constraints initial values are not changed in the solution process We discuss stages 6 to 1 and gt 0 Stage k 6 We have a scalar equation fs 0 which is the first equation in block 4 6 cf Figure and Figure 2 2 c Since it is nonlinear in x and y trial values for x and y are necessary and we need to determine values for them such that f3 0 is satisfied Stage k 5 We use the previously computed x y give trial values for x y and solve f 0 which is linear in z and y Stage k 4 We use computed x z y y to solve 0 fi fo fY for x y A no trial values are needed as this system is linear in z y A Stage k 3 We have two blocks 4 6 and 1 1 For block 4 6 we use computed x 2 x Y V V A to solve 0 fi fo f3 for x y X Then the e
44. ume SIP pe pv Output pe pv are two row n vectors such that S pe pv is block upper triangularized signa ture matrix where the entries below the diagonal blocks are all oo The th equation in the permuted system is pe i and the jth variable in the system is pv j Output cb an empty array a 2 x m positive integer matrix specifying the boundaries of the diagnostic see below block triangularization where m is the number of fine blocks minus one If B S pe pv the i j th block of the diagnostic block triangularization is B fb 1 i fb 1 i 1 1 fb 2 j fb 2 j 1 1 Note The number of output arguments should be either 2 3 or 4 Practically the diagnostic block triangularized form will look like e block B11 B12 is under determined By Bi Bi Bog e block B s is well determined and e block B34 B44 is over determined 4 3 Obtaining structural data 26 function cql getQLdata sadata function cql fql getQLdata sadata Assume SWP cql fal Output row vectors containing information about the quasi linearity of each diagonal block in the coarse and fine block triangularizations of a DAE structure respectively Let pe pv cb fb getBTF sadata For the coarse block triangularization if cq1 i true then the ith diagonal block comprising equations of indices pe cb i cb i 1 1 in variables of indices pv cb i cb i 1 1 is quasilinear in the leading derivatives of these variables otherwis
45. vu G 0 F u L ce 3 9 where G L c are constants Its signature matrix and system Jacobian are X y u V Hh ci xr y Alu v HL A 2 0 1 A 1 a B 2 0 1 B 1 y C 0 0 3 C e 24 IST 2 Par Tsp 1 u E 2 Ojo E 1 Ful 0 0 2 FUE 2u 2v dj 8 3 1 2 2 0 where 2c L cz A blank in X means c0 and in J means zero One can see just from X or from J without studying the equations themselves that the system splits into parts i e subsystems part 1 has equations A B C for variables x y A part 2 has equations D E F for variables u v u Part 1 influences part 2 by the F x entry in the lower left block of X but is uninfluenced by it since the top right block is blank thus giving a block lower triangular form BTF This coupling leaves the offsets of part 2 unchanged but increases by one those of part 1 from c 0 0 2 d 2 2 0 toc 1 1 3 dj 3 3 1 This seems to change the initial values part 1 requires paradoxical since part 1 is the uninfluenced one Namely the combined DAE is quasilinear and by considering the stages defined by 3 5 we find that IVs are needed for I Mh E Ue rri gt Y Y Y A UU au 3 5 Quasi linearity analysis 17 So part 1 now seems to need IVs for x y and A which as a stand alone system it did not see Table Of course this is false When one considers the offsets as defining a scheme for computing successive derivatives equiv
46. x2 in equations 0 fi fa For more information on the diagnostic block triangularization see the description of getBTF on p 5 3 Chemical Akzo Nobel We show DAESA s SA on the chemical Akzo Nobel problem 1 an index 1 DAE The function for evaluating the corresponding equations is in Figure function f akzonobel t y 8 7 k2 0 58 4 4 KIA 3 3 ki 1 K 3 Ks 1 ri r2 r3 r4 rb Fin f 1 2 3 f 4 f 5 6 end Figure 15 83 k3 0 09 k4 0 42 C02 0 9 H 737 k1xy 1 74 sqrt y 2 k2 y 3 y 4 k2 K y 1 y 5 k3 y 1 y 4 2 k4 y 6 2 sqrt y 2 k1A C02 H y 2 Dif y 1 1 Dif y 2 1 Dif y 3 1 Dif y 4 1 Dif y 5 1 Ks y 1 y 4 2 0 r1 r2 r3 r4 0 5xr1 r4 0 5 r5 Fin ri r2 13 r2 r3 2 0xr4 r2 r3 r5 y 6 5 5 DAESA function for evaluating the chemical Akzo Nobel DAE 5 3 Chemical Akzo Nobel 36 The script in Figure 5 6 produces the solution scheme in Figure 5 7 and the plot in Figure sadata daeSA akzonobel 6 showStruct sadata disptype fineblocks printSolScheme sadata outfile chemakzo txt varnames y Solution scheme for akzonobel problem Figure 5 6 DAESA script for analyzing AkzoNobel Initialization summary y3 y4 y5 yl y2 E ce k 0 f6 y6 f5 y5 f4 y4 f3 y3 2
Download Pdf Manuals
Related Search