Optimal Aerial Guidance in General Wind Fields

This paper presents an optimal guidance approach for a UAV point-to-point navigation in 2D, under wind perturbation, with a general desired airspeed proﬁle. A cost function weighting the travel time and the control eﬀort is minimized through the Pontryagin’s Minimum Principle, involving the derivatives of the airspeed velocities. Two iterative procedures for a guidance algorithm under general wind ﬁelds were developed, including an analytical solution for the optimal heading in minimum-time paths. Diﬀerent cases from the literature are compared and a hard wind scenario test is presented


I. INTRODUCTION
The last decade showed a growing interest in the development of planning and guidance for unmanned aerial vehicles (UAVs), with special emphasis on task and path planning, as well as optimal guidance [1], [2].In a top-down layer architecture, a UAV flight management system usually consists of the mission planning, the path planning and the guidance system [3].Upon receiving a set of waypoints to be reached, the guidance algorithm is responsible for providing the acceleration commands to the low level flight controllers in order to minimize a predefined performance criteria such as the traveled length, the traveled time, or the consumed energy [4].Additionally, to cope with the need for path replanning, it is desirable to have a computationally efficient algorithm in the onboard real-time embedded system [3].In this way, a closedform solution for the optimization procedure of the guidance algorithm is sometimes preferred instead of a numerical one.
In this paper, we consider the design of an optimal "timeeffort" trade-off guidance for a vehicle that flies between two given waypoints in 2D, under the influence of a spacedependent flow field (wind).To solve this problem, we leverage some of the classical approaches used in missile and spacecraft guidance [4], [5] by adding the effect of the wind field.In this way, the vehicle is considered as the "pursuer", while the target waypoint is the "fixed target" that, despite the wind, shall be reached at a given desired speed (rendez-vous problem).
The wind effect is an important nonlinear disturbance for aerial vehicles, with an impact that is proportional to the wind speed to the airspeed ratio.The effect of these perturbations E. C. de Paiva is with the School of Mechanical Engineering, University of Campinas, r.Mendeleyev, 200, Campinas, SP, Brazil, CEP 13083-860, email:elypaiva@fem.unicamp.br.Guilherme Pereira is with the Department of Mechanical and Aerospace Engineering, West Virgina University, 1306 Evansdale Drive, Morgantown, WV, USA, 26506-6106, email:guilherme.pereira@mail.wvu.edu on the flight path is especially relevant for the case of small unmanned aerial vehicles, as well as the lighter-thanair aircrafts case [6], [7].If on one hand, the wind can affect the UAV maneuvering capabilities (like the minimum turning radius), on the other hand, it may also help to save the energy consumption [8]- [10].This last point is one of the objectives of the proposed guidance algorithm, which minimizes a cost function weighting the control effort and the travel time, while benefiting from the energy of the flow field.
In the last years, a number of researchers have investigated the problem of guidance in the presence of wind [2], [9]- [16].However, most articles in the literature are focused on minimum-time problems with constant airspeed, as remarked in [16].Further, in many cases, the optimal trajectories are not feasible by not respecting the aircraft and the path constraints, yielding sharp turns that are not practical.And finally, many solutions rely on heavy computational numerical procedures for the guidance optimization, which may be a problem for the limited processing capacity of small UAVs.
The problem of optimal navigation in a flow field is closely related to the classical Zermelo navigation problem (ZNP) [17], [18].This problem was formulated by the german mathematician Ernst Zermelo, when the airship Graf Zeppelin circumnavigated the Earth, in August 1929.The objective was to find the minimum-time path for a vehicle in the presence of spatiotemporal drift flow field, when navigating between two points [19].The original ZNP problem evolved, along the years, under different extensions and variations, attracting the interest of the investigators [20], [21].One of the limitations of the original ZNP is that the motion of the vehicle is described by a "single integrator" kinematic model, where an ideal "instantaneous" control can be aplied.The inclusion of dynamics in the ZNP problem was first proposed in [22], when the authors included in the problem a self-propelled (or Newtonian) particle, under limited acceleration.Bakolas solved, in this work, the minimum-time problem for a spatially invariant wind [22], and later for an affine spatio-temporal wind [20], but for the intercept case, not for the rendez-vous.Also, in both cases, the numerical solution is computationally intense, the airspeed rate is limited (in norm-2) and the final velocity is free (intercept case), such that one can not impose both the terminal position and terminal velocity of the vehicle, as we do here.
Other relevant works over guidance with wind for aerial vehicles are the following.For the case of a long-distance airplane traveling horizontally between two points, an important research has been presented in [23], with an associated patent in [24].The authors used a neighboring optimal control (NOC) approach to solve the 2D Zermelo-like problem, assuming constant airspeed, finding about 1-3% reduction in the flight time w.r.t. the nominal straight flight.Solutions for the timeoptimal problem with free or constant airspeed flight are found in [9], [14], [25]- [27], as well as solutions based on Pursuit or Dubins [11], [28], and other guidance algorithms based on numerical procedures [10], [12], [15], [29], [30].The classic Zermelo-Markov-Dubins problem has been studied in [28], but again using sharp turns and assuming constant speed.None of these approaches provide a mechanism to extend the method to variable-speed trajectories.
Further, as far as we know, our work is the first analytical "time-energy" optimal solution for the rendez-vous problem including the wind.Recently, a generalized optimal cost solution was published in [31], but for the rendez-vous in space applications, disregarding the atmospheric drag (wind).In our methodology, we first minimize the generalized cost function, through the Pontryagin's Minimum Principle, assuming the case of a position-dependent wind, instead of a time-dependent wind, considered in our previous work [32].Notice that the methodology presented in this paper is more practical, since a time-dependent wind assumption is not useful, even when an iterative procedure is used, as it is necessary to know the future values of the wind speed.The new solution yields also analytical expressions for both the optimal control input and the optimal initial heading.In a second moment, the problem is generalized to a general wind case, through an iterative process that relies on a set of small path segments.Further, a desired airspeed profile may be assigned for the whole path, either by using acceleration correction/updates or by tuning the trade-off parameter of the cost function.
It is important to remark that our proposal is focused on the guidance optimization for navigation on wind between two given points, that can be further evolved to a path planning/guidance strategy.In the literature we find a lot of methods focused on high-level planning of optimal paths considering discretized domains.This is the case of KITE algorithm, used to control a Bell 206 autonomous helicopter [16], [30], [33].This work uses a decoupled trajectory optimization approach that computes feasible, time-optimal trajectories while accounting for wind disturbances.However the procedure requires a large computational effort as it is based on a recending horizon approach, under extensive numerical computation [16], [30], [33].Also, as far as we know, the airspeed profile in KITE is basically one of piecewise constant with smooth transitions between the path segments.
The rest of the paper is organized as follows.In Section II, we introduce the problem formulation.Section III presents the general solution of the optimization problem.Section IV generalizes the solution for the iterative guidance algorithm.Section V shows the simulation results, and section VI presents the final conclusions.

