{'fraction_non_alphanumeric': 0.04933706923943499, 'fraction_numerical': 0.016147202403304545, 'mean_word_length': 4.0072317037894125, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Trajectories of a UAV Base Station Using Hamilton-Jacobi Equations

We consider the problem of optimizing the trajectory of an Unmanned Aerial Vehicle (UAV). Assuming a traffic intensity map of users to be served, the UAV must travel from a given initial location to a final position within a given duration and serves the traffic on its way. The problem consists in finding the optimal trajectory that minimizes a certain cost depending on the velocity and on the amount of served traffic. We formulate the problem using the framework of Lagrangian mechanics. We derive closed-form formulas for the optimal trajectory when the traffic intensity is quadratic (single-phase) using Hamilton-Jacobi equations. When the traffic intensity is bi-phase, i.e. made of two quadratics, we provide necessary conditions of optimality that allow us to propose a gradient-based algorithm and a new algorithm based on the linear control properties of the quadratic model. These two solutions are of very low complexity because they rely on fast convergence numerical schemes and closed form formulas. These two approaches return a trajectory satisfying the necessary conditions of optimality. At last, we propose a data processing procedure based on a modified K-means algorithm to derive a bi-phase model and an optimal trajectory simulation from real traffic data.


I. INTRODUCTION
Unmanned Aerial Vehicles (UAV) are expected to play an increasing role in future wireless networks 1 [1]. UAVs may be deployed in an ad hoc manner when the traditional cellular infrastructure is missing. They can serve as relays to reach distant users outside the coverage of wireless networks.
They also may be used to disseminate data to ground stations or collect information from sensors. In this paper, we address one of the envisioned use cases for UAV-aided wireless communications, which 1 relates to cellular network offloading in highly crowded areas [1]. More specifically, we focus on the path planning problem or trajectory optimization problem that consists in finding an optimal path for a UAV Base Station (BS) that minimizes a certain cost depending on the velocity and on the amount of served traffic. Our approach is based on the Lagrangian mechanics framework and the use of Hamilton-Jacobi (HJ) equations.

A. Related Work
UAV trajectory optimization for networks has been tackled maybe for the first time in [2]. The model consists in a UAV flying over a sensor network from which it has to collect some data. The UAV can learn from previous experience, which is not assumed in our study. The problem of optimally deploying UAV BSs to serve traffic demand has been addressed in the literature by considering static UAVs BSs or relays, see e.g. [3], [4]. The goal is to optimally position the UAV so as to maximize the data rate with ground stations or the number of served users. In these works, the notion of trajectory is either ignored or restricted to be circular or linear. In robotics and autonomous systems, trajectory optimization is known as path planning. In this field, there are classical methods like Cell Decomposition, Potential Field Method, Probabilistic Road Map, or heuristic approaches, e.g. bio-inspired algorithms [5]. Authors of [6] have capitalized on this literature and proposed a path planning algorithm for drone BSs based on A* algorithm.
The main goal of these papers is to reach a destination while avoiding obstacles. In our work, we intend to minimize a certain cost function along the trajectory by controlling the velocity of the UAV. This goal is studied in optimal control theory [7] and is applied for example in aircraft trajectory planning [8]. Most numerical methods in control theory can be classified in direct and indirect methods. In direct methods, the problem is transformed in a non linear programming problem using discretized time, locations and controls. Direct methods are heavily applied in a series of recent publications in the field of UAV-aided communications, see e.g. [9]- [12]. Formulated problems are usually non-convex. The standard approach is hence to rely on Successive Convex Approximation (SCA), which iteratively minimizes a sequence of approximate convex functions. SCA is known to converge to a Karush-Kuhn-Tucker solution under mild conditions [13] but the quality of the solution may heavily depend on the initial guess. Here, simple heuristics or solutions to the Travelling Salesman Problem (TSP) or the Pickup-and-Deliver Problem (PDP) can be used for finding an initial trajectory [14]. With direct methods, because of the discretization, the differential equations and the constraints of the systems are satisfied only at discrete points. This can lead to less accurate solutions than indirect methods and the quality of the solution depends on the quantization step [15]. Although every iteration of SCA has a polynomial time complexity, practical resolution time may dramatically increase with the quantization grid and the dimension of the problem. On the other hand, indirect approaches relies on considering the Hamilton-Jacobi Partial Differential Equation associated to the optimal control problem (see e.g., [16], [17][chp. 10]). Several recent methods have been proposed to solve HJ Partial Differential Equations (PDE) in high dimensions. These include max-plus algebra methods [18], [19], dynamic programming and reinforcement learning [20], tensor decomposition techniques [21], sparse grids [22], model order reduction [23], polynomial approximation [24], optimization methods [25]- [28] and neural networks [29]- [32]. In this paper, we consider certain indirect methods that provide analytical solutions for certain classes of optimal control problem as we have shown in a preliminary study [33].

B. Contributions
Our contributions are the following: • Problem Formulation: To the best of our knowledge, this is the first time, after our preliminary study [33], that the UAV BS trajectory problem is formulated using the formalism of Lagrangian mechanics and solved using Hamilton-Jacobi equations. This approach provides closed-form equations when the potential is quadratic and thus very low complexity solutions compared to existing solutions in the literature.
• Closed-form expression of the optimal trajectory with single phase traffic intensity: When the traffic intensity map is made of a single hot spot or traffic hole, has a quadratic form (single phase), and is time-independent, closed form expressions for the optimal trajectory are derived. It follows a hyperbola for a hot spot and corresponds to a repulsor in mechanics. For a traffic hole, the trajectory is on an ellipse and corresponds to the case of an attractor in mechanics.
• Characterization of the optimal solution in multi-phase traffic intensity: When the traffic map has several hot spots or traffic holes (multi-phase) whose regions are separated by interfaces and is timeindependent, we derive necessary conditions to be fulfilled by the position and the instant at which the optimal trajectory crosses an interface (see Theorem 2).
• A gradient algorithm for bi-phase traffic: An in-depth analysis of convexity vs. non-convexity issues allows us to derive a gradient algorithm to solve the bi-phase problem (Algorithm 1). This algorithm finds a stationary point for the cost function. This algorithm has a complexity O(1) at every iteration, whereas iterations of the sequential convex optimization technique have polynomial time complexity.
• A new algorithm for the bi-phase optimization problem: A new algorithm, called the B-algorithm (Algorithm 2), is proposed based on the linear control properties of the quadratic model. This algorithm relies on a bisection scheme the complexity of which is proportional to the logarithm of the desired precision and closed form formulas.
• A data processing procedure: We propose a method to pre-process real measured traffic data in order to derive a bi-phase quadratic model. This procedure is based on smoothing steps followed by a modified K-means algorithm adapted to our quadratic model (Algorithm 3). In our numerical experiments, the optimal trajectory is computed in a region where real traffic data is available [34].
The paper is structured as follows. In Section II, we give the system model, its interpretation in terms of Lagrangian mechanics formulate the problem and give preliminary results. Section III is devoted to the characterization of the optimal trajectories for both the single-and bi-phase cases. Section IV presents our algorithms, Section V our data processing procedure and numerical experiments. Section VI concludes the paper.
Notations: The usual Euclidean scalar product between x ∈ R n and y ∈ R is denoted by x · y. The Euclidean norm x in R n of x ∈ R n is defined by x := √ x · x. The set of matrices with m rows, n columns and real entries is denoted by M m,n (R). The transpose of the A ∈ M m,n (R) is denoted by A † ∈ M n,m (R). We classically identify M m,1 (R) and M 1,n (R) as column vectors of R m and row vectors of R n , respectively. Let f : R n × R m → R defined by f (x, y) where x = (x 1 , . . . , x n ) ∈ R n and y = (y 1 , . . . , y m ) ∈ R m . Let a ∈ R n and b ∈ R m . We denote by ∂f ∂x i (a, b) the partial derivative of f with respect to the variable x i at (a, b) ∈ R n × R m . We also introduce the notations We also consider the following notation for partial Hessian matrices .
and Id n denotes the identity matrix of M n,n (R). We shall see that the value function S : R × R 2 × R × R 2 → R will play a fundamental role in this paper. We use the notations (T 1 , X 1 , T 2 , X 2 ) for S and therefore the partial derivatives of S at (t 1 , x 1 , t 2 , x 2 ) ∈ R × R 2 × R × R 2 are denoted as follows: x 2 ) and ∇ X 2 S(t 1 , x 1 , t 2 , x 2 ).

