Exact and Efficient Hamilton-Jacobi Reachability for Decoupled Systems

Exact and Efficient Hamilton-Jacobi Reachability for Decoupled Systems


Reachability analysis is important for studying optimal control problems and differential games, which are powerful theoretical tools for analyzing and modeling many practical problems in robotics, aircraft control, among other application areas. In reachability analysis, one is interested in computing the reachable set, defined as the set of states from which there exists a control, despite the worst disturbance, that can drive the system into a set of target states. The target states can be used to model either unsafe or desirable configurations, depending on the application. Many Hamilton-Jacobi formulations allow the computation of reachable sets; however, due to the exponential complexity scaling in computation time and space, problems involving approximately 5 dimensions become intractable. A number of methods that compute an approximate solution exist in the literature, but these methods trade off complexity for optimality. In this paper, we eliminate complexity-optimality trade-offs for time-invariant decoupled systems using a decoupled Hamilton-Jacobi formulation that enables the exact reconstruction of high dimensional solutions via low dimensional solutions of the decoupled subsystems. Our formulation is compatible with existing numerical tools, and we show the accuracy, computation benefits, and an application of our novel approach using two numerical examples.

I Introduction

Optimal control problems and differential games have been extensively studied [1, 2, 3, 4], and have received growing interest in the recent past. These powerful theoretical tools allow us to analyze a variety of real world problems, including path planning, collision avoidance, safety verification, among other applications in robotics, aircraft control, security, and other domains [5, 6, 7, 8].

In an optimal control problem, one aims to drive a controlled dynamical system into a set of states called the target set; depending on the application, the target set can model the set of either desirable or undesirable configurations. In a reachability framework, one aims to determine the backwards reachable set, defined as the set of states from which a control exists to drive the system into the target set. Differential games involve two adversarial players (Player 1 and Player 2). Player 2 seeks to drive a system to a target set, while Player 1 seeks to prevent Player 2 from doing so. One again aims to determine the backwards reachable set, which in this case is defined as the set of states from which a control from Player 2 exists to drive the system into the target set, despite the optimal adversarial control from Player 1.

Reachability is an effective way to analyze optimal control problems and differential games because it provides guarantees on system performance and safety. Reachability problems involving one player can be posed as a minimum (maximum) cost game where the player minimizes the minimum value over time of some cost function representing the proximity to the target set. In the case of a differential game, Player 1 maximizes the minimum cost over time, while Player 2 minimizes it. [1] has shown that the backwards reachable set can be obtained by solving a Hamilton-Jacobi Partial Differential Equation (HJ PDE) with a terminal condition specifying the target set. Many similar formulations of the backwards reachability problems also exist [9, 10, 11]. HJ reachability has been successfully used to solve problems such as aircraft collision avoidance [1], automated in-flight refueling [12], and reach-avoid games [13, 14].

The techniques for computing backwards reachable sets via solving an HJ PDE are very flexible and can be applied to a large variety of system dynamics when the problem dimensionality is low. Furthermore, many numerical tools have been developed to solve these equations, making the HJ approach practically appealing [15, 16, 17, 18]. For higher dimensional problems, various techniques such as those involving projections [19, 20], approximate dynamic programming [21], and occupation measure theory [22] have been proposed. While these approximation techniques alleviate the computation complexity, they give up optimality and sometimes give overly conservative results.

This paper resolves the complexity-optimality trade-off for time-invariant systems with decoupled dynamics. We present a decoupled formulation of HJ reachability for decoupled systems, defined in (3). By considering the decoupled component separately and solving lower dimensional HJ PDEs for each subsystem, we reduce the computation time and space complexity substantially. Our approach also exactly recovers the solution to the original, high dimensional PDE.

Ii Problem Formulation

Consider a differential game between two players described by the time-invariant system


