On Obtaining Energy-Optimal Trajectories for Landing of UAVs

: The optimization issues connected to a landing task of an unmanned aerial vehicle are discussed in the paper, based on a model of a mini-class drone. Three landing scenarios are considered, including minimum-time landing, landing with minimum energy consumption, and planned landing. With the use of classical dynamic programming techniques, including the minimum principle of Pontryagin, as well as the calculus of variations, the optimal altitude reference trajectories are found, to form the altitude control system in such a way as to mimic the proﬁle of the reference trajectory by the actual altitude of the UAV. The simulation results conducted with the use of the Simulink Support Package for Parrot Minidrones verify the correctness and effectiveness of the method, and open the research directions for further analysis, especially to tune altitude controller in a way, as to track the reference proﬁle. Up to this point, optimization tasks considered in the literature, with respect to the drones, were connected to swarm formation optimization, optimization of the take-off process or landing process limited to optimal path planning. This paper thus considers a new topic in the ﬁeld.


Introduction
Recent years have brought an increased popularity of unmanned aerial vehicles (UAVs), which, due to their increasing popularity in the market, have wider and wider applicability areas, ranging from pure recreation to scientific research. The latter implies the need to improve control algorithms which govern behaviour of UAVs to increase their safety and reliability. This can also be understood in the context of energy consumption efficiency.
A variety of tasks carried out by UAVs is a challenge for control engineers [1] responsible for controller tuning of drones. During in flight conditions, the UAVs are required to perform both agile, as well as precise maneuvers, also including maneuvers in formations of UAV, making the issues more complicated [2]. However, during the landing process, the situation is different, and the drones should closely abide some reference trajectories to achieve efficient landing with respect to some performance indices. In addition, at any stage of the flight, it might be necessary to land due to environmental issues, or low battery state, and to execute the flying scenario on the basis of, for example, some state machine.
In 2016, the Department and Aeronautics and Astronautics (AeroAstro) from Massachusetts Institute of Technology (MIT) designed the Simulink Support Package for Parrot Minidrones as Matlab's add-on [3]. The software enables one to design, simulate and test control strategies using real UAVs, and is dedicated to Rolling Spider and Mambo drones of Parrot. It can be also used, as in the current paper, to perform model-in-the-loop simulations, to verify the core idea of the paper, that is, generation of the optimal reference landing trajectory to achieve optimal energy consumption feature. This paper aims at finding optimal landing trajectories for different landing scenarios, as well as focusing on giving directions for modifications of the altitude control system, to tune their controllers in such a way as to achieve optimal performance, which will be a next research stage.
In the literature, one can find multiple approaches to landing trajectory optimization; however, they are usually connected either to landing on a moving target, or approaching the landing spot with a non-zero horizontal velocity. This is not the case here, as this approach can be easily thought of as a pick-and-place stage of the landing/take-off procedure.
Various papers give interesting applications to the drones' deployment, using optimization algorithms. For example, in Reference [4] a Hungarian algorithms is used on the basis of bipartie graphs to match drones to nesting stations, with the increasing interest for smart cities applications. In Reference [5], the authors have used reinforcement learning to obtain the solution, where the agent is able to learn the network topology and infer some information about the environment, in order to find the optimal trajectory that lets the UAV autonomously return to its landing spot within the flying time limit. Similarly, the authors of Reference [6] address trajectory planning issues for an aerial platform to identify and land on a moving car. In order to govern the behaviour of the UAV in real time, the model predictive control is adopted for a generated reference trajectory for the UAV, which are then tracked by the nonlinear feedback controller. The problem of flight path planning of the UAV landing on a moving vessel is discussed in Reference [7], where a method for calculating the optimal approach landing trajectory between an UAV and a small vessel is obtained on the basis of an iterative method. In Reference [8], an optimization task of landing swarms of UAVs is given. The paper is cited not in order to compare the approach presented therein, but to show the optimization has been used to find optimal position of the swarm, when reaching target position, not to generate a vertical take-off and landing-like (VTOL) trajectory of landing.
The author of Reference [9] considers the perched landing, on the contrary to a standard landing on a runaway, of a fixed-wing UAV, where the goal is to deliver the vehicle at point just above the runway surface with near-zero vertical velocity and finite horizontal velocity. During a perched landing, at a certain point in space it is required to have nominally zero vertical velocity and zero forward velocity. The authors of Reference [10] propose the optimization approach to obtain reference landing trajectories in an emergency situation for fixed-wing UAVs, and gives a good reference to the 2009 crash at Hudson River.
The authors of Reference [11] consider take-off trajectory optimization problems for a vertical takeoff and landing (VTOL) unmanned aerial vehicle. As in the case of this paper, the longitudinal kinematics equation and dynamics equation are established with atmospheric density, earth gravity, engine thrust and aerodynamic parameters taken into consideration. The objective function of the optimization problem is the fuel consumption and distance at the given time. Unfortunately, this study is performed for take-off procedure only. An increasing interest in VTOL UAVs, leading to the construction of hybrid structures, see Reference [12], might potentially lead to a further research in the topic at optimization in both take-off and landing stages.
As can be seen, a literature review, reveals that there is a need for the solution to the optimal landing task. The paper is structured as follows-Section 2 gives a simplified mathematical model of a quadrotor fitted to the vertical landing task, Section 3 provides details about model reduction to implement this task efficiently, in Section 4 some basic landing issues are presented, and Sections 5-7 describe minimum-time, minimum-energy, and velocity-penalised landing scenarios, respectively. In Section 8, a comparison of the obtained optimal reference trajectories is given. Section 9 gives an insight into the Simulink Support Package for Parrot Minidrones, to enable model-in-the-loop simulation, presenting the behaviour of the closed-loop system to track given reference trajectories, related to optimal performance. Finally, Section 10 gives a short summary of the work.

