International Journal of Aerospace Engineering
Volume 2019 (2019), Article ID 5158967, 15 pages
https://doi.org/10.1155/2019/5158967

Research Article

Multiaircraft Optimal 4D Trajectory Planning Using Logical Constraints

Dinesh B. Seenivasan, Alberto Olivares, and Ernesto Staffetti

Department of Signal Theory and Communications and Telematic Systems and Computing, Universidad Rey Juan Carlos, Camino del Molino s/n 28943 Fuenlabrada, Madrid, Spain

Correspondence should be addressed to Dinesh B. Seenivasan; seenivasan.dinesh@urjc.es

Received 15 January 2019; Accepted 19 September 2019; Published 3 December 2019

Academic Editor: Jose Carlos Páscoa

Copyright © 2019 Dinesh B. Seenivasan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 Section 2, 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

for , where describes the right-hand side of the differential equation

and describes the algebraic constraints

where and are the state and control sets, respectively, is a -dimensional state variable, is a -dimensional control input, is a vector of parameters, and represents time, in which and denote the initial time and final time for aircraft p, with .

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:

where

The objective function

is given in Bolza form. It is expressed as a combination of a Mayer term

and a Lagrange term

Functions

are assumed to be twice differentiable. Function is assumed to be piecewise Lipschitz continuous within the time interval , and the derivative of the algebraic right-hand side function with respect to , that is,

is assumed to be regular within the time interval , where denotes the -dimensional algebraic variable, that is, the state variables without time derivative. Vector represents the initial conditions given at the initial time and function

provides the terminal conditions at the final time , and it is assumed to be twice differentiable. The system must also satisfy algebraic path constraints within the time interval given by the vector function

with lower bound and upper bound . 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 such that the set of aircraft follow an admissible trajectory between the initial and final state that minimizes the performance index . The final time, , 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 has been subdivided into subintervals of equal length, whose endpoints are

with and . In each subinterval , , 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 , .

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:

The three dynamic equations in (14) are expressed in an aircraft-attached reference frame , and the three kinematic equations are expressed in a ground-based reference frame as shown in Figure 1.

(a) Top view

(b) Front view

(c) Lateral view

Figure 1: Aircraft state and forces.

The states of the system (14) are , and . Thus, . The state variables ,, and are the true airspeed, the heading angle, and the flight path angle, respectively. The state variables , and are the aircraft 3D position, the longitude, the latitude, and the altitude, respectively. The components of the aircraft position vector in two dimensions (, ) are approximated as and . Finally, the state variable is the aircraft mass. The control inputs are the bank angle , the engine thrust , and the lift coefficient . Thus, .

The parameter is the radius of Earth and is the speed-dependent fuel efficiency coefficient. Lift, , and drag, , are the components of the aerodynamic force. The parameter is the reference wing surface area and is the dynamic pressure. A parabolic drag polar and an International Standard Atmosphere (ISA) model are assumed. The lift coefficient 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]:

In (15), is the maximum reachable altitude and is the maximum operative altitude at a given mass (it increases as fuel is burned). is the Mach number and is the maximum operating Mach number. is the minimum speed coefficient, is the stall speed, is the maximum operating CAS, and and are, respectively, the maximum normal and longitudinal accelerations for civilian aircraft. Finally, and 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

where

is in disjunctive form, in which only one of the propositions must be satisfied and the proposition is either or . Term is a literal that can be either True or False and represents the negation or logical complement operator. The term is used to represent statements such as “longitude ”. Therefore, takes the form

where is assumed to be a 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:

Thus, taking into account (17), the logical expression (16) can be represented as follows:

To transcribe these disjunctions into a set of inequality constraints, a continuous variable is defined and related to each in (18). Thus, (20) can be expressed as follows:

It is immediate to check that if in the first term in (21), then constraint is not fulfilled. On the other hand, if , then is in fact , and thus, constraint is enforced. Finally, the last term in (21) guarantees that at least one of the propositions 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.

Let and be the unfixed time at the merging waypoint of aircraft and , respectively, and let be the safety time difference between two consecutive aircraft. Since in terms of the discretization (13), and , then multiaircraft time-based separation constraints can be expressed as follows:

where condition prevents unnecessary duplication of constraints. Constraints (22) can be rewritten as follows:

Following the technique described above to tackle constraints in disjunctive form, if we define new variables satisfying condition

Equation (23) can be transformed into

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 and 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 directions.

Let and be the positions of aircraft and at the endpoint of the discretization, respectively. If the safety distance-based separation in the and directions is denoted by and , respectively, then the distance-based separation constraints can be expressed as follows:

where the haversine formula has been considered for the distance in the direction with

As before, condition prevents unnecessary duplication of constraints, and the set of constraints (26) can be rewritten as follows:

Again, following the technique described above, new variables for , satisfying condition

are introduced, and Equation (28) can be transformed into

The last constraint in (30) ensures that at least one of the constraints in (28) is fulfilled, that is, aircraft and 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 and be the positions of opposite corners of the cuboid. Cuboid avoidance for a single aircraft involves that at every endpoint of the subintervals of the discretization(13), the position of the aircraft must remain outside it. In terms of logical constraints, this condition can be expressed as follows:

Note that constraints (31) are not enforced at the initial and final points of the discretization, and , 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 , for , satisfying condition

are introduced, and Equation (31) can be transformed into

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 Section 4.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 and 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 of the subintervals of the discretization (13), the position of the aircraft must remain inside it. In terms of logical constraints, this condition can be expressed as follows:

Note that, for one hand, constraints (34) are enforced at every point of the discretization except for the initial and final points, and , 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 , and be the minimum and maximum values of the state variables , and , respectively. In terms of logical constraints, the free flight mode condition can be expressed as follows:

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, , can be expressed as follows:

where and 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 satisfying condition

Equation (36) can be transformed into

The last constraint in (38) ensures that at least one of the conditions in (36) is fulfilled, which means that at each discretization instant , 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 instant . 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

where and are suitable constants determined by the continuation method. In the context of a minimization performance index such as (4a), a large enough value such that , 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:

where is a suitable constant. In the context of the problem considered in this work including the penalty term (39), a value of 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 Section 4, 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 waypoints

All 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 Figure 3. 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

m
7400
7000
7200

m
3350
3350
3350

deg
39.526
39.116
39.000

deg
40.575
40.575
40.575

deg
-5.327
-4.448
-3.325

deg
-3.422
-3.422
-3.422

m/s
130
130
130

m/s
110
110
110

deg
0
0
0

deg
0
0
0

deg
356
294
24

kg
65000
65000
65000

Table 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
64151

In Figure 4, 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 been  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 Figure 4, 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 been  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 Section 4.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 Table 3, 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

m
7400
7400
7400

m
3350
3350
3350

deg
39.000
39.000
38.701

deg
40.575
40.575
40.848

deg
-5.327
-3.325
-4.320

deg
-3.430
-5.215
-4.320

m/s
130
130
130

m/s
110
110
110

deg
0
0
0

deg
0
0
0

deg
356
294
294

kg
65000
65000
65000

Table 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

m
4176
4170
4194

deg
39.8414
39.8449
39.8458

deg
-4.3243
-4.3165
-4.3184

Table 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
63772