A. System Model
We consider a network area characterized by a traffic density at position z and time t. We intend to control the trajectory and the velocity of a UAV base station, which is located in z 0 z(t 0 ) at t 0 and shall reach a destination z T z(T ) at T with the aim of minimizing a cost determined by the velocity and the traffic, defined hereafter by (2). At (t, z), we assume that the UAV BS is able to cover an area, from which it can serve users (see Figure 1). The velocity of the UAV BS induces an energy cost. In this model, we control the velocity vector a of the UAV BS. The general form of the cost function is as follows where the first term is a cost related to the velocity of the vehicle (K is a positive constant), and · denotes the usual Euclidean norm. The higher is the speed, the higher is the energy cost. The second term is a user traffic intensity, i.e., the amount of traffic served by the UAV BS at (t, z). Note that a non-zero energy at null speed can be incorporated in the model by adding a constant. Without loss of generality, we assume that this constant is null.
Let S(t 0 , z 0 , T, z T ) be the minimal total cost along any trajectory between z 0 at t 0 and z T at T (also called the action in mechanics or value function in control theory). Let us define Ω(t 0 , T ) as the space of absolutely continuous functions from [t 0 ; T ] to R 2 . Our problem can now be formulated as follows where dz dt (t) = a(t), z(t 0 ) = z 0 , and J is the terminal cost defined by J(z) = 0 if z = z T and J(z) = +∞ otherwise. For simplicity reasons, we assume the existence and uniqueness of the optimal control a * (t) in (3) and denote the associated optimal trajectory z * (t). In a traffic map symmetric with respect to z 0 and z T , the reader can convince himself that the uniqueness is not guaranteed.

