Delay Differential-Algebraic Equations
vorgelegt von Phi Ha
geb. in NamDinh, Vietnam
von der Fakultät II - Mathematik und Naturwissenschaften der Technischen Universität Berlin
zur Erlangung des akademischen Grades Doktor der Naturwissenschaften
Dr. rer. nat.
-genehmigte Dissertation
Promotionsausschuss:
Vorsitzender: Prof. Dr. Yuri Suris
Gutachter: Prof. Dr. Volker Mehrmann Gutachter: Prof. Dr. Peter Kunkel Gutachterin: Prof. Dr. Caren Tischendorf
Tag der wissenschaftlichen Aussprache: 09. März 2015
Doing a PhD is a long hard work, and during the time working on this thesis, I have received the support from many persons.
First and foremost, I would like to express my deep gratitude to my advisor Prof. Dr. Volker Mehrmann for introducing me to the topic of Delay-DAEs, for his patience, support, advice and criticism. My thanks also go to all my colleagues in the Numerical group of the Technical University Berlin, where I have enjoyed the friendly atmosphere and the productive work environment. In addition, personal thanks are due to my former office mate Tobias Brüll, for being my weird and special friend, and also for his emotional suggestion to the behavior approach; to Helia Niroomandrad for three years sharing office with me; to Andreas Steinbrecher for many good and free lessons about DAEs; to Benjamin Unger for carefully reading this work and many valuable correc-tions and suggescorrec-tions. Assistance provided by Ma Vinh Tho was greatly appreciated.
I also wish to acknowledge the help and support in living and working provided by my Vietnamese friends, in particular Dr. Dam Quoc Phoi. With your friendship, you really lighted up my last one and half year in Berlin.
I’m grateful that this work is supported by the DFG project SFB 910 Control of self-organizing nonlinear systems: Theoretical methods and concepts of application during my four year PhD at the Technical University Berlin. As a student, I also want to thank the Berlin Mathematical School for their support, in particular the transition period during the end of Phase I and the beginning of Phase II.
Finally, this thesis is dedicated to my parents and my brother. Without my family’s encouragement, love and support, this work would never have been existent.
In dieser Arbeit legen wir die Grundlage für die Analyse der Lösbarkeit von Anfangs-wertproblemen (IVPs) für verzögerte differentiell-algebraische Gleichungen (DDAEs). Eine DDAE ist eine allgemeine Kombination von zwei mathematischen Objekten: dif-ferentiell-algebraischen Gleichungen (DAEs) und Differentialgleichungen mit Zeitver-zögerungen (DDEs). Diese Kombination führt zu vielen Schwierigkeiten in der Analyse von DDAEs, die weder für DAEs noch für DDEs auftreten. Zwei von diesen Schwie-rigkeiten, die große Auswirkungen auf die theoretische und numerischen Lösung von DDAEs haben, sind die Nicht-Kausalität und das Vorlaufen (advancedness) eines Sys-tems, die bisher für DDAEs nicht diskutiert wurden. Diese Schwierigkeiten zu überwin-den und neue Erkenntnisse über die Analyse der Lösbarkeit von DDAEs zu entwickeln sind das Thema der Arbeit, deren Inhalt aus drei wesentlichen Schwerpunkten besteht. Das Studium von linearen DDAE Systemen mit variablen Koeffizienten ist der erste Schwerpunkt der Arbeit. Im Gegensatz zu kausalen DDAEs, die sich ähnlich wie DAEs verhalten, muss (sogar für lineare) nichtkausale DAEs die Strukturtheorie von mehr als zwei Matrizen verwendet werden, und daher kann die Analyse der Lösbarkeit der ent-sprechenden IVPs durch die klassische Theorie von Matrizenpaaren nicht angewen-det werden. Ein wichtiger Beitrag des ersten Teils der Arbeit ist, den Zusammenhang zwischen Existenz und Eindeutigkeit von Lösungen von IVPs für DDAEs und der Regu-larität von entweder einem Matrizenpaar oder einem Matrizentripel zu untersuchen. Hierbei ist zu unterscheiden, ob das Zeitintervall beschränkt oder unbeschränkt ist. Das erzielte Ergebnis erlaubt uns, die Lösbarkeit einer linearen DDAE durch Spektra-leigenschaften der Koeffizientenmatrizen zu studieren. Wir wenden Matrix-Polynome an, um die Lösbarkeit einer großen Klasse von zeitverzögerten Systemen einschließ-lich unter- und überbestimmter Systeme zu analysieren. Zwei Klassen von DDAEs, die häufig in Anwendungen Verwendung finden sind DDAEs vom retardierten und vom neutralen Typ. Für diese beiden wichtigen Klassen untersuchen wir einen weiteren, insbesondere algebraischen Ansatz zur Analyse ihrer Lösbarkeit. Dieser Ansatz zeigt die Ähnlichkeit dieser zwei Typen von DDAEs zu DAEs, und ist geeignet für weitere Studien über Eigenschaften hinsichtlich Regelung und Steuerung von DDAEs.
Um die theoretische Grundlage für die numerische Lösung von IVPs für DDAEs zu erzielen, betrachten wir als zweiten Schwerpunkt allgemeine Systeme mit linea-ren zeitvarianten Koeffizienten. Es wird gezeigt, dass klassische schrittweise Metho-de (’Method of steps’) für nichtkausale Systeme nicht anwendbar ist, und es werMetho-den zwei neue Ansätze für die Regularisierung von allgemeinen linearen DDAEs vorgestellt. Der erste Ansatz modifiziert die ’method of steps’ so, dass die Lösung auf fortschrei-tenden Teilintervallen des Integrationsbereichs berechnet werden kann. Es wird ein
Shift-Index eingeführt, der ein Maß für die Nichtkausalität des System darstellt. Damit wird die Regularisierungstheorie von DAEs auf DDAEs verallgemeinert. Abhängig vom Typ von der DDAE, führt diese Regularisierungstechnik zu verschiedenen strangeness-freien Formulierungen, die sehr kompliziert sein können, falls eine DDAE von vorlau-fendem Typ ist. Die zweite Methode formuliert das AWP für die DDAE als ein Rand-wertproblem (RWP) für eine hochdimensionale DAE um, sodass die DAE Theorie für die Regularisierung angewendet werden kann. Weiter wird gezeigt, dass die so konstru-ierte DAE in diesem differentiell-algebraischen RWP einen beliebig großen Strangen-essindex haben kann.
Der letzte Teil der Arbeit ist die Entwicklung von zwei neuen Integrationsalgorith-men für allgemeine lineare DDAEs, welche auf den Ansätzen des zweiten Teils dieser Arbeit basieren. Die Hauptidee dieser Algorithmen ist es, bekannte numerische Me-thoden (für DDEs und für DAEs) auf die strangeness-freie Formulierung anzuwenden, die über oben genannte Regularisationstechniken erhalten wird. Der erste Algorith-mus, der auf einer Verallgemeinerung der ’Method of Steps’ basiert, behandelt DDAEs von retardierten und neutralen, aber nicht von vorlaufendem Typ. Der zweite Algo-rithmus, welcher differentiell-algebraische RWPe löst, ist für alle drei Typen geeignet. Abschließend zeigen wir die Effizienz und Robustheit der vorgestellten Algorithmen im Vergleich mit dem häufig benutzten numerischen Softwarepaket RADAR5 [60].
With this thesis, we aim to lay the foundation for the solvability analysis of initial value problems (IVPs) for delay differential-algebraic equations (DDAEs). In our context, a DDAE is a general combination of two important mathematical objects: differential-algebraic equations (DAEs) and delay differential equations (DDEs). This combination has led to many difficulties in the analysis of DDAEs, which occur neither for DAEs nor for DDEs. Two of these difficulties, which have strong influence on the theoretical and numerical solutions of DDAEs, are the noncausality and the advancedness of a system, previously not discussed for DDAEs. To overcome these difficulties in order to give new insights to the solvability analysis of DDAEs is the purpose of this work, whose main content is focused on three topics presented in the following.
The investigation of systems with time invariant coefficients is the first topic of in-terest. In contrast to causal DDAEs, which possess many similar properties like DAEs, already in the linear case, noncausal DDAEs take into account the structure of more than two matrices and consequently, the solvability analysis of the corresponding IVPs can not be completely analyzed by the classical theory of matrix pairs. An important contribution of the first part of this thesis is to point out the link between the existence and uniqueness of solutions of IVPs for DDAEs and the regularity of either a matrix pair or a matrix triple, depending on whether the time interval is bounded or not. This re-sult allows to study the solvability of a linear DDAE by investigating spectral properties of its matrix coefficients. In more details, we apply a matrix polynomial approach to study the solvability analysis for a much broader class of time delay systems including both under- and overdetermined systems. Another approach in the theory of DAEs, namely an algebraic method, is examined to study general DDAEs of retarded and neu-tral types, which usually occur in applications. This approach has shown the similarity between these two types of DDAEs and DAEs, and it is suitable for further investigation on control properties of DDAEs.
Aiming at the theoretical background for the solution procedure for DDAEs, the second topic involves the consideration of general systems with linear time varying coefficients. Observing the failure of the classical method of steps for noncausal sys-tems, we propose two new approaches for the regularization of general linear DDAEs. The first approach aims to modify the method of steps so that one can compute the solution on consecutive sub-intervals of the integration time interval. By introduc-ing the shift index concept to estimate the noncausality of a system, we generalize the regularization procedure of DAEs for DDAEs. Depending on the types of DDAEs, this regularization technique leads to different strangeness-free formulations, which can be very complicated if a DDAE is of advanced type. On the other hand, the second
ap-proach reformulates an IVP for a DDAE as a BVP for a high dimensional DAE so that one can make use of the DAE theory for the regularization. It is further shown that the con-structed DAE in this differential-algebraic BVP can have an arbitrarily high strangeness index.
The last part of this thesis is the development of two new integration algorithms for general linear DDAEs by using the two approaches introduced in the second part. The main core of these algorithms is to apply well-known numerical methods (for DDEs or for DAEs) to the strangeness-free formulation obtained by different regular-ization techniques discussed above. The first algorithm, based on a generalregular-ization of the method of steps, successfully handles DDAEs of retarded and neutral types but not systems of advanced type. On the other hand, the second algorithm, based on solving differential-algebraic BVPs, is suitable for dealing with all three types of DDAEs. Con-cluding, the efficiency and robustness of both algorithms are demonstrated in com-parison with the commonly used numerical solver RADAR5 [60].
Hiermit versichere ich an Eides statt, dass ich die vorliegende Arbeit selbstständig und ohne Benutzung anderer als der angegebenen Hilfsmittel angefertigt habe. Die aus fremden Quellen direkt oder indirekt übernommenen Gedanken sind als solche kennt-lich gemacht.
Teile der vorliegenden Arbeit sind in Kooperation mit Prof. Dr. Volker Mehrmann und Dr. Andreas Steinbrecher entstanden und sind als solche kenntlich gemacht. Kapitel2,
3und5zitieren Ergebnisse der Artikel
[63] Analysis and Reformulation of Linear Delay Differential-Algebraic Equations. P. Ha and V. Mehrmann. Electr. J. Lin. Alg., Vol. 23 (2013), pp. 703-730, 2013, [65] Analysis of Linear Variable Coefficient Delay Differential-Algebraic Equations.
P. Ha, V. Mehrmann and A. Steinbrecher. J. Dynam. Differential Equations, pp. 1-26, 2014.
Die, in den Artikeln [63,65], diskutieren Ideen und deren Umsetzung sowie die Lite-raturrecherche waren überwiegend bis außschließlich mein Anteil. Die Ausarbeitung dieser Artikel wurde mehrheitlich von mir ausgeführt. Prof. Dr. Volker Mehrmann und Dr. Andreas Steinbrecher standen mir dabei mit ihrem Wissen und ihrer Erfahrung zur Seite.
Hiermit erkläre ich, die vorliegende Arbeit weder im Inland noch im Ausland in glei-cher oder ähnliglei-cher Form als Dissertation, Diplom- oder ähnliche Prüfungsarbeit an-gemeldet oder eingereicht zu haben und weiterhin keine Promotionsabsicht an einer anderen Hochschule oder Fakultät beantragt zu haben.
Der wissenschaftliche Inhalt der Arbeit wurde weder in seiner Vollständigkeit noch in der vorliegenden Form veröffentlicht.
Kapitel5ist jeweils die überarbeitete und erweiterte Fassungen des Preprints [64] Analysis and numerical solution of linear delay differential-algebraic equations.
P. Ha and V. Mehrmann, in preparation. Kapitel2,3und5zitieren Ergebnisse der Artikel
[63] Analysis and Reformulation of Linear Delay Differential-Algebraic Equations. P. Ha and V. Mehrmann. Electr. J. Lin. Alg., Vol. 23 (2013), pp. 703-730, 2013, [65] Analysis of Linear Variable Coefficient Delay Differential-Algebraic Equations.
P. Ha and V. Mehrmann and A. Steinbrecher. J. Dynam. Differential Equations, pp. 1-26, 2014.
Nomenclature xix
Abbreviation xxi
1 Introduction 1
1.1 Applications . . . 4
1.1.1 Human balance control and multibody control systems . . . 4
1.1.2 Time-delayed electric circuits. . . 7
1.1.3 DDAEs in fluid dynamics . . . 8
1.1.4 DDAEs in chemical engineering . . . 9
2 Fundamentals of DAEs and DDEs 11 2.1 Time Invariant Differential-Algebraic Equations . . . 11
2.2 Time Varying Differential-Algebraic Equations . . . 13
2.3 High-order Differential-Algebraic Equations . . . 21
2.4 Classification of DDEs . . . 24
2.5 Solution concepts and the method of steps for DDEs . . . 25
3 Fundamentals of DDAEs 29 3.1 Basic concepts of DDAEs . . . 29
3.2 Prior work on DDAEs . . . 32
3.3 Characteristics of general linear DDAEs . . . 35
3.4 Transforming multiple delays into single delay. . . 38
4 Solvability Analysis of General Linear Time Invariant DDAEs 41 4.1 Basic properties . . . 41
Contents xiv
4.2 Systems on bounded time intervals . . . 43
4.3 Systems on unbounded time intervals. . . 47
4.3.1 Structures of matrix triples via constant equivalence transforma-tions . . . 47
4.3.2 A matrix polynomial approach . . . 54
4.3.3 Reformulation of non-advanced systems . . . 59
5 Solvability Analysis of General Linear Time Varying DDAEs 63 5.1 Generalization of the method of steps . . . 64
5.2 Solvability analysis via system classification . . . 70
5.3 Regularization of DDAEs by operational arrays . . . 73
5.4 Boundary Value Problem method . . . 79
6 Numerical solutions of IVPs for DDAEs 87 6.1 Application of the generalized method of steps for DDAEs . . . 87
6.2 Application of the BVP method to DDAEs . . . 89
6.3 Examples . . . 93
7 Solvability Analysis of General Nonlinear DDAEs 97 Conclusion and Outlook 101 Conclusion. . . 101
1.1 Inverted one-bar pendulum. . . 5
1.2 Time-delayed electric circuit, [97], p. 536.. . . 7
1.3 Williams-Otto process, [109], p. 304. . . 9
2.1 Classification of DDEs with discrete delays, compare [120]. . . 25
2.2 Discontinuity propagation of the solutions to (2.25) (left) and (2.26) (right). 26 3.1 Numerical solution and absolute error of the IVP (1.2) for (3.19) with con-stant stepsize h = 0.1.. . . 35
5.1 Absolute error vs. length of the time intervalIfor (5.45). . . 84
5.2 Numerical solution and absolute error of (5.46). . . 85
6.1 Numerical solution and absolute error of the IVP (6.7) with constant step-size h = 0.01. . . . 90
6.2 Numerical solution and absolute error of the IVP (6.14) with constant stepsize h = 0.01. . . . 93
6.3 Relative error of the DDAE (6.15) (left) and of the DDAE (6.16) (right) with constant stepsize h = 0.01. . . . 95
6.4 The solution and the absolute error of the DDAE (6.17) with constant stepsize h = 0.01. . . . 96
5.1 Strangeness-free formulation for different classes of DDAEs. . . 73
5.2 Consistency of initial conditions for different classes of DDAEs. . . 74
5.3 Theoretical comparison of the generalized method of steps and the BVP method. . . 85
m number of equations
n number of unknowns
I= [0, tf) time interval
Iτ= [−τ, tf) time interval including history
N0(resp.,N) set of non-negative (resp., positive) integer numbers Z (resp., Q) set of integer numbers (resp., rational numbers) R (resp., R+) set of real numbers (resp., positive real numbers)
C set of complex numbers
Rm,n(resp.,Cm,n) set of m by n matrices whose entries are inR (resp., C) Cm,n[ξ] set of m by n matrix polynomials of variableξ, p.55
Ck(I,Cn) vector space of all k-times continuously differentiable functions fromIintoCn
˙
x, ¨x, x(i ) total derivatives of x(t ) with respect to t σ(E, A) spectrum of a matrix pair (E , A), p.9
σ(E, A,B) spectrum of a matrix triple (E , A, B ), p.33
ν (or ν(E, A)) Kronecker-index of a regular matrix pair (E , A), p.12
µ (or µ(E, A)) strangeness-index of a DAE or of a matrix function pair (E , A), p.17
κ(i) shift-index of a linear DDAE with respect to i , p.58
I , In identity matrix (∈ Rn,n)
ED Drazin inverse of a matrix E
UT transpose of a matrix (or a matrix function) U
UH conjugate transpose of a matrix (or a matrix function) U ker(M ) (right) null space of a matrix M
corange(M ) (left) null space of a matrix M , i.e., ker(MT)
diag(A1, . . . , Ak) the block diagonal matrix where the blocks on the main diagonal are A1, . . . , Ak
⊗ Kronecker product of two matrices
vec the vec operator
IVP Initial Value Problem BVP Boundary Value Problem
ODE Ordinary Differential Equation DDE Delay-Differential Equation
DAE Differential-Algebraic Equation (without delay) DDAE Delay Differential-Algebraic Equation
gcd Greatest common divisor lcm Least common multiple
Introduction
In many real-life applications, the mathematical models are described by Delay Differe-ntial-Algebraic Equations (DDAEs) of the following form
F(t , x(t ), ˙x(t ), x(t − τ)) = 0, (1.1) on a given time intervalI= [0, tf) ⊂ R+, where ˙x denotes the first (time) derivative of the vector valued function x. Hereτ > 0 is a constant delay. The function x maps from
Iτ:= [−τ, tf) toCn. The function F takes values inCm. Furthermore, it is allowed that the time intervalIcan be unbounded, i. e., tf = ∞.
By linearizing general DDAEs of the form (1.1) around a trajectory, we obtain a linear time varying DDAE
E (t ) ˙x(t ) = A(t)x(t) + B(t)x(t − τ) + f (t), (1.2a) where the coefficients are matrix-valued functions E , A, B :I→ Cm,n, f :I→ Cm. In particular, the linearization at a stationary solution results in a linear time invariant DDAE where the functions E , A, B become constant matrices. The reason for calling (1.1) and (1.2a) delay differential-algebraic equations is that the partial derivative Fx˙ and the matrix E (t ) can be rank deficient, and consequently, systems (1.1) and (1.2a) can contain not only differential constraints but also some algebraic constraints. In general, to obtain an IVP, the following initial function is usually added to (1.1) (resp. (1.2a))
x|[−τ,0]= φ : [−τ, 0] → Cn. (1.2b) Talking about linear DDAEs, one could think about them as the general combination of two well-known mathematical objects: first, (non-delayed) Differential-Algebraic Equations (DAEs) of the form
E (t ) ˙x(t ) = A(t)x(t) + f (t), (1.3) and second, Delay Differential Equations (DDEs) of the form
˙
x(t ) = A(t)x(t) + B(t)x(t − τ) + f (t), (1.4) Differential-Algebraic Equations (DAEs), since the pioneering work of Gear [51], have been of interest in a tremendous amount of research, especially in the last two decades,
2
due to their vital role in automatic modeling and the rapid advancement of modern computers, which makes it possible to solve very large problems numerically. There now exists a collection of surveys and monographs in this field and we in particular refer to [6, 18, 22, 49, 72, 75, 82, 110]. The automatic modeling approach, making use of DAEs, is a convenient tool to model physical and engineering problems that involve constraints: for example, electronic networks (taking into account the Kirch-hoff’s laws), or mechanical systems (accounting to the positions constraints of a mov-ing mass on a surface), or chemical processes (involvmov-ing mass or energy conservation laws). Based on DAEs, implemented software packages such as Modelica (Dymola) [42], Matlab Simulink [91,92], or Spice [2,133], have enabled the possibility to model and simulate very complex physical phenomena in a very simple manner.
Compared with DAEs, time delay systems, in particular Delay Differential Equa-tions, have a much longer history, date back to 1874 with an equation describing elas-ticity effects proposed by Boltzman ([20]). Until now, this is still a very active research field with a great number of references. One of the main reasons for this long time interest is that mathematical models are expected to describe real-life processes as good as possible and many of them include time lag in their dynamics. This type of delay is often referred as intrinsic delays. Various well-known examples of real-life sys-tems involving intrinsic delays can be found in classical references [15,40,56,74,97] as well as in very recent monographs [46,95,126,134], ranging in numerous applica-tions in medicine and biology, population dynamics, chemistry, economics, viscoelas-ticity, physics, mechanics, engineering sciences, etc. In addition, another rich source of time delay systems comes from the control context. The cause is that any control system involving feedbacks will almost certainly involve time delay, due to the exis-tence of a finite time that is required to sense information and then to react to it, see e.g. [37,40,95,126,134]. Furthermore, delays can also be used as control parameters to achieve desired behaviors, for example stability or controllability, of a dynamical sys-tem, see [71,100,107,141]. Some other sources of delay, which are also indispensable in applications, can be found in [58,111,125].
The above historical introduction shows that delays must be taken into account to simulate or to control certain processes, and in situations where the processes are described by DAEs then certainly one has to analyze DDAEs. DDAEs, therefore, are of high practical relevance and they present a young field with increasing importance, see e.g. [5,10,11,19,25,29–31,41,60,61,63,65,83,84,94,105,122,131,142,143]. How-ever, the mixture of these two types of equations (DAEs and DDEs) leads to substantial mathematical difficulties and many open questions for DDAEs, even for fundamental problems such as the solvability analysis of linear time invariant systems.
In order to understand general nonlinear DDAEs of the form (1.1), the first and very important step is to understand linear DDAEs of the form (1.2a). This work mostly aims to lay the foundation for the solvability analysis to the corresponding Initial Value Problem (1.2), where we are interested in both theoretical and numerical solutions. The short outline of this work is as follows.
In Chapter2we deliver some important results on DAEs and on DDEs that will be useful later. While the first part of this chapter, consisting of the first three sections,
fo-cuses on DAEs, the second part of this chapter concentrates on DDEs. In particular, the first part presents the concept of strangeness index, the strangeness-free form and the reformulation algorithm to transform a given DAE into its strangeness-free form. The second part of this chapter discusses the solution concept and the classification for time delay systems. Finally, the method of steps, a standard tool for investigating the-oretical and numerical solutions of time delay systems is presented. On one side, the presented results here will guide our thinking and show us what to expect for DDAEs. On the other side, we will show that the study of linear DDAEs present many difficulties which occur in neither the theory of DAEs nor that of DDEs.
Basic concepts and properties of the linear DDAE (1.2a) are discussed in Chapter
3. Thereafter illustrating that implicit and more complicated DDE structures rather than the DDE (1.4) may occur, we propose in Section3.1the piecewise differentiable solution concept as well as the system classification for DDAEs. We briefly review prior work about DDAEs in Section3.2. In particular, the application of the method of steps to causal systems, as it has been presented in the literature [5,10,25,29,60, 83,122,
131, 142,143], is re-examined. We however demonstrate that the method of steps is not suitable for noncausal systems and consequently, should be modified for general DDAEs. Important characteristics of DDAEs presented in Section3.3show that DDAEs are neither DAEs nor DDEs, and actually they merit separate investigation in their own rights, as has been commented in [10,65]. Then, we propose some challenging prob-lems, which motivate the research in subsequent chapters. Finally, Section3.4 con-structs the transformation to convert a multiple delay DDAE into a single delay one, leaving the trajectory invariant. Using this transformation, all the results derived for single delay systems, in particular the solvability analysis, can at once be extended to multiple delay systems.
The foundation of the present work is mainly built in the next three chapters where we consecutively investigate Linear Time Invariant systems (Chapter4), Linear Time Varying systems (Chapter5) and the (numerical) solution procedures for Linear Time Varying systems (Chapter6).
Most pertinent in the analysis of linear time invariant DAEs is the well-known Kronecker-Weierstraß matrix pencil theory. One important result states that the ex-istence and uniqueness of a solution to the DAE (1.3) is closely related to the regularity of the matrix pair (E , A). This motivates our study in Chapter 4, where we address the solvability analysis of linear time invariant DDAEs. An important contribution of Chapter4 is to point out the relation between the existence and uniqueness of a solution to the DDAE (1.2a) and the regularity of either the matrix pair (E , A) or the matrix triple (E , A, B ), depending on whether the time interval Iis bounded or not. The proposed matrix polynomial approach considered in Section4.3.2even shows the applicability to a much broader class of equations including both underdetermined and overdetermined systems. Thereafter we present as third approach the algebraic method [75,127] to study DDAEs of retarded and neutral types with special emphasis on the similarity between these two types of DDAEs and DAEs.
1.1. Applications 4
Aiming at the theoretical background for the numerical solution of the IVP (1.2), Chapter5deals with general systems of linear time varying coefficients. Observing the failure of the method of steps, we propose two new approaches to handle general lin-ear DDAEs, which aim at the same goal is to solve the IVP (1.2) and also to obtain the necessary and sufficient conditions for a consistent initial function. The first approach, presented in Sections5.1-5.3, aims to modify the method of steps so that one can still compute the solution of the IVP (1.2) on consecutive intervals. More important, in Section5.3we present Algorithm5.2to construct the strangeness-free formulation of retarded and neutral DDAEs by using operational arrays. The sufficient conditions for the successful implementation of this algorithm is given in Hypothesis 5.22. On the other hand, the second approach aims to remove the delay by reformulating the IVP (1.2) as a BVP for a high dimensional DAE, and hence the solution of the IVP (1.2) is obtained by solving this differential-algebraic BVP. This approach is analyzed in Sec-tion 5.4. We further show that the constructed DAE in this BVP can have arbitrarily high strangeness index, which is proportional to the length of the time intervalI. To round up this chapter, the comparison of these two approaches and some illustrative examples are presented to confirm the theoretical results.
In Chapter 6, we consider the numerical solution of the IVP (1.2) using the two approaches introduced in Chapter5. To follow the first approach, under Hypothesis
5.22, we determine, via Algorithm5.1, the strangeness-free formulation of the DDAE (1.2a) pointwise. The numerical solution of the IVP (1.2) is obtained by implementing Radau collocation methods with Lagrange interpolation (for x(t − τ)) to the resulting strangeness-free DDAEs of Algorithm5.1. On the other hand, to follow the second ap-proach, we compute the numerical solution of the differential-algebraic BVP by Gauß-Lobatto collocation methods, as presented in the literature [79,80,130]. The existence, uniqueness, and the convergence results for the numerical solution of the IVP (1.2) are given in Theorems6.2,6.4.
For the sake of completeness, in Chapter 7we review some important results in prior work about the solvability analysis of IVPs for nonlinear DDAEs. Therein, we also discuss the limitation of prior studies and the motivation for further research. Finally, we give some conclusions and possible research problems found during this work.
1.1 Applications
Due to the broad range of applications of differential-algebraic systems and the natural physical meaning of the time delay, it is clear that DDAEs play a crucial role in the physical modeling and theoretical understanding of numerous applications in science and engineering. In this section we present several practical examples, which aim to give the readers a glimpse in how DDAEs occur in practice.
1.1.1 Human balance control and multibody control systems
Falling accidents for older peoples often occurs while walking and in many cases the immediate cause is not simply “slips and trips” but unknown [88,96]. In order to
min-1.1. Applications 6
taneously implemented based on the assumption that the balance control could be entirely depended on the biochemical properties of the joints, connective issues, etc. of the human body [136,137]. However, as demonstrated in subsequent experiments, not only these forces but also neural feedback control takes part in the control mech-anisms for balance [87,98]. This leads to an important consequence that time-delays must be included in the control force, due to the fact that there is a significant time lag, since the variables are measured until the force is applied (for human body, the neu-ral latencies are approximately from one to five hundred milliseconds, [96,128,138]). Consequently, the force applied to the cart becomes→−F (t −τ) and system (1.5) becomes a delay differential-algebraic system. It is worth to note that in [96,128,138] and the references therein, the time delay feedback control force−→F is chosen as
− →
F = k1X (t − τ) + k2θ(t − τ) + k3X (t − τ) + k4˙ θ(t − τ),˙
where the constants ki, i = 1,...,4, are chosen so that the upright position of the pen-dulum is stabilized. Furthermore, it is assumed that all the measurements of X , ˙X ,θ, ˙θ occur at the same time.
The inverted pendulum without time delay is one special case of multibody sys-tems, which are frequently modeled by differential-algebraic equations. A multibody system is a mechanical system described in terms of bodies and connections, where the mass is assumed to be concentrated entirely in rigid or elastic bodies. The bodies are coupled by massless connections, which are the source for constraints of the body motion. For example, let us consider the motion of np bodies described by npposition coordinates p(t ) and np velocity coordinates v(t ). Connections like joints cause holo-nomic constraints 0 = g (p, t), together with other differential constraints will lead to the system of mixed differential and algebraic equations for the dynamics of the multi-body system. By utilizing a Lagrange multiplierλ, one obtains Lagrange equations of the first kind
˙ p(t ) = v(t), M (p,t) ˙v(t) = →−F − ∂ ∂tg (p, t ) ¶T λ, (1.6) g (p, t ) = 0,
whereM (p,t) denotes the mass matrix and−→F stands for the external force. In the case that the external force is instantaneous, i.e.,→−F =−→F (p(t ), v(t ), t ), for standard works on the theoretical and numerical solutions of (1.6), we refer to [44,112,118,119,124,127] and the references therein.
However, because of the time delay in neural feedback control as in the case of hu-man balance control, or because of unavoidable time delays in both controllers and actuators, the force−→F in the presence of delays takes the form
− →
F =−→F (p(t − τ1), v(t − τ2), t )
whereτ1andτ2are the time delays in the paths of displacement and velocity feedback, respectively. In this case, the system (1.6) becomes the following delay
differential-algebraic system ˙ p(t ) = v(t), M (p,t) ˙v(t) = −→F (p(t − τ1), v(t − τ2), t ) − ∂t∂ g (p, t ) ¶T λ, (1.7) g (p, t ) = 0.
1.1.2 Time-delayed electric circuits
In conventional circuit analysis, the circuit equations are usually set up by standard ap-proaches for the systematic formulation like the modified nodal analysis [69]. It turns out that the resulting systems are differential-algebraic equations, see [7,47,48,62,82,
110,132,133]. On the other hand, early work of time delays in electric circuits was con-sidered by Minorski in [97] (1962). The circuit is given in Figure1.2. Here L, C are
con-Figure 1.2: Time-delayed electric circuit, [97], p. 536.
stants of the circuit that involves in series two resistors R1and R. The voltage across R1 is applied to a linear amplifier A and then applied through a special phase-shifting net-work P that creates a constant time delayτ between the input and the output. Under these conditions, the voltage across R in series with the output of P is VP= GR1i (t −τ), where G being the gain of the amplifier to R measured through the network P . In this manner, the following system is derived in [97]
uL(t ) = Ld i (t ) d t , ˙ uC = 1 Ci (t ), uL(t ) + uC(t ) + (R + R1)i (t ) +GR1i (t − τ) = 0,
1.1. Applications 8
which can be rewritten in the form of the linear DDAE 2 4 0 0 L 0 1 0 0 0 0 3 5 2 4 ˙ uL(t ) ˙ uC(t ) ˙ i (t ) 3 5= 2 4 1 0 0 0 0 C1 1 1 R + R1 3 5 2 4 uL(t ) uC(t ) i (t ) 3 5+ 2 4 0 0 0 0 0 0 0 0 GR1 3 5 2 4 uL(t − τ) uC(t − τ) i (t − τ) 3 5. (1.8)
Since the late 1960s, circuits which include delayed elements turn out to be popular in the engineering community due to the fact that these circuits have a very important role due to the increase of performance of very large scale integration (VLSI) systems, [12,21,114–116]. The two typical types of circuits where delays usually occur are cir-cuits with transmission lines (TL) [21], and partial element equivalent circuits (PEEC’s) [114,115]. These circuits result in neutral DDEs of the following form, see [12],
˙
y(t ) = Ly(t) + M y(t − τ) + N ˙y(t − τ), for all t ≥ t0, (1.9) together with an initial functionφ(t). The delay τ > 0 is constant and t0is the starting time. For the solution procedure, directly applying numerical methods such as Runge-Kutta methods or multi-step methods requires an approximation of the derivative ˙y(t − τ), which is a difficult task. The simple trick to overcome this difficulty is to introduce a new functionΨ(t) := y(t) − N y(t − τ) and to rewrite (1.9) as
˙
Ψ(t) = LΨ(t) + (M + LN)y(t − τ), y(t ) = Ψ(t) + N y(t − τ),
which in fact is the DDAE • I 0 0 0 ‚ •Ψ(t)˙ ˙ y(t ) ‚ =•L 0 I −I ‚ •Ψ(t) y(t ) ‚ +•0 M + LN 0 N ‚ •Ψ(t − τ) y(t − τ) ‚ .
Certainly, one may argue that this DDAE, which is reformulated from a neutral DDE, is quite artificial. However, the neutral DDE (1.9) is indeed a reduced system obtained by neglecting certain variables and certain algebraic equations, using the Kirchhoff’s voltage law or the Kirchhoff’s current law. In general, the unreduced system, which involves algebraic constraints obtained by these laws, is a delay differential-algebraic system.
1.1.3 DDAEs in fluid dynamics
The dynamical behavior of a system in fluid mechanics and turbulence modeling is often described by the incompressible Navier-Stokes equation of the form
∂u
∂t − ν∆u + ∇p + (u · ∇)u = f in (0, ∞) × Ω, ∇ · u = 0 in (0, ∞) × Ω,
whereν > 0 is the viscosity, u = u(t,ξ) is the velocity field which is a function of the time t and the position ξ, p is the pressure, f is the external force. Recently, there has been an increasing interest in the situation where the trajectories of some fluid particles have a delayτ to follow the fluid [85, 104]. Furthermore, from the control
perspective, it is favorable to control the system by another external force g = g (t,u(t − τ,ξ)) which involves some hereditary characteristics [34,50]. This leads to the following time-delayed version of the incompressible Navier-Stokes equation
∂u
∂t − ν∆u + ∇p + (u(t − τ, ξ) · ∇)u = f + g (t , u(t − τ, ξ)) in (0, ∞) × Ω,
∇ · u = 0 in (0, ∞) × Ω, (1.10a)
together with the following initial and boundary conditions u = 0 on (0, ∞) × ∂Ω, u(0, x) = u0(x) in Ω,
u(t , x) = φ(t, x) in (−τ,0) × Ω.
(1.10b)
To obtain the numerical solution to the initial-boundary value problem (1.10), the fre-quently used method is to discretize the space variable by finite difference or finite element methods [57] and consequently one obtains a delay differential-algebraic sys-tem.
1.1.4 DDAEs in chemical engineering
Chemical engineering is also a field where DAEs and DDEs occur in many applications, and within this section we present one frequently discussed example in the literature namely the Williams-Otto process, [113,120,135]. The schematic of the flow-sheet for this process is shown in Figure1.3. Two kinds of raw materials A and B upon entering
Figure 1.3: Williams-Otto process, [109], p. 304.
the chemical reactor take part in three chemical reactions which produce the desired product P , along with some by-products. The feed rates of the raw materials (in pounds per hour) are denoted by FAand FB, respectively. In order to settle an undesirable by-product out of the reactant mixture, a heat exchanger is used to cool the reactants to a certain temperature. This settling takes place in the decanter and after that the ma-terial, which contains the desired product, impurities, and a certain percentage of the raw material with some by-products of the chemical reaction, will enter a distillation
1.1. Applications 10
column. There the valuable product is separated from the impurities, whereas the raw material with the by-products is recycled to the chemical reactor to be reprocessed. The recycle loop, which ensures that useful products will not be discarded, introduces a significant time delay into the problem. In practical situations, it is not at all unusual for material to take ten minutes to travel from the chemical reactor through the cooler, the decanter, the distillation column, and then recycling to the reactor. This is because of the distance separating the various stages of the overall process and lengths of pip-ing in the stages themselves.
Even though the originally proposed differential equation governing this chemical pro-cess is nonlinear, see [135], for the determination of a proper correction in the feed rates at the desired operating point, a corresponding linearized model is useful, see [113]. For a recycle time of ten minutes, the linearized equation describing the reac-tion in the chemical reactor is the following DDE with the delayτ = 10
˙ x(t ) = A0x(t ) + B0x(t − τ) + f (t), (1.11) with A0 = 2 6 6 6 4 −4.93 −1.01 0 0 −3.20 −5.30 −12.8 0 6.40 0.347 −32.5 −1.04 0 0.833 11.0 −3.96 3 7 7 7 5 , B0= 2 6 6 6 4 1.92 0 0 0 0 1.92 0 0 0 0 1.87 0 0 0 0 0.724 3 7 7 7 5 , f (t ) = 1 6VR 2 6 6 6 4 1 0 0 1 0 0 0 0 3 7 7 7 5 • δFA δFB ‚ , x = 2 6 6 6 4 xA xB xC xP 3 7 7 7 5 ,
where VR is the volume of the chemical reactor (VR ≈ 2.628m3),δFA andδFB are the deviations in the feed rates of the raw materials A and B , respectively, from their nomi-nal values. The dimensionless components xA, xB, xC and xPrepresent the deviations from their nominal values in the weight compositions of the raw materials A and B , of the intermediate product C , and of the desired product P , respectively. Note that six variables appear in [135] but the other two variables xG1 and xG2correspond to the
waste product G, which does not react with the other chemicals, and so only the four variables given above occur in the DDE (1.11).
Furthermore, taking into account the mass conservation law, one obtains the following algebraic constraint
xA(t ) + xB(t ) + xG1(t ) + xG2(t ) + xC(t ) + xP(t )
= g (t , xG1(t − τ), xG2(t − τ), xC(t − τ), xP(t − τ)) + δFA+ δFB,
(1.12)
where the scalar function g represents the weight of the recycled raw material with the by-products. The combined system (1.11)-(1.12) therefore gives rise to a DDAE in the variable
x(t ) =£xA(t ) xB(t ) xG1(t ) xG2(t ) xC(t ) xP(t )
⁄T .
Fundamentals of DAEs and DDEs
Before studying DDAEs, it is worthwhile to spend some time to recall the essentials of the theory of differential algebraic equations (DAEs) and of delay differential equations (DDEs) as well.
The first part of this chapter, consisting of the first three sections, focuses on DAEs. In the first two sections, we recall the solvability analysis of first order systems, which follows by an extension to arbitrarily high order systems in the third section. In more detail, the existence and uniqueness of a solution to a DAE is studied via the spectral analysis (Section2.1), or via the structure of the matrix function coefficients (Sections
2.2and2.3).
The second part of this chapter collects results on DDEs. Section2.4recalls important facts concerning the solution concept and the classification of different classes of time-delay equations. These have important consequences in the discontinuity propagation and the numerical integration of DDEs. We then present in Section2.5the method of steps, a standard tool for investigating the theoretical as well as the numerical solutions to DDEs.
2.1 Time Invariant Differential-Algebraic Equations
In this section we consider time invariant differential-algebraic equations of the form
E ˙x(t ) = Ax(t) + f (t). (2.1)
This equation (2.1) is obviously a special case of (1.2a) where the matrix function B is identically zero, and the matrix functions E , A become constant matrices inCm,n.
Definition 2.1. The matrix pair (E , A) in (Cm,n)2 is called regular if m = n and the so-called characteristic polynomial p(λ) defined by p(λ) = det(λE − A) is the non-zero polynomial. Equivalently, the matrix pair (E , A) is regular if and only if the (finite) spec-trum σ(E, A) := {λ ∈ C | det(λE − A) = 0} is not the entire C.
Definition 2.2. A square matrix N ∈ Cn,nis called nilpotent of nilpotency indexν = ν(N) if Nν= 0 and Nν−16= 0.
We recall one basic fact that the nilpotency index does not exceed the size of a nilpotent matrix. Once the matrix pair (E , A) is regular, the solvability of (2.1) is char-acterized by using Kronecker-Weierstraß canonical form, see e.g., [26,38].
2.1. Time Invariant Differential-Algebraic Equations 12 Theorem 2.3. (Kronecker-Weierstraß canonical form) Let E , A ∈ Cn,n and (E , A) be regular. Then there exist nonsingular matrices W , T ∈ Cn,nsuch that
(W E T, W AT ) = •Id 0 0 N ‚ ,• J 0 0 Ia ‚¶ , (2.2)
where N is a nilpotent matrix. Moreover, it is allowed that one or the other blocks is not present.
Definition 2.4. Consider a pair (E , A) of square matrices that is regular and has a
canon-ical form as in (2.2). The quantityν defined by
ν = (
0 if the block N is absent in (2.2), ν(N) if N is present in (2.2),
is called the Kronecker index of the matrix pair (E , A) and is denoted byν(E, A).
Definition 2.5. Let E ∈ Cn,nand letν = ν(E, I) be the Kronecker index of the matrix pair (E , I ). A matrix X ∈ Cn,n satisfying
E X = X E , X E X = X , X Eν+1 = Eν, is called the Drazin inverse of E .
Clearly, the Kronecker indexν(E, I) exists for every square matrix E, and it is zero if and only if E is nonsingular. In this case, the Drazin inverse of E is nothing else than E−1. Otherwise, the existence and uniqueness of the Drazin inverse is given in the following theorem.
Theorem 2.6. Let E ∈ Cn,nand P be a nonsingular matrix such that
E = P•C 0
0 N
‚ P−1, where C is regular and N is nilpotent. Then,
ED = P•C −1 0
0 0
‚ P−1. Proof. For the proof, see [33].
Letλ0be an arbitrary number such that the inverse of matrix λ0E − A exists, we define
e
E := (λ0E − A)−1E , A := (λ0E − A)e −1A, ˜f := (λ0E − A)−1f .
Making use of the Drazin inverse, one obtains the explicit representation for the solu-tion of (2.1) in the following theorem.
Theorem 2.7. Let the matrix pair (E , A) as in (2.1) be regular. Furthermore, let f ∈ Cν(I,Cn) withν is the Kronecker index of the matrix pair (E, A). Then every solution x of the DAE (2.1) has the form
x(t ) = eEeDAte e EDE x(0) +e Z t 0 eEeDA(t −s)e e EDE ˜ef (s)d s − (I − eE D e E ) ν−1 X j =0 (EeAe D )jAe D ˜ f( j )(t ), for all t ∈I.
Proof. For the proof see Theorem 3.1.3, [26].
In the case that the matrix pair (E , A) is not regular, it is well-known that the cor-responding IVP for the DDAE (2.1) either has more than one solution or there are ar-bitrarily smooth inhomogeneities for which there is no solution at all, as presented in the following theorem.
Theorem 2.8. Let (E , A) be as in (2.1) and suppose that (E , A) is a singular matrix pair. 1. If rank(λE − A) < n for all λ ∈ C, then the homogeneous initial value problem
E ˙x(t ) = Ax(t), x(0) = 0, has a nontrivial solution.
2. If rank(λE − A) = n for some λ ∈ C and hence m > n, then there exist arbitrarily smooth inhomogeneities f for which the corresponding differential-algebraic equation is not solvable.
Proof. For the proof see Theorem 2.14, [75].
2.2 Time Varying Differential-Algebraic Equations
As we have seen in Section2.1, the solvability analysis of the time invariant DAE (2.1) is characterized by the structure of the matrix pair (E , A). The aim of this section is to study linear time varying coefficients DAEs of the form
E (t ) ˙x(t ) = A(t)x(t) + f (t), for all t ∈I, (2.3a)
x(0) = x0, (2.3b)
by analyzing the structure of the matrix function pair (E , A). To do that, we shall slightly modify the algebraic approach [75] by transforming system (2.3a) only from the left. The reason for choosing this approach is that it is suitable for not only uniquely solv-able DAEs, but also for over- and under-determined DAEs, which is an important ad-vantage for our study of DDAEs later.
The solution concept for DAEs of the form (2.3a) is stated in the next definition.
Definition 2.9. A function x :I→ Cnis called:
i) a (classical) solution of (2.3a) if x ∈ C1(I,Cn) and x satisfies (2.3a) pointwise, ii) a (classical) solution of the IVP (2.3) if x is a solution of (2.3a) and satisfies (2.3b).
2.2. Time Varying Differential-Algebraic Equations 14
An initial vector x0is called consistent to system (2.3a) if the IVP (2.3) has a solution. System (2.3a) is called solvable if it has at least one solution. It is called regular if in addition, for any consistent initial vector x0, the corresponding IVP (2.3) has a unique solution.
We will make frequent use of the following results, compare Theorems 3.9, 3.25 in [75].
Theorem 2.10. Let E ∈ C`(I,Cm,n),` ∈ N0∪{∞}, with constant rank E(t ) = r for all t ∈I. Then there exist pointwise unitary functions U ∈ C`(I,Cm,m) and V ∈ C`(I,Cn,n), such that UHEV = •Σ 0 0 0 ‚ , or UHE =•E10 ‚ ,
with pointwise nonsingularΣ ∈ C`(I,Cr,r), and E1has full row rank r .
Theorem 2.11. LetI⊂ R be a closed interval and M ∈ C (I,Cm,n). Then there exist open intervalsIj ⊂I, j ∈ N, with
[ j ∈N
Ij=I, Ii∩Ij= ; for i 6= j, and integers rj∈ N0, j ∈ N such that
rank M (t ) = rj for all t ∈Ij.
Lemma 2.12. For the pair (P,Q) with P ∈ C`(I,Cp,n), Q ∈ Ck(I,Cq,n),`, k ∈ N0∪{∞}, as-sume that there exist two integers rQ6r[P ;Q]such that rankQ(t ) = rQ and rank
•P (t ) Q(t ) ‚
= r[P ;Q]for all t ∈I. Then, there exists
• S 0 Z1 Z2
‚
∈ Cmin{`,k}(I,Cp,p+q) that satisfies the fol-lowing conditions. i) • S Z1 ‚ ∈ C (I,Cp,p) is pointwise unitary, ii) Z1P + Z2Q = 0,
iii) the function SP has pointwise full row rank, and the pair (SP,Q) satisfies
rank •SP Q
‚¶
= rank(SP ) + rank(Q).
Proof. Since Q has constant rank onI, one can apply Theorem2.10to factorize Q, and then partition P conformably to get
2 4 Ip 0 0 U11H 0 U12H 3 5 •P Q ‚ £V11 V12⁄ = 2 4 P1 P2 Σ 0 0 0 3 5 p rQ q − rQ , (2.4)
where U1=£U11 U12⁄ ∈ Ck(I,Cq,q), V1=£V11 V12⁄ ∈ Ck(I,Cn,n) are pointwise unitary functions, andΣ ∈ Ck(I,CrQ,rQ) is pointwise nonsingular. The sizes of the block rows in
(2.4) are p, rQ, q − rQ. Moreover, note that in (2.4), P2also has constant rank due to
rank(P2) = rank •Ip 0 0 U1H ‚ •P Q ‚ £V11 V12 ⁄ ¶ − rank(Σ) = r[P ;Q]− rQ.
Then, by Theorem 2.10, there exists a pointwise unitary function U2H=• S
Z1 ‚
∈ Cmin{`,k}(I,Cp,p) such that U2HP2= • S Z1 ‚ P2= •P12 0 ‚ , (2.5)
where P12∈ Cmin{`,k}(I,Cr[P ;Q]−rQ,n−rQ) has pointwise full row rank.
Combining (2.4) and (2.5), one obtains 2 6 6 6 4 S 0 Z1 0 0 U11H 0 U12H 3 7 7 7 5 •P Q ‚ £V11 V12⁄ = 2 6 6 6 4 P11 P12 P21 0 Σ 0 0 0 3 7 7 7 5 r[P ;Q]− rQ p − r[P ;Q]+ rQ rQ q − rQ ,
where P12has pointwise full row rank andΣ is pointwise nonsingular onI.
Consequently, SP =£P11 P12⁄ V1−1 has pointwise full row rank. Moreover, one sees that rank •SP Q ‚¶ = rank •P11Σ P12 0 ‚¶ = rank¡£0 P12⁄¢ + rank¡£Σ 0⁄¢ = rank(SP ) + rank(Q).
SinceΣ ∈ Ck(I,CrQ,rQ) is pointwise nonsingular, it implies thatΣ−1∈ Ck(I,CrQ,rQ).
Fi-nally, setting Z2:= −P21Σ−1U11H∈ Cmin{`,k}(I,Cp−r[P ;Q]+rQ,q), we obtain
Z1P + Z2Q =¡[P21 0] − P21Σ−1[Σ 0]¢V1−1= 0, which completes the proof.
Following the algebraic approach [75], we rewrite equation (2.3a) in the form E (t )d
dt− A(t ) ¶
x(t ) = f (t), (2.6)
for any t ∈I. For notational convenience, we will omit the time variable t in all matrix functions.
Making use of Theorem2.11and restricting ourselves if necessary to subintervals, we may assume that the following assumption holds.
Assumption 2.13. For the pair of matrix functions (E , A) of the DAE (2.3a), there exist integers r, a such that
rank(E ) = r, rank¡£E A⁄¢ = r + a for all t ∈I.
2.2. Time Varying Differential-Algebraic Equations 16 Lemma 2.14. Consider the DAE (2.3a) and suppose that Assumption2.13holds. Then, there exists a pointwise unitary function P1∈ C (I,Cm,m) such that by scaling system (2.6) with P1from the left one obtains a new system in the following form
2 4 M11dtd − M12 −M22 0 3 5x = 2 4 f1 f2 f3 3 5 r a v , (2.7)
where the functions M11∈ C (I,Cr,n), M22∈ C (I,Ca,n) have pointwise full row rank. Here the sizes of the block row equations are r, a and v = m − r − a.
Proof. First we determine a pointwise unitary function PE :I→ Cm,mvia Theorem2.10 or a smooth QR-decomposition, see [39], that compresses the matrix function E . This yields PE E d dt − A ¶ =•M11 d dt− M12 − ˜M22 ‚ r m − r ,
such that M11 has full row rank. Continuing, by compressing the block ˜M22 with a pointwise unitary function PA:I→ Cm−r,m−r, this yields
•Ir 0 0 PA ‚ PE E d dt− A ¶ =•Ir 0 0 PA ‚ •M11d dt− M12 − ˜M22 ‚ = 2 4 M11dtd − M12 −M22 0 3 5,
where M11and M22have pointwise full row rank. Setting P1:=
•Ir 0 0 PA
‚
PE, we arrive at (2.7).
The formula (2.7) in Lemma2.14clearly shows that the number of scalar (nontriv-ial) differential equations in system (2.3a) is r , while the number of scalar (nontrivial) algebraic constraints is a.
In the following we suppose that the function M22is continuously differentiable. Again, to be able to apply Lemma2.12, the following assumption is necessary.
Assumption 2.15. For the DAE (2.7), there existsm ∈ N such that the functions M11, M22b satisfy
rank •M11 M22
‚¶
=m, for all t ∈b I.
Under Assumption 2.15, applying Lemma 2.12to the pair (M11, M22) implies the existence of matrix functions S, Z1, Z2 of appropriate sizes that have the following properties
i) the function• S Z1
‚
∈ C (I,Cr,r) is pointwise unitary,
ii) the function SM11has pointwise full row rank and the following identities hold:
and
rank •SM11 M22
‚¶
= rank(SM11) + rank(M22). Define the operator
P2:= 2 6 6 6 4 S 0 0 Z1 Z2dtd 0 0 Ia 0 0 0 Iv 3 7 7 7 5 d s a v , (2.9)
where r = d + s, we see that P2has a left-inverse given by the formula
P2−1= 2 6 6 6 4 • S Z1 ‚−1 −• S Z1 ‚−1• 0 Z2dtd ‚ 0 0 Ia 0 0 0 Iv 3 7 7 7 5 d + s a v .
Applying the operator P2to system (2.7), we obtain 2 6 6 6 4 S 0 0 Z1 Z2dtd 0 0 Ia 0 0 0 Iv 3 7 7 7 5 2 4 M11dtd − M12 −M22 0 3 5x = 2 6 6 6 4 S 0 0 Z1 Z2dtd 0 0 Ia 0 0 0 Iv 3 7 7 7 5 2 4 f1 f2 f3 3 5 d s a v , or equivalently, 2 6 6 6 4 SM11dtd − SM12 Z1M11dtd − Z1M12− Z2dtd M22 −M22 0 3 7 7 7 5 x = 2 6 6 6 4 S f1 Z1f1+ Z2f˙2 f2 f3 3 7 7 7 5 d s a v . (2.10)
According to dtd M22= ˙M22+ M22dtd and (2.8), it follows that the second block equation of (2.10) becomes ¡−Z1M12− Z2M˙22¢ x = Z1f1+ Z2f˙2. As a result, (2.10) becomes 2 6 6 6 4 SM11dtd − SM12 −Z1M12− Z2M˙22 −M22 0 3 7 7 7 5 x = 2 6 6 6 4 S f1 Z1f1+ Z2f˙2 f2 f3 3 7 7 7 5 d s a v . (2.11)
It is worth to note that the existence of the left-inverse P2−1 of P2guarantees that the step of transforming system (2.6) via (2.7) to (2.11) does not alter the solution set of sys-tem (2.6). Furthermore, the number of scalar differential equations has been reduced from r to d . Continuing this reduction process leads us to the following algorithm.
2.2. Time Varying Differential-Algebraic Equations 18 Algorithm 2.1 Reformulation algorithm for the DAE (2.3a)
1: Set i = 0 and let E0= E , A0= A, f0= f , r0= r , a0= a.
2: Determine a pointwise unitary function P1 as in Lemma 2.14to bring the DAE ¡Ei d dt− A i¢ x = fito the form 2 4 M11dtd − M12 −M22 0 3 5x = 2 4 f1 f2 f3 3 5 ri ai vi , (2.12)
where the functions M11, M22have pointwise full row rank.
3: if rank£MT 11 M22T
⁄T
= ri+ aithen STOP with the resulting system (2.12),
4: else proceed to 5.
5: Determine the operator P2as in (2.9) and apply it to system (2.12) results in 2 6 6 6 4 SM11dtd − SM12 −Z1M12− Z2M˙22 −M22 0 3 7 7 7 5 x = 2 6 6 6 4 S f1 Z1f1+ Z2f˙2 f2 f3 3 7 7 7 5 di si ai vi . 6: Increase i by 1, set Ei:= 2 6 6 6 4 SM11 0 0 0 3 7 7 7 5 , Ai:= 2 6 6 6 4 SM12 Z1M12+ Z2M˙22 M22 0 3 7 7 7 5 , fi= 2 6 6 6 4 S f1 Z1f1+ Z2f2˙ f2 f3 3 7 7 7 5 ,
and repeat the process from 2.
7: end if
Since ri +1= ri− si, Algorithm2.1terminates after a finite number of iterations. This guarantees the existence of the so-called strangeness index of the DAE (2.3a) defined byµ = min{i ∈ N0, ri= ri +1}. We also callµ the strangeness index of the pair (E, A).
Theorem 2.16. Consider the DAE (2.3a) and assume that its strangeness indexµ is well defined. Then, the DAE (2.3a) has the same solution set as the resulting DAE, which we denote by 2 4 ˆ M11dtd − ˆM12 − ˆM22 0 3 5x = 2 6 4 ˆ f1 ˆ f2 ˆ f3 3 7 5, dµ aµ vµ (2.13) where • ˆ M11 ˆ M22 ‚
has pointwise full row rank. The functions ˆf2, ˆf3depend on f , ˙f , . . . , f(µ), while the function ˆf1depends only on f .
The quantities dµ, aµ, vµ and uµ:= n − dµ− aµare called the characteristic quan-tities of the DAE (2.3a). In particular, uµis the number of undetermined variables con-tained in the state vector-valued function x. Following the notation in [75], we also call (2.13) the strangeness-free formulation of the DAE (2.3a).
As a direct result, the solvability of the DAE (2.3a) is analyzed in the following corollary.
Corollary 2.17. Consider the DAE (2.3a) and assume that its strangeness indexµ is well defined. Then, the following assertions hold.
i) The DAE (2.3a) is solvable if and only if in (2.13), one has either ˆf3= 0 or vµ= 0. ii) The initial condition x0is consistent if and only if in addition, − ˆM22(0)x0= ˆf2(0). In this case, (2.3a) has the same solution set as the underlying linear system
• ˆ M11 − ˆM22 ‚ ˙ x =" ˆM12˙ˆ M22 # x + " ˆf1 ˙ˆf2 # . (2.14)
iii) Furthermore, (2.3a) is regular if and only if in addition uµ= 0. If this is the case, (2.14) becomes an ODE, which is often called an underlying ODE in the literature, see [6,22,82].
Remark 2.18. Whereas other index concepts, such as differentiation index [22], per-turbation index [66], or tractability index [82], aiming at the resulting system is ei-ther an underlying ODE or an inherent ODE, the goal of the strangeness index is the strangeness-free formulation (2.13) and the underlying linear system (2.14). Conse-quently, the strangeness index is suitable for general DAEs, which can be underdeter-mined or overdeterunderdeter-mined.
The validity of Algorithm2.1is illustrated by the next example.
Example 2.19. We apply Algorithm2.1to the following DAE •−t t2 −1 t ‚ ˙ x +•1 0 0 1 ‚ x =•0 0 ‚ . (2.15)
First scaling the system with P1=
• 0 −1 −1 t ‚ yields system (2.12) • d dt −t d dt− 1 −1 t ‚ x =•0 0 ‚ . (2.16)
Then, applying the operator
P2=•1 d dt 0 1 ‚
to (2.16) usingdtdt = 1 + tdtd, we obtain the strangeness-free formulation (2.13) • 0 0 −1 t ‚ x =•0 0 ‚ .
The strangeness-index is µ = 1, and the characteristic invariants are dµ = 0, aµ = 1, vµ= 1, uµ= 1.
In order to understand the effect of the reformulation procedure performed in Al-gorithm2.1to DDAEs, we apply it to a modification of the DAE (1.3) including a
func-2.2. Time Varying Differential-Algebraic Equations 20
tion parameter that reads
E (t ) ˙x(t ) = A(t)x(t) + T (t)λ(t) + f (t), (2.17a) together with an initial vector
x(0) = x0. (2.17b)
Hereλ :I→ Cnand T :I→ Cm,n. We further assume that E , A, T and f are sufficiently smooth.
The smoothness comparison between the function parameterλ and the state variable x gives rise to the system classification as follows.
Definition 2.20. The parameter dependent DAE (2.17a) is called:
i) retarded if for any continuous functionλ, there exists a solution x to the IVP (2.17). Since the solution x is continuously differentiable, formally, we will say x is smoother thanλ.
ii) neutral if for any continuously differentiable functionλ, there exists a solution x to the IVP (2.17). Formally, we will say x is at least as smooth asλ.
iii) The remaining case, whereλ must be at least two times continuously differen-tiable to guarantee the existence of a solution x to (2.17), is called advanced. This classification leads to different forms of the resulting DAE (2.13) as in the fol-lowing lemma.
Lemma 2.21. Suppose that the parameter dependent DAE (2.17a) is not advanced and the strangeness indexµ is well-defined for the function pair (E, A) of (2.17a). Consid-ering T (t )λ(t) + f (t) as a new inhomogeneity, then Algorithm2.1 applied to (2.17a) results in a system 2 4 ˆ E1(t ) 0 0 3 5x(t ) =˙ 2 4 ˆ A1(t ) ˆ A2(t ) 0 3 5x(t ) + 2 4 ˆ T10(t ) ˆ T20(t ) ˆ T30(t ) 3 5λ(t) + 2 4 0 0 ˆ T31(t ) 3 5 ˙λ(t) + 2 6 4 ˆ f1(t ) ˆ f2(t ) ˆ f3(t ) 3 7 5, (2.18)
where the matrix function •ˆ
E1 ˆ A2
‚
is of pointwise full row rank. In addition, if (2.17a) is of retarded type then ˆT20= 0 and ˆT31= 0.
Proof. When applying the strangeness-free formulation to the DAE (2.17a), the as-sumption that the system is not advanced ensures that all algebraic constraints of (2.17a) have the form
0 = eA2(t )x(t ) + eT2(t )λ(t) + ˜f2(t ),
for some matrix functions Ae2, Te2, ˜f2. On the other hand, all differential equations of (2.17a) have the form
e
E1(t ) ˙x(t ) = eA1(t )x(t ) + eT1(t )λ(t) + ˜f1(t ),
for some matrix functionsEe1, Ae1, Te1, ˜f1. Moreover, consistency conditions forλ and for the inhomogeneity of the DAE (2.17a) can only arise from one of the following three sources:
i) Adding an algebraic equation to another algebraic equation; ii) Adding a differential equation to another differential equation;
iii) Adding the derivative of some algebraic equation to a differential equation. Therefore, the consistency condition for the inhomogeneity of (2.17a) does not contain derivatives of λ of order bigger than one. This means that Algorithm2.1 applied to (2.17a) results in the DAE of the form (2.18).
2.3 High-order Differential-Algebraic Equations
As will be seen later in Example3.5, even though the DDAE (1.2a) is of first order, it is necessary to study high-order DAEs, since a first order DDAE may contain a hidden high-order DAE. Therefore, in this section we discuss high-order DAEs of the form
Ak(t )x(k)(t ) + ··· + A1(t ) ˙x(t ) + A0(t )x(t ) = f (t), (2.19a) on the time interval t ∈I= [0, tf) ⊂ R, where Ai :I→ Cm,n, for i = 1,...,k, Ak 6= 0, and f :I→ Cm. Similar to the first order case, we investigate the solvability of an IVP consisting of (2.19a) and the initial conditions
x(k−1)(0) = x0(k−1), . . . , ˙x(0) = ˙x0, x(0) = x0. (2.19b) The solution concept for high-order DAEs of the form (2.19a) is stated in the next defi-nition.
Definition 2.22. A function x :I→ Cnis called
i) a (classical) solution of (2.19a) if x ∈ Ck(I,Cn) and x satisfies (2.19a) pointwise, ii) a (classical) solution of the IVP (2.19) if x is a solution of (2.19a) and satisfies (2.19b). We introduce X0:= [x0(k−1),T· · · x0T]T ∈ Cknas an initial vector of the IVP (2.19).
iii) An initial vector X0is called consistent to system (2.19a) if the IVP (2.19) has a solu-tion.
iv) System (2.19a) is called solvable if it has at least one solution. It is called regular if it is solvable and for any consistent initial vector X0, the corresponding IVP (2.19) has a unique solution.
First order DAEs (k = 1) play a crucial role in various applications in science and en-gineering. Numerous investigations including both theoretical and numerical aspects are well known, see [6,22,75] and references therein. Another case of (2.19a) where k = 2 arise in various applications, ranging from constrained mechanical systems, see e.g. [43,53,99,108], to electrical and electro-mechanical systems [8,9], heterogeneous systems [102], and traveling waves [32,89], etc.
Even though special cases (k = 1, 2) of system (2.19a) are considered, there are not many available results about the general case, see [93]. More important, it has been observed, see e. g., [3,93,117,123,139] that the classical approach of ordinary differ-ential equations, which transforms system (2.19a) into a first order DAE by introducing new variables that represent derivatives of x(t ), may lead to a number of mathemat-ical difficulties, for example, unnecessary requirements on the smoothness of the in-homogeneity, or even failure of numerical methods applied to a resulting first order system. The direct treatment for (2.19a) therefore aims at a suitable reformulation that
