Autonomous parafoil precision landing using convex real-time optimized guidance and control

To overcome the limitations of current parafoil precision landing capabilities, an efficient real-time convex optimized guidance and control strategy is presented. Successive convexification of the parafoil guidance problem guarantees local optimality with polynomial convergence rate for efficient real-time implementation, where each iteration is dynamically feasible. Our approach shows reliable and fast numerical convergence through in-flight recalculation of time of flight and a new optimal trajectory to cope with time-varying dynamics. The efficiency of our strategy is demonstrated via a comparative analysis of the existing X-38 in-flight demonstrated guidance and control system. Exhaustive Monte-Carlo simulations show performance improvements of about one order of magnitude. The concept proposed is simple, yet general, as it scales to any atmospheric parafoil landing system and allows efficient implementation relying only on the turn rate information.


A. Background
The increased availability of onboard computation power and efficient optimization algorithms associated with automated code generators [1] have dramatically boosted the industrial propagation of embedded optimized guidance and control strategies. Initiated by the research of Açikmeşe, numerous aerospace problems can be solved in realtime [2]. With SpaceX and Blackmore, it was possible to land spaceships accurately [3,4] using that strategy. Convex optimization is widely used in real-time trajectory optimization problems thanks to the maturity of the technology [5,6] and the convergence speed to a global solution. The successive convexification (SCVx) method used in [7] for instance, has been proven to be superlinearly convergent [8].
The successes of the adoption of embedded onboard optimization relies on the numerical efficiencies of powerful algorithms and methods such as interior-point methods (IPM) for Second-Order Cone Programs (SOCP) and Semi-Definite Programming (SDP). For further overview, the reader is referenced to [Yabar2020].
Originally, parafoil guidance has been performed using a sequence of open-loop trajectories relying on the concept of energy management. The various phases of flight are triggered by a trade-off between potential and kinetic energy indicators guiding sequentially the parafoil to its landing site. While this method is reliable and fast, it fails in taking into consideration the uncertainties of the parafoil and wind model, leading to poor landing precision performances.
When the landing precision is a driving factor, especially in unknown and rugged terrain, such as for scientific payload delivery for exploration, or fairing recovery, the former technologies do not meet the requirements. To overcome these limitations, this paper demonstrates the potentials of using onboard optimized trajectory generation and control to increase drastically performance figures.
Driven by the need to embed the optimization onboard, series of authors have approached the problem from a different angle. To start with, the algorithm introduced by Slegers in [9] uses a Model Predictive Control (MPC) on a 2-Degrees Of Freedom (DOF) model for the terminal guidance phase of the landing trajectory. The method can update the reference with an high refresh rate mitigating the effect of the highly unknown wind, which is assumed constant in the optimization process.
In [10], the authors use Particle Swarm Optimisation (PSO) to generate a trajectory online with an update rate of 1Hz. It can generate smooth trajectories thanks to the penalty on the heading acceleration, it shows good performances when compared against other real-time methods. The idea exposed on [11] embodies the actuator limitation in the optimization problem which allows both a refresh frequency of 1Hz and a guarantee of non-saturation of actuators.
The innovation of the guidance proposed in [12] is the use of Euler's elastica functions to get the trajectory that minimises the maximum curvature, taking advantage of the similarity of the equations of motion with a structure problem. Sacrificing a bit of precision, the algorithm run in real-time. The authors of [13] decide to tackle the wind uncertainty problem thanks to a new frame called the WFF (wind fixed frame).
The parafoil model used in [14] takes into account that the horizontal and vertical velocity are influenced by the altitude through the air density, leading to more representative results. The paper brings a method to take into account non-convex and difficult terrain while it ignores the winch dynamics.
In terms of collision avoidance, [15] presents an Rapidly-exploring Random Tree (RRT) algorithm to avoid the especially difficult urban environment.
Finally, [16] uses a more realistic model compared to aforementionned methods. As a result, it takes into account the flared landing. However, the real-time implementation is not specified.

