Semidefinite Outer Approximation of the Backward Reachable Set of Discretetime Autonomous Polynomial Systems
Abstract
We approximate the backward reachable set of discretetime autonomous polynomial systems using the recently developed occupation measure approach. We formulate the problem as an infinitedimensional linear programming (LP) problem on measures and its dual on continuous functions. Then we approximate the LP by a hierarchy of finitedimensional semidefinite programming (SDP) programs on moments of measures and their duals on sumsofsquares polynomials. Finally we solve the SDP’s and obtain a sequence of outer approximations of the backward reachable set. We demonstrate our approach on three dynamical systems. As a special case, we also show how to approximate the preimage of a compact semialgebraic set under a polynomial map.
I Introduction
Given a discretetime autonomous polynomial dynamical system and a target set in state space, we are interested in the initial set from which the system evolves into the target set in finite time steps without leaving the state constraint set. We call this set the backward reachable set, also known as the region of attraction (ROA).
Estimating the backward reachable set is a fundamental problem in robotics and control. Henrion et al in 2012 formulated the computation of the ROA of continuoustime polynomial systems as an infinitedimensional linear programming (LP) problem in the cone of nonnegative Borel measures, and then solved the LP approximately using a hierarchy of finite dimensional semidefinite programming (SDP) relaxations [5]. Their approach, also known as the occupation measure approach (or the Lasserre hierarchy strategy on occupation measure [4]), originally proposed in [10], has the benefits of convexity and asymptotically vanishing conservatism, compared to other existing approaches; see [5] and the references therein. Since then, the occupation measure approach has been used widely, e.g., for inner approximation of the ROA of continuoustime polynomial systems [7], inner approximation of the ROA of continuoustime controlled polynomial systems [8], approximation of the maximum controlled invariant set of discretetime and continuoustime polynomial systems [9], controller synthesis for continuoustime polynomial systems [14], approximation of the backward reachable set of continuoustime controlled hybrid polynomial systems [18], and controller synthesis for continuoustime hybrid polynomial systems [22].
Despite the rich literature regarding continuoustime polynomial systems, studies on discretetime polynomial systems are sparse. We are interested in discretetime polynomial systems. The most important reason is that robot control is effectively always implemented as sampleddata control with the controller implemented on a digital computer. Additionally, when modeling the robot making and breaking contacts with the environment, the continuoustime formulation requires dealing with measure differential inclusions for impacts [17]. In contrast, timestepping methods [19] make it possible to capture the the complexity of the constrained hybrid dynamics and impact efficiently, without worrying about impulsive events and event detection [3, 15]. As another related example, in step capturability analysis, used to study balancing in legged robots, the decisionmaking is discrete on a footsteptofootstep level, and the entire problem formulation is asking about the viability kernel, also known as the backward reachable set [16].
The reason for the sparse literature in the discretetime world is probably that, in the occupation measure approach for dynamical systems, the LP formulation requires the Liouville equation, which is a partial differential equation describing the evolution of the state distribution over time, but there is no natural analogue of the Liouville equation in the discretetime world. Not until 2016 did Magron et al. devise the discretetime version of the Liouville equation for automonous polynomial systems [12]. They used it to approximate the forward reachable set for discretetime polynomial systems.
In this paper, we adopt the occupation measure approach to over approximate the backward reachable set of discretetime autonomous polynomial systems, a step towards our ultimate goal of controller synthesis for discretetime hybrid polynomial systems. We follow the Lasserre hierarchy strategy to first formulate the problem as an infinitedimensional LP on measures and its dual on continuous functions, and then approximate the LP by a hierarchy of finitedimensional SDP programs on moments of measures and their duals on sumsofsquares (SOS) polynomials. By solving the hierarchy of SDP’s we obtain a sequence of outer approximations of the backward reachable set. The benefit of this framework is that we solve convex problems instead of generally nonconvex problems, and theoretically the approximations to the real set are closer and closer with higher and higher degrees of approximation. ^{1}^{1}1Unlike the continuoustime case [5], as far as we know, no one has established the vanishing conservatism for any discretetime problems yet. See Remark 1 under Proposition 1. Besides the most natural optimization formulation, we introduce a few variants based on ideas in [8]. We also show that, as a special case, the “preimage problem”, the problem of computing the outer approximation of the preimage of a compact semialgebraic set under a polynomial map, can be solved easily with any of the optimization formulations. Our work can be viewed as the timereverse version of [12] and [13], combined with ideas in [8], and a pathway towards controller synthesis.
Ii Problem formulation
Iia Problem statement
Let . (resp. ) stands for the set of polynomials (resp. of degree at most ) in the variable . Given the discretetime autonomous polynomial system
where and is a polynomial map, a time step , and a target set , define the step backward reachable set
This is the set of points in that enter the target region within time steps and whose trajectories do not leave before entering . Once a point enters , what happens to it next is not the concern of this paper. We are interested in the backward reachable set
which is the union of all points in that enter in finite time.
Assume
is a compact semialgebraic set. Furthermore, assume is “simple”, in the sense that the moments of the Lebesgue measure on are available, e.g., can be a ball or a box.
Assume
is a compact semialgebraic set. In practice, might be a small region around the origin, so when the system enters , many control methods, such as linear quadratic regulator (LQR) and feedback linearization, can regulate the system to the origin.
While the approach described in the paper approximates the backward reachable set of autonomous systems, it can also approximate the backward controllable set of systems subject to polynomial statefeedback control inputs. This is immediately seen by plugging the polynomial control law into the control affine polynomial system , yielding a polynomial closedloop dynamical system. For example, suppose we linearize a nonlinear dynamical system around an unstable fixed point and compute the LQR controller around that point. Then we can approximate the backward controllable set of the LQR controller by applying the approach in this paper to some loworder Taylor expansion of the closedloop dynamical system. This is illustrated in the last example in Section V.
IiB Notations and assumptions
In this subsection, we introduce some notations in measure theory and functional analysis. We also make some assumptions on the systems we work with. For an introduction to real analysis, functional analysis, and polynomial optimization, please refer to [2], [1], and [11], respectively.
Let be a compact set. is the Banach space of continuous functions on equipped with the supnorm. denotes the topological dual of , i.e., the set of all continuous linear functionals on . is the Banach space of finite signed Radon measures on the Borel algebra equipped with the total variation norm. By Riesz Representation Theorem, is isometrically isomorphic to . (resp. ) denotes the cone of nonnegative elements of (resp. ). The topology in is the strong topology of uniform convergence while the topology in is the weakstar topology. For any , stands for the restriction of the Lebesgue measure on . For , we say is dominated by , denoted by , if .
Let , and , where ’s (resp. ’s) are the defining polynomials of (resp. ). (resp. ) denotes the cone of SOS polynomials (resp. SOS polynomials of degree up to ) in the variable . (resp. ) denotes the truncated quadratic module generated by the defining polynomials of (resp. ), assuming (resp. ):
Later when we apply the Lasserre hierarchy strategy to formulate SDP approximations on moments, we want the equivalence between positivity on a compact semialgebraic set and the existence of SOS representations, which was established by Putinar’s Positivstellensatz (Section 2.5 in [11]). So we make the following assumption:

For any , there exists (resp. ) such that
The assumption is easily satisfied by setting one of ’s (resp. ’s) to be (resp. ).
IiC Discretetime Liouville equation
Let and be a polynomial map. The pushforward measure is defined to be
for all .
The Liouville equation for discretetime autonomous polynomial systems proposed in [12] is
where . We can think of the initial measure as the distribution of the mass of the initial states of the system trajectories (not necessarily normalized to 1), the occupation measure as describing the volume occupied by the state trajectories, and the final measure as the distribution of the mass of the final states of the system trajectories. For example, , , and is a solution to the Liouville equation, describing the system trajectory , where is the Dirac measure centered on .
Iii Optimization formulation
Iiia The natural LP formulation and its variants
In order to approximate the backward reachable set, the natural optimization formulation is to maximize the mass of the initial measure as in [5], i.e., the following LP:
(1) The first constraint is the Liouville equation. The final measure is constrained to be supported on , because we want the system to end up in . The second constraint ensures that the initial measure is dominated by the Lebesgue measure on . The third constraint bounds the mass of the measure by times the volume of . The LP (1) admits an optimal solution that has the following properties: (i) is the restriction of the Lebesgue measure on a set containing the step backward reachable set ; (ii) the objective , i.e., the mass of , is the volume of the set .
Now we introduce some variants of the LP (1).

Replace the objective by , i.e., the mass of the final measure. The LP still admits the optimal solution with properties described above.

Use the trick in [8] to write the final measure as a summation of two measures , , and in the objective maximize the mass of while minimizing the mass of .

