# Multiaircraft Optimal 4D Trajectory Planning Using Logical Constraints

Authors: Dinesh B. Seenivasan; Alberto Olivares; Ernesto Staffetti
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

This paper studies the trajectory planning problem for multiple aircraft with logical constraints in disjunctive form which arise in modeling passage through waypoints, distance-based and time-based separation constraints, decision-making processes, conflict resolution policies, no-fly zones, or obstacle or storm avoidance. Enforcing separation between aircraft, passage through waypoints, and obstacle avoidance is especially demanding in terms of modeling efforts. Indeed, in general, separation constraints require the introduction of auxiliary integer variables in the model; for passage constraints, a multiphase optimal control approach is used, and for obstacle avoidance constraints, geometric approximations of the obstacles are introduced. Multiple phases increase model complexity, and the presence of integer variables in the model has the drawback of combinatorial complexity of the corresponding mixed-integer optimal control problem. In this paper, an embedding approach is employed to transform logical constraints in disjunctive form into inequality and equality constraints which involve only continuous auxiliary variables. In this way, the optimal control problem with logical constraints is converted into a smooth optimal control problem which is solved using traditional techniques, thereby reducing the computational complexity of finding the solution. The effectiveness of the approach is demonstrated through several numerical experiments by computing the optimal trajectories of multiple aircraft in converging and intersecting arrival routes with time-based separation constraints, distance-based separation constraints, and operational constraints.


---

## Body

## 1. Introduction

In air traffic management (ATM), the flight of several aircraft can be modeled as a hybrid dynamical system, which can be regarded as a set of interacting continuous dynamical systems. A number of frameworks have been proposed to model hybrid dynamical systems, in which, in general, differential equations describe the dynamics of each system, whereas logical constraints describe the behavior of the systems during the interactions among them and the interaction with the environment in which they operate. In the ATM context, logical constraints describe, for instance, policies to apply in conflict detection and resolution and operational constraints to be fulfilled during flight. The main operational constraints to be fulfilled during flight are separation constraints, keep-out constraints to avoid no-fly zones, and passage constraints through or by waypoints [1].Given a set of aircraft, separation constraints between them can be expressed as follows: pairwise, they must keep a vertical distance greater than a minimum vertical safety distance or a horizontal distance greater than a minimum horizontal safety distance. The minimum horizontal separation distance can be fixed or variable. In the latter case, it can be established based on the turbulence generated by the preceding aircraft and the ability of the following aircraft to resist turbulence [2]. Obstacles and no-fly zones are in general polyhedral regions of airspace. However, the corresponding keep-out constraints are usually introduced by bounding ellipsoids around obstacles. In some cases, this is a coarse approximation. Keep-out constraints from a polyhedral region of airspace can be expressed as follows: each aircraft must stay outside one of the half-spaces defined by the planes that supports the faces of the polyhedron. This method for modeling keep-out constraints from obstacles can be extended to model passage constraints through windows or waypoints in the airspace. In this manner, the multiphase modeling of the problem is avoided. This is an interesting possibility, since multiphase optimal control models imply the introduction of additional variables to the problem, such as the duration of the phases and additional constraints such as the linkage constraints between them to enforce continuity of the state variables between contiguous phases [3].It is easy to see that all the constraints mentioned above are expressed in disjunctive form. Standard modeling techniques are able to tackle constraints in disjunctive form using binary variables. The four-dimensional (4D) trajectory planning problem for multiple aircraft with logical constraints in disjunctive form can be solved as an optimal control problem (OCP) for a hybrid dynamical system and a common approach for solving this class of problems is to formulate them as a mixed-integer programming problem. In [4], the optimal cooperative three-dimensional (3D) conflict resolution problem among multiple aircraft has been solved in which separation constraints among aircraft expressed in disjunctive form have been included in the model using continuous auxiliary variables. In [5], the optimal path planning problem for multiple unmanned aerial vehicles in the horizontal plane with collision avoidance has been studied, in which constraints for collision avoidance with rectangular obstacles expressed in disjunctive form are included in the model using continuous auxiliary variables. In [6], the trajectory optimization problem for multiple aircraft landing on a single runway in the presence of constraints on the airspace has been treated. The constraints considered are passage constraints through windows in the airspace, and optimal trajectories have been determined by solving a nonsequential constrained multiple-phase optimal control problem.In this paper, the embedding technique proposed in [5] to model rectangular obstacle avoidance in the horizontal plane has been extended to model time-based and distance-based separation constraints and passage through waypoints constraints in trajectory optimization for multiple aircraft. The modeling of passage constraints through waypoints has been done by defining vertical walls in the airspace with a cuboidal window around the waypoint. In this way, introducing multiple phases in the model to enforce passage through waypoints is avoided. Moreover, the dimensions of the windows can be easily calibrated to induce a fly-by or a fly-through the waypoint.This study can be classified into the category of continuous descent operations (CDO) [7]. During CDO, aircraft descend from the cruise altitude to the final approach fix at or near idle thrust without level segments at low altitude minimizing the need for high thrust levels to remain at a constant altitude and reducing the environmental impact. Actually, the term CDO makes reference to the different techniques to maximize operational efficiency and, at the same time, fulfilling local airspace requirements and constraints. These operations are known as continuous descent arrivals, optimized profile descents (OPDs), tailored arrivals, 3D path arrival management, and continuous descent approaches (CDA). In particular, an OPD is a descent profile normally associated with a standard terminal arrival route (STAR) and designed to allow maximum use of a CDO. Planning CDOs is one of the functions of the so-called arrival managers (AMANs) whose purpose is to ensure an optimal sequencing and spacing of arrival traffic [8].Most of the previous research on CDOs based on optimal control theory focused on the trajectory optimization of a single aircraft. In [9], a multiphase optimal control method based on the pseudospectral technique has been employed to optimize vertical trajectories for individual aircraft in CDAs. Since the lateral path is assumed to be given by a STAR procedure, this work focused on optimizing vertical profile only using time and fuel consumption as performance indices. All the phases are formulated based on operational constraints and flap/gear schedules. The initial along track distance is free. Hence, it is possible to calculate both the optimal top-of-descent (TOD) and CDA trajectory. The optimal trajectories have been computed for two aircraft types: a Boeing 737-500 and a Boeing 767-400.In [10], the vertical trajectory optimization for the en route descent phase of an aircraft has been studied in the presence of both along-track and cross winds, which are both modeled as functions of altitude. Flight idle thrust was assumed during the entire descent phase. The problem is formulated as an optimal control problem. The flight range was specified from a point during the latter stages of the cruise to the meter fix. Calibrated air speed (CAS) and Mach constraints, which are the state path constraints, are considered, along with flight path angle constraints, and a maximum descent rate limit, which is a mixed input and state path constraint. The descent trajectory is optimized with respect to two cost functionals: fuel and emissions. The effects of wind speed, windshear, and cross-wind on the optimal trajectory have been analyzed using the models of two types of aircraft, Boeing 737-500 and Boeing 767-400.Less research efforts have been devoted to combined optimization of trajectories of multiple aircraft and sequencing for approaching a terminal manoeuvring area (TMA) in which all aircraft follow CDAs, while satisfying the operational requirements. This fact motivated the study presented in this paper.Two types of CDA exist depending on the lateral path followed, generally referred to as CDA under vectoring and advanced CDA. In the first case, the lateral path followed by the aircraft is assumed to be specified through instructions provided by the air traffic control (ATC). In the second case, the lateral path of the aircraft is predefined and is based on the STAR.Each of the problems that has been solved to validate the method proposed in this paper can be stated as follows. Given the dynamic models of a set of aircraft, their initial and final states, a set of operational constraints, find the optimal trajectories that steer the aircraft from the initial to the final states, fulfilling all the constraints and optimizing an objective functional.In particular, the optimal trajectories of multiple aircraft in converging arrival routes are computed taking into account time separation constraints, distance separation constraints, and their optimized profile descent along a STAR lateral profile. The problem has been solved using optimal control techniques. In particular, the OCP is transcribed using a Hermite-Simpson collocation method [11]. The resulting nonlinear programming (NLP) problem has been solved using the NLP solver IPOPT [12].The paper is structured as follows. In Section2, the general optimal control problem for multiple dynamical systems is stated and the direct collocation approach for its resolution is described. In Section 3, the aircraft equation of motion and the flight envelope constraints are stated. In Section 4, the general approach to model logical constraints is presented, which is then particularized to model time-based and distance-based separation constraints between aircraft, obstacle avoidance constraints, and waypoint constraints. In Section 5, the results of the application of the proposed method to solve several trajectory optimization problems for multiple aircraft with logical constraints in disjunctive form are reported and discussed. Finally, in Section 6, some conclusions are drawn.


## 2. Optimal Control Approach

### 2.1. Statement of the Optimal Control Problem

