First Order Hadamard Variation of the Harmonic Navigation Function on a Sphere World

,


Introduction
The aviation industry currently faces meaningful challenges to overcome the increasing traffic across the world.Safety solutions are a demand in all flight phases in order to cope with traffic capacity, all-weather conditions and efficiency.In order to address the 5% per year increase in air traffic, the Single European Sky ATM Research (SESAR) program in Europe and the NextGen program in the USA have been initiated to design new rules and tools for future air traffic management (ATM).
Each flight phase has specific operational constraints.Of interest here is the en-route flight phase.This phase is comprised from completion of initial climb through cruise altitude and completion of controlled descent to the initial approach fix.En-route air traffic is currently managed by subdividing airspace into sectors and air routes.In each sector, a team of air traffic controllers supervises the transiting aircraft and prevents conflicts by deviating some aircraft from their planned route.Each team communicates with the teams responsible for neighbouring sectors to transfer aircraft flying through several sectors.
As traffic grows, the workload in a sector may exceed air traffic controller capabilities.In this case, the sector is subdivided into smaller sectors, thus lowering the number of aircraft under the responsibility of an air traffic controller at any given time.However, this increases the time spent transferring aircraft between sectors up to the point where coordinating between sectors takes precedence over separating aircraft and solving conflicts.The SESAR program overcomes this pitfall by delegating most of the separation task to the aircraft themselves, thus inducing a higher level of automation.In this context, fixed air routes will no longer be necessary and trajectories will be planned in 4D in such a way that conflicts are avoided by design.
A major challenge associated with aircraft trajectory planning, in addition to operational constraints and flight dynamics, is the influence of wind.Much research has been conducted to plan such routes several days in advance.However, no algorithm is known to produce planning such that the resulting aircraft trajectories be robust to uncertain head-or tailwind.Artificial potential fields, as first used by [14], enable a mobile to reach a destination while avoiding obstacles in a complex environment.This method has been extended to dynamic environments, with multiple mobile trajectory planning as in [15], coordinating between mobiles in [31] and moving obstacles in [8].Navigation functions, a specific case of artificial potential fields, introduced originally for robots trajectory planning by [6] and [25], were extended to generate aircraft trajectories.Such trajectories meet ATM restrictions on speed and curvature.Furthermore, a proof of conflict avoidance can be obtained.Navigation functions have also been used to plan trajectories under a wind constraint in [23] for sailboat navigation.Recent developments in navigation functions also include their application to the case of real-time spacecraft guidance [26], or in uncertain environments [9].However, the navigation function works under the assumption that aircraft or other mobiles move in a deterministic way.We present in this paper a method to account for time uncertainty in aircraft trajectory planning.To account for this uncertainty, we consider the variation of the navigation function with respect to variations of the boundary of the configuration space.Of interest in this paper is the practical construction of Hadamard's variational formula, as presented in [10] for the Green's function of the Laplace equation on a 2-dimensional sphere world.
The purpose of this paper is to present a method for computing the Hadamard variational formula given a specific domain perturbation.This variational formula requires knowledge of the Green's function on the initial domain.Numerical evaluations of the Green's function for certain problems, such as the Helmholtz operator on periodic structures [13], heat diffusion in 1D [1], the elliptic problem [2], or the exterior Neumann problem [28].
An analytical expression of the Green's function for Laplace's equation has also been obtained in multiply connected domaines [3].However, this formulation is difficult to implement in practice, as it relies on the Schottky-Klein prime function [4].
The motivation for this paper is from designing schemes to plan trajectories for multiple mobiles under uncertainty using harmonic navigation functions.
After a state of the art on navigation functions in section 2.1, section 2.2 presents a method for numerically computing the Green's function of the Laplacian operator, and section 2.3 applies this method to the computation of the variation of the navigation function with respect to the domain boundary variation.The methods are illustrated and the results for this method are exhibited in section 3. Finally, future work is outlined in section 4.

Navigation functions
Navigation functions on an arbitrary sphere world were introduced by Rimon and Koditschek [16] in the context of robot navigation for generating trajectories that present a guarantee of obstacle avoidance.Navigation function are a family of potential functions, with their maximum on the boundary of the obstacles, their minimum at the destination and no interior local minimum or maximum.Formally, navigation functions are defined as follows: 1. Analytic on the interior of F 2. Polar on F with a minimum at an interior point q d ∈ •F 3. Admissible on F