B. Preliminary Results From Lagrangian Mechanics
We provide in this section important results from the Lagrangian mechanics for the convenience of the reader.
In the Newtonian classical framework that is used here (see (2)), the impulsion is the product of the particle mass by its velocity (hence the standard term "impulsion").

Definition 2.
The Hamiltonian function is defined as H(t, z, p) := max a∈R 2 p · a − L(t, z, a).
Proof. See Appendix A.
The Euler-Lagrange equation is the first-order necessary condition for optimality and holds for every point on the optimal trajectory. Lemma 2. If the Lagrangian L(t, z, a) is time-independent and α-homogeneous in z and a for α > 0, i.e., L(λz, λa) = |λ| α L(z, a) for all λ ∈ R, S given by (3) reads Proof. See Appendix B.
Proof. See Appendix C.
From now, we assume that the Lagrangian is time-independent, i.e., L(t, z, a) = L(z, a), and is an even function in a, i.e., L(z, −a) = L(z, a). A direct consequence is that H is time-independent and is an even function in p, i.e., we write H(t, z, p) = H(z, p) and H(z, −p) = H(z, p).

III. OPTIMAL TRAJECTORY
In this section, we characterize the optimal trajectory when the traffic intensity is a quadratic form and also when it is made of two regions of quadratic form separated by an interface. We call these two cases single-phase and bi-phase intensities respectively. Both cases satisfy our assumptions on the Lagrangian with α = 2.
A. Single-Phase Optimal Trajectory Assume that the traffic intensity is of the form u(z) = 1 2 u 0 ||z|| 2 . When u 0 > 0, this function models a traffic hole in z = 0. When u 0 < 0, it models a traffic hot spot at z = 0. We disregard the case u 0 = 0 because it corresponds to a constant traffic intensity that is not of interest in this paper. Thus the cost function has the following form Note that

