A stochastic approach for computing the domain of attraction without trajectory simulation

We present an approach for the numerical computation of the domain 
of attraction of some asymptotically stable set for continuous-time autonomous 
systems. It is based on a set-oriented approximation of the original 
dynamical system by a Markov jump process. The domain of attraction is 
extracted from absorption probabilities of the jump process. The method does 
not perform any trajectory simulation, integrals of the underlying vector eld 
on the boundary of partition elements are computed instead.


1.
Introduction. In numerous applications, one needs to compute domains of attraction (DOA) of some asymptotically stable set. A naive computation, e.g. creating a grid and simulating all grid points long enough to decide whether they converge to the steady state, can be computationally quite expensive. But even more sophisticated approaches [4,9,11,13] use (forward or backward) simulation of the underlying ordinary differential equation (ODE). Zubov's method [20] (see also in [10]) does not use trajectory simulation, however a partial differential equation has to be solved. Similar holds for Lyapunov function based techniques [7].
In this paper we present a numerical method for approximating the subset of the DOA of continuous-time autonomous systems lying in a previously defined set X, also without any trajectory simulation. Our method is based on a set-oriented approach to approximate properties of dynamical systems [1,3]. We partition the state space X into disjoint sets, which serve to discretize the original deterministic dynamical system, and relate it to a finite state non-deterministic system -a Markov jump process. Then we construct the infinitesimal generator (matrix) of the latter by integrating the vector field of the ODE on the boundary of the partition elements. From this matrix we compute absorption probabilities for the finite state system, which yield the desired information about the DOA.
For discrete-time systems, the idea goes back to [8]; an extensive study of cellto-cell mappings is given in [12].
An outline of this paper is as follows. Section 2 gives a brief introduction to deterministic and non-deterministic dynamical systems. Section 3 shows the mathematical instrument describing statistical properties of dynamical systems -the Frobenius-Perron operator. Section 4 contains the main tool we are working with INFINITESIMAL GENERATOR APPROACH 855 -the associated infinitesimal generator, and describes its numerical approximation. Finally, Section 5 assembles the knowledge gained in the previous sections to develop our approach for the approximation of the DOA for continuous-time systems; and we show a numerical example.

Background.
2.1. Dynamical systems. The temporal variation of states x(t) we are interested in is given by an (autonomous) We are only looking at a compact subset X of the whole state space R d . We call the pair (X, φ t ) a dynamical system, and we refer to X as state space. Objects which stay unchanged under the flow are of particular interest. A set A ⊂ X is called (forward) invariant, if φ t A = A for all t ∈ R ≥0 . The simplest example of an invariant set is a fixed point, i.e. a point x ∈ X such that v(x) = 0. If one already has located an invariant set A, the dynamical behavior in its vicinity may be analyzed -e.g. whether orbits tend to A as t → ∞.
is called the domain of attraction of A.

