Optimal Trajectories of a UAV Base Station Using Lagrangian Mechanics

TL;DR

This paper formulates the optimal UAV base station trajectory problem using Lagrangian mechanics, deriving closed-form solutions for quadratic traffic and proposing an iterative algorithm for more complex traffic patterns, with numerical validation.

Contribution

It introduces a novel application of Lagrangian mechanics to optimize UAV trajectories and develops an algorithm for bi-phase traffic scenarios, advancing UAV mobility optimization methods.

Findings

Closed-form solutions for quadratic traffic intensity.

An iterative algorithm for bi-phase traffic optimization.

Numerical results demonstrating improved trajectories over MPC.

Abstract

In this paper, we consider the problem of optimizing the trajectory of an Unmanned Aerial Vehicle (UAV) Base Station (BS). We consider a map characterized by a traffic intensity of users to be served. The UAV BS must travel from a given initial location at an initial time 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. When the traffic intensity is quadratic (single-phase), we derive closed-form formulas for the optimal trajectory. When the traffic intensity is bi-phase, we provide necessary conditions of optimality and propose an Alternating Optimization Algorithm that returns a trajectory satisfying these conditions. The Algorithm is initialized with a Model Predictive Control (MPC) online algorithm. Numerical results show how we improve the trajectory with respect to the MPC solution.

Equations

$$\mathcal{L}(t,z,a)=\frac{K}{2}||a|{|}^2-u(t,z)$$

$$S(t_0,z_0,T,z_T)$$

$$p(t,z,a):={\nabla }_a\mathcal{L}(t,z,a)$$

$$H(t,z,p):=\underset{a\in {\mathbb{R}}^2}{\max }p\cdot a-\mathcal{L}(t,z,a).$$

$$\frac{d}{dt}{\nabla }_a\mathcal{L}(t,z^{\ast }(t),a^{\ast }(t))={\nabla }_z\mathcal{L}(t,z^{\ast }(t),a^{\ast }(t))$$

$$\frac{dp}{dt}(t,z^{\ast }(t),a^{\ast }(t))={\nabla }_z\mathcal{L}(t,z^{\ast }(t),a^{\ast }(t))$$

$$S(t_0,z_0,T,z_T)=\frac{1}{\alpha }[z\cdot p{]}_{t_0}^T+J(z_T).$$

$$\frac{\partial S}{\partial t_0}(t,z^{\ast }(t),T,z_T)=H(t,z^{\ast }(t),-p^{\ast }(t))$$

$$\frac{\partial S}{\partial T}(t_0,z_0,t,z^{\ast }(t))=-H(t,z^{\ast }(t),p^{\ast }(t))$$

$$p^{\ast }(t)={\nabla }_a\mathcal{L}(t,z^{\ast }(t),a^{\ast }(t))={\nabla }_{z}S(t_{,}{z}^{\ast }(t),T,z_T)$$

$$\mathcal{L}(z,a)=\frac{1}{2}K||a|{|}^2-\frac{1}{2}{u}_0||z|{|}^2$$

$$p(z,a)={\nabla }_a\mathcal{L}(z,a)=Ka$$

$$S(t_0,z_0,T,z_T)=\frac{K\omega }{2\sinh \omega (T-t_0)}\left((|z_0{|}^2+|z_T{|}^2)\cosh \omega (T-t_0)-2z_0\cdot z_T\right)+J(z_T)$$

$$z^{\ast }(t)=\frac{z_T\sinh (\omega (t-t_0))+z_0\sinh (\omega (T-t))}{\sinh (\omega (T-t_0))}$$

$$a^{\ast }(t)=\omega \frac{z_T\cosh (\omega (t-T))-z_0\cosh (\omega (T-t))}{\sinh (\omega (T-t_0))}$$

$$S(t_0,z_0,T,z_T)=\frac{K\omega }{2\sin \omega (T-t_0)}\left((|z_0{|}^2+|z_T{|}^2)\cos \omega (T-t_0)-2z_0\cdot z_T\right)+J(z_T)$$

$$z^{\ast }(t)=\frac{z_T\sin (\omega (t-t_0))+z_0\sin (\omega (T-t))}{\sin (\omega (T-t_0))}$$

$$a^{\ast }(t)=\omega \frac{z_T\cos (\omega (t-t_0))-z_0\cos (\omega (T-t))}{\sin (\omega (T-t_0))}$$