1) Trajectory Equation:
In the single phase case, we have a closed form expression of the trajectory.
Theorem 1. If u 0 < 0, the cost function is given by (14), the optimal trajectory is and the control is given by where ω 2 = − u 0 K . If u 0 > 0, the cost function is given by (17), the optimal trajectory is and the control is given by where ω 2 = u 0 K .
Proof. See Appendix D.
Corollary 1. If the user traffic intensity is of the form u(t, z) = 1 2 u 0 ||z|| 2 + u 0 z · e + u 1 with u 1 ∈ R and e ∈ R 2 , then definez = z + e,z 0 = z 0 + e,z T = z T + e and trajectories given in Theorem 1 are valid by replacing z, z 0 , z T byz,z 0 ,z T , respectively. The cost function becomes: The system is thus equivalent to the one assumed in Theorem 1 by changing the origin of the locations to the barycentre z b of the z i .
2) Traffic Hot Spot, Traffic Hole: We assume that there is a hot spot or a traffic hole located in z h and that the traffic intensity is of the form u(t, z) We can apply Corollary 1 with e = −z h . Figure 2 shows optimal trajectories when z h is a hot spot, i.e., for u 0 < 0, and different values of K. The starting point is z 0 and the destination is z T . When K increases, the velocity cost increases and the trajectories tend to the straight line between z 0 and z T , which minimizes the speed. When K is small, the UAV can go fast to z h , reduces its speed in the vicinity of the hot spot and then goes fast to the destination (in order to decrease the cost function (2) by increasing its traffic contribution). Figure 3 shows optimal trajectories when z h is a traffic hole, i.e., for u 0 > 0. In Figure 3a, T is smaller than the period of the ellipse, i.e., 2π ω > T . When K decreases, the UAV can spend more time in the areas of higher traffic intensity. In Figure 3b, T is larger than the period. In this case, the trajectory follows one period of the ellipse whose equation is given by (18) plus a part of the same ellipse from z 0 to z T .

B. Multi-Phase Trajectory Characterization
We now consider a traffic intensity (or potential) consisting in two quadratic functions separated by an interface I of equal potentials delimiting two regions 1 and 2. The interface is defined by an equation where C is a constant and f is a differentiable function. We assume that the optimal trajectory   crosses the interface only once, at position ξ and time τ . The impulsion p * is defined everywhere on the optimal trajectory between (t 0 , z 0 ) and (T, z T ). The following notations will be used in the sequel: Theorem 2. The location and time (ξ, τ ) of interface crossing are characterized by the following equations for some Lagrange multiplier µ ∈ R Proof. See Appendix E.
Equation (22) expresses the fact the energy is conserved when crossing the interface. One can show that actually the energy is conserved along the whole trajectory. Equation (21) is related to the conservation of the tangential component of the impulsion at the interface. Equation (23) is the interface equation at ξ. One can show that under the assumption of equal potential on the interface, the kinetic energy, the impulsion, and the velocity vector are conserved across the interface.

