Document [original]

Regular solutions of DAE hybrid systems

and regularization techniques ∗

Peter Kunkel †Volker Mehrmann ‡

18.04.17

Abstract. The solvability and regularity of hybrid differential-algebraic systems (DAEs) is

studied, and classical stability estimates are extended to hybrid DAE systems. Different reasons

for non-regularity are discussed and appropriate regularization techniques are presented. This

includes a generalization of Filippov regularization in the case of so-called chattering. The results

are illustrated by several numerical examples.

Keywords. Differential-algebraic equation, hybrid system, switched system, index reduction,

existence and uniqueness of solutions, Filippov regularization, strangeness index, chattering.

AMS(MOS) subject classification. 34K34, 34K28, 65L80, 65L05, 34K32

1 Introduction

Differential-algebraic equations (DAEs) are widely used in the modeling process of complex dy-

namical systems, see, e. g., [3, 9]. In many applications the mathematical system equations switch

between different models. One then speaks of switched or hybrid systems, see, e. g., [4, 6, 22].

Classical application areas of switched systems are robot manipulators [17], mechatronics [16],

mechanical systems with dry friction [10, 20], or automatic gear-boxes [15]. Other applications

include electronic circuits where different device models are used for different frequency ranges,

or switching elements like diodes are used, see [22], process control in chemical engineering [1],

and in control systems [24, 26], where the value of a control switches between different operation

modes. See also [12, 18] for further examples. For an overview of modeling, analysis, simulation

and control of hybrid systems, see, e. g., [4, 23, 24, 25].

In this paper we study the solvability of general nonlinear hybrid systems of differential-algebraic

equations and extend the work of [14, 15, 24, 25]. We characterize the existence of regular solutions

and when no regular solutions exist. We also present stability estimates for general hybrid systems

of DAEs and derive regularization techniques that allow the numerical integration.

A hybrid system of DAEs on a (nontrivial) closed interval I= [t, t]⊆Rconsists of model

equations describing the different modes of the hybrid system, switching conditions describing

when a mode loses its validity, and transition functions describing the initiation of a new mode.

Hence, a hybrid system of DAEs with Jmodes is described as follows. The model equations are

given by nonlinear DAEs of the form

Fj(t, xj,˙xj)=0, j = 1, . . . , J,

Fj:I×Dj

x×Dj

˙x→Rnj,(1)

∗Supported through the Research-in-Pairs Program at Mathematisches Forschungsinstitut Oberwolfach.

†Mathematisches Institut, Universit¨at Leipzig, Augustusplatz 10, D-04109 Leipzig, Fed. Rep. Germany,

[email protected]. Supported by Deutsche Forschungsgemeinschaft under grant no. KU964/7-1.

‡Institut f¨ur Mathematik, MA 4-5, Technische Universit¨at Berlin, D-10623 Berlin, Fed. Rep. Germany,

[email protected]. Supported by European Research Council through the Advanced Grant MODSIM-

CONMP.

where Dj

x,Dj

˙x⊆Rnjdenote the open domains of definition of the variables xjand ˙xjin the

arguments of Fj.

We assume each of these mode systems to be regular in the sense that they satisfy Hypothesis 4.2

of [19], see Section 2 below. In particular, for each mode we have a set Lj

µjdescribing the possible

consistent initial conditions for this mode.

The switching conditions have the form

l(t, xj)=0, l = 1, . . . , Lj, j = 1, . . . , J,

l:I×Dj

x→R,(2)

together with the so-called feasible regions

Fj={(t, xj)∈I×Dj

x|Sj

l(t, xj)≥0, l = 1, . . . , Lj}(3)

in which the solution is required to lie if the system is in the jth mode.

The transition functions have the form

Tkj

l:I×Dj

x→Lk

µk, l = 1, . . . , Lj, j = 1, . . . , J, (4)

where k=σj(l), with σj:{1, . . . , Lj}→{1, . . . , J}, is the new mode when mode jwas terminated

due to a sign change in the switching function Sj

For our analysis we will assume that all functions are sufficiently smooth, i. e., sufficiently often

continuously differentiable. In particular, we assume that all occurring derivatives exist.

The paper is organized as follows. In Section 2 we recall the basics of differential-algebraic

equations and in Section 3 we discuss the regularity of hybrid DAE systems. The extension of

classical stability estimates to hybrid DAE systems is treated in Section 4. Possible non-regular

behavior and regularization techniques for hybrid DAEs, including a generalization of the Filippov

regularization, are presented in Section 5 and some numerical examples are presented in Section 6.

We conclude with a summary in Section 7.

2 Preliminaries

In order to introduce the concept of regularity for solutions of hybrid DAE systems, the DAEs

of the single modes must be regular. To recall the regularity of general DAEs, let Fdescribe a

nonlinear DAE F(t, x, ˙x)=0,

F:I×Dx×D˙x,Dx,D˙x⊆Rnopen, (5)

as it is given in each mode of (1). Following an idea of Campbell, see, e. g., [7], we consider

so-called derivative array equations

Fµ(t, x, y)=0, y = ( ˙x, . . . , x(µ+1)),

Fµ:I×Dx×Dy→R(µ+1)n,Dx⊆Rn,Dy⊆R(µ+1)nopen,

defined by stacking Fand its first µderivatives of Fon top of each other, thus

Fµ=





dt F

dt )µF







According to [19], regularity of the DAE (5) can be guaranteed by requiring the following hypoth-

esis, in which by Fµ;zwe denote the Jacobian of Fµwith respect to the variable z.

Hypothesis 1 Consider a DAE of the form (5). There exist integers µ,a, and dsuch that the

set

Lµ=F−1

µ({0})

associated with Fis nonempty and such that for every point (t0, x0, y0)∈Lµ, there exists a

(sufficiently small) neighborhood Uin which the following properties hold:

1. We have rank Fµ;y= (µ+1)n−aon Lµ∩Usuch that there exists a smooth matrix function Z2

of size (µ+ 1)n×aand pointwise maximal rank, satisfying ZT

2Fµ;y= 0 on Lµ∩U.

2. We have rank ZT

2Fµ;x=aon Usuch that there exists a smooth matrix function T2of size

n×d,d=n−a, and pointwise maximal rank, satisfying ZT

2Fµ;xT2= 0 on U.

3. We have rank F˙xT2=don Usuch that there exists a smooth matrix function Z1of size

n×dand pointwise maximal rank, satisfying rank ZT

1F˙xT2=don U.

If a DAE system satisfies Hypothesis 1, then the smallest µfor which this is the case is called

the strangeness index of Fand systems with µ= 0 are called strangeness-free.

We will utilize Hypothesis 1 in the following form. Given a solution x∗:I→Rnof the derivative

array system Fµ= 0, in the sense that there exists a continuous path (t, x∗(t),P(t)) ∈Lµfor all

t∈Iwith ˙x∗=P[In0 ]T, Hypothesis 1 guarantees the existence of (smooth) functions ˆ

F1,ˆ

defining a so-called reduced DAE

F1(t, x, ˙x)=0,(ddifferential equations)

F2(t, x) = 0,(aalgebraic equations) (6)

with x∗being a solution of (6). Moreover, the reduced DAE satisfies Hypothesis 1 with µ= 0

and, therefore, possesses solutions for sufficiently small perturbations of the initial value x∗(t0)

within the set {x0∈Rn|ˆ

F2(t0, x0) = 0}.

A (smooth) transformation of the variable xaccording to

x=Qx1

x2,(7)

where Q:I→Rn×nis pointwise orthogonal, then yields a decoupled DAE of the form

