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) = y 0 , (1) for which we want the solution y ( t ; y 0 ) given by numerical integration methods i.e. a sequence of pairs ( t i , y i ) such that y i ≈ y ( t i ; y 0 ) . 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 or 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 � f 1 ( t , y ) if g ( t , y ) � 0 y = ˙ with y (0) = y 0 , (2) f 2 ( t , y ) otherwise for which we want the solution y ( t ; y 0 ) given by numerical integration methods i.e. a sequence of pairs ( t i , y i ) such that y i ≈ y ( t i ; y 0 ) . 6 / 40
Example: zero-crossing detection A simple example � f 1 ( t , y ) if g ( y ) � 0 y = ˙ . f 2 ( t , y ) otherwise Legend 7 / 40
Zero-crossing event detection Main steps Detection of zero-crossing event Is one of the zero-crossing changed its sign between [ t n , t n + h n ]? 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) y n + ( τ 3 − 2 τ 2 + τ )( t 2 − t 1 ) f ( y n ) + ( − 2 τ 3 + 3 τ 2 ) y n +1 + ( τ 3 − τ 2 )( t 2 − t 1 ) f ( y n +1 ) with τ = t − t n h n 8 / 40
Zero-crossing event detection Main steps Detection of zero-crossing event Is one of the zero-crossing changed its sign between [ t n , t n + h n ]? 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: f 1 the dynamic, f 2 the dynamic, g the zero-crossing function, y 0 initial condition, t 0 starting time, t end end time, h integration step-size, tol t ← t 0 ; y ← y 0 ; f ← f 1 ; while t < t end do Print( t , y ); y 1 ← Euler( f , t , y , h ); y 2 ← Heun( f , t , y , h ); if ComputeError( y 1 , y 2 ) is smaller than tol then if g ( y ) · g ( y 1 ) < 0 then Compute p ( t ) from y , f ( y ), y 1 and f ( y 1 ); [ t − , t + ] = FindZero ( g ( p ( t ))); Print ( t + t − , p ( t − )); f ← f 2 ; y ← p ( t + ); t ← t + t + ; end y ← y 1 ; t ← t + h ; h ← ComputeNewH ( h , y 1 , y 2 ); 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 1 (˙ x ( t ) , x ( t ) , t ) F 2 (˙ x ( t ) , x ( t ) , t ) F (˙ x , x , t ) = = 0 . . . F n (˙ x ( t ) , x ( t ) , t ) x ( t ) , x ( t ) ∈ R n . with ˙ This system is a DAE if the Jacobian matrix ∂ F is singular ∂ ˙ x 12 / 40
Example of DAE The following system is a DAE x 1 − ˙ x 1 + 1 = 0 � � � � x 1 − ˙ x 1 + 1 x 1 ⇒ F (˙ x , x , t ) = with x = x 1 x 2 + 2 ˙ x 2 x 1 x 2 + 2 = 0 ˙ The Jacobian of F w.r.t. ˙ x is � ∂ F 1 ∂ F 1 � � � � ∂ F � ∂ F − 1 0 ∂ ˙ ∂ ˙ x 1 x 2 x = = ⇒ det = 0 0 ∂ ˙ ∂ F 2 ∂ F 2 x 2 ∂ ˙ x ∂ ˙ x 1 ∂ ˙ x 2 Note in this example ˙ x 2 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 ˙ x 1 the equation x 1 − ˙ x 1 + 1 = 0 ⇒ ˙ x 1 = x 1 + 1 Substitute ˙ x 1 in ˙ x 1 x 2 + 2 = 0 we get x 1 = x 1 + 1 ˙ Ordinary differential equation ( x 1 + 1) x 2 + 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 reactor 1 (CSTR) with elementary reactions: A ⇋ B → C assuming reactant A with a in-flow rate F a and concentration C A 0 Reversible reaction A ⇋ B is much faster that B → C , i.e., k 1 ≫ k 2 ˙ V = F a − F C A = F a ˙ V ( C A 0 − C A ) − R 1 C B = − F a ˙ V C B + R 1 − R 2 C C = − F a ˙ C C + R 2 0 = C A − C B K eq 0 = R 2 − k 2 C B 1 Control of Nonlinear DAE Systems with Applications to Chemical Processes 15 / 40
Engineering examples of DAE - Chemical reaction R 1 and R 2 rates of reactions F output flow C A , C B and C C are concentrations of A , B and C . Let x = ( V , C A , C B , C C ) z = ( R 1 , R 2 ) we get x = f ( x , z ) ˙ 0 = g ( x , z ) 16 / 40
Engineering examples of DAE - Mechanical system Pendulum second Newton’s law Mechanical energy conservation x 2 + y 2 = ℓ 2 x = − F m ¨ ℓ x y = mg F m ¨ ℓ y x 1 = x 3 ˙ x 2 = x 4 ˙ x 3 = − F ˙ ℓ x 1 x 4 = g F ˙ ℓ x 2 0 = x 2 1 + x 2 2 − ℓ 2 17 / 40
Engineering examples of DAE - Electrical system Ohm’s law C ˙ L ˙ V C = i C , V = i L , V R = Ri R Kirchoff’s voltage and current laws Conservation of current i E = i R , i R = i C , i C = i L Conservation of energy V R + V L + V C + V E = 0 And we get V C = 1 ˙ C i L V L = 1 ˙ Li L 0 = V R + Ri E 0 = V E + V R + V C + V L 0 = i L − i E 18 / 40
Engineering examples of DAE - Electrical system Let x = ( V C , V L , V R , i L , i E ) we have 1 0 0 0 0 C 1 0 0 0 0 L x = ˙ 0 0 0 0 0 x 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 R x + V E 0 = 1 1 1 0 0 1 0 0 0 1 − 1 0 which is of the form x = Ax ˙ 0 = Bx + Dz 19 / 40
Method of Lines for PDE Consider the linear PDE (diffusion equations) u ( x = x 0 , t ) = u b ∂ t ( x , t ) = D ∂ 2 u ∂ u ∂ x 2 ( x , t ) with ∂ u ∂ x ( x = x f , t ) = 0 and D a constant. Using method of lines, we have with an equally spaced grid for x (finite difference) ∂ 2 u ∂ x 2 ≈ u i +1 − 2 u i + u i − 1 ∆ x 2 � � + O ∆ x 2 Hence, we get du i dt = D u i +1 − 2 u i + u i − 1 for i = 1 , 2 , . . . , M ∆ x 2 20 / 40
Method of Lines for PDE In other words, we get the system u 1 = u b du 2 dt = D u 3 − 2 u 2 + u b ∆ x 2 du 3 dt = D u 4 − 2 u 3 + u 2 ∆ x 2 . . . du M = D u M +1 − 2 u M + u M − 1 ∆ x 2 dt u M +1 = u M Note u M +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 � x = z ˙ x = z ˙ F (˙ x , x , t ) = 0 ⇔ 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
Recommend
More recommend
