# Multiobjective Simulated Annealing for Collision Avoidance in ATM Accounting for Three Admissible Maneuvers

Authors: A. Mateos; A. Jiménez-Martín
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Technological advances are required to accommodate air traffic control systems for the future growth of air traffic. Particularly, detection and resolution of conflicts between aircrafts is a problem that has attracted much attention in the last decade becoming vital to improve the safety standards in free flight unstructured environments. We propose using the archive simulated annealing-based multiobjective optimization algorithm to deal with such a problem, accounting for three admissible maneuvers (velocity, turn, and altitude changes) in a multiobjective context. The minimization of the maneuver number and magnitude, time delays, or deviations in the leaving points are considered for analysis. The optimal values for the algorithm parameter set are identified in the more complex instance in which all aircrafts have conflicts between each other accounting for 5, 10, and 20 aircrafts. Moreover, the performance of the proposed approach is analyzed by means of a comparison with the Pareto front, computed using brute force for 5 aircrafts and the algorithm is also illustrated with a random instance with 20 aircrafts.


---

## Body

## 1. Introduction

Cargo and air traffic (AT) congestion has experienced a general exponential growth throughout the world over the last decade. Every minute of the day, both morning and afternoon, there are about 11,000 aircrafts in the air somewhere in the world, as can be seen in real time athttps://www.flightradar24.com/.2014 was the first year in which 100,000 flights per day were exceeded. Europe’s largest airports handle about 2,000 daily takeoffs and landings. This trend continues to increase gradually and estimates predict bending movements until 2030.With the systems currently available, the air traffic control agencies are not able to efficiently manage this large increase which is taking place due to several factors as follows:(1)
Efficient use of airspace: currently, the airspace is rigidly structured and with a large number of constraints that aircrafts have to comply with. They must fly along predetermined routes through certain waypoints, which are set by the agencies of air traffic control (ATC), something that usually fails to produce optimal results. Aircrafts are not allowed to fly directly to their final destination taking advantage of favorable winds without making changes to their trajectories causing unnecessary fuel costs, which can indirectly cause increases in ticket prices. This problem is particularly evident in transoceanic routes, which are experiencing the greatest growth in demand.
(2)
Increased ATC workload: AT controllers have, among other functions, to prevent collisions between aircrafts and redirect routes to avoid adverse conditions. In congested areas, such as regions near to airports calledterminal radar approach controls (TRACONs), AT controllers often simplify their high workloads making aircraft maintain default routes outside these regions, causing delays in landings and takeoffs.
(3)
Slow communication: communication is restricted to a tedious voice communication between aircraft and AT controllers, causing frequent bottleneck situations.In view of the problems described, the aviation community has been working in recent years on a concept calledfree flight. This innovative concept allows pilots to choose their own routes, altitude, and velocity to reduce delays and manage the use of aircraft fuel more efficiently. The preferences of the pilots will be restricted only in very congested airspace areas or to prevent unauthorized entry into military areas.Free flight is potentially possible due to the availability of technologies such asglobal positioning systems (GPS); communications data links such asautomatic dependence surveillance-broadcast (ADS-B); detection systems and collision avoidance; and powerful computing increasingly being implemented in aircrafts.In addition, there are several decision support tools that reduce the workload of both AT controllers and pilots and optimize their capacity, for example, the detection and resolution of conflicts in airspace sectors, landings and takeoffs management in airports, and organizational systems of the workload of AT controllers to better organize their tasks to increase productivity.These technological advances will also allow current ATC systems to accommodate the future growth of air traffic. Algorithms to detect and solve aircraft conflicts are vital to improve the safety standards in free flight unstructured environments. These systems can be used on land by ATC or by theflight management system (FMS) of each aircraft.In this paper, we focus on the development of algorithmic tools for aircraftconflict detection and resolution (CDR) problem. We assume that each aircraft is surrounded by cylinder representing a security virtual volume. Conflict between two aircrafts occurs when the respective aircraft security volumes overlap.Different approaches can be found in literature to deal with collision avoidance accounting for different number and types of admissible maneuvers for aircrafts and with different solution approaches, including the use of exact solvers, simulation techniques, and metaheuristics. The work in [1] present a survey with the most important of these up to the year 2000, whereas [2] focuses on approaches from 2000 up to 2012.One of the first approaches to deal with collision avoidance was [3]. A path planning problem among given waypoints avoiding all possible conflicts was considered aimed at minimizing the total flight time. Two mixed-integer linear programs were proposed accounting for velocity changes and angle changes as admissible maneuvers, respectively. The work in [4] proposes a three-dimensional formulation as a mixed-integer nonlinear program in which only velocity changes were admissible. CPLEX was used for the resolution in both approaches.Simulation techniques have also been used to handle CDR problems. For instance, [5] analyzes the economic performance of a specific conflict resolution strategy based on velocity change between two aircrafts in terms of extra time and fuel consumption. The work in [6] also considers a velocity regulation problem, but from a different perspective, distinguishing betweencrossing conflicts (the wider), in which the aircrafts intersect at some point and security cylinders overlap, andconflicts trail, caused when an aircraft pursues another, both with different velocities.Neural networks have been also used for performing velocity changes in CDR problems [5, 7, 8].More recently, [9] focuses on mixed-integer optimization models based on velocity regulation. They propose to accelerate or decelerate during a specified time interval, reverting back to the original velocity once the conflict is avoided. They propose a heuristic procedure where the problem is decomposed and locally exactly solved.Other less frequent proposals consider turn changes that lead to nonlinear optimization models. For instance, [10] proposes a two-step approach. First, a nonconvex mixed-integer nonlinear optimization is used to minimize the weighted aircraft angle variations. Then, a set of unconstrained quadratic optimization models are considered, where aircrafts are forced to return to their original flight plan as soon as possible once there is no aircraft in conflict with any other. Both an exact and an approximate resolution are proposed. In the second, the turn changes are discretized to reduce the search space.Different metaheuristics have been proposed for solving CDR models accounting for turn changes, such asant colony systems [11, 12],genetic algorithms [13],variable neighborhood search [14], andparticle swarm optimization [15], which uses a series of waypoints the aircrafts can pass through.A pretty realistic proposal is described in [16], wherein the acceleration variable is added to the model. It intends to solve conflicts discretizing the time remaining until it occurs at different intervals, optimizing acceleration, and velocity that should be assigned to each aircraft. A nonlinear mixed 0-1 model is used to solve the problem, which is iteratively linearized by using Taylor polynomials. This approach is then enhanced in [17], extending control to aircraft outside the aviation sector to manage, that is, taking into account those aircrafts leaving it or entering it. Moreover, they take into account the conflicts that may arise when an aircraft is climbing or descending to change altitude.The work in [18] improves the velocity change model by adding altitude changes when necessary, for example, in head-to-head conflict situations. A multiobjective perspective is considered including objectives such as velocity variation and total number of maneuvers and forcing to return to the original flight configuration when no aircrafts are in conflict. An exactly solved mixed 0-1 linear optimization model is used, with small computational time for the execution making it suitable for real-time use.In [19] an innovative point of view based on the choice of different strategies to avoid conflicts is proposed. An original trajectory model using B-splines is introduced together with a new semi-infinite programming formulation of the separation constraint involved in CDR problems.In this paper we propose usingsimulated annealing to deal with a CDR problem accounting for three admissible maneuvers (velocity, turn, and altitude changes) in a multiobjective context. Specifically, thearchive simulated annealing-based multiobjective optimization algorithm (AMOSA) has been adapted to the CDR problem accounting for objectives such as minimizing the maneuver number and magnitude, time delays, or deviations in the leaving points.Both the possibility of performing three types of maneuvers and the multiobjective context make this paper an original contribution regarding previous works on CDR problems.The paper is structured as follows. The mathematical modeling for the multiobjective problem under consideration is introduced in Section2, including the identification of parameters, decision variables, and constraints and the description and modelization of the candidate objective functions for analysis. AMOSA and its adaptation to the considered CDR problem are described in Section 3. Section 4 deals with the parameter setting and the performance analysis when 5 aircrafts are considered. The parameter setting for 10 and 20 aircrafts and an example illustrating the flexibility of the proposed algorithm are provided in Section 5. Finally, some conclusions are provided in Section 6.


## 2. Mathematical Modeling

We assume that in a particular moment there aren aircrafts in anaerial sector, a cubic volume in the space managed by an AT controller; see Figure 1. We have to decide which maneuvers to perform to avoid possible collisions between them. We also assume that the decision on maneuvers will be effective until n aircrafts leave the aerial sector or until a new aircraft enters the aerial sector. At the moment a new aircraft enters the aerial sector, the analysis we propose should be again carried out to identify new maneuvers to avoid new possible collisions caused by the entering aircraft.Figure 1
Aerial sector with several aircrafts.Moreover, we assume that the number of maneuvers performed by the aircrafts from the moment they entered the aerial sector until the moment the analysis is carried out is known. This information will be useful when analyzing the dispersion of maneuvers objective, as described afterwards.The collision avoidance problem can be mathematically modeled as follows.A conflict between two aircrafts occurs if the horizontal and/or vertical distances between them are smaller than some given security limits, making the security volume adopt a cylindrical form; see Figure2. Thus, the security cylinders corresponding to any couple of aircrafts should never intersect to avoid conflicts.Figure 2
Security volume of an aircraft.In the approach we propose, we account for three types of aircraft maneuvers:velocity change (V
C),altitude change (A
C), andturn change (T
C). A best type of maneuver does not exist since each maneuver has advantages and drawbacks. For example, the T
C maneuver is efficient regarding fuel consumption, but it makes the aircraft leave its original trajectory, which constitutes an inconvenience.Moreover, we consider a multiobjective perspective of the problem including the minimization of the maneuver number, magnitude and dispersion, time delays, and deviations in the leaving points.Thedispersion of the maneuvers objective aims at spreading the effort of all aircraft maneuvers in an attempt to avoid situations in which some aircrafts continuously perform maneuvers over time whereas other aircrafts do not do so. Thedeviations in the leaving points objective minimize the sum of the distances between the theoretical leaving points according to the initial trajectory when entering the aerial sector and the real leaving points after possible maneuver performances.


### 2.1. Parameters and Decision Variables

We considern aircrafts in an aerial sector at the time t. First, parameters accounting for thefeatures of the aircrafts are as follows:v
min
⁡
,
i: minimum velocity for the aircraft i.
v
max
⁡
,
i: maximum velocity for the aircraft i.
z
min
⁡
,
i: minimum altitude for the aircraft i.
z
max
⁡
,
i: maximum altitude for the aircraft i.
β
max
⁡
,
i: maximum variation of the angle for the aircraft i.Secondly, initial parameters for each aircraft when entering the aerial sector (denoted by subscriptini) are as follows:v
ini
,
i: velocity of the aircraft i upon entry.
z
ini
,
i: altitude of the aircraft i upon entry.
x
ini
,
i: abscise of the aircraft i upon entry.
y
ini
,
i: ordinate of the aircraft i upon entry.
α
ini
,
i: angle of the aircraft i with respect to the horizontal upon entry.Next, we consider final parameters for each aircraft when leaving the aerial sector assuming that no maneuvers have been performed, that is, according to its initial trajectory when entering the aerial sector, as follows:x
fin
,
i: estimated abscise of the aircraft i when leaving the aerial sector.
y
fin
,
i: estimated ordinate of the aircraft i when leaving the aerial sector.
z
fin
,
i: estimated altitude of the aircraft i when leaving the aerial sector.The parameters accounting for theconfiguration of the aircrafts at the time t are as follows:v
i
t: velocity of the aircraft i at the time t.
z
i
t: aircraft i altitude at the time t.
x
i
t: aircraft i abscise at the time t.
y
i
t: aircraft i ordinate at the time t.
t
i: time necessary to arrive to the bound of the aerial sector with a constant velocity v
i
t.
α
i
t: angle with respect to the horizontal of the aircraft i at the instant t.
m
a
n
i
t: number of maneuvers performed by the aircraft i at the time t since it entered in the aerial sector.Finally, parameters accounting forsecurity distances and the collision risk for aircrafts are as follows:s
hor
,
i: horizontal security distance for the aircraft i.
s
ver
,
i: vertical security distance for the aircraft i.
c
hor
,
ver: relative importance between vertical and horizontal risk.According to current standards, the horizontal security distance is usually 5 nautical miles, whereas the vertical security distance is 1000 feet (see Figure2), but other values could be used in the analysis.Regarding thedecision variables, as mentioned before, we propose three types of aircraft maneuvers:velocity change (VC),altitude change (AC
⁡), andturn change (TC). Thus, we consider three binary variables vc
i, ac
i, and t
c
i for each aircraft i, pointing out whether a velocity, altitude, or turn change is performed, respectively, and a continuous variable q
i representing the magnitude (proportion) of the change performed. Note that v
c
i
+
a
c
i
+
t
c
i
≤
1  ∀
i since each aircraft is allowed to perform, at most, one maneuver, and q
i
∈
[
-
1,1
] since the changes can be negative or positive.Consequently, a solution will consist of a vector with four elements per aircraft, identifying the maneuver performed and its magnitude.