The multiaircraft flight planning problem considered in this paper can be regarded as a multitrajectory optimization problem in which the motion of each aircraft has been modeled as a differential algebraic dynamic system(1)∑p=fp:Xp×Up×ℝnsp×tIp,tFp⟶ℝnxp,gp:Xp×Up×ℝnsp×tIp,tFp⟶Rnzp,for p=1,2,⋯,Np, where fp describes the right-hand side of the differential equation
(2)ẋpt=fpxpt,upt,sp,t,and gp describes the algebraic constraints
(3)0=gpxpt,upt,sp,t,where Xp⊆ℝnxp and Up⊆ℝnup are the state and control sets, respectively, xpt∈ℝnxp is a nxp-dimensional state variable, upt∈ℝnup is a nup-dimensional control input, sp∈ℝnsp is a vector of parameters, and t∈tIp,tFp represents time, in which tIp and tFp denote the initial time and final time for aircraft p, with p=1,2,⋯,Np.Since this multiaircraft flight planning problem also involves operative performances and flight envelope conditions for multiple aircraft, as well as the optimization of a specified performance index, the multitrajectory optimization problem can be formulated as an OCP of a set of dynamic systems in which the goal is to find the trajectories and the corresponding control inputs that steer the states of the systems between two configurations, satisfying a set of constraints on the state and/or control variables while minimizing an objective functional.Therefore, the optimal control problem considered in this work can be stated as follows:(4a)minJxt,ut,s,t=∑p=1NpΦtFp,xptFp+∑p=1Np∫tIptFpLpxpt,upt,sp,tdt(4b)subjecttoẋt=fxt,ut,s,t(4c)0=gxt,ut,s,t(4d)ϕl≤ϕxt,ut,s,t≤ϕu(4e)xtI=xI(4f)ψxtF=0,where
(5)xt=x1t,x2t,⋯,xNptT,ut=u1t,u2t,⋯,uNptT,s=s1,s2,⋯,sNpT.The objective function(6)J:ℝnx×ℝnu×ℝns×tI,tF⟶ℝis given in Bolza form. It is expressed as a combination of a Mayer term
(7)∑p=1NpΦtFp,xptFpand a Lagrange term
(8)∑p=1Np∫tIptFpLpxpt,upt,sp,tdt.Functions(9)Φp:tIp,tFp×ℝnxp⟶ℝ,Lp:ℝnxp×ℝnup×ℝnzp×tIp,tFp⟶ℝare assumed to be twice differentiable. Function f is assumed to be piecewise Lipschitz continuous within the time interval tI,tF, and the derivative of the algebraic right-hand side function g with respect to z, that is,
(10)∂g∂z∈ℝnz×nzis assumed to be regular within the time interval tI,tF, where zt∈ℝnz denotes the nz-dimensional algebraic variable, that is, the state variables without time derivative. Vector xI∈ℝnx represents the initial conditions given at the initial time tI and function
(11)ψ:ℝnx⟶ℝnψprovides the terminal conditions at the final time tF, and it is assumed to be twice differentiable. The system must also satisfy algebraic path constraints within the time interval tI,tF given by the vector function
(12)ϕ:ℝnx×ℝnu×ℝns×tI,tF⟶ℝnϕwith lower bound ϕl∈ℝnϕ and upper bound ϕu∈ℝnϕ. Function ϕ is assumed to be twice differentiable.In the objective function (4), the Lagrange term represents a running cost, whereas the Mayer terms represent a terminal cost. A usual Lagrange objective function is to minimize the total amount of fuel consumed during the flight. A typical Mayer objective function is to minimize the duration of the flight. Equations (4b) and (4c) represent the differential-algebraic equation system that governs the motion of the dynamical system, e.g., the aircraft. Equation (4d) models the physical performance limitations of the dynamical system, typically expressed as upper and lower bounds on both states and control variables. Equations (4e) and (4f) denote the boundary (initial and final, respectively) conditions of the aircraft. Note that Equations (4c) and (4d) will also include the logical constraints that model conflict detection and resolution, and operational constraints as described in Section 4, which are of special interest for the problem studied in this paper.Hence, the optimal control problem (4a), (4b), (4c), (4d), (4e) and (4f) consists in finding an admissible control u∗t such that the set of aircraft follow an admissible trajectory x∗t between the initial and final state that minimizes the performance index Jxt,ut,s,t. The final time, tF, may be fixed or free.


### 2.2. Direct Collocation Transcription of the Optimal Control Problem

