We report on recent work [16] on the reachability analysis of nonlinear ordinary diferential equations ( odes). Relying on Carleman linearization and Krylov projection, we describe a method that, given a nonlinear ode system, generates a small linear approximation of the original dynamics. The construction is independent of the initial condition. Used in conjunction with zonotopes, this yields CKR, an accurate reachability analysis algorithm.
for global approximation of nonlinear odes by linear odes, based on polyflows . These are systems that
are exactly linearizable via a change of variables. The technique in [28] is based on building polyflows
that approximate the original system, using as a basis the Lie derivatives up to some order and as
→ +∞ the approximation of [28] becomes exact. Note that this is an asymptotic result that does
not easily yield concrete bounds for a fixed . Systems that are exactly linearizable via polynomial
changes of variables are the subject of [
In [
In the field of reachability for continuous and hybrid systems, state-of-the-art tools like Flow * [23]
and cora [
The investigation reported here is part of a broader research agenda, aimed at integrating formal
methods for dynamical and safety-related systems [
We introduce a linearization method for system (1) that is strongly related to Carleman embedding [29] and then discuss an approach to reduce the dimension of the linearized system. For = (1, ..., ) a vector of state variables, we consider a system of odes
˙ = () where = (1, ..., ) is a vector field of locally Lipschitz analytic functions defined on some open subset Ω ⊆ R. For 0 ∈ Ω, we let (; 0) be the unique solution of the ode system with the initial condition (0) = 0: the unique solution exists and is real analytic (Picard-Lindelöf theorem).
For a real analytic function defined on some open subset of R that includes the trajectories (; 0) for 0 ∈ Ω, we will be interested in studying the observable of the system (2) via , that is the function ∘ (; 0) = ((; 0)). Recall that ℒ () := ⟨∇, ⟩ = ∑︀ =1 · is the Lie derivative of (w.r.t.
ℒ(+1)() := ℒ (ℒ()()). We shall ()() is the -th Lie derivative, defined inductively by ), and ℒ omit the subscript whenever it is understood from the context.
Let us fix a set = { 1, 2, ...} of functions : R → R. For instance might be all monomial functions. We assume that there are unique = ( 1, ..., ) ∈ R and basis vector := ( 1, ..., ) such that
= ∑︁ =
=1 where is any observable function. Otherwise, all we require from the functions in is that they are analytic1, and that the Lie derivative of each can in turn be expressed as a unique linear combination of elements from : ℒ( ) = ∑︀≥ 1 . We let denote the × matrix of the coeficients for 1 ≤ , ≤ , and be the × matrix of possibly nonzero elements , = ,+ ; that is, is chosen large enough to ensure that, for 1 ≤ ≤ , we have = 0 for each > + . We let =△ ( +1, ..., +) . The Carleman linearization (or embedding) of (2) is given by the following linear system in the variables = (1, ..., ) and initial condition
˙ = + ((; 0)) (4) 1This can be weakened to analyticity in some open set containing all the trajectories (; 0) for 0 ∈ Ω. (2) (3) (0) = (0) =: 0 . (5) The following result is an almost immediate consequence of the existence and uniqueness of the solution of odes (Picard-Lindelöf). For a detailed proof, see [8, Th.3].
Theorem 1 (Carleman linearization). Let 0 ∈ Ω. Then ((; 0)) is the unique solution of the system (4) with (0) as in (5).
Note that we cannot explicitly build the system (4), as the function ((; 0)) is in general not available. This leads us to consider an approximation where we neglect the “remainder” ((; 0)), the truncated Carleman linearization of dimension
˙ = (0) = 0(= (0)) . (6)
Now, we discuss a method to reduce the dimension of (6), while keeping certain, still local, accuracy guarantees. Fix = ⟨, ⟩ ( ∈ R ), an observable of interest, as in (3). We consider the -dimensional Krylov space2 generated by and , that is the subspace of R
:= span{, , ( )2, ..., ( )− 1} .
Let = [1| · · · | ] be an orthonormal basis of , represented as a × matrix. Consider the projection of onto and represent it w.r.t. the basis , that is the × matrix := .
Given a vector of distinct state variables = (1, ..., ) , we let the reduced linear system derived from (4) and the corresponding initial condition, derived from (5), be defined as:
˙ = (7) (0) = 0 =: 0 .
Note that the reduced equations (7) do not depend on 0 ∈ Ω. Informally speaking, the solution (; 0) of the reduced system describes the evolution of the vector ((; 0)), projected onto the subspace , in the coordinates of the basis . We note that there exists a well-known algorithm for the eficient, “on the fly” construction of the matrices , , the Arnoldi iteration [34]. Recalling that = ⟨, ⟩ it is natural to consider the following approximation of ((; 0)).
Definition 1 (reduced observable dynamics). For each 0 ∈ Ω and 0 = 0, we define the function: (; 0) := (; 0) .
̂︀ In fact, we will see that 1 = /||||2, while is orthogonal to for > 1. Hence (8) can be simplified to (; 0) = ||||2 1(; 0) .
̂︀ In order to study the quality of this approximation, we introduce the error function relative to (8) (9) (; 0) := ((; 0)) − ̂︀(; 0) .
The following result confirms that this error is small near = 0. Indeed, the Taylor expansions of (; 0) and ((; 0)) up to order − 1 coincide: ̂︀ Theorem 2. For each 0 ∈ Ω, the function (; 0) is () around = 0.
Explicit local bounds of the error function can be obtained from the Taylor theorem with remainder in Lagrange form, assuming we can construct validated enclosures and of ( ; ) and dd ( ; ), ̂︀ respectively, for (, ) ranging in a small set — which is possible by standard techniques, see e.g. [33] and references therein. Under suitable, stability conditions, also a global bound of the error function can be given. Let denote the projection of onto ⊥, the orthogonal complement of , and define the remainder function as ℎ() := () + (). Then, for any > 0 such that ( ; 0) is defined for ∈ [0, ] and assuming additionally is stable, we can prove ∫︁ | (; 0)| ≤ || ||2
|ℎ(( ; 0))| d 0 where > 0 is a constant independent of . Qualitatively speaking, (10) says that, for a stable , the behaviour of the global error is determined by |ℎ(( ; 0))|: if this function decays fast enough to be integrable over [0, +∞), then (; 0) will be globally bounded. (10) 2For an introduction to Krylov spaces, see e.g. [34].
We will apply the outlined linearization scheme to compute an approximation (; 0) of the flow ̂︀ (; 0), and then use it to compute an overapproximation of the reachable set of the nonlinear system (1) at fixed times: 1, 2, .... This goal will be achieved by applying the scheme of Section 2 to each of the observable functions = , for = 1, ..., in turn. Using the notation in that section, for each = 1, ..., , let () the coeficient vector of in the chosen basis , that is = () , and (), () the corresponding basis and reduced matrix. We define the approximate flow by ̂︀(; 0) := (̂︀1(; 0), ..., (; 0)) where, as an instance of (9), we have ̂︀
̂︀(; 0) := ||()||2 1()(; 0) ( = 1, ..., ) (11) with ()(; 0) the solution of the linear initial value problem (7) for = . As the solution of a linear system of odes, each component in (11) can be written explicitly as ̂︀(; 0) = ||()||2 ︀( e() ︀) 1 () (0) for = 1, ..., , where (e(· ))1 denotes the first row of the exponential matrix. Note that, as a function of 0, for a fixed , (; 0) is a linear combination of the components of the ̂︀ basis (0). It is also convenient to introduce the following error vector: (; 0) := (; 0) − ̂︀(; 0).
In what follows, we will consider the case where an initial set 0 is given, rather than an individual initial state 0. The general idea of the algorithm is to use ̂︀(Δ; · ) to propagate a reachset from one time point to the next, by a time-step of Δ. The propagated set needs then to be ‘inflated’ to compensate for approximation errors. The concrete way in which propagation and inflation are carried out will depend on the representation that will be adopted for sets of reachable states. We will give a generic description of the method, independent of the type of representation. Then we will instantiate it to a concrete method by considering a specific set representations, zonotopes [27].
Let 0 be a compact set of initial states and > 0 be such that ( ; 0) is well-defined for each 0 ∈ 0 and ∈ [0, ]. In the interval [0, ], choose time points 0 = 0, 1, · · · , = , with Δ := − − 1 > 0 for 1 ≤ ≤ . The algorithm, which we christen ckr for Carleman-Krylov Reachability, builds a sequence of compact sets 0, 1, ..., ⊆ R, the reachsets, s.t. (a) is an overapproximation of the set of reachable states at time : ⊇ (; 0) = {(, ) : ∈ 0}, in particular 0 = 0; (b) belongs to a pre-specified class of compact subsets of R, say (e.g. polytopes, zonotopes,...). The method essentially applies three set operations on R to build given − 1: • propagation: propagates linearly a reachset − 1 from − 1 to , thus obtaining ̃︀; • enclosure computation: generates a compact set that includes the error vector at ; • inflation : computes a bloated version of the propagated reachset ̃︀ that includes ̃︀ + , where + denotes here Minkowski sum3.
A formal description of these operations, as well as a proof of correctness of the general method, can be
found in the journal version [
Now, let us consider a specific set representation: zonotopes. Given a column vector ∈ R (center)
and a matrix formed by ≥ 1 column vectors = [1| · · · | ] ∈ R× (generators), the corresponding
3For , ⊆ R, + := { + : ∈ , ∈ }.
zonotope is the set Z(, ) := { + · : ∈ [
We quantify the accuracy of the returned sequence of reachsets using two diferent measures: (1) the
volume of the final reachset ; (2) only for cora and CKR, the average volume of the reachsets at times
0, ..., = . ckr is extremely accurate for all the models we have considered. In particular the final
reachsets produced by ckr are tighter than those produced by the other two tools in all cases, often
significantly so. With one exception, ckr is the most accurate algorithm, also when considering average
accuracy across diferent examples. In terms of execution time, Flow * is the most efective tool on the
considered examples. ckr times are in line with, or at least comparable with, Flow* ’s. As an example,
we report in Fig. 2 a graphical comparison of the reachsets produced by the three algorithms when
applied to the Van der Pol system. Numerical values and full details on experiments can be found in
[
10 8 6 4 2 0 -2 -4-3 3 2 1 y0 −1 −2 3 2 1 0 −1 −2 −3 −2 -2 -1 0 1 2 3 Acknowledgments Work partially supported by the project SERICS (PE00000014), EU funded NRRP MUR program - NextGenerationEU.
The authors have not employed any Generative AI tools. [21] X. Chen. Reachability analysis of non-linear hybrid systems using Taylor models. Ph.D. thesis, RWTH
Aachen University, 2015. [22] X. Chen, E. Abraham, and S. Sankaranarayanan. Taylor Model Flowpipe Construction for Non-linear Hybrid Systems. In Real-Time Systems Symposium (RTSS), pp. 183-192, IEEE, 2012. 10.1109/RTSS.2012.70. [23] X. Chen, E. Abraham, and S. Sankaranarayanan. Flow* : An Analyzer for Non-linear Hybrid Systems. In Computer Aided Verification (CAV) , LNCS 8044, pp. 258-263, Springer, 2013. 10.1007/9783-642-39799-8_18. [24] A. Chutinan, B. H. Krogh. Computational techniques for hybrid system verification. IEEE Transactions on Automatic Control 48(1): 64-75, 2003. 10.1109/TAC.2002.806655. [25] M. Forets, A. Pouly. Explicit Error Bounds for Carleman Linearization. arXiv:1711.02552, 2017. [26] M. Forets, C. Schilling. Reachability of Weakly Nonlinear Systems Using Carleman linearization.
In Proc. of Reachability Problems 2021, LNCS 13035:85–99, 2021. 10.1007/978-3-030-89716-1_6. [27] A. Girard. Reachability of Uncertain Linear Systems Using Zonotopes. In: Morari, M., Thiele, L. (eds) Hybrid Systems: Computation and Control. HSCC 2005. LNCS 3414: 291–305, Springer, 2005. 10.1007/978-3-540-31954-2_19. [28] R. M. Jungers, P. Tabuada. Non-local Linearization of Nonlinear Diferential Equations via Polyflows.
2019 American Control Conference (ACC), pp. 1-6, 2019. 10.23919/ACC.2019.8814337. [29] K. Kowalski, W.-H. Steeb. Nonlinear Dynamical Systems and Carleman Linearization. World Scientific, 1991. 10.1142/1347. [30] K. Makino and M. Berz. Rigorous integration of flows and ODEs using Taylor models. In Proceedings of Symbolic-Numeric Computation, pages 79–84, ACM, 2009. 10.1145/1577190.1577206. [31] H.D. Menéndez, M. Boreale, D. Gorla, D. Clark. Output Sampling for Output Diversity in Automatic Unit Test Generation. IEEE Transactions on Software Engineering 48(2): 295-308, 2022. 10.1109/TSE.2020.2987377. [32] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin. A time-dependent Hamilton–Jacobi formulation of reachable sets for continuous dynamic games. IEEE Transactions on Automatic Control, 50:947–957, 2005. 10.1109/TAC.2005.851439. [33] N. S. Nedialkov, K. R. Jackson, G. F. Corliss. Validated solutions of initial value problems for ordinary diferential equations. Applied Mathematics and Computation. 105(1): 21-68, Elsevier, 1999. /10.1016/S0096-3003(98)10083-8. [34] Y. Saad. Iterative methods for sparse linear systems. SIAM, 2003. 10.1137/1.9780898718003. [35] S. Sankaranarayanan. Automatic abstraction of non-linear systems using change of bases transformations. roceedings of the 14th international conference on Hybrid systems: computation and control (HSCC): 143-152, 2011. 10.1145/1967701.196772. [36] S. Sankaranarayanan. Change-Of-Bases Abstractions for Non-Linear Systems. CoRR abs/1204.4347, 2012. [37] A. Tiwari, G. Khanna. Nonlinear systems: Approximating reach sets. In Hybrid Systems: Computation and Control (HSCC):600-614, ACM, 2004. 10.1007/978-3-540-24743-2_40. [38] B. Xue, M. Fränzle, N. Zhan. Under-Approximating Reach Sets for Polynomial Continuous Systems.
Proceedings of the 21st International Conference on Hybrid Systems: Computation and Control (part of CPS Week) (HSCC), pp. 51–60, ACM, 2018. 10.1145/3178126.3178133.