• No results found

Computation of a Canonical Form for Linear Differential-Algebraic Equations

N/A
N/A
Protected

Academic year: 2021

Share "Computation of a Canonical Form for Linear Differential-Algebraic Equations"

Copied!
13
0
0

Loading.... (view fulltext now)

Full text

(1)

Computation of a canonical form for linear

differential-algebraic equations

Markus Gerdin

Division of Automatic Control

Department of Electrical Engineering

Link¨

opings universitet, SE-581 83 Link¨

oping, Sweden

WWW: http://www.control.isy.liu.se

E-mail: gerdin@isy.liu.se

April 7, 2004

AUTOMATIC CONTROL

COM

MUNICATION SYSTEMS

LINKÖPING

Report no.: LiTH-ISY-R-2602

Submitted to Reglerm¨

ote 2004

Technical reports from the Control & Communication group in Link¨oping are available at http://www.control.isy.liu.se/publications.

(2)
(3)

Computation of a canonical form for linear

differential-algebraic equations

Markus Gerdin

April 7, 2004

Abstract

This paper describes how a commonly used canonical form for linear differential-algebraic equations can be computed using numerical software from the linear algebra package LAPACK. This makes it possible to au-tomate for example observer construction and parameter estimation in linear models generated by a modeling language like Modelica.

1

Introduction

In recent years object-oriented modeling languages have become increasingly popular. An example of such a language is Modelica (Fritzson, 2004; Tiller, 2001). Modeling languages of this type make it possible to build models by connecting submodels in a manner that parallels the physical construction. A consequence of this viewpoint is that it is usually not possible to specify in advance what variables are inputs or outputs from a given submodel. A further consequence of this is that the resulting model is not in state space form. Instead the model is a collection of equations, some of which contain derivatives and some of which are static relations. A model of this form is sometimes referred to as a DAE (differential algebraic equations) model. It can be noted that these models are a special case of the so-called behavioral models discussed in, e.g., (Polderman and Willems, 1998). Consider a linear DAE,

E ˙ξ(t) = Jξ(t) + Ku(t) (1)

whereξ(t) is a vector of physical variables and u(t) is the input. E and J are constant square matrices and K is a constant matrix. Linear DAE are also known as descriptor systems, implicit systems and singular systems. For this case it is possible to make a transformation into a canonical form,

 I 0 0 N  Q−1ξ(t) =˙  A 0 0 I  Q−1ξ(t) +B D  u(t) (2)

whereN is a nilpotent matrix. This canonical form has been used extensively in the literature, for example to examine general system properties and obtain the solution (e.g., Ljung and Glad, 2004; Dai, 1989; Brenan et al., 1996), to develop control strategies (e.g., Dai, 1989), to estimate ξ(t) (e.g., Sch¨on et al., 2003), and to estimate parameters (e.g., Gerdin et al., 2003). The transformation itself was treated in (Gantmacher, 1960).

(4)

However, as far as we know it has not been thoroughly studied how this transformation can be computed numerically. The only reference we have found is (Varga, 1992), where a method for computation of the transformation is mentioned briefly. This contribution therefore discusses how the transformation can be computed. The approach here will include pointers to implementations of some algorithms in the linear algebra package LAPACK (Anderson et al., 1999). LAPACK is a is a free collection of routines written in Fortran77 that can be used for systems of linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems. LAPACK is more or less the standard way to solve this kind of problems, and is used by commercial software likeMatlab.

2

The Canonical Form

In this section we provide a proof for the canonical form. We also use the canonical form to show how the solution of a linear DAE can be calculated. The transformations discussed in this section can be found in for example (Dai, 1989), but the proofs which are presented here have been constructed so that the indicated calculations can be computed by numerical software in a reliable manner. How the different steps of the proofs can be computed numerically is studied in Section 4. It can be noted that the system must be regular for the transformation to exist. A linear DAE (1) is regular if det(sE − J) 6≡ 0, that is the determinant is not zero for alls. By Laplace transforming (1), it can be realized that regularity is equivalent to the existence of a unique solution. This is also discussed in for example (Dai, 1989).

The main result is presented in Lemma 3 and Theorem 1, but to derive these results we use a series of lemmas as described below. The first lemma describes how the system matricesE and J simultaneously can be written in triangular form with the zero eigenvalues ofE sorted to the lower right block.

Lemma 1 Consider a system

