{'fraction_non_alphanumeric': 0.043056567808487975, 'fraction_numerical': 0.010259572798116276, 'mean_word_length': 4.226486961617345, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatio-Temporal Decomposition of Sum-of-Squares Programs for the Region of Attraction and Reachability

This paper presents a method for calculating Region of Attraction of a target set (not necessarily an equilibrium) for controlled polynomial dynamical systems, using a hierarchy of semidefinite programming problems (SDPs). Our approach builds on previous work and addresses its main issue, the fast-growing memory demands for solving large-scale SDPs. The main idea in this work is in dissecting the original resource-demanding problem into multiple smaller, interconnected, and easier to solve problems. This is achieved by spatio-temporal splitting akin to methods based on partial differential equations. We show that the splitting procedure retains the convergence and outer-approximation guarantees of the previous work, while achieving higher precision in less time and with smaller memory footprint.


I. INTRODUCTION
This paper deals with stability and reachability analysis of nonlinear dynamical systems. A prominent method for quantifying the stability of a nonlinear system is based on calculating its Region Of Attraction (ROA) with respect to a given target set, which is the topic of this work. Using a time-reversal, one obtains the reachable set. Both of these objects are crucial in assessing stability and safety of control system; as a concrete example, let us mention the Flight Envelopes in the aerospace industry that directly coincide with the ROA [1]. A more recent application would be the study of closed-loop parameters variation, of both the controller and the controlled plant, and its effect on the ROA. For example, given a vehicle on the road, these parameters could be the distribution of mass of the vehicle or various road conditions [2].
The dominant methods for computing ROA are based on Lyapunov functions whose level sets provide inner approximations the ROA [3]. For polynomial systems, the Lyapunov function can be found by solving semidefinite programming problems (SDPs) [4] and [5]. These methods only work for autonomous systems and ROA's with respect to given equilibria.
This work expands the approach presented in [6], which is based on optimizing over trajectories of a polynomial dynamical systems modelled as occupation measures and is not limited to an autonomous system nor a stable equilibrium. The contribution proposed in this paper is dissection of the original problem into multiple interconnected problems of lower complexity. The original formulation is discretized in both the time and state variables, aiming to strike a trade-off between the original optimization-based approach [6] where a single polynomial is defined on the entire state-space and various methods based on extensive discretization such as the PDE-based approach [7], set oriented methods [8] or transfer operator approaches [9]. By doing so, our approach, based fully on convex optimization, preserves the outer approximation and convergence guarantees of [6], both of which are explicitly proven in this work, while significantly improving scalability of the method.
Structure of this paper Section II presents the problem statement and Section III introduces the splitting procedure with the proofs of outer approximation and convergence. Section IV states the practical sum-of-squares variant of the problem which is then demonstrated in Section V on numerical examples. The paper ends with a conclusion and future work discussion in sections VI and VII respectively.
Notation The symbol Z j denotes the set of consecutive integers {1, 2, . . . , j}. The Lebesgue measure (i.e., the volume) of a set is denoted by λ(A). The indicator function I A (x) of a set A is the function that takes the value 1 for x ∈ A, and 0 otherwise. The interior of a set A is denoted by A • . The symbols C(A) and C 1 (A) denote respectively the sets of continuous and continuously differentiable functions on A. The set of Borel measurable functions defined on a set A taking values in a set B is denoted by L(A; B) II. PROBLEM STATEMENT Let us consider the nonlinear system with control where x(t) ∈ R n is the state vector, u(t) ∈ R m is the control input vector, t is time, T > 0 is the final time and f is the vector field, which is assumed to be polynomial in variables x and u. The state and control input are constrained by basic semialgebraic sets where g U j (u), g X j (x), and g X T j (x) are polynomials. The region of attraction (ROA) is then defined as where "a.e." stands for "almost everywhere" with respect to the Lebesgue measure.
III. TIME AND STATE SPACE SPLITTING It was shown in [6], that the ROA can be characterized by an infinite-dimensional linear programming (LP) problem in the space of Borel measures or by its dual in the space of continuous functions. These LPs can then be approximated by SDPs, with guaranteed convergence. For space reason, the entire exposition in this work will be in the "dual" setting of continuous functions approximated by polynomials.
Let us restate the original problem from [6]: with variables w(x) ∈ C(X) and v(t, x) ∈ C([0, T ] × X). Any minimizing sequence (w k , v k ) for (4) satisfies w k ≥ I X 0 (x) and w k → I X 0 in L 1 as well as {x ∈ X : v k (0, x) ≥ 0} ⊃ X 0 with convergence in terms of the volume discrepancy tending to zero (see [6] for proof).
Let us now split the state space X into I closed subsets X i and the time interval where K is the number of time splits (meaning K − 1 intervals). It is assumed that The function w(x) will be split into I functions w i (x) and v(t, x) will be split into Assuming that neighbouring subsets X i share boundaries, let us define the set of indices of these neighbours as We are now ready to write the split version of (4) where h a,b is a normal vector of a shared boundary between two neighbouring sets X a and X b , and x a,b is a set of points on said boundary such that x a,b ∈ X a ∩ X b . For simplicity we assume that the normal vector h a,b is independent of x; the case of polynomial or rational dependence of h on x can also be handled [10,Section 4.3]. The optimization variables in (10) are the continuously differentiable functions of v i,k , each defined some neighborhood of [T k , T k+1 ] × X i and the continuous functions w i , each defined on X i . Let us now show that the modified problem (10) provides a guaranteed outer approximation of X 0 which can be defined asX Theorem 1: For any pair (v, w) feasible in (10), it holds that v(0, ·) ≥ 0 on X 0 andX v,0 ⊃ X 0 . Proof: We first need to show that the discontinuous function v(t, x(t)) is decreasing along the system trajectories. That is, given two time instants t α ≤ t β , we want to show that If we v were to be differentiable, this follows by simply integrating the first constraint of (4) along a trajectory.
Let us now show that (12) holds even for the discontinuous v(t, x(t)) as defined in (8). The first constraint of (10) ensures that the values of v decrease whenever the trajectory resides in the interior of one of the sets X i . We therefore need to argue only about what happens on the boundary of these sets.
a) Time splits: The result for time splits follows immediately from the fifth constraint of (10).
b) State-space splits: Let us first assume that the state space X is split into two parts, X α and X β by a hyper-plane with normal vector h, pointing from X α to X β , so that for x α ∈ X α and x β ∈ X β . The function v(t, x(t)), now split between X α and X β , is defined as Let x 0 ∈ X 0 ∩ X α , t α ∈ [0, T ] and u(·) be given. Let x(·|x 0 ) be the trajectory starting at t α generated by u(·) and suppose that x(t | x 0 ) ∈ X, u(t) ∈ U for t ∈ [t α , T ]. Assume further that this trajectory crosses from X α to X β at the crossing time and assume that this trajectory stays in The last constraint of (10) implies that whereas the first constraint implies and Let us now calculate the value of By inspecting (16), (17) and (17), we can see that the last three summands are nonpositive and we get the inequality v β (t β , x(t β )) ≤ v α (t α , x(t α )), which is equivalent to (12), recalling the definition of v(t, x(t)) in (14). The procedure for the negative trajectory direction is analogous. We note that this analysis encompasses the subtle case of the trajectory sliding on the boundary between the sets X α and X β . By induction, we can prove the inequality for arbitrary splitting of the state-space and time axis by considering a sequence of crossing times associated to a given trajectory. Therefore v(t β , x(t β )) ≤ v(t α , x(t α )) for any 0 ≤ t α ≤ t β ≤ T . By setting t α = 0, t β = T , and using the constraints of (10), we get v(T, for any x 0 ∈ X 0 as desired. This also implies that x 0 ∈X v,0 and henceX v,0 ⊃ X 0 . IV. SOS REPRESENTATION We can now obtain the SDP representation of (10) by applying Putinar's Positivstellensatz [11]. For example, given polynomials c(x) and g(x) the inequality is implied by where q(x) and s(x) are sum-of-squares polynomials. The condition that a polynomial s of degree 2d is sum-of-squares is in turn equivalent to is a basis of polynomials up degree d and hence this constraint is SDP representable. The SOS approximation of (10) reads where z = [t, x, u] , w i (x) and v i,k (t, x) are polynomials, w i is a vector of coefficients of w i (x) and l i is a vector of Lebesgue measure moments indexed with respect to the same basis as the coefficients of w i . The decision variables in the problem are the polynomials v i,k and w i as well as sum-of-squares multipliers q, s and s. The symbols g X i , g X T i , g U i and g τ k denote the column vectors of polynomials describing the sets X i , X T ∩ X i , U and [T k , T k+1 ] in that order. The degrees of all polynomial decision variables is chosen such that the degrees of all polynomials appearing in (23) do not exceed a given relaxation order d. This is a design parameter controlling the accuracy of the approximation.
Given the picewise polynomial functions (w d , v d ) of degree d constructed from a solution to (23) as in (7) and (8), the outer approximation to the ROA is defined by Convergence of the SDP approximations holds under the classical Archimedianity assumption, e.g., [6, Assumption 3].
Theorem 2: For each d ∈ N, we have X d ⊃ X 0 . If in addition the algebraic description of each element of the space-time partition used in (23) satisfies the Archimedianity condition, then lim d→∞ λ(X d \ X 0 ) = 0.
Proof: The result follows from [6, Theorem 6] by taking all functions defined on the partition equal, i.e., v i,k = v and w i = w for some polynomials v and w; this leads to the setting of [6, Theorem 6].

A. Practical implications
The ROA with splits is expected to improve accuracy of the original formulation by allowing one to trade off the degree of the polynomials for number of splits. By increasing the degree d, the size of the SDP will increase with the rate of the binomial coefficient n+d n . By fixing d and splitting the state-space into λ cells, the SDP size will grow only linearly with the number of the cells with rate λ n+d n . The numerical results in the following section suggest that even the computation time grows linearly. A full quantitative understanding of this tradeoff is beyond the scope of the present work; the first step in this direction would be the convergence rate bounds of [12].
V. NUMERICAL EXAMPLES This section presents numerical examples, showcasing the performance difference between the proposed method and the original approach from [6].
The first example in Section V-A shows the influence of the split positions on the resulting ROA. It is shown that by having the splits exactly at the boundaries of the ROA, we can retrieve the theoretical indicator function I X 0 (x).
The second example Section V-B benchmarks the algorithm on a Brockett integrator, which mimics a kinematic model of a nonholonomic system (it can be shown that three-dimensional nonholonomic vehicle with two inputs can be transformed into the Brockett integrator [13]).
Finally, Section V-C presents a comparison of computational demands of the proposed method and the original one from [6].
All the examples were implemented in MATLAB [14] with the use of YALMIP [15]. The YALMIP's sum-of-squares package [16] was used for rapid prototyping; for larger examples, the SDPs were assembled using a custom routine. All SDP's were solved by MOSEK [17].

A. Univariate cubic dynamics
This example shows that one can find the ROA with a very low degree polynomials by correctly positioning the splits.
The system in question is defined aṡ In Fig. 1, we can see a comparison between the original method (without splits) and multiple calculations with splits, going from the inside of the real ROA to the outside. The ROA estimates here are given by {x : w(x) ≥ 1}, which follows from (11) and (10). We can observe that the ROA estimates get more precise, the closer the splits are to the real ROA. Let us define the estimate of I X 0 (x) asĪ X 0 (x), which takes 1 on {x : w(x) ≥ 1} and 0 otherwise. We can observe from Fig. 1, that for the exact split it holds thatĪ X 0 = I X 0 and we obtain the theoretically optimal estimate. This example shows that our method can be used in an iterative manner with splits along an inner approximation of the ROA (such as the one in [18]) as a starting point.

B. Brockett integrator
The Brockett integrator is defined according to [19] aṡ With X = {x ∈ R 3 : ||x|| ∞ ≤ 1}, X T = {0}, U = {u ∈ R 2 : ||u|| 2 ≤ 1}, and T = 1. As was stated before, this system usually serves as a benchmark for nonholonomic control strategies, because it is the simplest system for which there exists no continuous control law which would make the origin asymptotically stable [19]. We shall use the system for calculation of the controlled ROA, which can be computed analytically [20] as The Fig. 2 shows that given a fixed time for the calculation, the proposed approach is always better than the original one and that the split-method approaches the real volume much faster, although neither of the two methods reached the real volume, due to memory constraints.
A visual example of the difference between the two methods can be seen in the Fig. 3 where two ROA's with the same computation time are compared, with one being calculated by the original method and the other by the proposed method. We can see that there is a notable difference between the two approximations, and that the better approximation is done by the lower-degree polynomials. C. Performance and scalability 1) Problem size: First, we shall investigate the accuracy of the algorithm with increasing size of the SDP. The problem size is measured as the number of nonzero elements in the A matrix of the SDP min c x for variable x ∈ R n , convex cone K and data A ∈ R m×n , b ∈ R m , and c ∈ R n . We can clearly see in Fig. 4 that the split versions are always more precise than the non-split version with the same memory footprint.
2) Computation time: The problem size increases linearly with the number of cells (as was explained in IV-A), but the computation time does not necessarily have to follow the same pattern. In this case, however, the computation time also showed linear growth as can be seen in the Figures 5 and 6 for the Brockett integrator and the Double integrator respectively. The Double integrator is defined asẋ [6, 9.3] for more details. In the Fig. 6, the variables to be split were chosen randomly and the splits were always halving the largest interval of the randomly selected variable. This was done in order to ensure that the linear growth is not simply a fortunate result of a particular split order.