B. Contribution
Onboard optimization for parafoil guidance and control was considered by Slegers in [17]. While the linear MPC offers convergence and recursive feasibility guarantees, the low precision of the linear model used leads to a largely suboptimal solution, and lower landing precision performances. On the other hand, the work of Rademacher in [18] considers a more accurate model, including more physical parameters of the parafoil. It shows better results in terms of landing precision. Besides, it considers a realistic wind profile. However, the performances depends largely on the precision of estimation of the aerodynamic coefficients of the parafoil model. A small error on each parameter adds up to large landing precision error as shown in the Monte-Carlo simulations. Besides, the convergence guarantees are not clear for an arbitrary initial guess. The contribution of this paper is to provide a method with fast convergence guarantees with a 4-DOF model, which is intermediately accurate model with respect to the models used in the literature. In [19] the authors use a sequential convex approach as well but relies on a 6-DOF model. While the resulting trajectory will be easier to track, each iterate is not guaranteed to be dynamically feasible. However, in our approach, since we use a 4-DOF model, the aerodynamic coefficients do not have to be known, and the method is fully transparent to parameters knowledge and estimation. The two main assumptions are the profile of air density and the level of saturation of the actuators. Overall, it enables the use of onboard optimization for a larger set of missions. It assumes the landing cone is free of any obstacle, ignoring any unevenness on the landing area. The landing cone is defined by the set of position where a dynamical feasible trajectory exist, leading to the landing site.

C. Outline
In the section II the parafoil landing problem is introduced. It contains the parafoil and wind model used and the optimization problem to be solved and it provides an estimation of the speed and time of flight. Then, in the section III, a resolution of the aforementioned parafoil landing problem is presented. The resolution has the possibility to be implemented in real-time since each iterate consist of a convex problem and yields feasible solution to the original nonconvex parafoil problem. The possibility to generate a new reference trajectory alleviates the effect of the uncertainties and unmodelled dynamics on the landing precision. Namely, the algorithm presented allows the continuous nonconvex problem to be discretized, convexified and linearized for real-time resolution. A simple controller tracking the generated reference trajectory is finally introduced. In the following section IV, the performances are assessed with a Monte-Carlo simulation, as well as compared to the guidance of the X-38 mission.

II. Problem formulation
A. Parafoil and environment modeling

Parafoil
Modeling the dynamical system, the parafoil, constitutes the first step of the problem formulation and there exists several ways to do it. First, we can use the first principles of physics and consider methodically all the forces and torques acting on the payload and parafoil, possibly considering highly nonlinear effects such as aeroelasticity. This classic modeling is accurate but contains many parameters and takes long to optimize over, which is redibitory for onboard optimization. It is commonly used as a benchmark model to verify the numerical performances of the guidance algorithm offline. Alternatively, we can base our model on the empirical observation that the parafoil rate of fall and the horizontal speed are roughly constant for a given air density, leading to a 4-DOF model. This kinematics model is illustrated on the Fig.1 and as in [20], it depends on the wind strength . This model is representative of the global parafoil motion and is therefore useful for guidance purposes. However, the 4-DOF model does not allow optimization for flare maneuver and generally ignores ground effect. Those dynamical limitations of the model will result in the sub-optimality of the solution of the parafoil landing problem. As it will be seen later, using a simpler model allows to run the optimization online. Both models, the one based on forces and torques (6-DOF) and and the overal motion (4-DOF) are compared in the Table 1.

Model 4-DOF 6-DOF Common use Trajectory optimization Simulation
Limitations Roll and pitch angle Flexible modes Nature Kinematics Dynamics Table 1 Comparison between 4-DOF and 6-DOF parafoil model Choosing and as respectively the position and the wind speed in the j-direction, as the heading angle and as the altitude, the parafoil model used for the optimization is the following The control input ( ) corresponds to an asymmetrical deflection of the parafoil trailing edges.