A direct numerical method has been employed to transcribe the OCP into a NLP problem. More specifically, a Hermite-Simpson direct collocation method [11] has been used. The time interval tI,tF has been subdivided into N subintervals of equal length, whose endpoints are
(13)t0,t1,⋯,tN,with t0=tI and tN=tF. In each subinterval ti,ti+1, i=0,⋯,N−1, the Hermite-Simpson numerical integration scheme has been used.The set of constraints of the resulting NLP problem includes the Hermite-Simpson system constraints that correspond to the differential constraint (4b) and the discretized versions of the other constraints of the optimal control problem. They include the algebraic constraints (4c), the state and control envelope constraints (4d), and the boundary conditions (4e) and (4f). The unknowns of the NLP problem are the values of the state and the control variables at the endpoints of each subinterval ti,ti+1, i=0,⋯,N−1.To solve the resulting NLP problem, the open source IPOPT solver [12] has been employed. It implements an interior point line search filter method and it is able to handle properly large-scale sparse nonconvex NLP problems, with a large number of equality and inequality constraints. Source and binary files are available at the Computational Infrastructure for Operations Research (COIN-OR) website (https://www.coin-or.org/).


## 2.1. Statement of the Optimal Control Problem

The multiaircraft flight planning problem considered in this paper can be regarded as a multitrajectory optimization problem in which the motion of each aircraft has been modeled as a differential algebraic dynamic system(1)∑p=fp:Xp×Up×ℝnsp×tIp,tFp⟶ℝnxp,gp:Xp×Up×ℝnsp×tIp,tFp⟶Rnzp,for p=1,2,⋯,Np, where fp describes the right-hand side of the differential equation
(2)ẋpt=fpxpt,upt,sp,t,and gp describes the algebraic constraints
(3)0=gpxpt,upt,sp,t,where Xp⊆ℝnxp and Up⊆ℝnup are the state and control sets, respectively, xpt∈ℝnxp is a nxp-dimensional state variable, upt∈ℝnup is a nup-dimensional control input, sp∈ℝnsp is a vector of parameters, and t∈tIp,tFp represents time, in which tIp and tFp denote the initial time and final time for aircraft p, with p=1,2,⋯,Np.Since this multiaircraft flight planning problem also involves operative performances and flight envelope conditions for multiple aircraft, as well as the optimization of a specified performance index, the multitrajectory optimization problem can be formulated as an OCP of a set of dynamic systems in which the goal is to find the trajectories and the corresponding control inputs that steer the states of the systems between two configurations, satisfying a set of constraints on the state and/or control variables while minimizing an objective functional.Therefore, the optimal control problem considered in this work can be stated as follows:(4a)minJxt,ut,s,t=∑p=1NpΦtFp,xptFp+∑p=1Np∫tIptFpLpxpt,upt,sp,tdt(4b)subjecttoẋt=fxt,ut,s,t(4c)0=gxt,ut,s,t(4d)ϕl≤ϕxt,ut,s,t≤ϕu(4e)xtI=xI(4f)ψxtF=0,where
(5)xt=x1t,x2t,⋯,xNptT,ut=u1t,u2t,⋯,uNptT,s=s1,s2,⋯,sNpT.The objective function(6)J:ℝnx×ℝnu×ℝns×tI,tF⟶ℝis given in Bolza form. It is expressed as a combination of a Mayer term
(7)∑p=1NpΦtFp,xptFpand a Lagrange term
(8)∑p=1Np∫tIptFpLpxpt,upt,sp,tdt.Functions(9)Φp:tIp,tFp×ℝnxp⟶ℝ,Lp:ℝnxp×ℝnup×ℝnzp×tIp,tFp⟶ℝare assumed to be twice differentiable. Function f is assumed to be piecewise Lipschitz continuous within the time interval tI,tF, and the derivative of the algebraic right-hand side function g with respect to z, that is,
(10)∂g∂z∈ℝnz×nzis assumed to be regular within the time interval tI,tF, where zt∈ℝnz denotes the nz-dimensional algebraic variable, that is, the state variables without time derivative. Vector xI∈ℝnx represents the initial conditions given at the initial time tI and function
(11)ψ:ℝnx⟶ℝnψprovides the terminal conditions at the final time tF, and it is assumed to be twice differentiable. The system must also satisfy algebraic path constraints within the time interval tI,tF given by the vector function
(12)ϕ:ℝnx×ℝnu×ℝns×tI,tF⟶ℝnϕwith lower bound ϕl∈ℝnϕ and upper bound ϕu∈ℝnϕ. Function ϕ is assumed to be twice differentiable.In the objective function (4), the Lagrange term represents a running cost, whereas the Mayer terms represent a terminal cost. A usual Lagrange objective function is to minimize the total amount of fuel consumed during the flight. A typical Mayer objective function is to minimize the duration of the flight. Equations (4b) and (4c) represent the differential-algebraic equation system that governs the motion of the dynamical system, e.g., the aircraft. Equation (4d) models the physical performance limitations of the dynamical system, typically expressed as upper and lower bounds on both states and control variables. Equations (4e) and (4f) denote the boundary (initial and final, respectively) conditions of the aircraft. Note that Equations (4c) and (4d) will also include the logical constraints that model conflict detection and resolution, and operational constraints as described in Section 4, which are of special interest for the problem studied in this paper.Hence, the optimal control problem (4a), (4b), (4c), (4d), (4e) and (4f) consists in finding an admissible control u∗t such that the set of aircraft follow an admissible trajectory x∗t between the initial and final state that minimizes the performance index Jxt,ut,s,t. The final time, tF, may be fixed or free.


## 2.2. Direct Collocation Transcription of the Optimal Control Problem

A direct numerical method has been employed to transcribe the OCP into a NLP problem. More specifically, a Hermite-Simpson direct collocation method [11] has been used. The time interval tI,tF has been subdivided into N subintervals of equal length, whose endpoints are
(13)t0,t1,⋯,tN,with t0=tI and tN=tF. In each subinterval ti,ti+1, i=0,⋯,N−1, the Hermite-Simpson numerical integration scheme has been used.The set of constraints of the resulting NLP problem includes the Hermite-Simpson system constraints that correspond to the differential constraint (4b) and the discretized versions of the other constraints of the optimal control problem. They include the algebraic constraints (4c), the state and control envelope constraints (4d), and the boundary conditions (4e) and (4f). The unknowns of the NLP problem are the values of the state and the control variables at the endpoints of each subinterval ti,ti+1, i=0,⋯,N−1.To solve the resulting NLP problem, the open source IPOPT solver [12] has been employed. It implements an interior point line search filter method and it is able to handle properly large-scale sparse nonconvex NLP problems, with a large number of equality and inequality constraints. Source and binary files are available at the Computational Infrastructure for Operations Research (COIN-OR) website (https://www.coin-or.org/).


## 3. Aircraft Model Description

Following [13], a common three-degree-of-freedom dynamic model has been used which describes the point variable-mass motion of the aircraft over a spherical Earth model. In particular, a symmetric flight has been considered. Thus, it has been assumed that there is no sideslip and all forces lie in the plane of symmetry of aircraft.


### 3.1. Equations of Motion

The following equations of motion of the aircraft have been considered:(14)V̇t=Tt−Dhet,Vt,CLt−mt·g·sinγtmt,χ̇t=Lhet,Vt,CLt·sinμtmt·Vt·cosγt,γ̇t=Lhet,Vt,CLt·cosμt−mt·g·cosγtmt·Vt,λėt=Vt·cosγt·cosχtR·cosθet,θėt=Vt·cosγt·sinχtR,hėt=Vt·sinγt,ṁt=−Tt·ηVt.The three dynamic equations in (14) are expressed in an aircraft-attached reference frame xw,yw,zw, and the three kinematic equations are expressed in a ground-based reference frame xe,ye,ze as shown in Figure 1.Figure 1
Aircraft state and forces.
(a)
Top view
(b)
Front view
(c)
Lateral viewThe states of the system (14) are V,χ,γ,λe,θe,he, and m. Thus, xt=Vt,χt,γt,λet,θet,het,mt. The state variables V,χ, and γ are the true airspeed, the heading angle, and the flight path angle, respectively. The state variables λe,θe, and he are the aircraft 3D position, the longitude, the latitude, and the altitude, respectively. The components of the aircraft position vector in two dimensions (xe, ye) are approximated as xe=λe·R+he·cosθe and ye=θe·R+he. Finally, the state variable m is the aircraft mass. The control inputs are the bank angle μ, the engine thrust T, and the lift coefficient CL. Thus, ut=Tt,μt,CLt.The parameterR is the radius of Earth and η is the speed-dependent fuel efficiency coefficient. Lift, L=CLSq^, and drag, D=CDSq^, are the components of the aerodynamic force. The parameter S is the reference wing surface area and q^=1/2ρV2 is the dynamic pressure. A parabolic drag polar CD=CD0+KCL2 and an International Standard Atmosphere (ISA) model are assumed. The lift coefficient CL is a known function of the angle of attack α and the Mach number.Note that differential equations in (14) take the form of (4b) of the continuous optimal control problem stated in Section 2.1.


### 3.2. Flight Envelope Constraints

Flight envelope constraints are derived from the geometry of the aircraft, structural limitations, engine power, and aerodynamic characteristics. The performance limitation model and the parameters have been obtained from the Base of Aircraft Data (BADA), version 3.6 [14]:
(15)0≤het≤minhM0,hut,γmin≤γt≤γmax,Mt≤MM0,mmin≤mt≤mmax,V̇t≤a¯l,CvVst≤Vt≤VMo,γ̇tVt≤a¯n,0.1≤CLt≤CLmax,Tmint≤Tt≤Tmaxt,μt≤μ̄.In (15), hM0 is the maximum reachable altitude and hut is the maximum operative altitude at a given mass (it increases as fuel is burned). Mt is the Mach number and MM0 is the maximum operating Mach number. Cv is the minimum speed coefficient, Vst is the stall speed, VM0 is the maximum operating CAS, and ān and āl are, respectively, the maximum normal and longitudinal accelerations for civilian aircraft. Finally, Tmint and Tmaxt correspond to the minimum and maximum available thrust, respectively, and μ̄ corresponds to the maximum bank angle due to structural limitations.Note that inequality constraints in (15) take the form of (4d) of the continuous optimal control problem stated in Section 2.1.


## 3.1. Equations of Motion

The following equations of motion of the aircraft have been considered:(14)V̇t=Tt−Dhet,Vt,CLt−mt·g·sinγtmt,χ̇t=Lhet,Vt,CLt·sinμtmt·Vt·cosγt,γ̇t=Lhet,Vt,CLt·cosμt−mt·g·cosγtmt·Vt,λėt=Vt·cosγt·cosχtR·cosθet,θėt=Vt·cosγt·sinχtR,hėt=Vt·sinγt,ṁt=−Tt·ηVt.The three dynamic equations in (14) are expressed in an aircraft-attached reference frame xw,yw,zw, and the three kinematic equations are expressed in a ground-based reference frame xe,ye,ze as shown in Figure 1.Figure 1
Aircraft state and forces.
(a)
Top view
(b)
Front view
(c)
Lateral viewThe states of the system (14) are V,χ,γ,λe,θe,he, and m. Thus, xt=Vt,χt,γt,λet,θet,het,mt. The state variables V,χ, and γ are the true airspeed, the heading angle, and the flight path angle, respectively. The state variables λe,θe, and he are the aircraft 3D position, the longitude, the latitude, and the altitude, respectively. The components of the aircraft position vector in two dimensions (xe, ye) are approximated as xe=λe·R+he·cosθe and ye=θe·R+he. Finally, the state variable m is the aircraft mass. The control inputs are the bank angle μ, the engine thrust T, and the lift coefficient CL. Thus, ut=Tt,μt,CLt.The parameterR is the radius of Earth and η is the speed-dependent fuel efficiency coefficient. Lift, L=CLSq^, and drag, D=CDSq^, are the components of the aerodynamic force. The parameter S is the reference wing surface area and q^=1/2ρV2 is the dynamic pressure. A parabolic drag polar CD=CD0+KCL2 and an International Standard Atmosphere (ISA) model are assumed. The lift coefficient CL is a known function of the angle of attack α and the Mach number.Note that differential equations in (14) take the form of (4b) of the continuous optimal control problem stated in Section 2.1.


## 3.2. Flight Envelope Constraints

Flight envelope constraints are derived from the geometry of the aircraft, structural limitations, engine power, and aerodynamic characteristics. The performance limitation model and the parameters have been obtained from the Base of Aircraft Data (BADA), version 3.6 [14]:
(15)0≤het≤minhM0,hut,γmin≤γt≤γmax,Mt≤MM0,mmin≤mt≤mmax,V̇t≤a¯l,CvVst≤Vt≤VMo,γ̇tVt≤a¯n,0.1≤CLt≤CLmax,Tmint≤Tt≤Tmaxt,μt≤μ̄.In (15), hM0 is the maximum reachable altitude and hut is the maximum operative altitude at a given mass (it increases as fuel is burned). Mt is the Mach number and MM0 is the maximum operating Mach number. Cv is the minimum speed coefficient, Vst is the stall speed, VM0 is the maximum operating CAS, and ān and āl are, respectively, the maximum normal and longitudinal accelerations for civilian aircraft. Finally, Tmint and Tmaxt correspond to the minimum and maximum available thrust, respectively, and μ̄ corresponds to the maximum bank angle due to structural limitations.Note that inequality constraints in (15) take the form of (4d) of the continuous optimal control problem stated in Section 2.1.


## 4. Logical Constraints Modeling

In this section, the approach proposed in [5] has been followed in which an extension of the embedding optimal control technique stated in [15] and developed in [16] was proposed. The embedding technique introduced in [15, 16], to transform hybrid optimal control problems into traditional smooth optimal control problems, in which the discrete aspect of the system arisen only from switches in the dynamic equations, was adapted in [5] to deal with the logical (discrete) components which also might appear as constraints.It was shown in [17] that every Boolean expression can be transformed into conjunctive normal form (CNF). Thus, it has been assumed that any logical constraint considered in this study can be written as a CNF expression
(16)Q1∧Q2∧⋯∧Qn,where
(17)Qi=Pi1∨Pi2∨⋯∨Pimiis in disjunctive form, in which only one of the propositions must be satisfied and the proposition Pij is either Xij or ¬Xij. Term Xij is a literal that can be either True or False and ¬ represents the negation or logical complement operator. The term Xij is used to represent statements such as “longitude λe≤40”. Therefore, Pij takes the form
(18)gijxt≤0,where gij:ℝnx⟶ℝ is assumed to be a C1 function.In order to include the logical constraint (16) in a smooth continuous optimal control problem formulation, it must be converted into a set of equality or inequality constraints in which binary variables are not employed. In this way, the combinatorial complexity of integer programming is eluded.Transcribing the conjunction in (16) is straightforward since it is equivalent to the following expression:
(19)∀i∈1,2,⋯,n:Qi.Thus, taking into account (17), the logical expression (16) can be represented as follows:
(20)∀i∈1,2,⋯,n:Pi1∨Pi2∨⋯∨Pimi.To transcribe these disjunctions into a set of inequality constraints, a continuous variableβij∈0,1 is defined and related to each Pij in (18). Thus, (20) can be expressed as follows:
(21)∀i∈1,2,⋯,n:βij·gijxt≤0,0≤βij≤1,∑j=1miβij=1.It is immediate to check that ifβij=0 in the first term in (21), then constraint gijxt≤0 is not fulfilled. On the other hand, if 0<βij≤1, then βij·gijxt≤0 is in fact gijxt≤0, and thus, constraint gijxt≤0 is enforced. Finally, the last term in (21) guarantees that at least one of the propositions Pij holds.Note that, as expected, equality and inequality constraints in (21) take the form of (4c) and (4d), respectively, of the traditional continuous optimal control problem stated in Section 2.1. In the following subsections, the application of this technique to several instances of interest in the ATM context will be presented in detail.


### 4.1. Collision Avoidance Constraints

#### 4.1.1. Time-Based Separation between Aircraft

One of the problems of most interest in ATM is the conflict detection and resolution problem [18]. Two different instances of this problem will be studied in this subsection and in the next one. In a first instance, a collision avoidance model among different aircraft along routes converging at the same waypoint has been considered, in which a safety time-based separation is guaranteed at the merging waypoint.LettpF and tqF be the unfixed time at the merging waypoint of aircraft p and q, respectively, and let dt be the safety time difference between two consecutive aircraft. Since in terms of the discretization (13), tpF=tpN and tqF=tqN, then multiaircraft time-based separation constraints can be expressed as follows:
(22)∀q>p:∣tpN−tqN∣≥dt,where condition q>p prevents unnecessary duplication of constraints. Constraints (22) can be rewritten as follows:
(23)∀q>p:tpN−tqN≥dtortqN−tpN≥dt.Following the technique described above to tackle constraints in disjunctive form, if we define new variablesδpN,qN1,δpN,qN2∈0,1 satisfying condition
(24)δpN,qN1+δpN,qN2=1,Equation (23) can be transformed into
(25)∀q>p:δpN,qN1tpN−tqN−dt≥0,δpN,qN2tqN−tpN−dt≥0,0≤δpN,qNj≤1,j=1,2,δpN,qN1+δpN,qN2=1.The last constraint in (25) ensures that at least one of the constraints in (23) is fulfilled, that is, the safety time-based separation between aircraft p and q is guaranteed.


#### 4.1.2. Distance-Based Separation between Aircraft

In the second instance, to model collision avoidance among different aircraft, at each endpoint of the subintervals of the discretization (13) and for every couple of aircraft, a safety distance-based separation is guaranteed in the horizontal or vertical directions, that is, λ−θ or h directions.Letλpi,θpi,hpi and λqi,θqi,hqi be the positions of aircraft p and q at the endpoint ti of the discretization, respectively. If the safety distance-based separation in the λ−θ and h directions is denoted by dλθ and dh, respectively, then the distance-based separation constraints can be expressed as follows:
(26)∀q>p,∀i∈1,2,⋯,N:2Ratan2ζi,1−ζi≥dλθor∣hpi−hqi∣≥dh,where the haversine formula has been considered for the distance in the λ−θ direction with
(27)ζi=sin2θqi−θpi2+cosθpicosθqisin2λqi−λpi2.As before, conditionq>p prevents unnecessary duplication of constraints, and the set of constraints (26) can be rewritten as follows:
(28)∀q>p,∀i∈1,2,⋯,N:2Ratan2ζi,1−ζi≥dλθorhpi−hqi≥dhorhqi−hpi≥dh.Again, following the technique described above, new variablesνpi,qij∈0,1 for i=1,2,⋯,N,j=1,2,3, satisfying condition
(29)∑j=13νpi,qij=1,are introduced, and Equation (28) can be transformed into
(30)∀q>p,∀i∈1,2,⋯,N:νpi,qi12Ratan2ζi,1−ζi−dλθ≥0,νpi,qi2hpi−hqi−dh≥0,νpi,qi3hqi−hpi−dh≥0,0≤νpi,qij≤1,j=1,2,3,∑j=13νpi,qij=1.The last constraint in (30) ensures that at least one of the constraints in (28) is fulfilled, that is, aircraft p and q are guaranteed to be a safe distance apart. Moreover, since this embedded logical constraint approach is quite general, any other distance separation model described in terms of Equation (17) can be considered.


#### 4.1.3. Obstacle Avoidance

In the third instance, an obstacle avoidance problem has been considered. Following [5], in which a set of stationary rectangular obstacles in a two-dimensional environment were assumed, stationary obstacles on the three-dimensional space have been enveloped by cuboids. Cuboids give computational advantages with respect to other models, such as ellipsoids, since they give rise to simple bound constraints instead of quadratic constraints, and mathematical programming solvers are able to deal with bound constraints more efficiently than with nonlinear constraints. Additionally, this simple approach allows us to specify each cuboid by giving only the coordinates of two opposite corners as shown in Figure 2.Figure 2
Two opposite corners of a cuboid in the 3D space.Letλl,θl,hl and λu,θu,hu be the positions of opposite corners of the cuboid. Cuboid avoidance for a single aircraft involves that at every endpoint ti of the subintervals of the discretization(13), the position of the aircraft λi,θi,hi must remain outside it. In terms of logical constraints, this condition can be expressed as follows:
(31)∀i∈1,2,⋯,N−1:λi−λl≤0orλu−λi≤0orθi−θl≤0orθu−θi≤0orhi−hl≤0orhu−hi≤0.Note that constraints (31) are not enforced at the initial and final points of the discretization, t0 and tN, since in practice both positions are fixed. Moreover, due to the nature of the collocation technique described in Section 2.2, constraints are not enforced between the endpoints of the discretized differential system of equations. Therefore, small intrusions inside the cuboids might occur. A smaller discretization size implies a larger intrusion, and vice versa. In practice, this drawback can be easily overcome by considering a sufficient safety margin in the modeling of the cuboids.Following the same technique used in the previous sections, new variablesξij∈0,1, for i=1,2,⋯,N−1,j=1,2,⋯,6, satisfying condition
(32)∑j=16ξij=1,are introduced, and Equation (31) can be transformed into
(33)∀i∈1,2,⋯,N−1:ξi1λi−λl≤0,ξi2λu−λi≤0,ξi3θi−θl≤0,ξi4θu−θi≤0,ξi5hi−hl≤0,ξi6hu−hi≤0,0≤ξij≤1,j=1,2,⋯,6,∑j=16ξij=1.Again, as mentioned above, the last constraint in (33) ensures that at least one of the constraints in (31) is fulfilled, which means that aircraft flies outside the cuboid. Also note that the set of constraints in (33) can be straightforward extended to multiobstacle and multiaircraft settings. As already pointed out in Section 4.1.2, since this embedded logical constraint approach is quite general, any other obstacle model described in terms of Equation (17) can be considered.


### 4.2. Waypoint Constraints

Finally, the modeling of an aircraft flying through a waypoint has been considered. The model has been based on the use of cuboids similar to those considered in Section4.1.3. In this case, a cuboid centered at each waypoint has been defined. However, unlike the obstacle avoidance constraints, in this setting, the aircraft has forced to pass through the cuboid instead of avoiding it.LetλWl,θWl,hWl and λWu,θWu,hWu be the positions of opposite corners of a cuboid surrounding a single waypoint. Flying by this waypoint (that is, passing through the corresponding cuboid) involves that at every endpoint ti of the subintervals of the discretization (13), the position of the aircraft λi,θi,hi must remain inside it. In terms of logical constraints, this condition can be expressed as follows:
(34)∀i∈1,2,⋯,N−1:λWl−λi≤0,λi−λWu≤0,θWl−θi≤0,θi−θWu≤0,hWl−hi≤0,hi−hWu≤0.Note that, for one hand, constraints (34) are enforced at every point of the discretization except for the initial and final points, t0 and tN, since there is no a priori knowledge about when the aircraft is going to fly by the waypoint. On the other hand, these constraints obviously make sense only when the aircraft is closed enough to the waypoint.To overcome this drawback, a second auxiliary cuboid is considered to modelled free flight mode of the aircraft. Letλmin,θmin,hmin,λmax,θmax, and hmax be the minimum and maximum values of the state variables λ,θ, and h, respectively. In terms of logical constraints, the free flight mode condition can be expressed as follows:
(35)∀i∈1,2,⋯,N−1:λmin−λi≤0,λi−λmax≤0,θmin−θi≤0,θi−θmax≤0,hmin−hi≤0,hi−hmax≤0.Then, the transcription into a logical disjunction which allows to select along the whole trajectory between flying by the waypoint mode (WM) or free flight mode (FM), namely,WM∨FM, can be expressed as follows:
(36)∀i∈1,2,⋯,N−1:WMiorFMi,where WMi and FMi denote if at discretization instant i, the aircraft is in waypoint mode or free flight mode, respectively.Once again, following the technique described above, if we define new variablesκi1,κi2∈0,1 satisfying condition
(37)κi1+κi2=1,foralli=1,2,⋯,N−1,Equation (36) can be transformed into
(38)∀i∈1,2,⋯,N−1:κi1λWl−λi≤0,κi2λmin−λi≤0,κi1λi−λWu≤0,κi2λi−λmax≤0,κi1θWl−θi≤0,κi2θmin−θi≤0,κi1θi−θWu≤0,κi2θi−θmax≤0,κi1hWl−hi≤0,κi2hmin−hi≤0,κi1hi−hWu≤0,κi2hi−hmax≤0,0≤κij≤1,j=1,2,∑j=12κij=1.The last constraint in (38) ensures that at least one of the conditions in (36) is fulfilled, which means that at each discretization instant i, the aircraft flies in FM or WM. Note that, depending on the performance index (4a) considered in the optimal control problem, two related undesired issues could potentially arise.On the one hand, the optimal solution could provide a trajectory in which the aircraft flies in FM for each discretization instanti. Therefore, in order to force the aircraft to actually fly by the waypoint, a penalty term must also be added to the numerical transcription of the performance index (4a). The use of this penalty term, which is set forth to encode the desired control objectives, implies that the numerical solution of the optimal control problem must combine the usual collocation technique describe in Section 2.2 with a penalty function methodology. In particular, in this work, the well-known continuation method has been implemented following a similar approach to [19]. In this specific case of waypoint modeling, the penalty term takes the form
(39)c1∑i=1N−1κi1+c2∑i=1N−1κi2,where c1 and c2 are suitable constants determined by the continuation method. In the context of a minimization performance index such as (4a), a large enough value c2 such that c2≫c1>0, which penalizes the free flight mode, guarantees that the aircraft actually flies by the waypoint.On the other hand, the optimal solution could provide a trajectory in which the aircraft flies by the waypoint more than once. This situation can be easily avoid introducing in the model the following simple constraint:(40)∑i=1N−1κi1≤c3,where c3∈1,2,⋯,N−1 is a suitable constant. In the context of the problem considered in this work including the penalty term (39), a value of c3 between 2 and 4 is enough to avoid this potential undesired issue.Note that this waypoint constraint setting can be straightforward extended to multiwaypoint and multiaircraft settings. Moreover, once again, since this embedded logical constraint approach is quite general, any other waypoint constraint described in terms of Equation (17) can be considered.


## 4.1. Collision Avoidance Constraints

### 4.1.1. Time-Based Separation between Aircraft

One of the problems of most interest in ATM is the conflict detection and resolution problem [18]. Two different instances of this problem will be studied in this subsection and in the next one. In a first instance, a collision avoidance model among different aircraft along routes converging at the same waypoint has been considered, in which a safety time-based separation is guaranteed at the merging waypoint.LettpF and tqF be the unfixed time at the merging waypoint of aircraft p and q, respectively, and let dt be the safety time difference between two consecutive aircraft. Since in terms of the discretization (13), tpF=tpN and tqF=tqN, then multiaircraft time-based separation constraints can be expressed as follows:
(22)∀q>p:∣tpN−tqN∣≥dt,where condition q>p prevents unnecessary duplication of constraints. Constraints (22) can be rewritten as follows:
(23)∀q>p:tpN−tqN≥dtortqN−tpN≥dt.Following the technique described above to tackle constraints in disjunctive form, if we define new variablesδpN,qN1,δpN,qN2∈0,1 satisfying condition
(24)δpN,qN1+δpN,qN2=1,Equation (23) can be transformed into
(25)∀q>p:δpN,qN1tpN−tqN−dt≥0,δpN,qN2tqN−tpN−dt≥0,0≤δpN,qNj≤1,j=1,2,δpN,qN1+δpN,qN2=1.The last constraint in (25) ensures that at least one of the constraints in (23) is fulfilled, that is, the safety time-based separation between aircraft p and q is guaranteed.


### 4.1.2. Distance-Based Separation between Aircraft

In the second instance, to model collision avoidance among different aircraft, at each endpoint of the subintervals of the discretization (13) and for every couple of aircraft, a safety distance-based separation is guaranteed in the horizontal or vertical directions, that is, λ−θ or h directions.Letλpi,θpi,hpi and λqi,θqi,hqi be the positions of aircraft p and q at the endpoint ti of the discretization, respectively. If the safety distance-based separation in the λ−θ and h directions is denoted by dλθ and dh, respectively, then the distance-based separation constraints can be expressed as follows:
(26)∀q>p,∀i∈1,2,⋯,N:2Ratan2ζi,1−ζi≥dλθor∣hpi−hqi∣≥dh,where the haversine formula has been considered for the distance in the λ−θ direction with
(27)ζi=sin2θqi−θpi2+cosθpicosθqisin2λqi−λpi2.As before, conditionq>p prevents unnecessary duplication of constraints, and the set of constraints (26) can be rewritten as follows:
(28)∀q>p,∀i∈1,2,⋯,N:2Ratan2ζi,1−ζi≥dλθorhpi−hqi≥dhorhqi−hpi≥dh.Again, following the technique described above, new variablesνpi,qij∈0,1 for i=1,2,⋯,N,j=1,2,3, satisfying condition
(29)∑j=13νpi,qij=1,are introduced, and Equation (28) can be transformed into
(30)∀q>p,∀i∈1,2,⋯,N:νpi,qi12Ratan2ζi,1−ζi−dλθ≥0,νpi,qi2hpi−hqi−dh≥0,νpi,qi3hqi−hpi−dh≥0,0≤νpi,qij≤1,j=1,2,3,∑j=13νpi,qij=1.The last constraint in (30) ensures that at least one of the constraints in (28) is fulfilled, that is, aircraft p and q are guaranteed to be a safe distance apart. Moreover, since this embedded logical constraint approach is quite general, any other distance separation model described in terms of Equation (17) can be considered.


### 4.1.3. Obstacle Avoidance

In the third instance, an obstacle avoidance problem has been considered. Following [5], in which a set of stationary rectangular obstacles in a two-dimensional environment were assumed, stationary obstacles on the three-dimensional space have been enveloped by cuboids. Cuboids give computational advantages with respect to other models, such as ellipsoids, since they give rise to simple bound constraints instead of quadratic constraints, and mathematical programming solvers are able to deal with bound constraints more efficiently than with nonlinear constraints. Additionally, this simple approach allows us to specify each cuboid by giving only the coordinates of two opposite corners as shown in Figure 2.Figure 2
Two opposite corners of a cuboid in the 3D space.Letλl,θl,hl and λu,θu,hu be the positions of opposite corners of the cuboid. Cuboid avoidance for a single aircraft involves that at every endpoint ti of the subintervals of the discretization(13), the position of the aircraft λi,θi,hi must remain outside it. In terms of logical constraints, this condition can be expressed as follows:
(31)∀i∈1,2,⋯,N−1:λi−λl≤0orλu−λi≤0orθi−θl≤0orθu−θi≤0orhi−hl≤0orhu−hi≤0.Note that constraints (31) are not enforced at the initial and final points of the discretization, t0 and tN, since in practice both positions are fixed. Moreover, due to the nature of the collocation technique described in Section 2.2, constraints are not enforced between the endpoints of the discretized differential system of equations. Therefore, small intrusions inside the cuboids might occur. A smaller discretization size implies a larger intrusion, and vice versa. In practice, this drawback can be easily overcome by considering a sufficient safety margin in the modeling of the cuboids.Following the same technique used in the previous sections, new variablesξij∈0,1, for i=1,2,⋯,N−1,j=1,2,⋯,6, satisfying condition
(32)∑j=16ξij=1,are introduced, and Equation (31) can be transformed into
(33)∀i∈1,2,⋯,N−1:ξi1λi−λl≤0,ξi2λu−λi≤0,ξi3θi−θl≤0,ξi4θu−θi≤0,ξi5hi−hl≤0,ξi6hu−hi≤0,0≤ξij≤1,j=1,2,⋯,6,∑j=16ξij=1.Again, as mentioned above, the last constraint in (33) ensures that at least one of the constraints in (31) is fulfilled, which means that aircraft flies outside the cuboid. Also note that the set of constraints in (33) can be straightforward extended to multiobstacle and multiaircraft settings. As already pointed out in Section 4.1.2, since this embedded logical constraint approach is quite general, any other obstacle model described in terms of Equation (17) can be considered.


## 4.1.1. Time-Based Separation between Aircraft

One of the problems of most interest in ATM is the conflict detection and resolution problem [18]. Two different instances of this problem will be studied in this subsection and in the next one. In a first instance, a collision avoidance model among different aircraft along routes converging at the same waypoint has been considered, in which a safety time-based separation is guaranteed at the merging waypoint.LettpF and tqF be the unfixed time at the merging waypoint of aircraft p and q, respectively, and let dt be the safety time difference between two consecutive aircraft. Since in terms of the discretization (13), tpF=tpN and tqF=tqN, then multiaircraft time-based separation constraints can be expressed as follows:
(22)∀q>p:∣tpN−tqN∣≥dt,where condition q>p prevents unnecessary duplication of constraints. Constraints (22) can be rewritten as follows:
(23)∀q>p:tpN−tqN≥dtortqN−tpN≥dt.Following the technique described above to tackle constraints in disjunctive form, if we define new variablesδpN,qN1,δpN,qN2∈0,1 satisfying condition
(24)δpN,qN1+δpN,qN2=1,Equation (23) can be transformed into
(25)∀q>p:δpN,qN1tpN−tqN−dt≥0,δpN,qN2tqN−tpN−dt≥0,0≤δpN,qNj≤1,j=1,2,δpN,qN1+δpN,qN2=1.The last constraint in (25) ensures that at least one of the constraints in (23) is fulfilled, that is, the safety time-based separation between aircraft p and q is guaranteed.


## 4.1.2. Distance-Based Separation between Aircraft

In the second instance, to model collision avoidance among different aircraft, at each endpoint of the subintervals of the discretization (13) and for every couple of aircraft, a safety distance-based separation is guaranteed in the horizontal or vertical directions, that is, λ−θ or h directions.Letλpi,θpi,hpi and λqi,θqi,hqi be the positions of aircraft p and q at the endpoint ti of the discretization, respectively. If the safety distance-based separation in the λ−θ and h directions is denoted by dλθ and dh, respectively, then the distance-based separation constraints can be expressed as follows:
(26)∀q>p,∀i∈1,2,⋯,N:2Ratan2ζi,1−ζi≥dλθor∣hpi−hqi∣≥dh,where the haversine formula has been considered for the distance in the λ−θ direction with
(27)ζi=sin2θqi−θpi2+cosθpicosθqisin2λqi−λpi2.As before, conditionq>p prevents unnecessary duplication of constraints, and the set of constraints (26) can be rewritten as follows:
(28)∀q>p,∀i∈1,2,⋯,N:2Ratan2ζi,1−ζi≥dλθorhpi−hqi≥dhorhqi−hpi≥dh.Again, following the technique described above, new variablesνpi,qij∈0,1 for i=1,2,⋯,N,j=1,2,3, satisfying condition
(29)∑j=13νpi,qij=1,are introduced, and Equation (28) can be transformed into
(30)∀q>p,∀i∈1,2,⋯,N:νpi,qi12Ratan2ζi,1−ζi−dλθ≥0,νpi,qi2hpi−hqi−dh≥0,νpi,qi3hqi−hpi−dh≥0,0≤νpi,qij≤1,j=1,2,3,∑j=13νpi,qij=1.The last constraint in (30) ensures that at least one of the constraints in (28) is fulfilled, that is, aircraft p and q are guaranteed to be a safe distance apart. Moreover, since this embedded logical constraint approach is quite general, any other distance separation model described in terms of Equation (17) can be considered.


## 4.1.3. Obstacle Avoidance

In the third instance, an obstacle avoidance problem has been considered. Following [5], in which a set of stationary rectangular obstacles in a two-dimensional environment were assumed, stationary obstacles on the three-dimensional space have been enveloped by cuboids. Cuboids give computational advantages with respect to other models, such as ellipsoids, since they give rise to simple bound constraints instead of quadratic constraints, and mathematical programming solvers are able to deal with bound constraints more efficiently than with nonlinear constraints. Additionally, this simple approach allows us to specify each cuboid by giving only the coordinates of two opposite corners as shown in Figure 2.Figure 2
Two opposite corners of a cuboid in the 3D space.Letλl,θl,hl and λu,θu,hu be the positions of opposite corners of the cuboid. Cuboid avoidance for a single aircraft involves that at every endpoint ti of the subintervals of the discretization(13), the position of the aircraft λi,θi,hi must remain outside it. In terms of logical constraints, this condition can be expressed as follows:
(31)∀i∈1,2,⋯,N−1:λi−λl≤0orλu−λi≤0orθi−θl≤0orθu−θi≤0orhi−hl≤0orhu−hi≤0.Note that constraints (31) are not enforced at the initial and final points of the discretization, t0 and tN, since in practice both positions are fixed. Moreover, due to the nature of the collocation technique described in Section 2.2, constraints are not enforced between the endpoints of the discretized differential system of equations. Therefore, small intrusions inside the cuboids might occur. A smaller discretization size implies a larger intrusion, and vice versa. In practice, this drawback can be easily overcome by considering a sufficient safety margin in the modeling of the cuboids.Following the same technique used in the previous sections, new variablesξij∈0,1, for i=1,2,⋯,N−1,j=1,2,⋯,6, satisfying condition
(32)∑j=16ξij=1,are introduced, and Equation (31) can be transformed into
(33)∀i∈1,2,⋯,N−1:ξi1λi−λl≤0,ξi2λu−λi≤0,ξi3θi−θl≤0,ξi4θu−θi≤0,ξi5hi−hl≤0,ξi6hu−hi≤0,0≤ξij≤1,j=1,2,⋯,6,∑j=16ξij=1.Again, as mentioned above, the last constraint in (33) ensures that at least one of the constraints in (31) is fulfilled, which means that aircraft flies outside the cuboid. Also note that the set of constraints in (33) can be straightforward extended to multiobstacle and multiaircraft settings. As already pointed out in Section 4.1.2, since this embedded logical constraint approach is quite general, any other obstacle model described in terms of Equation (17) can be considered.


## 4.2. Waypoint Constraints

Finally, the modeling of an aircraft flying through a waypoint has been considered. The model has been based on the use of cuboids similar to those considered in Section4.1.3. In this case, a cuboid centered at each waypoint has been defined. However, unlike the obstacle avoidance constraints, in this setting, the aircraft has forced to pass through the cuboid instead of avoiding it.LetλWl,θWl,hWl and λWu,θWu,hWu be the positions of opposite corners of a cuboid surrounding a single waypoint. Flying by this waypoint (that is, passing through the corresponding cuboid) involves that at every endpoint ti of the subintervals of the discretization (13), the position of the aircraft λi,θi,hi must remain inside it. In terms of logical constraints, this condition can be expressed as follows:
(34)∀i∈1,2,⋯,N−1:λWl−λi≤0,λi−λWu≤0,θWl−θi≤0,θi−θWu≤0,hWl−hi≤0,hi−hWu≤0.Note that, for one hand, constraints (34) are enforced at every point of the discretization except for the initial and final points, t0 and tN, since there is no a priori knowledge about when the aircraft is going to fly by the waypoint. On the other hand, these constraints obviously make sense only when the aircraft is closed enough to the waypoint.To overcome this drawback, a second auxiliary cuboid is considered to modelled free flight mode of the aircraft. Letλmin,θmin,hmin,λmax,θmax, and hmax be the minimum and maximum values of the state variables λ,θ, and h, respectively. In terms of logical constraints, the free flight mode condition can be expressed as follows:
(35)∀i∈1,2,⋯,N−1:λmin−λi≤0,λi−λmax≤0,θmin−θi≤0,θi−θmax≤0,hmin−hi≤0,hi−hmax≤0.Then, the transcription into a logical disjunction which allows to select along the whole trajectory between flying by the waypoint mode (WM) or free flight mode (FM), namely,WM∨FM, can be expressed as follows:
(36)∀i∈1,2,⋯,N−1:WMiorFMi,where WMi and FMi denote if at discretization instant i, the aircraft is in waypoint mode or free flight mode, respectively.Once again, following the technique described above, if we define new variablesκi1,κi2∈0,1 satisfying condition
(37)κi1+κi2=1,foralli=1,2,⋯,N−1,Equation (36) can be transformed into
(38)∀i∈1,2,⋯,N−1:κi1λWl−λi≤0,κi2λmin−λi≤0,κi1λi−λWu≤0,κi2λi−λmax≤0,κi1θWl−θi≤0,κi2θmin−θi≤0,κi1θi−θWu≤0,κi2θi−θmax≤0,κi1hWl−hi≤0,κi2hmin−hi≤0,κi1hi−hWu≤0,κi2hi−hmax≤0,0≤κij≤1,j=1,2,∑j=12κij=1.The last constraint in (38) ensures that at least one of the conditions in (36) is fulfilled, which means that at each discretization instant i, the aircraft flies in FM or WM. Note that, depending on the performance index (4a) considered in the optimal control problem, two related undesired issues could potentially arise.On the one hand, the optimal solution could provide a trajectory in which the aircraft flies in FM for each discretization instanti. Therefore, in order to force the aircraft to actually fly by the waypoint, a penalty term must also be added to the numerical transcription of the performance index (4a). The use of this penalty term, which is set forth to encode the desired control objectives, implies that the numerical solution of the optimal control problem must combine the usual collocation technique describe in Section 2.2 with a penalty function methodology. In particular, in this work, the well-known continuation method has been implemented following a similar approach to [19]. In this specific case of waypoint modeling, the penalty term takes the form
(39)c1∑i=1N−1κi1+c2∑i=1N−1κi2,where c1 and c2 are suitable constants determined by the continuation method. In the context of a minimization performance index such as (4a), a large enough value c2 such that c2≫c1>0, which penalizes the free flight mode, guarantees that the aircraft actually flies by the waypoint.On the other hand, the optimal solution could provide a trajectory in which the aircraft flies by the waypoint more than once. This situation can be easily avoid introducing in the model the following simple constraint:(40)∑i=1N−1κi1≤c3,where c3∈1,2,⋯,N−1 is a suitable constant. In the context of the problem considered in this work including the penalty term (39), a value of c3 between 2 and 4 is enough to avoid this potential undesired issue.Note that this waypoint constraint setting can be straightforward extended to multiwaypoint and multiaircraft settings. Moreover, once again, since this embedded logical constraint approach is quite general, any other waypoint constraint described in terms of Equation (17) can be considered.


## 5. Numerical Results

To show the effectiveness of the methodology describe in Section4, the following three numerical experiments have been carried out:(i) Experiment 1. Minimum-time continuous descent of three aircraft along routes converging at the same waypoint with time-based separation constraint among them(ii) Experiment 2. Minimum-time continuous descent of three aircraft along intersecting routes with distance-based separation constraint(iii) Experiment 3. Minimum-time STAR-based continuous descent of three aircraft along converging routes. The considered STAR procedure includes sequencing the aircraft at a merging point and passing through two other waypointsAll of these numerical experiments involve Airbus A-320 BADA 3.6 aircraft models in which the performance index is the sum of the duration of the flights of the three aircraft. In some of them, constraints derived from current flight regulation have been introduced, namely, time-based separation and distance-based separation operational constraints have been imposed. In particular, a minimum horizontal distance separation between aircraft of 5000 m and vertical separation between aircraft of 1000 m has been considered, whereas the considered minimum time separation between aircraft has been 200 s. These specific values have been chosen taking as a reference the aviation regulation [20] in which, in general, aircraft have to be separated by at least 3 NM or 3 min in the TMA.The numerical experiments have been conducted on an Intel Core i7 2.8 GHz CPU with 16 GB RAM.


### 5.1. Experiment 1: Minimum-Time Continuous Descent with Time-Based Separation Constraints

In this experiment, a CDA under vectoring has been considered, that is, the lateral path followed by the aircraft has been assumed to be specified through instructions provided by the ATC. In particular, the boundary conditions of the state variables have been selected from the chart of the Adolfo Suárez Madrid-Barajas (LEMD/MAD) TMA shown in Figure3. The initial position of Aircraft 1, Aircraft 2, and Aircraft 3 has been supposed to be coincident with the ROLDO, SOTUK, and MORAL waypoints, respectively. Their common final position has been assumed to be the LALPI waypoint. The three aircraft have been supposed to go directly from their initial positions to the final waypoint, that is, no constraints on the lateral path have been considered. Neither time-based nor distance-based separation constraints have been included in the model. The initial mass of the three aircraft has been assumed equal to the maximum landing weight of the aircraft. The specific boundary conditions of the state variables are given in Table 1 and the final mass and time given by the solution of the corresponding OCP are reported for each aircraft in the first row of Table 2.Figure 3
Chart of the Adolfo Suárez Madrid-Barajas (LEMD/MAD) STAR 10-2A1.Table 1
Boundary conditions for Experiment 1 and Experiment 3.
Symbol
Unit
Aircraft 1
Aircraft 2
Aircraft 3
hI
m
7400
7000
7200
hF
m
3350
3350
3350
θI
deg
39.526
39.116
39.000
θF
deg
40.575
40.575
40.575
λI
deg
-5.327
-4.448
-3.325
λF
deg
-3.422
-3.422
-3.422
VI
m/s
130
130
130
VF
m/s
110
110
110
μI
deg
0
0
0
γI
deg
0
0
0
χI
deg
356
294
24
mI
kg
65000
65000
65000Table 2
Results of Experiment 1.
# aircraft
Final time (s)
Final mass (kg)
Without time constraint
1
1539
64023
2
1419
64142
3
1350
64151
With time constraint
1
1750
64109
2
1550
64216
3
1350
64151In Figure4, the 3D view of the paths obtained in the solution is represented in thick lines, whereas in Figures 5 and 6, the horizontal and vertical profiles are represented in solid lines.Figure 4
Experiment 1. 3D view of the paths with (thin lines) and without (thick lines) time-based separation constraints.Figure 5
Experiment 1. Horizontal profiles with (TC, dashed lines) and without (FF, solid lines) time-based separation constraints.Figure 6
Experiment 1. Vertical profiles with (TC, dashed lines) and without (FF, solid lines) time-based separation constraints.The computation time to find the solution has been15.8098 s. The final time of the three aircraft is 1539 s, 1419 s, and 1350 s, respectively, corresponding to 22-25 min of flight. The difference between the final time of aircraft 2 and 3 is 69 s. The difference between the final time of aircraft 1 and 2 is 120 s. The difference between the final time of aircraft 1 and 3 is 189 s. These differences in time among the three aircraft at LALPI waypoint are below the required time separation of 200 s established above. Therefore, a new experiment has been conducted, which includes time-based separation logical constraint of at least 200 s between aircraft as described in Section 4.1.1. The final mass and time given by the solution are reported for each aircraft in the second row of Table 2.In Figure4, the 3D view of the paths is represented in thin lines, and in Figures 5 and 6, the horizontal and vertical profiles are represented in dashed lines. It can be seen that whereas there are slight changes in the horizontal profiles, the vertical profiles changed significantly.In this case, the computation time to find the solution has been20.8883 s, and the final time of the three aircraft is 1750 s, 1550 s, and 1350 s, respectively. Since the objective function to be minimized is the duration of the flights, aircraft maintain mutually the minimum required time separation of 200 s. The mass consumption is represented in Figure 7, in which the solution without time-based separation constraints is represented in solid lines, whereas the solution with time-based separation constraints is represented in dashed lines. It can be seen that introducing time separation constraints increases the mass consumption of Aircraft 1 and Aircraft 2, whereas that of Aircraft 3 remains unchanged.Figure 7
Experiment 1. Mass consumption with (TC, dashed lines) and without (FF, solid lines) time-based separation constraints.


### 5.2. Experiment 2: Minimum-Time Continuous Descent of Three Aircraft along Intersecting Routes with Distance-Based Separation Constraint

This experiment has been designed to test the effectiveness of the logical constraint formalism described in Section4.1.2, where three aircraft start at and have to reach three different positions flying along intersecting routes in which a mid-air conflict arises, in the sense that the experiment is intentionally designed in such a way that they reach almost the same point at the same time.Neither time-based nor distance-based separation constraints have been included in the model. The initial mass of the three aircraft has been assumed equal to the maximum landing weight of the aircraft. The specific boundary conditions of the state variables are given in Table3, the information about the positions of the aircraft at the conflict region is reported in Table 4, and the final mass and time given by the solution of the corresponding OCP are reported for each aircraft in the first row of Table 5.Table 3
Boundary conditions for Experiment 2.
Symbol
Unit
Aircraft 1
Aircraft 2
Aircraft 3
hI
m
7400
7400
7400
hF
m
3350
3350
3350
θI
deg
39.000
39.000
38.701
θF
deg
40.575
40.575
40.848
λI
deg
-5.327
-3.325
-4.320
λF
deg
-3.430
-5.215
-4.320
VI
m/s
130
130
130
VF
m/s
110
110
110
μI
deg
0
0
0
γI
deg
0
0
0
χI
deg
356
294
294
mI
kg
65000
65000
65000Table 4
Positions of the aircraft at conflict region (at time 980 s) without distance constraint for Experiment 2.
Symbol
Unit
Aircraft 1
Aircraft 2
Aircraft 3
hI
m
4176
4170
4194
θI
deg
39.8414
39.8449
39.8458
λI
deg
-4.3243
-4.3165
-4.3184Table 5
Results of Experiment 2.
# aircraft
Final time (s)
Final mass (kg)
Without distance constraint
1
1838
63787
2
1838
63824
3
1838
63751
With distance constraint
1
1841
63772
2
1839
63900
3
1843
63772In Figure8, the 3D view of the paths obtained in the solution is represented in thick lines, whereas Figures 9 and 10, the horizontal and vertical profiles are represented in solid lines. The computation time to find the solution has been 65.6747 s, and the final time of three aircraft is 1838 s.Figure 8
Experiment 2. 3D view of the paths with (thin lines) and without (thick lines) distance-based separation constraints.Figure 9
Experiment 2. Horizontal profiles with (DC, dashed lines) and without (FF, solid lines) distance-based separation constraints.Figure 10
Experiment 2. Vertical profiles with (DC, dashed lines) and without (FF, solid lines) distance-based separation constraints.To avoid the conflict, a new experiment has been conducted, which includes distance-based separation logical constraints of at least 5000 m (2.7 NM) in the horizontal profile and at least 1000 m (0.5 NM) in the vertical profile as described in Section 4.1.2. The final mass and time obtained in the solution are reported for each aircraft in the second row of Table 5. In Figure 8, the 3D view of the paths is represented in thin lines, whereas in Figures 9 and 10, the horizontal and vertical profiles are represented in dashed lines.It can be seen from Figures9 and 11 that the conflict is avoided by activating the logical constraints associated with the horizontal profile. In particular, from Figure 11, it can be observed that the logical distance constraint between Aircraft 1 and Aircraft 2 is active in the interval 982,1012 s, between Aircraft 2 and Aircraft 3 it is active in the interval 1011,1042 s, and between Aircraft 1 and Aircraft 3 it is active in the interval 952,983 s. The final time of Aircraft 1, Aircraft 2, and Aircraft 3 is similar, being 1841 s, 1838 s, and 1843 s, respectively. The computation time to find the solution has been 69.3823 s. The mass consumption is depicted in Figure 12, in which the solution without distance-based separation constraints is represented in solid lines, whereas the solution with distance-based separation constraints is represented in dashed lines. It can be seen that introducing distance-based separation constraints does not significantly change the mass consumption of the aircraft.Figure 11
Experiment 2. Distance among aircraft with distance-based separation constraints (DC).Figure 12
Experiment 2. Mass consumption with (DC, dashed lines) and without (FF, solid lines) distance-based separation constraints.


### 5.3. Experiment 3: Minimum-Time STAR-Based Continuous Descent along Converging Routes

In this experiment, a STAR-based CDA has been considered, that is, the lateral path followed by the aircraft has been assumed to be specified in a navigation chart. In particular, as in Experiment 1, the boundary conditions of the state variables have been selected from the chart of the Adolfo Suárez Madrid-Barajas (LEMD/MAD) TMA shown in Figure3. The initial position of Aircraft 1, Aircraft 2, and Aircraft 3 is supposed to be coincident with the ROLDO, SOTUK, and MORAL waypoints, respectively. Aircraft are constrained to pass through TODNO and RESBI waypoints and their common final position is assumed to be the LALPI waypoint. For the setting of the cuboids centered at the given waypoints, the modeling defined in Equation (34) has been considered. In particular, the cuboid centered at TODNO waypoint has been defined by the two corners 39.560∘,−4.24∘,5400m and 39.640∘,−4.160∘,6000m, whereas the cuboid centered at RESBI waypoint has been defined by the two corners 40.400∘,−4.150∘,4200m and 40.480∘,−4.070∘,4800m.The initial mass of the three aircraft has been assumed equal to the maximum landing weight of the aircraft. The specific boundary conditions of the state variables are the same as in Experiment 1 and are given in Table1, and the final mass and time given by the solution of the corresponding OCP are reported for each aircraft in the second row of Table 6. In spite of the fact that the three aircraft have to pass through two waypoints, neither distance-based conflict nor time-based conflict arises during the flight except for the last part of the time discretization. In particular, Aircraft 1 and Aircraft 3 reach the LALPI waypoint with a time separation of 6 s only. The computation time to find the solution has been 41.7981 s.Table 6
Results of Experiment 3.
# aircraft
Final time (s)
Final mass (kg)
Without logical constraints
1
1539
64023
2
1419
64142
3
1350
64151
With waypoint constraints
1
1971
63813
2
1626
64071
3
1965
63863
With waypoint and time constraints
1
2109
63908
2
1619
64127
3
1909
63956Therefore, a new experiment has been conducted, which, besides waypoint constraints, also includes time-based separation logical constraint of at least 200 s between aircraft as described in Section4.1.1. The final mass and time given by the solution are reported for each aircraft in the third row of Table 6. In Figure 13, the 3D view of the paths is represented in thin lines, whereas in Figures 14 and 15, the horizontal and vertical profiles are represented in dashed lines.Figure 13
Experiment 3. 3D view of the paths with (thin lines) and without (thick lines) waypoints and time-based separation constraints.Figure 14
Experiment 3. Horizontal profiles with (WT, dashed lines) and without (FF, solid lines) waypoints and time-based separation constraints.Figure 15
Experiment 3. Vertical profiles with (WT, dashed lines) and without (FF, solid lines) waypoints and time-based separation constraints.It can be seen that, as expected, there are significant changes in both the horizontal and vertical profiles. The final time of the three aircraft is 2109 s, 1619 s, and 1909 s, respectively. The difference between the final time of aircraft 1 and 3 is 200 s, whereas the difference between the final time of aircraft 3 and 2 is 290 s. Notice that, unlike Experiment 1, aircraft 3 and 2 do not maintain the minimum required time separation of 200 s. This is due to the fact that the performance index includes in this case, besides the minimization of the duration of the flights, the penalty term associated with the waypoints constraints. Therefore, the saturation of the time-based constraints may not happen in this setting. In this case, the computation time to find the solution has been49.0766 s.The mass consumption is depicted in Figure16, in which the solution without logical constraints is represented in solid lines, whereas the solution with waypoints constraints is represented in dashed lines. It can be seen that introducing waypoint constraints increases the mass consumption of Aircraft 1 and Aircraft 3, whereas that of Aircraft 2 remains almost unchanged.Figure 16
Experiment 3. Mass consumption with (WT, dashed lines) and without (FF, solid lines) waypoints and time-based separation constraints.


## 5.1. Experiment 1: Minimum-Time Continuous Descent with Time-Based Separation Constraints

In this experiment, a CDA under vectoring has been considered, that is, the lateral path followed by the aircraft has been assumed to be specified through instructions provided by the ATC. In particular, the boundary conditions of the state variables have been selected from the chart of the Adolfo Suárez Madrid-Barajas (LEMD/MAD) TMA shown in Figure3. The initial position of Aircraft 1, Aircraft 2, and Aircraft 3 has been supposed to be coincident with the ROLDO, SOTUK, and MORAL waypoints, respectively. Their common final position has been assumed to be the LALPI waypoint. The three aircraft have been supposed to go directly from their initial positions to the final waypoint, that is, no constraints on the lateral path have been considered. Neither time-based nor distance-based separation constraints have been included in the model. The initial mass of the three aircraft has been assumed equal to the maximum landing weight of the aircraft. The specific boundary conditions of the state variables are given in Table 1 and the final mass and time given by the solution of the corresponding OCP are reported for each aircraft in the first row of Table 2.Figure 3
Chart of the Adolfo Suárez Madrid-Barajas (LEMD/MAD) STAR 10-2A1.Table 1
Boundary conditions for Experiment 1 and Experiment 3.
Symbol
Unit
Aircraft 1
Aircraft 2
Aircraft 3
hI
m
7400
7000
7200
hF
m
3350
3350
3350
θI
deg
39.526
39.116
39.000
θF
deg
40.575
40.575
40.575
λI
deg
-5.327
-4.448
-3.325
λF
deg
-3.422
-3.422
-3.422
VI
m/s
130
130
130
VF
m/s
110
110
110
μI
deg
0
0
0
γI
deg
0
0
0
χI
deg
356
294
24
mI
kg
65000
65000
65000Table 2
Results of Experiment 1.
# aircraft
Final time (s)
Final mass (kg)
Without time constraint
1
1539
64023
2
1419
64142
3
1350
64151
With time constraint
1
1750
64109
2
1550
64216
3
1350
64151In Figure4, the 3D view of the paths obtained in the solution is represented in thick lines, whereas in Figures 5 and 6, the horizontal and vertical profiles are represented in solid lines.Figure 4
Experiment 1. 3D view of the paths with (thin lines) and without (thick lines) time-based separation constraints.Figure 5
Experiment 1. Horizontal profiles with (TC, dashed lines) and without (FF, solid lines) time-based separation constraints.Figure 6
Experiment 1. Vertical profiles with (TC, dashed lines) and without (FF, solid lines) time-based separation constraints.The computation time to find the solution has been15.8098 s. The final time of the three aircraft is 1539 s, 1419 s, and 1350 s, respectively, corresponding to 22-25 min of flight. The difference between the final time of aircraft 2 and 3 is 69 s. The difference between the final time of aircraft 1 and 2 is 120 s. The difference between the final time of aircraft 1 and 3 is 189 s. These differences in time among the three aircraft at LALPI waypoint are below the required time separation of 200 s established above. Therefore, a new experiment has been conducted, which includes time-based separation logical constraint of at least 200 s between aircraft as described in Section 4.1.1. The final mass and time given by the solution are reported for each aircraft in the second row of Table 2.In Figure4, the 3D view of the paths is represented in thin lines, and in Figures 5 and 6, the horizontal and vertical profiles are represented in dashed lines. It can be seen that whereas there are slight changes in the horizontal profiles, the vertical profiles changed significantly.In this case, the computation time to find the solution has been20.8883 s, and the final time of the three aircraft is 1750 s, 1550 s, and 1350 s, respectively. Since the objective function to be minimized is the duration of the flights, aircraft maintain mutually the minimum required time separation of 200 s. The mass consumption is represented in Figure 7, in which the solution without time-based separation constraints is represented in solid lines, whereas the solution with time-based separation constraints is represented in dashed lines. It can be seen that introducing time separation constraints increases the mass consumption of Aircraft 1 and Aircraft 2, whereas that of Aircraft 3 remains unchanged.Figure 7
Experiment 1. Mass consumption with (TC, dashed lines) and without (FF, solid lines) time-based separation constraints.


## 5.2. Experiment 2: Minimum-Time Continuous Descent of Three Aircraft along Intersecting Routes with Distance-Based Separation Constraint

This experiment has been designed to test the effectiveness of the logical constraint formalism described in Section4.1.2, where three aircraft start at and have to reach three different positions flying along intersecting routes in which a mid-air conflict arises, in the sense that the experiment is intentionally designed in such a way that they reach almost the same point at the same time.Neither time-based nor distance-based separation constraints have been included in the model. The initial mass of the three aircraft has been assumed equal to the maximum landing weight of the aircraft. The specific boundary conditions of the state variables are given in Table3, the information about the positions of the aircraft at the conflict region is reported in Table 4, and the final mass and time given by the solution of the corresponding OCP are reported for each aircraft in the first row of Table 5.Table 3
Boundary conditions for Experiment 2.
Symbol
Unit
Aircraft 1
Aircraft 2
Aircraft 3
hI
m
7400
7400
7400
hF
m
3350
3350
3350
θI
deg
39.000
39.000
38.701
θF
deg
40.575
40.575
40.848
λI
deg
-5.327
-3.325
-4.320
λF
deg
-3.430
-5.215
-4.320
VI
m/s
130
130
130
VF
m/s
110
110
110
μI
deg
0
0
0
γI
deg
0
0
0
χI
deg
356
294
294
mI
kg
65000
65000
65000Table 4
Positions of the aircraft at conflict region (at time 980 s) without distance constraint for Experiment 2.
Symbol
Unit
Aircraft 1
Aircraft 2
Aircraft 3
hI
m
4176
4170
4194
θI
deg
39.8414
39.8449
39.8458
λI
deg
-4.3243
-4.3165
-4.3184Table 5
Results of Experiment 2.
# aircraft
Final time (s)
Final mass (kg)
Without distance constraint
1
1838
63787
2
1838
63824
3
1838
63751
With distance constraint
1
1841
63772
2
1839
63900
3
1843
63772In Figure8, the 3D view of the paths obtained in the solution is represented in thick lines, whereas Figures 9 and 10, the horizontal and vertical profiles are represented in solid lines. The computation time to find the solution has been 65.6747 s, and the final time of three aircraft is 1838 s.Figure 8
Experiment 2. 3D view of the paths with (thin lines) and without (thick lines) distance-based separation constraints.Figure 9
Experiment 2. Horizontal profiles with (DC, dashed lines) and without (FF, solid lines) distance-based separation constraints.Figure 10
Experiment 2. Vertical profiles with (DC, dashed lines) and without (FF, solid lines) distance-based separation constraints.To avoid the conflict, a new experiment has been conducted, which includes distance-based separation logical constraints of at least 5000 m (2.7 NM) in the horizontal profile and at least 1000 m (0.5 NM) in the vertical profile as described in Section 4.1.2. The final mass and time obtained in the solution are reported for each aircraft in the second row of Table 5. In Figure 8, the 3D view of the paths is represented in thin lines, whereas in Figures 9 and 10, the horizontal and vertical profiles are represented in dashed lines.It can be seen from Figures9 and 11 that the conflict is avoided by activating the logical constraints associated with the horizontal profile. In particular, from Figure 11, it can be observed that the logical distance constraint between Aircraft 1 and Aircraft 2 is active in the interval 982,1012 s, between Aircraft 2 and Aircraft 3 it is active in the interval 1011,1042 s, and between Aircraft 1 and Aircraft 3 it is active in the interval 952,983 s. The final time of Aircraft 1, Aircraft 2, and Aircraft 3 is similar, being 1841 s, 1838 s, and 1843 s, respectively. The computation time to find the solution has been 69.3823 s. The mass consumption is depicted in Figure 12, in which the solution without distance-based separation constraints is represented in solid lines, whereas the solution with distance-based separation constraints is represented in dashed lines. It can be seen that introducing distance-based separation constraints does not significantly change the mass consumption of the aircraft.Figure 11
Experiment 2. Distance among aircraft with distance-based separation constraints (DC).Figure 12
Experiment 2. Mass consumption with (DC, dashed lines) and without (FF, solid lines) distance-based separation constraints.


## 5.3. Experiment 3: Minimum-Time STAR-Based Continuous Descent along Converging Routes

In this experiment, a STAR-based CDA has been considered, that is, the lateral path followed by the aircraft has been assumed to be specified in a navigation chart. In particular, as in Experiment 1, the boundary conditions of the state variables have been selected from the chart of the Adolfo Suárez Madrid-Barajas (LEMD/MAD) TMA shown in Figure3. The initial position of Aircraft 1, Aircraft 2, and Aircraft 3 is supposed to be coincident with the ROLDO, SOTUK, and MORAL waypoints, respectively. Aircraft are constrained to pass through TODNO and RESBI waypoints and their common final position is assumed to be the LALPI waypoint. For the setting of the cuboids centered at the given waypoints, the modeling defined in Equation (34) has been considered. In particular, the cuboid centered at TODNO waypoint has been defined by the two corners 39.560∘,−4.24∘,5400m and 39.640∘,−4.160∘,6000m, whereas the cuboid centered at RESBI waypoint has been defined by the two corners 40.400∘,−4.150∘,4200m and 40.480∘,−4.070∘,4800m.The initial mass of the three aircraft has been assumed equal to the maximum landing weight of the aircraft. The specific boundary conditions of the state variables are the same as in Experiment 1 and are given in Table1, and the final mass and time given by the solution of the corresponding OCP are reported for each aircraft in the second row of Table 6. In spite of the fact that the three aircraft have to pass through two waypoints, neither distance-based conflict nor time-based conflict arises during the flight except for the last part of the time discretization. In particular, Aircraft 1 and Aircraft 3 reach the LALPI waypoint with a time separation of 6 s only. The computation time to find the solution has been 41.7981 s.Table 6
Results of Experiment 3.
# aircraft
Final time (s)
Final mass (kg)
Without logical constraints
1
1539
64023
2
1419
64142
3
1350
64151
With waypoint constraints
1
1971
63813
2
1626
64071
3
1965
63863
With waypoint and time constraints
1
2109
63908
2
1619
64127
3
1909
63956Therefore, a new experiment has been conducted, which, besides waypoint constraints, also includes time-based separation logical constraint of at least 200 s between aircraft as described in Section4.1.1. The final mass and time given by the solution are reported for each aircraft in the third row of Table 6. In Figure 13, the 3D view of the paths is represented in thin lines, whereas in Figures 14 and 15, the horizontal and vertical profiles are represented in dashed lines.Figure 13
Experiment 3. 3D view of the paths with (thin lines) and without (thick lines) waypoints and time-based separation constraints.Figure 14
Experiment 3. Horizontal profiles with (WT, dashed lines) and without (FF, solid lines) waypoints and time-based separation constraints.Figure 15
Experiment 3. Vertical profiles with (WT, dashed lines) and without (FF, solid lines) waypoints and time-based separation constraints.It can be seen that, as expected, there are significant changes in both the horizontal and vertical profiles. The final time of the three aircraft is 2109 s, 1619 s, and 1909 s, respectively. The difference between the final time of aircraft 1 and 3 is 200 s, whereas the difference between the final time of aircraft 3 and 2 is 290 s. Notice that, unlike Experiment 1, aircraft 3 and 2 do not maintain the minimum required time separation of 200 s. This is due to the fact that the performance index includes in this case, besides the minimization of the duration of the flights, the penalty term associated with the waypoints constraints. Therefore, the saturation of the time-based constraints may not happen in this setting. In this case, the computation time to find the solution has been49.0766 s.The mass consumption is depicted in Figure16, in which the solution without logical constraints is represented in solid lines, whereas the solution with waypoints constraints is represented in dashed lines. It can be seen that introducing waypoint constraints increases the mass consumption of Aircraft 1 and Aircraft 3, whereas that of Aircraft 2 remains almost unchanged.Figure 16
Experiment 3. Mass consumption with (WT, dashed lines) and without (FF, solid lines) waypoints and time-based separation constraints.


## 6. Conclusions

In this paper, the trajectory planning problem for multiple aircraft has been studied in which logical constraints in disjunctive form are included in the model. The logical constraints in disjunctive form have been transformed into inequality and equality constraints which involves only continuous auxiliary variables. In this way, the optimal control problem with logical constraints has been converted into a smooth optimal control problem which has been solved using standard techniques. This approach has been applied to the computation of the optimized profile descent of multiple aircraft in converging and intersecting arrival routes within the Adolfo Suárez Madrid-Barajas (LEMD/MAD) TMA. The results show the effectiveness of the proposed technique.



---

*Source: 5158967-2019-12-03.xml*