A computational method for non-convex reachable sets using optimal control

A computational method for the approximation of reachable sets for non-linear dynamic systems is suggested. The method is based on a discretization of the interesting region and a projection onto grid points. The projections require to solve optimal control problems which are solved by a direct discretization approach. These optimal control problems allow a flexible formulation and it is possible to add non-linear state and/or control constraints and boundary conditions to the dynamic system. Numerical results for non-convex reachable sets are presented. Possible applications include robust optimal control problems.


I. INTRODUCTION
The subject of this paper is the description of an algorithm for the approximation of reachable sets of non-linear control problems. Reachable sets are interesting because they allow to study the future development of dynamic systems under the influence of control variations and variations in parameters. For instance, changes in climate can be studied for appropriate models using different environmental influence factors like carbondioxid concentrations or global mean temperature, see [9]. Further examples on fish harvesting, lake pollution, and spruce-budworm control can be found in [8]. The idea of our approach is to project grid points from an equidistant grid onto the reachable set. Each projection requires to solve an optimal control problem. The occurring optimal control problems are not solved theoretically by use of the Pontryagin's maximum principle as in [23] but numerically by suitable discretization methods. This approach turns out to be powerful in practice and allows to include control and/or state constraints and even boundary conditions. Results concerning the convergence of discretized optimal control problems can be found in [20], [10], [16] and the references stated therein.
Let t 0 < T be given and let U = ∅ be convex and compact. Moreover, let an initial state x 0 ∈ R n be given. Consider the following nonlinear control problem.
Problem 1: For a given u ∈ L ∞ ([t 0 , T ], R m ) find x ∈ W 1,∞ ([t 0 , T ], R n ) with a.e. in [t 0 , T ], a.e. in [t 0 , T ]. The task is to compute the reachable set at time T which is defined as follows: y ∈ R n | ∃u(·) control function and ∃x(·) corresponding solution of Problem 1 with x(T ) = y Many properties of the reachable set are known for linear control problems with f being linear in x and u. Most importantly, it can be shown that the reachable set for linear control problems is a convex set as a consequence of Aumann's integral. Various methods for the approximation of reachable sets for linear control problems have been suggested, among them are set-valued integration schemes [1], optimal control techniques [23], [2], external and inner ellipsoidal techniques [17], [18], [19], estimation methods [12].
However, in the non-linear case, very few methods are known. In general, the reachable set is non-convex. Chahma [9] used set-valued discretization methods for nonlinear problems with state constraints. The analysis of error estimates and the numerics are continued in [4].
A discrete counterpart of Problem 1 is constructed as follows. Consider a suitable one-step discretization scheme, e.g. an explicit Runge-Kutta method, with increment function Φ on an equidistant time grid with time points t i = t 0 + ih, i = 0, 1, . . . , N , and step size h = (T − t 0 )/N . Then, using an appropriate parameterization u h of the control on the grid, e.g. a piecewise constant control function, a discrete counterpart of the continuous control problem is defined as follows.
Problem 2: For a discretized control function u h (·) find a solution x h (·) with In the simplest case we have Φ(x, u, h) = f (x, u), that is Euler's method. Of course, higher order integration schemes can be used as well.
An approximation of the continuous reachable set R(T, t 0 , x 0 ) is given by the discrete reachable set defined by R h (t N , t 0 , x 0 ) := y ∈ R n | ∃u h (·) discretized control and ∃x h (·) corresponding solution of Problem 2 with x h (t N ) = y .