˙x1=L(t, x1),(ddifferential equations)

x2=R(t, x1),(aalgebraic equations) (8)

with appropriate (smooth) functions L,R. For details, see [19].

Since for given ˆ

t, every (t0, x0, y0)∈Lµin a neighborhood Uof (ˆ

t, x∗(ˆ

t),P(ˆ

t)) fixes a solution of

the reduced DAE, we have a flow Φt,t0:U→I×Rnfor all tin a sufficiently small neighborhood

of t0. The flow maps (t0, x0, y0) to (t, x(t)), where x(t) is the value of the corresponding solution

at t. Note that the existence of a (global) solution on the whole of Iguarantees the existence of

Φt,t0for all t∈Iprovided Uis chosen sufficiently small.

3 Regular solutions of hybrid systems of DAEs

To study the regularity of hybrid systems of DAEs, we assume that every mode in (1) is regular

in the sense that Fjsatisfies Hypothesis 1 with some characteristic values µj, aj, djand a set

µjof possible consistent initial conditions, as it was introduced in Hypothesis 1. Using then the

decoupled index reduced formulation (8) for each mode,

˙xj

1=Lj(t, xj

1), xj

2=Rj(t, xj

1), j = 1, . . . , J, (9)

one obtains transformed switching conditions (2) which have the form

l(t, xj

1)=0, l = 1, . . . , Lj, j = 1, . . . , J

after eliminating the variables x2with the help of the algebraic equations in (9). Accordingly, the

transition functions can be transformed to

Tkj

l(t, xj

1) = (t, xk

1), l = 1, . . . , Lj, j = 1, . . . , J, k ∈ {1, . . . , J}.

Remark 2 If a switching condition is given in the more general form

l(t, xj,˙xj)=0,

then it should reduce to the form

l(t, xj

1, xj

2,˙xj

1)=0

before the elimination of xj

2. In particular, there should be no argument ˙xj

2, since the presence

of ˙xj

2would result in the need of a further differentiation of Fj. This requirement corresponds

to the possibility in the ODE case to eliminate differentiated arguments with the help of the

differential equation.

In order to define a regular solution of a hybrid system on a given fixed interval [t0, T], we must

define regular solutions of the single modes (which is already obtained by requiring Hypothesis 1

for each mode) as well as a regular switching behavior. Thus, we first discuss the regularity of a

piece of the overall solution that corresponds to a specific mode. In the following, we call these

pieces of the solution regular arcs.

Definition 3 Consider a hybrid DAE (1) in the interval I= [t, t]and let j∈ {0, . . . , J}. A

function x∈C1(I,Rnj)is called a regular arc (of mode j)if

1. there exists z= (x,y)such that (t, z)∈Lj

µjand (t, x(t)) = Φj

t,t(t, zj)for all t∈I;

2. the initial value is feasible according to (t, x)∈Fjand there exists an εj>0and a neigh-

borhood Uof (t, z)such that

Φj

t,t0(t0, z0)∈

◦

for all (t0, z0)∈Uand all t∈(t0, t0+εj], where

◦

Fjdenotes the interior of the set Fj;

3. there exists ˆ

l∈ {1, . . . , Lj}such that

(a) the end point tis a regular solution of the one-dimensional problem (Sj

l◦Φj

t,t)(t, z)=0

in the sense that d

dt (Sj

l◦Φj

t,t)(t, z)|t=t6= 0,

(b) strict feasibility holds according to (Sj

l◦Φj

t,t)(t, z)>0for all t∈(t, t),

according to (Sj

l◦Φj

t,t)(t, z)>0for all t∈(t, t]and all l∈ {1, . . . , Lj}with l6=ˆ

Remark 4 The conditions in Definition 3 can be interpreted as follows. Condition 1. just states

that xis solution of the DAE of mode jbelonging to the (consistent) initial condition given by

x(t) = x. Condition 2. requires that the initial state is feasible and the corresponding solution

stays strictly feasible according to

(Sj

l◦Φj

t,t0)(t0, z0)>0, j = 1, . . . , J

for sufficiently small perturbations of the initial condition at least for some small time interval

whose length can be bounded away from zero. Condition 3. requires that the end point is a regular

solution of a unique switching function along the trajectory and that no other switching condition

was satisfied before.

To simplify the notation in the following definition of a regular solution of a hybrid system, we

assume that every mode contains the termination condition t−t= 0 as a switching condition.

Definition 5 Consider a hybrid DAE (1) in the interval I= [t, t]. Given a collection of N∈N

switching points

t= (t0, t1,· · · , tN), t =t0< t1<· · · < tN=t

and a collection of modes

j= (j0, j1,· · · , jN−1), ji∈ {1, . . . , J}for i= 0, . . . , N −1,

a finite collection

x= (xj0

0, xj1

1, . . . , xjN−1

N−1)

of functions xji∈C1([ti, ti+1],Rnji)is called a regular solution of the hybrid system (1) if xjiis

a regular arc of mode jifor every i= 0, . . . , N −1, and the following condition holds:

Let

z= (zj0

0, zj1

1, . . . , zjN−1

N−1)

be the corresponding initial values and let

l= (l0, l1, . . . , lN−1)

be the uniquely defined activated switching functions according to the definition of a regular arc.

Then the transitions satisfy

ji+1 =σji(li)

and

(ti+1, zji+1

i+1 )=(Tji+1,ji

li◦Φji

ti+1,ti)(ti, zji

i), i = 0, . . . , N −2.

Remark 6 A regular solution according to Definition 5 has the following properties. Since we

require a finite switching structure given by tand j, there is a common ε= minj∈{0,...,J}εj>0

for all arcs of a regular solution. By assumption, the activated switching condition in mode jN−1

of the last interval is given by t−t= 0 which leads to the termination of the integration process

for the hybrid system.

The notion of a regular solution is justified by the following result.

Theorem 7 Consider a hybrid system of the form (1) and let xwith corresponding tand jbe

a regular solution of the given hybrid system. Then (1) possesses a regular solution for every

initial condition xj0(t0) = ˜x0with (t0,˜xj0

0,˜yj0

0)∈Lj0

µj0from a sufficiently small neighborhood of

(t0, xj0

0, yj0

0)=(t0, zj0

0)and the final value ˜xjN−1(tN)depends smoothly on ˜x0.

Proof. Since there exists a solution of the hybrid system, there is a flow Φji

t,˜

t(˜

t, ˜zji

i) for all consistent

(˜

t, ˜zji

i) in a neighborhood of (ti, zji

i) and all tin a neighborhood of [ti, ti+1].

By assumption, the switching point ti+1 is a regular solution of

Sji

li(Φji

t,ti(ti, zji

i)) = 0.

Hence, by the implicit function theorem, see, e. g., [21], there exists a (smooth) function Rji

liwith

Sji

li(Φji

t,˜

t(˜

t, ˜zji

i)) ≡0 for ˆ

t=Rji

li(˜

t, ˜zji

i).

Since the sign conditions hold for all switching functions Sjl

l, the so obtained ˆ

tremains the first

switching point that occurs for all consistent (˜

t, ˜zji

i) in a neighborhood of (ti, zji

i). As a conse-

quence, at ˆ

tthe system changes from mode jito mode ji+1 and Φji

t,˜

t(˜

t, ˜zji

i) depends smoothly on

consistent initial values (˜

t, ˜zji

i) in an appropriate neighborhood of (ti, zji

i). Hence,

(˜

t, ˜zji+1

i+1 ) = Tji+1,ji

li(Φji

t,˜

t(˜

t, ˜zji

i)) with ˆ

t=Rji