Morse on Ω
Once the navigation function is computed, robot trajectories are determined by following −∇ f at non critical points or the steepest descent given by the Hessian at critical points.The trajectories hereby obtained reach the destination while avoiding obstacles.
Navigation functions have been extended to aircraft trajectory planning [27], taking into account ATM considerations such as limits on speed and curvature, as aircraft can only fly at [−3%, +6%] of their nominal speed [5].
Of present interest are harmonic navigation functions on a manifold of R 2 with non-overlapping spherical boundaries.Such a domain is referred to as a sphere world.Let D be a sphere world as described and z d the destination point in the interior of D.
Consider the Dirichlet problem ∆φ = 0 on the interior of D where is the Laplacian operator in 2 dimensions.As stated in [24], the function φ on L 2 (D) defined by this problem is a navigation function.Indeed, no interior point can be a local extrema so that φ is polar.φ is admissible by construction and is Morse on D However, navigation functions work under the assumption that trajectories are deterministic.
The kind of trajectory obtained by this method is not robust against uncertainties, and thus phenomenon such as wind may cause the separation norms to be violated.This prevents navigation functions from being deployed in an operational context.The purpose of this work is to provide tools for extending navigation functions to the case of aircraft trajectories under a wind uncertainty.
In the following sections, we will introduce the Green's function and a numerical method to obtain it on sphere worlds, and its application to navigation function under a domain boundary uncertainty.

2.2.
Semi-analytical approximation of the Green's function for the Laplacian operator

Green's function
Let Ω ∈ R 2 be a planar bounded domain with the boundary ∂Ω = i C i such that C i be non intersecting circles and C 1 be exterior to all other circles, as illustrated in Figure 1.The centers of circles C i are noted c i and the radii r i .In this configuration, ∂Ω is Lipschitz-continuous and Ω is a Lipschitz domain.
Let us consider the equation where f ∈ L 2 (Ω) is a given function.If the domain boundary ∂Ω is sufficiently smooth, then Equation where δ s is the Dirac distribution centered in s ∈ Ω. Contrary to Equation 3, the right hand term in Equation 4 is not L 2 .However, the L p -integrability of the Green's function for elliptical equations on Lipschitz domains is given by [19].Furthermore, as presented in [21], it follows from Equation 4that the solution u to Equation 3 has the form