II. THE ALGORITHM
We suggest a computational method which allows to approximate the reachable set of a nonlinear dynamic system at a fixed time point T by using optimal control techniques.
The algorithm works with a grid G h with step size h and projects each element in G h onto the reachable set of the dynamic system. Projecting a grid point w.r.t. the Euclidean norm leads to an optimal control problem and the following algorithm for the approximation of the reachable set.
Algorithm: (i) Choose a region G ⊆ R n and discretize G by a grid G h ⊂ G with step-size h such that each element of G can be approximated by a grid point with error h. (ii) For every g h ∈ G h solve the following optimal control problem OCP (g h ): Remark: There are alternative ways to define a reachable set approximation in (iii), for instance one could use the union of all grid points which yield distance (approximately) zero in the objective function of the respective optimal control problems. Or one could approximate the reachable set by the complement of the union of balls centered at the grid points whose radii are defined by the objective function values of the above optimal control problems. In the upcoming paper [3] theoretical approximation properties of the above algorithm are investigated showing that the reachable set will be approximated of at least order O(h) for Euler's method under suitable assumptions. Herein, it is important to synchronize the grid size to the time discretization. Higher order approximations are possible as well.
The same discretization scheme as defined earlier in Problem 2 leads to a discretized version of OCP(g h ) called DOCP (g h ): Let x h (·; g h ) and u h (·; g h ) denote the solution of DOCP(g h ).
In DOCP(g h ) u h is a suitable control discretization. For simplicity, u h will be a piecewise constant control approximation on the grid, i.e.
Obviously, DOCP(g h ) is an approximation of OCP(g h ) and any global solution of OCP(g h ) is an element of It remains to solve DOCP(g h ) (or ideally OCP(g h )). OCP(g h ) and its discrete counterpart DOCP(g h ) are in general non-convex problems and may exhibit all difficulties that may occur for general (discretized) optimal control problems like ill-conditioning, non-regularity, singular sub-arcs, etc. Particularly, they may possess local minima, which may cause problems as this may result in not detecting parts of the actual reachable set.
In order to make DOCP(g h ) accessible to numerical methods, we assume that the control set U is defined by box constraints, i.e.
In order to reduce the number of variables of DOCP(g h ) the equations can be eliminated recursively according to Herein, we identified u h with the control parameterization z.
Using these expressions, an equivalent optimization problem with variable z arises: Problem 3: This is a bound constraint nonlinear program and it can be solved, for instance, by a sequential quadratic programming (SQP) method or any other suitable nonlinear programming method. As all these methods are well-known and wellanalyzed, see for instance the book of Wright and Nocedal [22], we are not going into details here. All these methods have in common that they require the gradient of the objective function. Computing the gradient is the most costly operation during the numerical solution and it is important to exploit the structure of Problem 3.
There are basically four approaches for calculating derivatives: • The sensitivity ODE approach (also known as IND approach) is a sensitivity analysis of the integration scheme w.r.t. z. As the effort depends mainly on the number of variables and less on the number of constraints, it is particularly efficient if the number of nonlinear constraints exceeds the number of variables in a discretized optimal control problem. Details can be found in [5], [7], [21]. A discussion and comparison of several strategies can be found in [11]. • The adjoint ODE approach, see [6], is advantageous compared to the sensitivity ODE approach if the number of nonlinear constraints is less than the number of variables in the discretized optimal control problem. The effort mainly depends on the number of constraints and less on the number of variables. • Algorithmic differentiation, see [15], combines the sensitivity ODE approach (forward mode) and the adjoint ODE approach (backward mode). • Finite difference approximations are easy to implement but tend to be costly and it is difficult to control the accuracy of the computed gradients owing to round-off errors.
As Problem 3 only has box constraints, the adjoint approach for calculating gradients is the most efficient one. As we shall see, it only requires to integrate the differential equation from t 0 to T and the adjoint equation backwards from T to t 0 . We intend to calculate the gradient w.r.t. z of i = 0, 1, . . . , N − 1.

Consider the auxiliary functional
with multipliers λ i , i = 1, . . . , N . Differentiating J w.r.t. z, reordering terms, and neglecting arguments yields Herein, Φ x [t i ] is an abbreviation for Φ x (X i (z), z, h) and likewise for Φ z .
The calculation of the terms X i (z) is costly and shall be avoided. Hence, terms involving X i (z) have to be eliminated. This leads to the adjoint equation and the transversality condition λ N = −(X N (z) − g h ) . Notice, that the adjoint equation is solved backward in time. Exploiting these relations yields It is straightforward to show that G (z) = J (z) holds, see [14, Thm. 6.2.2], and thus we obtained a formula for the gradient of G.
Notice, that the derivatives Φ x and Φ z have to be computed. This is straightforward for Euler's method with Φ(x, u, h) = f (x, u), but it is more involved for more general Runge Kutta methods.
The overall integration scheme with (1) and (2) turns out to be a symplectic integration scheme, see [13]. Remarks: • The direct discretization approach outlined above can be easily extended to more complicated control constraints. Even state constraints and boundary conditions can be added. However, the calculation of gradients using the adjoint approach may not be the most suitable one if non-linear control and/or state constraints are present in the optimal control problem. In this case the sensitivity approach is preferable, details can be found in [14]. • The effort for solving the discretized optimal control problems and the number of grid points increases as the step-size h decreases. Please note that the computation of reachable sets is inherently a challenging task and this is reflected by the above comment. • Common nonlinear programming methods are only capable of finding a local minimizer of the above optimization problem. Global optimality is practically not achievable with a reasonable computational effort. However, approximation properties for the reachable set can only be guaranteed for global solutions. Local solutions may cause the approximate reachable set to miss out parts of the actual reachable set. The computational results in the next section however suggest that even local optimization methods usually lead to very good approximations. This effect seems to be compensated by considering many grid points.