II. PROBLEM FORMULATION
We present here the problem statement related to the optimal 2D trajectory generation of a UAV under wind perturbation, considering the minimization of a time/effort cost function, as shown in Figure 1.When travelling from point A to point B, the objective is to make the UAV reach the second point ("target"), in an optimal way, at a given constrained terminal velocity (rendez-vous problem), assuming a positiondependent wind distribution.For guidance purposes, let the motion of the UAV be described by the simplified point-mass dynamics below, which is a common assumption when the trajectories generated are to be used as reference signals for the low level flight controllers [4], [20]: is the inertial velocity vector and v a (t) = [v ax (t), v ay (t)] T ∈ R 2 is the relative speed, or airspeed vector, that can also be written as v a (t) = V a [cos(θ), sin(θ)] T .The wind speed vector is given by w = [w x , w y ] T ∈ R 2 , and it is supposed to be positiondependent, or w = w(r).The control input u(t) ∈ R 2 is the airspeed rate va (t), such that the optimal path will be driven by the optimal relative acceleration to be found.
Firstly, let us define the "control effort" as the integral of the squared norm of the airspeed rate (the control input), that is: Then, the optimal control problem for the trajectory generation between the two given points, with the generalized performance index J relating the control effort J u and the travel time t f , can be stated as: where C I > 0 is the trade-off coefficient between the costs associated to J u and t f .Thus, with a higher C I , the problem puts more weight on the minimum-time objective, and with a lower C I , the formulation tends to a "minimum control effort" problem.
It is interesting to note that this optimal control problem can also be interpreted as a generalization of the classical "Linear Quadratic Minimum Time" problem (LQMT) for the case of relative velocities.The development of LQMT comes from the decade of 1960, and it was the basis of one of the moon landing guidance approaches, known as E-guidance [34], aiming at minimum-fuel consumption.Although originated in the 1960's, the concepts of LQMT were presented in a full and analytical formulation by the first time in 1991 by Verriest and Lewis [35], and were used in a recent work [36] on UAV path planning to generate smooth trajectories with jerk/acceleration constraints.In this work of [36], using the classical LQMT (without wind), the authors call the control effort J u = u T u as the "smoothness of a trajectory".Such that we generalize this concept, with J u representing now the "smoothness in the airspeed", as the minimization of J u means the search for a solution with reduced variations in the airspeed along the time.
Meanwhile, minimizing the integral-square vT a va has also other possible advantages.First, the search for a less varying airspeed improves the vehicle maneuverability, as the aerodynamic actuation depends on the air velocity.Second, this integral-square criterion allows for the increase in the inertial velocity of the aircraft under increasing winds, as well as its decrease for decelerating winds, reducing airspeed and drag, and thus allowing for energy savings.

III. OPTIMAL SOLUTION FOR THE SIMPLIFIED WIND MODEL
In a first moment, in order to solve the proposed optimal control problem, we will suppose a simplified positiondependent wind profile, where the cross-track (or lateral) wind component is constant, or w y (t) = w yc , and the along-track (or longitudinal) wind component w x (t) is proportional to the lateral y-distance, or w x (t) = ky(t)+w xc , where k is the wind gradient and w xc is a wind bias or offset, as shown in Figure 2.This is a common simplification used in the classical version of the Zermelo problem [17], [18], [23].It could be interpreted as a model dominated by the longitudinal motion, where the cross-track wind w y (t) does not have a strong effect on the minimum-time trajectory solution, and the along-track wind w x (t) varies with the cross-track distance.We remark that, throughout this paper, we call "longitudinal" for the alongtrack distance and "lateral" for the cross-track distance, as our reference frame is the world frame.Further ahead, in this paper, the wind model will be generalized to one of piecewise linear components, through an iterative procedure, whose steps use the basic model solution, presented herein.
The simplified wind model can also be written in matricial form as: where K = 0 k 0 0 , w c = w xc w yc and r(t) = x(t) y(t) .
Such that the wind speed rate is given by ẇ(t) = Kv g (t).
Recalling the generalized performance index J = (J u + C I t f ), where C I > 0 is the trade-off coefficient between the control effort J u and the travel time t f , the optimal control problem for the trajectory generation, with the simplified wind model, can be restated as: This type of problem is closely related to the classical Zermelo problem, for the minimum time navigation with constant airspeed.An implicit analytical solution for a similar problem, has been derived in [17], [18] for the case where the x-direction wind speed is affine on the lateral y-distance (w xc = 0 and w yc = 0).However, for nonzero wind offsets, w xc = 0 and w yc = 0, as in our case, we find only numerical solutions in the literature.