In (b), replace the mass of in the objective by the mass of .
In the following, we take the variant (b) as an example, and work out its dual LP on continuous functions and its SDP relaxations.
IiiB Primaldual infinitedimensional LP
Given , we consider the LP:
(2) The supremum is taken over . As remarked in [8], it is unnecessary to restrict the support of to be by virtue of optimality in (2). The LP (2) admits an optimal solution that has the following properties: (i) the support of measure is empty; (ii) the mass of is the same as the mass of , i.e., ; (iii) is the restriction of the Lebesgue measure on a set containing the step backward reachable set .
The dual of the LP (2) is given by
(3) The infimum is taken over . Let be a feasible solution to the LP (3). Suppose is an admissible trajectory such that . Then and . So . Therefore, the 0superlevel set of the function , , is an outer approximation of . If, in addition , then the 1superlevel set of the function , , is an outer approximation of .
Iv Semidefinite relaxations
We have formulated the infinitedimensional LP on measures and its dual on continuous functions. However, we cannot solve them directly. The Lasserre hierarchy strategy is then to approximate the original infinitedimensional LP on measures by a hierarchy of finitedimensional SDP’s on moments of measures and their dual on SOS polynomials. By solving the hierarchy of relaxed SDP’s, we get a sequence of outer approximations to the backward reachable set.
Iva Preliminaries
Given , define .
For any set , denotes the indicator function on , i.e., if , and otherwise.
Any polynomial can be expressed in the monomial basis as
where , and can be identified with its vector of coefficients indexed by .
Given a sequence of real numbers , we define the linear functional by
If is a sequence of moments for some measure , i.e.,
then is called a representing measure for .
Define the moment matrix of order with entries indexed by multiindices (rows) and (columns)
If has a representing measure, then . However, the converse is generally not true.
Given a polynomial with coefficient vector , define the localizing matrix w.r.t. and to be the matrix indexed by multiindices (rows) and (columns)
If has a representing measure , then whenever the support of is contained in . Conversely, if is a compact semialgebraic set as defined in Section II, if Assumption 1 holds, and if , then has a finite Borel representing measure with support contained in (Theorem 3.8(b) in [11]).
IvB Primaldual finitedimensional SDP
For each , let , be the finite sequence of moments up to degree of the measure . Similarly, , and are finite sequences of moments up to degree associated with measures , , and , respectively. The infinitedimensional LP on measures (2) can be relaxed with the following semidefinite program on moments of measures:
(4) The dual of (4) is the following SDP on polynomials of degrees up to :
(5) 
Let (u_r,v_r,w_r)X_0.0r^T ⊇X_0^Tu_r = 0X_0r^T ⊇X_0^∞
Proof.
See the discussion under the LP (3). ∎
Remark 1. The proposition shows that we can over approximate the step backward reachable set by , the 0superlevel set of the function , and over approximate the backward reachable set simply by the 1superlevel set of if , which in practice is usually the case returned by the solver or can be set manually. Unlike the continuoustime case, the vanishing conservatism, e.g., a sequence of solutions whose component converges to in norm:
where is the closure of , has not been established yet, though we believe this should be the case when is an optimal or nearoptimal solution to SDP (5) for each . If the proof to the Theorem 4.1 item 2 in [12] were correct, then this can be proved analogously. Unfortunately, the proof there is incorrect, because the authors mixed up “pointwise convergence” with “uniform convergence”.
Remark 2. If we define
then the approximation by the sequence of sets is monotone.
V Special case: the preimage problem
In this section, we consider the preimage problem, a special case of the backward reachable set problem. We first define the problem, and then formulate the primal and dual LP’s. An example will be illustrated in the next section. This section can be considered as the reverse of part of [13].
Va Problem statement
Given a polynomial map , and a target set , where and are as defined in Section II.A, we are interested in approximating the preimage .
VB Primaldual LP’s
The LP (1) and its variants have degenerate forms that can handle this special case. Take LP (1) as an example, we formulate the LP as
(6) The dual is given by
(7) We can define the relaxed primaldual SDP’s accordingly, and the 1superlevel set of is an outer approximation of the preimage. What is interesting is that under certain assumptions the optimal solution to the relaxed dual SDP exists, and the sequence of outer approximations by the optimal solutions has vanishing conservatism. To be precise, for , the relaxed dual SDP is given by
(8) 
Let f^0.1(Z)X\f^1(Z)(v_r,w_r){w_r}_r1_f^1(Z)(x)L_1w_r(x) ≥1_f^1(Z)(x), ∀x∈X
Proof.
Analogous to the proof of Theorem 4.4 item 2 in [13]. ∎
Vi Examples
We approximate the preimage of a compact semialgebraic set and approximate the backward reachable set for three discretetime polynomial systems. All computations are done using MATLAB 2016b and the SDP solver MOSEK 8. In terms of the polynomial optimization toolbox, we used GloptiPoly 3 [6] for the first three examples, and turned to Spotless [21] for the last example, because GloptiPoly 3 was too slow on this example. In the last three examples, up to numerical precisions, so we set as in [12]. Therefore, in these examples, we are approximating .
Via Preimage of a semialgebraic set
Consider the polynomial map (Example 1 in [13]). Choose and .
We approximate the preimage by degree12 polynomials. We plotted the approximation set without and with the true preimage in figure 1. As shown in the figure, the gray area is the 1superlevel set of the component of the optimal solution, and therefore is the outer approximation of the preimage. The red area is the true preimage. Notice that we only consider the approximation inside , i.e., the area inside the dashed circle. So the approximation is accurate.
Fig. 1: Left: degree12 approximation of the preimage in gray. Right: degree12 approximation of the preimage in gray and the true preimage in red. ViB Cathala system
Consider the Cathala System (Example 5.2 in [12]):
and . We approximate the backward reachable set by degree10 and degree20 polynomials. As show in Figure 2, the gray areas are the 1superlevel set of the component of the optimal solution and hence approximate the backward reachable set. The red dots are the sampled points in the true backward reachable set. Since the number of samples is large, the red regions (the red disc in the middle plus the two small red regions just inside the boundary of the unit circle) in the figure represent the true backward reachable set. Notice that we only consider the region in , i.e., the region inside the dashed unit circle in the figures. The degree10 approximation is already very good and the degree20 approximation is even more precise.
Fig. 2: Left: degree10 approximation. Right: degree20 approximation. ViC Van der Pol oscillator
Consider the uncontrolled reversedtime Van der Pol oscillator (Example 9.2 in [5]) given by
Discretizing the model with the explicit Euler scheme with a sampling time , the discretetime system is
Choose and .
We approximate the backward reachable set by degree14, 16, 18, and 20 polynomials. As show in Figure 3, the gray areas are the approximate backward reachable sets. The areas enclosed by the red lines are the true backward reachable set, which was obtained analytically by integrating backwards in time.
Fig. 3: Degree14, 16, 18, and 20 approximations ViD Cartpole system under LQR control
Consider balancing the cartpole system [20] as shown in Figure 4. We can apply only horizontal force on the cart, so the system is underactuated. Let , and . We are interested in approximating the backward reachable set around its unstable fixed point when the system is under LQR control.
First, we linearize the system around the unstable fixed point and compute the LQR controller with . The physical parameters of the system are the same as in Example 1 of [3] and are shown in Table I. Second, we Taylor expand the closed loop dynamics about up to order 3. This gives us 3rd degree polynomial system dynamics. Finally, we approximate the backward reachable set with and .
The approximate backward reachable set is a subset in . We plot its 2dimensional section by setting and , as shown in Figure 5. The gray area is the section of the approximate backward reachable set, while the red dots are the sampled points in the true backward reachable set. Since the number of samples is large, the red regions in the figure represent the 2dimensional section of the true backward reachable set. The degree8 approximation is already good. The degree10 approximation is a little more closer. The computation time of the degree10 approximation is 3326 seconds.
Explanation Value cart mass 10 pole mass 1 pole length 0.5 gravitational acceleration 9.81 discretization time interval 0.05 TABLE I: Physical parameters of the cartpole system Fig. 5: Left: degree8 approximation. Right: degree10 approximation. Vii Conclusion
We have used recently developed occupation measure approach to over approximate the backward reachable set of discretetime autonomous polynomial systems. As shown in the last example, the approach can also be used to approximate the backward controllable set of discretetime controlled polynomial systems subject to alreadydesigned state feedback control laws. It is possible to inner approximate the backward reachable set by the strategy in [7]. However, what we are most interested in is the controller design for discretetime polynomial systems and for discretetime hybrid polynomial systems via the occupation measure approach, which could then be potentially used for humanoid push recovery by making and breaking multiple contacts with the environment [3, 15]. This may require devising a new form of discretetime Liouville equations, which shall be the focus of the future research.
Acknowledgment
This work was supported by Air Force/Lincoln Laboratory Award No. 7000374874 and Army Research Office Award No. W911NF1510166.
References
 [1] John B Conway. A course in functional analysis, volume 96. Springer Science & Business Media, 2013.
 [2] Gerald B Folland. Real analysis: modern techniques and their applications. John Wiley & Sons, 2013.
 [3] Weiqiao Han and Russ Tedrake. Feedback design for multicontact push recovery via lmi approximation of the piecewiseaffine quadratic regulator. In Humanoid Robotics (Humanoids), 2017 IEEERAS 17th International Conference on, pages 842–849. IEEE, 2017.
 [4] Didier Henrion. The lasserre hierarchy in robotics. http://webdav.tuebingen.mpg.de/robust_mpc_legged_robots/henrion_slides.pdf, May 2016. Accessed: 20180228.
 [5] Didier Henrion and Milan Korda. Convex computation of the region of attraction of polynomial control systems. IEEE Transactions on Automatic Control, 59(2):297–312, 2014.
 [6] Didier Henrion, JeanBernard Lasserre, and Johan Löfberg. Gloptipoly 3: moments, optimization and semidefinite programming. Optimization Methods & Software, 24(45):761–779, 2009.
 [7] Milan Korda, Didier Henrion, and Colin N Jones. Inner approximations of the region of attraction for polynomial dynamical systems. IFAC Proceedings Volumes, 46(23):534–539, 2013.
 [8] Milan Korda, Didier Henrion, and Colin N Jones. Controller design and region of attraction estimation for nonlinear dynamical systems. IFAC Proceedings Volumes, 47(3):2310–2316, 2014.
 [9] Milan Korda, Didier Henrion, and Colin N Jones. Convex computation of the maximum controlled invariant set for polynomial control systems. SIAM Journal on Control and Optimization, 52(5):2944–2969, 2014.
 [10] Jean B Lasserre, Didier Henrion, Christophe Prieur, and Emmanuel Trélat. Nonlinear optimal control via occupation measures and lmirelaxations. SIAM journal on control and optimization, 47(4):1643–1666, 2008.
 [11] JeanBernard Lasserre. Moments, positive polynomials and their applications, volume 1. World Scientific, 2010.
 [12] Victor Magron, PierreLoïc Garoche, Didier Henrion, and Xavier Thirioux. Semidefinite approximations of reachable sets for discretetime polynomial systems. arXiv preprint arXiv:1703.05085, 2017.
 [13] Victor Magron, Didier Henrion, and JeanBernard Lasserre. Semidefinite approximations of projections and polynomial images of semialgebraic sets. SIAM Journal on Optimization, 25(4):2143–2164, 2015.
 [14] Anirudha Majumdar, Ram Vasudevan, Mark M Tobenkin, and Russ Tedrake. Convex optimization of nonlinear feedback controllers via occupation measures. The International Journal of Robotics Research, 33(9):1209–1230, 2014.
 [15] Tobia Marcucci, Robin Deits, Marco Gabiccini, Antonio Biechi, and Russ Tedrake. Approximate hybrid model predictive control for multicontact push recovery in complex environments. In Humanoid Robotics (Humanoids), 2017 IEEERAS 17th International Conference on, pages 31–38. IEEE, 2017.
 [16] Michael Posa, Twan Koolen, and Russ Tedrake. Balancing and step recovery capturability via sumsofsquares optimization. In Robotics: Science and Systems, 2017.
 [17] Michael Posa, Mark Tobenkin, and Russ Tedrake. Stability analysis and control of rigidbody systems with impacts and friction. IEEE Transactions on Automatic Control, 61(6):1423–1437, 2016.
 [18] Victor Shia, Ram Vasudevan, Ruzena Bajcsy, and Russ Tedrake. Convex computation of the reachable set for controlled polynomial hybrid systems. In Decision and Control (CDC), 2014 IEEE 53rd Annual Conference on, pages 1499–1506. IEEE, 2014.
 [19] David E Stewart and Jeffrey C Trinkle. An implicit timestepping scheme for rigid body dynamics with inelastic collisions and coulomb friction. International Journal for Numerical Methods in Engineering, 39(15):2673–2691, 1996.
 [20] Russ Tedrake. Underactuated robotics: Algorithms for walking, running, swimming, flying, and manipulation (course notes for mit 6.832). Downloaded in Fall, 2014.
 [21] Mark M Tobenkin, Frank Permenter, and Alexandre Megretski. Spotless polynomial and conic optimization, 2013.
 [22] Pengcheng Zhao, Shankar Mohan, and Ram Vasudevan. Optimal control for nonlinear hybrid systems via convex relaxations. arXiv preprint arXiv:1702.04310, 2017.

