Abstract
We study stability for dynamical systems specified by autonomous stochastic differential equations of the form , with an -valued Itô process and an -valued Wiener process, and the functions and are Lipschitz and vanish at the origin, making it an equilibrium for the system. The concept of asymptotic stability in probability of the null solution is well known and implies that solutions started arbitrarily close to the origin remain close and converge to it. The concept therefore pertains exclusively to system properties local to the origin. We wish to address the matter in a more practical manner: Allowing for a (small) probability that solutions escape from the origin, how far away can they then be started? To this end we define a probabilistic version of the basin of attraction, the -BOA, with the property that any solution started within it stays close and converges to the origin with probability at least . We then develop a method using a local Lyapunov function and a nonlocal one to obtain rigid lower bounds on -BOA.
1. Introduction
Lyapunov theory of stability [1] for deterministic systems has been very successful as a tool to study qualitative system behavior. For mathematically oriented studies of the Lyapunov theory, compare, for example, [2–4] or [5–7] for a more modern treatment including complete Lyapunov functions. For more application oriented surveys, compare, for example, [8–10]. The centerpiece of the theory is the Lyapunov function, a real function of state-space which is nonincreasing along solution trajectories. A deterministic Lyapunov function gives rigid estimates on the basin of attraction (BOA) of an equilibrium through its sublevel sets. The equilibrium under consideration can be assumed to be the origin without loss of generality and one then often speaks of the stability of the null solution. Its usefulness has sparked the development of methods to find Lyapunov functions. Breaking the system domain up and addressing each part separately can sometimes work in the deterministic case. With , one starts by generating a local Lyapunov function on the small neighbourhood of the origin and then finding a nonlocal Lyapunov function , using entirely different methods suited away from the equilibrium.
In this paper (see Section 2), we shall develop the theory necessary to enable such a splitting for a stochastic system. Advanced numerical methods taking advantage of the splitting are in the final stages of development by the authors and their collaborators; see [11] for first results on the local part at the equilibrium. An example using analytical solution at the origin and the shooting method further away is presented in Section 3 in this paper.
To set the stage, we shall start by discussing the local and nonlocal Lyapunov functions in the deterministic setting. For definition of the notations used, see Section 1.1. Consider a deterministic system given by an autonomous ordinary differential equation (ODE) , where is sufficiently smooth and such that . To characterize asymptotic stability of the null solution it is convenient to resort to the function classes and (defined in Section 1.1). The null solution is called (locally) asymptotically stable (AS) if there exists a function and a neighbourhood of the origin such that for every and we have , where is the solution to the ODE started at at . In particular, is a subset of the equilibrium’s BOA defined as . If one can even take the origin is called globally asymptotically stable (GAS). Obviously (local) AS is too weak a condition and GAS is an unnecessary strict condition for most applications. One is interested in AS where the set is of reasonable size for the problem at hand. A rigorous lower bound on the BOA can be made if a Lyapunov function for the system is known, that is, a continuously differentiable function defined on an open neighbourhood of the origin and such that and for some and all . For a Lyapunov function all sublevel sets , , that are compact subsets of are subsets of the BOA.
For the general deterministic setting, constructing a Lyapunov function on a reasonably sized set is usually a very difficult problem. Often, however, the construction of a Lyapunov function on a small neighbourhood of the origin is simple. Indeed, if the spectrum of the Jacobian of at the origin does not have any purely imaginary points, then the origin is AS iff the equation has a positive definite solution for a positive definite . Indeed, then is a Lyapunov function for the ODE on and on some neighbourhood of the origin it is a Lyapunov function for the nonlinear system . The size of , however, depends strongly on the nonlinearities of ; compare, for example, [12, p. 659] for an explicit estimate. By using this or a different method to compute a local Lyapunov function for the system , one can compute a set that is a subset of the equilibrium’s BOA. To assert that a larger set is in the BOA, it now suffices to show that whenever , because the local Lyapunov function and the semigroup property of the flow deliver the missing part. For many methods to numerically compute Lyapunov functions for the system on a reasonably sized it is either advantageous or even necessary to use the fact that a (small) set is already known to be in the BOA; compare, for example, [12–18] or [19, Section ] and the references therein for an overview.
Now, consider a set of stochastic differential equations (SDEs) of the form: , and , with equilibrium at the origin; that is, and . As in the deterministic case, a Lyapunov function is a function defined on some neighbourhood of the origin, such that for some . The condition for deterministic systems that is decreasing along solution trajectories, that is, for an , is replaced by the condition that the expectation of the process is a decreasing function of time. Here is the stochastic process started at at time . By Itô’s formulawith where denotes the -th element of the matrix . The expectation of the second integral in the Itô formula is zero because the process is adapted and therefore the decrease of the expectation, that is, , is implied by . In order to find such a Lyapunov function for a stochastic system, one could go about as in the deterministic case and solve the partial differential equation (PDE) for an . However, in the deterministic case, the solution to the PDE will automatically have local minima at the AS equilibrium but this is not the case for a solution to . The PDE is markedly different, and the condition alone is no longer sufficient to ensure that the Lyapunov function has a local minimum at the equilibrium at the origin. It is therefore a natural first approach to try and solve the PDE imposing boundary conditions and on the boundary . The problem with this approach is that the PDE cannot be strictly elliptic at the equilibrium at the origin and the standard theory for existence and uniqueness of solutions to PDE does not apply. Making a cut-out of some small neighbourhood containing the origin from the domain and solving the PDE on instead can ensure strict ellipticity and hence existence and uniqueness of the solution; compare, for example, [20, Theorem 6.14], and better yet, solving so that solutions satisfy both a minimum and a maximum principle (of elliptic PDE) will guarantee that the minima and the maxima of are taken on the boundary and with on and on we automatically have on .
Even if the PDE is not strictly elliptic, as will happen if , the condition instead of means that one has more room to actually compute numerically a true Lyapunov function, rather than a numerical approximation of a function we know is a Lyapunov function, as done in [21] where solutions to the PDE are interpreted in the viscosity sense. This will be discussed briefly in Section 4.
Despite these advantages, the cut-out of course removes the equilibrium itself from the domain of the Lyapunov function and one must take care off that the null solution is actually stable. In Section 2, we describe the details of how this can be achieved by using a different (local) Lyapunov function on a superset of the cut-out set . In Section 3, we show an example of our approach and in Section 4, we draw a roadmap of how the results in Section 2 can be used to develop a general method to give rigid estimates of -BOA of reasonable size for practical applications.
1.1. Notation
We denote a (column)-vector in in boldface and its transpose by . We also write the matrix-valued function in boldface. We denote the Euclidian norm of a vector by and for a matrix denotes the induced matrix norm. We denote the scalar-product of two vectors by . We write sets in calligraphic and denote by its interior, by its closure, by its boundary, and by its complement. We define and denote by the set of the natural numbers larger than zero.
A function is said to be of class if it is continuous, monotonically increasing, , and . A function is said to be of class if it is continuous, monotonically decreasing, and . A function is said to be of class if is of class for every fixed and is of class for every fixed .
We write and for probability and expectation, respectively. The underlying probability spaces should always be clear from the context. Conditional probabilities and expectations are denoted by and , respectively. We denote the characteristic function of a set by . The abbreviation a.s. stands for almost surely, that is, with probability one, and means equal a.s. By and , we denote the set of all -valued and -valued adapted processes such that a.s. for all . Finally, .
2. Stochastic Differential Equations and Stability
Stochastic differential equations (SDEs) often come up when modeling a physical system described by an ODE and perturbing this system with noise, corresponding to either uncertainty, measurement error, or unknown complications in the system. The stochastic integral as presented by Itô is in some sense the natural way to approach the situation when the underlying deterministic system responds causally to the noise. This is the case for example when describing financial markets when it is instrumental that the noise has no autocorrelation. In cases when the noise is an effective model of some unknown and/or complicated dynamic subsystem, then it can be more natural to represent this via the so-called Stratonovich stochastic integral. In this paper, we shall restrict our attention to the Itô stochastic integral since a Stratonovich approach is easily represented by a slightly modified Itô system. We study the stability of the origin of the SDE of Itô typewhere and are (globally) Lipschitz continuous functions; that is, there exists a such thatand such that and .
We consider strong solutions to (3), that is, given a filtered probability space satisfying the usual conditions ( complete, right continuous, containing all null sets), a -valued Brownian motion defined on , and an measurable initial distribution fulfilling , an -valued stochastic process defined on is called a (unique) solution to (3) if it has the following properties (cf., e.g., [22, Section 2.3], [23, Section 21]):
Existence(i) is continuous a.s. and -adapted.(ii)The processes and belong to and , respectively.(iii)For every , the equation holds a.s. The second integral is interpreted in the Itô sense.
Uniqueness. Any solution to (3) is indistinguishable from ; that is, for every , we have .
For deterministic initial solutions, that is, a.s., we write for the solution.
The following theorem provides the standard solution theory, the backdrop to our study of stability later.
Theorem 1. The SDE (3) with initial distribution has a unique solution . Further, is a strong Markov process; that is, for every bounded Borel-measurable function , any -stopping-time a.s., and any , we have
For a proof of the last theorem, compare, for example, [22, Section and Theorem ]. Note that since we have and the Lipschitz condition implies the so-called linear growth condition.
A plethora of concepts are in use concerning the stability of SDEs. Here we will be concerned with the so-called asymptotic stability in probability of the zero solution [24, ], also referred to as stochastic asymptotic stability [22, Definition ]. For a more detailed discussion of the stability of SDEs, see the books by Khasminskii [24] or Mao [22]. For completeness, we recall the definitions of stability in probability and asymptotic stability in probability.
Definition 2 (stability in probability). The null solution to the SDE (3) is said to be stable in probability (SiP) if for every and every there exists a such that
Definition 3 (asymptotic stability in probability (ASiP)). The null solution to the SDE (3) is said to be asymptotically stable in probability (ASiP) if it is SiP and in addition for every there exists a such that
Our definitions of SiP and ASiP are equivalent to the more common The reason for our formulation is that we want to look at a more practical concept related to such stability, namely, stochastic analog of the BOA in the stability theory for deterministic systems. We therefore back from the limit and consider the following instead: Given some confidence how far from the origin the sample paths can start and still approach the equilibrium as with probability greater than or equal to . This is the motivation for the next definition.
Definition 4 (-basin of attraction (-BOA)). Consider system (3) and let . One refers to the set as the -basin of attraction or short -BOA.
Note that a 1-BOA corresponds to the usual BOA for deterministic systems.
For the SDE (3), the associated generator is given byfor some appropriately differentiable with . Notice that this is just the drift term in the expression for the stochastic differential of the process ; compare (1). As discussed in the Introduction, the generator for a stochastic system corresponds to the orbital derivative of a deterministic system.
Definition 5 (local Lyapunov function). Consider system (3). A function , where is a domain, is called a (local) Lyapunov function for system (3) if there are functions , such that fulfills the properties: (i) for all .(ii) for all .
Remark 6. It is of vital importance that is not necessarily differentiable at the equilibrium, because otherwise a large number of systems with an ASiP null solution do not possess a Lyapunov function; compare [24, Remark 5.5].
The following theorem provides the first centerpiece of Lyapunov stability theory for our application; compare [24, Theorem 5.5 and Corollary 5.1].
Theorem 7. If there exists a local Lyapunov function as in Definition 5 for system (3), then the null solution is ASiP. Further, let and assume that is a compact subset of . Then, for every , the set is a subset of the -BOA of the origin.
Proof. That the null solution is ASiP follows directly by [24, Theorem 5.5 and Corollary 5.1].
To prove the second assertion fix and an arbitrary . We must show that Let be the first exit-time of the process from . Then is a stopping-time and using identical argumentation as in the proof of [24, Theorem 5.5], we note that a.s. on we have . It thus suffices to show that .
By [24, Lemma 5.4] the stochastic process is a supermartingale; compare also Remark 8. Therefore for any . Since , we can take limits to get . Hence which proves the theorem.
Remark 8. If the condition is violated, even at only one point , one cannot conclude that is a supermartingale; compare [24, p. 149]. The reason why does not have to hold true at the equilibrium at the origin is because it is so-called inaccessible point of the SDE, compare again [24, p. 149], meaning that
Instead of trying to compute a local Lyapunov function for system (3) on a large domain , we suggest to compute it on a small , for example, by using linearization, and then to compute a nonlocal Lyapunov function, defined below, to compute a reasonably sized rigid estimate on the equilibrium’s -BOA. The advantages of this procedure were explained in the Introduction and are addressed again in Section 4, where we give a roadmap of our approach.
Definition 9 (nonlocal Lyapunov function). Let , , be simply connected compact neighbourhoods of the origin with boundaries and set . A function for system (3) such that (1) for all , , , and either(2a) for all or(2b) and is positive definite for all , is called a nonlocal Lyapunov function for system (3). We refer to as the outer boundary of and as the inner boundary of .
Lemma 10. Assume system (3) has a nonlocal Lyapunov function as in Definition 9 and let . Let be the exit-time from . Then
Proof. We first assume that fulfills (2a) in Definition 9 and set and note that since the solution process remains in for all , we get by [24, Lemma 3.2] and the nonnegativity of that Now and then Taking the limit delivers and then .
If the condition (2b) in Definition 9 is fulfilled (instead of condition (2a)), then, by [20, Theorem 6.14], it follows that the boundary value problem has a solution . Just as above with substituted for , it follows that .
For the second proposition of the lemma define the stopping-times and and observe that a.s. because is continuous a.s. Now and then that is, we have because , which completes the proof.
We now show how to combine a local Lyapunov function and a nonlocal one to get a rigorous lower bound on -BOA. For the schematic representation of the interrelation between the (sub)level sets in Theorem 11 see Figure 1.
Theorem 11. Consider system (3) and assume there exist a local Lyapunov function as in Definition 5 and a nonlocal Lyapunov function as in Definition 9 for the system. Assume further that the constants and from Theorem 7 and the constants and the set from Lemma 10 are such that Then the set is a subset of the -BOA of the origin, where .
Proof. Let be arbitrary but fixed. We must show that Define the stopping-times , , and . In the following, we suppress the argument of stopping-times (as usual) if it can be seen from the context. Define for every the eventNote that since is closed and the sample paths are continuous a.s., we have .
Let us first assume ; we deal with the case below. Then Because the events and are disjoint we haveThe strong Markov property of (cf. Theorem 1) delivers (cf., e.g., [25, Theorem ] or [24, Theorem ]) Define Now , a Borel set with respect to the relative topology of , defines a probability measure on and, for example, by [23, Lemma 8.7], we getFurther, because by Lemma 10, we have and we get from (28) the estimate which together with and the fact that impliesSimilarly, we get for every thatNow let be arbitrary and define the stopping-time . Then The assumption delivers by Lemma 10 that . By using the assumption we can show just as above that Thusand especiallyPlugging this estimate for into (33) deliversfor every .
That estimate (39) also holds true if can be seen from the fact that estimate (37) is clearly valid for any . Combining (37) and (38) and noting that we get We finish the proof by showing that is equal to a.s. Note first that cannot stay in for an infinitely long contiguous time interval a.s. by Lemma 10. Second, as in the proof of Theorem 7, we have a.s. if stays in for an infinitely long contiguous time interval. Thus implies the following event a.s.: There exists a sequence of increasing times , , such that and for all . Because and and by what we have showed above the probability of this event is no larger than , which completes the proof.
Remark 12. Note that in the proof of Theorem 11 the conditions on and in (3) need only to hold within , because we stop our solutions at the boundary , and therefore and can be altered outside of without changing the conclusions of the theorem. It thus suffices that and fulfill the Lipschitz condition (4) on , which is, for example, always the case if and are locally Lipschitz.
3. Worked Out Example
We consider an example to show how our approach works. Because we do not want to consider appropriate methods to solve the relevant PDEs in this paper we consider a 1-dimensional example, but a highly nonlinear one. The system is given by the SDEso that our coefficient functions in (3) are given by and and we have only a single Wiener process, denoted by because we reserve for the local Lyapunov function, driving the noise.
We first compute a local Lyapunov function for (41) as in Definition 5 using the ansatz with . Direct calculations with reveal so We fix , , and . This corresponds to the deterministic system with an unstable null solution. The noise, however, stabilizes the null solution, which can be seen from for all because then Thus with . To compute a rigid lower estimate on the -BOA of the equilibrium using we can, by Theorem 7 and because is an even function, solve . Thus . In the following we use which delivers .
To enlarge our estimate we solve the PDE (note that for ) with the boundary conditions and with ; compare Lemma 10 and Theorem 11. Because the system is 1-dimensional the PDE is also an ODE and can be conveniently solved with the shooting method. With the initial conditions and , we get using a very conservative step-size of . This corresponds to and . From the results (cf. Figure 2), we can estimate and . This corresponds to and . Note that by extending canonically for negative , we have and . Thus Theorem 11 delivers that is affirmed to be contained in the -BOA of the equilibrium at the origin with For comparison we get the rigid bound , that is, on the -BOA, if we only use the local Lyapunov function . Our combined method thus delivers a considerably larger set that is affirmed to be contained in the equilibrium’s -BOA, namely, instead of from the local Lyapunov function. For different the results are comparable. In Figure 2 we show nonlocal Lyapunov functions computed for the system, not only with but also for to give a better idea of what is happening.
4. Conclusions and Further Work
In this paper, we have worked out theory in Theorem 11 for a method to give rigid estimates on the -basin of attraction (-BOA) of the null solution as defined in Definition 4 for stochastic systems given by SDEs (3). We showed how our method works for a highly nonlinear 1-dimensional system in Section 3. For higher dimensional systems, one has to consider appropriate methods for our setting to solve PDEs.
Our roadmap for further advancing the method is as follows:(1)Given system (3), we start by finding a local Lyapunov function as in Definition 5 and an as large as possible, such that is a compact subset of . One method for doing this involves linearizing the system (i.e., the functions and ) via Taylor expansion around . For a wide class of systems, we expect the linearized system to admit locally a Lyapunov function of the form with and a positive definite matrix; compare [24, Theorem 5.12]. We have worked with our collaborators on a method using sum-of-squares programming for computing such Lyapunov functions for linear systems and the initial results are promising [11].(2)Next we determine a nonlocal Lyapunov function with bounded but as large as practical. This must be done in such a way as to ensure and for appropriate , ; compare Theorem 11. From these Lyapunov functions we know that for every , , the set is a subset of the -BOA of the null solution, where . We propose to address this by solving the following PDE boundary value problem using an appropriate numerical method: If the PDE is strictly elliptic, that is, the matrix is positive definite for all , then the PDE possesses a unique solution; compare, for example, [20, Theorem 6.14]. In this case, we can set and from a rigid error estimate of the numerical approximate solution to the unique solution, we can compute rigid sublevel sets of and thus give rigid lower bounds on -BOA. If the PDE is not strictly elliptic, then we suggest fixing . In this case, we must be sure that the solution we compute numerically really fulfills for all and that its sublevel sets are compatible with Theorem 11. That is, the computed function must be a true nonlocal Lyapunov function for the system and not just an approximation!
We are optimistic that our suggested method to compute rigid estimates to -BOA can be developed into a useful tool for the study of stochastic systems, both in practice and in theory.
Conflicts of Interest
There are no conflicts of interest related to this paper.
Acknowledgments
This work was supported by The Icelandic Research Fund, Grant no. 152429-051.