Environment
A precise model of the environment, such as the wind and the air density, will improve the representativeness of the problem and naturally improve the quality of the solution. In the optimization parafoil model of the Eq.(1), the horizontal and vertical velocities, depend on the air density modeled as a function of the altitude as where ℎ = 1.225 kg m −3 , = 2.256 × 10 −5 m −1 and = 4.2559 [20].
As for the wind, it is generated thanks to a Dryden filter, similarly to [21] as shown on terms; a steady-state profile , a low-frequency/high-amplitude gust , a high-frequency/low-amplitude gust . There are parametrized with the altitude and the ground speed . . The resulting wind model used is To recover the wind in both horizontal directions, we use two distinct Dryden filters with different seed for the white noise [21]. We consider the vertical wind speed to be negligible compared to the vertical speed of the parafoil.
However, it could be easily integrated inside the time of flight estimation.

B. Optimization problem
The parafoil landing problem consists in finding an optimal time-varying reference control input ( ) for the asymmetrical deflection, leading the parafoil towards the landing site with a desired heading angle, minimizing the quadratic control cost for higher control margins. The desired final state x depends on the mission. The objective where the parameters , and are derived in the following section (II.C). We define As we seek to minimize both the control energy and the missed final state x , the parafoil landing problem is multi-objective. The vector = [ 1 2 ] is a numerical artifact to select a single point of the pareto front of the multi-objective optimization problem.
This way, this selects a trajectory that reaches x at the lowest control cost. If the dynamics of the parafoil does not allow it to reach x , it will minimize the missed distance x( ) − x , at whatever control cost. It is worth highlighting * The . * is used to disambiguate the solution of the problem from its variables. that the optimization problem will always have a dynamically feasible solution, since the final state appears only in the objective function. However, it can be difficult to obtain because of the non-convexity of the constraint (4b), preventing the use of a convex solver. The procedure to convexify those constraints is outlined in the section III.A.1. Finally, the optimization problem is continuous and cannot be understood by a numerical solver. The discretization procedure is explained in the section III.A.2.

Speed
The 2 most important parameters of the model (1) are the horizontal and vertical speed . Because the density of the air decreases with the altitude, both and will be attitude dependent and considering this variation will improve the fidelity of the model and subsequently the landing precision. To determine their functions of the altitude, respectively ( ) and ( ), we assume the lift force constant and equal to the weight along the trajectory and the quantity .2 will remain unchanged [15] where . = √ 2 + 2 . Based on the measure of the speeds and density at the initial altitude 0 , respectively, ( 0 ), ( 0 ) and ( 0 ), we have the value of the horizontal and vertical speed at any altitude

Altitude evolution
Because the vertical speed is now known at any altitude, there is a direct relation between the altitude and the time elapsed. Thanks to this relation, we will be able to transform the control ( ) as a function of the altitude ( ).
This is an asset since the wind is considered known as a function of the altitude and it will contribute heavily to the dynamics of the system. We can define the following change of variables The relation can also be inverted to recover the corresponding altitude as a function of time. The idea to exchange the time and altitude variable is also used in [10].
Moreover, we can evaluate this function for the final altitude to get an estimation of the fixed final time of the trajectory.
The final time depends only on the current air density and vertical speed, which can both be measured in-flight.
Having a variable final time would constitute another non-convexity to get rid of, as shown in [22].

A. Approach
The parafoil landing problem as introduced in the Eq.(4)is non-convex because of the trigonometric functions. The general idea of the resolution will be to substitute the non-convex functions by constrained linear functions. As such, to exploit the convergence speed of the resolution of convex problems, we will solve the problem via a Sequential Convex Programming (SCP) algorithm similarly as in [23]. For every iteration, we solve a convexified version of the problem, linearized around the previous solution. We initialize with a dynamically feasible zero-input trajectory and the problem is then sequentially solved until convergence. At each iteration, the solution provided by the algorithm is dynamically feasible and this algorithm is guaranteed to converge to a numerical feasible solution [8]. Because each iteration is dynamically feasible and each iteration is numerically feasible, the output of the algorithm is always a dynamically feasible trajectory. Besides, since each iteration of the algorithm is dynamically feasible, we can choose to stop the convergence process at any given time and output a suboptimal solution to comply with the real-time constraints. Each iteration is a convex problem, which can be solved in polynomial time. The relative value of the cost function between two consecutive iterations is chosen as a stopping criteria for the iterative algorithm. However, the parafoil landing problem is continuous and its solution has an infinite dimension. To be computationally tractable, it is finally shown how to discretize the problem by constraining the input function ( ) to be piece-wise linear.