### 2.2. Constraints

As a first approach, the main constraint should be avoiding conflicts between aircrafts; that is, the intersection between security cylinders should always be empty for any two aircrafts in the aerial sector under consideration. A geometric construction for detecting conflict situations is considered [3, 10]; see Figure 3.Figure 3
Detecting conflict situations.The velocity vectors of two aircraftsi and j are (1)
u
→
i
=
v
new
,
i
×
cos
⁡
α
new
,
i
,
v
new
,
i
×
sin
⁡
α
new
,
i
,
u
→
j
=
v
new
,
j
×
cos
⁡
α
new
,
j
,
v
new
,
j
×
sin
⁡
α
new
,
j
,respectively, where v
new
,
i is the new velocity considering the VC maneuver and α
new
,
i is the new angle with the horizontal plane by adding the magnitude of the TC maneuver.The basic idea of the model comes from the construction of the relative velocity vectoru
→
i
-
u
→
j; see Figure 3. The two straight lines parallel to the relative velocity vector and tangent to the security circle of aircraft j (dotted lines in Figure 3) define a region where the intersection with the trajectory for aircraft i is a segment named theshadow segment.Ahorizontal conflict occurs if the security cylinder of aircraft i intersects theshadow segment generated by aircraft j or, on the contrary, since u
→
i
-
u
→
j and u
→
j
-
u
→
i are parallels.Considering now the cutting planes that are tangent to both cylinders (see Figure4) and the angles g
i
j and l
i
j, there is no conflict if one of the following two conditions holds: (2)
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≥
tan
⁡
l
i
j
,
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≤
tan
⁡
g
i
j
.Figure 4
Main angles and distances.Note that the functions at the left of the above expressions can cause a zero denominator. These cases are referred to asmodel pathological situations and produce unstable solutions since a conflict between two aircrafts may be erroneously determined due to the null denominator, forcing the aircrafts to crash in the worse case. This situation is detected when |
x
i
-
x
j
|
<
s
hor
,
i
+
s
hor
,
j. Therefore, variables α
new
,
i, α
new
,
j and parameters l
i
j and g
i
j, which represent angles, are rotated π
/
2 radians when computing expressions for horizontal risk detection to overcome such pathological situation.Vertical conflicts are detected more easily considering the security cylinders. Computing the vertical distance between two aircrafts can detect these conflicts. Figure 5 shows the modeling of the vertical conflicts.Figure 5
Vertical security distance and turn vector.However, very restrictive situations may occur in an aerial sector at certain times. For example, it could be very difficult to find a feasible solution without (either horizontal or vertical) conflicts when there is a high density of nearby aircrafts. Therefore, we have relaxed the collision avoidance constraint, which becomes an additional objective function from now on, as described in the next section. This allows us to better explore the solution space but, in contrast, it complicates reaching the optimal solution.Finally, the following constraints check whether the new velocity, altitude, and turn satisfy the features of the corresponding aircraft:(3)
v
min
⁡
,
i
≤
v
new
,
i
≤
v
max
⁡
,
i
,
z
min
⁡
,
i
≤
z
new
,
i
≤
z
max
⁡
,
i
,
α
new
,
i
-
α
i
≤
β
max
⁡
,
i
,where v
new
,
i,  z
new
,
i, and α
new
,
i are the new aircraft configuration.


### 2.3. Objective Functions

Six different objective functions will be considered for analysis: specifically, minimizing the magnitude of maneuvers, collision risks, number of maneuvers, time delays, deviations in the leaving points, and the maneuver dispersion. Further information and the mathematical modeling of such objective functions are provided below.


#### 2.3.1. Objective 1: Minimizing Maneuver Magnitudes

It makes sense to claim that aircrafts perform maneuvers as smoothly as possible to avoid conflicts; that is, abrupt maneuvers that may disturb passengers or even be dangerous should be avoided. According to the above, the dispersion magnitude maneuvers performed by aircrafts should be incorporated into the analysis. Moreover, a high dispersion should be penalized, that is, situations in which maneuver magnitudes are very high for some aircrafts and very low for others.The first objective function,f
1, can be then modeled by the sum of the average maneuver magnitude and a dispersion term: (4)
m
i
n
f
1
=
a
1
+
w
1
n
∑
i
=
1
n
a
1
-
q
i
,where (5)
a
1
=
1
n
∑
i
=
1
n
q
i
,and w
1 represents the relative importance of the dispersion term regarding the average maneuver magnitude.


#### 2.3.2. Objective 2: Minimizing Collision Risks

This objective is the constraint we decided to relax. As cited before, this allows a better exploration of the solution space. Moreover, if we consider two solutions in which there are no conflicts we can compute a conflict risk measure for each of them accounting for the distances between all pairs of aircrafts.There is a conflict between a couple of aircrafts if their collision risk is positiver
i
j
>
0; otherwise (r
i
j
≤
0) there is no conflict.We differentiate the situation in which there is at least one conflict between a couple of aircrafts and the one in which there are no conflicts to compute an average collision risk (a
2). In the first case, a
2 is computed taking into account only those r
i
j
>
0, avoiding that negative values nullify the positive. In the second case (r
i
j
≤
0, ∀
i
,
j), all r
i
j are used to compute a
2, which is a negative value: (6)
a
2
=
1
k
1
∑
i
=
1
n
∑
j
=
i
+
1
n
r
i
j
,
i
f
r
i
j
≤
0
∀
i
,
j
1
k
2
∑
i
=
1
n
∑
j
=
i
+
1
n
φ
i
j
,
i
f
∃
a
t
l
e
a
s
t
o
n
e
r
i
j
>
0
,where (7)
φ
i
j
=
r
i
j
,
i
f
r
i
j
>
0
0
,
o
t
h
e
r
w
i
s
e
,and k
1 and k
2 are the numbers of elements considered in the average computations, that is, k
1
=
n
!
/
(
n
-
2
)
!
4, since all pairs of aircrafts are considered, whereas k
2 is the number of couples with conflict.The above two situations are also considered in the modelization of objective function,f
2. The average collision risk is minimized in both cases, but in the second case (a
2
≤
0), the dispersion on the collision risk values is also considered: (8)
m
i
n
f
2
=
a
2
,
i
f
a
2
>
0
a
2
+
w
2
k
1
∑
i
=
1
n
∑
j
=
i
+
1
n
a
2
-
r
i
j
,
e
l
s
e
.Next, we clarify how the collision risk for a couple of aircrafts is computed. The collision riskr
i
j between two aircrafts i and j is computed in four different ways, depending on the combinations of vertical and horizontal conflicts:(9)
r
i
j
=
r
hor
,
i
j
+
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
≤
0
∧
r
ver
,
i
j
≤
0
r
hor
,
i
j
+
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
>
0
∧
r
ver
,
i
j
>
0
r
hor
,
i
j
,
i
f
r
hor
,
i
j
≤
0
∧
r
ver
,
i
j
>
0
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
>
0
∧
r
ver
,
i
j
≤
0
.To computehorizontal collision risk, we take into account the analysis of the horizontal conflict shown in Section 2.2, which concludes that there is no horizontal conflict if one of the following two conditions holds: (10)
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≥
tan
⁡
l
i
j
,
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≤
tan
⁡
g
i
j
.Thehorizontal collision risk (r
hor
,
i
j) is computed as (11)
r
hor
,
i
j
=
min
⁡
γ
i
j
,
δ
i
j
,
γ
i
j
=
tan
⁡
l
i
j
-
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
,
δ
i
j
=
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
-
tan
⁡
g
i
j
.Thus, it is possible to assess how far two aircrafts are to horizontally invade their respective security cylinders or, in case of conflict, to measure the intensity of invasion.The new aircraft configuration changes or not depending on the maneuver being performed:(12)
v
new
,
i
=
v
i
,
i
f
v
c
i
=
0
v
i
+
q
i
v
i
-
v
min
⁡
,
i
,
i
f
v
c
i
=
1
∧
q
i
≤
0
v
i
+
q
i
v
max
⁡
,
i
-
v
i
,
i
f
v
c
i
=
1
∧
q
i
>
0
,
α
new
,
i
=
α
i
,
i
f
t
c
i
=
0
α
i
+
q
i
β
max
⁡
,
i
,
i
f
t
c
i
=
1
.Parametersl
i
j and g
i
j are computed on the basis of the angles and distances shown in Figure 4: (13)
l
i
j
=
ω
i
j
+
θ
i
j
,
g
i
j
=
ω
i
j
-
θ
i
j
,
ω
i
j
=
arctan
⁡
y
i
-
y
j
x
i
-
x
j
,
θ
i
j
=
arcsin
⁡
s
hor
,
i
+
s
hor
,
j
/
2
d
i
j
/
2
.Thevertical collision risk (r
ver
,
i
j) just consists of computing the vertical distance between aircrafts added to the vertical security distance, so that if there is no conflict, then r
ver
,
i
j
≤
0 and vice versa:(14)
r
ver
,
i
j
=
max
⁡
z
new
,
j
-
z
new
,
i
+
s
ver
,
j
,
z
new
,
j
-
z
new
,
i
+
s
ver
,
i
,
i
f
z
new
,
i
>
z
new
,
j
max
⁡
z
new
,
i
-
z
new
,
j
+
s
ver
,
i
,
z
new
,
i
-
z
new
,
j
+
s
ver
,
j
,
i
f
z
new
,
i
≤
z
new
,
j
,where z
new
,
i takes different values depending on whether or not this maneuver is performed and on whether the altitude is increased or decreased:(15)
z
new
,
i
=
0
,
i
f
a
c
i
=
0
z
i
+
q
i
z
i
-
z
min
⁡
,
i
,
i
f
a
c
i
=
1
∧
q
i
≤
0
z
i
+
q
i
z
max
⁡
,
i
-
z
i
,
i
f
a
c
i
=
1
∧
q
i
>
0
.


#### 2.3.3. Objective 3: Minimizing the Numbers of Maneuvers

The objective is to minimize the number of maneuvers performed byn aircrafts. This is directly connected with the AT controllers workload, since they communicate the corresponding maneuvers to the pilots. Therefore, this goal is equivalent to minimizing the controllers’ workload. To model this objective function we just have to sum up the three binary variables associated to possible maneuvers of the aircrafts in the solution under consideration(16)
m
i
n
f
3
=
∑
i
=
1
n
v
c
i
+
a
c
i
+
t
c
i
.Note again that v
c
i
+
a
c
i
+
t
c
i
≤
1  ∀
i since each aircraft can perform at most one type of maneuver.


#### 2.3.4. Objective 4: Minimizing Time Delays

The time an aircraft will leave the aerial sector may differ from the expected time according to its initial trajectory (when entering the aerial sector) if aVC or TC maneuver is performed. Then, the aim is now minimizing the sum of the delays for the aircrafts. Moreover, we will again take into account the dispersion of such delays.This objective function,f
4, is computed by adding the average delay (a
4) and the dispersion of the delays: (17)
m
i
n
f
4
=
a
4
+
w
4
n
∑
i
=
1
n
a
4
-
t
new
,
i
-
t
i
,where a
4
=
(
1
/
n
)
∑
i
=
1
n
|
t
new
,
i
-
t
i
| and (18)
t
new
,
i
=
x
fin
,
i
-
x
i
2
+
y
fin
,
i
-
y
i
2
v
new
,
i
,is the time aircraft i leaves the aerial sector once the maneuvers associated to the solution under consideration are performed. Note that t
new
,
i
=
t
i if an AC
⁡ maneuver is performed.Inf
4, w
4
∈
[
0,1
] represents the relative importance of the dispersion term regarding the average delay, analogously to w
1 and w
2 in objective functions f
1 and f
2, respectively.