A. Optimal Control Solution
To solve the optimal control problem (P2) above, which has free terminal time and fixed terminal states, we use the Pontryagin's Minimum Principle (PMP), that furnishes a necessary condition for a control input to be optimal.Indeed, as our problem works with linear dynamics and convex input sets, in our case PMP yields a condition that is both necessary and sufficient for the optimality, thus allowing to completely characterize the time-optimal control law.
The Hamiltonian function for our problem is written as where p r and p v are the costate vectors associated with the position and velocity vectors, respectively.A necessary condition for the optimality is that H u = ∂H ∂u = 0.This relation allows to find the optimal control input that minimizes the Hamiltonian, which in our case becomes: Also, recalling that the relation ẇ = Kv g holds for our wind model, we can rewrite the Hamiltonian equation as Such that, using the fact that u = vg − ẇ = vg − Kv g , according to Pontryagin Minimum's Principle (PMP) the costates derivatives for this problem are given by which implies that p r is constant and p v is linear with time.Thus, from ( 6), the velocity costate vector yields the optimal control signal, that is the airspeed rate vector, or whose integration over time yields the expressions for the airspeed velocities as: Such that the expressions for the inertial velocities become: And the expressions for the whole set of positions/velocities components are

B. Evaluation of the Optimal Costates
Now that we have the polynomial expressions for the optimal velocities/positions, we need to evaluate the coefficients of these polynomials, that are related to the optimal costates.With this purpose, we consider the following initial/final conditions: Which are equivalent to the bounding conditions for the inertial velocities as Note that for the equations above, as we assume here a position-dependent wind model, it is possible to evaluate the wind values at the end of the path (t = t f ).As the target position is at the origin, we have w x (t f ) = ky(t f ) + w xc = w xc .And with the help of the iterative procedure, presented in the next section, this will allow for an efficient guidance algorithm.This was not possible with our previous approach [32], where the wind was a pure function of time, unrelated to the position.
Let us define two intermediate variables δ vax , δ vay as: Thus, from the four equations in (11) evaluated for the terminal time t f , we have four equations and four unknowns (p rx , p ry , p vx0 , p vy0 ), whose solution leads to the final analytical expressions for the optimal costates as: Such that the optimal control input u * (t) is finally written, from (8), as a function of these costates as: Note that the coefficients of this linear acceleration command, given by the costates, depend on the value of the terminal travel time (t f ), such that the acceleration increases for small values of t f .This is the reason of the huge increase in acceleration commands when approaching the target, a common problem found in missile guidance, which usually requires the command adaptation for the close near-target maneuvering.We will get back to this point further ahead.
Finally, the particular value of the terminal relative acceleration for t = t f is written as: C. Derivation of the General Time-Optimal Polynomial The optimal terminal time t f can be found by equating the terminal Hamiltonian H(t f ) to zero, which is the so-called transversality condition.Thus, using equation ( 7), we have: The expansion of this last expression will lead to the quartic time-optimal polynomial, whose roots yield the optimal travel time t f .Calling v g (t f ) = [v gxf , v gyf ] T , and using ( 15) in ( 16), after a tedious algebraic manipulation, we finally come to the quartic time-optimal polynomial in t f as: where with z 1 = 6v gx0 + 4δ vax − 3ky 0 z 2 = 6v gy0 + 4δ vay (18) It is important to mention that, in the case of a constant wind speed (k = 0, or null wind gradient), when the coefficient a 3 becomes zero, this optimal polynomial turns out to be the same depressed quartic polynomial appearing in the previous works of [35], [36], [37] and [38], when the wind was not considered.

IV. SUBOPTIMAL ITERATIVE SOLUTION FOR THE GENERAL WIND PROFILE
The optimal solution for the guidance problem derived in the previous section yields the optimal time and optimal accelerations that minimizes a generalized cost function, considering the simplified wind model presented in Figure 2. If a general position-dependent wind profile was to be considered in the previous optimal control problem (P2), it would become extremely complex, and the optimal solution would be hard to be derived.However, if we consider an iterative procedure, where inside each iteration, corresponding to a path segment (Figure 3), we approximate the wind profile by the simplified model (using the results of previous section), then it is possible to solve the guidance problem in a suboptimal way for a general wind configuration.
We analyse in this section two different types of iterative procedures, that will generate the two guidance approaches proposed in this paper.The first one tries to impose constant longitudinal airspeed ( vax = 0) to the vehicle inside the given path segment or iteration.This will lead to LARG-"Longitudinal Airspeed Rate Guidance" approach.The second one tries to impose null airspeed rate ( Va = 0) to the vehicle inside the segment, leading to ARG-"Airspeed Rate Guidance" approach.In the first case (LARG), we also derive an analytical solution for the optimal initial heading inside each segment as to minimize the travel time, yielding the approach called here LARG-h.