E ˙ξ(t) = Jξ(t) + Ku(t) (3)

If (3) is regular, then there exist non-singular matricesP1 andQ1 such that P1EQ1=  E1 E2 0 E3  andP1JQ1=  J1 J2 0 J3  (4) whereE1is non-singular,E3 is upper triangular with all diagonal elements zero andJ3 is non-singular and upper triangular.

Note that either the first or the second block row in (4) may be of size zero. Proof. The Kronecker canonical form of a regular matrix pencil discussed in, e.g., (Kailath, 1980, Chapter 6) directly shows that it is possible to perform the transformation (4).

In the case when the matrix pencil is regular, the Kronecker canonical form is also called the Weierstrass canonical form. The Kronecker and Weier-strass canonical forms are also discussed by (Gantmacher, 1960, Chapter 12). The original works by Weierstrass and Kronecker are (Weierstrass, 1867) and (Kronecker, 1890).

(5)

Note that the full Kronecker form is not computed by the numerical software discussed in Section 4. The Kronecker form is here just a convenient way of showing that the transformation (4) is possible.

The next two lemmas describe how the internal variables of the system can be separated into two parts by making the systems matrices block diagonal. Lemma 2 Consider (4). There exist matrices L and R such that

 I L 0 I   E1 E2 0 E3   I R 0 I  =  E1 0 0 E3  (5) and  I L 0 I   J1 J2 0 J3   I R 0 I  =  J1 0 0 J3  . (6)

See (K˚agstr¨om, 1994) and references therein for a proof of this lemma. Lemma 3 Consider a system

E ˙ξ(t) = Jξ(t) + Ku(t) (7)

If (7) is regular, there exist non-singular matricesP and Q such that the trans-formation

P EQQ−1ξ(t) = P JQQ˙ −1ξ(t) + P Ku(t) (8)

gives the system  I 0 0 N  Q−1ξ(t) =˙  A 0 0 I  Q−1ξ(t) +B D  u(t) (9)

whereN is a nilpotent matrix.

Proof. LetP1andQ1 be the matrices in Lemma 1 and define P2=  I L 0 I  (10a) Q2=  I R 0 I  (10b) P3=  E−1 1 0 0 J3−1  (10c) whereL and R are from Lemma 2. Also let

P = P3P2P1 (11a) Q = Q1Q2. (11b) Then P EQ =  I 0 0 J3−1E3  (12) and P JQ =  E−1 1 J1 0 0 I  (13)

(6)

Here N = J3−1E3 is nilpotent since E3 is upper triangular with zero diagonal elements andJ3−1 is upper triangular. J3−1 is upper triangular sinceJ3is.

DefiningA = E1−1J1 finally gives us the desired form (9).

We are now ready to present the main result in this section, which shows how a solution of the system equations can be obtained. We get this result by observing that the first block row of (9) just is a normal state-space description and showing that the solution of the second block row is a sum of the input and some of its derivatives.

Theorem 1 Consider a system

E ˙ξ(t) = Jξ(t) + Ku(t) (14)

If (14) is regular, its solution can be described by ˙

w1(t) = Aw1(t) + Bu(t) (15a)

w2(t) = −Du(t) − m−1X i=1 NiDu(i)(t) (15b)  w1(t) w2(t)  =Q−1ξ(t). (15c)

Proof. According to Lemma 3 we can without loss of generality assume that the system is in the form

 I 0 0 N   ˙ w1(t) ˙ w2(t)  =  A 0 0 I   w1(t) w2(t)  +  B D  u(t) (16a)  w1(t) w2(t)  =Q−1ξ(t). (16b) where w(t) =  w1(t) w2(t)  (17) is partitioned according to the matrices.

Now, ifN = 0 we have that

w2(t) = −Du(t) (18)

and we are done. IfN 6= 0 we can multiply the second row of (16a) with N and get

N2w2˙ (t) = Nw2(t) + NDu(t). (19)

We now differentiate (19) and insert the second row of (16a). This gives w2(t) = −Du(t) − ND ˙u(t) + N2w¨2(t) (20) If N2 = 0 we are done, otherwise we just continue until Nm = 0 (this is true for somem since N is nilpotent). We would then arrive at an expression like

w2(t) = −Du(t) − m−1X

i=1

(7)

u(t) I1(t)

I2(t) I3(t)

R L

Figure 1: A small electrical circuit.

and the proof is complete.