Simplified Mathematical Model of a Quadrotor
Let us describe a model of a Parrot Rolling Spider quadrotor driven by DC motors. In order to define position and orientation of this drone, global and local coordinate systems need to be introduced.
The global system is described in NED convention (North-East-Down), and the local system has its origin in the geometrical center of the drone. Its axes, that is, x and y point towards the two adjacent motors, whereas the axis z points downwards. The orientation of the drone can be described using the Euler angles [13] as listed below. The coordinate system of a UAV and the Euler angles are depicted in Figure 1, with: • φ -angle w.r.t. the x axis (roll angle), • θ -angle w.r.t. the y axis (pitch angle), • ψ -angle w.r.t. the z axis (yaw angle). The transformation between global and local coordinates requires one to multiply three rotation matrices R yaw , R pitch and R roll , respectively, to obtain the final transformation matrix R W2B . In the transformations given below, the following notation has been adopted s(α) := sin(α), c(α) := cos(α), and: The backward transformation from local to global coordinates requires the matrix R W2B to be inverted, Similarly, the transformation from local rotational speeds Ω = [p, q, r] T to Euler angle derivativeṡ Θ = [φ,θ,ψ] T is possible when the inverse of Wronskian matrix W is used [13]: According to the second Newton's principle, the acceleration exerted by the unbalanced force F w on a body of a mass m is proportional to this force, and inversely proportional to its mass, that is, Having assumed that F w = G − T, where G is a gravity force vector, and T denotes the thrust vector generated by the propellers, Equation (5) can be extended to the forṁ For the sake of brevity, the arrows above symbols will be henceforth omitted, as only numerical values will be considered. By substituting the matrix R B2W from Equation (2) to Equation (6), it is possible to obtain the acceleration vector in three axes of the local coordinate system of the UAV The thrust T generated by the propellers can be calculated on the basis of formulae below: where: The considered quadrotor (n = 4) has six degrees of freedom, thus, in order to describe its state, 12 state variables are necessary [13], In Equation (11), the symbol x refers to a full state vector, whereas x andẋ are position and velocity of the UAV, respectively, in axis x, and so forth.

