Approximation of Reachable Sets using Optimal Control Algorithms

We investigate and analyze a computational method for the approximation of reachable sets for 
 nonlinear dynamic systems. The method uses grids to cover the region of 
 interest and the distance function to the reachable set 
 evaluated at grid points. A convergence analysis is provided and shows the convergence of three 
 different types of discrete set approximations to the reachable set. 
 The distance functions can be computed numerically by suitable optimal 
 control problems in combination with direct discretization techniques 
 which allows adaptive calculations of reachable sets. 
 Several numerical examples with nonconvex reachable sets are presented.


Introduction
The subject of this paper is the analysis and extension of an algorithm in [5] for the approximation of reachable sets of non-linear control problems.
Reachable sets appear in various applications, e.g. in the presence of disturbances of parameters in control problems, in estimates in terminal points of all solutions of a control problem, differential inclusions, and differential games.Specific applications can be found in population models, fish harvesting, collision avoidance, weather forecasts, climate models, space orbit calculations, and many others, see [12], [17], [37].
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.Various methods for the approximation of reachable sets for linear control problems have been suggested, among them are set-valued integration schemes [2], optimal control techniques [52,47,4], external and inner ellipsoidal techniques [36,37,38], estimation methods [23,34,21], see also references listed therein.
However, in the non-linear case less methods are known (for an overview see [12]) as the reachable set is non-convex.Häckl in [29] used time discretization combined with ε-grids in state space (see also [46]), Chahma [12] used set-valued Runge-Kutta methods for nonlinear problems with state constraints, in [6] the Euler's method is studied with a detailed analysis of the discretization error in state space.
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, where the optimal value yields the distance of the grid point to the reachable set.The corresponding optimal control problems are not solved theoretically by use of the Pontryagin's maximum principle as in [52] but numerically by suitable discretization methods.The resulting DFOG method (DFOG = distance field on grids) 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 [42,19,30] and the references stated therein.
A similar approach using optimal control techniques was discussed in [4] in the special case of linear control problems.Herein, the optimal value of the optimal control problem provides the support function exploiting the convexity of the reachable set.In the nonlinear case this approach is not applicable anymore, as the reachable set is in general non-convex.Hence, we switch to the distance function which allow the description of non-convex sets, see [13,Sect. 1].
We also would like to mention an interesting similar approach in [9] which uses distance functions in the objective function of optimal control problems and a path-following idea for boundary points of the viability kernel.The distance function penalizes infeasible states or states outside an a priori chosen bounding box.In the DFOG method the objective function is the squared Euclidean distance and the optimal value function is the distance function.The approach in [9] creates nonsmooth optimization problems which are solved by simulated annealing techniques.The use of distance functions allows a quicker computation of viability kernels in comparison to other approaches.
In the calculation of the value function, which plays an important role in level-set methods, a similar approach with optimal control problems is applied in [37,39] to derive the Hamilton-Jacobi equations.Level sets of the value function then determine the reachable set.
Further approaches can be found in [17] (simulation of trajectories with piecewise constant control functions by Runge-Kutta methods), approximation schemes using Volterra series in [33] or methods based on zonotopes in [26] for hybrid linear systems.
The paper is organized as follows.In Section 1 the problem of computing the reachable set of a non-linear control problem is introduced.Section 2 provides three discrete approximation strategies for closed sets which are rather elementary.Their approximation accuracy depends on the gridwidth in state space discretization and is also discussed.The strategies depend Lipschitz w.r.t. the set which is discretized.Further properties and connections are also included in this part.In Section 3 properties and the approximation of reachable sets are recalled to lay the basis for the overall estimates of the DFOG method.Section 4 discusses the numerical implementation of the three strategies from Section 2. Here, the grid point is the parameter entering the optimal control problems which leads to the calculation of projections to the (unknown) discrete reachable set.Using the results of Section 3, error estimates w.r.t.stepsize and state-space grids are presented.In the DFOG method the popular direct discretization of optimal control problems is combined with an approximation strategy for sets based on distance functions and best approximations.Numerical examples for various non-linear control problems are presented in Section 5. Finally, concluding remarks and an outlook is given in Section 6.
Let t 0 < T be given and let U = ∅ be a convex and compact subset of R m .Moreover, let an initial state x 0 ∈ R n be given.Consider the following nonlinear control problem.
The task is to compute the reachable set at time T which is defined as follows: and ∃ corresponding solution x(•) of Control Problem 1.1 with x(T ) = y Consider a suitable one-step discretization scheme, e.g. an explicit Runge-Kutta method, with increment function Φ on a time grid with time points t i = t 0 + ih, i = 0, 1, . . ., N, and stepsize h = (T − t 0 )/N.For simplicity the grid is chosen equidistantly.Then, using the discretization scheme, a discrete counterpart of the continuous control problem is defined as follows.
Problem 1.2 For a discretized control function u h (•) ∈ U h , i.e. u h : I → U and U h being a finite dimensional approximation of U, find a solution x h (•) with In the simplest case we choose Φ(t, x, u, h) = f (t, x, u) and U h as the piecewise constant functions with values in U, that is Euler's method.This method is highly studied in the approximation of control problems, see [19] for an overview.
An approximation of the continuous reachable set R(T, t 0 , x 0 ) is given by the discrete reachable set defined by and 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 based on appropriate models, changes in climate can be studied for different environmental influence factors, like carbondioxid concentrations or temperature, see [12,Sect. 5.3].
Another field of applications are robust control problems.Herein, a control problem with uncertain dynamics is considered: where P denotes an appropriate parameter region.Let u * be a given control law (e.g. an optimal control for a fixed p * (•) ∈ P ) and let x(u * , p)(•) denote the solution for any p(•) ∈ P .The task is to decide whether u * robustly obeys given constraints, e.g.whether c(t, x(u * , p)(t), u * (t)) ≤ 0 ∀p(•) ∈ P holds.This constraint can be checked, if the reachable set of x for a fixed u * and for varying p(•) ∈ P can be approximated.The following lemma leads to a characterization of (parts of) a set by the complement of open balls whose radii are given by the distance function.
Lemma 2.2 Let S ⊂ R n be closed and nonempty.Then, where r(x) = dist(x, S).
Proof: The representation follows immediately from [13, § 1, Corollary 6.2].Please notice that the case ŝ = x, i.e. x ∈ S, can be safely ignored in the union above, since r(x) = x − ŝ = 0 and the interior of the ball is empty.
The Hausdorff distance can be expressed not only by the Minkowski duality, but more general for closed sets by the distance function, which immediately results in regularity properties for the distance function.
Proposition 2.3 Let S, S ⊂ R n be closed and nonempty.Then,