Markov processes.
We use the theory of Markov processes to describe nondeterministic dynamical behavior. From now on let (X, B, m) be a measure space with the compact state space X and Lebesgue measure m. We also work with a given probability space (Ω, F, P). Definition 2.2 (Stochastic process). Let T = N or T = R ≥0 , and let X be the set of possible states. Then a family {Z t } t∈T , where Z t : Ω → X for all t ∈ T , of X-valued random variables is called a stochastic process. Definition 2.3 (Discrete time Markov chain). We call a stochastic process {Z t } t∈T a discrete time Markov chain, if T = N, and there is a q : x ∈ X and t ∈ T.
A stochastic process {Z t } t∈T is called a finite state discrete time Markov chain, if T = N, X = {x 1 , . . . , x n } is a finite state space, and P ∈ R n×n , with (i) P ij ≥ 0 and (ii) n i=1 P ij = 1, describes the transition probabilities, i.e. P ji = P(Z t+1 = x j |Z t = x i ), for all t ∈ T and i, j ∈ {1, . . . , n}.
The matrix P is called the transition matrix. A matrix satisfying (i) and (ii) is called stochastic.
Definition 2.4 (Continuous time Markov chain). Let T = R ≥0 and state space X = {x 1 , . . . , x n }. Further, let G ∈ R n×n be such that (i) G ij ≥ 0 for i = j, and (ii) Then P t is a stochastic matrix for every t ≥ 0. 1 The process {Z t } t∈T with P t ji = P(Z t+s = x j |Z s = x i ) for all t, s ≥ 0, and i, j ∈ {1, . . . , n} is called a finite state continuous time Markov process. It is often called also a Markov jump process (MJP). The matrix G is called the (infinitesimal) generator of the process.
One may think of such MJPs as follows (cf. Section 2.6 in [18]). Being in an arbitrary state x i at an arbitrary time s, we remain at the given state until some random time s + t, when we jump at random to another state x j . The jump time t has exponential distribution with parameter G ii , and the probability of jumping to state Consider now a finite state discrete time Markov chain with state space X = {x 0 , x 1 , . . . , x n } and transition matrix P . Let x 0 be an absorbing state, i.e. P j0 = δ j0 , where δ ij is the Kronecker delta, and there is another state from which the probability of moving to x 0 is nonzero. If the process enters state x 0 , it will stay there forever. Erasing the state x 0 from X and the corresponding row and column from P , we get a new state space X and transition matrix P , which can be associated with a new process. Since there is a j ∈ {1, . . . , n} with n i=1 P ij < 1, the process may exit the state space X and terminate. We call such processes leaky. For the state space X we think of the fictive state x 0 as everything being outside of X, hence outside the scope of interest.
Analogous considerations can be done with a continuous time process. Its generator will satisfy n i=1 G ij < 0 for some j ∈ {1, . . . , n}.

Absorption probabilities.
Throughout this section we consider MJPs on a finite state space X = {x 1 , . . . , x n }. Furthermore, we assume that there is one absorbing state x 1 , and that the process is possibly leaky. Then, under a reachability assumption the process will either end in the absorbing state x 1 , or "leak out" (i.e. terminate in the fictive state x 0 ; cf. above). 2 Let p ∈ [0, 1] n denote the absorption probabilities for the absorbing state, i.e.
We wish to compute p from the generator G of the process.
(a) The eigenvalue problem method. Fix some τ > 0 and define the discrete time process {Y k } k∈N by Y k := Z kτ . The new process has the transition matrix P = e τ G , and the absorption probabilities are the same as for {Z t } t≥0 . Now, the absorption probabilities for {Y k } k∈N satisfy the well-known one-step relation which can be written as p P = p (cf. also Theorem 1.3.2 in [18]). Hence p is a left eigenvector of P at eigenvalue 1. By P = e τ G it is easy to see that then p G = 0 (cf. also Theorem 3.3.1 in [18]).
(b) The partition method. If x 1 is an absorbing state, the first column of the generator G is the zero vector, and we can partition G as with q ∈ R n−1 and Q ∈ R (n−1)×(n−1) . From this, the transition matrix after time t can be written as Note, that the jth entry of the first row of e tG yields the probability that the MJP starting in x j is in x 1 at time t. Once the process is absorbed in state x 1 it stays there forever, hence for t = ∞ the first row contains the absorption probabilities p j . Trivially, p 1 = 1. Let p = (p 2 , . . . , p n ) . Now, from (2) we obtain where I denotes the identity matrix. Assume that lim t→∞ e tQ = 0 (since this is exactly the transitivity property of the states {x 2 , . . . , x n }, the same reachability assumption as above ensure this). Then, it holds also that Q is non-singular. 3 Solving the equation above, we proved