#### 2.3.5. Objective 5: Minimizing Deviations in the Leaving Points

The point at which an aircraft will leave the aerial sector may differ from the expected point according to its initial trajectory (when entering the aerial sector) as a consequence of the maneuvers performed. As such, the aim is now minimizing the sum of the distances between both leaving points for each aircraft. Moreover, we will again take into account the dispersion of such distances.This objective function,f
5, is computed by adding the average (a
5) and the dispersion of the distances between the original leaving point and the new leaving point once the maneuvers associated to the solution under consideration are performed: (19)
m
i
n
f
5
=
a
5
+
w
5
n
∑
i
=
1
n
a
5
-
d
i
w
i
t
h
a
5
=
1
n
∑
i
=
1
n
d
i
,where the distance is computed as follows: (20)
d
i
=
x
fin
,
i
-
x
fin
,
new
,
i
2
+
y
fin
,
i
-
y
fin
,
new
,
i
2
+
z
fin
,
i
-
z
fin
,
new
,
i
2
,where x
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
x
i
,
α
new
,
i,  y
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
y
i
,
α
new
,
i,  and z
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
z
i
,
α
new
,
i are the coordinates of the new leaving point, with α
new
,
i being the new angle with respect to the horizontal of the aircraft i once the maneuvers associated to the solution under consideration are performed.w
5
∈
[
0,1
] again represents the relative importance of the dispersion term regarding the average value (a
5).


#### 2.3.6. Objective 6: Minimizing Maneuver Dispersion

The last objective under consideration is related to the dispersion of maneuvers over time, that is, over multiple executions of the analysis. The aim is to share the maneuver effort among the aircrafts to attempt to avoid situations in which some aircrafts continuously perform maneuvers over time, whereas other aircrafts scarcely do it.For this, as mentioned before, we assume that a vector including the number of maneuvers performed by each aircraft from the moments they entered the aerial sector until the moment the analysis is carried out is available (m
a
n
1
,
…
,
m
a
n
i
,
…
,
m
a
n
n). Thus, a dispersion measure can be computed as follows: (21)
m
i
n
f
6
=
1
n
∑
i
=
1
n
a
6
-
m
a
n
n
e
w
i
,with a
6
=
(
∑
i
=
1
n
‍
m
a
n
new
,
i
)
/
n being the average number of maneuvers and (22)
m
a
n
new
,
i
=
m
a
n
i
,
i
f
v
c
i
=
0
∧
a
c
i
=
0
∧
t
c
i
=
0
m
a
n
i
+
1
,
i
f
v
c
i
=
1
∨
a
c
i
=
1
∨
t
c
i
=
1
,being the new number of maneuvers accumulated by aircraft i if the maneuvers associated to the solution under consideration are performed.


## 2.1. Parameters and Decision Variables

We considern aircrafts in an aerial sector at the time t. First, parameters accounting for thefeatures of the aircrafts are as follows:v
min
⁡
,
i: minimum velocity for the aircraft i.
v
max
⁡
,
i: maximum velocity for the aircraft i.
z
min
⁡
,
i: minimum altitude for the aircraft i.
z
max
⁡
,
i: maximum altitude for the aircraft i.
β
max
⁡
,
i: maximum variation of the angle for the aircraft i.Secondly, initial parameters for each aircraft when entering the aerial sector (denoted by subscriptini) are as follows:v
ini
,
i: velocity of the aircraft i upon entry.
z
ini
,
i: altitude of the aircraft i upon entry.
x
ini
,
i: abscise of the aircraft i upon entry.
y
ini
,
i: ordinate of the aircraft i upon entry.
α
ini
,
i: angle of the aircraft i with respect to the horizontal upon entry.Next, we consider final parameters for each aircraft when leaving the aerial sector assuming that no maneuvers have been performed, that is, according to its initial trajectory when entering the aerial sector, as follows:x
fin
,
i: estimated abscise of the aircraft i when leaving the aerial sector.
y
fin
,
i: estimated ordinate of the aircraft i when leaving the aerial sector.
z
fin
,
i: estimated altitude of the aircraft i when leaving the aerial sector.The parameters accounting for theconfiguration of the aircrafts at the time t are as follows:v
i
t: velocity of the aircraft i at the time t.
z
i
t: aircraft i altitude at the time t.
x
i
t: aircraft i abscise at the time t.
y
i
t: aircraft i ordinate at the time t.
t
i: time necessary to arrive to the bound of the aerial sector with a constant velocity v
i
t.
α
i
t: angle with respect to the horizontal of the aircraft i at the instant t.
m
a
n
i
t: number of maneuvers performed by the aircraft i at the time t since it entered in the aerial sector.Finally, parameters accounting forsecurity distances and the collision risk for aircrafts are as follows:s
hor
,
i: horizontal security distance for the aircraft i.
s
ver
,
i: vertical security distance for the aircraft i.
c
hor
,
ver: relative importance between vertical and horizontal risk.According to current standards, the horizontal security distance is usually 5 nautical miles, whereas the vertical security distance is 1000 feet (see Figure2), but other values could be used in the analysis.Regarding thedecision variables, as mentioned before, we propose three types of aircraft maneuvers:velocity change (VC),altitude change (AC
⁡), andturn change (TC). Thus, we consider three binary variables vc
i, ac
i, and t
c
i for each aircraft i, pointing out whether a velocity, altitude, or turn change is performed, respectively, and a continuous variable q
i representing the magnitude (proportion) of the change performed. Note that v
c
i
+
a
c
i
+
t
c
i
≤
1  ∀
i since each aircraft is allowed to perform, at most, one maneuver, and q
i
∈
[
-
1,1
] since the changes can be negative or positive.Consequently, a solution will consist of a vector with four elements per aircraft, identifying the maneuver performed and its magnitude.


## 2.2. Constraints

As a first approach, the main constraint should be avoiding conflicts between aircrafts; that is, the intersection between security cylinders should always be empty for any two aircrafts in the aerial sector under consideration. A geometric construction for detecting conflict situations is considered [3, 10]; see Figure 3.Figure 3
Detecting conflict situations.The velocity vectors of two aircraftsi and j are (1)
u
→
i
=
v
new
,
i
×
cos
⁡
α
new
,
i
,
v
new
,
i
×
sin
⁡
α
new
,
i
,
u
→
j
=
v
new
,
j
×
cos
⁡
α
new
,
j
,
v
new
,
j
×
sin
⁡
α
new
,
j
,respectively, where v
new
,
i is the new velocity considering the VC maneuver and α
new
,
i is the new angle with the horizontal plane by adding the magnitude of the TC maneuver.The basic idea of the model comes from the construction of the relative velocity vectoru
→
i
-
u
→
j; see Figure 3. The two straight lines parallel to the relative velocity vector and tangent to the security circle of aircraft j (dotted lines in Figure 3) define a region where the intersection with the trajectory for aircraft i is a segment named theshadow segment.Ahorizontal conflict occurs if the security cylinder of aircraft i intersects theshadow segment generated by aircraft j or, on the contrary, since u
→
i
-
u
→
j and u
→
j
-
u
→
i are parallels.Considering now the cutting planes that are tangent to both cylinders (see Figure4) and the angles g
i
j and l
i
j, there is no conflict if one of the following two conditions holds: (2)
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≥
tan
⁡
l
i
j
,
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≤
tan
⁡
g
i
j
.Figure 4
Main angles and distances.Note that the functions at the left of the above expressions can cause a zero denominator. These cases are referred to asmodel pathological situations and produce unstable solutions since a conflict between two aircrafts may be erroneously determined due to the null denominator, forcing the aircrafts to crash in the worse case. This situation is detected when |
x
i
-
x
j
|
<
s
hor
,
i
+
s
hor
,
j. Therefore, variables α
new
,
i, α
new
,
j and parameters l
i
j and g
i
j, which represent angles, are rotated π
/
2 radians when computing expressions for horizontal risk detection to overcome such pathological situation.Vertical conflicts are detected more easily considering the security cylinders. Computing the vertical distance between two aircrafts can detect these conflicts. Figure 5 shows the modeling of the vertical conflicts.Figure 5
Vertical security distance and turn vector.However, very restrictive situations may occur in an aerial sector at certain times. For example, it could be very difficult to find a feasible solution without (either horizontal or vertical) conflicts when there is a high density of nearby aircrafts. Therefore, we have relaxed the collision avoidance constraint, which becomes an additional objective function from now on, as described in the next section. This allows us to better explore the solution space but, in contrast, it complicates reaching the optimal solution.Finally, the following constraints check whether the new velocity, altitude, and turn satisfy the features of the corresponding aircraft:(3)
v
min
⁡
,
i
≤
v
new
,
i
≤
v
max
⁡
,
i
,
z
min
⁡
,
i
≤
z
new
,
i
≤
z
max
⁡
,
i
,
α
new
,
i
-
α
i
≤
β
max
⁡
,
i
,where v
new
,
i,  z
new
,
i, and α
new
,
i are the new aircraft configuration.


## 2.3. Objective Functions

Six different objective functions will be considered for analysis: specifically, minimizing the magnitude of maneuvers, collision risks, number of maneuvers, time delays, deviations in the leaving points, and the maneuver dispersion. Further information and the mathematical modeling of such objective functions are provided below.


### 2.3.1. Objective 1: Minimizing Maneuver Magnitudes

It makes sense to claim that aircrafts perform maneuvers as smoothly as possible to avoid conflicts; that is, abrupt maneuvers that may disturb passengers or even be dangerous should be avoided. According to the above, the dispersion magnitude maneuvers performed by aircrafts should be incorporated into the analysis. Moreover, a high dispersion should be penalized, that is, situations in which maneuver magnitudes are very high for some aircrafts and very low for others.The first objective function,f
1, can be then modeled by the sum of the average maneuver magnitude and a dispersion term: (4)
m
i
n
f
1
=
a
1
+
w
1
n
∑
i
=
1
n
a
1
-
q
i
,where (5)
a
1
=
1
n
∑
i
=
1
n
q
i
,and w
1 represents the relative importance of the dispersion term regarding the average maneuver magnitude.


### 2.3.2. Objective 2: Minimizing Collision Risks