Inner/Outer Approximation of Sets
If not otherwise specified, let ρ • Z n be a grid with grid size ρ > 0. For a given set S ⊂ R n we denote by S ρ the slightly enlarged set S restricted to the grid, i.e.
The strategy (2) to extend a set depending on the distance of grid points in one coordinate before restricting it to the grid is commonly used in set-valued analysis, especially in viability theory, see e.g.[50,48], as well as [16,12,6].We explain the need to extend sets by an example that Lars Grüne provided to us, compare also [28,Section 2.3].In Figure 1 the dark colored set S (a small ball with an attached horizontal line segment with irrational x 2 -coordinate) is not well approximated by its grid projection S ∩ ρZ n for rational figures ρ.Only one element of S, the origin, is hit by grid points (marked with black diamonds) for the values ρ = 2, 1, 1 2 .Thus, the Hausdorff distance remains constant and no point of the line segment with irrational coordinate is ever hit by grid points.In contrary, Figure 2 shows that the extension S ρ of S (dark colored set) yields a better approximation of S. Now, 2, 3 or 5 grid points lie in the dark colored shaded extension for ρ = 2, 1, 1  2 respectively.Here, the Hausdorff distance between S and its (shrinking) extension S ρ is proportional to ρ.
The following lemma suggests that this extension provides a first strategy to approximate a closed set via grid points.Although the result is wellknown, we provide the proof for the convenience of the reader.Lemma 2.4 Assume that S is closed and nonempty, ρ > 0, and ρ • Z n be the (infinite) grid.Then, the distance of a point s in S to a grid point can be estimated by Denote the first strategy by then S ρ = M 1 ρ (S) and the following estimate holds: We define the grid point g ρ ∈ ρZ n componentwise with the Gauss bracket ⌊•⌋ for the biggest integer not exceeding a real number: Since dist(g ρ , S) can be estimated by g ρ − s , we found a bound for the one-sided Hausdorff distance d(S, M 1 ρ (S)).On the other hand, choosing so that the estimate follows.The equality of the first strategy with the extension S ρ is clear by the equivalence of dist(g ρ , S) A second elementary strategy is described in the following lemma and is based on best approximations of grid points.
Lemma 2.5 Assume that S is closed and nonempty, ρ > 0 and ρ • Z n be the (infinite) grid.Let us denote the second strategy by Then, Proof: Clearly, d(M 2 ρ (S), S) = 0, since all best approximations are elements of S. On the other hand, for any g ρ ∈ M 1 ρ (S) there exists s ∈ S with by Lemma 2.4.For this grid point g ρ let us choose ŝρ ∈ Π S (g ρ ).Then, The following lemma provides an approximate set representation via the complement using only finitely many open balls instead of infinitely many ones as in Lemma 2.2.Lemma 2.6 Let S ⊂ R n be closed and nonempty, ρ > 0 and ρ • Z n be the (infinite) grid.Define the third strategy by with r(x) = dist(x, S).Then, obviously holds.Each v ∈ M 3 ρ (S) cannot lie in the interior of the ball B r(gρ) (g ρ ) for all g ρ ∈ ρZ n due to (5).Hence, the inequality v − g ρ ≥ dist( g ρ , S) holds for is valid by Lemma 2.4.Hence, there exists s ∈ S with Surprisingly, all above estimates of the three discretization strategies M i ρ (S), i = 1, 2, 3, do not depend on the set S and the regularity of its boundary.
Let us specialize these results to compact sets S to be able to restrict our discretization to a finite set G ρ of grid points.Although we change the set of grid points slightly, we stick to the old notation M i ρ (S) to avoid additional notation.
Corollary 2.7 Let S ⊂ R n be compact and nonempty with a (compact) bounding box G ⊂ R n , i.e. S ⊂ G.For a given grid size ρ > 0 we use the extension G ρ ⊂ ρZ n in (2) with finitely many grid points.If we replace ρ • Z n by G ρ in the three strategies, the following holds: Proof: Let us check the three strategies separately, since only small modifications of the proofs of Lemmas 2.4-2.6 are necessary.
(i) If we take s ∈ S, it is obvious that the constructed grid point 2 ρ.On the other hand, each grid point g ρ in G ρ with distance to S not exceeding 2 ρB 1 (0) so that the rest of the proof remains unchanged.(ii) We only need to adapt the second part of the proof slightly.The best approximation g ρ of an element s ∈ S to the set M1 ρ (S) has distance smaller or equal to √ n 2 ρ.Now, the same reasoning as in (i) shows that 2 ρB 1 (0) which means that g ρ is an element of G ρ .The proof can be finished as in Lemma 2.5.
(iii) Similar to (ii) one shows that g ρ ∈ G ρ for a chosen v ∈ M 3 ρ (S) in Lemma 2.6.The rest of the proof remains unchanged.
The next proposition compares the three strategies.The first strategy (restricted to grid points in S) is included in the second one 1 , the latter is a subset of the third strategy2 and this strategy is contained in the original set S. The first two strategies consist only of finitely many points, the third one will later allow an adaptive modification of the approximation.Proposition 2.8 Let S ⊂ R n be closed and nonempty and let ρ > 0 be the grid size.Then, If S is compact and G is a compact bounding box and we substitute ρ • Z n by G ρ in all three strategies, we will have the same inclusions.
Proof: A grid point g ρ that lies in S ∩ M 1 ρ (S) coincides with its best approximation in S and thus is contained in the second strategy.
Take a best approximation ŝρ ∈ S of a grid point g ρ ∈ ρ • Z n .Let us assume that there exists g ρ ∈ ρ • Z n with ŝρ − g ρ < dist( g ρ , S).This grid point g ρ could not lie in S, since we would get the contradiction ŝρ − g ρ < 0. Its best approximation s ∈ Π S ( g ρ ) fulfills which would be a contradiction to the optimality of ŝρ .Hence, such a grid point g ρ cannot exist and ŝρ is an element of the set of the third strategy.
The last inclusion is already shown in the proof of Lemma 2.6, while the reasoning is the same for the three strategies using finitely many points in G ρ instead of infinitely many in ρ • Z n .
We end this section with a perturbation result for the three discretization strategies.Except a fixed error of term O(ρ) each of the three strategies behaves Lipschitz continuous w.r.t. the change in the set.Furthermore, all three strategies differ only in Hausdorff distance by O(ρ).Lemma 2.9 Let S, S ⊂ R n be closed and nonempty and ρ > 0 be the grid size.Then, for all strategies i, j = 1, 2, 3, where c 1 = 1 and c 2 = c 3 = 2.
Proof: The first estimate follows immediately from Lemmas 2.4-2.6,since ) .The second one is clear by the estimate The application of these three strategies to reachable sets is obvious.Since we know that the discrete reachable sets are close to the continuous one under appropriate assumptions, we will apply the three discretization strategies to the discrete reachable set and we will still have a good approximation of the reachable set in Problem 1.1.
Then, the closure of the reachable set for the right-hand side F equals the reachable set for the convexified (relaxed) right-hand side co F .
In our problem setting we need to ensure that the reachable set is closed and nonempty, for the numerical computations the compactness is essential.Both is guaranteed by the following proposition.
be continuous w.r.t.all arguments, continuously differentiable w.r.t.x and U : I ⇒ R m be continuous with compact, nonempty images.Consider F : I × R n ⇒ R n defined as Assume that F has convex images and is of linear growth, i.e.
Then, the reachable set is compact and nonempty.
Proof: To apply [31, Theorem 20.1] we only have to show that But this follows from the linear growth condition: The previous result refers to a parametrization of the right-hand side of the differential inclusion.This formulation is convenient for the verification of the assumptions for the examples in Section 5.Alternatively, it would be possible to state a similar result with the set-valued right-hand side.In [31,Corollary 20.2], the relaxation theorem as in Proposition 3.1 is formulated for a parametrization.
For our analysis we formulate a central result establishing the convergence of reachable sets with slightly more restrictive assumptions on the right-hand side based on the results in [18], see also [55].
Then, the following convergence result holds for all N ∈ N, h = T −t 0 N : Proof: cf.[18, Section 1, Theorem]