li(˜

t, ˜zji

depends smoothly on (˜

t, ˜zji

i) in an appropriate neighborhood of (ti, zji

i) and is consistent for the

mode ji+1 due to the required properties of Tji+1,ji

li.

In this way we have shown that small changes in the initial conditions of a mode do not change

the switching that is activated and only yield small changes in the initial condition for the next

mode.

Putting all modes together, we have shown that there exists a neighborhood of (t0, zj0

0) such that

for all consistent (t0,˜zj0

0) from this neighborhood the corresponding initial value problem for (1)

possesses a solution with the same collection of modes and the same activated switching functions

given both by the sequence j. The switching points ˜

tiand initial values ˜zji

idepend smoothly on

the initial value (t0,˜zj0

0). In particular, the final value ˜xjN−1(tN) is a smooth function of the initial

value (t0,˜zj0

0).

Remark 8 It follows from Theorem 7 that we have

(tN, xjN−1) = Ψ(t0, zj0

with a smooth function

Ψ=ΦjN−1

tN,tN−1◦N−2



i=0

(Tji+1,ji

li◦Φji

ti+1,ti),

where

ti+1 =Rji

li(ti, zji

i) = Rji

li◦i−1



m=0

(Tjm+1,jm

lm◦Φjm

tm+1,tm)(t0, zj0

0).

Remark 9 If suitable spaces for a Banach space formulation of an initial value problem for hybrid

DAE systems are needed, as for example in the treatment of optimal control problems for hybrid

systems, then it is necessary to transform the switching structure to a fixed grid. This can for

example be done via a linear transformation

θi: [ti, ti+1]→[i, i + 1], θi(t) = t−ti

ti+1 −ti

+k.

To summarize, we have introduced the concept of regular solutions for hybrid systems of DAEs

and we have shown that such a solution stays regular in a small neighborhood with the same

switching structure and that we obtain a flow of the system for the whole interval Iunder consid-

eration. In the next section we derive a stability estimate for such regular solutions.

4 Stability estimates

For ordinary differential equations

˙x=f(t, x) (10)

there is a well-known stability estimate which states that (under suitable assumptions) the solution

depends Lipschitz-continuously on perturbations in the initial value and in the evaluation of the

right hand side of the ODE, see, e. g., [13, Theorem I.10.3].

Theorem 10 Let x∈C1([t0, t1],Rn)be a solution of (10) in a real interval [t0, t1]and let ˜x∈

C1([t0, t1],Rn)satisfy

(a) k˜x(t0)−x(t0)k ≤ δ,

(b) k˙

˜x(t)−f(t, ˜x(t))k ≤ βfor all t∈[t0, t1],

with some constant L > 0. Then

k˜x(t)−x(t)k ≤ eL(t−t0)δ+1

L(eL(t−t0)−1)βfor all t∈[t0, t1].

In this section we derive a corresponding stability estimate for hybrid systems of DAEs, which

includes possible perturbations in all parts of the involved computations thus generalizing Theo-

rem 7.

We first consider hybrid ODEs and consider a single mode given by (10) together with a regular

arc x∈C1([tj, tj+1],Rn). This arc is governed by a unique switching condition denoted here by

S(t, x)=0.

In particular, we assume that fand Sare defined on a compact neighborhood of

{(t, x(t)) |t∈[tj, tj+1]} ⊆ R×Rn

and that (a) S(tj+1, x(tj+1)) = 0,d

dt S(tj+1, x(tj+1)) 6= 0,

(b) S(t, x(t)) >0 for all t∈(tj, tj+1).

In this way, the terminal point tj+1 is fixed by being the first zero larger than tjof the function

u(t) = S(t, x(t)).(11)

Suppose that ˜x∈C1([˜

tj,˜

tj+1],Rn) is a perturbed solution which satisfies

(a) k˜x(˜

tj)−x(tj)k ≤ δ,

(b) |˜

tj−tj| ≤ τ,

˜x(t)−f(t, ˜x(t))k ≤ βfor all t∈[tj, tj+1],

(d) |S(˜

tj+1,˜x(˜

tj+1))| ≤ σ,

(12)

for sufficiently small δ, τ, β, σ > 0. Our goal is to estimate k˜x(˜

tj+1)−x(tj+1)kand |˜

tj+1 −tj+1|.

Note, however, that the correct terminal point ˜

tj+1 under perturbations may not be the first zero

of the switching function, since under perturbations we cannot guarantee the feasibility of the

initial value.

For the correct choice of the terminal point, we first observe that if τis sufficiently small, then

the solution xcan be extended to ˜

tjeven if ˜

tjis not contained in [tj, tj+1]. Since fis bounded on

its compact domain of definition, say by a constant M, then it holds that

kx(˜

tj)−x(tj)k=kx(tj+s(˜

tj−tj))|1

0k=kR1

0˙x(tj+s(˜

tj−tj))(˜

tj−tj)dsk

≤maxs∈[0,1] kf(tj+s(˜

tj−tj), x(tj+s(˜

tj−tj)))k |˜

tj−tj|ds ≤Mτ. (13)

It follows from this inequality and (12a) that

k˜x(˜

tj)−˜x(tj)k=k˜x(˜

tj)−x(tj)k+kx(˜

tj)−x(tj)k ≤ δ+Mτ, (14)

and Theorem 10 applied to the interval [tj, tj+1] implies that

k˜x(t)−x(t)k ≤ eL(t−tj)(δ+Mτ) + 1

L(eL(t−tj)−1)β, (15)

as long as both functions xand ˜xare defined and satisfy

kf(t, ˜x(t)) −f(t, x(t))k ≤ Lk˜x(t)−x(t)k

with some constant L > 0.

Having studied the influence of a perturbation in the initial state and during the integration

of the system, we now must discuss the influence of the perturbations on the accepted switching

point ˜

tj+1.

Since we only detect switching points by the change of sign of the switching function along the

computed solution, it is natural to assume that both xand ˜xare still defined in a neighborhood

of their final points. This then implies that there exists a sufficiently small γ > 0 such that tj+1 is

the unique zero of udefined by (11) in the interval (tj+1 −γ, tj+1 +γ). Due to (15), the function

˜u(t) = S(t, ˜x(t)) (16)

will also change sign in (tj+1 −γ, tj+1 +γ) for sufficiently small δ, τ, β > 0.

For sufficiently small γ > 0, then

dt S(t, x(t)) 6= 0 for all t∈(tj+1 −γ, tj+1 +γ)

and the estimate

dt S(t, ˜x(t)) −d

dt S(t, x(t))k=kSt(t, ˜x(t)) + Sx(t, ˜x(t))˙

˜x(t)−St(t, x(t)) −Sx(t, x(t)) ˙x(t)k

≤ kSt(t, ˜x(t)) −St(t, x(t))k+kSx(t, ˜x(t))(˙

˜x(t)−f(t, ˜x(t)))k

+kSx(t, ˜x(t))f(t, ˜x(t)) −Sx(t, x(t)) ˙x(t)k

≤D(k˜x(t)−x(t)k+β)

holds for some constant D > 0. This shows that

dt S(t, ˜x(t)) 6= 0 for all t∈(tj+1 −γ, tj+1 +γ)

for sufficiently small δ, τ, β > 0. Hence, the function ˜ufrom (16) has a unique zero in the interval

(tj+1 −γ, tj+1 +γ).

We now assume that ˜

tj+1 ∈(tj+1 −γ, tj+1 +γ) and that σ > 0 is sufficiently small. Then we

consider the nonlinear system of equations

S(t, ˆ

tj+1,ˆxj+1,ˆσ)=0,

defined by ˆ

S(t, ˆ

tj+1,ˆxj+1,ˆσ) = S(t, x(t, ˆ

tj+1,ˆxj+1)) −ˆσ,

where x(t, ˆ

tj+1,ˆxj+1) denotes the (local) solution (10) satisfying the initial condition x(ˆ

tj+1) =

ˆxj+1.

Obviously, the function ˆ

Sis defined in a neighborhood of (tj+1,˜

tj+1, x(˜

tj+1),0) and satisfies

(a) ˆ

S(t1,˜

tj+1, x(˜

tj+1),0)=S(tj+1, x(tj+1,˜

tj+1, x1(˜

tj+1))) = S(tj+1, x(tj+1)) = 0,

(b) ∂

∂t ˆ

S(tj+1,˜

tj+1, x(˜

tj+1),0)=St(tj+1, x(tj+1)) + Sx(tj+1, x(tj+1)) ˙x(tj+1)6= 0.

Thus, we can apply the implicit function theorem and we obtain that locally there exists a func-

tion Zwhich is as smooth as ˆ

Sand satisfies

(a) tj+1 =Z(˜

tj+1, x(˜

tj+1),0),

(b) ˆ

S(Z(ˆ

tj+1,ˆxj+1,ˆσ),ˆ

tj+1,ˆxj+1,ˆσ) = 0 for all (ˆ

tj+1,ˆxj+1,ˆσ).

Since ˆ

S(˜

tj+1,˜

tj+1,˜x(˜

tj+1), S(˜

tj+1,˜x(˜

tj+1)))

=S(˜

tj+1, x(˜

tj+1,˜

tj+1,˜x(˜

tj+1)) −S(˜

tj+1,˜x(˜

tj+1)) = 0

within the validity of the implicit function theorem for sufficiently small δ, τ, β, σ > 0, we conclude

that

tj+1 =Z(˜

tj+1,˜x(˜

tj+1), S(˜

tj+1,˜x(˜

tj+1))).

Hence

|˜

tj+1 −tj+1|=|Z(˜

tj+1,˜x(˜

tj+1), S(˜

tj+1,˜x(˜

tj+1))) −Z(˜

tj+1, x(˜

tj+1),0)|

≤K(k˜x(˜

tj+1)−x(˜

tj+1)k+σ),

with some constant K > 0. Finally, in the same way as in (13), we get

k˜x(˜

tj+1)−x(tj+1)k≤k˜x(˜

tj+1)−x(˜

tj+1)k+kx(˜

tj+1)−x(tj+1)k

≤ k˜x(˜

tj+1)−x(˜

tj+1)k+M|˜

tj+1 −tj+1|.

We summarize this analysis in the following stability estimate for a single mode in a switched

system of ODEs.

Lemma 11 Let all assumptions in the above construction for a single mode of a hybrid system of

ODEs be satisfied. Then, for sufficiently small δ, τ, β, σ > 0there exist constants κδ, κτ, κβ, κσ>0

such that

k˜x(˜

tj+1)−x(tj+1)k+|˜

tj+1 −tj+1| ≤ κδδ+κττ+κββ+κσσ. (17)

Remark 12 The unique zero in (tj+1 −γ, tj+1 +γ) of ˜ufrom (16), however, does not need to

be the first zero of ˜ubeyond ˜

tj. Due to the perturbations it is possible that there exists a zero

immediately after ˜

tj, typically together with the observation that the given initial state seems

to be infeasible. If this is the case, then the implementation of the zero finder for the switching

points must take care to skip such an artificial switching point in order to obtain a robust solver

for hybrid systems.

To extend Lemma 11 from a single regular arc to several regular arcs, we now assume that

x= (xj0

0, xj1

1, . . . , xjN−1

N−1) represents a regular solution of a hybrid system of ODEs and that

x= (˜xj0

0,˜xj1

1,...,˜xjN−1

N−1) is an approximate solution. In particular, we assume that the same

switching structure is valid, i. e., that we have the same modes for xand ˜

x, but instead of

switching points t= (t0, t1,· · · , tN) for xwe may have (possibly) perturbed switching points

t= (˜

t0,˜

t1,· · · ,˜

tN) for ˜

For the transitions, we assume that

(a) (ti+1, xji+1

i+1 (ti+1)) = T(ti+1, xji

i(ti+1)),

(b) k(˜

ti+1,˜xji+1

i+1 (˜

ti+1)) −T(˜

ti+1,˜xji

i(˜

ti+1))k ≤ %, i= 0, . . . , N −2

where for simplicity we have omitted the (uniquely fixed) lower index of the transition functions.

If the transition functions are Lipschitz continuous, then we obtain that

k(˜

ti+1,˜xji+1

i+1 (˜

ti+1)) −(ti+1, xji+1

i+1 (ti+1))k

=k(˜

ti+1,˜xji+1

i+1 (˜

ti+1)) −T(˜

ti+1,˜xji

i(˜

ti+1))k

+kT(˜

ti+1,˜xji

i(˜

ti+1)) −T(ti+1, xji

i(ti+1))k

≤%+Uk(˜

ti+1,˜xji

i(˜

ti+1)) −(ti+1, xji

i(ti+1))k

(18)

with some constant U > 0. Combining inductively the estimates (18) for every transition with

the estimate (17) for every regular arc, we obtain

k(˜

tN,˜xjN−1

N−1(˜

tN)) −(tN, xjN−1

N−1(tN))k ≤ κδδ+κττ+κββ+κσσ+κ%%.

for sufficiently small δ, τ, β, σ, % > 0 with suitable constants κδ, κτ, κβ, κσ, κ%>0.

Assuming that the initial and final position are not subject to perturbations, as they typically

are given as data, we are allowed to set τ= 0 and it remains to estimate the perturbation of the

final state.

Lemma 13 Let all assumptions in the above construction for a hybrid system of ODEs be satisfied.

Then for sufficiently small δ, β, σ, % > 0, there exist constants κδ, κβ, κσ, κ%>0such that

k˜xjN−1(tN)−xjN−1(tN)k ≤ κδδ+κββ+κσσ+κ%%.

Having obtained a stability estimate for hybrid ODEs, we can easily extend these estimates to

regular hybrid systems of DAEs.

For strangeness-free DAEs in the reduced formulation (8), we use (12c) and assume that

(a) k˙

˜xj

1(t)− L(t, ˜xj

1(t))k ≤ β,

(b) k˜xj

2(t)− R(t, ˜xj

1(t))k ≤ η, for all t∈[˜

tj,˜

tj+1], (19)

for every regular arc, together with all above assumptions for the (decoupled) ordinary differential

equation for xj

1including Lipschitz continuity for all involved functions.

Since the transformation Qin (7) is only used for a theoretical reformulation of the problem

to state the arising perturbations as in (19), it is not subject to computational errors. Since the

computation of the algebraic part, given here as xj

2, is not to subject to a propagation of the error,

we only need to transform the initial value to get the initial value of the differential part, apply

the above results for hybrid systems of ODEs to the differential part, and then consider (19b) at

the end point in order to determine the final algebraic part and final solution by transforming

back according to (7). Since the involved transformations (7) are linear with bounded norm and

since

k˜xjN−1

2,N−1(tN)−xjN−1

2,N−1(tN)k

≤ k˜xjN−1

2,N−1(tN)− R(tN,˜x1,N−1(t))k+kR(tN,˜x1,N−1(tN)) − R(tN, x1,N−1(tN))k

≤η+Rk˜x1,N−1(tN)−x1,N−1(tN)k

with some constant R > 0, Lipschitz continuity is maintained and we obtain the following result.

Theorem 14 Let all assumptions in the above construction for a hybrid system of DAEs be sat-

isfied. Then, for sufficiently small δ, β, η, σ, %, η > 0, there exist constants κδ, κβ, κη, κσ, κ%, κη>0

such that

k˜xjN−1(tN)−xjN−1(tN)k ≤ κδδ+κττ+κββ+κηη+κσσ+κ%%+κηη. (20)

Remark 15 If we assume that there are no computational errors, i. e., if β=η=σ=%=η= 0,

then the estimate (20) reduces to

k˜xjN−1(tN)−xjN−1(tN)k ≤ κδδ,

for sufficiently small δ > 0, cp. Theorem 7.

In this section we have obtained stability estimates for hybrid ODEs and extended the results to

regular strangeness-free hybrid DAEs. There is, however, the immediate question what happens

in the case that a hybrid system of DAEs loses regularity. In this case regularization techniques

are necessary. We discuss this issue in the next section.

5 Non-regular behavior

In this section we discuss hybrid DAEs which have a lack of regularity. In view of Definition 3,

such a lack of regularity can have several reasons which include at least the following cases:

1. more than one switching condition is satisfied;

2. the switching condition has a non-regular solution;

3. there is no finite collection of switching points;

4. the flow does not lead to (Sj

l◦Φj

t,t)(t, zj)>0 for t∈(tj, tj+ε] for all l∈ {1, . . . , Lj}, since

(a) the flow is in at least one of the switching surfaces, i. e., we have that (Sj

l◦Φj

t,t)(t, zj)=0

for t∈(tj, tj+ε] and some l∈ {1, . . . , Lj},

(b) we have instantaneous further switching, i. e., we have infeasibility due to (Sj

l◦Φj

t,t)(t, zj)<

0 for t∈(tj, tj+ε] and some l∈ {1, . . . , Lj}.

Since a numerical treatment of hybrid DAEs with a non-regular solution makes no sense, it is

necessary to perform a regularization by remodeling the given problem. Some of the discussed cases

of non-regularity allow for ad hoc regularization techniques, but observe that these techniques may

be non-physical, since they are not based on the physical background, but they are needed for the

sensible numerical integration and they may also give an indication for a possible reformulation

on physical grounds.

5.1 Regularization techniques

In the following, we separately discuss various regularization techniques, but it should be noted

that these cases may occur simultaneously, so that it may be necessary to apply several of these

together in order to obtain a regularization.

1. Non-unique switching condition. If the end point (tj+1, xj(tj+1)) satisfies more than

one switching condition, then we can simply replace the current mode by selecting one of the

activated switching conditions and omitting the others. The particular choice may be viewed as a

dedicated hierarchy in the switching conditions.

2. Non-regular solution of the switching condition. If the solution of the equation for

the switching condition is non-regular, then small perturbations can lead to a different change of

modes. A well-known example of this case is the modeling of billiard [6, 12], when a moving ball

just touches another ball. For this example, we know that the overall solution depends smoothly on

the initial condition but that the condition number of the Jacobian that determines the switching

point may be very large, especially when there are several hits. In the general case, however, there

seems to be no good ad hoc strategy for regularizing the problem.

3. Infinite number of switching points. It may happen that an infinite number of switching

points occurs in the given integration interval. Since we have assumed this interval to be compact,

an infinite number of switching points implies that there is at least one accumulation point in

the set of switching points. A prominent example for this case is the bouncing ball [6, 12], when

the ball loses a certain percentage of its momentum with every bounce. A simple regularization

strategy in this case is to characterize the behavior of the physical system after the accumulation

point by a mode and to switch to this mode a sufficiently small amount of time before we reach

the accumulation point.

4.(a) Flow in switching surface. In the case that the flow is in the switching surface, we

have that

(Sj

l◦Φj

t,tj)(tj, zj

0) = 0 for all t∈[tj, tj+ε]

holds for one or more l∈ {1, . . . , Lj}. But this means nothing else than that these Sj

lare first

integrals of the DAE of mode j. Thus, we actually have an overdetermined but consistent system

of DAEs of the form

Fj(t, xj,˙xj)=0,

Sj(t, xj)=0,(21)

where Sjgathers all the first integrals among the switching functions. To study the properties

of this system, we go over to the formulation (9) so that (21) takes the form

˙xj

1=Lj(t, xj

1),

2=Rj(t, xj

1),

0 = Sj(t, xj

1, xj

2).

Thus, we have additional constraints of the form

Sj(t, xj

1,Rj(t, xj

1)) = 0,(22)

which may in parts be trivial equations 0 = 0. Of course, these redundant equations can be

simply omitted. Assuming then that d

dxj

Sjhas full row rank, we can solve (22) for some of the

components of x1. For these components, the corresponding differential equations contained in

˙xj

1=Lj(t, xj

1) are redundant and can be omitted. If we denote by Pjthe projector which selects

the components of xj

1that are not fixed by (22), then we get a DAE of the form

Pj˙xj

1=PjLj(t, xj

1),

2=Rj(t, xj

1),

0 = Sj(t, xj

1, xj

2).(23)

The original mode jthen should be replaced by the regular DAE (23), together with the switching

functions which were not first integrals of the original DAE. Of course, the actual construction

should and can be performed on the original (possibly higher-index) formulation of the mode. For

corresponding details, see [19].

4.(b) Instantaneous further switching. The case of instantaneous further switching is

characterized by the observation that

(Sj

l◦Φj

t,tj)(tj, xj(tj)) <0 for all t∈(tj, tj+ε]

holds for one or more l∈ {1, . . . , Lj}. Similar to the discussion in Case 1., we then must select one

of the activated switching conditions, say ˆ

lj, in order to fix a unique new mode for the system. If

we denote the activated switching condition of mode kwhich has led to the current mode jby ˆ

lk,

then we can replace mode kby a new mode consisting of the same DAE and switching functions,

but replacing σk(ˆ

lk) = jby σk(ˆ

lk) = m, where m=σj(ˆ

lj), and the transition function Tj,k

lkby

Tm,k

lk◦Φj

tj,tj◦Tj,k

lk.

The discussed regularization techniques may have to be iterated to deal with systems where

several of the described non-regularities arise. In particular, there may be several instantaneous

further switchings one after the other and one must hope that this terminates after a finite number

of times. It is, however, possible that this iteration actually leads to an infinite loop. One such

case is chattering which is discussed in the next subsection.

5.2 Chattering

Chattering is the phenomenon that we just switch back to the previous mode and the previous

state requires to switch again to the current system. This frequently happens in actual physical

systems, see e. g., [15, 22, 24], and in this case it is necessary to modify the model to incorporate

the chattering effect. On the other hand, even correctly modeled chattering typically leads to

systems with highly oscillatory solutions. For such systems the numerical solution usually is very

difficult and costly. Since one is often not really interested in resolving the chattering but in the

more global behavior of the systems, it is also of interest to develop regularization techniques

especially for this case.

Since it is sufficient for the discussion of chattering to consider only two modes and the specific

switching conditions which are responsible for the chattering, we assume that the hybrid system

is given by the two DAEs

F1(t, x1,˙x1)=0, F2(t, x2,˙x2)=0,

which we may assume to be already strangeness-free, the switching functions S1and S2, and

transition functions

T21 :M1→L2

0, T12 :M2→L1

T12

T21

S1= (S1)−1({0})

S2= (S2)−1({0})

Figure 1: Occurence of chattering

where M1= (S1)−1({0}) and M2= (S2)−1({0}).

Let ˆ

T21 :M1→M2and ˆ

T12 :M2→M1be defined by the following implications

T21(t, x1)=(t, x2,˙x2) =⇒ˆ

T21(t, x1) = (t, x2),

T12(t, x2)=(t, x1,˙x1) =⇒ˆ

T12(t, x2)=(t, x1).

Then the effect of chattering is characterized by the conditions

T21(M1)⊆M2,ˆ

T12(M2)⊆M1

∂

∂t S1+∂

∂x1S1˙x1<0,∂

∂t S2+∂

∂x2S2˙x2<0,(24)

where the sign conditions hold on the switching surfaces in L1

0and L2

0belonging to the switching

conditions that are responsible for the chattering, and by

T21 ◦ˆ

T12 ◦ˆ

T21 =ˆ

T21,ˆ

T12 ◦ˆ

T21 ◦ˆ

T12 =ˆ

T12.(25)

This means that ˆ

T21 is an inner and outer inverse of ˆ

T12 and vice versa, [5, 8].

A very simple method to regularize such a hybrid system is to introduce hysteresis [2] into the

system by moving the switching surface in such a way that the transition functions lead to the

interior of the feasible region, for example by choosing some small δ > 0, replacing the switching

conditions by

S1(t, x1) + δ= 0, S2(t, x2) + δ= 0,(26)

and taking care that the transition functions still map to M1and M2, see the numerical experiments

below.

A second method to deal with chattering is the so-called Filippov regularization, see [11], in the

case that the two modes and the switching surfaces coincide. The idea of Filippov regularization

is to take the unique convex combination of the two flows which yields a flow that is restricted to

the common switching surface. Actually, this corresponds to the introduction of a new mode into

the original hybrid system together with suitable switching conditions and transfer functions.

Being originally developed for hybrid ODE models and generalized to hybrid DAEs in [15, 24, 25],

we generalize this concept further to the case that the variables of the involved models do not

coincide.

For the construction, we first discuss the (autonomous) decoupled case, with d1differential and

a1algebraic equations in the first mode and d2differential and a2algebraic equations in the second

mode,

(d1eq.) ˙x1

1=L1(x1

1),(d2eq.) ˙x2

1=L2(x2

1),

(a1eq.) x1

2=R1(x1

1),(a2eq.) x2

2=R2(x2

1),

cp. Figure 5.2, and we assume without loss of generality that d1≥d2.

Chattering occurs at the time point tif

(T21 ◦ T 12)(x2

1(t)) = x2

1(t) (27)

and the conditions

(S2◦ T 12)(x2

1(t)) = 0,S2(x2

1(t)) = 0,

1(T12(x2

1(t)))L1(T12(x2

1(t))) <0,S2

1(x2

1(t))L2(x2

1(t)) <0(28)

are satisfied, where the lower index as in S1

1means the derivative with respect to x1

A possible regularization would be to restrict the flows to the switching surface. For the second

mode this leads to

(d2eq.) ˙x2

1=L2(x2

1),

(a2eq.) x2

2=R2(x2

1),

(1 eq.) 0 = S2(x2

1),(29)

which is an over-determined and inconsistent system. In particular, it is not directly evident what

would be a meaningful flow on the manifold defined by the algebraic constraints in (29).

Observing that the differentiation of the transition relation

1=T12(x1

yields

˙x2

1=T21

1(x1

1) ˙x1

1=T21

1(T12(x2

1))L1(T12(x2

1)),(30)

the idea of Filippov regularization in the ODE case is to combine the two proposed flows from (29)

and (30) to a possible flow on the constraint manifold of (29).

In this way, we obtain the DAE

1=T21

1(T12(x2

1))L1(T12(x2

1)),

1=L2(x2

1),

(d2eq.) ˙x2

1=αw1

1+ (1 −α)w2

(a2eq.) x2

2=R2(x2

1),

(1 eq.) 0 = S2(x2

(31)

for the unknowns (x2

1, x2

2, α), if we assume that the quantities w1

1, w2

1are eliminated. In particular,

we only need to consider the DAE

(d2eq.) ˙x2

1=αT21

1(T12(x2

1))L1(T12(x2

1)) + (1 −α)L2(x2

1),

(1 eq.) 0 = S2(x2

1)(32)

for (x2

1, α).

Note that (32) is a semi-explicit DAE of the form

˙x=f(x, y),0 = g(x).

It is well-known that such a DAE cannot be strangeness-free, but it satisfies Hypothesis 1 with

µ= 1 if gx(x)fy(x, y) is nonsingular along the solution.

To deal with this situation we assume that regularization with the help of hysteresis as described

above is possible. In particular, we assume that extending the solution a little bit beyond the

switching surface in one mode and then switching to the other mode yields a point which is

feasible with respect to the switching condition in this mode, i. e.,

1(T21(x1

1(t)))T21

1(x1

1(t))L1(x1

1(t)) >0,

1(T12(x2

1(t)))T12

1(x2

1(t))L2(x2

1(t)) >0.(33)

In (32) the quantity gx(x)fy(x, y) is a scalar given by

γ=S2

1(x2

1)[T21

1(T12(x2

1))L1(T12(x2

1)) − L2(x2

1)].

Setting x1

1=T12(x2

1) and using (27), which gives T21(x1

1) = x2

1, we obtain that

γ=S2

1(T21(x1

1))T21

1(x1

1)L1(x1

1)− S2

1(x2

1)L2(x2

and this gives

γ > 0 (34)

due to (28) and (33). Hence, the DAE (31) and so the DAE (32) satisfy Hypothesis 1 with µ= 1.

The corresponding hidden constraint is given by

1(x2

1)[αT21

1(T12(x2

1))L1(T12(x2

1)) + (1 −α)L2(x2

1)] = 0.(35)

Utilizing (34), we can solve (35) for αto obtain

α=−S2

1(x2

1)L2(x2

(x2

1)T21

(T12(x2

1)) − S2

(x2

1)L2(x2

1).

Due to the sign conditions (28) and (33), we see that

α∈(0,1).

Thus we actually use an appropriate convex combination of the involved flows to get a flow on the

constraint manifold of (31).

In the non-autonomous case of the form

(d1eq.) ˙x1

1=L1(t, x1

1),(d2eq.) ˙x2

1=L2(t, x2

1),

(a1eq.) x1

2=R1(t, x1

1),(a2eq.) x2

2=R2(t, x2

1),(36)

we must replace the relation (27) by

(T21 ◦ T 12)(t, x2

1(t)) = (t, x2

1(t))

and, utilizing

T12(t, x2

1(t)) = (t, x1

1(t)),d

dt (T12(t, x2

1(t))) = (1,˙x1

1(t)),

the relations (28) by

(S2◦ T 12)(t, x2

1(t)) = 0,

t(T12(t, x2

1(t))) + S1

1(T12(t, x2

1(t)))L1(T12(t, x2

1(t))) <0,

S2(t, x2

1(t)) = 0,

t(t, x2

1(t)) + S2

1(t, x2

1(t))L2(t, x2

1(t)) <0.

(37)

The relation for x2

1becomes

˙x2

1= Π2

2(T21

t(T12(t, x2

1)) + T21

1(T12(t, x2

1))L1(T12(t, x2

1))),

where Π2

2denotes the projection onto the second argument of (t, x2

1). The regularized DAE then

reads w1

1= Π2

2(T21

t(T12(t, x2

1)) + T21

1(T12(t, x2

1))L1(T12(t, x2

1)))),

1=L2(t, x2

1),

(d2eq.) ˙x2

1=αw1

1+ (1 −α)w2

(a2eq.) x2

2=R2(t, x2

1),

(1 eq.) 0 = S2(t, x2

1).

Eliminating w1

1, w2

1yields the relevant part

(d2eq.) ˙x2

1=αPi2

2(T21

t(T12(t, x2

1)) + T21

1(T12(t, x2

1))L1(T12(t, x2

1))))

+(1 −α)L2(t, x2

1),

(1 eq.) 0 = S2(t, x2

1).

(38)

The assumptions (33) take the form

t(T21(t, x1

1(t))) + S2

t,x2

1(T21(t, x1

1(t)))T21

1(t, x1

1(t))L1(t, x1

1(t)) >0,

t(T12(t, x2

1(t)))S1

1(T12(t, x2

1(t)))T12

1(t, x2

1(t))L2(t, x2

1(t)) >0.(39)

The hidden constraint in (38) is given by

t(t, x2

1) + S2

1(t, x2

·[αΠ2

2(T21

t(T12(t, x2

1)) + T21

1(T12(t, x2

1))L1(T12(t, x2

1))) + (1 −α)L2(t, x2

1)]

=α[S2

t(t, x2

1) + S2

1(t, x2

1)Π2

2(T21

t(T12(t, x2

1)) + T21

1(T12(t, x2

1))L1(T12(t, x2

1)))]

(1 −α)[S2

t(t, x2

1) + S2

1(t, x2

1)L2(t, x2

1)].

The term in the first bracket is positive due to (39), and the term in the second bracket is negative

due to (37), such that we are in the same situation as in the autonomous case leading to a well-

defined α∈(0,1).

In order to lift (36) to the general formulation

(d1eq.) ˆ

1(t, x1,˙x1) = 0,(d2eq.) ˆ

1(t, x2,˙x2)=0,

(a1eq.) ˆ

2(t, x1) = 0,(a2eq.) ˆ

2(t, x2)=0,

we consider the combined system

(d1eq.) ˆ

1(t, x1,˙w1) = 0,(d2eq.) ˆ

1(t, x2,˙w2)=0,

[ˆ

2(t, x1)=0,] [ ˆ

2(t, x2)=0,]

(a1eq.) ˆ

2;t(t, x1) + ˆ

2;x1(t, x1)w1= 0,(a2eq.) ˆ

2;t(t, x2) + ˆ

2;x2(t, x2)w2= 0,

[S1(t, x1)=0,] [ S2(t, x2)=0,]

(n1eq.) (t, x1) = T12(t, x3),(n2eq.) (t, x2) = (t, x3),

(n3eq.) ˙x3=αΠ2

2(T21

t(T12(t, x3)) + T21

x1(T12(t, x3))w1) + (1 −α)w2,

(a2eq.) ˆ

2(t, x3)=0,

(1 eq.) S2(t, x3)=0,

(40)

where the bracketed relations are redundant and can therefore be omitted. Since the matrices

"ˆ

1,˙x1

2,x1#,"ˆ

1,˙x2

2,x2#

are nonsingular due to Hypothesis 1, the first part of (40) fixes w1∈Rn1and w2∈Rn2. The

second part of (40) fixes x1∈Rn1and x2∈Rn2in terms of x3∈Rn2. The last three relations of

(40) finally form an over-determined DAE for x3∈Rn2.

If we can show that the flow ˙x3is consistent with ˆ

2(t, x3) = 0, then there are actually only

d2relevant relations within the n2equations defining ˙x3, which would lead to a square system for

the unknowns (x3, α).

Starting from

2,t(t, x3) + ˆ

2,x2˙x3(t, x3)

=ˆ

2,t(t, x3) + ˆ

2,x2(t, x3)

·[αΠ2

2(T21

t(T12(t, x3)) + T21

x1(T12(t, x3))w1) + (1 −α)w2]

=α[ˆ

2,t(t, x3) + ˆ

2,x2(t, x3)Π2

2(T21

t(T12(t, x3)) + T21

x1(T12(t, x3))w1)]

+ (1 −α)[ ˆ

2,t(t, x3) + ˆ

2,x2(t, x3)w2],

(41)

we see that the second bracket in the last term of (41) vanishes due to (40). Furthermore, since

2(T21(t, x1)) ≡0

yields

F2;t(T21(t, x1)) + ˆ

F2;x2(T21(t, x1))T21

t(t, x1)≡0,

F2;x2(T21(t, x1))Π2

2T21

x1(t, x1)≡0,

the first bracket in the last term of (41) vanishes as well. Hence, the desired consistency of the

above DAE for x3is shown.

In this section we have discussed regularization techniques for various cases of non-regularity.

In particular, we have generalized the well-known Filippov regularization to the DAE case.

6 Examples and numerical experiments

To illustrate the analysis in the last sections, in this section we present some examples and numer-

ical experiments. We use a C++ implementation of a DAE class for arbitrary index, see [19], using

automatic differentiation together with a root finder, extended by an interface for mode switching

along the lines described in the previous sections.

Example 16 A two-dimensional pendulum can be modeled via the Hamilton formalism starting

with the Hamilton function

H(p, q) = 1

2m−1p2−mgL cos q,

with mthe mass, Lthe length of the pendulum, gthe gravity constant, qthe generalized position

(here the angle between the negative vertical axis and the pendulum), and pthe corresponding

generalized momentum. In particular, Hdescribes the total energy consisting of kinetic and

potential energy. By the Hamilton formalism, the equations of motion are given by

˙p=−∇qH(p, q) = −mgL sin q,

˙q=∇pH(p, q) = m−1p.

It is well-known that such a system is conservative, i. e., that the total energy is preserved, due to

dt H(p, q) = ∇pH(p, q)T˙p+∇qH(p, q)T˙q

=−∇pH(p, q)T∇qH(p, q) + ∇qH(p, q)T∇pH(p, q)=0.

Let H0denote the total energy fixed by the initial values for p, q. Then, the switching function

S(t, p, q) = H(p, q)−H0

has the property that the flow connected with the equations of motions lies in the switching surface.

Consequently, we are in Case 4.(a) of Section 5. In this example, the only algebraic constraint in

the arising overdetermined problem

˙p=−mgL sin q, ˙q=m−1p, H(p, q)−H0= 0

has the Jacobian

J=m−1p mgL sin q,

corresponding to ZT

2Fµ;xof Hypothesis 1. Hence, the variables not fixed by the algebraic con-

straints are described by

T2=−mgL sin q

m−1p,

again in the notation of Hypothesis 1. Taking T2to project onto the relevant part of the original

equations of motion, we end up with the regularized system

−mgL sin q( ˙p+mgL sin q) + m−1p( ˙q−m−1p)=0,

2m−1p2−mgL cos q−H0= 0.

-0.02

-0.015

-0.01

-0.005

0.005

0.01

0.015

0.02

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

Figure 2: Regularized chattering by generalized Filippov (black) and hysteresis (red)

Example 17 Chattering may occur in the hybrid system

˙x1

1=−1,˙x1

2= 0,˙x2

1= 1 −t2

1(t, x1

1, x1

2) = x1

1, S2

1(t, x2

1) = −x2

T21

1(t, x1

1, x1

2)=(t, x1

1,1−t2), T12

1(t, x2

1)=(t, x2

1,0,−1,0)

for t∈(−1,1). The generalized Filippov regularization as described in the previous section is

given by the new mode

˙x1

2= 0,

1= 0,

1=x2

1(t, x1

1, x1

2, x2

1)=1,

2(t, x1

1, x1

2, x2

1)=1−t2,

T13

1(t, x1

1, x1

2, x2

1) = (t, x1

1, x1

2,−1,0),

T23

1(t, x1

1, x1

2, x2

1) = (t, x2

1,1−t2),

together with new switching and transfer functions

1(t, x1

1, x1

2) = x1

1, S2

1(t, x2

1) = −x2

T31

1(t, x1

1, x1

2)=(t, x1

1, x1

2, x1

1,0,0,0), T32

1(t, x2

1) = (t, x2

1,0, x2

1,0,0,0)

for the original modes for t∈(−1,1). Figure 6 shows the solution profile of the first component of

all modes for the initial value (t0, z1

0) = (−1,1,1,−1,0) in Mode 1 using the generalized Filippov

regularization and hysteresis as described in (26) for δ= 0.01.

Example 18 An evaporator vessel can be modeled by the hybrid system

I˙p1=−Rpp1, I ˙p2=−Rpp2+L2,

C˙

L1=−L1/Rb+fin, C ˙

L2=−p2/I −L2/Rb+fin

S1(t, p1, L1) = L1−Lth, S2(t, p2, L2) = −L2+Lth,

T21(t, p1, L1)=(t, p1, L1,˙p2,˙

L2), T12(t, p2, L2)=(t, p2, L2,˙p1,˙

L1),

0.078

0.079

0.08

0.081

0.082

0.083

0.084

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14

Figure 3: Regularized chattering by generalized Filippov (black) and hysteresis (red)

where the derivatives in the transition functions denote the values obtained by the differential

equations of the new mode, see [24] and references therein. The two modes of this example are

separated by the line L=Lth. With the choice

Rb= 1, Rp= 0.5, I = 0.5, C = 15, fin = 0.25, Lth = 0.08

of parameters we have that ˙

L > 0 in the first mode while in the second mode the line L=Rb(fin −

p/I) separates the regions with ˙

L > 0 for smaller pand ˙

L < 0 for larger p. Hence chattering occurs

when the solution hits L=Lth in the region to the right of the line L=Rb(fin −p/I), cp. Figure 6.

Here, Filippov regularization yields the new mode

I˙p3=−αRpp3−(1 −α)(Rpp3−L3),

C˙

L3=−αL1/Rb−(1 −α)(p3/I +L3/Rb) + fin,

L3−Lth = 0,

1(t, p3, L3) = p3/I +L3/Rb−fin,

2(t, p3, L3) = −L3/Rb+fin,

T31(t, p3, L3)=(t, p3, L3,˙p3,˙

L3),

T32(t, p3, L3)=(t, p3, L3,˙p3,˙

L3).

Figure 6 shows the solution in the phase plane (p, L) for the initial value (t0, p0, L0) = (0.0,0.05,0.06)

in Mode 1 using the generalized Filippov regularization and hysteresis as described in (26) for

δ= 0.0004.

7 Conclusions

We have discussed the analysis and numerical solution of hybrid systems of DAEs. It has been

shown that it is possible to characterize regular arcs of a solution and how to estimate the stability

of a global solution. For systems which do not satisfy the regularity assumptions, we have described

different regularization techniques, in particular the generalization of Filippov regularization to

hybrid systems of DAEs where we even allow that the modes are formulated in different variables.

We have illustrated the results via several numerical examples. It is an open problem whether the

described cases of non-regularity are a complete set.

References

[1] J. Agrawal, K.M. Moudgalya, and A.K. Pani. Sliding motion of discontinuous dynamical

systems described by semi-implicit index one differential algebraic equations. Chemical Engin.

Science, 61:4722–4731, 2006.

[2] J.C. Alexander and T.I. Seidman. Sliding modes in intersecting switching surfaces II: hys-

teresis. Houston J. Math., 25:185–211, 1999.

[3] Modelica Association. Modelica language specification, version 3.0. Modelica Association.

[4] P.I. Barton and C.K. Lee. Modeling, simulation, sensitivity analysis, and optimization of

hybrid systems. ACM Trans. Modeling Computer Simulation, 12:256–289, 2002.

[5] A. Ben-Israel and T. N. E. Greville. Generalized Inverses: Theory and Applications. John

Wiley and Sons, New York, NY, 1973.

[6] B. Brogliato. Nonsmooth Mechanics. Springer-Verlag, New York, 1999.

[7] S. L. Campbell. A general form for solvable linear time varying singular systems of differential

equations. SIAM J. Math. Anal., 18:1101–1115, 1987.

[8] S. L. Campbell and C. D. Meyer. Generalized Inverses of Linear Transformations. Pitman,

San Francisco, CA, 1979.

[9] Dynasim AB, Ideon Research Park - SE-223 70 Lund - Sweden. Dymola,Multi-Engineering

Modelling and Simulation, 2006.

[10] E. Eich-Soellner and C. F¨uhrer. Numerical Methods in Multibody Systems. Teubner Verlag,

Stuttgart, Germany, 1998.

[11] A.F. Filippov. Differential Equations with Discontinuous Right-hand Sides. Kluwer, Dor-

drecht, 1998.

[12] R. Goebel, R.G. Sanfelice, and A. Teel. Hybrid dynamical systems. Control Systems, IEEE,

29:28–93, 2009.

[13] E. Hairer, S. P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I: Nonstiff

Problems. Springer-Verlag, Berlin, Germany, 1st edition, 1987.

[14] P. Hamann. Modellierung und Simulation von realen Planetengetrieben. Diplomarbeit, In-

stitut f¨ur Mathematik, TU Berlin, Berlin, Germany, 2003.

[15] P. Hamann and V. Mehrmann. Numerical solution of hybrid differential-algebraic equations.

Comp. Meth. Appl. Mech. Eng., 197:693–705, 2008.

[16] R. Isermann. Mechatronic systems: concepts and applications. Trans. Institute of Measure-

ment and Control, 22:29–55, 2000.

[17] D. Jeon and M. Tomizuka. Learning hybrid force and position control of robot manipulators.

IEEE Trans. Robotic. Automat., 9:423–431, 1996.

[18] T. Krilavicius. Hybrid techniques for hybrid systems. PhD thesis, University of Twente,

Enschede, 2006.

[19] P. Kunkel and V. Mehrmann. Differential-Algebraic Equations. Analysis and Numerical So-

lution. EMS Publishing House, Z¨urich, Switzerland, 2006.

[20] C. Lacoursi`ere. Ghosts and Machines: Regularized Variational Methods for Interactive Sim-

ulation of Multibodies with Dry Frictional Contacts. Phd thesis, Ume˚a University, Sweden,

2007.

[21] S. Lang. Analysis I. Addison-Wesley, Reading, MA, 3rd edition, 1973.

[22] J. Lunze and F. Lamnabhi-Laguarrigue. Handbook of Hybrid Systems Control. Theory, Tools

and Applications. Cambridge University Press, Cambridge, UK, 2009.

[23] J. Lygeros, G.J. Pappas, and S. Sastry. An Introduction to Hybrid System Modeling, Analysis,

and Control. Preprints of the First Nonlinear Control Network Pedagogical School, Athens,

Greece, pages 307–329, 1999.

[24] V. Mehrmann and L. Wunderlich. Hybrid systems of differential-algebraic equations – analysis

and numerical solution. J. Process Control, 19:1218–1228, 2009.

[25] L. Wunderlich. Analysis and Numerical Solution of Structured and Switched Differential-

Algebraic Systems. Phd thesis, TU Berlin, 2008.

[26] M. Yu, L. Wang, T. Chu, and G. Xie. Stabilization of networked control systems with data

packet dropout and network delays via switching system approach. In IEEE Conference on

Decision and Control, pages 3539–3544, 2004.