Zubov's equation for state-constrained perturbed nonlinear systems

The paper gives a characterization of the uniform robust domain of attraction for a nite non-linear controlled system subject to perturbations and state constraints. We extend the Zubov approach to characterize this domain by means of the value function of a suitable in nite horizon state-constrained control problem which at the same time is a Lyapunov function for the system. We provide associated Hamilton-Jacobi-Bellman equations and prove existence and uniqueness of the solutions of these generalized Zubov equations.


Introduction
The domain of attraction of a locally asymptotically stable fixed point is an important object in the analysis of the long term behavior of dynamical systems. A seminal result was made by V.I. Zubov [21], who proved that for dynamical systems induced by ordinary differential equations the domain of attraction can be represented as a sublevel set of a solution of a first order partial differential equation, the Zubov equation. Beyond characterizing the domain of attraction, the solution to Zubov's equation also yields a Lyapunov function. Zubov's idea has been used in different contexts [1,18,14] and was extended in various ways [3,9,11,6,10,17]. Particularly the latter references have established a connection of Zubov's method to optimal control, as they identify the solutions of (generalized versions of) Zubov's equation as optimal value functions of suitable optimal control problems. This connection is established by the observation that Zubov's equation is a Hamilton-Jacobi-Bellman type equation and is thus naturally linked to an optimal control problem. Typically, the generalized Zubov equations derived in these references have to be treated in the framework of viscosity solutions, since classical solutions are in general not guaranteed to exist.
In many practical applications of ordinary differential equations the range of physically meaningful or desired states is restricted to a subset of the state space. This fact has long since been recognized in optimal control in which the treatment of state constraints is a classical yet still active research topic. A typical feature of state constraints is that they may affect the continuity of the optimal value function as well as structural properties of this function. However, when introduced in an appropriate way, continuity can be maintained as demonstrated by the construction proposed in [5] for a reachability problem and extended to a more general setting of optimal control problems in [2,13]. This approach allows to incorporate state constraints by means of appropriately chosen penalization functions into the Hamilton-Jacobi framework without loosing the continuity of the optimal value function while providing an exact representation of the state constraints.
In this paper we show that the approach of incorporating state constraints from [5,2] can be successfully carried over to the Zubov equation. We will do this for ordinary differential equations subject to deterministic perturbations, a setting first considered in [9], noting that the extension to the controlled [10] or game theoretic setting [17] is straightforward 1 . Following this approach, we are able to show that besides continuity all of the main features of Zubov's method can be maintained under state constraints. More precisely, we prove that for the state constrained Zubov problem the optimal value function • is continuous, cf. Theorem 3.1 and Remark 3.2 • characterizes the domain of attraction as a sub-level set, cf. Theorem 3.1, and is a Lyapunov function for the system, cf. Remark 4.3 • is the unique viscosity solution of an appropriate Hamilton-Jacobi-Bellman equation, i.e., of a generalized Zubov equation, cf. Theorem 4.8.
Moreover, the fact that the continuity of the optimal value function is maintained implies that the approach can be used as the basis for numerical computations, which we illustrate in the last section of this paper. As usual in the Zubov literature, we will formulate our results in parallel for two variants of the Zubov equation which are linked through the Kruzhkov transform, since both variants have their specific advantages.
The paper is organized as follows. In Section 2, we specify the problem and prove useful properties of domains of attractions under state constraints. In Section 3, we provide optimal control characterizations of the uniform domain of attraction for which in Section 4 we define the related generalized Zubov equations and prove existence and uniqueness of the corresponding viscosity solutions. Finally, Section 5 concludes our paper with a numerical example. Throughout the paper, in order to avoid repetitions we provide detailed proofs only for those results whose proofs significantly differ from the unconstrained case.

