Semidefinite Outer Approximation of the Backward Reachable Set of Discrete-time Autonomous Polynomial Systems

Semidefinite Outer Approximation of the Backward Reachable Set of Discrete-time Autonomous Polynomial Systems

Weiqiao Han and Russ Tedrake Computer Science and Artificial Intelligence Laboratory, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA. weiqiaoh,russt@mit.edu
Abstract

We approximate the backward reachable set of discrete-time autonomous polynomial systems using the recently developed occupation measure approach. We formulate the problem as an infinite-dimensional linear programming (LP) problem on measures and its dual on continuous functions. Then we approximate the LP by a hierarchy of finite-dimensional semidefinite programming (SDP) programs on moments of measures and their duals on sums-of-squares 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 semi-algebraic set under a polynomial map.

I Introduction

Given a discrete-time 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 continuous-time polynomial systems as an infinite-dimensional 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 continuous-time polynomial systems [7], inner approximation of the ROA of continuous-time controlled polynomial systems [8], approximation of the maximum controlled invariant set of discrete-time and continuous-time polynomial systems [9], controller synthesis for continuous-time polynomial systems [14], approximation of the backward reachable set of continuous-time controlled hybrid polynomial systems [18], and controller synthesis for continuous-time hybrid polynomial systems [22].

Despite the rich literature regarding continuous-time polynomial systems, studies on discrete-time polynomial systems are sparse. We are interested in discrete-time polynomial systems. The most important reason is that robot control is effectively always implemented as sampled-data control with the controller implemented on a digital computer. Additionally, when modeling the robot making and breaking contacts with the environment, the continuous-time formulation requires dealing with measure differential inclusions for impacts [17]. In contrast, time-stepping 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 decision-making is discrete on a footstep-to-footstep 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 discrete-time 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 discrete-time world. Not until 2016 did Magron et al. devise the discrete-time version of the Liouville equation for automonous polynomial systems [12]. They used it to approximate the forward reachable set for discrete-time polynomial systems.