III. NUMERICAL EXAMPLES
The true reachable sets in this section are either known analytically (Kenderov example) or reference solutions using exhaustive calculations are available, see [9]. Using these reference solutions, it turns out that the above optimal control strategy actually finds correct approximations to the reachable set without missing out parts of it, despite the fact that only local optimization algorithms were used. In addition to the standard algorithm outlined above, an adaptive version was created. The basic idea of the adaptive algorithm is based on the following lemma.
Lemma 1: Let g h be a grid point and x (T ; g h ) an optimal solution of DOCP(g h ). Then no grid point within the ball B r (g h ) and radius r = x h (T ; g h ) − g h is reachable.
Proof: As r is supposed to be minimal, no grid point within the ball B r (g h ) is reachable.
Lemma 1 implies that the grid points within the ball B r (g h ) need not to be projected. This simple strategy turned out to be useful in order to reduce the computational effort.

A. Example 1: Brachistochrone
The first example is the well-known Brachistochrone example posed in the 17th century by Johann Bernoulli. The corresponding control problem reads as follows.
x (t) = 2gy(t) cos(u(t)) y (t) = 2gy(t) sin(u(t)) x(0) = 0 The numerical computations reveal that the reachable set is convex although the control problem itself is non-convex. Figure 1 shows the numerical results for a discretization of N = 5, 10, 20, 40: The numerical solution of this problem turned out to be nasty because the root function in the differential equations is not defined for negative arguments. To avoid trajectories with negative y component, an additional state constraint y(t) ≥ 0 was introduced. The resulting optimal control problems became more difficult to solve and more advanced methods than previously mentioned are required.

B. Example 2: Rayleigh Problem
The second example originates from an optimal control problem of an electric circuit and is known as the Rayleigh problem. The nonlinear control problem reads as follows.
CPU times for the Rayleigh problem can be found in Table I. The third example was suggested by Petar Kenderov. It is constructed in such a way that the reachable set is a circle, that is the reachable set is a set of measure zero. The nonlinear control problem reads as follows.
x (t) = 8 (a 11 x(t) + a 12 y(t) − 2a 12 y(t)u(t)) y (t) = 8 (−a 12 x(t) + a 11 y(t) + 2a 12 x(t)u(t)) Herein, a 11 = σ 2 − 1, a 12 = σ √ 1 − σ 2 , and σ = 0.9. The numerical computations reveal that the reachable set is non-convex and the approximations apparently converge to a circle. Figure 3 shows the numerical results for a discretization of N = 20, 40, 80, 160, 320: CPU times for Kenderov's problem are summarized in Table II. The adaptive algorithm turns out to be very efficient in view of CPU time because the low dimension of the reachable set. However, for this example, the adaptive algorithm causes some regions of the reachable set to be cut-off because the optimization algorithm only found local minimizers instead of global ones.

A. Domains of Attraction
There is also a close relationship between reachable sets and domains of attraction for dynamical systems, that is the domain of attraction of a given point x T at time t 0 is just the reachable set R(t 0 , T, x T ) of the backward dynamic system.

B. Robust Control
Another field of applications are robust control problems. Consider a control problem with uncertain dynamics where P denotes an appropriate parameter region within a suitable function space. Let u * be a given control law, for instance an optimal control for a fixed p * (·) ∈ P . Let x(u * , p)(·) denote the solution for any p(·) ∈ P . The task in robust control is to decide whether u * robustly obeys given constraints, e.g. whether c(x(u * , p)(t), u * (t)) ≤ 0 holds for any p(·) ∈ P . If the reachable set of x for a fixed u * and for varying p(·) ∈ P can be approximated for any time t, then the validity of the above constraint can be checked by inserting the reachable set, i.e.
A similar task can be defined for feedback systems. Suppose the control u is fixed by a nonlinear feedback law u = g(x). Introducing this relation into the dynamic system yields x (t) = f (x(t), g(x(t)), p(t)), p(·) ∈ P.
The solution operator x(p)(·) and the reachable set w.r.t. p allows to verify whether inequality constraints of typẽ c(x(p)(t)) ≤ 0 are fulfilled or not.

V. CONCLUSIONS AND FUTURE WORKS
A computational method for the approximation of reachable sets for non-linear control systems was presented. For a number of examples the method showed its capability of computing the reachable set. The results are comparable with the ones computed by a parallelized version of set-valued Euler's method in [9]. Differently to the approach in [9], which needs intermediate state grids of step-size O(h 2 ) at each of the N Euler steps, this approach only needs one grid of step-size O(h). Another advantage of the method is the flexibility with respect to adding control and state constraints and boundary conditions.
A thorough analysis of theoretical approximation properties of the proposed algorithm has to be completed. Moreover, the effect of local optimization methods has to be studied in more detail, in particular in combination with adaptive algorithms as local solutions of the projection problems may cause parts of the reachable set to be cut-off. Furthermore, it would be desirable to not only approximate the reachable set at a given time point T but instead for a whole time interval [t 0 , T ]. It has to be investigated whether the solutions of the optimal control problems, which live on [t 0 , T ], can be used to at least approximate the reachable sets at intermediate time points in [t 0 , T ].