$$H_1(\xi (\tau ),p^{\ast }({\tau }^{-}))-H_2(\xi (\tau ),p^{\ast }({\tau }^{+}))$$

$$p^{\ast }({\tau }^{-})-p^{\ast }({\tau }^{+})-\mu {\nabla }_{z}f(\xi )$$

$$f(\xi )$$

$${\nabla }_{z_T}{S}_1(t_0,z_0,\tau ,\xi )+{\nabla }_{z_0}{S}_2(\tau ,\xi ,T,z_T)=Kh(\xi -B)$$

$$h={\omega }_1\coth ({\omega }_1(\tau -t_0))+{\omega }_2\coth ({\omega }_2(T-\tau ))$$

$$B$$

$$Kh(\xi -B)-\mu {\nabla }_{\xi }f(\xi )=0$$

$$\delta S$$

$${\int }_{t_0}^T{\nabla }_a\mathcal{L}(t,z,a)\cdot \frac{d(\delta z)}{dt}dt$$

$$0={\int }_{t_0}^T\left[{\nabla }_z\mathcal{L}(t,z,a)-\frac{d}{dt}{\nabla }_a\mathcal{L}(t,z,a)\right]\cdot \delta z(t)dt.$$

$$\delta S$$

$$z\cdot \frac{\partial \mathcal{L}(z,\dot{z})}{\partial z}+\dot{z}\cdot \frac{\partial \mathcal{L}(z,\dot{z})}{\partial z}$$

$$S(t_0,z_0,T,z_T)$$

Full text

I Introduction

Unmanned Aerial Vehicles (UAV) are expected to play an increasing role in future wireless networks111J. Darbon is supported by NSF DMS-1820821. [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 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.

I-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 problem consists in optimizing the trajectory length of the UAV under the constraint that it collects the required amount of data from every sensor. Authors use a reinforcement learning approach where improved trajectories are sequentially learned over several tour iterations. This model is different from ours as it allows the UAV to learn the optimal trajectory from previous experience. 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 a very recent work [5], a data rate-energy trade-off is studied. 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 [6]. In this aim, there are classical methods like Cell Decomposition, Potential Field Method or Probabilistic Road Map and there are heuristic approaches, e.g. bio-inspired algorithms. Authors of [7] 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, and in [7] the speed cannot be controlled. 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 [8] and is applied for example in aircraft trajectory planning [9]. Most numerical methods in control theory can be classified in direct and indirect methods. Indirect methods provide analytical solutions from the calculus of variations and use first order necessary conditions for a trajectory to be optimal. 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 very recent publications in the field of UAV-aided communications. In [10] for example, a UAV relay assists the communication between a source and a destination. As the resulting problem is non-convex, it is first approximated and then solved by successive convex optimization. In [11], the objective is to maximize the energy efficiency of a UAV-to-ground station communication by taking into account the propulsion energy consumption and by optimizing the trajectory. Again, sequential convex optimization is applied to an approximated problem. In the same vein, [12] considers multiple-UAV BSs used to serve fixed users. The quality of the solution to the nonlinear program may heavily depend on the initial guess. Authors thus propose an heuristic based on circular trajectories to initialize their algorithm. 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 [13]. Although every iteration of the sequential convex optimization technique has a polynomial time complexity, practical resolution time may dramatically increase with the quantization grid and the dimension of the problem. We thus propose in this paper an indirect approach based on Lagrangian mechanics that has the advantage to provide closed-form expressions for the optimal trajectories when the potential is quadratic (we say single-phase). When the potential is quadratic by region (or multi-phase) the optimization process consists in finding the right crossing time and location on the interface of the regions. This question is an active field of research in control theory, see e.g. [14]. As explained in [15, 16], a trajectory optimization problem can be decomposed in different phases or arcs. Phases are sequential in time, i.e., they partition the time domain. Differential equations describing the system dynamics cannot change during a phase. This point of view allows us to consider the multi-phase problem.

I-B Contributions

Our contributions are the following:

• •

Problem Formulation: To the best of our knowledge, this is the first time that the UAV BS trajectory problem is formulated using the formalism of Lagrangian mechanics. 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 consists in a part of hyperbole 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 time-independent, we derive necessary conditions to be fulfilled by the position and the instant at which the optimal trajectory crosses an interface (see Theorem 2).

• •

An online algorithm for multi-phase time-varying traffic intensity: When the traffic map is multi-phase and is time-varying, we propose an online algorithm based on MPC.

• •

An Alternating Optimization Algorithm for bi-phase time-independent traffic intensity: When the traffic intensity is made of two hot spots separated by an interface (bi-phase) and is time-independent, we propose an Alternating Optimization Algorithm that 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 (see Algorithm 1).

The paper is structured as follows. In Section II we give the system model and its interpretation in terms of Lagrangian mechanics. In Section III, we formulate the problem and give preliminary results. Section IV is devoted to the characterization of the optimal trajectories. Section V presents our algorithms and Section VI concludes the paper.

Notations: Let $f:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ defined by $f(x,y)$ where $x=(x_{1},\dots,x_{n})\in\mathbb{R}^{n}$ and $y=(y_{1},\dots,y_{m})\in\mathbb{R}^{m}$. Let $a\in\mathbb{R}^{n}$ and $b\in\mathbb{R}^{m}$. We denote by $\frac{\partial f}{\partial x_{i}}(a,b)$ the partial derivative of $f$ with respect to the variable $x_{i}$ at $(a,b)\in\mathbb{R}^{n}\times\mathbb{R}^{m}$. We also introduce the notations $\nabla_{x}f(a,b)=(\frac{\partial f}{\partial x_{1}}(a,b),\dots,\frac{\partial f}{\partial x_{n}}(a,b))\in\mathbb{R}^{n}$ and $\nabla_{y}f(a,b)=(\frac{\partial f}{\partial y_{1}}(a,b),\dots,\frac{\partial f}{\partial y_{m}}(a,b))\in\mathbb{R}^{m}$.

II System Model and Interpretation

II-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}\triangleq z(t_{0})$ at $t_{0}$ and shall reach a destination $z_{T}\triangleq z(T)$ at $T$ with the aim of minimizing a cost determined by the velocity and the traffic, defined hereafter by (1). 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