A. Null Longitudinal Airspeed Rate (LARG guidance)
Firstly, we analyse the guidance approach where the commanded longitudinal airspeed acceleration is zero, or vax = 0, which will shown to be useful when the path is dominated by a longitudinal movement and longitudinal airspeed.
From the first component of the acceleration vector in equation (8), and using (12) and (13), in order to assure vax = 0, or constant longitudinal airspeed (v ax = cnst), we need to impose: Thus, to obtain a null acceleration in the longitudinal axis, it is sufficient to make p vx0 = 0, as with constant longitudinal airspeed we have δv ax = (v axf − v ax0 ) = 0, or: which is equivalent to the simple second order polynomial in t f : Polynomial P 1 (t f ) in (19) implies that there are at most two possible optimal times with corresponding optimal costates p x = 0, p vx0 = 0, such that the longitudinal airspeed is null along the given iteration or segment, for the simplified wind profile given by w x (t) = ky(t) + w xc and w y (t) = w yc .
The idea here is to find the given optimal travel time t f 0 from the quadratic polynomial in (19), and subtitute this value of t f 0 in the general quartic polynomial (18) to find the corresponding trade-off parameter C I0 for the given condition.This parameter C I0 could then be increased or decreased to find faster or slower alternative solutions.
Finding the optimal initial heading θ * 0 for the minimum time flight with constant v ax (LARG-h) Still considering the constant longitudinal airspeed case (v ax = cnst), it is interesting to investigate the influence of the initial heading angle θ 0 on the airspeed profile and on the total flight time.To illustrate this, Figure 4   We show below that it is possible to find the optimal initial heading θ * 0 for the minimum time flight under constant v ax , through the roots of a special 4th order polynomial in θ 0 .
Let us recall the quadratic polynomial P 1 (t f ), from (19), rewritten as: ) First, if we write the airspeed components as functions of the heading angle and the true airspeed, we have: And as we assume constant v ax along the segment, the following relation holds And supposing that the final heading angle has opposite signal w.r.t. the initial heading, assuming v ay0 = −v af , we have a second relation below (note that v ay0 = v af is also possible, but less useful): Such that the quadratic polynomial above can be rewritten as a function of the initial heading angle, or P 1 (θ 0 ), as ) Thus, if we obtain the partial derivative of P 1 (θ 0 ) with respect to θ 0 , equating it to zero, we have This relation between the optimal terminal time and the optimal heading angle can be substituted back in (24) to yield a full polynomial in θ 0 , that is ) If we use a simplified notation of S = sin(θ 0 ), C = cos(θ 0 ), and denoting α = (3ky 0 + 6w xc ), then after a tedious algebraic manipulation, we come to the quartic polynomial in S = sin(θ 0 ) below ) Thus, given the initial position (x 0 , y 0 ) and airspeed V a0 , as well as the longitudinal wind gradient k and wind offset w xc , the real roots of P 1 (S) above will give the optimal initial heading angle θ * 0 , such that the path from point A to point B, under constant longitudinal airspeed v ax0 = V a0 cos(θ * 0 ) shows the minimum travel time t * f = 6 k tan(θ * 0 ).In the case that P 1 (S) yields more than one real root, the optimal θ * 0 could be distinguished by the root that verifies the relation tan(θ * 0 ).Note that, in [18], the optimal minimum-time solution for a similar problem, but with constant V a instead of constant v ax , and also with y 0 = 0, w xc = w yc = 0, is given by the relation

B. Null Airspeed Rate (ARG guidance)
Here we will consider the problem of searching the solution for a travel path with constant airspeed V a , which is a typical case in the optimal guidance in flow fields, like in the Zermelotype problems, where usually one wants to find the minimumtime path under a given constant relative velocity.
Indeed, the solution will be obtained deriving the optimal time/accelerations, among all the solutions of the optimization problem (P2) in ( 4), such that Va (0) = 0.This means that we will find, from the general quartic time-polynomial in (18), the terminal time t f that will generate the optimal accelerations vax (t), vay (t) such that Va (0) = 0.It can also be interpreted as the search for the direction with which the airspeed tends to be constant at the beginning of the movement.We will show below that the optimal time solution t f for this condition is again given by a simple quadratic time-polynomial.This result will then be used in an iterative procedure to assure constant airspeed during the flight of the vehicle from point A to point B, provided that we use it iteratively, with small segments.
Another important consequence of this special solution is that the path could be further adjusted, through the tuning of C I , during the iterative procedure, to generate alternative paths with increased/decreased airspeeds.As far as we know, optimized paths with wind, under variable airspeeds, are only found in the literature using numerical solutions, and not with an analytical-based formulation as presented here.
The expansion of this last equation yields the polynomial (30) Thus, the positive real root of this 2nd order polynomial in t f , if it exists, will assure that Va (0) = 0.
Note that when we have v ay0 = 0, then P 2 (t f ) will be equivalent to the previous polynomial P 1 (tf ), yielding the same roots.This suggests the remark below.
Remark 4.1 When a given small segment path is dominated by the longitudinal motion, such that v ax v ay , the timeoptimal solution that yields constant longitudinal airspeed v ax for this segment is close to the solution that exhibits constant airspeed V a , as in this case V a ≈ v ax .This means that the optimal heading derivation from (27) can also be used, with reasonable approximation, to guide the vehicle in the search for the minimum-time solution under constant airspeed.