This objective is the constraint we decided to relax. As cited before, this allows a better exploration of the solution space. Moreover, if we consider two solutions in which there are no conflicts we can compute a conflict risk measure for each of them accounting for the distances between all pairs of aircrafts.There is a conflict between a couple of aircrafts if their collision risk is positiver
i
j
>
0; otherwise (r
i
j
≤
0) there is no conflict.We differentiate the situation in which there is at least one conflict between a couple of aircrafts and the one in which there are no conflicts to compute an average collision risk (a
2). In the first case, a
2 is computed taking into account only those r
i
j
>
0, avoiding that negative values nullify the positive. In the second case (r
i
j
≤
0, ∀
i
,
j), all r
i
j are used to compute a
2, which is a negative value: (6)
a
2
=
1
k
1
∑
i
=
1
n
∑
j
=
i
+
1
n
r
i
j
,
i
f
r
i
j
≤
0
∀
i
,
j
1
k
2
∑
i
=
1
n
∑
j
=
i
+
1
n
φ
i
j
,
i
f
∃
a
t
l
e
a
s
t
o
n
e
r
i
j
>
0
,where (7)
φ
i
j
=
r
i
j
,
i
f
r
i
j
>
0
0
,
o
t
h
e
r
w
i
s
e
,and k
1 and k
2 are the numbers of elements considered in the average computations, that is, k
1
=
n
!
/
(
n
-
2
)
!
4, since all pairs of aircrafts are considered, whereas k
2 is the number of couples with conflict.The above two situations are also considered in the modelization of objective function,f
2. The average collision risk is minimized in both cases, but in the second case (a
2
≤
0), the dispersion on the collision risk values is also considered: (8)
m
i
n
f
2
=
a
2
,
i
f
a
2
>
0
a
2
+
w
2
k
1
∑
i
=
1
n
∑
j
=
i
+
1
n
a
2
-
r
i
j
,
e
l
s
e
.Next, we clarify how the collision risk for a couple of aircrafts is computed. The collision riskr
i
j between two aircrafts i and j is computed in four different ways, depending on the combinations of vertical and horizontal conflicts:(9)
r
i
j
=
r
hor
,
i
j
+
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
≤
0
∧
r
ver
,
i
j
≤
0
r
hor
,
i
j
+
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
>
0
∧
r
ver
,
i
j
>
0
r
hor
,
i
j
,
i
f
r
hor
,
i
j
≤
0
∧
r
ver
,
i
j
>
0
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
>
0
∧
r
ver
,
i
j
≤
0
.To computehorizontal collision risk, we take into account the analysis of the horizontal conflict shown in Section 2.2, which concludes that there is no horizontal conflict if one of the following two conditions holds: (10)
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≥
tan
⁡
l
i
j
,
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≤
tan
⁡
g
i
j
.Thehorizontal collision risk (r
hor
,
i
j) is computed as (11)
r
hor
,
i
j
=
min
⁡
γ
i
j
,
δ
i
j
,
γ
i
j
=
tan
⁡
l
i
j
-
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
,
δ
i
j
=
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
-
tan
⁡
g
i
j
.Thus, it is possible to assess how far two aircrafts are to horizontally invade their respective security cylinders or, in case of conflict, to measure the intensity of invasion.The new aircraft configuration changes or not depending on the maneuver being performed:(12)
v
new
,
i
=
v
i
,
i
f
v
c
i
=
0
v
i
+
q
i
v
i
-
v
min
⁡
,
i
,
i
f
v
c
i
=
1
∧
q
i
≤
0
v
i
+
q
i
v
max
⁡
,
i
-
v
i
,
i
f
v
c
i
=
1
∧
q
i
>
0
,
α
new
,
i
=
α
i
,
i
f
t
c
i
=
0
α
i
+
q
i
β
max
⁡
,
i
,
i
f
t
c
i
=
1
.Parametersl
i
j and g
i
j are computed on the basis of the angles and distances shown in Figure 4: (13)
l
i
j
=
ω
i
j
+
θ
i
j
,
g
i
j
=
ω
i
j
-
θ
i
j
,
ω
i
j
=
arctan
⁡
y
i
-
y
j
x
i
-
x
j
,
θ
i
j
=
arcsin
⁡
s
hor
,
i
+
s
hor
,
j
/
2
d
i
j
/
2
.Thevertical collision risk (r
ver
,
i
j) just consists of computing the vertical distance between aircrafts added to the vertical security distance, so that if there is no conflict, then r
ver
,
i
j
≤
0 and vice versa:(14)
r
ver
,
i
j
=
max
⁡
z
new
,
j
-
z
new
,
i
+
s
ver
,
j
,
z
new
,
j
-
z
new
,
i
+
s
ver
,
i
,
i
f
z
new
,
i
>
z
new
,
j
max
⁡
z
new
,
i
-
z
new
,
j
+
s
ver
,
i
,
z
new
,
i
-
z
new
,
j
+
s
ver
,
j
,
i
f
z
new
,
i
≤
z
new
,
j
,where z
new
,
i takes different values depending on whether or not this maneuver is performed and on whether the altitude is increased or decreased:(15)
z
new
,
i
=
0
,
i
f
a
c
i
=
0
z
i
+
q
i
z
i
-
z
min
⁡
,
i
,
i
f
a
c
i
=
1
∧
q
i
≤
0
z
i
+
q
i
z
max
⁡
,
i
-
z
i
,
i
f
a
c
i
=
1
∧
q
i
>
0
.


### 2.3.3. Objective 3: Minimizing the Numbers of Maneuvers

The objective is to minimize the number of maneuvers performed byn aircrafts. This is directly connected with the AT controllers workload, since they communicate the corresponding maneuvers to the pilots. Therefore, this goal is equivalent to minimizing the controllers’ workload. To model this objective function we just have to sum up the three binary variables associated to possible maneuvers of the aircrafts in the solution under consideration(16)
m
i
n
f
3
=
∑
i
=
1
n
v
c
i
+
a
c
i
+
t
c
i
.Note again that v
c
i
+
a
c
i
+
t
c
i
≤
1  ∀
i since each aircraft can perform at most one type of maneuver.


### 2.3.4. Objective 4: Minimizing Time Delays

The time an aircraft will leave the aerial sector may differ from the expected time according to its initial trajectory (when entering the aerial sector) if aVC or TC maneuver is performed. Then, the aim is now minimizing the sum of the delays for the aircrafts. Moreover, we will again take into account the dispersion of such delays.This objective function,f
4, is computed by adding the average delay (a
4) and the dispersion of the delays: (17)
m
i
n
f
4
=
a
4
+
w
4
n
∑
i
=
1
n
a
4
-
t
new
,
i
-
t
i
,where a
4
=
(
1
/
n
)
∑
i
=
1
n
|
t
new
,
i
-
t
i
| and (18)
t
new
,
i
=
x
fin
,
i
-
x
i
2
+
y
fin
,
i
-
y
i
2
v
new
,
i
,is the time aircraft i leaves the aerial sector once the maneuvers associated to the solution under consideration are performed. Note that t
new
,
i
=
t
i if an AC
⁡ maneuver is performed.Inf
4, w
4
∈
[
0,1
] represents the relative importance of the dispersion term regarding the average delay, analogously to w
1 and w
2 in objective functions f
1 and f
2, respectively.


### 2.3.5. Objective 5: Minimizing Deviations in the Leaving Points

The point at which an aircraft will leave the aerial sector may differ from the expected point according to its initial trajectory (when entering the aerial sector) as a consequence of the maneuvers performed. As such, the aim is now minimizing the sum of the distances between both leaving points for each aircraft. Moreover, we will again take into account the dispersion of such distances.This objective function,f
5, is computed by adding the average (a
5) and the dispersion of the distances between the original leaving point and the new leaving point once the maneuvers associated to the solution under consideration are performed: (19)
m
i
n
f
5
=
a
5
+
w
5
n
∑
i
=
1
n
a
5
-
d
i
w
i
t
h
a
5
=
1
n
∑
i
=
1
n
d
i
,where the distance is computed as follows: (20)
d
i
=
x
fin
,
i
-
x
fin
,
new
,
i
2
+
y
fin
,
i
-
y
fin
,
new
,
i
2
+
z
fin
,
i
-
z
fin
,
new
,
i
2
,where x
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
x
i
,
α
new
,
i,  y
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
y
i
,
α
new
,
i,  and z
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
z
i
,
α
new
,
i are the coordinates of the new leaving point, with α
new
,
i being the new angle with respect to the horizontal of the aircraft i once the maneuvers associated to the solution under consideration are performed.w
5
∈
[
0,1
] again represents the relative importance of the dispersion term regarding the average value (a
5).


### 2.3.6. Objective 6: Minimizing Maneuver Dispersion

The last objective under consideration is related to the dispersion of maneuvers over time, that is, over multiple executions of the analysis. The aim is to share the maneuver effort among the aircrafts to attempt to avoid situations in which some aircrafts continuously perform maneuvers over time, whereas other aircrafts scarcely do it.For this, as mentioned before, we assume that a vector including the number of maneuvers performed by each aircraft from the moments they entered the aerial sector until the moment the analysis is carried out is available (m
a
n
1
,
…
,
m
a
n
i
,
…
,
m
a
n
n). Thus, a dispersion measure can be computed as follows: (21)
m
i
n
f
6
=
1
n
∑
i
=
1
n
a
6
-
m
a
n
n
e
w
i
,with a
6
=
(
∑
i
=
1
n
‍
m
a
n
new
,
i
)
/
n being the average number of maneuvers and (22)
m
a
n
new
,
i
=
m
a
n
i
,
i
f
v
c
i
=
0
∧
a
c
i
=
0
∧
t
c
i
=
0
m
a
n
i
+
1
,
i
f
v
c
i
=
1
∨
a
c
i
=
1
∨
t
c
i
=
1
,being the new number of maneuvers accumulated by aircraft i if the maneuvers associated to the solution under consideration are performed.


## 2.3.1. Objective 1: Minimizing Maneuver Magnitudes

It makes sense to claim that aircrafts perform maneuvers as smoothly as possible to avoid conflicts; that is, abrupt maneuvers that may disturb passengers or even be dangerous should be avoided. According to the above, the dispersion magnitude maneuvers performed by aircrafts should be incorporated into the analysis. Moreover, a high dispersion should be penalized, that is, situations in which maneuver magnitudes are very high for some aircrafts and very low for others.The first objective function,f
1, can be then modeled by the sum of the average maneuver magnitude and a dispersion term: (4)
m
i
n
f
1
=
a
1
+
w
1
n
∑
i
=
1
n
a
1
-
q
i
,where (5)
a
1
=
1
n
∑
i
=
1
n
q
i
,and w
1 represents the relative importance of the dispersion term regarding the average maneuver magnitude.


## 2.3.2. Objective 2: Minimizing Collision Risks

This objective is the constraint we decided to relax. As cited before, this allows a better exploration of the solution space. Moreover, if we consider two solutions in which there are no conflicts we can compute a conflict risk measure for each of them accounting for the distances between all pairs of aircrafts.There is a conflict between a couple of aircrafts if their collision risk is positiver
i
j
>
0; otherwise (r
i
j
≤
0) there is no conflict.We differentiate the situation in which there is at least one conflict between a couple of aircrafts and the one in which there are no conflicts to compute an average collision risk (a
2). In the first case, a
2 is computed taking into account only those r
i
j
>
0, avoiding that negative values nullify the positive. In the second case (r
i
j
≤
0, ∀
i
,
j), all r
i
j are used to compute a
2, which is a negative value: (6)
a
2
=
1
k
1
∑
i
=
1
n
∑
j
=
i
+
1
n
r
i
j
,
i
f
r
i
j
≤
0
∀
i
,
j
1
k
2
∑
i
=
1
n
∑
j
=
i
+
1
n
φ
i
j
,
i
f
∃
a
t
l
e
a
s
t
o
n
e
r
i
j
>
0
,where (7)
φ
i
j
=
r
i
j
,
i
f
r
i
j
>
0
0
,
o
t
h
e
r
w
i
s
e
,and k
1 and k
2 are the numbers of elements considered in the average computations, that is, k
1
=
n
!
/
(
n
-
2
)
!
4, since all pairs of aircrafts are considered, whereas k
2 is the number of couples with conflict.The above two situations are also considered in the modelization of objective function,f
2. The average collision risk is minimized in both cases, but in the second case (a
2
≤
0), the dispersion on the collision risk values is also considered: (8)
m
i
n
f
2
=
a
2
,
i
f
a
2
>
0
a
2
+
w
2
k
1
∑
i
=
1
n
∑
j
=
i
+
1
n
a
2
-
r
i
j
,
e
l
s
e
.Next, we clarify how the collision risk for a couple of aircrafts is computed. The collision riskr
i
j between two aircrafts i and j is computed in four different ways, depending on the combinations of vertical and horizontal conflicts:(9)
r
i
j
=
r
hor
,
i
j
+
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
≤
0
∧
r
ver
,
i
j
≤
0
r
hor
,
i
j
+
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
>
0
∧
r
ver
,
i
j
>
0
r
hor
,
i
j
,
i
f
r
hor
,
i
j
≤
0
∧
r
ver
,
i
j
>
0
c
hor
,
ver
×
r
ver
,
i
j
,
i
f
r
hor
,
i
j
>
0
∧
r
ver
,
i
j
≤
0
.To computehorizontal collision risk, we take into account the analysis of the horizontal conflict shown in Section 2.2, which concludes that there is no horizontal conflict if one of the following two conditions holds: (10)
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≥
tan
⁡
l
i
j
,
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
≤
tan
⁡
g
i
j
.Thehorizontal collision risk (r
hor
,
i
j) is computed as (11)
r
hor
,
i
j
=
min
⁡
γ
i
j
,
δ
i
j
,
γ
i
j
=
tan
⁡
l
i
j
-
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
,
δ
i
j
=
v
new
,
i
×
sin
⁡
α
new
,
i
-
v
new
,
j
×
sin
⁡
α
new
,
j
v
new
,
i
×
cos
⁡
α
new
,
i
-
v
new
,
j
×
cos
⁡
α
new
,
j
-
tan
⁡
g
i
j
.Thus, it is possible to assess how far two aircrafts are to horizontally invade their respective security cylinders or, in case of conflict, to measure the intensity of invasion.The new aircraft configuration changes or not depending on the maneuver being performed:(12)
v
new
,
i
=
v
i
,
i
f
v
c
i
=
0
v
i
+
q
i
v
i
-
v
min
⁡
,
i
,
i
f
v
c
i
=
1
∧
q
i
≤
0
v
i
+
q
i
v
max
⁡
,
i
-
v
i
,
i
f
v
c
i
=
1
∧
q
i
>
0
,
α
new
,
i
=
α
i
,
i
f
t
c
i
=
0
α
i
+
q
i
β
max
⁡
,
i
,
i
f
t
c
i
=
1
.Parametersl
i
j and g
i
j are computed on the basis of the angles and distances shown in Figure 4: (13)
l
i
j
=
ω
i
j
+
θ
i
j
,
g
i
j
=
ω
i
j
-
θ
i
j
,
ω
i
j
=
arctan
⁡
y
i
-
y
j
x
i
-
x
j
,
θ
i
j
=
arcsin
⁡
s
hor
,
i
+
s
hor
,
j
/
2
d
i
j
/
2
.Thevertical collision risk (r
ver
,
i
j) just consists of computing the vertical distance between aircrafts added to the vertical security distance, so that if there is no conflict, then r
ver
,
i
j
≤
0 and vice versa:(14)
r
ver
,
i
j
=
max
⁡
z
new
,
j
-
z
new
,
i
+
s
ver
,
j
,
z
new
,
j
-
z
new
,
i
+
s
ver
,
i
,
i
f
z
new
,
i
>
z
new
,
j
max
⁡
z
new
,
i
-
z
new
,
j
+
s
ver
,
i
,
z
new
,
i
-
z
new
,
j
+
s
ver
,
j
,
i
f
z
new
,
i
≤
z
new
,
j
,where z
new
,
i takes different values depending on whether or not this maneuver is performed and on whether the altitude is increased or decreased:(15)
z
new
,
i
=
0
,
i
f
a
c
i
=
0
z
i
+
q
i
z
i
-
z
min
⁡
,
i
,
i
f
a
c
i
=
1
∧
q
i
≤
0
z
i
+
q
i
z
max
⁡
,
i
-
z
i
,
i
f
a
c
i
=
1
∧
q
i
>
0
.