Note that the internal variables of the system may depend directly on deriva-tives of the input. However, it can be noted that the internal variables of physical systems seldom depend directly on derivatives of the input since this would for example lead to the internal variables taking infinite values for a step input. In the common case of no dependence on the derivative of the input, we will have

ND = 0. (22)

We conclude the section with an example which shows what the form (15) is for a simple electrical system.

Example 1 (Canonical form) Consider the electrical circuit in Figure 1. With u(t) as the input, the equations describing the systems are

 00 00 L0 0 0 0     ˙ I1(t) ˙ I2(t) ˙ I3(t)   =  01 −1 −10 0 0 −R 0    II12(t)(t) I3(t)   +  10 1   u(t) (23) Transforming the system into the form (9) gives

 10 00 00 0 0 0    00 01 10 1 0 −1     ˙ I1(t) ˙ I2(t) ˙ I3(t)   =  00 01 00 0 0 1    00 01 10 1 0 −1    II12(t)(t) I3(t)   +   1 L 1 R 1 R u(t) (24) Further transformation into the form (15) gives

˙ w1(t) = L1u(t) (25a) w2(t) = −  1 R 1 R  u(t) (25b)  w1(t) w2(t)  =  00 01 10 1 0 −1    II21(t)(t) I3(t)   . (25c)

(8)

We can here see how the state-space part has been singled out by the transfor-mation. In (25c) we can see that the state-space variablew1(t) is equal to I3(t). This is natural, since the only dynamic element in the circuit is the inductor. The two variables in w2(t) are I2(t) and I1(t) − I3(t). These variables depend directly on the input.

3

Generalized Eigenvalues

The computation of the canonical forms will be performed with tools that nor-mally are used for computation of generalized eigenvalues. Therefore, some theory for generalized eigenvalues will be presented in this section. The theory presented here about generalized eigenvalues can be found in for example (Bai et al., 2000). Another reference is (Golub and van Loan, 1996, Section 7.7).

Consider a matrix pencil

λE − J (26)

where the matrices E and J are n × n with constant real elements and λ is a scalar variable. We will assume that the pencil is regular, that is

det(λE − J) 6≡ 0 with respect to λ. (27) The generalized eigenvalues are defined as thoseλ for which

det(λE − J) = 0. (28)

If the degreep of the polynomial det(λE − J) is less than n, the pencil also has n − p infinite generalized eigenvalues. This happens when rank E < n (Golub and van Loan, 1996, Section 7.7). We illustrate the concepts with an example. Example 2 (Generalized eigenvalues) Consider the matrix pencil

λ  1 0 0 0   −1 0 1 −1  . (29) We have that det  λ  1 0 0 0   −1 0 1 −1  = 1 +λ (30)

so the matrix pencil has two generalized eigenvalues,∞ and −1.

Generalized eigenvectors will not be discussed here, the interested reader is instead referred to for example (Bai et al., 2000).

Since it may be difficult to solve Equation (28) for the generalized eigenval-ues, different transformations of the matrices that simplifies computation of the generalized eigenvalues exist. The transformations are of the form

P (λE − J)Q (31)

with invertible matrices P and Q. Such transformations do not change the eigenvalues since

(9)

One such form is the Kronecker canonical form. However, this form cannot in general be computed numerically in a reliable manner (Bai et al., 2000). For example it may change discontinuously with the elements of the matricesE and J. The transformation which we will use here is therefore instead the generalized Schur form which requires fewer operations and is more stable to compute (Bai et al., 2000).

The generalized Schur form of a real matrix pencil is a transformation such that

P (λE − J)Q (33)

is upper quasi-triangular, that is it is upper triangular with some 2 by 2 blocks, corresponding to complex eigenvalues, on the diagonal. P and Q are orthogonal matrices. The generalized Schur form can be computed with the LAPACK commands dgges or sgges. These commands also give the possibility to sort certain generalized eigenvalues to the lower right. An algorithm for ordering of the generalized eigenvalues is also discussed by (Sima, 1996). Here we will use the possibility to sort the infinite generalized eigenvalues to the lower right.

The generalized Schur form discussed here is also called the generalized real Schur form, since the original and transformed matrices only contain real ele-ments.

4

Computation of the Canonical Forms

The discussion in this section is based on the steps of the proof of the form in Theorem 1. We therefore begin by examining how the diagonalization in Lemma 1 can be performed numerically.

