Reachability Analysis for High-Index Linear Differential Algebraic - - PowerPoint PPT Presentation

Institute for Software Integrated Systems

Vanderbilt University

Reachability Analysis for High-Index Linear Differential Algebraic Equations (DAEs)

17th International Conference on Formal Modeling and Analysis of Timed Systems (FORMATS’19), August 27, 2019 Hoang-Dung Tran, Luan Viet Nguyen, Nathaniel Hamilton, Weiming Xiang & Taylor T. Johnson VeriVITAL-The Verification and Validation for Intelligent and Trustworthy Autonomy Laboratory (http://www.verivital.com) Electrical Engineering and Computer Science (EECS)

https://github.com/verivital/daev/

Motivation: Mass Dampers

2

[Intro to Structural Motion Control, Connor 2003]

3

Motivation

▪ Most existing cyber-physical systems (CPS) verification techniques only focus

• n physical behaviors as ordinary differential equations (ODEs), or hybrid

variants thereof (hybrid automata, etc.) ▪ Many CPS domains naturally model systems as DAEs instead of ODEs ▪ Mechatronics, robotics, electrical circuits, earthquake engineering, water distribution networks / fluid dynamics (certain problems), process/chemical engineering, …

Index-3 DAE system electrical generator (power) Index-2 interconnected rotating masses (IRM) system (automotive) Index-2 semi-discretized Stoke System (fluids) Index-3 damped mass-spring system (earthquake)

DAE Modeling Intuition

▪ Consider an RLC (resistor, inductor, capacitor) circuit ▪ Kirchhoff's current law (KCL) and voltage law (KVL) => algebraic constraints + ODEs for transient behavior

▪ KCL: conservation of current: 𝑗𝐹 = 𝑗𝑆 = 𝑗𝐷 = 𝑗𝑀 ▪ KVL: conservation of energy: 𝑊

𝑆 + 𝑊 𝐷 + 𝑊 𝑀 + 𝑊 𝐹 = 0

▪ Ohm’s laws: C ሶ 𝑊

𝐷 = 𝑗𝑑

L ሶ 𝑊

𝑀 = 𝑗𝑀

𝑊

𝑆 = 𝑆 𝑗𝑆

4

DAE Modeling Intuition

▪ Replace equal currents (𝑗𝑆 to 𝑗𝐹, 𝑗𝐷 to 𝑗𝑀), don’t have to, but reduces dimensionality for fewer state variables

ሶ 𝑊

𝐷 = 1 𝐷 𝑗𝑀

ሶ 𝑊

𝑀 = 1 𝑀 𝑗𝐹

0 = 𝑊

𝑆 + 𝑆𝑗𝐹

0 = 𝑊

𝐹 + 𝑊 𝑆 + 𝑊 𝐷 + 𝑊 𝑀

0 = 𝑗𝑀 − 𝑗𝐹

▪ Now a DAE system with:

5

𝑦 𝑢 = 𝑊

𝐷(𝑢)

𝑊

𝑀(𝑢)

𝑊

𝑆(𝑢)

𝑗𝑀(𝑢) 𝑗𝐹(𝑢)

DAE Modeling Intuition

▪ Linear DAE system: 𝑒𝑦 𝑒𝑢 = ሶ 𝑦 = 𝐵𝑦 0 = 𝐶𝑦 + 𝐸𝑨

6

𝑦 𝑢 = 𝑊

𝐷(𝑢)

𝑊

𝑀(𝑢)

𝑊

𝑆(𝑢)

𝑗𝑀(𝑢) 𝑗𝐹(𝑢) , 𝑨 𝑢 = 𝑊

𝐹(𝑢)

𝐶 = 1 𝑆 1 1 1 1 −1 𝐸 = 1

ሶ 𝑊

𝐷 = 1 𝐷 𝑗𝑀

ሶ 𝑊

𝑀 = 1 𝑀 𝑗𝐹

0 = 𝑊

𝑆 + 𝑆𝑗𝐹

0 = 𝑊

𝐹 + 𝑊 𝑆 + 𝑊 𝐷 + 𝑊 𝑀

0 = 𝑗𝑀 − 𝑗𝐹

𝐵 = 1 𝐷 1 𝑀

7

Motivation

▪ Most existing cyber-physical systems (CPS) verification techniques only focus

• n ODE dynamics, or hybrid variants thereof (hybrid automata, etc.)

▪ Verifying DAE systems is more complex than ODE systems ▪ No existing works (to our knowledge) on verifying high-index (>1) DAEs ▪ Scalability: state-space explosion / “curse of dimensionality” ▪ How to verify safety of systems with DAE dynamics?

Index-3 DAE system electrical generator (power) Index-2 interconnected rotating masses (IRM) system (automotive) Index-2 semi-discretized Stoke System (fluids) Index-3 damped mass-spring system (earthquake)

8

Linear DAE Systems

▪ Linear DAE System: 𝑭 ሶ

𝑦 𝒖 = 𝑩𝒚 𝒖 + 𝑪𝒗 𝒖

▪ 𝑦 𝑢 ∈ R𝑜 is the state vector ▪ 𝑣 𝑢 ∈ R𝑛 is the s input vector ▪ 𝐹, 𝐵 ∈ R𝑜×𝑜 and 𝐶 ∈ R𝑜×𝑛 are the DAEs matrices, where 𝐹 is singular (non- invertible) ▪ Index of a DAE: typically (can depend on initial conditions) the minimum number of times to differentiate DAEs wrt 𝑢 to get ODEs (“index reduction”), where ODEs are called index-0, can typically evaluate rank(E) to check

▪ Example: Index-2 interconnected rotating masses (IRM) system

Where 𝐾1 = 1, 𝐾2 = 2, 𝑁2 𝑢 + 𝑁3 𝑢 = 0, 𝑨1 𝑢 = 𝑨2(𝑢)

9

Linear DAE Systems

▪ Index-2 interconnected rotating masses (IRM) system Reachable sets computed using daev: https://github.com/verivital/daev

• 1. Decoupling

2. Consistency Checking

▪ Define a consistent space for the initial state and input ▪ Guarantee a solution for the DAE system

• 3. Construct reachable set for the decoupled system

▪ Using Star-sets and Simulation

• 4. Construct reachable set for original DAE system
• 5. Perform safety verification & falsification using computed reachable set

10

Our Approach

DAE AEs 𝐹 ሶ 𝑦 = 𝐵𝑦 + 𝐶𝑣 ODE ODEs ሶ 𝑦1 = 𝑂1𝑦1 + 𝐶𝑣 AC: : Alg Algebraic Con Const straints ts ሶ 𝑦𝑗 = 𝑂𝑗𝑦𝑗 + 𝑁𝑗𝑣

+ =

Mar arz Dec ecoupli ling

11

Index-1 Decoupling

▪ Definition (Tractability index). Assume that the DAE system 𝐹 ሶ 𝑦 𝑢 = 𝐵𝑦 𝑢 + 𝐶𝑣(𝑢) is solvable, i.e., the matrix pair (𝐹, 𝐵) is regular. A matrix chain is defined by: 𝐹0 = 𝐹, 𝐵0 = 𝐵 𝐹

𝑘+1 = 𝐹 𝑘 − 𝐵𝑘𝑅𝑘, 𝐵𝑘+1 = 𝐵𝑘𝑃 𝑘, 𝑘 ≥ 0, where 𝐹 𝑘𝑅𝑘 = 0, 𝑅𝑘 2 = 𝑅𝑘, 𝑄 𝑘 = 𝐽𝑜 − 𝑅𝑘

Where ∃ index 𝜈 s.t. 𝐹𝜈 is non-singular and ∀𝑘 ∈ 0, 𝜈 − 1 , 𝐹𝑘 is singular 𝜈 is called the tractability index A matrix pair (𝐹, 𝐵) is regular if det 𝑡𝐹 − 𝐵 ≠ 0 ▪ Lemma 1 (Index-1 DAE decoupling). An index-1 DAE system can be decoupled using the matrix chain defined as follows: Δ1: ሶ 𝑦1 𝑢 = 𝑂1𝑦1(𝑢) + 𝑁1𝑣(𝑢), ODE subsystems Δ2: ሶ 𝑦2 𝑢 = 𝑂2𝑦1(𝑢) + 𝑁2𝑣(𝑢), AC subsystems 𝑦 𝑢 = 𝑦1 𝑢 + 𝑦2(𝑢) 𝑦1 𝑢 = 𝑄0𝑦 𝑢 , 𝑂1 = 𝑄0𝐹1

−1𝐵0, 𝑁1 = 𝑄0𝐹1 −1𝐶

𝑦2 𝑢 = 𝑅0𝑦 𝑢 , 𝑂2 = 𝑅0𝐹1

−1𝐵0, 𝑁2 = 𝑅0𝐹1 −1𝐶

12

Index-2 Decoupling

▪ Lemma 2 (Index-2 DAE decoupling). An index-2 DAE system can be decoupled using the matrix chain defined as follows:

▪ Intuition: basically taking derivatives wrt 𝑢 of the algebraic constraint subsystems to get ODEs ▪ Scalability issue: increasing dimensionality, more state variables being introduced

Δ1: ሶ 𝑦1 𝑢 = 𝑂1𝑦1(𝑢) + 𝑁1𝑣(𝑢), ODE subsystems Δ2: ሶ 𝑦2 𝑢 = 𝑂2𝑦1(𝑢) + 𝑁2𝑣(𝑢), AC subsystems 1 𝑦 𝑢 = 𝑦1 𝑢 + 𝑦2 𝑢 + 𝑦3 𝑢 𝑦1 𝑢 = 𝑄0𝑄

1𝑦 𝑢 , 𝑂1 = 𝑄0𝑄 1𝐹2 −1𝐵2, 𝑁1 = 𝑄0𝑄 1𝐹2 −1𝐶

𝑦3 𝑢 = 𝑅0𝑦 𝑢 , 𝑂3 = 𝑅0𝑄

1𝐹2 −1𝐵2, 𝑁3 = 𝑅0𝑄 1𝐹2 −1𝐶, 𝑀3 = 𝑅0𝑅1

Δ3: ሶ 𝑦3 𝑢 = 𝑂3𝑦1 𝑢 + 𝑁3𝑣 𝑢 + 𝑀3 ሶ 𝑦2 𝑢 , AC subsystems 2 𝑦2 𝑢 = 𝑄0𝑅1𝑦 𝑢 , 𝑂2 = 𝑄0𝑅1𝐹2

−1𝐵2, 𝑁2 = 𝑄0𝑅1𝐹2 −1𝐶

13

Index-3 Decoupling

▪ Lemma 3 (Index-3 DAE decoupling). An index-3 DAE system can be decoupled using the matrix chain defined as follows: Δ4: ሶ 𝑦4 𝑢 = 𝑂4𝑦1 𝑢 + 𝑁4𝑣 𝑢 + 𝑀4 ሶ 𝑦3 𝑢 + 𝑎4 ሶ 𝑦2 𝑢 , AC subsystems 3 Δ1: ሶ 𝑦1 𝑢 = 𝑂1𝑦1(𝑢) + 𝑁1𝑣(𝑢), ODE subsystems Δ2: ሶ 𝑦2 𝑢 = 𝑂2𝑦1(𝑢) + 𝑁2𝑣(𝑢), AC subsystems 1 𝑦 𝑢 = 𝑦1 𝑢 + 𝑦2 𝑢 + 𝑦3 𝑢 + 𝑦4 𝑢 𝑦1 𝑢 = 𝑄0𝑄

1𝑄2𝑦 𝑢 , 𝑂1 = 𝑄0𝑄 1𝑄2𝐹3 −1𝐵3, 𝑁1 = 𝑄0𝑄 1𝑄2𝐹3 −1𝐶

𝑦3 𝑢 = 𝑄0𝑅1𝑦 𝑢 , 𝑂3 = 𝑄0𝑅1𝑄2𝐹3

−1𝐵3, 𝑁3 = 𝑄0𝑅1𝑄2𝐹3 −1𝐶, 𝑀3 = 𝑄0𝑅1𝑅2

Δ3: ሶ 𝑦3 𝑢 = 𝑂3𝑦1 𝑢 + 𝑁3𝑣 𝑢 + 𝑀3 ሶ 𝑦2 𝑢 , AC subsystems 2 𝑦2 𝑢 = 𝑄0𝑄

1𝑅2𝑦 𝑢 , 𝑂2 = 𝑄0𝑄 1𝑅2𝐹3 −1𝐵3, 𝑁2 = 𝑄0𝑄 1𝑅2𝐹3 −1𝐶

𝑦4 𝑢 = 𝑅0𝑦 𝑢 , 𝑂3 = 𝑅0𝑄

1𝑄2𝐹3 −1𝐵3, 𝑁4 = 𝑅0𝑄 1𝑄2𝐹3 −1𝐶, 𝑀4 = 𝑅0𝑅1, 𝑎4 = 𝑅0 𝑄 1𝑅2

14

Admissible Projectors

▪ Why is it needed?

15

Example: Decoupling for IRM System

▪ Consistent initial set of states ▪ IRM can be decoupled into one ODE and two AC subsystems

16

Consistency Checking

▪ To guarantee a solution for the DAE system, the initial states and inputs must satisfy the following conditions ▪ Where input 𝑣(𝑢) is smooth such that: ሶ 𝑣 𝑢 = 𝐵𝑣𝑣 𝑢 , 𝑣 0 = 𝑣0 ∈ U0

▪ 𝐵𝑣 ∈ R𝑛×𝑜: user-defined input matrix ▪ 𝑉0: the set of initial inputs

Index-1 DAE: 𝑦2 0 = 𝑂2𝑦1(0) + 𝑁2𝑣(0) Index-2 DAE: 𝑦2 0 = 𝑂2𝑦1 0 + 𝑁2𝑣 0 𝑦3 0 = 𝑂3𝑦1 0 + 𝑁3𝑣 0 + 𝑀3 ሶ 𝑦2 0 Index-3 DAE: 𝑦2 0 = 𝑂2𝑦1 0 + 𝑁2𝑣 0 𝑦3 0 = 𝑂3𝑦1 0 + 𝑁3𝑣 0 + 𝑀3 ሶ 𝑦2 0 𝑦4 0 = 𝑂4𝑦1 0 + 𝑁4𝑣 0 + 𝑀4 ሶ 𝑦3 0 + 𝑎4 ሶ 𝑦2 0

17

Consistency Checking

▪ Definition (Consistent space). Consider the DAE system Δ: 𝐹 ሶ 𝑦 𝑢 = 𝐵𝑦 𝑢 + 𝐶𝑣(𝑢), by letting 𝑣 𝑢 = 0, we define a consistent matrix Γ as: ▪ An initial state 𝑦0 is consistent if it is in the consistent space, i.e., Γ𝑦0 = 0 Index-1 Δ : Γ = 𝑅0 − 𝑂2𝑄 Index-2 Δ : 𝑄0𝑅1 − 𝑂2𝑄0𝑄

1

𝑅0 − 𝑂3 + 𝑀3𝑂2𝑂1 𝑄0𝑄

1

Index-2 Δ : 𝑄0𝑄

1𝑅2 − 𝑂2𝑄 0𝑄 1𝑄2

𝑄0𝑅1 − 𝑂3 + 𝑀3𝑂2𝑂

1 𝑄0𝑄 1𝑄2

𝑅0 − 𝑂4 + 𝑀4 𝑂3𝑂1 + 𝑀3𝑂2𝑂1

2 + 𝑎4𝑂2𝑂 1 𝑄0𝑄 1𝑄2

Then, 𝐿𝑓𝑠 Γ is the consistent space of the system Δ, also denotes null space of the matrix Γ

18

Reachability Analysis

▪ Definition (Modified Star-Set). A modified star set Θ is a tuple 〈𝑊, 𝑄〉, where 𝑊 = 𝑤1, 𝑤2, … , 𝑤𝑙 ∈ R𝑜×𝑙 is a star basis matrix and 𝑄is a linear

• predicate. The set of states represented by the star is given by:

Θ = {𝑦|𝑦 = Σ𝑗=1

𝑙

𝛽𝑗𝑤𝑗 = 𝑊 × 𝛽, 𝑄 𝛽 ≙ 𝐷𝛽 ≤ 𝑒} where, 𝛽 = [𝛽1= 1, 𝛽2, … , 𝛽𝑙]𝑈, 𝐷 ∈ R𝑞×𝑙, 𝑄 ∈ R𝑞 , and 𝑞 is the number of linear constraints. 𝑊 = 0 1 1 C= 1 −1 1 −1 1 −1 𝑒 = 1 −1 2 1 1 1

[Stanley Bak, Hoang-Dung Tran, Taylor T. Johnson, "Numerical Verification of Affine Systems with Up to a Billion Dimensions“, HSCC’19] [Hoang-Dung Tran, Patrick Musau, Diego Manzanas Lopez, Xiaodong Yang, Luan Viet Nguyen, Weiming Xiang, Taylor T. Johnson, "Star-Based Reachability Analysis for Deep Neural Networks", FM’19]

19

Reachability Analysis

▪ Lemma 4 (Reachable Set Construction). Given an autonomous DAE system 𝐹 ሶ 𝑦 𝑢 = 𝐵𝑦 𝑢 + 𝐶𝑣(𝑢) where 𝑣 𝑢 = 0 and a consistent initial set of states Θ(0) = 〈𝑊(0), 𝑄〉, let Θ1(t) be the reachable set at time 𝑢 of the corresponding ODE subsystem after decoupling. Then, the reachable set at time 𝒖 of the system is given by Θ t = 〈𝑊 𝑢 = 𝛺𝑊

1 𝑢 , 𝑄〉, where 𝛺 is a reachable set projector

defined as ▪ Recall 𝑂𝑗, 𝑀𝑘, 𝑎𝑙 are from Marz decoupling discussed earlier Index-1: 𝛺 = 𝐽𝑜 + 𝑂2 Index-2: 𝛺 = 𝐽𝑜 + 𝑂2 + 𝑂3 +𝑀3𝑂2𝑂1 Index-3: 𝛺 = 𝐽𝑜 + 𝑂2 + 𝑂3 +𝑀3𝑂2𝑂1 + 𝑀4𝑂3𝑂1 + 𝑀4𝑀3𝑂2𝑂1

2 + 𝑎4𝑂2𝑂1

20

Reachability Analysis

21

Bounded-time safety verification/falsification

22

Reachability Analysis for IRM System

▪ Sinusoid input ▪ A consistent initial set of states ▪ Safety verification w.r.t unsafe specification 𝑁2 𝑢 ≤ −0.8

Reachable set

Violation

An unsafe trace

23

Scalability Performance

Takeaways:

• Daev is scalable in verifying

large DAE systems (≥ 1K state variables) where other over- approximation approaches not applicable

• Daev can produce unsafe traces
• Available:

https://github.com/verivital/daev https://github.com/verivital/daev /releases/tag/formats2019

Benchmark details: ARCH’18 paper, “Linear Differential- Algebraic Equations”

RLC Circuit

24

Damped Mass Spring

25

Damped Mass Spring

26

Partial Element Equivalent Circuit (PEEC)

27

▪ Electromagnetics application: RF engineering

Stokes

28

Stokes

29

30

Scalability Analysis

▪ Stokes-equation PDE D-T: decoupling time, RSC-T: reachable set computation time CS-T: checking safety time V-T: verification time (overall total time sum)

Takeaway:

• Decoupling and reachable set computation

times dominate the time for verification process

• Time for checking safety is almost

unchanged and very small

• vode, dopri5, and dop853 solvers should

be used for large DAE systems Boundary conditions => algebraic constraints (Finite-difference method based

• n marker-and-cell [MAC])

31

Conclusion and future works

▪ Conclusion ✓ A simulation-based reachability analysis for high-index, linear DAE systems ✓ Based on the effective combination of a decoupling method and a reachable set computation using star-sets ✓ Design and implementation of the approach in a Python toolbox, called Daev: https://github.com/verivital/daev/ ✓ Applied to verify/falsify high-index linear DAE systems ✓ Approach can deal with DAE systems with up to thousands of state variables ▪ Future Work ✓ Enhance the time performance and the scalability of our approach ✓ Apply to verify million-dimensional DAE systems ✓ DAEs with hybrid/switching behavior (time or state-dependent)

32

Thank You

32

33

Thank You! Questions?

• Students

– VU EECS: Hoang-Dung Tran (PhD), Nate Hamilton (PhD), Ayana Wild (PhD), Patrick Musau (PhD), Xiaodong Yang (PhD), Ran Hao (PhD), Tianshu Bao (PhD), Diego Manzanas (PhD), Weiming Xiang (Postdoc), Joel Rosenfeld (Postdoc) – UTA CSE:Shafiul Chowdhury (PhD) – UTA Alumni: Luan Viet Nguyen (PhD), Omar Beg (PhD), Nathan Hervey (MS), Ruoshi Zhang (MS), Shweta Hardas (MS), Randy Long (MS), Rahul (MS), Amol (MS)

• Recent Collaborators

– Vanderbilt: Gabor Karsai, Xenofon Koutsoukos, Janos Sztipanovits, … – UTA: Ali Davoudi, Christoph Csallner, Matt Wright, Steve Mattingly, Colleen Casey – Illinois: Sayan Mitra, Marco Caccamo Lui Sha, Amy LaViers – AFRL: Stanley Bak and Steven Drager – Toyota: Jim Kapinski, Xiaoqing Jin, Jyo Deshmukh, Ken Butts, Issac Ito – Waterloo: Sebastian Fischmeister – Toronto: Andreas Veneris – ANU: Sergiy Bogomolov – UTSW: Ian White, Victor Salinas, Rama Ranganathan

Taylor T. Johnson http://www.TaylorTJohnson.com Taylor.Johnson@vanderbilt.edu http://www.verivital.com