C. Multi-phase optimal trajectory uniqueness and value function convexity issues
Looking for uniqueness/non-uniqueness of optimal multi-phase trajectories leads naturally to study the convexity of the multi-phase total cost. As stated in Appendix E this cost is additive and satisfies the dynamic programming principle which reads where S 1 and S 2 are themselves minimal. Hamilton-Jacobi equations applied to each cost component allows to study their first-an second-order differential properties such as convexity. For instance in the chosen quadratic model, each single-phase (minimal) total cost S 1 (t 0 , z 0 , τ, ξ), S 2 (τ, ξ, T, z T ) given by Corollary 1 and (14) is convex (since it is a positive quadratic form with respect to the spatial coordinates) with respect to the spatial position ξ, but not necessarily so with respect to the interface crossing time τ . We consider in what follows a general form of the single-phase value function between (t 1 , x 1 ) and (t 2 , x 2 ) that is denoted by We also note by ω the pulsation, φ = ω (t 2 − t 1 ) is the temporal phase and p 1 , p 2 are the the initial and final impulsions as derived from formula (16).
where p i is the impulsion at time T i , i.e., at extremity X i (i = 1, 2 and j = 3 − i).
ii) At most one eigenvalue of H(ψ i ) can be negative and α < 0 is a sufficient condition for this to hold, namely implying local non-convexity of the (single-phase) value function.
iii) The double-phase Hessian enjoys a similar structure as (25) and writes with respect to the variable and may thus be non-convex as well.
Proof. See Appendix F subsections F1 and F2.
In the single-phase case, the sufficient condition ii) i.e., α < 0, supports first a physical interpretation (see Fig. 4) and also a geometrical one (see next Theorem 4 and Fig. 5).  and v = x 2 −x 1 2 . (Note that u is the distance from the hotspot to the average of x 1 and x 2 , and the v is the half-distance between x 1 and x 2 .) Then the following sufficient condition Proof: We combine the formulas (53) and (54) in (58) to obtain which can be rewritten as This non-convexity condition for the single-phase value function interprets as: long phase (tanh φ 2 > 1 2 ) and long distances between the hotspot z h to both initial and final positions x 1 , x 2 , relatively to their mutual distance ( u v ≫ 1). Several preliminary simulations show that when this case happens for two phases, the total value function is indeed non-convex and that several (local) optimal solutions may indeed exist.

IV. ALGORITHMS
Previous propositions allow us to propose several numerical algorithms seeking optimal trajectories. In this section, we present two algorithms aiming this goal: a gradient descent algorithm and a bisection search method based on the linear control property (16).

A. A Gradient Descent Algorithm GRAD-ALGO
In this section, we propose GRAD-ALGO (Algorithm 1) which is an alternated optimization-based algorithm relying on the following procedures. a) Procedure for seeking an optimal ξ given a fixed τ : from Hamilton-Jacobi equations (Lemma 3 and Appendix C) the gradient of the total cost function with respect to ξ is p * (τ − )−p * (τ +). Equation (13) says that in the Newtonian framework the impulsion is proportional to the control variable a. Equation (16) says that in the quadratic model the velocity vector is, at any time a linear combination of centered initial and final positions. Then the searched gradient appears to be an affine function of ξ which reads Scalar h and point B (where the spatial gradient cancels at fixed τ i.e., p * (τ − ) = p * (τ +)) verify .
The equation involving the Lagrange multiplier (21) now reads and shows that the optimal location ξ * is the orthogonal projection of B on the interface I. b) Procedure for seeking an optimal τ given a fixed ξ: we use the result of Theorem 2. As also shown by Hamilton-Jacobi equations (Lemma 3 and Appendix C), the gradient of S with respect to τ is given by H 2 (ξ, p * (τ + )) − H 1 (ξ, p * (τ − )). The Hamiltonians are easily computed in each phase by applying the classical Newton formula at given location (z, t) We then update τ by using a simple gradient descent scheme (see Algorithm 1).
c) Stop criteria: the search procedure is performed until the impulsion and the Hamiltonian have converged with a given accuracy. For this purpose we consider the function g : which measures a relative "discrepancy" between x and y.