## 2.3.3. Objective 3: Minimizing the Numbers of Maneuvers

The objective is to minimize the number of maneuvers performed byn aircrafts. This is directly connected with the AT controllers workload, since they communicate the corresponding maneuvers to the pilots. Therefore, this goal is equivalent to minimizing the controllers’ workload. To model this objective function we just have to sum up the three binary variables associated to possible maneuvers of the aircrafts in the solution under consideration(16)
m
i
n
f
3
=
∑
i
=
1
n
v
c
i
+
a
c
i
+
t
c
i
.Note again that v
c
i
+
a
c
i
+
t
c
i
≤
1  ∀
i since each aircraft can perform at most one type of maneuver.


## 2.3.4. Objective 4: Minimizing Time Delays

The time an aircraft will leave the aerial sector may differ from the expected time according to its initial trajectory (when entering the aerial sector) if aVC or TC maneuver is performed. Then, the aim is now minimizing the sum of the delays for the aircrafts. Moreover, we will again take into account the dispersion of such delays.This objective function,f
4, is computed by adding the average delay (a
4) and the dispersion of the delays: (17)
m
i
n
f
4
=
a
4
+
w
4
n
∑
i
=
1
n
a
4
-
t
new
,
i
-
t
i
,where a
4
=
(
1
/
n
)
∑
i
=
1
n
|
t
new
,
i
-
t
i
| and (18)
t
new
,
i
=
x
fin
,
i
-
x
i
2
+
y
fin
,
i
-
y
i
2
v
new
,
i
,is the time aircraft i leaves the aerial sector once the maneuvers associated to the solution under consideration are performed. Note that t
new
,
i
=
t
i if an AC
⁡ maneuver is performed.Inf
4, w
4
∈
[
0,1
] represents the relative importance of the dispersion term regarding the average delay, analogously to w
1 and w
2 in objective functions f
1 and f
2, respectively.


## 2.3.5. Objective 5: Minimizing Deviations in the Leaving Points

The point at which an aircraft will leave the aerial sector may differ from the expected point according to its initial trajectory (when entering the aerial sector) as a consequence of the maneuvers performed. As such, the aim is now minimizing the sum of the distances between both leaving points for each aircraft. Moreover, we will again take into account the dispersion of such distances.This objective function,f
5, is computed by adding the average (a
5) and the dispersion of the distances between the original leaving point and the new leaving point once the maneuvers associated to the solution under consideration are performed: (19)
m
i
n
f
5
=
a
5
+
w
5
n
∑
i
=
1
n
a
5
-
d
i
w
i
t
h
a
5
=
1
n
∑
i
=
1
n
d
i
,where the distance is computed as follows: (20)
d
i
=
x
fin
,
i
-
x
fin
,
new
,
i
2
+
y
fin
,
i
-
y
fin
,
new
,
i
2
+
z
fin
,
i
-
z
fin
,
new
,
i
2
,where x
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
x
i
,
α
new
,
i,  y
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
y
i
,
α
new
,
i,  and z
fin
,
new
,
i
=
i
n
t
e
r
s
e
c
t
i
o
n
z
i
,
α
new
,
i are the coordinates of the new leaving point, with α
new
,
i being the new angle with respect to the horizontal of the aircraft i once the maneuvers associated to the solution under consideration are performed.w
5
∈
[
0,1
] again represents the relative importance of the dispersion term regarding the average value (a
5).


## 2.3.6. Objective 6: Minimizing Maneuver Dispersion

The last objective under consideration is related to the dispersion of maneuvers over time, that is, over multiple executions of the analysis. The aim is to share the maneuver effort among the aircrafts to attempt to avoid situations in which some aircrafts continuously perform maneuvers over time, whereas other aircrafts scarcely do it.For this, as mentioned before, we assume that a vector including the number of maneuvers performed by each aircraft from the moments they entered the aerial sector until the moment the analysis is carried out is available (m
a
n
1
,
…
,
m
a
n
i
,
…
,
m
a
n
n). Thus, a dispersion measure can be computed as follows: (21)
m
i
n
f
6
=
1
n
∑
i
=
1
n
a
6
-
m
a
n
n
e
w
i
,with a
6
=
(
∑
i
=
1
n
‍
m
a
n
new
,
i
)
/
n being the average number of maneuvers and (22)
m
a
n
new
,
i
=
m
a
n
i
,
i
f
v
c
i
=
0
∧
a
c
i
=
0
∧
t
c
i
=
0
m
a
n
i
+
1
,
i
f
v
c
i
=
1
∨
a
c
i
=
1
∨
t
c
i
=
1
,being the new number of maneuvers accumulated by aircraft i if the maneuvers associated to the solution under consideration are performed.


## 3. Multiobjective Simulated Annealing

Simulated annealing (SA) [20, 21] is a trajectorial metaheuristic which is named for and inspired by annealing in metallurgy.An initial feasible solution is randomly generated. In each iteration a new solutiony is randomly generated from the neighborhood of the current solution, y
∈
N
(
x
i
). If the new solution is better than the current one, then the algorithm moves to that solution (x
i
+
1
=
y); otherwise the movement to the worst solution is performed with certain probability. Note that accepting worse solutions allows for a more extensive search for the optimal solution and avoids trapping in local optima in early iterations. The probability of accepting a worse movement is a function of both a temperature factor and the change in the cost function as follows: (23)
p
=
e
-
f
y
-
f
x
i
/
T
i
,where T
i is the temperature in the ith iteration.The initial value of temperature (T
0) is high, which leads to a diversified search, since practically all movements are allowed. As the temperature decreases, the probability of accepting a worse movement falls. If the temperature is zero, then only better movements will be accepted, which makes simulated annealing work similar to hill climbing.Temperature is usually kept constant forL iterations and is then decreased after multiplying by a cooling rate (α
<
1). The algorithm stops when there has been a maximum number of iterations without accepting solutions.Metaheuristics have recently become very popular for multiobjective optimization. The aim is now to derive a good approximation of the efficient or Pareto set or, alternatively, take advantage of a decision-maker’s preferences to identify asatisficing efficient solution. There are several reasons that explain the increasing acceptance of SA and other metaheuristics; for instance, they converge speedily to Pareto-optimal solutions, handle both discrete and continuous problems with ease, and are less susceptible to the shape of the Pareto front.Multiobjective simulated annealing (MSA) was first proposed in [22]. The algorithm is analogous to the classical SA but now based on the concept of archiving the Pareto-optimal solutions and introducing a modification in the acceptance criteria of solutions, for which different approaches or a combination of them can be found in the literature aimed at increasing the probability of accepting nondominated solutions; see, for example, the method of Ulungu et al. (UMOSA) [23], the method of Suppapitnarm et al. (SMOSA) [24], or thePareto simulated annealing (PSA) [25]. A comparison of the above methods and other methods can be found in [26].


### 3.1. AMOSA Method

In this paper we consider thearchive simulated annealing-based multiobjective optimization algorithm (AMOSA) [27], which incorporates the concept ofarchive where the nondominated solutions generated are stored and determine the acceptance probability of a new solution taking into account thedomination status of the new solution (y) with the current one (x
i), as well as those in thearchive. For this purpose, theamount of domination measure is used [26] and defined as follows: given two solutions x
i and y, (24)
D
x
i
,
y
=
∏
j
=
1
m
f
j
x
i
-
f
j
y
R
j
,where m is the number of objectives and R
j is the range of the jth objective.Based on the domination status between current solutionx
i and new solution y, we can face the following situations:(1)
Ifx
i
dominates
y and k points from thearchive also dominate y, then x
i
+
1
=
y with probability (25)
p
=
e
-
D
a
v
g
/
T
,where (26)
D
a
v
g
=
∑
l
=
1
k
D
a
r
c
h
i
v
e
l
,
y
+
D
x
i
,
y
k
+
1denotes the average amount of domination of y by (
k
+
1
) points, namely, the current solution (x
i) and k points of thearchive.
(2)
Ify and x
i
are nondominating to each other, then we check the domination status of y and points in thearchive.
(a)
Ify is dominated by k solutions inarchive, then x
i
+
1
=
y with probability p
=
e
-
(
D
a
v
g
+
E
)
/
T, but now D
a
v
g
=
∑
l
=
1
k
D
a
r
c
h
i
v
e
l
,
y
/
k. E is a new element accounting for the possible redundancy associated to the incorporation of y to thearchive, as explained later.
(b)
Ify is nondominated by all solutions in thearchive, then x
i
+
1
=
y and add y to thearchive.
(c)
Ify dominates k solutions in thearchive, then x
i
+
1
=
y, add y to thearchive, and remove k dominated solutions from it.
(3)
Ify
dominates
x
i, then we check the domination status of y and points in thearchive.
(a)
Ify is dominated by k solutions in thearchive, we compute the minimum of the difference of domination amounts between y and the k points (Δ
d
o
m
min
⁡). Then, x
i
+
1 is the solution in thearchive corresponding to Δ
d
o
m
min
⁡ with probability (27)
p
=
1
1
+
e
-
Δ
d
o
m
min
⁡
.Else, x
i
+
1
=
y.
(b)
Ify does not dominate any solution in thearchive, then x
i
+
1
=
y, and add y to thearchive. If x
i is inarchive, then remove it from thearchive.
(c)
Ify dominates k solutions in thearchive, then x
i
+
1
=
y, add y to thearchive, and remove the k dominated solutions from it.ElementE in the above algorithm accounting for the redundancy would produce the inclusion of y in the archive, which is computed as follows: (28)
E
=
∑
j
=
1
m
C
j
R
j
2
,where C
j is the distance in each coordinate (objective) to the nearest solution to y and R
j is the range of each objective.As a result of the incorporation of the elementE, as the temperature decreases, solutions that increase the diversity in thearchive will be accepted. This makes the algorithm tend to converge faster, whereas in the exploitation phase it tends only to improve solutions and add diversity to the Pareto front.It is important to note that nondominated solutions are stored in thearchive up to a maximum, M
A. If the number of nondominated solutions exceeds M
A, clustering is applied to reduce the size to M
A [27].


#### 3.1.1. AMOSA Adaptation to the CDR Problem