Discrete Approximation of Reachable Sets
Under the assumptions of Proposition 3.2, the reachable set is compact and nonempty, hence we may apply the results of Section 2 to derive discrete approximations of the reachable sets.Nevertheless, the convergence result stated in Proposition 3.3 discusses only the discretization in time.To implement set-valued Euler's method we need the additional discretization in space, see [12, Section 4.2] and [6].In contrast to the cited articles, the choice of the grid points in the three discretization strategies in Lemmas 2.4-2.6 depend on values of the distance function or on closest points.Summarizing we obtain the following main result, where we formulate the abstract convergence result not only for the Euler method, but more generally for set-valued Runge-Kutta methods, see [54,53,12,3,4].Actually, any discretization method of order p could be treated by the theorem.Theorem 3.4 Let U ⊂ R m be convex, compact, nonempty and let f : Let R = R(T, t 0 , X 0 ) be the compact, nonempty reachable set of the differential inclusion, h = T −t 0 N be a stepsize with N ∈ N, let R h = R h (T, t 0 , X 0 ) be the discrete reachable set of the set-valued Runge-Kutta method which is assumed to be closed and nonempty.Assume that the choice of the set-valued method and U as well as the regularity of the parametrization guarantee that Then, taking ρ = C 2 • h p for the three discretization strategies in (3)-( 5) Proof: We just apply the triangle inequality for the Hausdorff distance and combine the convergence assumption (7) with the results of Corollary 2.7, since the reachable sets R h are bounded and the assumption yields compactness and nonemptiness.In [54,3] regularity assumptions for Runge-Kutta methods of order 2 applied to linear differential inclusions are formulated, in [53] also for the nonlinear ones with strongly convex right-hand sides.Typical conditions involve smoothness of the parametrizing function f (•, •, •) but also for the support functions (t, x) → δ * (l, F (t, x)) uniformly in normed directions l.
In the next corollary we restrict ourselves to the set-valued Euler's method.
Corollary 3.5 Consider the notations for U, R as in Theorem 3.4 and let f : Then, the reachable sets R and R h are compact, nonempty and the results of Theorem 3.4 hold with convergence order p = 1.
Proof: Everything follows from Theorem 3.4, if we can show that the assumption of the parametrization fits to the set-valued convergence results.Under the given assumptions F (•) is Lipschitz w.r.t.(t, x) and we can apply Propositions 3.1-3.2.In particular, we can drop the assumptions of the convex images in Proposition 3.2, see also [51].Due to this relaxation theorem, the Hausdorff distance remains unchanged and this proposition shows that the reachable set R is compact, nonempty.It is not difficult to conclude that this is valid for the discrete reachable sets R h , too.Proposition 3.3 clarifies that the assumption (7) holds.Hence, we can apply the theorem above.