B. The B-curve algorithm B-ALGO
The B-curve algorithm aims to overcome the non-convexity issues developed previously. It proceeds as follows. Consider the B point defined in formula (30). We notice that for each τ ∈ [t 0 , T ] this point is defined univocally knowing all parameters (t 0 , T, z 0 , z T , z h ) and all (quadratic) traffic profiles so that we can see B in (30)  u i (z) = (time-stationary) traffic u(z) from hotspot z hi at point z (i = 1, 2) return (u 2 (z) > u 1 (z)) end procedure repeat Then by formula (28) the spatial gradient of S at point ξ is This implies that the kinetic components of both Hamiltonians are equal at the interface. Now since ξ belongs to the interface, both phase potentials (traffics) are equal by definition, so that the total Hamiltonian is also conserved at the interface: this is the optimality condition required with respect to time τ . To summarize, in these conditions, both the Hamiltonian and the total impulsion are conserved at interface I, i.e., local optimality conditions hold for the total value function. It is also worth noticing that the related  In this section, we explain how from measured or estimated traffic load, we can derive a quadratic model that will allow us to apply our framework. To illustrate the procedure, we extract data from the open data set presented in [34]. Traffic data (in number of bytes) has been collected from an operational cellular network in a medium-size city in China and is recorded for every base station and every hour.  Figure 6a. We assume that this traffic is representative of the traffic intensity when the drone is launched. In order to fit this raw data to our model, we follow the pre-processing steps shown in Algorithm 3. The first smoothing consists in aggregating the traffic data on a grid of 50 steps in both longitude and latitude directions. The resulting elementary regions should correspond approximately to the drone coverage. The result is shown in Figure 6b. The data exhibits a very high variability with very high peaks around few locations. The second smoothing is a Locally Weighted Scatterplot Smoothing (LOWESS) [36]. We use here the Matlab function fit with the option "Lowess". The choice of the smoothing parameter α, i.e., the proportion of data points used for every local regression, has a decisive impact on the result. Increasing α has the effect of averaging out the different peaks. In our specific scenario, α = 0.25 yields Figure 7a with two local maxima. With α = 0.5, we obtain a single maximum.
In step 3 of the pre-processing, the traffic is normalized between 0 and 1 with no influence on the optimal trajectory.
The final pre-processing step is an adaptation of the classical K-means algorithm (see Algorithm 4) to fit to quadratic models. Inputs are the data points obtained after the normalization, K c , the number of clusters (or hot spots), K n the number of nearest neighbors and M the number of iterations. Every cluster is associated to a quadratic function (in our case, we have K c = 2). Every data point j is associated to a cluster and has a related label L j in {1, ..., K}. For every data point, a list of nearest neighbors is built (step 4). An arbitrary initial labelization is chosen (step 6). The algorithm then proceeds by iterations (steps [7][8][9][10][11][12][13][14][15][16]. At every iteration, if a point j has some neighbor with a different label (step 9), a best new label is found for j (step [10][11][12][13] in terms of quadratic error e k . The error e k measures the difference between the data points and the K c -quadratic model, which fits a quadratic function to every cluster assuming that j has label k (steps [18][19][20][21][22][23][24][25][26][27][28][29][30][31]. In Algorithm 4, KNN(K n , X, Y ) is a procedure that finds the K n nearest neighbors of (X, Y ) with respect to the Euclidian distance. We use the Matlab implementation knnsearch. NLLS(L) is a non-linear least square method that fits data points in L to a quadratic function of the form 1 2 u 0 |z − z h | 2 + u 1 . The Matlab implementation is based on a trust-region approach of the Levenberg-Marquardt Algorithm.  Figure 8 shows the quadratic error as a function of the number of iterations of the K-means algorithm for K c = 1 and K c = 2. The error is constant for K c = 1 as there is only one iteration, which performs the non-linear least square fitting for the single cluster. For K c = 2, we distinguish two cases: K n = 5 and K n = ∞. In the former case, only the 5 nearest neighbors of a data point are inspected to decide if a relabelization should be performed. In the later case, relabelization is systematically considered. Increasing K n increases the complexity of every iteration but provides a faster convergence.

B. Trajectory optimization results
1) Estimation procedure of parameters K and T: as in most similar algorithms, a 'prior' estimation of parameters is necessary since the optimal results strongly depend on them. Traffic parameters have been estimated just above, so the required parameters are first the mass K representing trade-off between the kinetic term and the traffic term in Lagrangian (2) and then the total duration time T which has also Algorithm 4 K-means with quadratic models 1: Input: K c (number of clusters), J (number of measurement points) (X j , Y j , Z j ), j = 1, . . . , J (coordinates and estimated traffic load Z j in (X j , Y j )), K n (number of nearest neighbors), M (number of iterations) 2: Output: z hk , u 0k , u 1k k = 1, . . . , K c (hot spot characteristics), e (quadratic error) 3: for j = 1, . . . , J do 4:  return e k , z hl , u 0l , u 1l , l = 1, . . . , K c 31: end procedure a significant influence on the optimal trajectory (by convention t 0 = 0). This procedure is developed in Appendix G and yields in our case: T ≈ 850, K = 50 × r 2 with r = 40 × 10 6 360 × 100 and we select: T = 1200, 1800s and K = 30, 60 × r 2 .
2) Results: optimal trajectories are shown in Figs. 9 and 10 and related numerical results in Table I.
Length of various trajectories and related average velocities are given in Table II.
The following main comments can be made on these results:     Iso-traffic levels Hot spot Single-phase Algorithm T=1200 Single-phase Algorithm T=1800 (b) K = 60, T = 1200 s and 1800 s, 1 phase. Fig. 9: Influence of T on the optimal trajectory.
• As already noted for the single-phase case (paragraph III-A2), decreasing mass K enables the drone to collect more traffic. In fact, during the allowed time interval, the drone will get closer to the hotspot with maximal traffic in order to decrease the total value function, i.e., zone 1 here.
• Increasing time interval (T ) will produce the same tendancy, allowing the drone to spend more time near the hotspot with higher traffic z h1 , see Figs. 9a and 10a.
• Since the drone has enough time to pick up traffic located close to hotspots, this is reflected in an average drone velocity about half its nominal value (see Table II).   and r = 40×10 6 360×100 ≈ 1111.1111 (corresponding to scaled spatial unity of 100m). Also, S 1ph denotes the total single-phase cost (i.e., assuming only one, effective hotspot).

VI. CONCLUSION
In this paper, we propose a Lagrangian approach to solve the UAV base station optimal trajectory problem. When the traffic intensity exhibits a single phase, closed-form expressions for the trajectory and speed are derived from Hamilton-Jacobi equations. When the traffic intensity exhibits multiple phases, we characterize the crossing time and location at the interface. We propose two low-complexity algorithms for  the bi-phase time-stationary traffic case that provide optimal crossing time and location on the interface and fulfill the necessary conditions of optimality. At last, we present a data processing procedure based on a modified K-means algorithm that derives a single-phase or bi-phase quadratic model from real traffic data. Further extensions of this work are envisioned to generalize the approach to three or more hot-spots and to consider multi-drone coordinated trajectories.

A. Proof of Lemma 1
In a neighborhood of the optimal trajectory, the first order variation of S is null δS = [∇ z L(t, z, a) · δz(t) + ∇ a L(t, z, a) · δa(t)] dt.
We now note that δa = δ dz dt = d(δz) dt . Integrating by part the second term in the integral of δS, we obtain Note that [δz ∂L ∂a ] T t 0 = 0 because z 0 and z T are fixed. Equating δS to zero gives As this should be true for every δz, L, z 0 and z T , we obtain the first result.
Note that z 0 is fixed and so δz in z 0 is null. We thus obtain the second result of the lemma.
Using (7) and (39), we have: z · dp dt +ż · p = αL or equivalently d(p·z) dt = αL. We can now integrate the cost function (3) along the optimal trajectory as follows C. Proof of Lemma 3 We assume that an optimal trajectory exists and we apply the principle of optimality on the optimal trajectory between (t, z * (t)) and (t + h, z * (t) + ah), where h > 0, to opbtain This implies that By using the same approach between t − h and t, we deduce in the same way equation (10)  From (6) and (12), we obtain the following ordinary differential equation of second degree:z = − u 0 K z. If u 0 K > 0, we define ω 2 = u 0 K and we look for an optimal trajectory of the form: z(t) = A cos(ωt)+B sin(ωt). If u 0 K < 0, we look for an optimal trajectory of the form: z(t) = A cosh(ωt) + B sinh(ωt) with ω 2 = − u 0 K . Let us denote z 0 = z(t 0 ) and a 0 = a(t 0 ) the initial conditions for z andż.
Take the case u 0 K < 0. Using the derivative of z(t) and identifying terms, we obtain: when u 0 < 0. In a similar way, we have when u 0 > 0. Injecting a(t 0 ) = a 0 in the equation of the trajectory provides the result.
For the computation of S, we now use the result of Lemma 2 as our cost function is 2-homogeneous.
From equation (8), we see that only initial and final conditions are required to compute the cost function.
Recall now that p = Ka. Injecting the equations of a(t 0 ) and a(T ) in (8), we obtain the result for the cost function.

E. Proof of Theorem 2
We assume that the location and time (ξ, τ ) of interface crossing is known and unique. The optimal trajectory between (z 0 , t 0 ) and (z T , T ) can be decomposed in two sub-trajectories that are themselves optimal between (z 0 , t 0 ) and (ξ, τ ) on the one hand and between (ξ, τ ) and (z T , T ) on the other hand, by the principle of optimality. In region 1, the optimal cost up to τ is Using Hamilton-Jacobi, we obtain In the same way, the optimal cost in region 2 is Using again Hamilton-Jacobi, we obtain A necessary condition for the optimality of τ is thus that is A necessary condition for the optimality of ξ in the total cost under the constraint f (ξ) = C is also where µ is a Lagrange multiplier associated to the constraint and where the second line comes from equation (11) of Hamilton-Jacobi. Thus we obtain precisely equation (21).
F. The structure and the diagonalization of single-and two-phase Hessian of the value fonction 1) Proof of Theorem 3 : Recall that we study the second-order differentiability properties of singlephase value function S(t 1 , x 1 , t 2 , x 2 ). Recall that from Eq. (1) we study Hessians of the form  where σ stands for S(t 1 , x 1 , t 2 , x 2 ) and i = 1, 2 It is clear by inspecting the symmetries of (14) with respect to the spatial coordinates that where g = ω coth φ and φ = ω (t 2 − t 1 ) is the temporal phase.
Concerning symmetry with respect to time variables, one also finds easily that The Hessian with respect to the variable χ i = (T i , X i ) i=1,2 ∈ R × R 2 has thus the following structure: where the vector Π i and real scalar α remain to be determined. Now recall that initial and final impulsions p 1 and p 2 as given by (41) write: In fact the expression of Π i and α appearing in Hessian (52) result both from the computation of ∂p i ∂T i for i = 1, 2. Using (51), the differentiation of (53) and (54) wrt. T 1 resp. T 2 leads to the following simple result Now, as far as the second partial derivative of value function S with respect to time (α) is concerned Recalling that the Hamiltonian for any Newtonian model has the form H(z, p, t) = and differentiating it with respect to the temporal phase φ at fixed extremities x 1 and x 2 leads to Now from (24) the two-phase Hessian is simply the sum of the two single-phase Hessians.
2) Diagonalizing the Hessian of single-and two-phase value function : since the Hessians of singleand two-phase value function do have the same structure, we address a general Hessian of the form (52).
However the condition α Kg < Π 2 (e.g. induced by the sufficient condition α < 0) implies that one of eigenvalues ν 1 , ν 2 < 0 namely the local non-convexity of the Hessian matrix H.
This property holds thus for each single-phase value function as well as for the total two-phase case.

G. Estimation of parameters
It is based on the following remarks. First, the interface location I is independent of K and T (since only depends on traffic variables). Then, final time T remains always estimated in seconds (t 0 = 0). Also, according to [38], a typical drone cell flights at speedV = 20 m/s and has an autonomy of about 28 min.
Last, we test our procedure for spatially-scaled data, where the convenient spatial unity is 100 m, so that the scaling ratio is Passing back into the original frame therefore implies that K ← K r 2 .