To adapt AMOSA for the resolution of the CDR problem considered in this paper, we must first identify the way solutions are represented and how the neighborhood of a solution is defined.As mentioned in the mathematical modeling of the CDR problem, each solution consists on a vector with four elements per aircraft, three of them being binary elements representing the maneuver performed (VC, AC
⁡, or TC) and another element representing its magnitude (q
i). However, when the model was implemented, it was decided to replace maneuver elements by a single element with 4 possible values depending on the maneuver performed (VC
=
1, AC
⁡
=
2, TC
=
3, and NM
=
0
), where NM means no maneuver performed. We denote by t
i
j
∈
{
0,1
,
2,3
} and q
i
j
∈
[
-
1,1
] the type and magnitude of the maneuver performed by the aircraft j in the solution x
i, respectively.Besides, given solutionx
i, a new solution y is randomly generated from its neighborhood as follows. First, an aircraft j is uniformly selected from the solution, then we randomly decide if the type of maneuver for this aircraft is changed:(29)
t
i
+
1
j
=
mod
⁡
r
2
,
4
,
i
f
r
1
<
φ
t
i
j
,
e
l
s
e
,where φ is an algorithm parameter and r
1
,
r
2
~
U
[
0,1
].If value 1, 2, or 3 is randomly selected, then we randomly generate the magnitude of the maneuver:(30)
q
i
+
1
j
=
q
i
j
+
q
,
i
f
r
4
<
0.5
q
i
j
-
q
,
e
l
s
e
,where q
,
r
4
~
U
[
0,1
].The above expression works whent
i
+
1
j
=
t
i
j; that is, the maneuver type remains. Otherwise, (31)
q
i
+
1
j
=
q
,
i
f
r
4
<
0.5
-
q
,
e
l
s
e
.We must check thatq
i
+
1
j
∈
[
-
1,1
]  ∀
j.Besides, AMOSA parameters must be fixed by analyzing the algorithm performance for a representative instance set accounting for different sizes of the problem, that is, for different numbers of aircrafts. For this, we have carried out different tests accounting for different combinations for a subset of the algorithm parameters, analyzing the algorithm performance for each combination, as described in the following sections. We must distinguish between the algorithm parameters and those associated with the nature of the CDR problem, that is, parameters related to air traffic management.We have used a standard instance to fix parameters, also used by other approaches in the literature. The instance consists ofn of aircrafts with the same altitude and speed and equidistant from the center of a circle and with direction toward the center of that circle; see Figure 6. This implies that all aircrafts have conflicts between each other, since if no change is performed they will meet in the center of the circle. This instance is the most difficult one we can face and is often used to measure the performance of new approaches.Figure 6
Standard instance of the CDR problem.Aircraft parameters were fixed values as close to reality as possible using as a reference the features of anairbus A320, that is, minimum velocity = 440 km/h, maximum velocity = 870 km/h, initial velocity = 800 km/h, minimum altitude = 10 km, maximum altitude = 15 km, initial altitude = 12 km, and maximum angle = 45 grades. Horizontal and vertical security distances are 4.0234 km and 0.1524 km, respectively. The dimension of a Spanish aerial sector was used, that is,aerial sector radio = 250 km. The number of maneuvers performed by the aircrafts from the moment they entered the aerial sector until the moment the analysis is carried out (m
a
n
i
t) will consist on values 0 and 1 for even and odd aircrafts; that is, one maneuver at most has been performed by aircrafts before the execution of the algorithm in time t.The parameter values whose combinations were analyzed are initial temperature,T
0 (for 5 aircrafts it is 100, 1000, 1
0
4, 1
0
5, and 1
0
6, whereas for 10 and 20 aircrafts it is 1
0
4, 1
0
5, 1
0
6, 1
0
7, and 1
0
8); number of iterations the temperature is kept constant (L
=
5,10,50,100,200); cooling rate (α
=
0.8,0.9,0.95,0.98,0.99); probability of maneuver change (φ
=
1
/
2, 1/3, 1/4, 1/5, and 1/6); magnitude change (q
=
0.05,0.1,0.2,0.3,0.4); and maximum number of iterations without accepting solutions (convergence criterion) {
10,20,50,100,200
}.The following indicators were used to measure the algorithm performance:(1)
Minimum value for each of objective functionsf
1 and f
2: although the aim of a multiobjective optimization algorithm is to find compromise solutions equidistantly spread around the Pareto front, this indicator can give an idea of the exploitation capacity of the algorithm.
(2)
Number of nondominated solutions found: as the previous indicator, it is not a good measure of the algorithm performance, since many solutions could be found but which are not well spread around the Pareto front. However, it may provide a rough idea of the operating capability of the algorithm.
(3)
Dispersion of solutions: this measure provides information about the good or bad distribution of solutions around the Pareto front. However, it could occur that the algorithm derives a set of solutions which are very well spread around the Pareto front but that set includes very few solutions. Thus, this measure should be interpreted together with the previous one. Similarly, solutions could be well distributed but far from the Pareto front.
A redundancy measure (R) is used, which consists of taking the difference between the average distance of each solution and its nearest one: (32)
R
=
1
p
∑
i
=
1
p
m
-
d
i
,
min
⁡
,where m
=
(
1
/
p
)
∑
i
=
1
p
d
i
,
min
⁡ and d
i
,
min
⁡ is the minimum distance of the solution i to any other solution.
(4)
Dominance of solutions: it is computed as the mean of the hypervolume dominated by each solution,H
-. Combined with the previous two metrics, it is possible to compare the performance of different executions of the algorithm with different parameters.If we simultaneously consider the previous four indicators, it is not possible to derive the best parameter values. However, this problem can be simplified if we aggregate the metrics that measure the algorithm performance. If we consider the relative importance of indicators by means of (empirically obtained) weights, transform maximizing into minimizing objectives, normalize them, and discard atypical values when necessary, then we have the following metric:(33)
m
1
,
s
=
0.05
f
1
,
s
m
a
x
s
f
1
,
s
+
0.1
f
2
,
s
-
m
i
n
s
f
2
,
s
-
m
i
n
s
f
2
,
s
+
0.1
m
a
x
s
N
s
-
N
s
m
a
x
s
N
s
+
0.3
R
s
m
a
x
s
R
s
+
0.45
m
a
x
s
H
-
s
-
H
-
s
m
a
x
s
H
-
s
,where f
1
,
s and f
2
,
s are the values in f
1 and f
2 for the sth execution, N
s is the number of nondominated solutions derived in the sth execution, R
s is the redundancy measure, and H
-
s is the mean hypervolume dominated, respectively.A second metric represents the computation time:(34)
m
2
,
s
=
t
e
j
s
m
a
x
s
t
e
j
s
.


## 3.1. AMOSA Method

In this paper we consider thearchive simulated annealing-based multiobjective optimization algorithm (AMOSA) [27], which incorporates the concept ofarchive where the nondominated solutions generated are stored and determine the acceptance probability of a new solution taking into account thedomination status of the new solution (y) with the current one (x
i), as well as those in thearchive. For this purpose, theamount of domination measure is used [26] and defined as follows: given two solutions x
i and y, (24)
D
x
i
,
y
=
∏
j
=
1
m
f
j
x
i
-
f
j
y
R
j
,where m is the number of objectives and R
j is the range of the jth objective.Based on the domination status between current solutionx
i and new solution y, we can face the following situations:(1)
Ifx
i
dominates
y and k points from thearchive also dominate y, then x
i
+
1
=
y with probability (25)
p
=
e
-
D
a
v
g
/
T
,where (26)
D
a
v
g
=
∑
l
=
1
k
D
a
r
c
h
i
v
e
l
,
y
+
D
x
i
,
y
k
+
1denotes the average amount of domination of y by (
k
+
1
) points, namely, the current solution (x
i) and k points of thearchive.
(2)
Ify and x
i
are nondominating to each other, then we check the domination status of y and points in thearchive.
(a)
Ify is dominated by k solutions inarchive, then x
i
+
1
=
y with probability p
=
e
-
(
D
a
v
g
+
E
)
/
T, but now D
a
v
g
=
∑
l
=
1
k
D
a
r
c
h
i
v
e
l
,
y
/
k. E is a new element accounting for the possible redundancy associated to the incorporation of y to thearchive, as explained later.
(b)
Ify is nondominated by all solutions in thearchive, then x
i
+
1
=
y and add y to thearchive.
(c)
Ify dominates k solutions in thearchive, then x
i
+
1
=
y, add y to thearchive, and remove k dominated solutions from it.
(3)
Ify
dominates
x
i, then we check the domination status of y and points in thearchive.
(a)
Ify is dominated by k solutions in thearchive, we compute the minimum of the difference of domination amounts between y and the k points (Δ
d
o
m
min
⁡). Then, x
i
+
1 is the solution in thearchive corresponding to Δ
d
o
m
min
⁡ with probability (27)
p
=
1
1
+
e
-
Δ
d
o
m
min
⁡
.Else, x
i
+
1
=
y.
(b)
Ify does not dominate any solution in thearchive, then x
i
+
1
=
y, and add y to thearchive. If x
i is inarchive, then remove it from thearchive.
(c)
Ify dominates k solutions in thearchive, then x
i
+
1
=
y, add y to thearchive, and remove the k dominated solutions from it.ElementE in the above algorithm accounting for the redundancy would produce the inclusion of y in the archive, which is computed as follows: (28)
E
=
∑
j
=
1
m
C
j
R
j
2
,where C
j is the distance in each coordinate (objective) to the nearest solution to y and R
j is the range of each objective.As a result of the incorporation of the elementE, as the temperature decreases, solutions that increase the diversity in thearchive will be accepted. This makes the algorithm tend to converge faster, whereas in the exploitation phase it tends only to improve solutions and add diversity to the Pareto front.It is important to note that nondominated solutions are stored in thearchive up to a maximum, M
A. If the number of nondominated solutions exceeds M
A, clustering is applied to reduce the size to M
A [27].


### 3.1.1. AMOSA Adaptation to the CDR Problem

