Analyse et r esolution num erique d equations alg ebro-diff - - PowerPoint PPT Presentation

Analyse et r´ esolution num´ erique d’´ equations alg´ ebro-diff´ erentielles

Alexandre Chapoutot

ENSTA ParisTech

2018-2019

Discontinuous simulation

1

Discontinuous simulation

2

Differential Algebraic Equations

2 / 40

IVP

Recall our starting point is the IVP of ODE defined by ˙ y = f (t, y) with y(0) = y0 , (1) for which we want the solution y(t; y0) given by numerical integration methods i.e. a sequence of pairs (ti, yi) such that yi ≈ y(ti; y0) .

3 / 40

Why do we consider discontinuities?

Need to model non-smooth behaviors, e.g., solid body in contact with each other interaction between computer and physics, e.g., control-command systems constraints on the system, e.g., robotic arm with limited space

4 / 40

Simulation with discontinuous systems

There are two kinds of events: time event: only depending on time as sampling state event: depending on a particular value of the solution of ODE

• r DAE.

To handle these events we need to adapt the simulation algorithm. Time events are known before the simulation starting. Hence we can use the step-size control to handle this. State event should be detect and handle on the fly. New algorithms are needed.

5 / 40

IVP with discontinuities

An IVP for ODE with discontinuities is defined by ˙ y =

• f1(t, y)

if g(t, y) 0 f2(t, y)

• therwise

with y(0) = y0 , (2) for which we want the solution y(t; y0) given by numerical integration methods i.e. a sequence of pairs (ti, yi) such that yi ≈ y(ti; y0) .

6 / 40

Example: zero-crossing detection

A simple example ˙ y =

• f1(t, y)

if g(y) 0 f2(t, y)

• therwise

. Legend

7 / 40

Zero-crossing event detection

Main steps

Detection of zero-crossing event Is one of the zero-crossing changed its sign between [tn, tn + hn]? Localization: if detection is true Bracket the most recent zero-crossing time using bisection method. Pass through the zero-crossing event in two steps:

Set the next major output to the left bound of the bracket time. Reset the solver with the state estimate at the right bound of bracket time.

Ingredients for the localization – 1

Continuous extension (method dependent) to easily estimate state. For example, ode23 uses Hermite interpolation

p(t) = (2τ 3 − 3τ 2 + 1)yn + (τ 3 − 2τ 2 + τ)(t2 − t1)f (yn) + (−2τ 3 + 3τ 2)yn+1 + (τ 3 − τ 2)(t2 − t1)f (yn+1)

with τ = t−tn

hn

8 / 40

Zero-crossing event detection

Main steps

Detection of zero-crossing event Is one of the zero-crossing changed its sign between [tn, tn + hn]? Localization: if detection is true Bracket the most recent zero-crossing time using bisection method. Pass through the zero-crossing event in two steps:

Set the next major output to the left bound of the bracket time. Reset the solver with the state estimate at the right bound of bracket time.

Ingredients for the localization – 2

The solve the equation g(t, p(t)) = 0 instead of g(t, y(t)) = 0 Note: as this equation is 1D then algorithm as bisection or Brent’s method can be used instead of Newton’s iteration.

8 / 40

Simulation algorithm

Data: f1 the dynamic, f2 the dynamic, g the zero-crossing function, y0 initial condition, t0 starting time, tend end time, h integration step-size, tol t ← t0; y ← y0; f ← f1; while t < tend do Print(t, y); y1 ← Euler(f ,t,y,h); y2 ← Heun(f ,t,y,h); if ComputeError(y1, y2) is smaller than tol then if g(y) · g(y1) < 0 then Compute p(t) from y, f (y), y1 and f (y1); [t−, t+] = FindZero (g(p(t))); Print (t + t−, p(t−)); f ← f2; y ← p(t+); t ← t + t+; end y ← y1; t ← t + h; h ← ComputeNewH (h, y1, y2); end h ← h/2 end

9 / 40

Discontinuities and numerical integration methods

One-step methods are more robust than multi-step in case of discontinuities

10 / 40

Differential Algebraic Equations

1

Discontinuous simulation

2

Differential Algebraic Equations

11 / 40

Definition of Differential Algebraic Equations (DAE)

We consider a differential system of equation F(˙ x, x, t) =      F1(˙ x(t), x(t), t) F2(˙ x(t), x(t), t) . . . Fn(˙ x(t), x(t), t)      = 0 with ˙ x(t), x(t) ∈ Rn. This system is a DAE if the Jacobian matrix ∂F ∂ ˙ x is singular