where is the system state, is the control of Player 1, and is the control of Player 2. We assume is uniformly continuous, bounded, and Lipschitz continuous in for fixed , and the control functions are drawn from the set of measurable functions1. As in [1, 23, 24], we allow Player 2 to only use nonanticipative strategies , defined by


We further assume that the system is a decoupled system.

Definition 1

Decoupled system. A system (1) is a decoupled system if it can be split into components, denoted where , that satisfy the following:


where is th component of the full state, is th component of the control of Player 1, and is th component of the control of Player 2. Based on this assumption and the assumptions on , we have that is uniformly continuous, bounded, and Lipschitz continuous in for fixed and are measurable. Note that .

Denote system trajectories, which are solutions to (1), as


satisfies initial conditions and the following differential equation almost everywhere


In our differential game, the goal of Player 2 is to drive the system into some target set , and the goal of Player 1 is to drive the system away from it. The set is represented as the zero sublevel set of a bounded, Lipschitz continuous function , .

Such a function always exists, since we can choose to be a signed distance function; we call the implicit surface function representing the set . In accordance with our decoupled dynamics, we assume that can be represented as a maximum of bounded, Lipschitz continuous functions , where are implicit surface functions representing so that . Note that with the definition of and , we have that


Given the decoupled system (3) and the target set in the form (6) represented by , our goal in this paper is to compute the backwards reachable set, , in the low-dimensional space of each of the decoupled components as opposed to in the full system state space . is defined as

Remark 1

One may have noticed that if , one would be able to simply find by solving (8) with as the target set, and then obtain . However, it is crucial to observe that we are interested in the case where , in which a simple union of would not yield the correct reachable set .

Iii Solution

Iii-a HJ Reachability: Full Formulation

In [1], the authors showed that the backwards reachable set can be obtained as the zero sublevel set of the viscosity solution [25] of the following terminal value HJ PDE:


from which we obtain from the bounded, Lipschitz function that is also continuous in both and [23].

[1] and similar approaches, such as [9, 10, 11], are compatible with well-established numerical methods [15, 16, 17, 18]. However, these approaches become intractable quickly as the dimensionality of the problem increases. Numerically, the solution is computed on a grid, and the number of grid points increases exponentially with the number of dimensions.

Decoupled dynamics allow for tractable or faster computation of reachable sets in the individual decoupled components. Some authors [19, 20] have proposed methods for combining or stitching together these reachable set components into the full reachable set. These methods work reasonably well, but introduce conservatism in various ways. In the next subsections, we will provide a method for combining solutions to the lower dimensional HJ PDEs to construct the exact full solution in the original HJ PDE.

Iii-B HJ Reachability: Decoupled Formulation

Observe that (8) can be viewed as an equation involving two cases. Depending on which of the arguments in the outer-most minimum is active, becomes one of (9) or (10):


This motivates us to define the following sets which characterize which of the outer-most minimum operation is active in (8).


Note that is the complement of , , for all time. At a given , in , satisfies (9); in , satisfies (10). We now show an important property of and in the Lemma and Corollary below. These will be used to show that our proposed decoupled formulation allows exact computation of , by computation of , value functions of lower dimensional spaces.

Lemma 1

for some such that .


Suppose , then by (9) we have the following:

  • . Thus, becomes independent of at .

  • since is an open set, there exists a neighborhood around that is contained in . Thus is also independent of in a neighborhood of .

  • By (11), we have .

Let and suppose at . Then, by (11), . Since is independent of , this necessarily means that .

This implies in a neighborhood of such that .


However, by (9), . In particular, then, we have for any


This means that , such that


which is a contradiction since is Lipschitz continuous. Therefore, since , we have that by (11).

Corollary 1



Suppose but . Since is the complement of , this implies . By Lemma 1, we must have that since , a contradiction.

We can now state our main theorem.

Theorem 1

The solution to (8) for a decoupled system with dynamics (3) and terminal condition is given by


where are the viscosity solutions to


and is the smallest time such that , i.e.