In this paper, we adopt the occupation measure approach to over approximate the backward reachable set of discrete-time autonomous polynomial systems, a step towards our ultimate goal of controller synthesis for discrete-time hybrid polynomial systems. We follow the Lasserre hierarchy strategy to first formulate the problem as an infinite-dimensional LP on measures and its dual on continuous functions, and then approximate the LP by a hierarchy of finite-dimensional SDP programs on moments of measures and their duals on sums-of-squares (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 non-convex problems, and theoretically the approximations to the real set are closer and closer with higher and higher degrees of approximation. 111Unlike the continuous-time case [5], as far as we know, no one has established the vanishing conservatism for any discrete-time 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 time-reverse version of [12] and [13], combined with ideas in [8], and a pathway towards controller synthesis.

Ii Problem formulation

Ii-a Problem statement

Let . (resp. ) stands for the set of polynomials (resp. of degree at most ) in the variable . Given the discrete-time 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 semi-algebraic 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 semi-algebraic 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 state-feedback control inputs. This is immediately seen by plugging the polynomial control law into the control affine polynomial system , yielding a polynomial closed-loop 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 low-order Taylor expansion of the closed-loop dynamical system. This is illustrated in the last example in Section V.

Ii-B 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 sup-norm. 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 non-negative elements of (resp. ). The topology in is the strong topology of uniform convergence while the topology in is the weak-star 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 semi-algebraic 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. ).

    Ii-C Discrete-time Liouville equation

    Let and be a polynomial map. The pushforward measure is defined to be

    for all .

    The Liouville equation for discrete-time 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

    Iii-a 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).

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

    2. 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 .

    3. 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.

    Iii-B Primal-dual infinite-dimensional 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 0-superlevel set of the function , , is an outer approximation of . If, in addition , then the 1-superlevel set of the function , , is an outer approximation of .

    Iv Semidefinite relaxations

    We have formulated the infinite-dimensional 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 infinite-dimensional LP on measures by a hierarchy of finite-dimensional 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.

    Iv-a 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 multi-indices (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 multi-indices (rows) and (columns)

    If has a representing measure , then whenever the support of is contained in . Conversely, if is a compact semi-algebraic 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]).

    Iv-B Primal-dual finite-dimensional 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 infinite-dimensional 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 0-superlevel set of the function , and over approximate the backward reachable set simply by the 1-superlevel set of if , which in practice is usually the case returned by the solver or can be set manually. Unlike the continuous-time 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 near-optimal 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].

      V-a 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 .

      V-B Primal-dual 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 primal-dual SDP’s accordingly, and the 1-superlevel 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 semi-algebraic set and approximate the backward reachable set for three discrete-time 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 .

        Vi-a Preimage of a semi-algebraic set

        Consider the polynomial map (Example 1 in [13]). Choose and .

        We approximate the preimage by degree-12 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 1-superlevel 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: degree-12 approximation of the preimage in gray. Right: degree-12 approximation of the preimage in gray and the true preimage in red.

        Vi-B Cathala system

        Consider the Cathala System (Example 5.2 in [12]):

        and . We approximate the backward reachable set by degree-10 and degree-20 polynomials. As show in Figure 2, the gray areas are the 1-superlevel 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 degree-10 approximation is already very good and the degree-20 approximation is even more precise.

        Fig. 2: Left: degree-10 approximation. Right: degree-20 approximation.

        Vi-C Van der Pol oscillator

        Consider the uncontrolled reversed-time Van der Pol oscillator (Example 9.2 in [5]) given by

        Discretizing the model with the explicit Euler scheme with a sampling time , the discrete-time system is

        Choose and .

        We approximate the backward reachable set by degree-14, 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: Degree-14, 16, 18, and 20 approximations

        Vi-D Cart-pole system under LQR control

        Fig. 4: The cart-pole system.

        Consider balancing the cart-pole 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 2-dimensional 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 2-dimensional section of the true backward reachable set. The degree-8 approximation is already good. The degree-10 approximation is a little more closer. The computation time of the degree-10 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 cart-pole system
        Fig. 5: Left: degree-8 approximation. Right: degree-10 approximation.

        Vii Conclusion

        We have used recently developed occupation measure approach to over approximate the backward reachable set of discrete-time autonomous polynomial systems. As shown in the last example, the approach can also be used to approximate the backward controllable set of discrete-time controlled polynomial systems subject to already-designed 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 discrete-time polynomial systems and for discrete-time 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 discrete-time 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. W911NF-15-1-0166.

        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 multi-contact push recovery via lmi approximation of the piecewise-affine quadratic regulator. In Humanoid Robotics (Humanoids), 2017 IEEE-RAS 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: 2018-02-28.
        • [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, Jean-Bernard Lasserre, and Johan Löfberg. Gloptipoly 3: moments, optimization and semidefinite programming. Optimization Methods & Software, 24(4-5):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 lmi-relaxations. SIAM journal on control and optimization, 47(4):1643–1666, 2008.
        • [11] Jean-Bernard Lasserre. Moments, positive polynomials and their applications, volume 1. World Scientific, 2010.
        • [12] Victor Magron, Pierre-Loïc Garoche, Didier Henrion, and Xavier Thirioux. Semidefinite approximations of reachable sets for discrete-time polynomial systems. arXiv preprint arXiv:1703.05085, 2017.
        • [13] Victor Magron, Didier Henrion, and Jean-Bernard 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 multi-contact push recovery in complex environments. In Humanoid Robotics (Humanoids), 2017 IEEE-RAS 17th International Conference on, pages 31–38. IEEE, 2017.
        • [16] Michael Posa, Twan Koolen, and Russ Tedrake. Balancing and step recovery capturability via sums-of-squares optimization. In Robotics: Science and Systems, 2017.
        • [17] Michael Posa, Mark Tobenkin, and Russ Tedrake. Stability analysis and control of rigid-body 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 time-stepping 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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
131607
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description