The goal is to find matrices P and Q such that P (λE − J)Q = λ  E1 E2 0 E3  +  J1 J2 0 J3  (34) where E1 is non-singular, E3 is upper triangular with all diagonal elements zero and J3 is non-singular and upper triangular. This is exactly the form we get if we compute the generalized Schur form with the infinite generalized eigenvalues sorted to the lower right. This computation can be performed with the LAPACK commands dgges or sgges. E1 corresponds to finite generalized eigenvalues and is non-singular since it is upper quasi-triangular with non-zero diagonal elements andE3 corresponds to infinite generalized eigenvalues and is upper triangular with zero diagonal elements. J3 is non-singular, otherwise the pencil would not be regular.

The next step is to compute the matricesR and L in Lemma 2, that is we want to solve the system

 I L 0 I   E1 E2 0 E3   I R 0 I  =  E1 0 0 E3  (35a)  I L 0 I   J1 J2 0 J3   I R 0 I  =  J1 0 0 J3  . (35b)

(10)

yields  E1 E1R + E2+LE3 0 E3  =  E1 0 0 E3  (36a)  J1 J1R + J2+LJ3 0 J3  =  J1 0 0 J3  (36b) which is equivalent to the system

E1R + LE3=−E2 (37a)

J1R + LJ3=−J2. (37b)

Equation (37) is a generalized Sylvester equation (K˚agstr¨om, 1994). The gener-alized Sylvester equation (37) can be solved from the linear system of equations (K˚agstr¨om, 1994)  In⊗ E1 ET 3 ⊗ Im In⊗ J1 JT 3 ⊗ Im   vec(R) vec(L)  =  − vec(E2) − vec(J2)  . (38)

HereInis an identity matrix with the same size asE3andJ3,Imis an identity matrix with the same size as E1 and J1, ⊗ represents the Kronecker product and vec(X) denotes an ordered stack of the columns of a matrix X from left to right starting with the first column.

One way to solve the generalized Sylvester equation (37) is thus to use the linear system of equations (38). This system can be quite large, so it may be a better choice to use specialized software such as the LAPACK routines stgsyl or dtgsyl.

The steps in the proof of Lemma 3 and Theorem 1 only contain standard matrix manipulations, such as multiplication and inversion. They are straight-forward to implement, and will not be discussed further here.

5

Summary of the computations

In this section a summary of the steps to compute the canonical forms is pro-vided. It can be used to implement the computations without studying Section 4 in detail. The summary is provided as a numbered list with the necessary com-putations.

1. Start with a system

E ˙ξ(t) = Jξ(t) + Ku(t) (39)

that should be transformed into the form  I 0 0 N  Q−1ξ(t) =˙  A 0 0 I  Q−1ξ(t) +B D  u(t) (40) or ˙ w1(t) = Aw1(t) + Bu(t) (41a) w2(t) = −Du(t) − m−1X i=1 NiDu(i)(t) (41b)  w1(t) w2(t)  =Q−1ξ(t). (41c)

(11)

2. Compute the generalized Schur form of the matrix pencilλE − J so that P1(λE − J)Q1=λ  E1 E2 0 E3  +  J1 J2 0 J3  . (42)

The generalized eigenvalues should be sorted so that diagonal elements of E1 contain only non-zero elements and the diagonal elements of E3 are zero. This computation can be made with one of the LAPACK commands dgges and sgges.

3. Solve the generalized Sylvester equation (43) to get the matricesL and R.

E1R + LE3=−E2 (43a)

J1R + LJ3=−J2. (43b)

The generalized Sylvester equation (43) can be solved from the linear equation system (44) or with the LAPACK commands stgsyl or dtgsyl.

 In⊗ E1 ET 3 ⊗ Im In⊗ J1 JT 3 ⊗ Im   vec(R) vec(L)  =  − vec(E2) − vec(J2)  . (44)

Here In is an identity matrix with the same size as E3 and J3, Im is an identity matrix with the same size as E1 and J1, ⊗ represents the Kronecker product and vec(X) denotes an ordered stack of the columns of a matrixX from left to right starting with the first column.

4. We now get the form (40) and (41) according to P =  E−1 1 0 0 J3−1   I L 0 I  P1 (45a) Q = Q1  I R 0 I  (45b) N = J−1 3 E3 (45c) A = E−1 1 J1 (45d)  B D  =P K. (45e)

6

Conclusions