Numerical Realization
We suggest a numerical method which allows to approximate reachable sets for non-linear problems using optimal control techniques.The theoretical approximation properties using a combination of distance functions and statespace grids have been derived in Theorem 3.4.Although it is possible to consider other set-valued Runge-Kutta methods, we restrict ourselves to the simplest case, the set-valued Euler's method in [18,55,12,6].

DFOG Method
The algorithm is presented only for the set-valued Euler's method, but can easily be adapted to set-valued Runge-Kutta methods of higher order in the light of Theorem 3.4.
The algorithm works with a grid G h with stepsize h and projects each element in G h onto the reachable set of the dynamic system.For simplicity, we choose C 2 = 1 in adapting the grid size ρ = C 2 h to the stepsize h in time.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.The algorithm is called DFOG method, since it applies distance fields3 on grids.G by its grid extension G h in (2) with stepsize ρ = h (or ρ = h p for higher-order methods).
(ii) For every g h ∈ G h solve the following optimal control problem: Let x ⋆ (•; g h ) and u ⋆ (•; g h ) denote the solution of OCP(g h ).
(iii) Define the reachable set approximation (relative to G h ) according to one of the three strategies (3)-( 5) stated in Lemmas 2.4-2.6.
Please notice the following relations in optimal control problem OCP (g h ): • The admissible set is the set of solutions of Problem 1.1.
• The set of all endpoints x(T ) of admissible solutions forms the reachable set R = R(T, t 0 , X 0 ).
• The optimal value coincides with 1 2 • dist(g h , R) 2 .• The endpoint x ⋆ (T ; g h ) is a best approximation of the grid point g h in the reachable set R.
Similar relations hold also for the following discretized version of OCP(g h ) and Problem 1.2, provided the same discretization scheme is used.
In DOCP(g h ), u h is a suitable control discretization.For simplicity, we restrict ourselves to Euler's method so that Φ(t, x, u, h) = f (t, x, u) and u h will be the following piecewise constant control approximation on the grid Obviously, DOCP(g h ) is an approximation of OCP(g h ) and any global solution of OCP(g h ) is an element of Π R (g h ) with R = R(T, t 0 , x 0 ) and any global solution of DOCP(g h ) computes an element of Π R h (g h ) with R h = R h (T, t 0 , x 0 ).The convergence results in the previous section guarantee that the reachable set is approximated with at least order h for Euler's method.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 subarcs, etc. Particularly, they may possess local minima, which as we shall see later may cause problems in combination with an adaptive strategy.
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.
Let z := (u 0 , u 1 , . . ., u N −1 ) ⊤ .Then the constraints u h (t i ) ∈ U read as In order to reduce the number of variables of DOCP(g h ) the equations can be eliminated recursively according to Herein, we identified the grid function u h (•) with the control parameterization z.Using these expressions, an equivalent optimization problem with variable z arises: 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 well-analyzed, see for instance the book of Wright and Nocedal [45], we are not going into details here.All these methods have in common that they require the gradient of the objective function, which is the most costly operation during the numerical solution.Hence, it is important to exploit the structure of Problem 4.2.
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 the discretized optimal control problem.Details can be found in Bock [7], Caracotsios and Stewart [11], Maly and Petzold [43].A discussion and comparison of several strategies can be found in Feehery et al. [20].
• The adjoint ODE approach, see Cao et al. [10], 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.
• Finite difference approximations are easy to implement but tend to be costly and it is difficult to control the accuracy of the computed gradients.
Since Problem 4.2 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 where Following [25] let us 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 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 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 Φ(t, x, u, h) = f (t, x, u), but it is more involved for more general Runge Kutta methods.