Model Reduction
The analysis of the landing process requires the model to be simplified. It is assumed that the UAV does not move in XY plane and keeps a horizontal alignment, thus without any loss of generality it can be assumed that: This assumption is realistic, and present in multiple pick-and-place tasks, as eventually, the control system of the UAV should keep this alignment. When environmental disturbances take place, the UAV is likely to be moved away from the landing position, but when emergency or planned landing is executed, the precision of landing in a spot is not of prime importance. It is a trajectory of descend the most important factor, thus the precise location of landing is not. Abiding Equations (12) and (13) strictly, would stipulate no exogenous disturbances occur. This simplification is somewhat naturel for tilt-hex constructions, see Reference [14] where the accumulative thrust force can be exerted in non-vertical direction, without tilting the structure, to compensate for the disturbances.
It has also been assumed that all four motors have the same rotational speed, leading to simplifying Equations (7) and (10) to the form of: In order to enable any further analysis, the third coordinate system has been introduced, with the origin at the ground level, and the z axis pointing upwards. As the landing takes place while the altitude is not extremally high, it has been assumed that the air density does not vary, then: The symbolic expressions α and u have been defined as (:=): leading to the introduction of a simplified model of dynamics in the z axiṡ in the third coordinate system, with the output equation (state-space equations) y = x 2 . The dynamical model of the UAV can be put in a matrix form (with the time indices omitted for brevity)ẋ with the output equation As can be seen, it is a typical case of a system which is structurally unstable (as g = 0 at all times), resembling many other dynamical structures, as a bicycle robot, subject to the same force, see Reference [15].

Landing Issues
The landing procedure is different from any control task of the UAV during the main phase of flight. The major difference is that reducing the current altitude to the ground level accepts no overshoots, which are typical in any transient processes during regulation tasks. The other issue is that different landing scenarios should be taken into consideration, as the UAV might undergo emergency landing or might land in a planned way. The first one can be given rise due to, for example, low battery levels or due to environmental issues. All the cases imply the landing trajectory h start → h ground , to vary, where h start denotes the initial altitude of the UAV at the time instant t 0 , and h ground denotes the ground level.
The following have been henceforth assumed: where z(t) denotes the position of the UAV in z axis in the third coordinate system. Among the others, it is important to classify landing procedures in the following pattern: velocity-penalized landing, as they arise from typical and emergency scenarios which can be potentially considered. It is to be stressed that the generated reference trajectories should, in principle, be further fed to the altitude control loop to preserve desired behaviour of the UAV.

Minimum-Time Landing
The first scenario considered is the landing in the minimum time, which requires the optimization techniques to be adopted, and is connected to the following cost function: where t f is the duration of the landing procedure (final time), with the assumption that the thrust can be generated upwards only, and the propellers have a limited rotational speed Following a standard routine for minimum-time problems, the following need to be calculated: to verify that the determinants of G 1 , G 2 are non zero for α = 0, and whether the minimum-time control is unique [16]. As known, such a scenario leads to a two-phase control: free falling phase, and braking phase with maximal u = u max [16]. This is a common solution to minimum-time problems arising in control engineering practice, using for example, the minimum principle. The minimum-time control combined of two phases, makes use of the properties of the system, to allow maximum acceleration increase with no energy lost on rotating the propellers (motors stop), and the braking stage. In this sense, it resembles other approaches identified in the field, to use internal system properties to carry on control tasks, see Reference [17]. Below, t s denotes the switching time, t f is the grounding time, and T max = αu max .
First landing phase: Second landing phase: Having assumed that assumptions Equations (25) and (26) hold andż(t s ) = −gt s , the velocity at the first phase of the landing procedure satisfies: and for ζ : The switching t s and grounding t f times satisfy: On the basis of Equations (34) and (36), it is possible to provide the analytical formula for the landing trajectory Sample landing trajectories for the Craig model have been presented in Figure 2 with the parameters of the Parrot Rolling Spider UAV and various initial altitudes:

Minimum-Energy Landing (False-zero Landing)
Minimum-energy landing is connected to the minimization of the performance index taking both altitude and control u effort (energy) into account, with the weight β and leading to the following optimization task: The task has been solved using the minimum principle of Pontryagin by minimizing the following Hamiltonian function [16,18]: By transforming Equation (21) one gets: Using the terminal condition H * (t) = H * (t f ) = 0, and subject to no constraints imposed on u(t): x 1 = e rt , 0 = Ce rt (4) + e rt , 0 = Cr 4 e rt + e rt | : e rt , 0 = Cr 4 + 1, A general formula describing x 1 (t), x 2 (t) is given as: x 1 (t) =e ψt (C 1 c(ψt) + C 2 s(ψt)) + e −ψt (C 3 c(ψt) + C 4 s(ψt)), x 2 (t) = (x 1 (t))dt, x 2 (t) =C 1 e ψt 2ψ 2 (ψc(ψt) + ψs(ψt)) + C 2 e ψt 2ψ 2 (ψs(ψt) − ψc(ψt)) + + C 3 e −ψt 2ψ 2 (−ψc(ψt) + ψs(ψt)) + C 4 e −ψt 2ψ 2 (−ψs(ψt) − ψc(ψt)) , Using assumptions Equations (25) and (26), the integral constants can be found: Finally, the formula describing trajectories depicted in Figure 3 has been obtained. As can be seen, the altitude while landing, falls below zero due to the overshoot phenomenon, what is not acceptable, though definitely shortens the regulation time (typical regulation time vs. damping ratio interplay). In order to compensate for the latter, it has been established that x 2 (t) reaches the ground level at t = π ψ , and the landing phase can be splitted into two phases: where f (t) results from Equation (43) and reduces the altitude between h start → 0. In order to eliminate the impact of the overshoot on the landing phase, the ground level should be virtually raised so that the minimum of Equation (43) equals always zero, what explains the 'false zero' name of the approach. The latter has been shown in Figure 4. The latter implies that the conditions Equations (25) and (26) must be satisfied at all times where: From relations Equation (46) one gets: , ), (47)  It is to be borne in mind that in order to use the trajectory evaluated, the correct value of β must be stipulated first, and subsequently used in the performance index Equation (38). The evaluation of this parameter can be done iteratively by creating, at first, two vectors, which: include a set of initial altitudes, say ranging between 0.5 m up to 5 m with the step of 0.05 m, and the vector of 4000 equally-spaced samples of candidate β values in the range of 0; 3 × 10 −12 . The choice of the final value of the parameter is based on evaluation of various trajectories for every pair of the above-mentioned parameters, according to relation Equation (38). Secondly, the samples with negative thrust force need to be removed, and finally, the value of β is based on the minimal value of J for consecutive initial altitude. The results of such iterative calculations is presented in Figure 5. The function implementing the calculations in Matlab has been presented in Appendix A (Algorithm A1).
Once the final formula for the landing trajectory x 2 (t) is found, and accompanied by the optimal β, it is possible to calculate the trajectories for a specific test conditions. Sample trajectories have been presented in Figure 6.

Velocity-Penalized Landing
The trajectories presented in Sections 5 and 6 can be used during emergency landing, when either time or energy considerations are of prime importance. In these cases, the landing process is fast, and coupled with a huge jerk or snap. However, during the planned landing phase it is useful to take the performance index penalizing the velocity of a drone in the z axis and taking energy considerations into account. During the calculations, a variational calculus approach with constraints is adopted [16], and a penalty taking control signal energy is included into the performance index, to shorten the landing time, It is necessary to transform the cost function Equation (48) using Equations (21) and (22) into: Having assumed that γ := β α , In order to find the extremum of the functional Equation (49), it is necessary to satisfy the following Euler-Poisson equations for all the functions which are presented in Reference [18]: where x i (n) stands for the n-th order derivative of x i with respect to time. Abiding standard Euler notation, the sub-derivative expressions take the following form: for n = 1 one gets δ f δẋ i , for n = 2, δ f δẍ i , and so forth. On the basis of the sub-integral function in Equation (49), a Lagrangean must be defined at first. The term γg 2 has not been included in the performance index, as it is a constant value and has no impact on the precise localisation of the extremum of the functional. In addition, the Lagrange multiplier λ is included into the formula and results from the assumption Equation (27) x 2 ≥ 0 ⇒ −x 2 + κ 2 = 0: In order to formulate the Euler-Poisson equation, it is necessary to calculate the following derivatives of L: On the basis of the latter, and Equation (50), the final set the necessary conditions of optimality is formulated: Two cases need to be considered here, namely λ = 0 and λ = 0: Case I -constraint Equation (27) is active (z = 0): Case II -constraint Equation (27) is inactive (z > 0): γ , By using Equation (22), a time derivative of Equation (54) can be calculated: ) , On the basis of the expressions from Equations (25), (26), (54) and (56), the integral constants C 1÷4 can be obtained: As in the case of minimum-energy landing, β should also be optimized using, preferably, the same method, though in this case β ∈ 0; 2 × 10 −6 , and the cost function takes the form of Equation (48). The listing of a function implementing this idea is given in Appendix A (Algorithm A1). Figure 7 presents velocity-penalized landing trajectories for various initial altitudes.