In Figure 8, 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  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 ( NM) in the horizontal profile and at least 1000 m ( 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 Figures 9 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  s, between Aircraft 2 and Aircraft 3 it is active in the interval  s, and between Aircraft 1 and Aircraft 3 it is active in the interval  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  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 Figure 3. 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 and , whereas the cuboid centered at RESBI waypoint has been defined by the two corners and .

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 Table 1, 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  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
63956

Therefore, 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 Section 4.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 been  s.

The mass consumption is depicted in Figure 16, 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.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Funding

This work has been supported by the Spanish grants TRA2013-47619-C2-2-R and TRA2017-91203-EXP.

References

A. Gardi, R. Sabatini, and S. Ramasamy, “Multi-objective optimisation of aircraft flight trajectories in the ATM and avionics context,” Progress in Aerospace Sciences, vol. 83, pp. 1–36, 2016.

A. Andreeva-Mori, S. Suzuki, and E. Itoh, “Rule derivation for arrival aircraft sequencing,” Aerospace Science and Technology, vol. 30, no. 1, pp. 200–209, 2013.

P. Bonami, A. Olivares, M. Soler, and E. Staffetti, “Multiphase mixed-integer optimal control approach to aircraft trajectory optimization,” Journal of Guidance, Control, and Dynamics, vol. 36, no. 5, pp. 1267–1277, 2013.

A. U. Raghunathan, V. Gopal, D. Subramanian, L. T. Biegler, and T. Samad, “Dynamic optimization strategies for three-dimensional conflict resolution of multiple aircraft,” Journal of Aircraft, vol. 27, no. 4, pp. 586–594, 2004.

S. Wei, M. Zěfran, and R. A. DeCarlo, “Optimal control of robotic systems with logical constraints: Application to UAV path planning,” in 2008 IEEE International Conference on Robotics and Automation, pp. 176–181, Pasadena, CA, USA, 2008.

K. Mohan, M. A. Patterson, and A. V. Rao, “Optimal trajectory and control generation for landing of multiple aircraft in the presence of obstacles,” in AIAA Guidance, Navigation, and Control Conference, pp. 1–6, Minneapolis, MN, USA, 2012.

ICAO, “Continuous descent operations (CDO) manual,” Tech. Rep., Tech. Rep. Document 9931, International Civil Aviation Organization, 2010.

E. Itoh, M. Brown, A. Senoguchi, N. Wickramasinghe, and S. Fukushima, “Future arrival management collaborating with trajectory-based operations,” in Air Traffic Management and Systems II, Electronic Navigation Research Institute, Ed., vol. 420 of Lecture Notes in Electrical Engineering, pp. 137–156, Springer, 2017.

S. G. Park and J.-P. Clarke, “Optimal control based vertical trajectory determination for continuous descent arrival procedures,” Journal of Aircraft, vol. 52, no. 5, pp. 1469–1480, 2015.

S. G. Park and J. P. Clarke, “Vertical trajectory optimization to minimize environmental impact in the presence of wind,” Journal of Aircraft, vol. 53, no. 3, pp. 725–737, 2016.

C. R. Hargraves and S. W. Paris, “Direct trajectory optimization using nonlinear programming and collocation,” Journal of Guidance, Control, and Dynamics, vol. 10, no. 4, pp. 338–342, 1987.

A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical Programming, vol. 106, no. 1, pp. 25–57, 2006.

M. Soler, A. Olivares, and E. Staffetti, “Multiphase optimal control framework for commercial aircraft four-dimensional flight-planning problems,” Journal of Aircraft, vol. 52, no. 1, pp. 274–286, 2015.

A. Nuic, User Manual for the Base of Aircraft Data (BADA) Revision 3.6, EUROCONTROL Experimental Center, Bretigny, France, 2005.

S. C. Bengea and R. A. DeCarlo, “Optimal control of switching systems,” Automatica, vol. 41, no. 1, pp. 11–27, 2005.

S. Wei, K. Uthaichana, M. Žefran, R. A. DeCarlo, and S. Bengea, “Applications of numerical optimal control to nonlinear hybrid systems,” Nonlinear Analysis: Hybrid Systems, vol. 1, no. 2, pp. 264–279, 2007.

T. M. Cavalier, P. M. Pardalos, and A. L. Soyster, “Modeling and integer programming techniques applied to propositional calculus,” Computers and Operations Research, vol. 17, no. 6, pp. 561–570, 1990.

S. Hao, S. Cheng, and Y. Zhang, “A multi-aircraft conflict detection and resolution method for 4-dimensional trajectory-based operation,” Chinese Journal of Aeronautics, vol. 31, no. 7, pp. 1579–1593, 2018.

S. Wei, K. Uthaichana, M. Zefran, and R. DeCarlo, “Hybrid model predictive control for the stabilization of wheeled mobile robots subject to wheel slippage,” IEEE Transactions on Control Systems Technology, vol. 21, no. 6, pp. 2181–2193, 2013.

ICAO, “Procedures for Air Navigation Services — Air Traffic Management,” Tech. Rep., Tech. Rep. 4444, International Civil Aviation Organization, 2016.