VI. CONCLUSION
We showed that the convex SOS-based calculation of ROA can be extended by splitting the time and state space to achieve better accuracy with lower degree polynomials, while keeping the outer-approximation guarantees. The method is faster and provides more accurate estimates of the ROA. It was also empirically demonstrated that the computation time grows linearly with increasing number of cells. Blue approximation is by single degree 10 polynomial (754s) and the red is by sixteen polynomials of degree 8 (753s). The red, lower-degree, approximation is visibly closer to the real ROA.  VII. FUTURE WORK The sparse problem structure can be exploited by ADMM-like algorithm, where the ADMM could alternate between solving small independent SDPs, and enforcing the flow constraints between them. This would allow for parallelization and significantly decrease the required computation time.
The splits can be found in an optimization-based fashion by differentiating the SDP [21] with respect to the location (or more general parametrization of the splits), thereby increasing the accuracy while keeping the memory requirements constant. This also leads to the option of creating the splits in a less naive fashion than the one presented here. Using inner-approximation of the ROA lends itself as an obvious candidate, due to the behaviour shown in V-A.
The same splitting approach can also be applied to other SOS-based algorithms, e.g., to the optimal control problem [22], where the increased precision would lead to more precise estimates of the optimal control sequences and also allow for control of higher dimensional systems. When applied to the problem of the maximum control invariant set estimation [23], temporal splitting is no longer possible but could be replaced by introducing multiple discount factors.
Finally, a natural way to decrease the computational complexity further is to combine the approach with the sparsity-based decomposition method of [24], where the time and state-space splitting would be applied to the subsystems that cannot decomposed using the method of [24].