In this chapter, we consider general autonomous nonlinear systems formalized as a set of ODEs
˙
x(t) =f(x(t)), ∀x(t)∈D⊆Rn, (4.1) wherexis the vector of state variables,D is the domain andf :D7→Rn is locally Lipschitz continuous.
The time dependency is often omitted for simplicity of notation.
4.2.1 Problem Formulation and Objective
Suppose the system under investigation has a position of rest denoted byxs∈D, i.e. f(xs) =0. First we will establish stability ofxsfor the system (4.1) in the sense of Lyapunov. In fact, there are different forms of stability such as exponential, uniform or bounded-input bounded-output (BIBO), see [76]; in this chapter, however, we particularly consider the so-called Lyapunov asymptotic stability, since it is directly associated with the estimation of the ROA. Lyapunov asymptotic stability establishes whether the initial state x(0) of (4.1) will always remain confined within a small neighbour surrounding xs, if they are located near the equilibrium point. Formally, asymptotic stability is defined precisely as:
Definition 4.1. Lyapunov asymptotic stability: The equilibrium pointxs∈D of (4.1) is said to be Lyapunov stable att= 0 if, and only if, for any>0there exists aδ>0such that
|x(0)−xs|<δ ⇒ |χ(x(0);t)−xs|<,∀t≥0,
and asymptotical stability is established if xs is stable and locally attractive; that is the region δ have to be chosen such that
|x(0)−xs|<δ ⇒ lim
t→∞χ(x(0);t) =xs.
Here δ and correspond to two arbitrary regions in the state-space, and χ(x(0);t) specifies a system trajectory of (4.1) starting from the initial conditionx(0).
Assuming that the equilibrium pointxsof (4.1) is stable according to Def. 4.1, then one can clearly see thatδdefines a region in the domainD from which any initial state is guaranteed to be attracted byxs. In this chapter we seek to find this stability region, which is formally defined as follows:
Definition 4.2. Region of attraction: Given the nonlinear system (4.1) whose equilibrium point xs is locally asymptotically stable in the sense of Lyapunov according to Def. 4.1, then the ROA surroundingxsis expressed via the set
Se(xs) =n
x0∈D : lim
t→∞χ(x(0);t) =xso .
Here the subscript e specifies the exact stability region which is hard, or even impossible, to express analytically except for a specific class of nonlinear systems; hence, similarly to computation of reachable sets, our objective is to find an estimate of this region such that S(xs)⊆Se(xs) withS(xs) as the set associated with the under-approximation of the exact ROA.
4.2.2 Existing Techniques
As stated earlier and illustrated in Fig. 4.1, the existing methods to estimate the ROA can be broadly cat-egorized into Lyapunov-based and Lyapunov-free techniques. Here we shall briefly illustrate a technique from each category in order to compare them with our proposed contribution later in Sec. 4.4.
4.2.2.1 Level Set Method
The first class of techniques is based on an Eulerian scheme using level set methods. Generally speaking, LSMs are a collection of numerical algorithms solving a specific class of PDEs often encountered in simulation of dynamic implicit surfaces for a variety of applications, such as fluids, image processing, and computer vision [97]. Instead of explicitly representing a surface via its vertices, edges or faces, the LSM describes surfaces implicity via a level set function, often denoted byφ(t,x) :Rn7→R. Using LSMs, one can compute the backward reachable set starting from a small neighboured surroundingxs.
Here finding an estimate of the ROA is reformulated into a reachability problem, which in turn is casted into solving the following first-order time-dependent HJI PDE:
∂φ(t,x)
where H is a basic Hamiltonian function. For further details about Hamiltonian dynamics, HJI formu-lation and viscosity solutions, the reader is referred to [21, Ch. 2]. It is proven in [99] that the viscosity solution of (4.2) provides a surface representationφ(t,x)≤0 corresponding to the continuous reachable set of differential state variables; that is the backward reachable set is expressed via:
BackwardReach(R(0), t) :=n
x∈D : φ(t,x)≤0o
=Se(xs),
withBackwardReachdenoting the backward reachable set starting from the initial set of statesR(0).
4.2.2.2 Lyapunov method
The second class of techniques is based on Lyapunov direct method. Without loss of generality, we place the equilibrium pointxsat the origin 0n. Let V(x) :D7→Rdenote a scalar positive definite function, that isV(x)>0, andV(0n) = 0. This function is often refereed to as the so-called Lyapunov function of (4.1), where its partial derivatives are continuous functions within the domain D surrounding xs. Furthermore, the derivative with respect to time has to be negatively definite inD, that is:
V˙ =
Here the problem of providing an estimate of the ROA is reformulated into finding an optimal Lyapunov function, whose sub-level set provides the largest estimate of the ROA. Formally, the aforementioned task is equivalent to solving the following optimization problem [59]:
sup
with the constantcv∈R+ andρ(·) specifying a pre-definable measure of the sub-level set V(cv).
Example of reachable set not attracted by the target set Examples of reachable sets attracted by the target set Recursive partitioning of the grid (black boxes)
Enlargement of the target set (blue boxes) Stablity region (union of light gray boxes) Initial reachable set selected for chosen reachable set examples (dark gray boxes)
Figure 4.2: Step-by-step computation of the stability region using forward reachable sets. The blue boxes represent the iterative enlargement of the target setTg constructed around an equilibrium point (step 1).
The boxes show the recursive partitioning of the grid to provide an accurate estimate of the ROA (step 2), see remark4.2. The dark gray areas are the cells whose reachable set is attracted by the enlarged target set, whereas the white areas correspond to the cells whose reachable set does not converge. Three random cells were chosen to illustrate the evolution of the reachable set (dark red), where the computation stops (highlighted by the red area) when the reachable set of thei-th cellRi(t), t≥0 is either a subset of the target setRi(t)⊆Tg or if the Lagrangian remainder is not a subset of the maximum linearization errors L*Lmax, see (4.8).
4.2.3 Proposed Approach
The problem is approached in this chapter using forward reachable sets based on the algorithms presented earlier in Sec. 2.4. Hereafter, we introduce the target set denoted by Tg. Basically, the target set is a small region surrounding the stable equilibrium point, based on Lyapunov asymptotic stability, see Def. 4.1. With the introduction of the target set, our problem of findingS(xs) is reformulated into
S(xs) :=R1(0) ∪ R2(0)∪ . . . ∪ RNcells(0)⊆Se(xs) (4.4)
s.t.
∀i∈ {1, . . . , Ncells}:
∃t≥0 : Ri(t)Def.4.1⊆ Tg, Ri(t)∈reach(Ri(0)) :=n
x(t)∈D : x(t) satisfies (4.1) within t≥0 forx(0)∈Ri(0)o . In (4.4), the stability region is specified via the union of initial reachable setsRi(0), i∈ {1, . . . , Ncells}, whereNcells ∈N specifies the number of cells resulting from a partitioning of the domainD. The main idea is to discretize the working domain into smaller regions and examine whether each cell belongs to the ROA; that is for eachRi, we check at each time instant whether the reachable set of differential state variables is confined within a target set denoted byTg, see Fig. 4.2. As stated earlier throughout Ch. 2, it is proven that the exact reachable set for the class of nonlinear systems are not computable [111];
thus an over-approximation, which includes all behaviours of (4.1) is performed as tightly as possible, see Fig. 2.7.
4.2.4 Organization
The remainder of the chapter is organized as follows: In Sec. 4.3 we briefly recall reachability computa-tions for the class of systems described via nonlinear ODEs, then we present the estimation algorithm.
In Sec. 4.4we illustrate the applicability of our algorithm and compare it with existing techniques. The chapter is summarized in Sec. 4.5.