12 / 40

Example of DAE

The following system is a DAE x1 − ˙ x1 + 1 = 0 ˙ x1x2 + 2 = 0 ⇒ F(˙ x, x, t) =

• x1 − ˙

x1 + 1 ˙ x1x2 + 2

• with

x =

• x1

x2

• The Jacobian of F w.r.t. ˙

x is ∂F ∂ ˙ x = ∂F1

∂ ˙ x1 ∂F1 ∂ ˙ x2 ∂F2 ∂ ˙ x1 ∂F2 ∂ ˙ x2

• =
• −1

x2

• ⇒

det ∂F ∂ ˙ x

• = 0

Note in this example ˙ x2 is not explicitly defined.

13 / 40

Example of DAE continued

Solving DAE is a hard challenge either symbolically or numerically. Special DAE forms are usually considered: linear, Hessenberg form, etc.

Example, we rewrite the previous system

solving for ˙ x1 the equation x1 − ˙ x1 + 1 = 0 ⇒ ˙ x1 = x1 + 1 Substitute ˙ x1 in ˙ x1x2 + 2 = 0 we get ˙ x1 = x1 + 1

Ordinary differential equation

(x1 + 1) x2 + 2 = 0

Algebraic equation

Note: this form of DAE is used in many engineering applications. mechanical engineering, process engineering, electrical engineering, etc. Usually: dynamics of the process + laws of conservation

14 / 40

Engineering examples of DAE - Chemical reaction

An isothermal continuous flow stirred-tank reactor1 (CSTR) with elementary reactions: A ⇋ B → C assuming reactant A with a in-flow rate Fa and concentration CA0 Reversible reaction A ⇋ B is much faster that B → C, i.e., k1 ≫ k2

˙ V = Fa − F ˙ CA = Fa V (CA0 − CA) − R1 ˙ CB = −Fa V CB + R1 − R2 ˙ CC = −Fa C C + R2 0 = CA − CB Keq 0 = R2 − k2CB

1Control of Nonlinear DAE Systems with Applications to Chemical Processes

15 / 40

Engineering examples of DAE - Chemical reaction

R1 and R2 rates of reactions F output flow CA, CB and CC are concentrations of A, B and C. Let x = (V , CA, CB, CC) z = (R1, R2) we get ˙ x = f (x, z) 0 = g(x, z)

16 / 40

Engineering examples of DAE - Mechanical system

Pendulum

second Newton’s law m¨ x = −F ℓ x m¨ y = mg F ℓ y Mechanical energy conservation x2 + y 2 = ℓ2

˙ x1 = x3 ˙ x2 = x4 ˙ x3 = −F ℓ x1 ˙ x4 = g F ℓ x2 0 = x 2

1 + x 2 2 − ℓ2

17 / 40

Engineering examples of DAE - Electrical system

Ohm’s law

C ˙ VC = iC, L ˙ V = iL, VR = RiR

Kirchoff’s voltage and current laws Conservation of current

iE = iR, iR = iC, iC = iL

Conservation of energy

VR + VL + VC + VE = 0

And we get

˙ VC = 1 C iL ˙ VL = 1 LiL 0 = VR + RiE 0 = VE + VR + VC + VL 0 = iL − iE

18 / 40

Engineering examples of DAE - Electrical system

Let x = (VC, VL, VR, iL, iE) we have ˙ x =      

1 C 1 L

      x 0 =   1 R 1 1 1 1 −1   x +   1   VE which is of the form ˙ x = Ax 0 = Bx + Dz

19 / 40

Method of Lines for PDE

Consider the linear PDE (diffusion equations) ∂u ∂t (x, t) = D ∂2u ∂x2 (x, t) with    u(x = x0, t) = ub ∂u ∂x (x = xf , t) = 0 and D a constant. Using method of lines, we have with an equally spaced grid for x (finite difference) ∂2u ∂x2 ≈ ui+1 − 2ui + ui−1 ∆x2 + O

• ∆x2

Hence, we get dui dt = D ui+1 − 2ui + ui−1 ∆x2 for i = 1, 2, . . . , M

20 / 40

Method of Lines for PDE

In other words, we get the system u1 = ub du2 dt = D u3 − 2u2 + ub ∆x2 du3 dt = D u4 − 2u3 + u2 ∆x2 . . . duM dt = D uM+1 − 2uM + uM−1 ∆x2 uM+1 = uM Note uM+1 is outside of the grid so we add an extra constraints. Hence we get a DAE

21 / 40

Method of Lines for PDE