Input substitution
To form the convexified problem, we introduce a change of variables and a quadratic constraint to remove the inherently non-convex trigonometric functions. A similar method was used in [13].
Likewise, we define a new state vector x ( ) = ( ) ( ) and the heading angle disappeared from the initial state vector. After a few algebraic operations, the parafoil landing problem (4) is equivalent to where the parameters u 0 , x 0 , x and u are inferred from the problem (4).
The solution of the Problem 4 { * ( ), x * ( )} can be recovered from the solution of the Problem (12) {u * ( ), x * ( )} The estimation of the wind profile w( ) goes beyond the scope of this article. As emphasized in [13], it is important to consider the wind proactively rather than reactively to achieve better landing precision performances.

Discretization
Standard numerical tools such as [24] require the variables of the optimization problem to be real numbers. Indeed, searching on the set of all possible functions would be intractable. To convert the continuous problem into a discrete counterpart, we define a temporal mesh and impose the constraints only on its nodes. Motivated by [25], we choose to constrain the input function to the following format which interpolates linearly the value between u ( ) and u ( +1 ) with Although Δ was selected constant for the simulations, we could select a more dense temporal mesh closer to the landing point. The drastic reduction of the dimension of the problem introduced by this discretization results in an approximated solution. Given the nature of the problem, its solution is smooth and well fitted by a piece-wise linear function provided that the mesh is dense enough. In practice, increasing iteratively the number of nodes, the quality of the solution did not improve further beyond 40 nodes, indicating the convergence to the continuous problem's solution for a fixed final time.
In the following, for the sake of simplicity, the subscript will refer to the corresponding time value, such as, x = x ( ) and u = u ( ). Using this input format, we can derive an approximation of the objective as a quadratic function from Eq.(12a).
. Likewise, we have an approximation of the rate of change of the heading angle as a quadratic function using the substituted input variables. The constraints (12d) becomes We discretize the dynamics without loss of precision using the state-transition matrix [23]. In the time-invariant case, we can recover its explicit formulation. With piece-wise affine input from (14), the dynamics constraint (12b) can be expressed as where represents the integrated influence of the wind between and +1 and the state transition matrices are

Linearization of constraints
For the problem to be convex, the last step is to linearize the equality constraint (11). Using a first-order truncated Taylor series expansion around a known previous iterationū( ) gives an approximated linear constraint Since this approximation is valid only whenū( ) ≈ u ( ), we have to impose u −ū 2 ≤ for each iteration and collocation point with a small enough .

B. Algorithm and convex programming
Based on all the previous transformations, we can transform the Problem (12) into a discrete convex problem solved for each iteration. At the iteration (n+1), we solve The equality constraint (20) have been replaced by two inequalities (21d) and (21d) with ℎ close to 0, to give the solver more flexibility to converge. The parameter constrain the gap between two consecutive iterations, trading off between the quality of the linearization and the convergence speed. The cost function is evaluated between each iteration and the convergence of the first SCP is assumed when its variation reaches a chosen threshold. At that stage, we switch to the second SCP of the algorithm where the parameter ℎ becomes a variable to be minimized and is added to the cost function. At the iteration (n+1), we have Overall, the method used stems from [2], where we have two Sequential Convex Programming, the first one giving an initial solution to the second, as proposed in [26]. The first SCP will yield a solution with a small dynamical infeasibility of ℎ . The second SCP will include that value inside the cost function, so that the dynamical infeasibility is minimized. Implementing the second SCP is not essential and is only meant to improve the solution slightly. In practice, the selection of ℎ is driven by the following trade-off. When ℎ is small, the convergence of the SCP is slower but the resulting solution is more accurate. Oppositely, when ℎ is larger, less iterations will be needed and the resulting solution will be more dynamically accurate. In practice, the second SCP converges usually after a single iteration. In order to accelerate the convergence, we can also infer an initial trajectory based on a previous computed optimal trajectory, where the initial state is updated. The algorithm is summarised below. To start the iterative process, we choose a constant initial guess for u ( ), propagating from the initial state, such as ] and we derive x accordingly with Eq. (18). Each iteration, including the initial guess will yield a feasible trajectory, oppositely to [23], where an artificial input is included to deal with the artificial infeasabilities.
The method presented here avoid using the artificial input from [23], by tailoring the sCVX algorithm to the more specific parafoil landing problem. To reduce the number of iterations required to converge, further development could focus on a better feasible initial guess. On Fig.3, we see that each iteration get the resulting trajectory closer to the targeted final state in x = [0, 0] with the heading angle aligned with the X-axis. The variation of state and input and also converges. The shape of the resulting trajectory is not obvious since it takes into account the wind perturbation. Because of the choice of discretization, the asymmetrical deflection to be applied on the parafoil * ( ) is piecewise-constant. The Fig.4 shows an optimal trajectory with the associated constraints and converging total cost.