Remark 1.
There is an analogous derivation of absorption probabilities for discrete time Markov chains, see Section 10.3 in [12], and also [14].
3. Statistical properties of dynamical systems.
3.1. Transfer operators. Transfer operators are used for describing the evolution of statistical properties of the underlying system. For an X-valued random variable Z we write Z ∼ f ("Z is distributed according to f ") for a density f ∈ L 1 (X) if, 4 Given a stochastic process {Z t } t≥0 , we are interested in the distribution of Z t provided Z 0 ∼ f . For the systems considered here, we give the following abstract definition.
Definition 3.1 (Frobenius-Perron operator). Given a Markov process {Z t } t∈T , the bounded linear operator P t : is called the Frobenius-Perron operator (FPO).
For the dynamical systems in consideration, there are explicit representations of the FPO (for (a) and (b) see [18]; for (c) and (d) see [17]): (a) Finite state discrete time Markov chain: One may see easily, that the FPO is the transition matrix of the chain, i.e. P 1 = P . Also, it holds P t = P t for every t ∈ N. (b) Finite state continuous time Markov chain: On analogy to (a), we have P t = e tG for t ≥ 0. (c) General discrete time Markov chain: It holds (d) Continuous time deterministic case: For the special deterministic case, where Z t = φ t x for Z 0 = x ∈ X, and X is an invariant set, one has where Dφ −t is the Jacobian matrix of φ −t . In general, the set X we are looking at will not be invariant, and we do not wish to keep track of any orbit which leaves the set X. Hence, we consider an orbit leaving X as being "lost", and not entering X any more. 5 This change in the dynamics causes a change in the FPO as well. It turns to

Discretization of the FPO as Markov chain.
This section shows how a discretization of a deterministic system can be connected to a non-deterministic system. Partition the state space X into connected disjoint subsets B 1 , . . . , B n ; i.e. n i=1 B i = X. In all our considerations, the B i will be hyperrectangles ("boxes"). Let χ i , i = 1, . . . , n, denote the characteristic function on B i . We define an approximation space V n := span(χ 1 , . . . , χ n ), and represent functions/operators in/on V n w.r.t. the basis (χ 1 1 , . . . , χ 1 n ), where χ 1 i := χ i /m(B i ). Also, we define a projection π n : L 1 (X) → V n onto V n by Let P t n : L 1 (X) → V n be the Galerkin projection of P t on V n , i.e. P t n = π n P t π n . It can be shown that, with the given basis, P t n Vn has a matrix representation P t n with which is the probability of {φ t x ∈ B i }, provided x ∈ B j is chosen at random according to a uniform distribution. The discretization of the FPO yields an approximation of the deterministic process Z t = φ t (with probability 1). It is approximated by a process {Y t n } t≥0 , such that if Y 0 n = y ∈ X, then where j y is the (almost everywhere unique) index with y ∈ B jy [5]. In general, if Y 0 n ∼ f , then Y t n ∼ P t n f . Since π n converges pointwise in L 1 (X) as n → ∞, and P t is a bounded operator, we have pointwise convergence of P t n to P t as n → ∞. This means Y t n converges to Z t in distribution as n → ∞ for every fixed t ≥ 0, whenever the initial distribution is a density.
Note that if no orbits leave X, the transition matrix P t n is stochastic; otherwise the process {Y t n } t≥0 may be leaky, i.e. n i=1 P t ij < 1 for some j. Remark 2. Even though the process {Y t n } t≥0 has state space X, note that it has "essentially" a finite state space, since it can be defined solely by a finite number of box-to-box transition rates -the entries of the transition matrix P t n .
4. The approximate dynamical system.

4.1.
The discrete generator. In this section we follow the aim of designing (by discretization) a continuous time "essentially" finite state Markov process (see Remark 2) which is a suitable approximation of the deterministic dynamical system given byẋ = v(x), in the sense of convergence in distribution, cf. Section 4.2.
We have seen in Section 2.2 that if we have a finite state continuous time Markov chain, then its generator G and time-t-transition matrix P t satisfy the relations P t = e tG , and G = lim where I is the identity matrix. Approximating the dynamics of a dynamical system on a continuous space by a finite state continuous time Markov process is equivalent to approximating the transfer operator P t by a finite dimensional operator of the form e tGn , where G n is a finite rank generator. Since the discretization from Section 3.2 already yields a finite rank approximation P t n of P t , our candidate for the generator will be its time derivative at t = 0; although there is no G ∈ R n×n , in general, such that P t n = e tG would hold. 6 By applying the same notation as in Section 3.2, we define G n : L 1 (X) → V n by G n f := lim t→0 P t n f − π n f t .
Lemma 4.1 shows that G n is well-defined. Note that π n G n = G n = G n π n .
Lemma 4.1 (Lemma 5.11 [16]). For x ∈ ∂B j , define n j (x) to be the unit normal vector pointing out of B j . The boxes B j satisfy that n j exists almost everywhere on ∂B j (measured by the d − 1 dimensional Lebesgue measure m d−1 on ∂B j ). The matrix representation of G n Vn : V n → V n w.r.t. the basis (χ 1 1 , . . . , χ 1 n ) is given by where x · y denotes the dot product of the vectors x and y, and f + denotes the positive part of the function f .
Equation (5) shows that G n,ij is the flux of probability mass from B j to B i if i = j, and −G n,jj is the flux of probability mass leaving B j -provided mass is equidistributed in B j . Thus G n,ij = 0 whenever B i and B j have no common d − 1 dimensional face, hence G n is a sparse matrix.