Case 1: By Lemma 1, we have . Therefore, satisfies (9) , so . Case 2 of this proof would then imply .

Case 2: Consider a target set represented by the zero sublevel set of the function , where . By (8), we have that . Define functions such that , then at , we have , and


where is the indicator function that is when its argument is true and otherwise, and is an matrix in of all zeros except for in the rows corresponding to the component where it is the identity matrix in .

Now, consider all points , in which satisfies (10). Substituting into (10), we have


Equation (19) states that in the region where is the maximum among , we have , where satisfies (16).

Consider auxiliary functions which satisfy, for all and all ,


By Corollary 1, we have that and both satisfy the same PDE with the same terminal conditions, . Therefore, .

Iii-C Decoupled Formulation Algorithm

Algorithmically, Theorem 1 states the following:

  1. for some .

  2. .

  3. where satisfies (16).

This gives us an efficient way to computed by computing , transforming the original -dimensional problem of computing into the -dimensional problems of computing . Based on the conclusions we drew, the following algorithm exactly computes , which satisfies (8), with the above-mentioned computation benefits:

  1. Initialize .

  2. Compute , by solving (16).

  3. Initialize .

  4. Decrement from 0 to ; for each time step :

    1. Set the auxiliary variable

      This step correctly computes to be for all .

    2. Update the value function

      This step correctly computes to satisfy for all .

Iii-D Computation Time and Space Complexity Comparison

For a state space discretization of grid points in each dimension, the computation time complexity decreases from for the original problem in , to for the subproblems in . This is a computation speed improvement of many orders of magnitude.

Directly solving (8) on a computational domain has a space complexity of , where is the number of time steps of being stored, since we need to store an -dimensional grid for each of the time steps. The algorithm presented in III-C involves computing on computation domains . Each thus has a space complexity of , making the overall space complexity .

From , we can then reconstruct in any domain . Thus, by choosing to be a small subset of , we can always avoid additional space complexity. This allows us to access . In practice, one would choose to be in a small region around a state of interest (eg. the current system state), and access the value function as well as its gradient at ; this allows one to determine of whether is in the reachable set based on the sign of , and compute the optimal controls for both Player 1 and Player 2 respectively based on as follows:


Note that we do not need to store for all . In fact, in many situations, storing in the entire is infeasible since the space complexity is exponential with the dimension of . With that caveat, we restate our algorithm for computing from , explicitly noting memory allocation, to show that one only needs to store for , where is a very small subset of :

  1. Initialize for .

  2. Compute in by solving 16.

  3. Initialize in a small computation domain .

  4. Decrement from 0 to ; for each time step , perform the following computations in :

Iii-E Numerical Implementation

Our proposed decoupled formulation involves solving (16) for each of the subsystems. As already mentioned, many numerical tools already exist for solving (16); we will use the implementation in [15]. For the examples in this paper, we used the numerical schemes below.

For the numerical Hamiltonian in (16), we used the Lax-Friedrich approximation [26]. For the numerical spatial derivatives , we used a fifth-order accurate weighted essentially non-oscillatory scheme [26, 27]. For numerical time derivatives , we used a third-order accurate total variation diminishing Runge-Kutta scheme [27, 28].

Computations were done on a computer with an Intel Core i7-2600K CPU running at 3.4 GHz, with 16 GB of memory.

Iv 4D Quadrotor Collision Avoidance

Consider a simple quadrotor model consisting of two decoupled double-integrators:


denote the - and -position of the quadrotor, and denote the - and -velocity. The control signals are the - and -acceleration of the quadrotor, constrained to be between and .

Now, consider two quadrotors in a pursuit-evasion game, in which the evader (Player 1) aims to avoid collision, while the pursuer (Player 2) aims to cause a collision. The relative coordinates of the two quadrotors are given by the following state variables:


Given the above relative state variables, the relative dynamics of the two quadrotors are given by