[TABLE]

where the first term is a cost related to the velocity of the vehicle ($K$ is a positive constant), and $\|\cdot\|$ 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.

III Lagrangian Mechanics Formulation

III-A Problem Formulation

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 $\Omega(t_{0},T)$ as the space of absolutely continuous functions from $[t_{0};T]$ to $\mathbb{R}^{2}$. Our problem can now be formulated as follows

[TABLE]

where $\frac{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)=+\infty$ otherwise. For simplicity reasons, we assume the existence and uniqueness of the optimal control $a^{*}(t)$ in (2) 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.

III-B Preliminary Results From Lagrangian Mechanics

We provide in this section important results from the Lagrangian mechanics for the convenience of the reader.

Definition 1** (Impulsion).**

The impulsion function is defined as

[TABLE]

In the Newtonian classical framework that is used here (see (1)), 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

[TABLE]

Lemma 1** (Euler-Lagrange Equations).**

Along the optimal trajectory $z^{*}(t)$ that starts from $z_{0}$ at $t_{0}$ and ends at $z_{T}$ at $T$, we have

[TABLE]

or equivalently

[TABLE]

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 $\mathcal{L}(t,z,a)$ is time-independent and $\alpha$-homogeneous in $z$ and $a$ for $\alpha>0$, i.e., $\mathcal{L}(\lambda z,\lambda a)=|\lambda|^{\alpha}\mathcal{L}(z,a)$ for all $\lambda\in\mathbb{R}$, $S$ given by (2) reads

[TABLE]

Proof.

See Appendix -B. ∎

Lemma 3** (Hamilton-Jacobi).**

Along the optimal trajectory, we have for $t\in(t_{0};T)$

[TABLE]

[TABLE]

where

[TABLE]

Proof.

See Appendix -C. ∎

From now, we assume that the Lagrangian is time-independent, i.e., $\mathcal{L}(t,z,a)=\mathcal{L}(z,a)$, and is an even function in $a$, i.e., $\mathcal{L}(z,-a)=\mathcal{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)$.

IV 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 interface222We leave for further work the way to approximate a realistic traffic intensity map by a set of regions with intensities of quadratic form.. We call these two cases single-phase and multiple-phase intensities respectively. Both cases satisfy our assumptions on the Lagrangian with $\alpha=2$.

IV-A Single-Phase Optimal Trajectory