Figure 1.1. MOL solution of Eq. (1.1) illustrating the origin of the method of lines

22 / 40

Classification of DAE

Nonlinear DAE if it is of the form F(˙ x, x, t) = 0 and it is nonlinear w.r.t. any one of ˙ x, x, or t Linear DAE if it is of the form A(t)˙ x + B(t)x = c(t) If A(t) ≡ A and B(t) ≡ B then the DAE is time-invariant Semi-explicit DAE it is of the form ˙ x = f (t, x, z) 0 = g(t, x, z) z is the algebraic variable and x is a differential/state variable Fully implicit DAE it is of the form F(˙ x, x, t) = 0

23 / 40

Classification of DAE - cont

Note any DAE can be written in a semi-explicit form.

Conversion of fully implicit form

F(˙ x, x, t) = 0

˙ x=z

⇔

• ˙

x = z 0 = F(z, x, t) Remark this transformation does not make the solution more easier to get But useful in case of linear DAE, see next.

24 / 40

Classification of DAE - cont

Consider a linear time-invariant DAE A˙ x + Bx + b(t) = 0 assuming that λA + B (matrix pencil) is not singular for some scalar λ. Then it exists non-singular matrices G and H of size n × n such that: GAH = Im N

• and

GBH = J In−m

• Im is the identity matrix of size m × m (m ≤ n)

In−m is the identity matrix of size (n − m) × (n − m) N is a nilpotent matrix, i.e., ∃p ∈ N+, Np = 0 J ∈ Rm×m

25 / 40

Classification of DAE - cont

Hence A˙ x + Bx + b(t) = 0 ⇔ (GAH)(H−1)˙ x + (GBH)(H−1)x + Gb(t) = 0 ⇔ Im N

• H−1 ˙

x + J In−m

• H−1x + Gb(t) = 0

⇔ with w(t) = H−1x Im N

• ˙

w + J In−m

• w + Gb(t) = 0

Let w = (w1, w2)T with w1 ∈ Rm and w2 ∈ Rn−m, b = (b1, b2)T we get ˙ w1 + Jw1 + b1(t) = 0 Nw1 + w2 + b2(t) = 0

From Nilpotency property, we get

˙ w1 = −Jw1 − b1(t) 0 = −(Np)−1w2 − (Np)−1b2(t)

26 / 40

Index of DAE

Remark

There are several definitions of an index. Each measure a different aspect of the DAE. Differential index (δ) measure the degree of singularity. Perturbation index (π) measure the influence of numerical approximation. etc.

Definition of differential index

The index of a DAE system F (˙ x, x, t) = 0 is the minimum number of times certain equations in the DAE must be differentiated w.r.t. t, in

• rder to transform the problem into an ODE.

Remark: (differential) index can be seen as a measure of the distance between the DAE and the corresponding ODE. Remark: mathematical properties are lost with differentiation!

27 / 40

DAE and index

Definition of index

The differential index k of a sufficiently smooth DAE is the smallest k such that: F(˙ x, x, t) = 0 ∂F ∂t (˙ x, x, t) = 0 . . . ∂kF ∂tk (˙ x, x, t) = 0 uniquely determines ˙ x as a continuous function of (x, t).

28 / 40

Differential index and DAE – example

Let ˙ x1 = x1 + 1 (x1 + 1) x2 + 2 = 0 with x2 the algebraic variable. Differentiation of g w.r.t. t, d dt g(x1, x2) = 0 ⇒ ˙ x1x2 + (x1 + 1)˙ x2 = 0 ⇒ ˙ x2 = − ˙ x1x2 x11 = −x2 Only one differentiation is needed to define ˙ x2, this DAE is index 1 Other examples, CSTR is index 2 Pendulum is index 3 There are higher index DAEs (index > 1) Index reduction is used to go from higher index to lower index DAE (cf Khalil Ghorbal’s lecture)

29 / 40

DAE family and differential index

Index 0

ODE system ˙ x = f (t, x(t))

Index 1

Algebraic equation y = q(t)

Index 1

DAE in Hessenberg form of index 1 ˙ x = f (t, x, y) 0 = g (x, y) with ∂g ∂y is non-singular

30 / 40

Examples of differential index - cont.

Index 2

DAE in Hessenberg form of index 2 ˙ x = f (t, x, y) 0 = g (t, x) with ∂g ∂x ∂f ∂y is non-singular

Index 3

DAE in Hessenberg form of index 3 ˙ x = f (t, x, y, z) ˙ y = g (t, x, y) 0 = h(t, y) with ∂h ∂y ∂g ∂x ∂f ∂z is non-singular e.g., mechanical systems