Setting of the problem
We consider a non-linear system subject to time varying perturbations given by the ordinary differential equationẏ where f : R N × U → R N is the vector field, u is the perturbation input that is supposed to take values in the compact set U ⊂ R m (with m ≥ 1). We assume that the system has a locally robustly asymptotically stable equilibrium (for details on this assumption see (A2), below); without loss of generality we assume that this equilibrium is located in the origin.
We define the set of admissible input functions as and assume that the vector field f in (2.1) satisfies the following regularity assumptions (where · is the Euclidean norm and B(y, R) the ball around y with radius R > 0 in R N ): Note that (ii) may be assumed without loss of generality replacing f by f /(1 + f ) if necessary.
Under assumption (A1), for any u ∈ U and for any x ∈ R N , there exists a unique absolutely continuous trajectory y = y u x (·) satisfying (2.1) for almost all t ∈ R and y(0) = x. For any T > 0, the set of all feasible trajectories on [0, T ] starting in x will be denoted as is absolutely continuous, satisfies (2.1) for some u ∈ U and almost all t ∈ (0, T ), and y(0) = x .
Again under assumption (A1), for any T > 0 and For every t ≥ 0 and x ∈ R N , we denote by R(x, t) the reachable set, at time t, starting from the position x, that is It is clear that under assumption (A1) the set R(x, t) is bounded for every x ∈ R N and every t ≥ 0. Moreover, we will use the set R ∞ (x) := t≥0 R(x, t).
Throughout the paper, we assume that the origin x = 0 is an equilibrium for the system (2.1) and that we have the following local uniform asymptotic stability property f (0, 0) = 0; (ii ) There exists a KL-function 2 β andr > 0 such that Note that (ii) implies the existence ofε > 0 such that Our goal in this paper is to characterize the robust domain of attraction of the origin under state constraints. Here, "robust" is to be understood in the sense "for all perturbations u ∈ U". In order to formally define this object, we consider an open set Ω ad ⊂ R N of admissible states, whose complement Ω c ad = R N \ Ω ad defines the state constraints or obstacles we impose on the system. We assume that the state constraints do not touch the equilibrium x = 0, i.e., that there existsr > 0 such that B(0,r) ⊂ Ω ad . We may without loss of generality assume that thisr > 0 is the same as in (A2).
Using this admissible set, we define the set of admissible input functions as and the robust domain of attraction of x = 0 w.r.t. Ω ad as D := x ∈ R N U ad (x) = U and y u x (t) → 0 for all u ∈ U .
In words, D contains all those initial values x for which all corresponding solutions y u x (·) stay in Ω ad and converge to 0.
In order to relate the domain of attraction to a Zubov type PDE, similar to [9] we need to consider a uniform version of D. In extension to the "uniform attraction rate" imposed in [9], here we also need to consider a "uniform distance" to the obstacle set Ω c ad . To this end, we define the distance between a point x ∈ R N and a set A ⊂ R N by dist(x, A) := inf y∈A x − y . Then, for δ ≥ 0, we define the set of δ-admissible input functions as , Ω c ad ) > δ for all t ≥ 0}. Note that U ad,0 (x) = U ad (x). The robust domain of uniform attraction of x = 0 w.r.t. Ω ad is then defined by The set D 0 is a uniform version of D in the sense that for each x ∈ D 0 the trajectories have a positive distance δ to Ω c ad and a speed of convergence β(t) towards 0 which both do not depend on u ∈ U. The relation between D 0 and D will be analysed below.
An important object for our analysis is the first hitting time with respect to B(0,ε), that is the first time a given trajectory reaches the ball B(0,ε) withε from (2.2): This observation in particular implies B(0,ε) ⊂ D 0 , with δ =r/2 and β(t) = β(t,ε).
The following lemma establishes further links between the sets D 0 , D and the hitting time t(x, u).
Proof: (i) Let x ∈ D 0 , pick the corresponding δ > 0 and β(t) and let T > 0 be such that Hence, t(x, u) < T for all u ∈ U and thus sup u∈U t(x, u) ≤ T < ∞. Conversely, let x ∈ R N be such that T := sup u∈U t(x, u) < ∞ and U ad,δ (x) = U for some δ > 0. Then, for t ≥ T we have y u x (t) ≤ β(r, t). Hence, for t ≥ T we can choose β(t) := β(r, t). For t ∈ [0, T ], since the vector field is bounded, the solutions are uniformly bounded, say, by a constant M > 0. Setting β(t) := M for t ∈ [0, T ] then yields the function β(t) with the desired properties.
(ii) Let x ∈ D 0 and consider the corresponding δ > 0 and β(·). Let T > 0 be such that β(t) <ε/2 for all t ≥ T . By continuity of the trajectories in the initial value x there exists a neighbourhood B(x, ε) of x such that y u To see the converse inclusion int D ⊂ D 0 , let x ∈ (int D) \ D 0 . Since x ∈ D 0 , by (i) either sup u∈U t(x, u) = ∞ or U ad,δ = U for all δ > 0 must hold. If the former holds, then as in [9, Proof of Proposition 2.3(iv)] we obtain x ∈ ∂D, contradicting x ∈ int D.
Hence assume T := sup u∈U t(x, u) < ∞. Then U ad,δ = U for all δ > 0 must hold and thus we can conclude the existence of u n ∈ U and t n > 0 with dist(y un x (t n ), Ω c ad ) → 0 as n → ∞. Since (2.3) and t(x, u n ) ≤ T imply t n ≤ T and f is bounded by (A1)(ii), both t n and y n := y un x (t n ) must be bounded. Hence, the Lipschitz assumption (A1)(i) yields that for any ε > 0 the set {y un ξ (t n ) | ξ ∈ B(x, ε)} contains a ball B(y n , ρ) with ρ > 0 independent of n. For sufficiently large n this implies B(y n , ρ) ⊂ Ω ad . This means that u n ∈ U ad,0 (x n ) for some x n ∈ B(x, ε) and consequently x n ∈ D. Since ε > 0 was arbitrary, Remark 2.2 In words, (iii) implies that D 0 and D coincide except for sets with void interior. Moreover, D 0 is the largest subset of D which can be expressed as a strict sub- which is the core of Zubov's method. Indeed, the fact that Zubov's method naturally leads to a characterization of an open set is also the reason for choosing the set of admissible states Ω ad as an open set and not -as it is more common in the control literature -as a closed set.
3 Characterization of the set D 0 As an important step towards characterizing D 0 via a Zubov type equation, in this section we introduce an optimal control problem whose optimal value function characterizes D 0 . The idea used here has been developed in many papers for the case of differential systems without state constraints, see, e.g., [6,9,10,11]. The main novelty here is the extension to the case when the trajectories are constrained to an open set Ω ad . To define the adequate optimal control problem, we consider a running cost g : R N × U → R satisfying: (A3) (i) g is continuous and for each R > 0 there exists L g,R ≥ 0 such that: (ii) g(0, u) = 0 for any u ∈ U , and for every c > 0 we have: Existence of such a function follows, e.g., from the construction in [10, Section 3]. Note that the converse implication to (iii), i.e., that the integral value is Like in the other Zubov papers cited above, the cost g measures the convergence of the solutions to 0. In order to incorporate state constraints into our setting, we need another cost that indicates when a given trajectory violates the constraints in Ω ad . For this, we consider a function h : R N → [0, +∞] satisfying: The function h can be chosen for instance as h( . This general expression always satisfies the requirements of (A4), however, for particular domains Ω ad , simpler expressions may exist and could be preferred. In any case, the exact expression of the function h is not important as long as h satisfies (A4). The same remark holds for the function g satisfying (A3).
Consider also the Kruzhkov transformed optimal value function v : The following theorem shows the relation between V , v and D 0 .
Theorem 3.1 Assume (A1)-(A2) and let g and h satisfy (A3) and (A4), respectively. Then, Since y u x (t) ≤ β(t) for any u ∈ U, the reachable set R ∞ (x) is bounded, hence R ∞ (x) is compact. Moreover, since U = U ad,δ (x) for some δ > 0, R ∞ (x) ⊂ Ω ad follows and since h is continuous on Ω ad it will attain a (finite) maximum on R ∞ (x) which shows the claim. Let x ∈ D 0 . Then either the existence of β(t) or the existence of δ in the definition of D 0 is not satisfied. In the first case, we can proceed as in [9, Proof of Proposition 3.1(i)] in order to conclude V (x) = ∞. In the second case, the existence of β again implies that R ∞ (x) is bounded. Moreover, the non-existence of δ implies the existence of controls u n and times t n with dist(y un x (t n ), Ω c ad ) → 0. Since x n := y un x (t n ) lies in the bounded set R ∞ (x), we can find a subsequence x n k converging to some x 0 ∈ Ω ad . This implies h(x n k ) → ∞ and since V (x) ≥ h(x n k ) we obtain V (x) = ∞.
(ii) In order to prove continuity, let x, y ∈ D 0 . Then we obtain For (3.6), we can follow the proof of [9, Proposition 3.1(iii)] in order to see that this term vanishes as x 2 → x 1 . In order to show the same for (3.6), let u ∈ U be arbitrary. Then, as x 1 ∈ D 0 , the convergence h(y u x 1 (t)) → 0 holds as t → ∞, and due to the uniformity of this convergence in u, the continuous dependence of the trajectories on the initial value and the local asymptotic stability, this convergence is uniform for all u ∈ U and all initial values near x 1 . This means that we can find a neighborhood B(x 1 , δ) and a function γ(t) → 0 such that |h(x u x 2 (t))| ≤ γ(t) holds for all x 2 ∈ B(x 1 , δ). This implies that the supremum over t in (3.7) is attained on a finite interval [0, T ]. On a compact time interval, however, the map x → h(y u x (t)) is continuous in x uniformly in u and t, implying that (3.7) also tends to 0 as x 2 → x 1 . This shows the desired continuity.
For the second assertion, first consider x n → x ∈ D 0 . Since x ∈ D 0 , either the existence of β(t) or the existence of δ must fail for this point. In the first case, V (x n ) → ∞ follows by the same arguments as in the proof of [9, Proposition 3.1(iv)]. In the second case, as in the proof of Lemma 2.1(iii) we find a bounded sequence of times t n ∈ [0, ∞) and a sequence of functions u n ∈ U with dist(y un x (t n ), Ω c ad ) → 0 as n → ∞. Since t n is bounded and x n → x we can conclude dist(y un xn (t n ), Ω c ad ) → 0 and thus V (x n ) ≥ h(y un xn (t n )) → ∞. In case x n → ∞ the assertion follows similar to [9, Proof of Proposition 3.1(iv)].
(iii) Follows immediately from (i) and (ii) using the formula v(x) = 1 − e −V (x) . Remark 3.2 Both the original V and the transformed v have advantages and disadvantages. The original function V has an easier definition and is thus easier to handle in proofs. However, it is unbounded and only defined on the subset D 0 which makes it inconvenient for constructing numerical approximations and turns the determination of D 0 via V into a free boundary value problem. In contrast to this, the function v is defined, bounded and continuous on the whole R N . Thus, it is much better suited for numerical computations. However, this advantage comes at the expense of a more complicated definition.

Zubov's equation
is a viscosity solution of (4.8) if it is sub-and supersolution.
The first step to proving that V and v satisfy Zubov-type equations is to establish suitable dynamic programming principles for these functions. The arguments for establishing these principles are standard (see, e.g., [4]) and the corresponding proofs are thus omitted.
For the formulation of these principles, we define the set Since

Moreover, for any open set
Then, the following dynamic programming principle holds:

Proposition 4.2 Abbreviate
G(x, t, u) := t 0 g(y u x (τ ), u(τ )) dτ (4.9) Then, under (A1), the following assertions are satisfied: (i) For all x ∈ D 0 and all T ≥ 0 we have: (4.10) (ii) For all x ∈ R N and all T ≥ 0, we have: . (4.13) From (4.10), standard viscosity solution arguments yield that for all x ∈ Ω 0 the function V satisfies the Hamilton-Jacobi-Bellman equation in the viscosity sense. Likewise, (4.12) implies that the transformed function for all x ∈ Ω 0 . Moreover, for all x ∈ D 0 and all t ≥ 0, the dynamic programming principle (4.11) implies which yields that V is a viscosity super-solution of (4.14) on the whole set D 0 . Similarly, (4.13) implies that v is a viscosity super-solution of (4.15) on the whole R N .

Remark 4.3 Inequality (4.16) implies
for all u ∈ U and all x = 0. This implies that the function V is a Lyapunov function for the system for all perturbation inputs u ∈ U. Thus, the modifications we have made in order to incorporate the obstacles Ω c ad have not affected the Lyapunov function property of Zubov's original approach.
The following theorem introduces the two central equations of this paper, i.e., the two variants of the Zubov equation with state constraints. Its proof follows again by standard viscosity solution arguments from the respective dynamic programming principles.
Likewise, the function v is viscosity solution of the equation While the fact that V and v solve (4.17) and (4.18), respectively, follows by standard viscosity solution techniques, establishing uniqueness of the solutions to these equations requires more arguments which we provide next.
In order to incorporate the state constraints into our analysis, we start with a localized version of [9, Proposition 3.5]. Proposition 4.5 (i) Let W be a l.s.c. super-solution of (4.14) on an open set O ⊂ R N . Then W satisfies a super-optimality principle, that is for any x ∈ O, we have (ii) Let U be a u.s.c. sub-solution of (4.14) on an open set O ⊂ R N and let U : O → R be a continuous function with U ≤ U . Then U satisfies a sub-optimality principle, that is for any x ∈ O we have .
(iv) Letũ be a u.s.c. sub-solution of (4.15) on an open set O ⊂ R N . Thenũ satisfies a sub-optimality principle, that is for any x ∈ O, we havẽ with G from (4.9).
The proof of this proposition is omitted as it follows along the same lines as similar results in the literature: The link between the notion of sub-solution (resp. super-solution) and sub-optimality (super-optimality) has been analysed in several papers, we refer to [19,20] or [4, Chapter III] for a proof using viscosity techniques. Another proof based on non smooth analysis is given in [12].
The main consequence of the above proposition are the following comparison principle for the equations (4.17) and (4.18). If, moreover, W (x n ) → ∞ holds whenever x n → x ∈ ∂Ω, then Ω = D 0 follows.
Proof: (i) We first show the inequality W (x) ≥ V (x) on D 0 . Let x ∈ D 0 , fix some ε > 0 and pick u ε and t ε such that the supremum on the right hand side of (3.4) is attained up to this ε. This implies x (τ ), u ε (τ ))dτ + h(y uε x (t ε )) + ε.
Since W is a super-solution of (4.17) on Ω, it is also a super-solution of (4.14) on Ω. Hence, (4.20) holds with O = Ω. Since any solution starting in D 0 stays in D 0 for all future times, we can use (4.19) in order to obtain the inequality Combining the two inequalities and using that W being a super-solution of (4.17) implies which shows the claim since ε > 0 was arbitrary.
In order to show the identity Ω = D 0 we again use that W is a super-solution of (4.14) on Ω and thus satisfies (4.19) with O = Ω. We show that this property implies that for each x ∈ D 0 the value W (x) cannot be finite. To this end, note that for every u ∈ U holds for all t ∈ [0, t 0 (x, u, Ω)]. Now, we distinguish three cases.
(a) t 0 (x, u, Ω) = ∞ and y u x (t) ∈ D 0 for all t ≥ 0. In this case the integral grows unboundedly and W (x) = ∞ follows.
(b) t 0 (x, u, Ω) = ∞ and y u x (t) ∈ D 0 for some t ≥ 0. In this case by continuity of the solution there exist times t n with x n := y u x (t n ) ∈ D 0 and dist(x n , ∂D 0 ) → 0 as n → ∞. Then, the inequality W ≥ V on D 0 and the fact that V (x n ) → ∞ for dist(x n , ∂D 0 ) → 0 implies W (x n ) → ∞ for n → ∞ and hence again W (x) = ∞ follows.
Summarizing, in all cases we obtain W (x) = ∞, which shows that W cannot attain finite values outside D 0 and thus Ω = D 0 must hold.
(ii) Consider the set which is open since both U and h are continuous. For x ∈ Ω 1 we have which proves the desired inequality. In order to prove the inequality for x ∈ Ω 1 , observe that since U is a sub-solution of (4.17), it is a sub-solution of (4.14) on Ω 1 , hence (4.20) holds with U = U on O = Ω 1 . Hence, for each ε > 0 we can find a control u ε satisfying U (x) ≤ t 0 g(y uε x (τ ), u ε (τ ))dτ + U (y uε x (t)) + ε for all t ∈ [0, t 0 ] with t 0 := t 0 (x, u ε , Ω 1 ). Now, from (4.16) we get for all t ≥ 0. Hence, for all t ∈ [0, t 0 ] we get Now, in case t 0 is finite, by definition of Ω 1 we get In case t 0 = ∞, since x ∈ D 0 we get y uε x (t) → 0, when t tends to +∞, and thus by continuity of U and V for each sufficiently large t we get Hence, in both cases we get U (x) − V (x) ≤ 2ε which shows the claim.
The next proposition yields the analogous comparison result for the scaled optimal value function v. (ii) Let u : Ω → R, D 0 ⊆ Ω, be a continuous sub-solution of (4.18) with u(0) ≤ 0. Then u(x) ≤ v(x) holds for all x ∈ R N with v from (3.5).
As a direct consequence of Propositions 4.6 and 4.7, we get the following uniqueness result. Proof. (i) In order to prove the inclusion D 0 ⊆ Ω we proceed by contradiction. Assuming that this inclusion does not hold, we can find x ∈ ∂Ω ∩ D 0 . Let x n → x be a sequence in Ω. Since W is a viscosity solution of (4.17) on D 0 ∩ Ω, by Proposition 4.6(ii) we obtain The assertion (ii) follows immediately from Proposition 4.6, while (iii) is an immediate consequence of Proposition 4.7.