Remark 4.3
• The direct discretization approach outlined above can be easily extended to more complicated control constraints.Even state constraints and boundary conditions can be added, see [5].However, the calculcation 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 Gerdts [25].
• The effort for solving the discretized optimal control problems obviously increases as the stepsize h decreases.
• 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.Thus, all calculations in Section 5 may contain inaccuracies owing to local minimality.Nevertheless, the obtained numerical results are in good correspondence to the results in [12].

Numerical Examples
In all following pictures, grid points g h ∈ h•Z n in red color indicate a negative status of the optimizer for the discrete optimization problem DOCP(g h ).In this case, the corresponding grid point (or its best approximation, depending on the strategy that is used) will be colored in red.Furthermore, we use Euler's method as direct discretization method, the grid width in state space equal to the stepsize in time and an optimizer OCPID-DAE1 [24].It is surprising that the algorithm produces nice results, although it only returns a local minimum.Note that local minimizers still define reachable points.We observed that the lack of global optimality is often cured by considering many grid points.

Kenderov's Example
This example was suggested by Petar Kenderov (see [12,Example 5.2.1]).It is constructed in such a way that the reachable set is a sphere, 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)) 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 sphere.Figure 3 shows the numerical results for a discretization of N = 20, 40, 80, 160, 320 and the linear convergence w.r.t.ρ = h.
Please notice that the axes are different for the first two pictures in the first line of Figure 3.
For this example, the results of the three strategies differ in an obvious way, see Figure 4. Without the extension in the first strategy (a), only two grid points lie on the reachable set, whereas a major part of the reachable set appears with extension (b).The best result of the four pictures is (c), the second strategy.Since the best approximations of grid points do not need to be grid points themselves, we can drop the extension here.Due to the onedimensionality, the third strategy (d) produces reasonable theoretical results, but the visualization is not helpful.