To adapt AMOSA for the resolution of the CDR problem considered in this paper, we must first identify the way solutions are represented and how the neighborhood of a solution is defined.As mentioned in the mathematical modeling of the CDR problem, each solution consists on a vector with four elements per aircraft, three of them being binary elements representing the maneuver performed (VC, AC
⁡, or TC) and another element representing its magnitude (q
i). However, when the model was implemented, it was decided to replace maneuver elements by a single element with 4 possible values depending on the maneuver performed (VC
=
1, AC
⁡
=
2, TC
=
3, and NM
=
0
), where NM means no maneuver performed. We denote by t
i
j
∈
{
0,1
,
2,3
} and q
i
j
∈
[
-
1,1
] the type and magnitude of the maneuver performed by the aircraft j in the solution x
i, respectively.Besides, given solutionx
i, a new solution y is randomly generated from its neighborhood as follows. First, an aircraft j is uniformly selected from the solution, then we randomly decide if the type of maneuver for this aircraft is changed:(29)
t
i
+
1
j
=
mod
⁡
r
2
,
4
,
i
f
r
1
<
φ
t
i
j
,
e
l
s
e
,where φ is an algorithm parameter and r
1
,
r
2
~
U
[
0,1
].If value 1, 2, or 3 is randomly selected, then we randomly generate the magnitude of the maneuver:(30)
q
i
+
1
j
=
q
i
j
+
q
,
i
f
r
4
<
0.5
q
i
j
-
q
,
e
l
s
e
,where q
,
r
4
~
U
[
0,1
].The above expression works whent
i
+
1
j
=
t
i
j; that is, the maneuver type remains. Otherwise, (31)
q
i
+
1
j
=
q
,
i
f
r
4
<
0.5
-
q
,
e
l
s
e
.We must check thatq
i
+
1
j
∈
[
-
1,1
]  ∀
j.Besides, AMOSA parameters must be fixed by analyzing the algorithm performance for a representative instance set accounting for different sizes of the problem, that is, for different numbers of aircrafts. For this, we have carried out different tests accounting for different combinations for a subset of the algorithm parameters, analyzing the algorithm performance for each combination, as described in the following sections. We must distinguish between the algorithm parameters and those associated with the nature of the CDR problem, that is, parameters related to air traffic management.We have used a standard instance to fix parameters, also used by other approaches in the literature. The instance consists ofn of aircrafts with the same altitude and speed and equidistant from the center of a circle and with direction toward the center of that circle; see Figure 6. This implies that all aircrafts have conflicts between each other, since if no change is performed they will meet in the center of the circle. This instance is the most difficult one we can face and is often used to measure the performance of new approaches.Figure 6
Standard instance of the CDR problem.Aircraft parameters were fixed values as close to reality as possible using as a reference the features of anairbus A320, that is, minimum velocity = 440 km/h, maximum velocity = 870 km/h, initial velocity = 800 km/h, minimum altitude = 10 km, maximum altitude = 15 km, initial altitude = 12 km, and maximum angle = 45 grades. Horizontal and vertical security distances are 4.0234 km and 0.1524 km, respectively. The dimension of a Spanish aerial sector was used, that is,aerial sector radio = 250 km. The number of maneuvers performed by the aircrafts from the moment they entered the aerial sector until the moment the analysis is carried out (m
a
n
i
t) will consist on values 0 and 1 for even and odd aircrafts; that is, one maneuver at most has been performed by aircrafts before the execution of the algorithm in time t.The parameter values whose combinations were analyzed are initial temperature,T
0 (for 5 aircrafts it is 100, 1000, 1
0
4, 1
0
5, and 1
0
6, whereas for 10 and 20 aircrafts it is 1
0
4, 1
0
5, 1
0
6, 1
0
7, and 1
0
8); number of iterations the temperature is kept constant (L
=
5,10,50,100,200); cooling rate (α
=
0.8,0.9,0.95,0.98,0.99); probability of maneuver change (φ
=
1
/
2, 1/3, 1/4, 1/5, and 1/6); magnitude change (q
=
0.05,0.1,0.2,0.3,0.4); and maximum number of iterations without accepting solutions (convergence criterion) {
10,20,50,100,200
}.The following indicators were used to measure the algorithm performance:(1)
Minimum value for each of objective functionsf
1 and f
2: although the aim of a multiobjective optimization algorithm is to find compromise solutions equidistantly spread around the Pareto front, this indicator can give an idea of the exploitation capacity of the algorithm.
(2)
Number of nondominated solutions found: as the previous indicator, it is not a good measure of the algorithm performance, since many solutions could be found but which are not well spread around the Pareto front. However, it may provide a rough idea of the operating capability of the algorithm.
(3)
Dispersion of solutions: this measure provides information about the good or bad distribution of solutions around the Pareto front. However, it could occur that the algorithm derives a set of solutions which are very well spread around the Pareto front but that set includes very few solutions. Thus, this measure should be interpreted together with the previous one. Similarly, solutions could be well distributed but far from the Pareto front.
A redundancy measure (R) is used, which consists of taking the difference between the average distance of each solution and its nearest one: (32)
R
=
1
p
∑
i
=
1
p
m
-
d
i
,
min
⁡
,where m
=
(
1
/
p
)
∑
i
=
1
p
d
i
,
min
⁡ and d
i
,
min
⁡ is the minimum distance of the solution i to any other solution.
(4)
Dominance of solutions: it is computed as the mean of the hypervolume dominated by each solution,H
-. Combined with the previous two metrics, it is possible to compare the performance of different executions of the algorithm with different parameters.If we simultaneously consider the previous four indicators, it is not possible to derive the best parameter values. However, this problem can be simplified if we aggregate the metrics that measure the algorithm performance. If we consider the relative importance of indicators by means of (empirically obtained) weights, transform maximizing into minimizing objectives, normalize them, and discard atypical values when necessary, then we have the following metric:(33)
m
1
,
s
=
0.05
f
1
,
s
m
a
x
s
f
1
,
s
+
0.1
f
2
,
s
-
m
i
n
s
f
2
,
s
-
m
i
n
s
f
2
,
s
+
0.1
m
a
x
s
N
s
-
N
s
m
a
x
s
N
s
+
0.3
R
s
m
a
x
s
R
s
+
0.45
m
a
x
s
H
-
s
-
H
-
s
m
a
x
s
H
-
s
,where f
1
,
s and f
2
,
s are the values in f
1 and f
2 for the sth execution, N
s is the number of nondominated solutions derived in the sth execution, R
s is the redundancy measure, and H
-
s is the mean hypervolume dominated, respectively.A second metric represents the computation time:(34)
m
2
,
s
=
t
e
j
s
m
a
x
s
t
e
j
s
.


## 3.1.1. AMOSA Adaptation to the CDR Problem

To adapt AMOSA for the resolution of the CDR problem considered in this paper, we must first identify the way solutions are represented and how the neighborhood of a solution is defined.As mentioned in the mathematical modeling of the CDR problem, each solution consists on a vector with four elements per aircraft, three of them being binary elements representing the maneuver performed (VC, AC
⁡, or TC) and another element representing its magnitude (q
i). However, when the model was implemented, it was decided to replace maneuver elements by a single element with 4 possible values depending on the maneuver performed (VC
=
1, AC
⁡
=
2, TC
=
3, and NM
=
0
), where NM means no maneuver performed. We denote by t
i
j
∈
{
0,1
,
2,3
} and q
i
j
∈
[
-
1,1
] the type and magnitude of the maneuver performed by the aircraft j in the solution x
i, respectively.Besides, given solutionx
i, a new solution y is randomly generated from its neighborhood as follows. First, an aircraft j is uniformly selected from the solution, then we randomly decide if the type of maneuver for this aircraft is changed:(29)
t
i
+
1
j
=
mod
⁡
r
2
,
4
,
i
f
r
1
<
φ
t
i
j
,
e
l
s
e
,where φ is an algorithm parameter and r
1
,
r
2
~
U
[
0,1
].If value 1, 2, or 3 is randomly selected, then we randomly generate the magnitude of the maneuver:(30)
q
i
+
1
j
=
q
i
j
+
q
,
i
f
r
4
<
0.5
q
i
j
-
q
,
e
l
s
e
,where q
,
r
4
~
U
[
0,1
].The above expression works whent
i
+
1
j
=
t
i
j; that is, the maneuver type remains. Otherwise, (31)
q
i
+
1
j
=
q
,
i
f
r
4
<
0.5
-
q
,
e
l
s
e
.We must check thatq
i
+
1
j
∈
[
-
1,1
]  ∀
j.Besides, AMOSA parameters must be fixed by analyzing the algorithm performance for a representative instance set accounting for different sizes of the problem, that is, for different numbers of aircrafts. For this, we have carried out different tests accounting for different combinations for a subset of the algorithm parameters, analyzing the algorithm performance for each combination, as described in the following sections. We must distinguish between the algorithm parameters and those associated with the nature of the CDR problem, that is, parameters related to air traffic management.We have used a standard instance to fix parameters, also used by other approaches in the literature. The instance consists ofn of aircrafts with the same altitude and speed and equidistant from the center of a circle and with direction toward the center of that circle; see Figure 6. This implies that all aircrafts have conflicts between each other, since if no change is performed they will meet in the center of the circle. This instance is the most difficult one we can face and is often used to measure the performance of new approaches.Figure 6
Standard instance of the CDR problem.Aircraft parameters were fixed values as close to reality as possible using as a reference the features of anairbus A320, that is, minimum velocity = 440 km/h, maximum velocity = 870 km/h, initial velocity = 800 km/h, minimum altitude = 10 km, maximum altitude = 15 km, initial altitude = 12 km, and maximum angle = 45 grades. Horizontal and vertical security distances are 4.0234 km and 0.1524 km, respectively. The dimension of a Spanish aerial sector was used, that is,aerial sector radio = 250 km. The number of maneuvers performed by the aircrafts from the moment they entered the aerial sector until the moment the analysis is carried out (m
a
n
i
t) will consist on values 0 and 1 for even and odd aircrafts; that is, one maneuver at most has been performed by aircrafts before the execution of the algorithm in time t.The parameter values whose combinations were analyzed are initial temperature,T
0 (for 5 aircrafts it is 100, 1000, 1
0
4, 1
0
5, and 1
0
6, whereas for 10 and 20 aircrafts it is 1
0
4, 1
0
5, 1
0
6, 1
0
7, and 1
0
8); number of iterations the temperature is kept constant (L
=
5,10,50,100,200); cooling rate (α
=
0.8,0.9,0.95,0.98,0.99); probability of maneuver change (φ
=
1
/
2, 1/3, 1/4, 1/5, and 1/6); magnitude change (q
=
0.05,0.1,0.2,0.3,0.4); and maximum number of iterations without accepting solutions (convergence criterion) {
10,20,50,100,200
}.The following indicators were used to measure the algorithm performance:(1)
Minimum value for each of objective functionsf
1 and f
2: although the aim of a multiobjective optimization algorithm is to find compromise solutions equidistantly spread around the Pareto front, this indicator can give an idea of the exploitation capacity of the algorithm.
(2)
Number of nondominated solutions found: as the previous indicator, it is not a good measure of the algorithm performance, since many solutions could be found but which are not well spread around the Pareto front. However, it may provide a rough idea of the operating capability of the algorithm.
(3)
Dispersion of solutions: this measure provides information about the good or bad distribution of solutions around the Pareto front. However, it could occur that the algorithm derives a set of solutions which are very well spread around the Pareto front but that set includes very few solutions. Thus, this measure should be interpreted together with the previous one. Similarly, solutions could be well distributed but far from the Pareto front.
A redundancy measure (R) is used, which consists of taking the difference between the average distance of each solution and its nearest one: (32)
R
=
1
p
∑
i
=
1
p
m
-
d
i
,
min
⁡
,where m
=
(
1
/
p
)
∑
i
=
1
p
d
i
,
min
⁡ and d
i
,
min
⁡ is the minimum distance of the solution i to any other solution.
(4)
Dominance of solutions: it is computed as the mean of the hypervolume dominated by each solution,H
-. Combined with the previous two metrics, it is possible to compare the performance of different executions of the algorithm with different parameters.If we simultaneously consider the previous four indicators, it is not possible to derive the best parameter values. However, this problem can be simplified if we aggregate the metrics that measure the algorithm performance. If we consider the relative importance of indicators by means of (empirically obtained) weights, transform maximizing into minimizing objectives, normalize them, and discard atypical values when necessary, then we have the following metric:(33)
m
1
,
s
=
0.05
f
1
,
s
m
a
x
s
f
1
,
s
+
0.1
f
2
,
s
-
m
i
n
s
f
2
,
s
-
m
i
n
s
f
2
,
s
+
0.1
m
a
x
s
N
s
-
N
s
m
a
x
s
N
s
+
0.3
R
s
m
a
x
s
R
s
+
0.45
m
a
x
s
H
-
s
-
H
-
s
m
a
x
s
H
-
s
,where f
1
,
s and f
2
,
s are the values in f
1 and f
2 for the sth execution, N
s is the number of nondominated solutions derived in the sth execution, R
s is the redundancy measure, and H
-
s is the mean hypervolume dominated, respectively.A second metric represents the computation time:(34)
m
2
,
s
=
t
e
j
s
m
a
x
s
t
e
j
s
.


## 4. Parameter Setting and Performance Analysis with 5 Aircrafts

