Computer-Aided stability analysis of differential-algebraic equations

Copyright © IFAC Nonlinear Control Systems, Stuttgart, Germany, 2004

ELSEVIER

IFAC PUBLICATIONS www.elsevier.com/locate/ifac

COMPUTER-AIDED STABILITY ANALYSIS OF DIFFERENTIAL-ALGEBRAIC EQUATIONS Christian Ebenbauer * Frank Allgower *

* Institute for Systems Theory in Engineering (1ST),

University of Stuttgart, Germany {ce,allgower}~ist.uni-stuttgart.de

Abstract: The purpose of this paper is twofold. In the first part, a simple stability result for time-invariant nonlinear differential-algebraic equations (DAE) with higher index is derived in terms of a new Lyapunov dissipation inequality. In the second part, it is shown that this stability result can be used for a computer-aided stability analysis of polynomial DAEs. Copyright © 2004 IFAC I
1. INTRODUCTION

Le., it allows a computer-aided stability analysis based on the sum of squares decomposition and on semidefinite programming (Parrilo, 2000). This is investigated in the second part of this paper. The remaining paper is organized as follows. In Section 2, the problem formulation is given as well as the new Lyapunov dissipation inequality is introduced. In Section 3, the computational aspect of the Lyapunov dissipation inequality for polynomial DAEs is investigated. Section 4 contains two examples, and finally, conclusions are given in Section 5.

Differential-algebraic equations (DAEs) (or descriptor systems, singular systems), Le., differential equations coupled with algebraic equations, have many applications in various fields of systems theory and control. For example in electrical circuits simulation, nonholonomic mechanical systems, constrained dynamical systems, or in switched systems theory and optimal control. Although there is a rapidly increasing interest in DAEs, there are only a few stability results about DAEs available (Hanke and Marz, 1996; Miiller, 1998; Rabier and Rheinboldt, 1994). The stability problem of DAEs is strongly related with the problem of stability on a constraint manifold. By restricting the stability problem to this manifold, the stability of an equilibrium point can be defined in the sense of Lyapunov by restricting the perturbed solutions to this manifold. Apart from this standard approach, a simple stability result in terms of a new Lyapunov dissipation inequality is derived in the first part of this paper. This approach allows to check stability of timeinvariant nonlinear DAEs, possibly with higher index. 1-1oreover, this condition can be used to prove stability of a polynomial DAE on a computer,

Notations. A function V : !Rn ~ !R is called positive semidefinite, if V(x) ~ 0 Vx and positive definite, if V(O) = 0, V(x) > 0 Vx =I O. FUrthermore, V is called radially unbounded, if V (x ) ~ 00 whenever Ilxll ~ 00. If V is differentiable, then the row vector ~~ (x) = \7V(x) denotes the derivative of V with respect to x. U(x) denotes a neighborhood of a point x E IRn and for x = (Xl, X2), a xl-neighborhood U X1 (X) is defined as the set of all points {(Xl, X2) E !Rn I Xl E U(xd}· Finally, Ilxll denotes the Euclidian norm of x E IR n and the set of polynomials in the variables x = [Xl, ..., Xn] with coefficients in IR is denoted by IR[x].

757

2. A SIMPLE STABILITY RESULT

In the general case, the DAE (1) is not given in semi-explicit form and may be of higher index, i.e., g(Xl' X2) = is not uniquely solvable for X2· However, this situation can be also handled by the idea given above. Before the stability result is derived, some regularity properties have to be assumed:

°

In the first part of the paper, the following problem is considered: Given a time-invariant nonlinear DAE

!(x,x) =0,

(1) Assumption 1. It is assumed that all solutions

is the state and the function f : }Rn X ]Rn ~ ]Rn is assumed to be sufficiently smooth with 1(0,0) = O. where x E

]Rn

x = x(t) of the DAE (1) satisfies f(x(t), x(t)) = 0 for all t ~ 0, and that the solution is sufficiently smooth, i.e., there are no impulses in the solution and the initial conditions are consistent.