Note that this system is decoupled, with as the first decoupled component, and as the second decoupled component. In the relative coordinates of the two quadrotors, we define the collision set of size , representing the configurations in which the two quadrotors are considered to have collided, as the following set:


with the corresponding implicit surface function where . Since we have a decoupled system, let be the following sets:


with corresponding implicit surface functions . Then, we have and .

We will set as the target set in our reachability problem, and compute the backwards reachable set from using three methods:

  • Solve (8) directly in to obtain , whose zero sublevel set represents .

  • Solve (16) in , , to obtain using our proposed decoupled formulation described in Section III-B.

  • Compute the analytic boundary of the reachable set.

For comparison purposes, for the first two methods we will compute on the computation domain . However, it is important to recall that for our proposed method described in Section III-B, we can significantly reduce space complexity by only storing a small part of .

Iv-a Reachable Set

Since the state space of our system is 4D, we visualize various slices of the reachable set, whose boundary is given by . Figures 1 shows these slices. The reachable set boundary computed using our proposed decoupled method is very close to the reachable set boundary computed by solving the full PDE (8) in and to analytic reachable set boundary. Figures 2 zooms in on the plots for a closer look.

Fig. 1: Various slices of the reachable set.
Fig. 2: Various slices of the reachable set, zoomed in.

Iv-B Performance

In order to quantify the computation error, we converted into a signed distance function from the boundary . This operation was first proposed in [29] and can be done by solving the reinitialization PDE formulated in [30]; for this operation, we use the implementation in [15]. We then evaluated approximately 24 million analytically-computed reachable set boundary points on the ; the resulting values represent how far each of the analytically-computed points are from the numerically-computed boundary. The values of on analytic boundary points are defined as the computation error.

Figure 3 shows the error as a function of grid spacing. In terms of the maximum error (red curve), the decoupled formulation results in a numerically-computed reachable set boundary that is accurate within the size of the grid spacing (black line). On average, the error is approximately an order of magnitude smaller than the size of grid spacing (blue curve). Furthermore, we can see numerical convergence to the analytic solution as the grid spacing size decreases.

Figure 4 shows the computation time as a function of the number of grid points in each dimension. Here, we can see that the decoupled formulation is orders of magnitude faster than the full formulation, and can be done with many more grid points in each dimension. Lastly, the slopes of curves in the log-log plot show an time complexity for the full formulation, and only for the decoupled formulation.

For the decoupled formulation, when we reconstruct the full value function in 4D (blue curve), the computation time hardly increases compared to when we do not perform full reconstruction (green curve). However, in general, we recommend that the value function in only a region near a state of interest should be computed. Without full reconstruction of the value function, we are able to obtain results with many more grid points (green curve), improving the accuracy of the numerical computation.

Fig. 3: Mean and maximum error of the reachable set computed using the decoupled formulation, as a function of the grid spacing.
Fig. 4: Computation time as a function of the number of grid points in each dimension.

V 6D Quadrotor Collision Avoidance

Consider relative dynamics augmented by the velocity of evader quadrotor, given in Equation (27). These dynamics are needed to impose a velocity limit on the quadrotor.


For this system, we consider a collision between the two quadrotors, as defined previously, to be unsafe configurations. Here, we denote this set


with corresponding implicit surface function .

We also consider configurations in which Player 1 is exceeding a velocity limit of in the - or - directions to be unsafe. This set of configurations is denoted :


We define the target set to be the union of the above configurations; represents all unsafe configurations. We will compute the reachable set from the target set as follows:

  1. Compute , the reachable set from target set .

  2. Compute , the reachable set from target set .

  3. Take the union to obtain .

V-a Reachable set