C. Controller
The controller is needed to compensate for the modeling errors and approximations and unmodelled perturbations and noise. The algorithm above gives a reference for the input * ( ) and for the states x * ( ), used as a feedforward to the controller illustrated on the Fig.5. It operates independently on both the lateral and longitudinal channel to ensure that the parafoil has respectively a minimum lateral and longitudinal error. The symmetrical and asymmetrical deflection of the parafoil controls the rate of fall and the rate of turn of the parafoil. The longitudinal control compensates for the error in the estimation of the time of flight and rate of fall. We assume that the parafoil navigation providex. Inspired by [18], for the lateral control, we have the following asymmetrical deflection and for the longitudinal control, we have the symmetrical deflection

Fig. 5 Architecture of the guidance and control algorithm
The control algorithm is implemented only for illustrative purposes and further development could focus on an optimally tuned controller. 1 and 2 are controller respectively on the cross-track error and the along-track error.

IV. Results
To assess the performances of the algorithm, a 6-DOF model developed by Rademacher [18] was used as a benchmark of the physical reality, including the saturations on the actuators. The 4-DOF model used in Eq.(1) was developed independently.

A. Monte-Carlo
The 600 runs of the Monte-Carlo simulation are used to show the good convergence properties of the algorithm exposed for a large set of initial conditions. On the Fig.6, we can see that all trajectories have led the parafoil close to the landing site. In this case, a single reference trajectory have been computed and tracked with the aforementioned controller. In practice, better results can be achieved when a new trajectory is computed when the parafoil is too far from the reference. Additionally, we can see on the right hand side that the algorithm is consistent. The mean missed distance and mean missed heading angle reaches a steady value after only a few simulations, indicating that the algorithm does not produce outliers. On the Fig. 7, we can see that regardless of the initial conditions in heading angle and position, a trajectory satisfying the final condition is found in less than 35 iterations. We can also point out that each iteration takes on average 0.004 s to be computed on a regular laptop, i7 16Go RAM.

B. Comparison
The online trajectory generation methods is compared against a more traditional method; the energy management.
The energy management method is tested for a 4-DOF X-38 parafoil model calibrated with flight data while the aforementioned algorithm is used on a comparable 6-DOF parafoil model, with the same wind model in both cases.
With similar initial conditions, we can observe the difference in landing precision on the Fig. 8. The final heading angle constraint has been removed for higher consistency and comparability, since the EM method does not target a specific heading angle. The dispersion of the landing precision is much smaller for the online optimization method.

V. Conclusion
The parafoil guidance problem have been solved by the application of a sequential convex programming. Using a substitution for the trigonometric function and linearization, it is possible to convert the non-convex problem into a sequence of convex problems, suitable for onboard implementation and real-time use. Especially, each iteration yields a feasible solution of the physical problem in polynomial time and can be used as a sub-optimal feedforward trajectory, to comply with real-time constraints. The algorithm takes into account the influence of the altitude-varying density of the air on the aerodynamics and requires only the control saturation level for the parafoil model. It also considers the wind distribution and find a locally optimal solution, minimizing the control energy or the missed landing or both when possible. The algorithm is tested on a 6-DOF parafoil model. The navigation error is not included in the simulation.
Further test will include hardware-in-the-loop simulation.