Bilinear Example
This example contains a bilinear term in the dynamics which are given below: Linear convergence and the numerical approximations of the reachable set for T = 1 and N = 10, 20, 40, 80, 160 are visualized in Figure 7.For this example, the results of the three strategies are rather similar, see Figure 8.The extension in the first strategy (b) does not improve significantly the reachable set and is very similar to the non-extended version (a).In Figures 8(c   It turns out that the extra effort for checking the balls in the third strategy does not lead to improvements.
Figure 9 shows an impression of the funnel of trajectories by plotting reachable sets for t =

Adaptive Version
The idea of an adaptive algorithm for the third strategy is as follows.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 and the ball could not contain an element of the reachable set as r is supposed to be minimal.Following a heuristic approach, all grid points within the ball B r (g h ) will not to be projected.
In Figure 11(a) the bounding box is shown in which all grid points are chosen.For each grid point the best approximation to the (unknown) reachable set (painted in blue color) is computed.The reachable set belongs to the Rayleigh problem in [12, Ex. 5.2.5], [5,Example 2].Since neighboring grid points create similar balls, a first idea of an adaptive algorithm is rather obvious.We do not compute projections of grid points that already lie inside an open ball calculated for another grid point, in contrary to Figure 11(b).In this way, the calculation time can be drastically reduced, but the approximation in state space is coarser.

Figure 11: non-adaptive algorithm: grid points lie in many balls
There is one pitfall.The above reasoning assumes that the numerical method is able to find a global minimizer of DOCP(g h ).However, in practise this cannot be guaranteed, if a local optimization method is used.If the method finds a local minimizer but not a global minimizer, then the radius r of the ball B r (g h ) is too large and too many grid points are cut off.A remedy for this problem would be to use global optimization methods, which unfortunately is not practical owing to the large dimension of the discretized optimal control problems.In our numerical examples we observe that the above mentioned shortcoming is cured by considering many grid points and storing all end-points of locally optimal trajectories, which might not be globally optimal but still are reachable.This pitfall will be demonstrated for Kenderov's example.This example is particularly well-suited for demonstration, because the reachable set is a set of measure zero.
The benefit in view of CPU time of the adaptive strategy depends on the initial region and on the dimension of the reachable set.Thus, the adaptive strategy is particularly efficient for low dimensional reachable sets, e.g.Kenderovs example below.
CPU times for Kenderov's problem are summarized in Table 1.The adaptive algorithm turns out to be very efficient in view of CPU, because of the low dimension of the reachable set.However, for this example, special measures have to be taken to avoid that the adaptive algorithm cuts off some regions of the reachable set, because the optimization algorithm only found local minimizers instead of global ones.Let us point out that the computation times in Tables 1 and 2 use the strategy ρ = h N .This is required by the convergence result in Theorem 3.4, but the discretization in space in Figures 7,12 is much finer than optically necessary.An obvious way to drastically reduce computational time is the choice of a rougher state space grid.
In Figure 12 we will compare the non-adaptive and the adaptive version of the third strategy.In the latter one, only those grid points, which do not lie in other open balls from other grid points, determine open balls.2. The adaptive algorithm turns out not to be so efficient in view of CPU time as for Kenderov's example, because the reachable set is two dimensional and the initial region already localizes the reachable set sufficiently well.

Outline
The proposed algorithm has the following potential advantages and extensions: • The approximation of reachable sets with higher order methods (Runge-Kutta methods of second order in [54,53]) is possible with adapted stepsize ρ = C 2 h 2 .
• Zooming into interesting sub-regions of the reachable set can be easily realized.
• DFOG requires to solve one optimal control problem for each point in the O(h)-grid in state space.In common implementations of the set-valued Euler method X h (t j+1 ) = x h (t j )∈X h (t j ) {x h (t j ) + hF (t j , x h (t j ))} (j = 0, 1, . . ., N − 1) as in [12] and [6], the sets within the union have to be approximated with an O(h 2 )-grid in state space to keep the overall order O(h) for subsequent times t k , k > j.
The effort for solving the finite dimensional optimization problem for O(N n ) grid points and exhaustive calculation of set-valued Euler steps from neighboring grid points with approximately the same images has to be compared.Our experience shows advantages of the DFOG method in computational time and memory consumption.
• The grid points do not need to be chosen in an equidistant fashion and adaptivity based on the third strategy is possible.
• The algorithm can be easily parallelized.
• Higher-dimensional systems can be treated as long as the discretized optimal control problems can be solved.
• State constraints, mixed control state constraints, and terminal conditions can be introduced in the algorithm.
• It remains a future a task to adapt the method to the regularity of the boundary of the reachable set, see [41].

Figure 1 :
Figure 1: bad approximation of the dark colored set without extension

Figure 2 :
Figure 2: better approximation of the light colored set with extension (dark colored set)

1
Properties and Approximations of Reachable Sets Proposition 3.1 Let the set-valued mapping F : I × R n ⇒ R n be locally Lipschitz w.r.t.x with integrable Lipschitz bounds L R (•) on balls B R (0) and have closed, nonempty images.Assume that F is of linear growth with integrable bound C(•), i.e.

Algorithm 4 . 1 (
DFOG method)(i) Choose a bounding region G ⊆ R n for the reachable set and approximate

Figure 5 :
Figure 5: Kenderov example: two zooms for third strategy for N = 20

Figure 6 :
Figure 6: Kenderov example: unchecked and checked balls in third strategy for N = 20

Figure 9 :
Figure 9: bilinear problem: evolution of reachable sets in time.

Figure 10 :
Figure 10: bilinear problem: evolution of reachable sets in time.

(a) N = 10 (Figure 12 :
Figure 12: non-adaptive and adaptive second strategy for the bilinear problem with T = 0.5 and N = 10, 40, 160 If needed, we use the notation π S (x) for one element of Π S (x).

Table 1 :
CPU times for Kenderov's problem: Comparison of non-adaptive and adaptive algorithm.

Table 2 :
CPU times for the bilinear problem: Comparison of non-adaptive and adaptive algorithm.