Taking the union , we obtain the 6D reachable set. To visualize , we compute 2D slices of the 6D reachable set at various values. This is done without computing the entire 6D reachable set by setting the computation domain to be in a large portion of the plane, at a small range of values. The 2D slices at shown in Figure 5. Each subplot shows two different pairs for a particular . The red boundary represents the slice with , while the blue boundary represents a slice with the evader velocity almost exceeding the limit of . The blue boundaries contain the red ones, because if the evader is already near the velocity limit, it would have more limited capability to avoid collisions.

Fig. 5: 2D slices of the 6D reachable set for the augmented relative dynamics of two quadrotors.

Figure 6 shows a simulation of the collision avoidance maneuver resulting from the reachable set. The red and blue quadrotors are initially traveling in opposite directions. Consider the situation in which the red quadrotor insists on staying on its intended path. In this case, the blue quadrotor must perform the optimal avoidance control whenever it reaches the boundary of the reachable set (), shown as the blue dotted boundary. When the blue quadrotor is no longer at the boundary of the reachable set (), it is free to perform any control.

Fig. 6: A simulation of collision avoidance between the two quadrotors using the 6D reachable set.

Vi Conclusions and Future Work

We have presented a decoupled formulation of the HJ reachability, enabling us to efficiently compute the viscosity solution of the HJ PDE (8), which gives the reachable set for optimal control problems and differential games. When the system dynamics are decoupled, our decoupled formulation allows the exact reconstruction of the solution to the full dimensional HJ PDE by solving the HJ PDEs corresponding to each of the decoupled components. Our novel approach enables the analysis of otherwise intractable problems. In addition, our formulation achieves this without sacrificing any optimality compared to the full formulation. We demonstrated the benefits of our approach using a 4D and 6D quadrotor system.


  1. A function between two measurable spaces and is said to be measurable if the preimage of a measurable set in is a measurable set in , that is: , with -algebras on ,.


  1. I. Mitchell, A. Bayen, and C. Tomlin, “A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games,” IEEE Transactions on Automatic Control, vol. 50, no. 7, pp. 947–957, July 2005.
  2. C. Tomlin, J. Lygeros, and S. Shankar Sastry, “A game theoretic approach to controller design for hybrid systems,” Proceedings of the IEEE, vol. 88, no. 7, pp. 949 –970, Jul 2000.
  3. T. Basar and G. Olsder, Dynamic Noncooperative Game Theory, 2nd ed.   Philadelphia, PA: SIAM, 1999.
  4. J. Lygeros, C. Tomlin, and S. Sastry, “Controllers for reachability specifications for hybrid systems,” Automatica, vol. 35, no. 3, pp. 349 – 370, 1999.
  5. Department of the Air Force, “United states air force unmanned aircraft systems flight plan 2009-2047,” oai.dtic.mil, Jan 2009.
  6. A. Madrigal, “Autonomous robots invade retail warehouses,” http://www.wired.com/wiredscience/2009/01/retailrobots/.
  7. H. Erzberger, “Automated conflict resolution for air traffic control,” 25 th International Congress of the Aeronautical Sciences, Jan 2006.
  8. B. P. Tice, “Unmanned aerial vehicles – the force multiplier of the 1990s,” Airpower Journal, Spring 1991.
  9. J. F. Fisac, M. Chen, C. J. Tomlin, and S. S. Sastry, “Reach-Avoid Problems with Time-Varying Dynamics, Targets and Constraints.” in 18th International Conference on Hybrid Systems: Computation and Controls, 2015.
  10. O. Bokanowski, N. Forcadel, and H. Zidani, “Reachability and minimal times for state constrained nonlinear problems without any controllability assumption,” SIAM Journal on Control and …, pp. 1–24, 2010.
  11. E. Barron, “Differential Games with Maximum Cost,” Nonlinear analysis: Theory, methods & applications, pp. 971–989, 1990.
  12. J. Ding, J. Sprinkle, S. S. Sastry, and C. J. Tomlin, “Reachability calculations for automated aerial refueling,” in Proceedings of the 47th IEEE Conference on Decision and Control, Cancun, Mexico, 2008.
  13. M. Chen, Z. Zhou, and C. Tomlin, “Multiplayer reach-avoid games,” in American Control Conference 2014, 2014.
  14. H. Huang, J. Ding, W. Zhang, and C. Tomlin, “A differential game approach to planning in adversarial scenarios: A case study on capture-the-flag,” in Robotics and Automation (ICRA), 2011 IEEE International Conference on, 2011, pp. 1451–1456.
  15. I. Mitchell, A Toolbox of Level Set Methods, 2009, http://people.cs.ubc.ca/~mitchell/ToolboxLS/index.html.
  16. I. M. Mitchell, “The flexible, extensible and efficient toolbox of level set methods,” Journal of Scientific Computing, vol. 35, no. 2-3.
  17. S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces.   Springer-Verlag, 2002, ISBN: 978-0-387-95482-0.
  18. J. A. Sethian, “A fast marching level set method for monotonically advancing fronts,” Proceedings of the National Academy of Sciences, vol. 93, no. 4, pp. 1591–1595, 1996. [Online]. Available: http://www.pnas.org/content/93/4/1591.abstract
  19. I. M. Mitchell, “Scalable calculation of reach sets and tubes for nonlinear systems with terminal integrators: A mixed implicit explicit formulation,” in Proceedings of the 14th International Conference on Hybrid Systems: Computation and Control, ser. HSCC ’11.   New York, NY, USA: ACM, 2011, pp. 103–112. [Online]. Available: http://doi.acm.org/10.1145/1967701.1967718
  20. I. M. Mitchell and C. J. Tomlin, “Overapproximating reachable sets by hamilton-jacobi projections,” Journal of Scientific Computing, vol. 19, no. 1-3, pp. 323–346, 2003.
  21. J. McGrew, J. How, L. Bush, B. Williams, and N. Roy, “Air combat strategy using approximate dynamic programming,” AIAA Guidance, Navigation, and Control Conference, Aug 2008.
  22. J. B. Lasserre, D. Henrion, C. Prieur, and E. Trélat, “Nonlinear optimal control via occupation measures and lmi-relaxations,” SIAM J. Control Optim., vol. 47, no. 4, pp. 1643–1666, Jun. 2008. [Online]. Available: http://dx.doi.org/10.1137/070685051
  23. L. C. Evans and P. E. Souganidis, “Differential games and representation formulas for solutions of Hamilton-Jacobi-Isaacs equations,” Indiana Univ. Math. J., vol. 33, no. 5, pp. 773–797, 1984.
  24. P. Varaiya, “On the existence of solutions to a differential game,” SIAM Journal on Control, vol. 5, no. 1, pp. 153–162, 1967.
  25. M. G. Crandall and P.-L. Lions, “Viscosity solutions of Hamilton-Jacobi equations,” Transactions of the American Mathematical Society, vol. 277, no. 1, pp. 1–42, 1983.
  26. S. Osher and C.-W. Shu, “High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations,” SIAM Journal on Numerical Analysis, vol. 28, no. 4, pp. pp. 907–922, 1991. [Online]. Available: http://www.jstor.org/stable/2157779
  27. S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces.   Springer Verlag, 2003.
  28. C.-W. Shu and S. Osher, “Efficient implementation of essentially non-oscillatory shock-capturing schemes,” Journal of Computational Physics, vol. 77, no. 2, pp. 439 – 471, 1988. [Online]. Available: http://www.sciencedirect.com/science/article/pii/0021999188901775
  29. D. L. Chopp, “Computing minimal surfaces via level set curvature flow,” Journal of Computational Physics, vol. 106, no. 1, pp. 77 – 91, 1993.
  30. M. Sussman, E. Fatemi, P. Smereka, and S. Osher, “An improved level set method for incompressible two-phase flows,” Computers & Fluids, vol. 27, no. 5–6, pp. 663 – 680, 1998.
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 minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

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 description