C. Iterative Guidance Algorithms
Remark 4.1 can be illustrated in the following example, shown in Figure 5, even if it does not verify exactly v ax v ay .The dataset for this simulation was: airspeed reference V a = 10, θ 0 = −180 deg, w x (t) = 0.3y(t) − 3 , w y (t) = 0.5, x 0 = 30, y 0 = 20, iteration step time of 0.25 sec.This example also shows how the two iterative procedures (ARG, LARG), presented in the previous subsections, can be used in a guidance algorithm for a general wind profile.Figure 5(a) uses LARG approach in a single step solution (no iterations).Note that the single step solution for ARG is not possible as it assures only "initial" null airspeed rate, forcing it to be used iteratively.Figure 5(b) uses LARG and Figure 5(c) uses ARG, both iteratively.Note that the relative velocities for these last two cases are very similar as the longitudinal motion is dominant here.Look also that v ax is piecewise constant inside each segment of Figure 5(b).In a general wind case this will not be true, as inside each iterative segment usually the real values of the wind offsets w xc and w yc are not constant.The idea for the general guidance algorithm is presented in Algorithm 1 below.This same algorithm structure can be used for all the guidance variants, except for the heading optimization part at the end, used only for LARG-h.Important time steps are δ t , ∆ t1 , ∆ t2 .The first one (δ t ) is the numerical simulation step, the second one is the time interval for the iterative segment (∆ t1 ), after which the wind data is updated, as well as the commanded acceleration for the segment.The last one (∆ t2 ) is the frequency at which the heading corrections are made (for LARG-h), upon the optimal heading computation for the segment.
An important point is the forgetting factor α used at the beginning of the loop.It is important for the guidance variants with heading optimization, LARG-h (α = 0.05 to 0.4) and LARG-h0 (α = 1).Such that LARG-h0 would use the full optimized heading command, while LARG-h would weight this value with the current heading value.This is important to avoid strong heading changes, which in practical applications could eventually saturate the available vehicle actuation or generate non-feasible sharp turns.
Another important point is the use of an airspeed profile (V a prof ile ), when updating the commanded airspeed.Indeed, we adopt this strategy here in order to generate lower/higher airspeed profiles, instead of adjusting the trade-off parameter C I , that may be somehow equivalent.
Algorithm 1 ARG (α = 0), LARG (α = 0), LARG-h0 (α = 1), LARG-h (α = 0.05 to 0.4) Initialise x 0 , y 0 , θ 0 , V a0 , V af , w x0 , w y0 , v gx0 c , v gy0 c , α, δ t , ∆ t1 , ∆ t2 while t f0 > ∆ t1 do v gx0 c ← αv gx0 c + (1 − α)v x old and also ) as to satisfy V a profile Find optimal time t f0 for the segment (from eq. ( 19) for LARG's and eq.( 30) for ARG) if t f0 < ∆ t1 ∆ t1 = t f0 (preparing for the last iteration) end if Find optimal costates and accelerations va (t), vg (t) Simulate segment iteration for interval t = 0 until ∆ t1 , to get v gx (t), v gy (t), x(t), y(t) Update wind, V a0 , θ 0 and variables for the next segment: if t is multiple of ∆ t2 (LARG-h and LARG-h0 ) Find optimal initial heading θ 0 = f (x 0 , y 0 , w x0 , k, V a0 ) (from eq. ( 27)) end if Update v ax0 c = V a0 cos(θ 0 ), v ay0 c = V a0 sin(θ 0 ), and update v gx0 c , v gy0 c end while Other practical aspects of implementation are commented in the following.The requirements of minimum turning radius to avoid sharp turns in the path due to aircraft limitations can be treated indirectly by decreasing the forgetting factor α, such as to yield a solution "more close to ARG", preventing strong heading changes, as commented above.To avoid possible obstacles in the path, or even to enforce flight corridors, "virtual zones" of strong heading winds could be considered.This could be done, for example, using artificial vector fields [39] to create "virtual winds" in specific regions of the path.When using these guidance approaches in a waypoint navigation, careful should be taken during the transition of waypoints (around the target) due to the increase in commanded accelerations with the decrease of t f , as the final distance to the target decreases.This is a typical problem found in missile and spacecraft guidance, as commented previously [4], [31].
Another important practical concern is to constraint the vehicle velocity, acceleration and jerk during the path, due to the operational specifications of the aircraft [16], [30].Possible solutions, for this problem, can be found in [36] that uses a similar approach as ours for the planning/guidance of quadcopters (without wind).The first possibility is to include, in the original optimization problem, a constraint for the terminal time as t f > t v where t v := p 0 v max is the minimum-time to reach the target (at distance p 0 ) under a given maximum velocity v max .Another possibility, less conservative, is to reshape the problem during the iterative procedure, since the low computational cost of the solution allows for that.In this way, since the optimal velocitiesaccelerations-jerks are simple polynomial functions of time, their extrema within the time interval [0, t f ] could be easily calculated and ajusted through the tuning of C I .Finally, in the case that the optimal time-polynomial used does not present a positive real root (which was seldomly observed by us), the terminal time t f used in the previous segment iteration could be adopted.Or even an alternate t f with k = 0 in the quartic polynomial could be calculated, as for the "null wind case" previous results show that there is always a positive real root [37], [38].

V. SIMULATION RESULTS
In this section we present and analyse 6 special cases, where six of them are benchmark problems from four different authors of the literature, and the last one is a navigation flight problem in a hard wind scenario including strong windshifts of up to 90% of the current airspeed.
The first and second cases come from the approach of [23], also related to Patent US6600991B1, from [24], where the minimum-time problem for an airplane traveling horizontally between two points in a wind field is examined.It was formulated as a neighboring optimal control (NOC) problem with a constant cross-track wind shear (gradient) that is generalized for other wind models through an iterative procedure.Note that the wind profile at a given point of the path is calculated by the linear combination of the profile in two reference points (from initial, middle, final points), as can be seen in Figures (6(a), 7(a)).The NOC solution, assuming constant airspeed and using pre-flight nominal data, showed about 1 to 3% reduction in flight time, compared with the nominal straight flight, as can be seen from Table 1.The results for these two cases are shown in Figures 6 and  7.In both figures 6(b,c) and 7(b,c), the simulation dataset used is: airspeed reference of 480 knots, LARG-h with step time of ∆ t2 = 5 min for heading updates, α = 0.05, and ∆ t1 = 5  1, where we used ARG and LARG-h approches (in blue).In both cases, the travel times obtained in [23] have similar values with ARG, which are a little better than LARG-h.However, we should remember that the small perturbation correction of NOC is based on the pre-computed nominal wind/speed of the path, while we use only the local wind data and the destiny point position.The final trajectories are very similar to those found in the cited work.Figure 7 also shows the response for the LARG approach in case 2, with a zooming on two portions of the corresponding longitudinal airspeed curve (v ax ).We would like to remark that zoom 1 shows a constant v ax value for LARG (except for the dynamics transient at the beginning of each segment ), while for zoom 2 this airspeed component is no longer constant.The reason is that in the first case (gray-shaded regions) the wind offsets w xc and w yc are constant, and thus the solution is indeed optimal.Still regarding case 2 with LARG, we show the time evolution of the Hamiltonian H(t) and the trade-off parameter C I in Figure 8.The value of C I is updated at each iteration of ∆ t1 = 5 minutes, according to equation ( 18), to yield the corresponding terminal time t f .Hamiltonian H(t) is numeri-cally close to zero, as expected.When approaching the target, t f decreases and acceleration u increases in Hamiltonian expression (5), and since the wind offsets are not constant in this part of the path, the wind model approximation degrades for this near-target region.
Figure 8. Hamiltonian H(t) and trade-off C I , for case 2 with LARG.When "near-target", t f decreases and acceleration u increases, and since wind offset is not constant at the end, the wind model approximation degrades at this point.
It is worthwhile to say that in addition to the requirement of pre-computed nominal wind data for the path, the NOC approach also needs to perform, at each iteration, coordinate transformations to match the current flight direction with the coordinate system of the nominal NOC gains.Furthermore, their approach works only for constant airspeed, while we can impose a given desired airspeed profile.Because of this, we can even answer a reverse question for this problem, that is: "would there be a reduced quadratic airspeed profile such that the total travel time is the same as in the nominal straight flight?"The answer is shown in Figure 7(d), for case 2 above, using ARG with a quadratic airspeed profile.This figure compares the straight nominal flight (constant airspeed) and the optimized path (ARG) with reduced airspeed, with the same travel time (t f = 831.4min) for both.In the middle of the path, the reduced airspeed case reaches about V a = 456.5 knots, instead of the nominal V a = 480 knots, thus possibly saving energy.
The third case comes from the work of [25] which is an extension of [23].The piecewise-constant solution is examined in the limit as the wind model becomes continuous, leading to an integral feedback law for the neighboring optimal guidance (NOC) [25].One of the examples treated in his work is presented in Figure 9, where the wind speed direction shows an inversion in the middle of the 10 "units long" segment.All the variables here are dimensionless.The simulation dataset is: airspeed of V a = 1, θ 0 = 170 deg, α = 0.05, ∆ t1 = 0.05, LARG-h with step time of ∆ t2 = 0.05, and gradient k varying from −1 to +1 inside wind window.The optimal path solution from [25] is shown in Figure 9(a), and it is very similar to the red path in Figure 9(b), that was obtained with ARG approach.The solution with LARG-h is shown in the same figure, in blue, with a travel time that is slightly greater than that of ARG (t f = 10.58 instead of t f = 10.46).Note that LARGh guides the vehicle toward greater lateral shifts following the better wind gradient at this region.But just after the wind speed shifts direction, the vehicle is then guided in the opposite side, searching for lower heading winds.The UAV stamps in Figure 9 are printed at each 3 seconds.
We were surprised with the excellent behaviour of ARG in Figure 9(b) with good antecipatory turning movement, as "if it knew" what was beyond the barrier in the middle of the region.In [25] this could be explained by the use of nominal wind values for all the segment, but in our case, as ARG uses only local information and does not search for optimal headings like LARG-h, we argued if this could be due to the proximity of the target.With this purpose, we simulated this case in a new configuration where the 10 units long window, with a wind inversion in the middle, is now located far away from the target, as shown in Figure 9(c).The result for LARG-h is somehow similar to the first case, but it is very different now for ARG.In this case, when the vehicle crosses the middle of the region, it only "slows down" the longitudinal motion to cope with the increasing heading wind, and it will only start turning (slightly) after leaving the window, which confirms that the good results of ARG in Figure 9(b) may be explained by the approach maneuver, when reaching the target.
The fourth case is a classic benchmark used to illustrate the Zermelo problem, commonly found in books of Optimal Control like [17] and [18], and also in many papers like [21].Again, all the variables here are dimensionless.As the solution for the general flow field is only obtained numerically, this problem usually appears in the literature with a simple wind profile like w x = ky, w y = 0, which for our wind model corresponds to null wind offsets w xc = w yc = 0, as in Figure 10(a,c).The red arrows in these figures are the heading direction vectors for the ARG solution.For the no wind offset case, a kind of look-up table or graphical solution to obtain the ideal heading along the minimum-time path is given in [18].The particular example shown in Figure 10(a) comes from [21].The dataset used in our simulations was V a = 2, k = −1, ∆ t1 = 0.02, ∆ t2 = 0.05, α = 0.3, x 0 = 7.32, y 0 = −3.727.In part (a) of this figure, we used no wind offsets and θ 0 = 110 deg, like in the original work of [21].In part (b), we used wind offsets and θ 0 = 110 deg.And in part (c) no offsets, but with different initial headings (θ 0 = 110 deg and θ 0 = 210 deg).The optimal minimumtime solution to reach the target, in part (a), is given by t f = 5.46439, according to [21], while we obtained a close value with ARG, that was t f = 5.61.The solution with LARG-h exhibits some oscilatory behaviour at the end due to the target approaching, as commented previously.
Figure 10(b) shows the ARG solution for the same problem but now with a wind offset configuration (w xc = −2.5, w yc = 1.5).In this case, a wind offset blowing up and left helps to reach the target faster.Although this is suboptimal, the solution found in the literature is only obtained numerically, like in previous cited works [18], [21].Figure 10(c) shows the same case of Figure 10(a), without offsets, but using different initial headings.In this case, ARG solution yields a long travel, with a bigger t f , when compared to LARG-h, showing its dependency with respect to the initial heading angle.
The fifth case comes from [26].In this paper, the author proves that in a point-symmetric flow, such as inside vortices or regions of eddie-driven upwelling/downwelling, to obtain a minimum-time solution, the rate of the steering angle (heading in our case) has to be equal to one-half of the instantaneous vertical vorticity [26].We simulated the particular case shown in Figure 11, with a purely rotational and time-invariant flow, with angular rate equal to ω = 0.4.The variables are dimensionless.Figure 11( ≈ ω, which is the condition for the characterization of the optimal solution in this case [26].In this simulation we obtained a travel time of t f = 2.46 and an angular rate of ∆θ ∆t = 1.01 2.46 = 0.41 ≈ 0.4 = ω.This intuitive result for a ship navigation is cited in their work as: "if the fluid rotates as a rigid body, then the fastest approach is to helm the ship at the rate given by the instantaneous fluid angular velocity" [26].Figure 11(b), shows the results of the same problem using LARG-h, for the same 3 initial heading conditions.The optimal terminal time for θ 0 = 144 deg is now t f = 2.54 instead of t f = 2.46.However, the resulting path for the 3 different initial conditions seems to converge to a common similar pattern.The last case is a complex scenario proposed here (Figure 12), including nonuniform wind gradient, strong windshifts of up to 90% of the airspeed and sensor/actuators delays and dynamics.Different tests were made aiming real applications, where we considered that the wind data (wind speed and gradient) are available at every 10 seconds.The wind gradient along the flight is calculated as k ≈ ∆w x ∆y , as if obtained from two different measures close to the vehicle, distant 10 m one from the other (∆y = 10 m).One possibility is to use a double parallel UAV flight, or even a wind forecast.The simulation dataset used here is: airspeed reference of V a =10 m/s, θ 0 = 170 deg, ∆ t1 = 10 sec, ∆ t2 = 10 sec, α = 0.4.We also consider a time constant of τ = 5 seconds in the actuation (acceleration) response, yielding approximately a settling time of t = 20 seconds.In the first test, ARG and LARG-h were compared for a constant airspeed flight (V a = 10 m/s) with a constant cross-track wind of w yc = 2 m/s, as shown in Figure 12(a).Figures 13(a,b) show the velocities responses for LARG-h and Figures 13(c,d) show the velocities responses for ARG, both relative and inertial.The green-shaded regions (I,II,III) in the path correspond to strong windshift transitions, where it is hard to hold the airspeed constant.The UAV stamps in Figure 12 are printed at each 10 sec, which is the rate of wind/acceleration update.The retard in responding to the abrupt transitions (light green areas) is due to the sparse sample of wind update and the actuation delays.This leads to peaks of airspeeds in the relative velocity profiles of Figures 13(a,c).Note the slow response of ARG in the first 100 seconds of simulation (red window in the XY-Path plot), due to an intentionally "bad choice" of the initial heading angle (θ 0 = 170 deg).Instead, LARG-h is able to quickly correct the direction in the first wind update, turning the vehicle left, in the search for better winds.Note also the fast heading corrections of LARG-h along the path, against a smoother path of ARG.The total travel time for ARG here is t f = 433.6 sec against t f = 304.6 sec for LARG-h.However, if we have chosen a better initial heading angle, the difference in the travel time between ARG and LARG-h would not be so large here.
ARG may present other types of shortcomings, as shown in the next three tests for this same scenario.The first test implements a decreased airspeed profile with LARG-h (Figure 14(a)), forcing a travel time close to the t f that would be obtained in the straight flight with ARG, that is t f = 445 sec.The dataset here is ∆ t1 = ∆ t2 = 10, α = 0.4.Figure 14(a) shows the resulting relative velocities for LARG-h in this case.Note that the airspeed in the middle of the path is about half the value of the airspeed in the extremes.The ARG solution for this test yielded a very irregular airspeed profile.
The second test uses a constant airspeed profile, but for a different initial heading (θ 0 = −150 deg), as shown in Figure 14(b).The ARG solution has difficulty to keep the constant airspeed of V a = 10 m/s in the first 100 seconds of travel.Note the interesting behaviour of the lateral velocity v ay , which for ARG seems to be a "smoother", or filtered version of the v ay of LARG-h.In this sense, it is interesting to see that by adapting accelerations to keep constant airspeed naturally helps ARG to find better favour winds.
The third test is a comparison between ARG and LARG-h (Figure 12(b)) for a strong lateral wind perturbation with 8 senoidal cycles of w yc along the path, ranging from -4 to +4 m/s, and using θ 0 = −180 deg.While the travel time for LARG-h is almost the same for the constant w yc case, the travel time with ARG is now around 30 seconds longer, and it has difficulty to keep the constant airspeed in the first 200 seconds of travel, as shown in the same figure (at the lower left part).However, despite the lateral perturbation, LARG-h is able to change heading and to adapt longitudinal acceleration accordingly.The two airspeed peaks for LARG-h around t = 200 sec are due to the strong windshifts in the lower part of the path (green areas), which was not present in the ARG path in this case.

VI. CONCLUSIONS
This paper derived an optimal guidance approach for a UAV point-to-point navigation in 2D under wind perturbation and considering a general desired airspeed profile.A cost function is minimized, through the Pontryagin's Minimum Principle, involving the derivatives of the airspeed velocities, weighting the travel time and the control effort.The problem can be seen as a generalization of the Zermelo problem, as the vehicle is able to track variable airspeeds, while minimizing a generalized performance index, instead of only minimumtime.Two iterative procedures were derived for a guidance algorithm under general wind fields.LARG-h0 guidance is an approach that makes the vehicle to move in the direction of the wind's gradient.ARG guidance is an approach that makes the vehicle adjust its acceleration as to minimize airspeed  variations.The mix of them leads to LARG-h, which combines the advantages of both, and seems to be efficient and robust even in hard wind scenarios.
For the LARG-h case, we also derived an analytical solution, through the roots of a quartic polynomial, to obtain the optimal heading for minimum-time paths, allowing for quick reshaping along the path, as well as the option for paths with increased/decreased accelerations (variable airspeed profile), useful for flight planning/replanning.In the simulations section, we applied the iterative algorithm to 5 benchmark problems found in the literature, from 4 different authors.Also, a challenging scenario case was proposed and solved, with strong windshifts of up to 90% of the current airspeed.
Future works include the development of a waypoint planning algorithm, that uses the present guidance proposal with the inclusion of a mechanism for the safe transition between the waypoints.Also the development of an obstacle avoidance technique, and the generalization for 3D cases are possible.465755/2014-3, FAPESP 2014/50851-0) and Fapesp Auto-VERDE (p.2018/04905-1).

Figure 1 .
Figure 1.Problem formulation for the optimal guidance between two points under general wind perturbation.Wind vector relates airspeed vector to inertial velocity vector.Airspeed heading θ shown w.r.t.world frame.

Figure 2 .
Figure 2. Simplified wind model representation, where the cross-track wind wy(t) is constant, and the along-track wind wx(t) varies linearly with the cross-track distance.

Figure 3 .
Figure 3. Idea of the iterative solution for the general wind profile.The whole path is divided in small segments, inside which the wind behaviour is approximated by the simplified wind model.
shows a particular problem, simulated with 3 different initial heading angle θ 0 for the vehicle (θ 0 = −100 deg, −130 deg, −160 deg), indicated through the arrows of V a0 at the starting point.The dataset for this simulation was: V a0 = V af = 0.5, w xc = w yc = 0, k = 0.5, x 0 = 4.9, y 0 = 1.66.Look in this figure, that although the longitudinal airspeed v ax is constant for the 3 solutions, case (II) yields the minimum-time path, and case (I) could be the "minimum-energy" path due to the reduced airspeed.

Figure 5 .
Figure 5. Illustration of LARG and ARG procedures for the path shown in (a).Relative velocities for the single step solution in (b), and iterative solutions for constant airspeed in (c) with LARG and (d) with ARG.

Figure 7 .
Figure 7. Case 2. Example from [23].(a) Original path results with NOC; (b) path results with ARG, LARG and LARG-h; (c) relative airspeed components; (d) constant x reduced airspeed for ARG and the same t f .

Figure 9 .
Figure 9. Case 3. Example from [25], with wind inversion in the middle of the path.(a) Original path results with NOC-integral; (b) path results with ARG and LARG-h close to the target; (c) path results far from target.
a) shows the resulting path obtained with ARG approach for three different initial heading angles, θ 0 = 120, 144, 190 deg.The dataset here is V a = 2, α = 0.1, ∆ t1 = 0.01, ∆ t2 = 0.01.For θ 0 = 144 deg, the resulting path is approximately an arc of cirle (middle path of Figure 11(a)) and the heading θ(t) is nearly a linear function of time (Figure 11(c)), with an angular rate close to the flow angular rate, or ∆θ ∆t

Table I TRAVEL
TIMES (IN MINUTES) FOR CASES 1 AND 2 FROM [23] VS PROPOSED HEREIN, TOGETHER WITH THE PERCENTAGE REDUCTION IN TRAVEL TIME W.R.T. STRAIGHT/NOMINAL FLIGHT.For case 1, the initial heading is θ 0 = −160 deg, and for case 2, the initial heading is θ 0 = +160 deg.The comparative travel times and the percentage reduction (w.r.t.straight flight) are presented in Table