Assume that the traffic intensity is of the form $u(z)=\frac{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

[TABLE]

Note that

[TABLE]

IV-A1 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 (13),

*the optimal trajectory is *

[TABLE]

and the control is given by

[TABLE]

*where $\omega^{2}=-\frac{u_{0}}{K}$.

If $u_{0}>0$, the cost function is given by (16),*

*the optimal trajectory is *

[TABLE]

and the control is given by

[TABLE]

where $\omega^{2}=\frac{u_{0}}{K}$.

Proof.

See Appendix -D. ∎

Corollary 1**.**

If the user traffic intensity is of the form $u(t,z)=\frac{1}{2}u_{0}||z||^{2}+u_{0}z\cdot e+u_{1}$ with $u_{1}\in\mathbb{R}$ and $e\in\mathbb{R}^{2}$, then define $\tilde{z}=z+e$, $\tilde{z}_{0}=z_{0}+e$, $\tilde{z}_{T}=z_{T}+e$ and trajectories given in Theorem 1 are valid by replacing $z$, $z_{0}$, $z_{T}$ by $\tilde{z}$, $\tilde{z}_{0}$, $\tilde{z}_{T}$, respectively. The cost function becomes: $S(t_{0},z_{0},T,z_{T})=\frac{1}{\alpha}[z\cdot p]_{t_{0}}^{T}+J(z_{T})-u_{1}(T-t_{0})$.

Corollary 2**.**

If the user traffic intensity is of the form $u(t,z)=\sum_{i}u_{i}||z-z_{i}||^{2}$ with $\sum_{i}u_{i}\neq 0$, then $u(t,z)=(\sum_{i}u_{i})||z-z_{b}||^{2}+\sum_{i}u_{i}||z_{i}-z_{b}||^{2}$ with $z_{b}=\frac{\sum_{i}u_{i}z_{i}}{\sum_{i}u_{i}}$. Define $\tilde{z}=z+z_{b}$, $\tilde{z}_{0}=z_{0}+z_{b}$, $\tilde{z}_{T}=z_{T}+z_{b}$, $\tilde{u}_{0}=\sum_{i}u_{i}$ and trajectories given in Theorem 1 are valid by replacing $z$, $z_{0}$, $z_{T}$, $u_{0}$ by $\tilde{z}$, $\tilde{z}_{0}$, $\tilde{z}_{T}$, $\tilde{u}_{0}$ respectively.

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}$.

IV-A2 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)=\frac{1}{2}u_{0}||z-z_{h}||^{2}+u_{1}=\frac{1}{2}u_{0}||z||^{2}-u_{0}z\cdot z_{h}+\frac{1}{2}u_{0}||z_{h}||^{2}+u_{1}$. 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.

Figure 3 shows optimal trajectories when $z_{h}$ is a traffic hole, i.e., for $u_{0}>0$. In Figure 3(a), $T$ is smaller than the period of the ellipse, i.e., $\frac{2\pi}{\omega}>T$. When $K$ decreases, the UAV can spend more time in the areas of higher traffic intensity. In Figure 3(b), $T$ is larger than the period. In this case, the trajectory follows one period of the ellipse whose equation is given by (17) plus a part of the same ellipse from $z_{0}$ to $z_{T}$.

IV-B Multi-Phase Trajectory Characterization

We now consider a traffic intensity (or potential) consisting in two quadratic functions separated by an interface $\mathcal{I}$ of equal potentials delimiting two regions $1$ and $2$. The interface is defined by an equation $f(z)=C$, where $C$ is a constant and $f$ is a differentiable function. We assume that the optimal trajectory crosses only once the interface at position $\xi$ at $\tau$.

Theorem 2**.**

The location and time $(\xi,\tau)$ of interface crossing are characterized by the following equations

[TABLE]

for some Lagrange multiplier $\mu\in\mathbb{R}$, where we recall that $p^{*}$ is defined with respect to the optimal trajectory between $(t_{0},z_{0})$ and $(T,z_{T})$, and where $p^{*}(\tau^{-})=\lim_{s\to\tau,s<\tau}p^{*}(s)$ and $p^{*}(\tau^{+})=\lim_{s\to\tau,s>\tau}p^{*}(s)$.

Proof.

See Appendix -E. ∎

Equation (19) 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 (20) is related to the conservation of the tangential component of the impulsion at the interface. Equation (21) is the interface equation at $\xi$. 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.