The results of several executions with different combination of parameters regarding both metrics are shown in Figure7, whereas Figure 8 shows nondominated points (executions) in Figure 7.Figure 7
Performance metric (see (33)) and computation time metric (see (34)).Figure 8
Nondominated and the selected (in red) solution regarding both metrics. Note that we want to minimize performance and computation time metrics, which are shown in Figure7.Twelve nondominated points were found. The combination of parameters we chose (pointed out in red) due to its very good performance and very low computation time (0.1156 seconds) isT
0
=
100, L
=
5, α
=
0.99, q
=
0.4, φ
=
1
/
4, and convergence iterations = 20.The valuations in the indicators for the selected combination aref
1
,
s
=
0.1947, f
2
,
s
=
-
960.7386 (0.0012 normalized), number of nondominated solutions = 115.2 (0.0938 normalized), R
s
=
78.7474 (0.0629 normalized), H
-
s
=
60.9578 (0.1343 normalized), and t
e
j
s
=
0.1156 seconds (m
2
,
s
=
0.00073), with 3583.8 iterations carried out (0.0064 normalized).We will now provide different figures showing the evolution of some interesting elements and the different individual objectives along the algorithm execution with the selected parameter combination. Note that each objective function will be individually considered but in all cases we ensure no collision (r
i
j
≤
0,  ∀
i
,
j).Figure9 shows the evolution of the acceptance probability of new solutions. At the beginning of the algorithm the acceptance probability is very close to 1. This means that in the first iterations the algorithm has a great capacity for exploration in the search space, accepting nondominated solutions that dominate the current one with any value of redundancy as well as solutions dominated by the current one. However, as the number of iterations progresses, the acceptance probability gradually decreases, leading the algorithm to an exploitation or intensification state, in which only solutions that dominate the current one or nondominated solutions with a redundancy value very small are accepted.Figure 9
Evolution of acceptance probability of new solutions.Figure10 shows the evolution of the number of nondominated solutions. It is possible to see that the number of nondominated solutions increases along the execution of the algorithm. An even higher number of iterations would imply a greater number of nondominated solutions. However, not surprisingly, there is little practical difference between, for example, 100 and 5000 solutions if we assume a similar distribution of both set points around the Pareto front.Figure 10
Evolution of the number of nondominated solutions.Figure11 shows the evolution of the minimum value for f
1 that ensures no collision. We can see that it improves from 0.4 to 0.19. This means that aircrafts in the solution with minimum value for f
1 would be forced to perform an average magnitude change maneuver of 20%, assuming that the variance between them is not taken into account (w
1
=
0) in this case.Figure 11
Evolution of the minimum value for maneuver magnitudes(
f
1
).Figure12 shows the evolution of the minimum values for f
2. The average collision risk (w
2
=
0) improves from a value close to -
1 to value -
30. This means that we go from a situation in which though the average risk collision is negative, there are some collisions, to another safer situation in which there are no collisions and slack in the security distances. This low risk has been reached on a solution in which several turn change (TC) maneuvers have been performed, since this maneuver more effectively reduces the risk.Figure 12
Evolution of the minimum value for collision risks(
f
2
).Regarding the evolution of the minimum values forf
3, the minimum and maximum number of maneuvers throughout the execution are 4 and 5, respectively; that is, less than 4 maneuvers cannot be performed to avoid collisions. Note that in the CDR instance under consideration the aircrafts have the same altitude and speed and are equidistant from the center of a circle and with direction toward the center of that circle and that minimum number of maneuvers which ensure no collision is 4.If we now pay attention to the evolution of the minimum values forf
4 that ensures no collision, the minimum value achieved for time delay is 0, since a solution involving only altitude change maneuvers was found. Note also that we have assumed that maneuvers are instantly performed.Figure13 shows the evolution of the minimum values for f
5 that ensures no collision. We can see that the distances between the points at which aircrafts will leave the aerial sector and the expected points according to their initial trajectory (when entering the aerial sector) decrease from 10 km to almost 0. This means that a solution in which the aircraft maneuvers are only velocity changes (VC) has not been found, since this maneuver ensures that the real and estimated leaving points are the same.Figure 13
Evolution of the minimum value for deviations in the leaving points(
f
5
).Finally, Figure14 shows the evolution of the minimum values for f
6 that ensures no collision. Maneuver dispersion begins at 0.65 and decreases to a value close to 0.17. Note that in the parameters setting we assumed that the number of maneuvers performed by the 5 aircrafts from the moment they entered the aerial sector until the moment the analysis is carried out was {
1,0
,
1,0
,
1
}; that is, aircrafts 1, 3, and 5 had previously performed one maneuver.Figure 14
Evolution of the minimum values for maneuver dispersion(
f
6
).As mentioned before, when we analyzedf
3, at least 4 maneuvers have to be performed to avoid collisions, leading to three possible dispersion values, the lowest being 0.17.Next, we will analyze the performance of the proposed algorithm by comparing the set of nondominated solutions derived from it with the real Pareto front, which can be computed using brute force for this instance with 5 aircrafts having conflicts between each other.Of course, we discretized maneuver magnitude to 200 values in the range[
-
1,1
]; that is, a step of 0.01 is considered. Besides, we considered all the possible combinations of maneuvers, 4
5. Only those solutions for which there was no conflict were stored, leading to 258719616 solutions (14,326 Gbs). This set was iteratively reduced, leading to 6641 nondominated solutions. Clearly, we obtained an approximation of the real Pareto front since a discretization in the maneuver magnitudes was done.The metric for comparison is the relationship of the hypervolume dominated by each nondominated solution derived by the algorithm with hypervolume dominated by the solution derived using brute force closest to the first, among those that dominate it:(35)
m
=
∑
i
=
1
m
h
a
l
g
,
i
h
o
p
t
,
min
⁡
,where h
o
p
t
,
min
⁡
=
m
i
n
j
=
1
n
{
h
o
p
t
,
j
,
h
a
l
g
,
i
} and h
o
p
t
,
j are the solutions that dominate h
a
l
g
,
i.We carried out 100 executions of the algorithm with the optimal parameter values identified above for 5 aircrafts. The algorithm performance (on average) is 0.25 worse regarding Pareto front. If a higher number of values were considered in the discretization of the maneuver magnitude a higher number of nondominated would be obtained and the performance would approach Pareto front. Figures15–19 show an example of nondominated solutions derived by the algorithm (in red) and using brute force (in blue) for the different objectives with respect to f
2, the collision risks. In these figures, we can see how the nondominated solutions derived by the algorithm are spread around the Pareto front.Figure 15
Pareto front and its approximation when the objectives shown on figure are maneuver magnitudes(
f
1
) and collision risks (
f
2
).Figure 16
Pareto front and its approximation when the objectives shown on figure are numbers of maneuvers(
f
3
) and collision risks (
f
2
).Figure 17
Pareto front and its approximation when the objectives shown on figure are time delays(
f
4
) and collision risks (
f
2
).Figure 18
Pareto front and its approximation when the objectives shown on figure are deviations in the leaving points(
f
5
) and collision risks (
f
2
).Figure 19
Pareto front and its approximation when the objectives shown on figure are maneuver dispersion(
f
6
) and collision risks (
f
2
).


## 5. Parameter Setting for 10 and 20 Aircrafts and Illustrative Example

In this section we first identify the best combination of parameter values for 10 and 20 aircrafts, analogously as in the previous section for 5. Then, we use a randomly generated instance with 20 aircrafts to illustrate the algorithm and its flexibility.First, we consider 10 aircrafts and the possible values forT
0: {
1
0
4
,
1
0
5
,
1
0
6
,
1
0
7
,
1
0
8
}. Figure 20 shows the values regarding the performance and computation time metrics for several executions, whereas nondominated points (executions) are shown in Figure 21.Figure 20
Performance and computation time metrics for 10 aircrafts.Figure 21
Nondominated and the selected (in red) solution for 10 aircrafts.The two points located at the right were discarded since they have a very high execution time for a tiny variation of performance with respect to the others. The following combination of parameters was selected:T
0
=
10000, L
=
5, α
=
0.8, t
=
0.05, φ
=
1
/
6, and convergence iterations = 100.The valuations in the different indicators by the selected combination aref
1
,
s
=
0.2776, f
2
,
s
=
-
154.9622 (0.0092 normalized), number of nondominated solutions = 48.8000 (0.1025 normalized), R
s
=
32.1758 (0.0498 normalized), H
-
s
=
607.0908 (0.0329 normalized), and t
e
j
s
=
0.0964 seconds (m
2
,
s
=
0.0019), with 1912 iterations carried out (0.0019 normalized).Regarding 20 aircrafts, the values regarding the performance and computation time metrics for several executions are shown in Figure22, whereas nondominated points (executions) are shown in Figure 23.Figure 22
Performance and computation time metrics for 20 aircrafts.Figure 23
Nondominated and selected (in red) solutions for 20 aircrafts.The combination of parameters chosen, see Figure23, was the following: T
0
=
10000, L
=
5, α
=
0.8, t
=
0.1, φ
=
1
/
5, and convergence iterations = 200.As mentioned before, the most difficult instance has been used to identify the combination of parameter values that yields the best performance of the algorithm. As such, once that optimal parameter set has been identified we will illustrate the algorithm for another randomly generated instance with 20 aircrafts; see Figures24 and 25.Figure 24
Random instance with 20 aircrafts, PlaneXY.Figure 25
Random instance with 20 aircrafts, PlaneXZ.The algorithm is then executed with the optimal parameters set, leading to a set of nondominated solutions. If only one solution has to be implemented we need to incorporate the preferences of a decision-maker (DM) to derive acompromise (nondominated) solution. Specifically, we assume that the DM only wants to minimize the delay (f
4) and the deviations in the leaving points (f
5), with weights 1000 and 1
/
4, respectively; that is, we consider the following expression: (36)
f
=
1000
×
f
4
+
1
4
×
f
5
.The compromise solution for the previous weighted function is (VC, 0.43, NM, −, VC, 0.26, TC, 0.07, NM, −, VC, -
0.08, AC
⁡, 0.95, VC, 0.07, VC, -
0.27, VC, -
0.49, VC, -
0.13, TC, 0.39, VC, 0.47, NM, −, TC, -
0.24, TC, 0.83, TC, -
0.06, VC, -
0.48, AC
⁡, 0.88, VC, -
0.38) with objective values f
1
=
0.381293, f
2
=
-
152.273331, f
3
=
17.0, f
4
=
0.013386, f
5
=
81.067757, and f
6
=
2.75.Let us pay attention to the aircraft in blue in Figures24 and 25, whose initial configuration is position = (
116.67
,
-
47.70,12.15
), entry point = (
255.00
,
-
90.91,12
), estimated leaving point = (
-
255.00,68.41,12
), current velocity = 833.84 km/h, initial velocity = 900.00 km/h, current angle = 2.83 grades, initial angle = 3.21 grades, maneuvers performed since it entered the aerial sector = 4, elapsed time since it entered the aerial sector = 0.16, and estimated time traveling in the aerial sector = 0.59.The compromise solution implies that this aircraft must perform a velocity change maneuver (VC) of magnitude 0.43. That is, its velocity will increase by 43%. This is consistent, since the current velocity before applying the algorithm was 833.84 km/h, and the initial velocity when it entered the aerial sector was 900 km/h. It could be interpreted as the aircraft increases its velocity to try to counteract the delay it would cause if we maintain the current velocity. Therefore, the aircraft configuration changes after running the algorithm and now the current velocity is 904.91 km/h and maneuvers performed so far are 5.The developed algorithm is very flexible since objectives can individually or in combination be optimized. For example, if we now want to optimizef
6 (minimizing maneuver dispersion), the following compromise solution is reached: (NM, −, NM, −, NM, −, NM, −, NM, −, NM, −, AC
⁡, -
0.20, AC
⁡, 0.09, TC, 0.57, NM, −, TC, 0.74, AC
⁡, -
1.00, NM, −, AC
⁡, 0.07, VC, -
0.25, NM, −, AC
⁡, 0.68, VC, -
0.06, TC, 0.78, NM, −), with objective values f
1
=
0.444926, f
2
=
-
151.998688, f
3
=
10.0, f
4
=
0.013348, f
5
=
86.590012, and f
6
=
2.36.We realize how, regardless of the maneuver performed, the algorithm tends to homogenize the number of maneuvers performed by the aircrafts. The minimum value forf
6 is 2.36 in comparison with value 2.75 obtained where f
4 and f
5 were optimized. The maneuvers performed before running the algorithm were (4, 9, 9, 10, 7, 8,3,4,3, 9,1,2, 6,1,3, 9,8,4,4, 10), whereas maneuvers performed later are (4, 9, 9, 10, 7, 8,4,5,4, 9,2,3, 6,2,4, 9,9,5,5, 10). Highlighted in bold are those aircrafts that have performed a maneuver.


## 6. Conclusions

We have dealt with the conflict detection and resolution (CDR) problem in air traffic management accounting for three admissible maneuvers (velocity, turn, and altitude changes). Moreover, a multiobjective context has been considered, accounting for the minimization of the maneuver number and magnitude, collision risks, time delays, deviations in the leaving points, and maneuver dispersion.Both the possibility of performing three types of maneuvers and the multiobjective context make this paper an original contribution regarding previous works on CDR problems.We have adapted thearchive simulated annealing-based multiobjective optimization algorithm to solve the CDR problem and conducted a thorough analysis of the effect of parameters in the resolution based on some instances of the problem, which led to a subset of optimal values for such parameters. The results derived by a discretized version of the algorithm for 5 aircrafts were compared with the real Pareto front computed using brute force, showing a good approximation. Moreover, a more complex random instance with 20 aircrafts has been considered to illustrate the algorithm and show its flexibility.As future research work, we propose improving and extending some of the model features and its resolution. We have assumed that the altitude changes are instantaneous. In general, the planes are far away enough from each other so that the time required to perform an altitude change is negligible. However, in some cases the algorithm could fail to detect a conflict.Another improvement could be allowing aircrafts to perform more than one maneuver simultaneously (e.g., VC and TC at once). This would mean better performances in some objectives. Moreover, the evaluation objectives, the most time-consuming part in the algorithm, could be parallelized to improve execution times to come near to output real-time solutions.Finally, the good performance of the approach proposed has been proven. However, its performance should be compared with the application of other metaheuristics for the specific CDR problem under consideration.



---

*Source: 8738014-2016-08-02.xml*