We have in this paper examined how a commonly used canonical form for linear differential-algebraic equations can be computed using numerical software. As discussed in the introduction, it is possible to for example estimate the states or unknown parameters using this canonical form. Models generated by a modeling language like Modelica are described as differential-algebraic equations. For linear Modelica models it is therefore possible to automate the procedures of for example observer construction and parameter estimation.

7

Acknowledgments

This work was supported by the Swedish Research Council and by the Founda-tion for Strategic Research (SSF).

(12)

References

Anderson, E., Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney and D. Sorensen (1999). LAPACK Users’ Guide. 3. ed. Society for Industrial and Applied Mathematics. Philadelphia.

Bai, Z., J. Demmel, J. Dongarra, A. Ruhe and H. van der Vorst (2000). Templates for the Solution of Algebraic Eigenvalue Problems, A Practical Guide. SIAM. Philadelphia.

Brenan, K.E., S.L. Campbell and L.R. Petzold (1996). Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. Classics In Ap-plied Mathematics. SIAM. Philadelphia.

Dai, L. (1989). Singular Control Systems. Lecture Notes in Control and Infor-mation Sciences. Springer-Verlag. Berlin, New York.

Fritzson, Peter (2004). Principles of Object-Oriented Modeling and Simulation with Modelica 2.1. Wiley-IEEE. New York.

Gantmacher, F.R. (1960). The Theory of Matrices. Vol. 2. Chelsea Publishing Company. New York.

Gerdin, M., T. Glad and L. Ljung (2003). Parameter estimation in linear differential-algebraic equations. In: Proceedings of the 13th IFAC sympo-sium on system identification. Rotterdam, the Netherlands. pp. 1530–1535. Golub, G.H. and C.F. van Loan (1996). Matrix Computations. 3 ed. The John

Hopkins University Press. Baltimore and London.

K˚agstr¨om, B. (1994). A perturbation analysis of the generalized sylvester equa-tion. Siam Journal on Matrix Analysis and Applications 15(4), 1045–1060. Kailath, T. (1980). Linear Systems. Information and Systems Sciences Series.

Prentice Hall. Englewood Cliffs, N.J.

Kronecker, L. (1890). Algebraische reduction der schaaren bilinearer formen. S.-B. Akad. Berlin pp. 763–776.

Ljung, L. and T. Glad (2004). Modellbygge och simulering. Studentlitteratur. In Swedish.

Polderman, J.W. and J.C. Willems (1998). Introduction to Mathematical Sys-tems Theory: a behavioral approach. Number 26 in: Texts in Applied Math-ematics. Springer-Verlag. New York.

Sch¨on, T., M. Gerdin, T. Glad and F. Gustafsson (2003). A modeling and filter-ing framework for linear differential-algebraic equations. In: Proceedfilter-ings of the 42nd IEEE Conference on Decision and Control. Maui, Hawaii, USA. pp. 892–897.

Sima, V. (1996). Algorithms for linear-quadratic optimization. Dekker. New York.

(13)

Tiller, M. (2001). Introduction to Physical Modeling with Modelica. Kluwer. Boston, Mass.

Varga, A. (1992). Numerical algorithms and software tools for analysis and modelling of descriptor systems. In: Prepr. of 2nd IFAC Workshop on System Structure and Control, Prague, Czechoslovakia. pp. 392–395. Weierstrass, K. (1867). Zur theorie der bilinearen und quadratischen formen.

References

Related documents

Much of the theory of weighted Bergman spaces have been developed with weights locally bounded from above and below, and the problem of removable singularities have the same solution

svarar på enkäten två gånger eller inte alls jobbar med barn eftersom enkäter lades i fikarum, bara 2% var män. Högsta ålder på deltagarna framkommer inte.

As has been explained above, linear diophantine equation in two variables can be solved by euclidean algorithm and extended euclidean algorithm. When the solution has been searched

Department of Electrical Engineering Linköping University.. SE-581

The principal findings of the present study were that: (1) prevalences of symptoms of CMD among male European professional footballers ranged from 3 % for smoking to 37 %

Syftet med studien är att undersöka hur lärare upplever att man bör bemöta elever i läs- och skrivsvårigheter/dyslexi i undervisningen i ämnet engelska, för en

Is it feasible to replace human bug assignment and fault localization with machine learning techniques in large scale industry SW development projects.. RQ2: Which machine

Group classification of linear Schrödinger equations.. by the algebraic method