V Algorithms

V-A An Online Algorithm: MPC

In this section, we first present an online algorithm based on MPC [17] (we omit the pseudo-code for space reasons). In a traffic intensity landscape made of multiple phases, the idea is to assume at every $t$ that the current phase won’t change from $t$ to $T$. Using this assumption, we compute the optimal trajectory as in the single phase case and take the next decision based on this trajectory. This algorithm has the advantage of being online, of low complexity and can be used in multiphase time-dependent traffic maps. We have however no guarantee of optimality.

V-B An Alternating Optimization Algorithm

We now study a time-independent bi-phase scenario, in which a trajectory from $z_{0}$ to $z_{T}$ crosses the interface at time $\tau$ and location $\xi$. We present an Alternating Optimization Algorithm (Algorithm 1) that provides a stationary trajectory in the sense of Theorem 2. The algorithm consists in alternatively optimizing $\tau$ (steps 9-17) and $\xi$ (steps 18-26) by using the results of Theorem 2. For every fixed $\tau$ and $\xi$, the current trajectory is the concatenation of the optimal trajectory between $(t_{0},z_{0})$ and $(\tau,\xi)$ and the optimal trajectory between $(\tau,\xi)$ and $(T,z_{T})$ (step 27). Every iteration of the algorithm only requires the evaluation of two Hamiltonians or the computation of a point $B$, see (23), and its projection on the interface. Therefore the complexity of an iteration is $O(1)$. In simulations, MPC is used to produce an initial trajectory.

V-B1 Procedure for seeking an optimal $\tau$ given a fixed $\xi$

We use the result of Theorem 2. As shown in its proof [18], the gradient of $S$ with respect to $\tau$ is given by $H_{2}(\xi,p^{*}(\tau^{+}))-H_{1}(\xi,p^{*}(\tau^{-}))$. We can thus compute the Hamiltonians in every region by differentiating the cost function (13) with respect to the final time in region 1 (see (9)) and with respect to the initial time in region 2 (see (8)). We then update $\tau$ by using a simple gradient descent scheme in step 11.

V-B2 Procedure for seeking an optimal $\xi$ given a fixed $\tau$

From Hamilton-Jacobi, the gradient of the total cost function with respect to $\xi$ is $p^{*}(\tau^{-})-p^{*}(\tau+)$ (see proof of Theorem 2 in [18]). Since in the Newtonian framework the impulsion is proportional to the control variable $a$ (see (12)) and since in a quadratic model the velocity vector is, at any time a linear combination of centered initial and final positions (15), this gradient appears to be an affine function of $\xi$ which reads

[TABLE]

Scalar Hessian $h$ and position $B$, where the spatial gradient cancels i.e., $p^{*}(\tau^{-})=p^{*}(\tau+)$ at fixed $\tau$ are given by:

[TABLE]

[TABLE]

The equation involving the Lagrange multiplier (20) new reads

[TABLE]

and shows that the optimal location $\xi^{*}$ is the orthogonal projection of $B$ on the interface ${\cal I}$. This projection is performed in steps 19-20 of the algorithm.

Figure 4 shows the MPC trajectory and the trajectory obtained from Algorithm 1 after 60 iterations in a bi-phase landscape. The traffic intensity is shown in three dimensions in Figure 1: It is a bi-phase landscape made of two hot-spots, where the peak of traffic in $z_{h2}$ is higher than in $z_{h1}$. The Alternating Optimization Algorithm has gradually moved the interface crossing time and location in order to spend more time in the second hot-spot and to go closer to $z_{h2}$. Figure 5 shows how the cost function has decreased along the iterations and thus how our algorithm has improved over the MPC solution. From iterations 1 to 45, $\tau$ has been gradually updated; at iteration 46, $\xi$ is updated once; $\xi$ is again updated once at iteration 59.

VI Conclusion

In this paper, we have proposed 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 given. When the traffic intensity exhibits multiple phases, we characterize the crossing time and location at the interface. In a first approach, we propose an online algorithm based on MPC for multi-phase and time-dependent traffic intensity, which allows to take into account the impact of each phase. We then propose an offline Alternating Optimization Algorithm for bi-phase time-independent traffic intensities that provides a stationary trajectory with respect to the crossing time and location on the interface and fulfills the necessary conditions of optimality. Numerical results show that we improve the trajectory obtained with MPC.