Numerical simulations
Theorem 4.8 in conjunction with Theorem 3.1 shows that by computing a solution to the Hamilton-Jacobi equations (4.17) or (4.18) we can compute the uniform domain of attraction D 0 . Due to the non-linearity of the state system, it is in general hopeless to get an analytical expression of the value function V or of v. However, the equations can be approximately solved numerically, at least in low space dimensions. To this end, it is advantageous to consider an approximation of the Kruzhkov transformed optimal value function v via (4.18). This is because v is bounded and because the Zubov equation that characterizes v can be solved on the whole space R N and thus in particular on a large bounded closed domain K that contains D 0 .
Consider a simple grid G with edges x i and a uniform time step h > 0. For simplicity, we assume that 0 is a node of the grid G. An approximation of v can be obtained by using a first order Semi-Lagrangian scheme (see [4, Appendix A]) leading to the problem of solving the following equation.
where v h is piecewise affine on each simplex of the grid G, and satisfies v h (0) = 0. The Semi-Lagrangian scheme is taken here just an example of possible approximation scheme. Other numerical schemes can also be considered like finite difference, Markov Chain methods and in order to improve the accuracy of the method the approximations hg(x i , u) (rectangle rule) and x i + hf (x i , u) (forward Euler approximation) can be replaced by more accurate quadrature rules and ODE solvers, respectively. In the numerical example, below, we used a rectangular quadrature rule and the forward Euler with step size h/50 for this purpose.
Unfortunately, (5.1) has a singularity in 0 since g(0, u) vanishes. This fact was already pointed out in [8] for the case without obstacle. As a consequence of this singularity, the usual fixed-point arguments do not apply and the convergence is not guaranteed. Following the ideas introduced in [8], we use a regularization of (4.18) by introducing an approximation of g: where ε > 0 is a small parameter. Then the regularized Zubov equation that is considered is (5.2) A good feature of this regularization is that v ε converges to v when ε tends to 0. An even better feature is that for all sufficiently small ε > 0 the sub-level characterization D 0 = {x ∈ R N | v ε (x) < 1} still holds; hence in terms of computing D 0 the regularization (5.2) does not introduce any errors. The arguments used in [8] for proving this property can be carried over in a straightforward way to the setting in this paper. The discretization of (5.2) leads to the equation: with again v ε h being a piecewise affine function on each simplex of the grid G and v ε h (0) = 0. Now a straightforward application of fixed-point arguments allow to conclude the existence and uniqueness of the solution v ε h of (5.3). In order to illustrate our findings we provide a numerical example for a system considered in [16, Section C.2] whose unperturbed version was presented in [14]. Its dynamics is given by the two-dimensional systeṁ In order to meet the global boundedness requirement (A1)(ii), as in [16] we perform the transformationf (x, u) = 5f (x, u) which does not change the domain of attraction. Using the perturbation range u ∈ [−0.03, 0.03] the origin is a locally asymptotically stable equilibrium. In fact, it is even exponentially stable, hence the running cost g(x, u) = 10 −3 x 2 satisfies all requirements in (A3). Using this function, we get the numerical approximation for v depicted in In the next computation we introduce a ball with radius r = 2 around the point (−12, 47.5) T as an obstacle. Since we are solving the Kruzhkov transformed Zubov equation The numerical computations were performed using the time step h = 0.2 and piecewise bilinear interpolation on an adaptively generated [15] rectangular finite element grid with about 30000 vertices. The resulting discretized equation (5.3) was solved using the Gauss-Seidel type increasing coordinate iteration described in [15,Section 3], the maximization was performed by discretizing U with 3 equidistant perturbation values.