31 / 40

Perturbation index

The DAE has the perturbation index k along a solution x if k is the smallest integer such that, for all functions x(t) having the defect f (˙ xδ, xδ, t) = δ(t) there exists an estimate x(t)−xδ(t) ≤ C

• x(t0) − xδ(t0) + max

t

δ(t) +maxt δ′(t) + · · · + max

t

δ(k−1)(t)

• for a constant C > 0, if δ is small enough.

Property:

δ ≤ π ≤ δ + 1

32 / 40

Issues of index reduction

Issues

Consistent initial conditions finding initial value for differential and algebraic variables may be very difficult. For F(˙ x, x, t) = 0 x0 is a consistent initial value, if there exists a smooth solution that fulfills x(0) = x0 and this solution is defined for all t. E.g., semi-explicit DAE with only x(0) = x0 what about the algebraic variable? Drift off effect when applying index reduction the solution of the lower index DAE may not be of the original index. In consequence, tools/methods to solve DAE should provide automatic index reduction be able to find consistent initial values e.g., Dymola/Modelica

33 / 40

Example of consistent initial value

Let ˙ u = −0.5(u + v) + q1(t) 0 = 0.5(u − v) − q2(t) If u(0) is given we can determine v(0) = u(0) − 2q2(0) and so ˙ u(0). Set u = y1 + y2 and v = y1 − y2 we get ˙ y1 + ˙ y2 = −y1 + q1(t) 0 = y2 − q2(t) For consistency we must have y2(0) = q2(0) but we can choose y1(0) arbitrarily but we cannot determine ˙ y1(0) without using ˙ y2(0) = ˙ q2(0).

34 / 40

Example of drift off effect

Going from index 3 pendulum to index 2 by differentiating the constraint x2

1 + x2 2 − ℓ2 = 0 leads to

˙ x1 = x3 (3) ˙ x2 = x4 (4) ˙ x3 = −F ℓ x1 (5) ˙ x4 = g F ℓ x2 (6) 0 = x1x3 + x2x4 (7)

m l

3.2. Expanding set of solutions due to index

Comments: solid line curve is the result of index 3 pendulum problem Constraint (7) says the velocity should orthogonal to the position. Index reduction increase the space of solution with dashed line curves

35 / 40

A small theory of DAE

For ODE, we have a theorem applying on a large class of problem proving the existence and unicity of the solution No such theorem exists for DAE Instead we have some theorems of solvability of different kinds of DAE Linear constant coefficient DAE Linear time varying coefficient DAE Non-linear DAE

36 / 40

Solvability of DAE

Definition

Let I be an open sub-interval of R, Ω a connected open subset of R2m+1, and F a differentiable function from Ω to Rm. Then the DAE F(˙ x, x, t) = 0 is solvable on I in Ω if there is an r-dimensional family of solutions φ(t, c) defined on a connected open set I × ˜ Ω, ˜ Ω ⊂ Rr, such that

1

φ(t, c) is defined on all of I for each c ∈ ˜ Ω

2

(φ′(t, c), φ(t, c), t) ∈ Ω for (t, c) ∈ I × ˜ Ω

3

If ψ(t) is any other solution with (ψ′(t, c), ψ(t, c), t) ∈ Ω then ψ(t) = φ(t, c) for some c ∈ ˜ Ω

4

The graph of φ as a function of (t, c) is an r + 1-dimensional manifold.

37 / 40

Solvability of linear constant constant DAE

Let A˙ x + Bx = f And consider the matrix pencil λA + B A matrix pencil is regular if det(λA + B) is not identically zero as a function of λ.

Theorem

The linear constant coefficient DAE is solvable if and only if λA + B is regular pencil. Note: the degree of nilpotency of the matrix N used in the decomposition is also the index number of the DAE.

38 / 40

Relation between DAE and stiff ODE

Singularly perturbed ODE systems are of the form ˙ x = f (t, x, y) (8) ε˙ y = g(t, x, y) (9) When ε = 0 then we get a DAE but Eq. (8) is usually stiff. DAE can be seen as infinitely stiff. Consequence not all numerical method to solve ODE can be used to solve DAE! Usually, we consider stiffly stable methods

39 / 40

Conclusion

DAE are a generalisation of ODE but there is no general theorem to prove existence of the solution of DAE differentiation used to index reduction can introduce singularities the class of numerical methods used to solve DAE is rather small compare to ODE.

40 / 40