-A Proof of Lemma 1

Around the optimal trajectory, the first order variation of $S$ is null.

[TABLE]

We now note that $\delta a=\delta\frac{dz}{dt}=\frac{d(\delta z)}{dt}$. Integrating by part the second term in the integral of $\delta S$, we have

[TABLE]

Note that $[\delta z\frac{\partial\mathcal{L}}{\partial a}]_{t_{0}}^{T}=0$ because $z_{0}$ and $z_{T}$ are fixed. Equating $\delta S$ to zero gives now

[TABLE]

As this should be true for every $\delta z$, $\mathcal{L}$, $z_{0}$ and $z_{T}$, we obtain the first result.

Assume that we have the optimal $a(t)$, the condition for $z(T)$ to be the optimal final position is

[TABLE]

Note that $z_{0}$ is fixed and so $\delta z$ in $z_{0}$ is null. We thus obtain the second result of the lemma.

-B Proof of Lemma 2

As $\mathcal{L}(z,a)$ is an homogeneous function of $z$ and $a$, we have: $\mathcal{L}(\lambda z,\lambda a)=|\lambda|^{\alpha}\mathcal{L}(z,a)$ for all $\lambda$ (in our case with $\alpha=2$). Deriving this expression with respect to $\lambda$, setting $\lambda=1$, and noting that $a=\dot{z}$ we obtain

[TABLE]

Using (6) and (29), we have: $z\cdot\frac{dp}{dt}+\dot{z}\cdot p=\alpha\mathcal{L}$ or equivalently $\frac{d(p\cdot z)}{dt}=\alpha\mathcal{L}$. We can now integrate the cost function (2) along the optimal trajectory as follows

[TABLE]

-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$. For simplicity, we omit variables $T$ and $z_{T}$.

[TABLE]

By using the same approach between $t-h$ and $t$, we deduce in the same way equation (9) when the final time $T$ is varying.

-D Proof of Theorem 1

From (5) and (11), we obtain the following ordinary differential equation of second degree: $\ddot{z}=-\frac{u_{0}}{K}z$. If $\frac{u_{0}}{K}>0$, we define $\omega^{2}=\frac{u_{0}}{K}$ and we look for an optimal trajectory of the form: $z(t)=A\cos(\omega t)+B\sin(\omega t)$. If $\frac{u_{0}}{K}<0$, we look for an optimal trajectory of the form: $z(t)=A\cosh(\omega t)+B\sinh(\omega t)$ with $\omega^{2}=-\frac{u_{0}}{K}$. Let us denote $z_{0}=z(t_{0})$ and $a_{0}=a(t_{0})$ the initial conditions for $z$ and $\dot{z}$.

Take the case $\frac{u_{0}}{K}<0$. Using the derivative of $z(t)$ and identifying terms, we obtain: $z(t)=z_{0}\cosh\omega(t-t_{0})+\frac{a_{0}}{\omega}\sinh\omega(t-t_{0})$. At $t=T$, we have also: $z_{T}=z_{0}\cosh\omega(T-t_{0})+\frac{a_{0}}{\omega}\sinh\omega(T-t_{0})$, from which we deduce

[TABLE]

when $u_{0}<0$. In a similar way, we have:

[TABLE]

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 (7), 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 (7), we obtain the result for the cost function.

-E Proof of Theorem 2

We assume that the location and time $(\xi,\tau)$ 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 $(\xi,\tau)$ on the one hand and between $(\xi,\tau)$ and $(z_{T},T)$ on the other hand, by the principle of optimality. In region 1, the optimal cost up to $\tau$ is

[TABLE]

Using Hamilton-Jacobi, we obtain

[TABLE]

In the same way, the optimal cost in region 2 is

[TABLE]

Using again Hamilton-Jacobi, we obtain

[TABLE]

The total cost along the optimal trajectory is the sum of the cost over the two regions

[TABLE]

A necessary condition for the optimality of $\tau$ is thus

[TABLE]

i.e.,

[TABLE]

A necessary condition for the optimality of $\xi$ in (38) under the constraint $f(\xi)=C$ is characterized by

[TABLE]

where $\mu$ is a Lagrange multiplier associated to the constraint and where the second line comes from equation (10) of Hamilton-Jacobi. Thus we obtain precisely equation (20).