Comparison of the Obtained Trajectories
The trajectories obtained in Sections 5-7 have been presented in Figure 8. As can be seen, the time of the landing phase and velocities differ as the function of the selected landing scenario. Having assumed that landing is terminated whenever the difference between the ground level and the current altitude does not exceed 0.01 m, the results have been given in Table 1. The value 0.01 m is arbitrary, and is reflected by the way the real Rolling Spider UAV is affected by disturbances, resulting from the air flow, and poor precision of its inertial measurement unit (IMU).
The methodology behind selecting the gain for the integral term in the PID controller is depicted in Figure 9. By comparing the experiments for various values of the I term gain, the value of 0.01 was selected to mimic relatively fast response time in hover state, while keeping good damping ratio of the control loop. It is a common practice in UAV that controllers actually act like PD controllers, when the I term has low gains.

Simulink Support Package for Parrot Minidrones
By using the Simulink Support Package for Parrot Minidrones in Matlab is was possible to simulate the behaviour of the drone in safe conditions, and to verify the obtained optimization results. The software enables one to obtain access to all the signals that would be available from the real machine, with the block diagram in the Simulink given in Figure 10.  The basic block diagram has been modified so that it is possible to provide the reference trajectory, see Figure 11. Its block diagram depicted using standard control blocks, is presented in Figure 10c. The position signal is connected here to the altitude component which is fed back from the estimator, see Reference [19]. This model, as well as many other mathematical models have been made available to the public by Mathematical Model Database, see Reference [20]. As the Simulink model is implemented in NED convention, the sign of the z coordinate in the Simulink block diagram has to be changed to the opposite.  . Building a signal generator to provide reference to the drone.
The control system provided by the developer contains the altitude, XY position, yaw, roll and pitch controllers. According to Equation (7) the translation in z axis in independent from the movement in the horizontal plane, thus the modification of the altitude controller is possible without any changes in the remaining controllers.
In order to eliminate the steady-state error, an integral action has been added to the original PD controller, with the integral gain set as K I = 0.01, see Figure 12, where the controller is given by (e(t) = z ref (t) − z(t), as the tracking error, and z ref (t) as the reference signal), u(t) = k P e(t) + k I e(t)dt + k D de dt . In this configuration of the controller, P and D gains of the altitude controller have been found by a systematic search within the admissible values' ranges, repeating the experiment to ascend to 2 m, hover, and descend the drone starting at 6 s from the altitude of 2 m to 0.5 m, to observe possible overshoot effects. The gain values have been given in Table 2. The corresponding trajectories from the closed-loop system are presented in Figure 13.

Summary
In this paper, we demonstrated the applicability of the idea of optimization of the landing process of the UAV to the model-in-the-loop simulation of the Rolling Spider drone. The encouraging performance of the model of the system for various optimization criteria is shown through simulations. The obtained results imply that it is possible to design the control loop in such a way, so that various frequency responses of the closed-loop systems can be obtained, to shape the system's response to a specific altitude reference signal. In order to do the latter, the altitude controller parameters need to be identified to ensure optimal-like behaviour.
In the future work, the authors plan to move from model-in-the-loop to hardware-in-the-loop experiments, using a motion capture system to precisely measure the altitude of the Rolling Spider drone, in comparison to the reference altitude. Such a knowledge enables one conduct further experiments with more complex approaches, including auto-tuning for real UAV experiment. Carrying out such experiments would make the method portable in the sense that it could also work for other types of UAVs, as just a reference profile and the tuning method need to be used.
On the other hand, it is to be stressed that for unstable structures, as presented in Reference [21], it is possible to build a composite feedback linearization control scheme where using an inverse dynamic model a reference trajectory generator is built, and the polynomial approximation of the trajectory is utilized, what calls for multivariable nonlinear programming problems to arise, to fit the polynomial into the waypoints. This might also be an interesting step, requiring good identification of the model, at first. As a wider and wider spectrum of the UAVs being constructed and successfully deployed, Reference [22], including a spectrum of control approaches, potentially enables researchers to identify new, optimal, control approaches to increase the performance of the flying structures at various stages of flight, including optimal landing procedure.