5.
Approximating the domain of attraction.
5.1. The setting and the approach. Let a dynamical system (X, φ t ) be given by the ODEẋ = v(x). Denote the attractor relative to X with A. Our aim is to approximate the DOA D of A. For this, we use a box-partition {B 1 , . . . , B Nn } of the state space X, and compute the discrete generator G n on this partition, as discussed in Section 4.1. For simplicity, assume that a known subset of the box partition is a box covering A n of the attractor A, such that there is no outward flow on the boundary of this covering. 7 We treat A n as one big state. This means collapsing the states (boxes) in A n into one state, and collapsing (adding up and erasing) the corresponding rows and columns of G n . Thus, we obtain a matrix representation G n of the discrete generator, having only one absorbing state. Finally, we compute the absorption probabilities of the other boxes for this state. Here we use the methods presented in Section 2.3 on G n .
The discrete generator yields a stochastic process, which approximates the deterministic dynamics in the sense of distributions (cf. Section 4.2), hence we expect the absorption probabilities to be high for boxes in the interior of D, we expect them to be low outside of D, and to decrease/increase as we approach the boundary of D from inside/outside.
We can extract an approximate DOA D n from the absorption probability vector p by defining D n to be the union of all boxes B i such that the absorption probability is at least c, for a given probability threshold 0 < c ≤ 1:
The Jacobian of the vector field at the origin has two different negative eigenvalues, hence the origin is asymptotically stable. We consider X = [−1, 1] 2 and search for the DOA of the origin in X. If the partition of X is fine enough, the box containing the origin will not yield outward flow on its boundary, 8 hence the corresponding discrete generator G n will have one absorbing state. Denote the index of this box containing the origin by i 0 . We choose a partition such that this box is unique, i.e. the origin is in its interior.
Denote by p the absorption probabilities for B i0 . We have that p is a solution of p G n = 0 with p i0 = 1 (cf. Section 2.3 method (a)); and also p can be obtained by method (b) in Section 2.3. Both methods give the same results, but method (b) is significantly faster than (a). Table 1  type, where P t n is assembled according to (4), and p is approximated by solving the eigenvalue problem (EVP) P t n p = p. Note, that this latter approach [8] uses trajectory simulation. 9 In Figure 1 we show a phase portrait of the system; and we plot the absorption probabilities computed with the two approaches above. Since we can visualize the absorption probabilities for a 2d system, there is no need to extract an approximate DOA by thresholding the probabilities.
Note the "blurred" edges on the right hand side in Figure 1. These are due to a numerical diffusion introduced by the method (see [16] Chapter 5); the price we pay for being able to omit trajectory simulations. Remark 4. In this example the computation time assembling the discrete generator G n (22 sec) is much larger than the time needed to compute the absorption probabilities p from it (1.2/0.2 sec). This fact falsely suggests that it does not make any big difference whether we use method (a) or (b) to compute p. The generator G n is a sparse matrix, and its computation scales linearly with the number of boxes N n , while computation of p by method (a) uses inverse vector iteration which has numerical costs O N 3 n -it solves one system of linear equations in each iteration step. In contrast, method (b) solves merely one system of linear equations.