Find a Lyapunov dissipation inequality that proves stability of the equilibrium point x = o.

°

Since 1(x (t), X(t)) = is satisfied identically in time, also -9t!(x(t), x(t)) = 0 for all t ~ 0 must be satisfied and also for all higher order derivatives must hold: tti~/(x(t),x(t)) = 0 for all t ~ O. Note that Assumption 1 is not easy to satisfy, but for DAEs obtained from physical systems it is often satisfied. However, define now the following stacked vector containing the first J..L derivatives of ! with respect of time:

The idea behind the Lyapunov dissipation inequality, which is introduced here, is very basic and is illustrated on the following (semi-explicit) DAE:

Xl = f(Xl' X2)

(2)

0= g(Xl' X2),

(3)

with I, 9 appropriately defined. Assume for the moment that for each x 1 there exists a unique X2 such that g(Xl' X2) = 0 and X2 can be expressed as X2 = h(xd, i.e., the system is of (differentiation) index one. Then, there are two ways to check stability. The standard way is the following: First, eliminate X2 from the first equation, i.e., Xl = I(Xl, h(xd), and then check stability in the usual way via finding a Lyapunov function V for Xl = /(Xl, h(Xl)). An alternative way is the following: Instead of substituting X2 from the first equation, one can search for a Lyapunov function V V(xd and a function P = P(Xl' X2) such that

FJL(~) = [ iJ~~~l) ],

(5)

t;i I(x, x)

with c; = (x, X, ... , x(JL+l)), -9t/(x, x) = %x!(x, x)x+ %x/(x, x)x aso. Hence, every solution x = x(t) of (1) satisfies FJL(c;(t)) = 0 for all t ~ o. The following theorem gives a simple stability result for the DAE (1):

Theorem 1. If there exists a radially unbounded differentiable positive definite function V : ]Rn ~ ]R, a function P : ]R(JL+2).n ~ ]R, and a number J..L such that

holds for all (x 1, X2). The basic idea behind inequality (4) is just to check negative semidefiniteness of V on the set of point (Xl, X2) where g(Xl,X2) = O. Otherwise, if g(Xl,X2) =I- 0, one can always find a (penalty) function P such that inequality (4) is satisfied, by just making p(Xl,X2) big enough. The idea of dominating the left side by the right side if g(Xl' X2) =I- 0, is well-known in optimization theory with constraints. There, p plays the role of a Lagrange multiplier. In nonlinear systems and control theory, the use of Lagrange multipliers is quite old. Not exactly but in a similar fashion, such an idea was used, at least implicitly, by Zubov to characterize the region of attraction via a partial differential equation (Hahn, 1967). Recently, this idea was used in (Ebenbauer and Allgower, 2004) to establish a dissipation inequality for nonlinear minimum-phase systems.

\7V(x)x ~ IIFJL(~)112p(~)

(6)

is satisfied for some x-neighborhood Ux(~ = 0), then the equilibrium point x = 0 of system (1) is stable.

Proof. Proof by contradiction. Assume that the dissipation inequality (6) is satisfied but the DAE (1) is not stable. Since any solution x = x (t) satisfies FJJ.(~(t)) = 0 for all t ~ 0, we have \7V(x(t))x(t) ~ 0 in some x-neighborhood. But from this Lyapunov inequality follows that the equilibrium point x = of the DAE (1) is stable.

°

Remark 1. Notice that the variables

~

=

(x, X, ..., x(JL+l)) in the dissipation inequality (6) have to be considered as independent variables which are of purely algebraic nature, i.e., instead of xCi) one could write Yi.

758