Green's function from conformal mappings
Definition 2. Let Ω ⊂ C be a domain.A fundamental solution of the Laplace equation is a mapping A representation theorem for such a G is given in [11]: Let G be a fundamental solution of Laplace equation in a domain Ω.Then it exists a function h : Ω × Ω → R, harmonic in its first argument, and a real constant λ > 0 such that for any (z, For a given λ, h, hence G, is uniquely defined by its values at the boundary ∂Ω.The special choice λ = 1/2π and vanishing boundary conditions yields the so-called Green's function of the Laplace equation.It can be obtained by solving for h the system: Theorem 4. Let G be Green's function for the Laplace equation in the domain Ω.The following relation holds in the distributional sense for any z 0 ∈ Ω: Proof.Using the standard operators ∂, ∂, it comes: ∆ z G(., z 0 ) = 4∂∂G(., z 0 ).For an arbitrary function φ holomorphic in Ω, continuous on ∂Ω: By Stokes theorem: Since dx ∧ dy = −2idz ∧ dz, the equation 10 yields: and finally using the relation: the right hand side term equals: Using the system defining h to obtain an approximate solution is not an easy task from a numerical point of view as it requires solving a Laplace equation that becomes increasingly ill-conditioned as z 0 approaches the boundary of the domain.This is mainly due to the fact that the singularity lies within the domain of interest.
In order to transform this inconvenient geometry into a more convenient one, we choose a conformal mapping that eliminates the singularity at z 0 .
Consider the following conformal mapping Provided that z 0 is interior to Ω, the map of Ω by T z 0 noted Ω T , as illustrated in Figure 2, is the unbounded region exterior to a set of non-overlapping circles of center c T j and of radii r T j defined by . Given that T z 0 is a conformal mapping, H is the harmonic solution exterior to non-overlapping disks that vanishes on the disk boundaries and such that H(z) ∼ log |z| as |z| → ∞.Furthermore, contrary to G(•, z 0 ), function H has no singularity on its domain.Function G(•, z 0 ) can then be determined given the knowledge of function In the following subsection we present a method for approximating H

Harmonic function exterior to non-overlapping disks
We seek a function H on domain Ω T as defined prior.We use the following Laurent expansion [17] for Mapping of the sphere world Ω by T z 0 .Ω T is the space exterior to three circles. with and where α, β j , γ jk and δ jk are scalars to be determined.This function is harmonic with H(z) ∼ log |z| as |z| → ∞ and vanishing values on the boundary of Ω T .
This series ∑ ∞ k=1 (γ jk − iδ jk )(z − c T j ) −k can be approximated by its N first terms.Thus we consider the approximation h of function The Green's function G(•, z 0 ) can then be approximated by the function g = h • T −1 z 0 .This function g has the additional benefit of being easily derived to any order as it is the composition of standard functions.Given Equation 18 the derivatives of h are ∀z ∈ ΩT , ∂h ∂x The derivatives of g are thus obtained by combining Equations 19, 15 and 12. Once the coefficients for h are determined, these derivatives come for free, as no additional coefficients need to be computed.
This property is of interest to us for the purpose of the Hadamard variation covered in section 2.3 as the normal derivatives on the domain boundary are required.
The values of coefficients of the expansion of h are obtained by a numerical method derived from that of [30], using collocation points equally spaced on the set of boundary circles, to approach the coefficients of 18 given a domain Ω.A least squares method can be used to set the boundary conditions.To ensure that the vanishing condition on the boundary is satisfied, a number n of collocation points denoted by z i , i ∈ [ [1, n]] are set on each circle.Values for the coefficients are chosen by minimizing the sum of squares of h(z i ).Finding the coefficients of h is then equivalent to solving the optimization problem min ||AX|| 2 s.t.v T X = 1 (20) where X and v are two vectors of length (1 and A is an n by (1 + J + 2N J) matrix. 129 The problem set in Equation 20 can be turned into a set of linear equations by using Lagrange multipliers [18].Given enough collocation points, this problem has a unique solution X as defined in Equation 21, from which we obtain the approximation h from Equation 18of the harmonic function on Ω T and the approximation g of the Green's function, with function g defined for a fixed point z 0 as

First order variation of the harmonic solution 130
The purpose of this section is to present the main results relating the behavior of the solution of an elliptic problem with Dirichlet boundary conditions to the variation of its border.Let (M, g) be a connected orientable smooth Riemannian manifold of dimension d and let γ i , i = 1 . . .N be a finite set of disjoint embedded smooth submanifolds of dimension d − 1 that partition M into two disjoint components M + , M − that are manifolds of dimension d with common boundary γ = ∪ i=1...N γ i .M + is assumed to be relatively compact.Finally, let X be a smooth vector field on M with flow φ(t, x) that is transversal to γ (except perhaps at a finite number of points).The flow of X will be denoted by φ : ] − , [×M → M, with as usual ∂ t φ(0, x) = X(x).For a fixed t, φ t will denote the function φ t : x → φ(t, x).The flow φ will model an admissible deformation of the manifold with boundary M + ∪ γ.Following [22], the deformed boundary at t ∈] − , [ will be defined as γ t = φ (t, γ). and the perturbed problem as: Since the primary goal is to obtain a formula for the Green's function variation, the problem stated in Equation 24 will be reformulated to accommodate a distribution in the right hand term of Equation 24.Likewise, boundary condition 25 will be simplified later to u t | γ t = 0.The final problem in distributional form is then: where ω is the riemanian volume form and T is the right hand term distribution.Recalling that ∆u t = div∇u t , it comes: where the subscript N denotes the component of the vector normal to γ and σ = i N ω, with N the normal vector to γ.Any test function v is the image by φ * t of a test function φ −1 * t v.The perturbed problem 27 can thus be rewritten so as to involve only integrals on M + : Taking the derivative of the first integral in the left hand side of 28 with respect to t at t = 0 yields: By Cartan's formula, div(X)ω = d(i X ω), thus: Putting it in 29, the derivative becomes: The derivative with respect to t at t = 0 of the boundary term: can be expressed as: with the boundary condition u| γ = 0, the first term vanishes.However, it will be kept in the sequel to obtain a more general variation formula.Gathering things together yields: Considering the problem: and putting back its solution as v in 33 and using: gives for the derivative of the solution: Finally, if u is taken to be the green function of the original problem, the well known Hadamard formula is obtained [29]:

Numerical validation of the Green's function approximation
In this section, we discuss the results obtained by implementing the above method and verify that the resulting functions present certain expected properties.We also justify the decisions regarding the choice of parameters for Equation 18in the light of the relative error and the computation time.We finally discuss the relevance of this method for the application at hand, namely computing navigation functions.
Python 2.7 was used to implement the method described in section 2.2 on an office computer with an Intel Core i7-4710MQ CPU , 2.50GHz, 8 cores.

Convergence of the method
Reconstruction of the harmonic solution for constant boundary conditions Given the Green's function G on a domain Ω × Ω, the harmonic solution can be found for any Dirichlet boundary conditions from Equation 5.In particular, one important property of the Green function is that if the function u in Equation 5 is constant and equal to 1, then we find that for any x ∈ Ω ∂Ω ∂G ∂n (x, y) dy = 1 (37) We build function g to estimate the Green's function as defined in Equation 23 on a domain featuring two inner circle boundaries of centers 2i and −1 − 2i and of radii 1. and 0.5 respectively, and bounded by an outer circle of radius 5. centered at the origin.We then evaluate the quality of the estimation by comparing the integral of its normal derivative to 1 : The integral is computed using the quad function from the scipy Python library [12], and based on a technique from the Fortran library QUADPACK.
In Figure 3 is a measure of the quality of the estimation from Equation 37, given by the logarithm of ε, relative to the number of collocation points at a given point x = −2 interior to the domain and for a given number of terms in the Laurent expansion N = 2, N = 5 or N = 10.Overall, this value follows a decreasing curve, thus the distance decreases as the number of collocation points increases.
There is a steep increase of the curve for values of n c around 2N.For higher fixed values of N, the curve starts at a higher point, and converges for higher values of n c , but in either case, the distance from ∂Ω ∂g ∂n (x, y) dy to 1 plateau around 10 −14 for n c > 2N.
This indicates that the obtained values are constrained by round-off errors.
There is no improvement for n c > 2N, since the problem is over-constrained, as can be seen in the linear problem stated in Equation 20, where 1 + J + 2(N − 1)J coefficients are to be determined from Jn c equations.
In the following section, we will further evaluate the quality of the approximation of the Green's function derived from our method by integrating over the domain boundary to find the harmonic solution for non-constant boundary conditions, similar to those of the foreseen application of our method.

Reconstruction for non-constant boundary conditions
On certain sphere worlds, an analytical solution can be found to Equation 3. Here, we consider the crown defined by the unit circle as an outer boundary, and the circle of center 2/5 and of radius 2/5 as an inner boundary as pictured in Figure 4.
For a domain bounded by two concentric circles, with constant Dirichlet conditions on each circle, the harmonic solution to this problem is a function of the form z → A + B log z with constants A and B chosen to fit the boundary condition.A Möbius mapping [20] can be found that transforms this symmetrical problem into that described in Figure 4. From here, and given Dirichlet boundary conditions such that the value on the inner boundary be constant equal to scalar a ∈ R and the value on the outer boundary be constant equal to b ∈ R, then the harmonic potential on this domain is function We can thus compare the solution reconstructed with the Green's function as described in Equation 5, noted u g The z i in this equation are placed on a fine enough Cartesian grid, such that the result has converged.
In Figure 5a, the logarithm of the distance between both solutions is presented.As the number of terms used in the Laurent series increases, the difference 1 decreases for high numbers of collocation points.As discussed previously, the error converges for values of n c greater than 2N.As such, as the number of terms N in the Laurent expansion increases, the error 1 converges for larger numbers n c of collocation points, but converges towards a lower value.
In Figure 5b, log 1 is represented for the comparison between the analytical solution and that obtained with the finite elements method using Lagrange elements of order 2. Here, we observe that the error for the solution computed with the Green's function is comparable to that for the finite elements method.
Although we compare the results obtained using the reconstruction from the Green's function to those of the finite element method, it is to be noted that reconstructing the harmonic solution from the Green's function is more time-consuming than solving the finite elements problem, as the Green's function simultaneously gives the harmonic solution for all boundary conditions.As such, our method is not appropriate for solving single navigation functions.However, in the context of solving boundary variation problems, the hereby presented method computing the Green's function is superior.We further illustrate this in the following paragraphs by exhibiting the time complexity.

Time complexity
In this section, we consider the execution time to compute the approximation of the Green's function at one point, noted g in section 2.2.We normalize the resulting times relative to the execution time in the test case where N = 2 and n c = 2.The results are presented in Figure 6 for different values of N and of n c .
As can be seen in Figure 6, the time required to compute the Green's function rapidly increases as the number of terms in the Laurent expansion increases.However, as seen in sections 3.1.1and 3.1.1,the solution rapidly converges and good results can be obtained with small values of N.
The order of magnitude for the base case in our configuration, with n c = 2 and N = 2, is a millisecond.This relatively small amount of time must be put in perspective with the fact that the approximation g of the Green's function may be computed hundreds of times for different points of the domain in order to estimate the value of the harmonic solution in one point.
In the following section, we will present a practical example of our method, relevant to the application at hand.In each of the Laurent expansion, ten terms are used.Thirty collocation points are set on the domain boundaries at regular intervals.
Sphere worlds with more obstacles can be tackled with no additional difficulty.

Computing navigation functions
In order to demonstrate the usability of our method in the context of aircraft trajectory planning, we have carried out the computation of navigation functions with the Green's function and traced some resulting trajectories.These are exhibited in Figure 8.The domain is defined with dimensions similar to those found in operational contexts.The outer boundary, figuring the limit of the airspace is 100NM in diameter ; the obstacle discs, representing the other aircraft are 5NM in diameter.The destination is near the center of the domain.
The trajectories are plotted using a gradient descent method [7] with a constant step.We find that the generated trajectories all converge towards the destination, while avoiding all obstacles with a wide margin, as is illustrated in Figure 8.

Numerical validation of the Hadamard variation of the Green's function
In the present section, we will exhibit the results obtained with the aforementioned method for variations of the domain boundary.We consider here a sphere world Ω bounded by one destination disc (c 0 , r 0 ), one obstacle disc (c 1 , r 1 ) and an outer boundary disc (c 2 , r 2 ), such as the one presented in Figure 7.We hereby study two manners of disturbing the domain boundary • by varying the center of the obstacle disc, • and by varying the radius of the obstacle disc.where n(w) is the normal vector to the domain boundary at point w ∈ ∂Ω.For a pure translation of the disc, the term δr is null.For a pure variation of the obstacle disc radius, the term δc is null.
In practice, for determining aircraft trajectories, only the center of the obstacle disc is submitted to uncertainty, since the radius obstacle disc is defined by the aircraft separation norms.
In the following paragraphs, we will consider the variation of the Green's function for small variations of the obstacle radius.
We first introduce notations.Let G : Ω × Ω be the Green's function computed at each η ∈ Ω from Equation 23.Given the study of parameters in section 3.1, the values of N and n c in computing the Green's function are set respectively to 10 and 30.δG denotes the variation of the Green's function as defined in Equation 36 and the Green's function estimated using the Hadamard variation is noted G = G + δG.The Green's function on the perturbed domain Ω * is noted G * : Ω * × Ω * and is similarly computed at each η ∈ Ω * from Equation 23.Finally, δG, the Hadamard variation of the Green's function is compared to G − G * , the actual variation of the Green's function on the entire domain Ω × Ω using the relative error for the variation as a metric where the norm || • || 1 is computed on Ω × Ω by discretizing the domain along a Cartesian grid as was done for Equation 40.This error, which is expected to approach 0 for small values of δρ, is given in Figure 9 for perturbations of decreasing amplitudes, for a dilation of the radius of the inner disc.

Conclusion
In this article, a method for effectively computing the navigation function when the domain boundary is uncertain was presented.It is a two-fold method, where section 2.2 aims at determining the Green's function for the Laplacian operator and section 2.3 uses the resulting Green's function for determining the Hadamard variation.This promising method can be extended to the case of stochastic partial differential equations for aircraft trajectory planning.

c 1 c 4 c 2 c 3 Figure 1 .
Figure 1.A sphere world bounded by an outer circle of center c 1 and three non overlapping inner circles of centers c 2 , c 3 and c 4 respectively.

10 ( 1 N = 10 (
a) log 1 as a function of n c for a fixed number of terms in the Laurent expansion N = 2, N = 5 or N = b) log 1 as a function of the number of mesh vertices for the finite element method with Lagrangian elements of order 2.

Figure 5 . 10 Figure 6 .
Figure 5. log 1 with our method based on the Green's function compared to that obtained with the finite elements method

Figure 7
Figure 7 depicts examples of sphere worlds with one obstacle and the singularity placed at different points in the domain.The resulting function vanishes except in the vicinity of the singularity.

( a )Figure 7 .
Figure 7. Approximation g of the Green's function on unit disc with obstacle of radius 0.1 centered at −0.2 − 0.5i and destination disc of radius 0.2 centered at 0.5i