Remark 2. In the theory of DAEs, the notion of index plays a fundamental role in theoretical as well as numerical questions. Although no index theory is (explicitly) used here, there is a strong connection between the differentiation index and the stacked vector Fp., which is outlined in the following: The manifold defined by the set of points x such that Fp. (x, X, ... , X(P.+l)) = 0 is solvable is called constraint manifold which contains all solutions x = x(t) of the DAE (1). In (Griepentrog, 1992; Griepentrog et al., 1992), this manifold was used to define the differentiation index. More precisely, the differentiation index J1- is the smallest integer such that for every pair (Yo, 'T/I, Y2, ... , Yp.+d, (Yo, 'T/2, Y2, ... ,Yp.+d willch satisfies Fp.(YO,'T/I,Y2, ... ,Yp.+d = 0 and FJ.L(YO,'T/2,Y2, ···'YJ.L+d = 0 follows that 'T/I = 'T/2, i.e., x is unique. Tills allows a very intuitive geometric interpretation of the differentiation index and establish the connection to the point of view of understanding a DAE as a implicitly defined vector field on a implicitly defined manifold. In particular, under certain regularity assumptions, a unique vector field v on the constraint manifold can be determined such that with x = v(x) follows I(x(t), v(x(t))) = 0, for all t ~ 0 (cf. (Griepentrog, 1992; Rabier and Rheinboldt, 1994)). In contrast to such an explicit characterization of tills vector field v, one can consider the stacked vector FJ.L as an implicit definition of this vector field v by setting FJ.L to zero. The implicit stability analysis via the dissipation inequality (6) corresponds then exactly to an explicit stability analysis of the vector field v on the constrained manifold.

3. POLYNOMIAL DAES AND SUM OF SQUARES In the following, DAEs of form

I(x,x)

=

0

(7)

are considered, where it is assumed that I : IRn x ---+ IRn is a polynomial map, i.e., I E IRn[x).

IR n

The algorithmic problem formulation of proving stability for the DAE (7) is as follows: Given I, J1-, and FJ.L. Find a radially unbounded differentiable positive definite function V and a function p such that the dissipation inequality (6) is satisfied. Of course, to search over all functions V, p is impossible. To implement the problem on a computer, a finite parametrization is necessary. Here, it is assumed that V and p are polynomials up to certain degrees. For example, if global stability has to be proven, following polynomial inequalities must be satisfied:

V(x»O

(8)

\7V(x)x ~ IIFJ.L(~)112p(~), for all ~ = (x, x, ..., X(J.L+I)) E IR(J.L+2).n. Note that V polynomial and positive definite, implies that V is radially unbounded. It is a (NP-)hard computational problem to verify such polynomial inequalities. But with the help of the recently established sum of squares decomposition, it is possible to verify such polynomial inequalities very efficiently. In the following, a short summary of the sum of squares decomposition is given. NIore details can be found in (Parrilo, 2000).

Remark 3. Theorem 1 is also applicable to DAEs of higher-index. In case that the differentiation index of the DAE (1) is known, it is reasonable to choose J1- to be equal to the differentiation index. Theorem 1 can be also applied if the DAE (1) has not an unique solution or if the J1- is not equal to the differentiation index or even if the differentiation index is not well-defined. Finally, notice that one has to satisfy Lyapunov inequality (6) in an x-neighborhood Ux(~ = 0). But under some smoothness assumptions, i.e., ~ = ~(t) becomes small whenever x = x(t) becomes small, one can replace Ux(~ = 0) by U(~ = 0).

A polynomial P E IR[x) of degree d = Ql +... +Qn, E No with real coefficients Ca , i.e., P(x) = ~ al a' f ·f L..." a l +... +a n
P;

Remark 4. Often it is of interest to establish asymptotic stability instead of stability. Then instead of (6), the strict inequality \7V(x)x < IIFJ.L(~),,2p(~) must be satisfied for all nonzero x around some x-neighborhood Ux(~ = 0).

Theorem 2. (Choi et al., 1995). A polynomial E IR[x) of degree 2d has a sum of squares decomposition if and only if there exists a positive semidefinite matrix Q such that

P

(9)

Remark 5. For implicitly defined ODEs, that is I(x,x) = 0, where ~ has full rank, the dissipation inequality (6) turns into \7V(x)x ~ 11/(x, x) 11 2 p(x, x), willch establish stability for implicitly defined ODEs.

where m is the vector of all monomials in of degree less or equal to d, Le., m = [1, Xl, X2, ... , Xn , XIX2, ... , x~). There exists (n~d) such monomials. Xl, ... , X n

759

This representation theorem tells us that all sum of squares polynomials can be parameterized by the (convex) set of positive semidefinite matrices. Obviously, all sum of squares polynomial are positive semidefinite, but the converse is not true. Furthermore, if one has to prove local stability, it is sufficient to check positive semidefiniteness of a polynomial P on a subset in lRn constrained by c(x) 2 0, c E lR[x], i.e., the set {x E lRn : c(x) 2 O}. Then, the following lemma is helpful (Papachristodoulou and Prajna, 2002).

Consider the following DAE in semi-explicit form:

Lemma 1. Let c E lR[x). A polynomial P E lR[x) is positive semidefinite on the set {x E lRn c(x) 2 O} if there exists a positive semidefinite polynomial q E lR[x) such that P(x)c(x )q(x) is positive semidefinite.

Substituting Xl, X2 into the last equation yields

+ Xl + X2 - x~ = 0 X2 + 2xr + 2xr X2 = 0

Xl

X2 -

2 3

3

-Xl

=0.

(10) (11 )

(12)

Instead of differentiating the whole equations, only the algebraic part is differentiated:

X2 - 2XrXI

=

O.

xix~ = O.

(13)

Differentiating again this equation leads to (14)

What makes these results especially interesting froln an engineering point of view is the fact that there exists efficient numerical algorithms, namely semidefinite programming algorithms, which allow to test the conditions of these theorems on a computer. The main advantage of a semidefinite program is that it is a convex optimization problem. Hence, the algorithms are often much more efficient than algorithms from algebraic geometry, which are often used to solve polynomial problems. In (Parrilo, 2000), the connection between the sum of squares decomposition in Theorem 2 and semidefinite programming was established, by showing that the existence of a sum of squares decomposition can be decided by solving a semidefinite program. Finally, the sum of squares decomposition was successfully applied in various tasks in systems and control theory. In particular, in (EI-Samad et al., 2003) a stability analysis under polynomial constraint were performed, but for explicit systems and without considering higher order derivatives of the constraints (so-called hidden constraints), as considered here and which are crucial for DAEs of higher index.

Combining equation (10) to (14) into a reduced stacked vector F2 = F2(XI, X2, X3, X3), following polynomial inequalities are obtained:

av (-Xl -a Xl

- x2

2

+ x3) -

av (2x 3 + 2X 2 x 2) + -a av X3 -a I l X2 x3

:::; IIF2 (XI, X2, X3, x3)11 2 P(XI, X2, X3, X3), for a positive definite V = V(XI, x2, X3) and for all (XI,X2,X3,X3) E lR 4 . With the help of the free available toolbox SOSTOOLS (Prajna et al., 2002), a Lyapunov function of degree four was found, which proves stability of the given polynomial DAE. A Linear Example. A special case of polynomial DAEs are linear DAEs. For linear DAEs, the dissipation inequality (6) with a quadratic Lyapunov function leads directly to a linear matrix inequality, which can be efficiently solved via semidefinite programming.

Consider the linear DAE in semi-explicit form

x = A ll x + A l2 Z 4. EXAMPLES

(15)

0= A 21 X + A 22 Z

A Polynomial Example. Although the polynomial inequalities (8) can be directly formulated as a SDP via the SOS decomposition, it is important to notice that the number of variables grows very fast with the number of differentiations, i.e. J.-L + 1 differentiations leads to (1l+2)·n variables. Such a big number of variables may cause computational problems already for small examples. Hence it is important to reduce the number of unknown variables as much as possible. This can be done by substituting variables, for example if the DAE is given in semi-explicit form, or by differentiating only (the algebraic) parts of the equations. This is demonstrated in the following simple example.

with

which is of index 2. Differentiating the linear DAE (15) leads to

x=

Ailx + A u A 12 Z

0= A 2l A ll x

760

+ A 12 Z

+ A 2l A 12 Z + A 22 Z.

(16)

(17)

Differentiating again these two equations leads to

'x' =

AilX + A u A 12 Z + A 12 Z

(18)

0= A21AilX + A21AllA12Z

(19)

REFERENCES Campell, S. L., C. T. Kelley and K. D. Yeomans (1996). Consistent initial conditions for unstructured higher index DAEs: A computational study. In: Conference on Computational Engineering in Systems Applications (CESA '96), Lille, France. pp. 416 - 421. Choi, M.D., T.Y. Lam and B. Reznick (1995). Sums of squares of real polynomials. In: Proceedings of Symposia in Pure Mathematics 58(2). pp. 103-126. Ebenbauer, C. and F. Allgower (2004). Minimurnphase property of nonlinear systems in terms of a dissipation inequality. In: American Control Conference (A CC), Boston, USA. pp. 1737-1742. El-Samad, H., S. Prajna, A. Papachristodoulou, M. Khammash and J.C. Doyle (2003). Model validation and robust stability analysis of the bacterial heat shock response using sostools. In: Proceedings of the 42nd IEEE Conference on Decision and Control, Maui, Hawaii USA. pp. 3766-3771. Griepentrog, E. (1992). Index reduction methods for differential-algebraic equations. In: Berliner Seminar on Differential-Algebraic Equations. Griepentrog, E., M. Hanke and R. Marz (1992). Toward a better understanding of differentialalgebraic equations. In: Berliner Seminar on Differential- A 1gebraic Equations. Hahn, W. (1967). Stability of Motion. Springer. Hanke, M. and R. Marz (1996). On asymptotics in case of DAE's.. Zeitschrift der angewandten Mathematik und Mechanik (ZAMM) 76, 99102. Muller, P.C. (1998). Stability and optimal control of nonlinear descriptor systems: A survey. Appl. Math. and Comp. Sci. 8, 269-286. Papachristodoulou, A. and S. Prajna (2002). On the construction of Lyapunov functions using the sum of squares decomposition. In: IEEE Conference on Decision and Control, Las Vegas, Nevada. pp. 3482-3487. Parrilo, P. (2000). Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization. PhD thesis. California Institute of Technology. Prajna, S., A. Papachristodoulou and P. A. Parrilo (2002). SOSTOOLS for MATLAB. Available from http://www.cds.caltech.edu/sostools. Rabier, P.J. and W. Rheinboldt (1994). A geometric treatment of implicit differentialalgebraic equations. J. of Differential Equations. 109, 110-146. Sturm, J.F. (2001). SeDuMi 1.05. Available from http://fewcal.kub.nl/ sturm/software/ sedumi.html.

+ A 21 A 12 Z + A 22 Z. Since x and 'x' appear only in the equations (16) and (18) respectively, it makes no sense to add this two equations in the stacked vector, since for any (x, X, z, z) there is a x and 'x such that the equations (16) and (18) are satisfied. Therefore, the following reduced stacked vector is considered:

P2(x, x, z, z, z) = A 12 0 -1 l1 o A 0 A 0 22 [ A21 o A 21 A 12 A 22 o] A 21 A u 0 A 21 Atl 0 A21AuA12 A 21 A 12 A 22

r~~

X

~

The dissipation inequality (6) with V(x) = x T Px, where P has to be a positive definite matrix, can now be written as a linear matrix inequality in the unknowns P, p: x T P±

+ xT Px S 2p . PiP2

(20)

for all (x, X, z, z, z) E }R12. With the help of the free available semidefinite programming toolbox SeDuIvIi (Sturm, 2001), a positive definite matrix P was found, which proves stability of the given linear DAE.

5. CONCLUSIONS In this paper, a simple stability result for timeinvariant nonlinear DAEs with higher index is established. This stability result is stated in terms of a new Lyapunov dissipation inequality. Such a formulation is especially helpful for polynomial DAEs in combination with computational tools like the sum of squares decomposition and semidefinite programming. From a theoretical point of view, the Lyapunov dissipation inequality may also deal as a starting point for controller design formulations. Analogous to the ODE case, synthesis dissipation inequalities in the spirit of passivity or control Lyapunov functions may be formulated.

761