Hindawi Publishing Corporation Mathematical Problems in Engineering Volume 2016, Article ID 8738014, 16 pages http://dx.doi.org/10.1155/2016/8738014 Research Article Multiobjective Simulated Annealing for Collision Avoidance in ATM Accounting for Three Admissible Maneuvers A.MateosandA.Jiménez-Martín DecisionAnalysisandStatisticsGroup,DepartamentodeInteligenciaArtificial,UniversidadPolite´cnicadeMadrid, CampusdeMontegancedos/n,BoadilladelMonte,28660Madrid,Spain CorrespondenceshouldbeaddressedtoA.Mateos;[email protected] Received8January2016;Accepted6June2016 AcademicEditor:MartinJ.Geiger Copyright©2016A.MateosandA.Jime´nez-Mart´ın. This is an open access article distributed under the Creative Commons AttributionLicense,whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkis properlycited. Technologicaladvancesarerequiredtoaccommodateairtrafficcontrolsystemsforthefuturegrowthofairtraffic.Particularly, detectionandresolutionofconflictsbetweenaircraftsisaproblemthathasattractedmuchattentioninthelastdecadebecoming vitaltoimprovethesafetystandardsinfreeflightunstructuredenvironments.Weproposeusingthearchivesimulatedannealing- basedmultiobjectiveoptimizationalgorithmtodealwithsuchaproblem,accountingforthreeadmissiblemaneuvers(velocity, turn,andaltitudechanges)inamultiobjectivecontext.Theminimizationofthemaneuvernumberandmagnitude,timedelays,or deviationsintheleavingpointsareconsideredforanalysis.Theoptimalvaluesforthealgorithmparametersetareidentifiedinthe morecomplexinstanceinwhichallaircraftshaveconflictsbetweeneachotheraccountingfor5,10,and20aircrafts.Moreover,the performanceoftheproposedapproachisanalyzedbymeansofacomparisonwiththeParetofront,computedusingbruteforcefor 5aircraftsandthealgorithmisalsoillustratedwitharandominstancewith20aircrafts. 1.Introduction mustflyalongpredeterminedroutesthroughcertain waypoints,whicharesetbytheagenciesofairtraffic Cargo and air traffic (AT) congestion has experienced a control(ATC),somethingthatusuallyfailstoproduce general exponential growth throughout the world over the optimalresults.Aircraftsarenotallowedtoflydirectly lastdecade.Everyminuteoftheday,bothmorningandafter- totheirfinaldestinationtakingadvantageoffavorable noon,thereareabout11,000aircraftsintheairsomewherein winds without making changes to their trajectories theworld,ascanbeseeninrealtimeathttps://www.flight- causingunnecessaryfuelcosts,whichcanindirectly radar24.com/. cause increases in ticket prices. This problem is 2014 was the first year in which 100,000 flights per day particularlyevidentintransoceanicroutes,whichare wereexceeded.Europe’slargestairportshandleabout2,000 experiencingthegreatestgrowthindemand. dailytakeoffsandlandings.Thistrendcontinuestoincrease gradually and estimates predict bending movements until (2)IncreasedATCworkload:ATcontrollershave,among 2030. other functions, to prevent collisions between air- Withthesystemscurrentlyavailable,theairtrafficcontrol craftsandredirectroutestoavoidadverseconditions. agenciesarenotabletoefficientlymanagethislargeincrease In congested areas, such as regions near to airports whichistakingplaceduetoseveralfactorsasfollows: calledterminalradarapproachcontrols(TRACONs), (1)Efficient use of airspace: currently, the airspace is AT controllers often simplify their high workloads rigidly structured and with a large number of con- makingaircraftmaintaindefaultroutesoutsidethese straints that aircrafts have to comply with. They regions,causingdelaysinlandingsandtakeoffs. 2 MathematicalProblemsinEngineering (3)Slowcommunication:communicationisrestrictedto Simulation techniques have also been used to handle atediousvoicecommunicationbetweenaircraftand CDRproblems.Forinstance,[5]analyzestheeconomicper- AT controllers, causing frequent bottleneck situa- formance of a specific conflict resolution strategy based on tions. velocitychangebetweentwoaircraftsintermsofextratime andfuelconsumption.Theworkin[6]alsoconsidersaveloc- ity regulation problem, but from a different perspective, In view of the problems described, the aviation com- distinguishingbetweencrossingconflicts(thewider),inwhich munity has been working in recent years on a concept the aircrafts intersect at some point and security cylinders called free flight. This innovative concept allows pilots to overlap, and conflicts trail, caused when an aircraft pursues choose their own routes, altitude, and velocity to reduce another,bothwithdifferentvelocities. delays and manage the use of aircraft fuel more efficiently. Neural networks have been also used for performing The preferences of the pilots will be restricted only in very velocitychangesinCDRproblems[5,7,8]. congested airspace areas or to prevent unauthorized entry Morerecently,[9]focusesonmixed-integeroptimization intomilitaryareas. modelsbasedonvelocityregulation.Theyproposetoaccel- Free flight is potentially possible due to the availability erateordecelerateduringaspecifiedtimeinterval,reverting of technologies such as global positioning systems (GPS); back to the original velocity once the conflict is avoided. communications data links such as automatic dependence They propose a heuristic procedure where the problem is surveillance-broadcast(ADS-B);detectionsystemsandcolli- decomposedandlocallyexactlysolved. sionavoidance;andpowerfulcomputingincreasinglybeing Otherlessfrequentproposalsconsiderturnchangesthat implementedinaircrafts. leadtononlinearoptimizationmodels.Forinstance,[10]pro- Inaddition,thereareseveraldecisionsupporttoolsthat posesatwo-stepapproach.First,anonconvexmixed-integer reduce the workload of both AT controllers and pilots and nonlinearoptimizationisusedtominimizetheweightedair- optimizetheircapacity,forexample,thedetectionandreso- craftanglevariations.Then,asetofunconstrainedquadratic lution of conflicts in airspace sectors, landings and takeoffs optimization models are considered, where aircrafts are management in airports, and organizational systems of the forcedtoreturntotheiroriginalflightplanassoonaspossible workload of AT controllers to better organize their tasks to once there is no aircraft in conflict with any other. Both an increaseproductivity. exact and an approximate resolution are proposed. In the These technological advances will also allow current second,theturnchangesarediscretizedtoreducethesearch ATC systems to accommodate the future growth of air space. traffic. Algorithms to detect and solve aircraft conflicts are Different metaheuristics have been proposed for solv- vital to improve the safety standards in free flight unstruc- ing CDR models accounting for turn changes, such as tured environments. These systems can be used on land by ant colony systems [11, 12], genetic algorithms [13], variable ATC or by the flight management system (FMS) of each neighborhood search [14], and particle swarm optimization aircraft. [15], which uses a series of waypoints the aircrafts can pass Inthispaper,wefocusonthedevelopmentofalgorithmic through. toolsforaircraftconflictdetectionandresolution(CDR)prob- A pretty realistic proposal is described in [16], wherein lem.Weassumethateachaircraftissurroundedbycylinder theaccelerationvariableisaddedtothemodel.Itintendsto representingasecurityvirtualvolume.Conflictbetweentwo solveconflictsdiscretizingthetimeremaininguntilitoccurs aircraftsoccurswhentherespectiveaircraftsecurityvolumes at different intervals, optimizing acceleration, and velocity overlap. thatshouldbeassignedtoeachaircraft.Anonlinearmixed Different approaches can be found in literature to deal 0-1 model is used to solve the problem, which is iteratively withcollisionavoidanceaccountingfordifferentnumberand linearized by using Taylor polynomials. This approach is typesofadmissiblemaneuversforaircraftsandwithdifferent then enhanced in [17], extending control to aircraft outside solution approaches, including the use of exact solvers, the aviation sector to manage, that is, taking into account simulation techniques, and metaheuristics. The work in [1] those aircrafts leaving it or entering it. Moreover, they take presentasurveywiththemostimportantoftheseuptothe intoaccounttheconflictsthatmayarisewhenanaircraftis year2000,whereas[2]focusesonapproachesfrom2000up climbingordescendingtochangealtitude. to2012. The work in [18] improves the velocity change model Oneofthefirstapproachestodealwithcollisionavoid- by adding altitude changes when necessary, for example, in ance was [3]. A path planning problem among given way- head-to-head conflict situations. A multiobjective perspec- points avoiding all possible conflicts was considered aimed tiveisconsideredincludingobjectivessuchasvelocityvari- atminimizingthetotalflighttime.Twomixed-integerlinear ationandtotalnumberofmaneuversandforcingtoreturnto programs were proposed accounting for velocity changes theoriginalflightconfigurationwhennoaircraftsareincon- andanglechangesasadmissiblemaneuvers,respectively.The flict.Anexactlysolvedmixed0-1linearoptimizationmodel work in [4] proposes a three-dimensional formulation as is used, with small computational time for the execution a mixed-integer nonlinear program in which only velocity makingitsuitableforreal-timeuse. changeswereadmissible.CPLEXwasusedfortheresolution In[19]aninnovativepointofviewbasedonthechoiceof inbothapproaches. differentstrategiestoavoidconflictsisproposed.Anoriginal MathematicalProblemsinEngineering 3 trajectorymodelusingB-splinesisintroducedtogetherwith anewsemi-infiniteprogrammingformulationofthesepara- tionconstraintinvolvedinCDRproblems. In this paper we propose using simulated annealing to deal with a CDR problem accounting for three admis- sible maneuvers (velocity, turn, and altitude changes) in a multiobjective context. Specifically, the archive simu- lated annealing-based multiobjective optimization algorithm (AMOSA)hasbeenadaptedtotheCDRproblemaccounting for objectives such as minimizing the maneuver number and magnitude, time delays, or deviations in the leaving points. Boththepossibilityofperformingthreetypesofmaneu- Figure1:Aerialsectorwithseveralaircrafts. vers and the multiobjective context make this paper an original contribution regarding previous works on CDR problems. The paper is structured as follows. The mathematical Security modeling for the multiobjective problem under considera- volume 1000feet tion is introduced in Section 2, including the identification of parameters, decision variables, and constraints and the description and modelization of the candidate objective 5 miles functionsforanalysis.AMOSAanditsadaptationtothecon- sidered CDR problem are described in Section 3. Section 4 Aircraft dealswiththeparametersettingandtheperformanceanalysis when5aircraftsareconsidered.Theparametersettingfor10 and20aircraftsandanexampleillustratingtheflexibilityof the proposed algorithm are provided in Section 5. Finally, someconclusionsareprovidedinSection6. Figure2:Securityvolumeofanaircraft. 2.MathematicalModeling Forexample,theTCmaneuverisefficientregardingfuelcon- Weassumethatinaparticularmomentthereare𝑛aircrafts sumption,butitmakestheaircraftleaveitsoriginaltrajectory, in an aerial sector, a cubic volume in the space managed whichconstitutesaninconvenience. by an AT controller; see Figure 1. We have to decide which Moreover,weconsideramultiobjectiveperspectiveofthe maneuvers to perform to avoid possible collisions between problemincludingtheminimizationofthemaneuvernum- them. We also assume that the decision on maneuvers will ber,magnitudeanddispersion,timedelays,anddeviationsin beeffectiveuntil𝑛aircraftsleavetheaerialsectororuntila theleavingpoints. new aircraft enters the aerial sector. At the moment a new Thedispersionofthemaneuversobjectiveaimsatspread- aircraftenterstheaerialsector,theanalysisweproposeshould ingtheeffortofallaircraftmaneuversinanattempttoavoid beagaincarriedouttoidentifynewmaneuverstoavoidnew situations in which some aircrafts continuously perform possiblecollisionscausedbytheenteringaircraft. maneuvers over time whereas other aircrafts do not do so. Moreover, we assume that the number of maneuvers The deviations in the leaving points objective minimize the performed by the aircrafts from the moment they entered sum of the distances between the theoretical leaving points the aerial sector until the moment the analysis is carried according to the initial trajectory when entering the aerial out is known. This information will be useful when ana- sector and the real leaving points after possible maneuver lyzing the dispersion of maneuvers objective, as described performances. afterwards. The collision avoidance problem can be mathematically modeledasfollows. 2.1. Parameters and Decision Variables. We consider 𝑛 air- Aconflictbetweentwoaircraftsoccursifthehorizontal crafts in an aerial sector at the time 𝑡. First, parameters and/orverticaldistancesbetweenthemaresmallerthansome accountingforthefeaturesoftheaircraftsareasfollows: gcyivlienndrsieccaulrfiotyrmli;msietes,Fmigaukrieng2.thTheusse,ctuhreitysevcoulruitmyecyaldinodpetras Vmin,𝑖:minimumvelocityfortheaircraft𝑖. correspondingtoanycoupleofaircraftsshouldneverinter- Vmax,𝑖:maximumvelocityfortheaircraft𝑖. secttoavoidconflicts. 𝑧min,𝑖:minimumaltitudefortheaircraft𝑖. Intheapproachwepropose,weaccountforthreetypes of aircraft maneuvers: velocity change (VC), altitude change 𝑧max,𝑖:maximumaltitudefortheaircraft𝑖. (AC), and turn change (TC). A best type of maneuver does 𝛽max,𝑖:maximumvariationoftheanglefortheaircraft notexistsinceeachmaneuverhasadvantagesanddrawbacks. 𝑖. 4 MathematicalProblemsinEngineering Secondly,initialparametersforeachaircraftwhenenter- ui−uj ingtheaerialsector(denotedbysubscriptini)areasfollows: Vini,𝑖:velocityoftheaircraft𝑖uponentry. 𝑧ini,𝑖:altitudeoftheaircraft𝑖uponentry. 𝑥ini,𝑖:absciseoftheaircraft𝑖uponentry. 𝑦ini,𝑖:ordinateoftheaircraft𝑖uponentry. 𝛼ini,𝑖: angle of the aircraft 𝑖 with respect to the u horizontaluponentry. j Shadow Next,weconsiderfinalparametersforeachaircraftwhen leaving the aerial sector assuming that no maneuvers have been performed, that is, according to its initial trajectory whenenteringtheaerialsector,asfollows: 𝑥fin,𝑖:estimatedabsciseoftheaircraft𝑖whenleaving ui theaerialsector. 𝑦fin,𝑖:estimatedordinateoftheaircraft𝑖whenleaving theaerialsector. 𝑧fin,𝑖:estimatedaltitudeoftheaircraft𝑖whenleaving Figure3:Detectingconflictsituations. theaerialsector. The parameters accounting for the configuration of the aircraftsatthetime𝑡areasfollows: V𝑖𝑡:velocityoftheaircraft𝑖atthetime𝑡. allowedtoperform,atmost,onemaneuver,and𝑞𝑖 ∈ [−1,1] 𝑧𝑡:aircraft𝑖altitudeatthetime𝑡. sincethechangescanbenegativeorpositive. 𝑖 Consequently,asolutionwillconsistofavectorwithfour 𝑥𝑖𝑡:aircraft𝑖absciseatthetime𝑡. elements per aircraft, identifying the maneuver performed 𝑦𝑡:aircraft𝑖ordinateatthetime𝑡. anditsmagnitude. 𝑖 𝑡𝑖:timenecessarytoarrivetotheboundoftheaerial sectorwithaconstantvelocityV𝑡. 2.2. Constraints. As a first approach, the main constraint 𝑖 should be avoiding conflicts between aircrafts; that is, the 𝛼𝑖𝑡:anglewithrespecttothehorizontaloftheaircraft intersection between security cylinders should always be 𝑖attheinstant𝑡. empty for any two aircrafts in the aerial sector under con- 𝑡 sideration. A geometric construction for detecting conflict man:numberofmaneuversperformedbytheaircraft 𝑖 𝑖atthetime𝑡sinceitenteredintheaerialsector. situationsisconsidered[3,10];seeFigure3. Thevelocityvectorsoftwoaircrafts𝑖and𝑗are Finally,parametersaccountingforsecuritydistancesand thecollisionriskforaircraftsareasfollows: 𝑢⃗𝑖 =(Vnew,𝑖×cos(𝛼new,𝑖),Vnew,𝑖×sin(𝛼new,𝑖)), 𝑠hor,𝑖:horizontalsecuritydistancefortheaircraft𝑖. (1) 𝑠ver,𝑖:verticalsecuritydistancefortheaircraft𝑖. 𝑢⃗𝑗 =(Vnew,𝑗×cos(𝛼new,𝑗),Vnew,𝑗×sin(𝛼new,𝑗)), 𝑐hor,ver:relativeimportancebetweenverticalandhor- izontalrisk. respectively,whereVnew,𝑖 isthenewvelocityconsideringthe According to current standards, the horizontal security VCmaneuverand𝛼new,𝑖isthenewanglewiththehorizontal distanceisusually5nauticalmiles,whereastheverticalsecu- planebyaddingthemagnitudeoftheTCmaneuver. ritydistanceis1000feet(seeFigure2),butothervaluescould Thebasicideaofthemodelcomesfromtheconstruction beusedintheanalysis. of the relative velocity vector 𝑢⃗𝑖 − 𝑢⃗𝑗; see Figure 3. The Regardingthedecisionvariables,asmentionedbefore,we twostraightlinesparalleltotherelativevelocityvectorand propose three types of aircraft maneuvers: velocity change tangent to the security circle of aircraft 𝑗 (dotted lines in (VC), altitude change (AC), and turn change (TC). Thus, Figure 3) define a region where the intersection with the we consider three binary variables vc𝑖, ac𝑖, and tc𝑖 for each trajectory for aircraft 𝑖 is a segment named the shadow aircraft 𝑖, pointing out whether a velocity, altitude, or turn segment. changeisperformed,respectively,andacontinuousvariable A horizontal conflict occurs if the security cylinder of 𝑞𝑖 representing the magnitude (proportion) of the change aircraft𝑖intersectstheshadowsegmentgeneratedbyaircraft performed.Notethatvc𝑖+ac𝑖+tc𝑖 ≤1 ∀𝑖sinceeachaircraftis 𝑗or,onthecontrary,since𝑢⃗𝑖−𝑢⃗𝑗and𝑢⃗𝑗−𝑢⃗𝑖areparallels. MathematicalProblemsinEngineering 5 Z Security r j distance d ij Aircraft 𝛼 ij Security distance g wij lij ij Figure5:Verticalsecuritydistanceandturnvector. Figure4:Mainanglesanddistances. Finally,thefollowingconstraintscheckwhetherthenew velocity, altitude, and turn satisfy the features of the corre- spondingaircraft: V ≤V ≤V , Considering now the cutting planes that are tangent to min,𝑖 new,𝑖 max,𝑖 bothcylinders(seeFigure4)andtheangles𝑔𝑖𝑗 and𝑙𝑖𝑗,there 𝑧min,𝑖 ≤𝑧new,𝑖 ≤𝑧max,𝑖, (3) isnoconflictifoneofthefollowingtwoconditionsholds: 𝛼 −𝛼 ≤𝛽 , new,𝑖 𝑖 max,𝑖 Vnew,𝑖×sin(𝛼new,𝑖)−Vnew,𝑗×sin(𝛼new,𝑗) ≥tan(𝑙𝑖𝑗), wtiohne.reVnew,𝑖, 𝑧new,𝑖,and𝛼new,𝑖 arethenewaircraftconfigura- Vnew,𝑖×cos(𝛼new,𝑖)−Vnew,𝑗×cos(𝛼new,𝑗) (2) 2.3.ObjectiveFunctions. Sixdifferentobjectivefunctionswill VVnneeww,,𝑖𝑖××csoins((𝛼𝛼nneeww,,𝑖𝑖))−−VVnneeww,,𝑗𝑗××scions((𝛼𝛼nneeww,,𝑗𝑗)) ≤tan(𝑔𝑖𝑗). bntieimtuceoddneesoliadfyemsr,eaddneefvouiravteairnosna,lscyiosnilslti:shsiepoenlecarifivisicknasgl,lypn,oumimnintbsie,mraniozdfintmghaetnhmeeuamnveearugs--, ver dispersion. Further information and the mathematical modelingofsuchobjectivefunctionsareprovidedbelow. Notethatthefunctionsattheleftoftheaboveexpressions can cause a zero denominator. These cases are referred 2.3.1.Objective1:MinimizingManeuverMagnitudes. Itmakes to as model pathological situations and produce unstable sensetoclaimthataircraftsperformmaneuversassmoothly solutions since a conflict between two aircrafts may be aspossibletoavoidconflicts;thatis,abruptmaneuversthat erroneouslydeterminedduetothenulldenominator,forcing may disturb passengers or even be dangerous should be the aircrafts to crash in the worse case. This situation is avoided. According to the above, the dispersion magnitude detectedwhen|𝑥𝑖 −𝑥𝑗| < 𝑠hor,𝑖 +𝑠hor,𝑗.Therefore,variables maneuvers performed by aircrafts should be incorporated 𝛼new,𝑖, 𝛼new,𝑗 and parameters 𝑙𝑖𝑗 and 𝑔𝑖𝑗, which represent into the analysis. Moreover, a high dispersion should be angles,arerotated𝜋/2radianswhencomputingexpressions penalized,thatis,situationsinwhichmaneuvermagnitudes forhorizontalriskdetection to overcomesuch pathological areveryhighforsomeaircraftsandverylowforothers. situation. Thefirstobjectivefunction,𝑓1,canbethenmodeledby Verticalconflictsaredetectedmoreeasilyconsideringthe thesumoftheaveragemaneuvermagnitudeandadispersion securitycylinders.Computingtheverticaldistancebetween term: two aircrafts can detect these conflicts. Figure 5 shows the modHeoliwngevoefr,thveevreyrtriecsatlrcicotnivfleictssit.uations may occur in an min 𝑓1 =𝑎1+ 𝑤𝑛1∑𝑛 󵄨󵄨󵄨󵄨𝑎1−󵄨󵄨󵄨󵄨𝑞𝑖󵄨󵄨󵄨󵄨󵄨󵄨󵄨󵄨, (4) 𝑖=1 aerial sector at certain times. For example, it could be very difficulttofindafeasiblesolutionwithout(eitherhorizontal where or vertical) conflictswhen there is a high density of nearby acoirncsrtarftasin.tTh, werheifcohreb,ewcoemheasveanrealadxdeidtiothnealcooblljiesciotinveavfuoindcatniocne 𝑎1 = 𝑛1∑𝑛 󵄨󵄨󵄨󵄨𝑞𝑖󵄨󵄨󵄨󵄨, (5) 𝑖=1 from now on, as described in the next section. This allows us to better explore the solution space but, in contrast, it and 𝑤1 represents the relative importance of the dispersion complicatesreachingtheoptimalsolution. termregardingtheaveragemaneuvermagnitude. 6 MathematicalProblemsinEngineering 2.3.2.Objective2:MinimizingCollisionRisks. Thisobjective whichconcludesthatthereisnohorizontalconflictifoneof is the constraint we decided to relax. As cited before, this thefollowingtwoconditionsholds: allowsabetterexplorationofthesolutionspace.Moreover,if wcaenccoonmsipduetretawcoonsoflliucttiroinsksimnewashuicrehftohreeraecahroefntohecmonaflciccotusnwte- Vnew,𝑖×sin(𝛼new,𝑖)−Vnew,𝑗×sin(𝛼new,𝑗) ≥tan(𝑙𝑖𝑗), ingforthedistancesbetweenallpairsofaircrafts. Vnew,𝑖×cos(𝛼new,𝑖)−Vnew,𝑗×cos(𝛼new,𝑗) There is a conflict between a couple of aircrafts if their collisionriskispositive𝑟𝑖𝑗 >0;otherwise(𝑟𝑖𝑗 ≤0)thereisno Vnew,𝑖×sin(𝛼new,𝑖)−Vnew,𝑗×sin(𝛼new,𝑗) (10) conflict. Vnew,𝑖×cos(𝛼new,𝑖)−Vnew,𝑗×cos(𝛼new,𝑗) Wedifferentiatethesituationinwhichthereisatleastone conflict between a couple of aircrafts and the one in which ≤tan(𝑔𝑖𝑗). there are no conflicts to compute an average collision risk (𝑎2).Inthefirstcase,𝑎2iscomputedtakingintoaccountonly Thehorizontalcollisionrisk(𝑟hor,𝑖𝑗)iscomputedas those𝑟𝑖𝑗 >0,avoidingthatnegativevaluesnullifythepositive. I𝑎n2,twhehiscehcoisnadnceagsaet(iv𝑟𝑖e𝑗v≤al0u,e∀: 𝑖,𝑗),all𝑟𝑖𝑗 areusedtocompute 𝑟hor,𝑖𝑗 =min{𝛾𝑖𝑗,𝛿𝑖𝑗}, 1 𝑛 𝑛 𝛾𝑖𝑗 =tan(𝑙𝑖𝑗) 𝑎2 ={{{{{{{{{{𝑘𝑘11∑∑𝑖=𝑛1𝑗=∑∑𝑛𝑖+1𝑟𝜑𝑖𝑗𝑖𝑗,, iiff ∃𝑟𝑖𝑗a≤t l0ea∀s𝑖t,𝑗one 𝑟𝑖𝑗 >0, (6) − VVnneeww,,𝑖𝑖××csoins((𝛼𝛼nneeww,,𝑖𝑖))−−VVnneeww,,𝑗𝑗××scions((𝛼𝛼nneeww,,𝑗𝑗)), (11) { 2𝑖=1𝑗=𝑖+1 𝛿 = Vnew,𝑖×sin(𝛼new,𝑖)−Vnew,𝑗×sin(𝛼new,𝑗) 𝑖𝑗 where Vnew,𝑖×cos(𝛼new,𝑖)−Vnew,𝑗×cos(𝛼new,𝑗) 𝜑𝑖𝑗 ={{𝑟𝑖𝑗, if 𝑟𝑖𝑗 >0 (7) −tan(𝑔𝑖𝑗). {0, otherwise, Thus, it is possible to assess how far two aircrafts are to and𝑘1and𝑘2arethenumbersofelementsconsideredinthe horizontallyinvadetheirrespectivesecuritycylindersor,in average computations, that is, 𝑘1 = 𝑛!/(𝑛 − 2)!4, since all caseofconflict,tomeasuretheintensityofinvasion. pairsofaircraftsareconsidered,whereas𝑘2isthenumberof Thenewaircraftconfigurationchangesornotdepending onthemaneuverbeingperformed: coupleswithconflict. The above two situations are also considered in the modelizationofobjectivefunction,𝑓2.Theaveragecollision {{{{V𝑖, if vc𝑖 =0 r(𝑎is2k≤is0m),inthime diziesdperinsiobnotohncthaseesc,olbliusitoninritshkevsaelucoensdiscaalssoe Vnew,𝑖 ={{{{V𝑖+𝑞𝑖(V𝑖−Vmin,𝑖), if vc𝑖 =1∧𝑞𝑖 ≤0 considered: {V𝑖+𝑞𝑖(Vmax,𝑖−V𝑖), if vc𝑖 =1∧𝑞𝑖 >0, (12) {{{𝑎2, if 𝑎2 >0 𝛼 ={𝛼𝑖, if tc𝑖 =0 min 𝑓2 ={{{𝑎2+ 𝑤𝑘2∑𝑛 ∑𝑛 󵄨󵄨󵄨󵄨󵄨𝑎2−𝑟𝑖𝑗󵄨󵄨󵄨󵄨󵄨, else. (8) new,𝑖 {{𝛼𝑖+𝑞𝑖𝛽max,𝑖, if tc𝑖 =1. { 1𝑖=1𝑗=𝑖+1 Parameters 𝑙𝑖𝑗 and 𝑔𝑖𝑗 are computed on the basis of the Next, we clarify how the collision risk for a couple of anglesanddistancesshowninFigure4: aircrafts is computed. The collision risk 𝑟𝑖𝑗 between two 𝑙 =𝜔 +𝜃 , aircrafts𝑖and𝑗iscomputedinfourdifferentways,depending 𝑖𝑗 𝑖𝑗 𝑖𝑗 onthecombinationsofverticalandhorizontalconflicts: 𝑔 =𝜔 −𝜃 , 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑟 𝑖𝑗 𝑦 −𝑦 ={{{{{{{{{{{{{{𝑟𝑟𝑟hhhooorrr,,,𝑖𝑖𝑖𝑗𝑗𝑗,++𝑐𝑐hhoorr,,vveerr××𝑟𝑟vveerr,,𝑖𝑖𝑗𝑗,, iiifff 𝑟𝑟𝑟hhhooorrr,,,𝑖𝑖𝑖𝑗𝑗𝑗 ≤>≤000∧∧∧𝑟𝑟𝑟vvveeerrr,,,𝑖𝑖𝑖𝑗𝑗𝑗 ≤>>000 (9) 𝜔𝜃𝑖𝑖𝑗𝑗 ==aarrccstainn((𝑠𝑥h𝑖𝑖or−,𝑖𝑑𝑥+𝑖𝑗𝑗𝑗𝑠/)h2o,r,𝑗/2). (13) {𝑐hor,ver×𝑟ver,𝑖𝑗, if 𝑟hor,𝑖𝑗 >0∧𝑟ver,𝑖𝑗 ≤0. Theverticalcollisionrisk(𝑟ver,𝑖𝑗)justconsistsofcomputing the vertical distance between aircrafts added to the vertical Tocomputehorizontalcollisionrisk,wetakeintoaccount securitydistance,sothatifthereisnoconflict,then𝑟ver,𝑖𝑗 ≤0 the analysis of the horizontalconflict shown in Section 2.2, andviceversa: MathematicalProblemsinEngineering 7 𝑟ver,𝑖𝑗 ={{{mmaaxx{{𝑧𝑧nneeww,,𝑗𝑖−−𝑧𝑧nneeww,,𝑗𝑖++𝑠𝑠vveerr,,𝑗𝑖,,𝑧𝑧nneeww,,𝑖𝑗−−𝑧𝑧nneeww,𝑗,𝑖++𝑠𝑠vveerr,,𝑗𝑖}},, iiff 𝑧𝑧nneeww,,𝑖𝑖 >≤𝑧𝑧nneeww,,𝑗𝑗, (14) where𝑧new,𝑖 takesdifferentvaluesdependingonwhetheror 2.3.5.Objective5:MinimizingDeviationsintheLeavingPoints. notthismaneuverisperformedandonwhetherthealtitude The point at which an aircraft will leave the aerial sector isincreasedordecreased: may differ from the expected point according to its initial trajectory(whenenteringtheaerialsector)asaconsequence {{{{0, if ac𝑖 =0 of the maneuvers performed. As such, the aim is now minimizing the sum of the distances between both leaving 𝑧new,𝑖 ={{{{𝑧𝑖+𝑞𝑖(𝑧𝑖−𝑧min,𝑖), if ac𝑖 =1∧𝑞𝑖 ≤0 (15) points for each aircraft. Moreover, we will again take into {𝑧𝑖+𝑞𝑖(𝑧max,𝑖−𝑧𝑖), if ac𝑖 =1∧𝑞𝑖 >0. accoThunisttohbejedcistipveersfiuonncotifosnu,c𝑓h5d,iisstacnocmesp.uted by adding the average(𝑎5)andthedispersionofthedistancesbetweenthe 2.3.3. Objective 3: Minimizing the Numbers of Maneuvers. original leaving point and the new leaving point once the The objective is to minimize the number of maneuvers maneuversassociatedtothesolutionunderconsiderationare performed by 𝑛 aircrafts. This is directly connected with performed: the AT controllers workload, since they communicate the correspondingmaneuverstothepilots.Therefore,thisgoalis ethqiusiovbaljeencttitvoemfuinncitmioinzinwgetjhuestchoanvterotollesrusm’wuoprkthloeatdh.rTeoebminoadreyl min 𝑓5 =𝑎5+ 𝑤𝑛5∑𝑛 󵄨󵄨󵄨󵄨𝑎5−𝑑𝑖󵄨󵄨󵄨󵄨 with 𝑎5 = 𝑛1∑𝑛 𝑑𝑖, (19) 𝑖=1 𝑖=1 variablesassociatedtopossiblemaneuversoftheaircraftsin thesolutionunderconsideration wherethedistanceiscomputedasfollows: 𝑛 min 𝑓3 =∑(vc𝑖+ac𝑖+tc𝑖). (16) 𝑑 𝑖=1 𝑖 (20) Noteagainthatvc𝑖 +ac𝑖 +tc𝑖 ≤ 1 ∀𝑖sinceeachaircraftcan =√(𝑥fin,𝑖−𝑥fin,new,𝑖)2+(𝑦fin,𝑖−𝑦fin,new,𝑖)2+(𝑧fin,𝑖−𝑧fin,new,𝑖)2, performatmostonetypeofmaneuver. where 𝑥fin,new,𝑖 = intersection(𝑥𝑖,𝛼new,𝑖), 𝑦fin,new,𝑖 = 2.3.4.Objective4:MinimizingTimeDelays. Thetimeanair- intersection(𝑦𝑖,𝛼new,𝑖), and𝑧fin,new,𝑖 = intersection(𝑧𝑖,𝛼new,𝑖) craftwillleavetheaerialsectormaydifferfromtheexpected are the coordinates of the new leaving point, with 𝛼new,𝑖 time according to its initial trajectory (when entering the being the new angle with respect to the horizontal of the aerialsector)ifaVCorTCmaneuverisperformed.Then,the aircraft𝑖oncethemaneuversassociatedtothesolutionunder aimisnowminimizingthesumofthedelaysfortheaircrafts. considerationareperformed. Moreover,wewillagaintakeintoaccountthedispersionof 𝑤5 ∈[0,1]againrepresentstherelativeimportanceofthe suchdelays. dispersiontermregardingtheaveragevalue(𝑎5). This objective function, 𝑓4, is computed by adding the averagedelay(𝑎4)andthedispersionofthedelays: 2.3.6.Objective6:MinimizingManeuverDispersion. Thelast objectiveunderconsiderationisrelatedtothedispersionof min 𝑓4 =𝑎4+ 𝑤𝑛4∑𝑛 󵄨󵄨󵄨󵄨𝑎4−󵄨󵄨󵄨󵄨𝑡new,𝑖−𝑡𝑖󵄨󵄨󵄨󵄨󵄨󵄨󵄨󵄨, (17) manaanlyesuivs.erThsoevaeirmtimiset,othshaatries,tohveemrmanueltuivpelereexffeocruttaiomnosnogftthhee 𝑖=1 aircraftstoattempttoavoidsituationsinwhichsomeaircrafts where𝑎4 =(1/𝑛)∑𝑛𝑖=1|𝑡new,𝑖−𝑡𝑖|and continuously perform maneuvers over time, whereas other aircraftsscarcelydoit. For this, as mentioned before, we assume that a vector √(𝑥 −𝑥)2+(𝑦 −𝑦)2 𝑡 = fin,𝑖 𝑖 fin,𝑖 𝑖 , (18) includingthenumberofmaneuversperformedbyeachair- new,𝑖 V craft from the moments they entered the aerial sector until new,𝑖 the moment the analysis is carried out is available (man1, isthetimeaircraft𝑖leavestheaerialsectoroncethemaneu- ...,man𝑖,...,man𝑛). Thus, a dispersion measure can be versassociatedtothesolutionunderconsiderationareper- computedasfollows: formed.Notethat𝑡new,𝑖 =𝑡𝑖ifanACmaneuverisperformed. ttohe𝑤Id1nisap𝑓ne4d,rs𝑤𝑤io24nin∈teor[bm0j,e1rce]tgivraeerpdfruiennsgcetntihtosentsahv𝑓ee1rraaegnlaedtdi𝑓ve2el,airyme,sappnoearclttoaigvnoeclueys.olyf min 𝑓6 = 𝑛1∑𝑖=𝑛1󵄨󵄨󵄨󵄨󵄨𝑎6−mannew𝑖󵄨󵄨󵄨󵄨󵄨, (21) 8 MathematicalProblemsinEngineering with 𝑎6 = (∑𝑛𝑖=1mannew,𝑖)/𝑛 being the average number of [23],themethodofSuppapitnarmetal.(SMOSA)[24],orthe maneuversand Paretosimulatedannealing (PSA)[25].Acomparisonofthe abovemethodsandothermethodscanbefoundin[26]. mannew,𝑖 {man𝑖, if vc𝑖 =0∧ac𝑖 =0∧tc𝑖 =0 (22) 3.1. AMOSA Method. In this paper we consider the archive = { simulated annealing-based multiobjective optimization algo- {man𝑖+1, if vc𝑖 =1∨ac𝑖 =1∨tc𝑖 =1, rithm (AMOSA) [27], which incorporates the concept of archive where the nondominated solutions generated are beingthenewnumberofmaneuversaccumulatedbyaircraft 𝑖ifthemaneuversassociatedtothesolutionunderconsider- stored and determine the acceptance probability of a new solutiontakingintoaccountthedominationstatusofthenew ationareperformed. solution(𝑦)withthecurrentone(𝑥𝑖),aswellasthoseinthe archive.Forthispurpose,theamountofdominationmeasure 3.MultiobjectiveSimulatedAnnealing isused[26]anddefinedasfollows:giventwosolutions𝑥𝑖and 𝑦, Simulated annealing (SA) [20, 21] is a trajectorial meta- heuristic which is named for and inspired by annealing in 𝑚 [𝑓 (𝑥)−𝑓 (𝑦)] 𝑗 𝑖 𝑗 metAalnluirngiyt.ialfeasiblesolutionisrandomlygenerated.Ineach 𝐷𝑥𝑖,𝑦 =∏𝑗=1 𝑅𝑗 , (24) iteration a new solution 𝑦 is randomly generated from the neighborhoodofthecurrentsolution,𝑦 ∈ 𝑁(𝑥𝑖).Ifthenew where𝑚isthenumberofobjectivesand𝑅𝑗istherangeofthe solution is better than the current one, then the algorithm 𝑗thobjective. movestothatsolution(𝑥𝑖+1 =𝑦);otherwisethemovementto Basedonthedominationstatusbetweencurrentsolution theworstsolutionisperformedwithcertainprobability.Note 𝑥𝑖andnewsolution𝑦,wecanfacethefollowingsituations: that accepting worse solutions allows for a more extensive searchfortheoptimalsolutionandavoidstrappinginlocal (1)If𝑥𝑖 dominates𝑦and𝑘pointsfromthearchivealso optima in early iterations. The probability of accepting a dominate𝑦,then𝑥𝑖+1 =𝑦withprobability worsemovementisafunctionofbothatemperaturefactor andthechangeinthecostfunctionasfollows: 𝑝=𝑒−𝐷avg/𝑇, (25) 𝑝=𝑒−(𝑓(𝑦)−𝑓(𝑥𝑖))/𝑇𝑖, (23) where wheThree𝑇𝑖inisittihalevtaelmuepeorfatteumrepienratthuer𝑖eth(𝑇i0te)riastihoingh.,whichleads 𝐷avg = ∑𝑘𝑙=1𝐷arc𝑘hi+ve𝑙1,𝑦+𝐷𝑥𝑖,𝑦 (26) to a diversified search, since practically all movements are allowed. As the temperature decreases, the probability of denotes the average amount of domination of 𝑦 by acceptingaworsemovementfalls.Ifthetemperatureiszero, (𝑘+1)points,namely,thecurrentsolution(𝑥𝑖)and then only better movements will be accepted, which makes 𝑘pointsofthearchive. simulatedannealingworksimilartohillclimbing. Temperatureisusuallykeptconstantfor𝐿iterationsand (2)If𝑦and𝑥𝑖 arenondominatingtoeachother,thenwe is then decreased after multiplying by a cooling rate (𝛼 < check the domination status of 𝑦 and points in the 1). The algorithm stops when there has been a maximum archive. numberofiterationswithoutacceptingsolutions. (a)If𝑦isdominatedby𝑘solutionsinarchive,then Metaheuristics have recently become very popular for multiobjectiveoptimization.Theaimisnowtoderiveagood 𝑥𝑖+1 = 𝑦withprobability𝑝 = 𝑒−(𝐷avg+𝐸)/𝑇,but approximationoftheefficientorParetosetor,alternatively, now𝐷avg =∑𝑘𝑙=1𝐷archive𝑙,𝑦/𝑘.𝐸isanewelement takeadvantageofadecision-maker’spreferencestoidentify accountingforthepossibleredundancyassoci- asatisficing efficientsolution.Thereareseveralreasonsthat atedtotheincorporationof𝑦tothearchive,as explain the increasing acceptance of SA and other meta- explainedlater. heuristics; for instance, they converge speedily to Pareto- (b)If 𝑦 is nondominated by all solutions in the optimalsolutions,handlebothdiscreteandcontinuousprob- archive,then𝑥𝑖+1 =𝑦andadd𝑦tothearchive. lems with ease, and are less susceptible to the shape of the (c)If 𝑦 dominates 𝑘 solutions in the archive, then Paretofront. 𝑥𝑖+1 = 𝑦, add 𝑦 to the archive, and remove 𝑘 Multiobjectivesimulatedannealing(MSA)wasfirstpro- dominatedsolutionsfromit. posedin[22].ThealgorithmisanalogoustotheclassicalSA butnowbasedontheconceptofarchivingthePareto-optimal (3)If𝑦dominates𝑥𝑖,thenwecheckthedominationstatus solutions and introducing a modification in the acceptance of𝑦andpointsinthearchive. criteria of solutions, for which different approaches or a combinationofthemcanbefoundintheliteratureaimedat (a)If𝑦isdominatedby𝑘solutionsinthearchive, increasing the probability of accepting nondominated solu- we compute the minimum of the difference tions;see,forexample,themethodofUlunguetal.(UMOSA) of domination amounts between 𝑦 and the 𝑘 MathematicalProblemsinEngineering 9 points (Δdommin). Then, 𝑥𝑖+1 is the solution 500 in the archive corresponding to Δdom with min 400 probability 300 1 𝑝= . (27) 200 1+𝑒−Δdom min 100 Else,𝑥𝑖+1 =𝑦. 0 (b)If 𝑦 does not dominate any solution in the −100 archive,then𝑥𝑖+1 =𝑦,andadd𝑦tothearchive. −200 If 𝑥𝑖 is in archive, then remove it from the archive. −300 (c)If 𝑦 dominates 𝑘 solutions in the archive, then −400 𝑥𝑖+1 =𝑦,add𝑦tothearchive,andremovethe𝑘 −500 dominatedsolutionsfromit. −500 −400 −300 −200 −100 0 100 200 300 400 500 Figure6:StandardinstanceoftheCDRproblem. Element 𝐸 in the above algorithm accounting for the redundancywouldproducetheinclusionof𝑦inthearchive, 𝑗isuniformlyselectedfromthesolution,thenwerandomly whichiscomputedasfollows: decideifthetypeofmaneuverforthisaircraftischanged: 𝐸=√𝑗∑𝑚=1(𝐶𝑅𝑗𝑗)2, (28) 𝑡𝑖𝑗+1 ={{{𝑡m𝑖𝑗,od(𝑟2,4), eifls𝑟e1, <𝜑 (29) where𝐶𝑗isthedistanceineachcoordinate(objective)tothe where𝜑isanalgorithmparameterand𝑟1,𝑟2 ∼𝑈[0,1]. nearestsolutionto𝑦and𝑅𝑗istherangeofeachobjective. Ifvalue1,2,or3israndomlyselected,thenwerandomly Asaresultoftheincorporationoftheelement𝐸,asthe generatethemagnitudeofthemaneuver: temperaturedecreases,solutionsthatincreasethediversityin ttoohnelcyoanrtcovheiirvmgeepwfraoilsvlteebre,sowalcuhcteeiropentaessdia.nnThdthieasdemdxpadkloievisteatrthsiioetnyalptgohoartsihteheimtPtaetrenendtdos 𝑞𝑖𝑗+1 ={{{𝑞𝑞𝑖𝑖𝑗𝑗+−𝑞𝑞,, iefls𝑟e4, <0.5 (30) front. where𝑞,𝑟4 ∼𝑈[0,1]. It is importanttonotethatnondominatedsolutionsare The above expression works when 𝑡𝑗 = 𝑡𝑗; that is, the storedinthearchiveuptoamaximum,MA.Ifthenumberof 𝑖+1 𝑖 maneuvertyperemains.Otherwise, nondominatedsolutionsexceedsMA,clusteringisappliedto reducethesizetoMA[27]. 𝑞𝑖𝑗+1 ={{𝑞, if 𝑟4 <0.5 (31) 3.1.1. AMOSA Adaptation to the CDR Problem. To adapt {−𝑞, else. AMOSAfortheresolutionoftheCDRproblemconsidered in this paper, we must first identify the way solutions are Wemustcheckthat𝑞𝑗 ∈[−1,1] ∀𝑗. 𝑖+1 represented and how the neighborhood of a solution is Besides,AMOSAparametersmustbefixedbyanalyzing defined. the algorithm performance for a representative instance set As mentioned in the mathematical modeling of the accounting for different sizes of the problem, that is, for CDR problem, each solution consists on a vector with four different numbers of aircrafts. For this, we have carried out elements per aircraft, three of them being binary elements different tests accounting for different combinations for a representingthemaneuverperformed(VC,AC,orTC)and subsetofthealgorithmparameters,analyzingthealgorithm another element representing its magnitude (𝑞𝑖). However, performance for each combination, as described in the fol- whenthemodelwasimplemented,itwasdecidedtoreplace lowingsections.Wemustdistinguishbetweenthealgorithm maneuverelementsbyasingleelementwith4possiblevalues parametersandthoseassociatedwiththenatureoftheCDR depending on the maneuver performed(VC = 1, AC = 2, problem,thatis,parametersrelatedtoairtrafficmanagement. TC = 3, and NM = 0), where NM means no maneuver We have used a standard instance to fix parameters, performed. We denote by 𝑡𝑗 ∈ {0,1,2,3} and 𝑞𝑗 ∈ [−1,1] alsousedbyotherapproachesintheliterature.Theinstance 𝑖 𝑖 the type and magnitude of the maneuver performed by the consists of 𝑛 of aircrafts with the same altitude and speed aircraft𝑗inthesolution𝑥𝑖,respectively. andequidistantfromthecenterofacircleandwithdirection Besides,givensolution𝑥𝑖,anewsolution𝑦israndomly toward the center of that circle; see Figure 6. This implies generatedfromitsneighborhoodasfollows.First,anaircraft thatallaircraftshaveconflictsbetweeneachother,sinceifno 10 MathematicalProblemsinEngineering changeisperformedtheywillmeetinthecenterofthecircle. where 𝑚 = (1/𝑝)∑𝑝𝑖=1𝑑𝑖,min and 𝑑𝑖,min is the mini- Thisinstanceisthemostdifficultonewecanfaceandisoften mumdistanceofthesolution𝑖toanyothersolution. usedtomeasuretheperformanceofnewapproaches. (4)Dominanceofsolutions:itiscomputedasthemean Aircraftparameterswerefixedvaluesasclosetorealityas ofthehypervolumedominatedbyeachsolution,𝐻. possibleusingasareferencethefeaturesofanairbusA320, Combinedwiththeprevioustwometrics,itispossible thatis,minimumvelocity=440km/h,maximumvelocity= tocomparetheperformanceofdifferentexecutionsof 870km/h, initial velocity = 800km/h, minimum altitude = thealgorithmwithdifferentparameters. 10km, maximum altitude = 15km, initial altitude = 12km, and maximum angle = 45 grades. Horizontal and vertical If we simultaneously consider the previous four indica- securitydistancesare4.0234kmand0.1524km,respectively. tors, it is not possible to derive the best parameter values. The dimension of a Spanish aerial sector was used, that However, this problem can be simplified if we aggregate is, aerial sector radio = 250km. The number of maneuvers the metrics that measure the algorithm performance. If we performed by the aircrafts from the moment they entered consider the relative importance of indicators by means of the aerial sector until the moment the analysis is carried (empirically obtained) weights, transform maximizing into 𝑡 out (man𝑖) will consist on values 0 and 1 for even and odd minimizingobjectives,normalizethem,anddiscardatypical aircrafts;thatis,onemaneuveratmosthasbeenperformed valueswhennecessary,thenwehavethefollowingmetric: byaircraftsbeforetheexecutionofthealgorithmintime𝑡. a1r0e4,Thi1n0iet5ip,aalarntaedmm1ep0tee6rr,avwtauhlrueeer,se𝑇aws0hfo(ofsroer1c0o5maanbidricn2raa0tftioasnirisctwraiesftrse10ai0tn,iasl1y01z00e04d,, 𝑚1,𝑠 =0.05max𝑓𝑠1(,𝑠𝑓1,𝑠) +0.1𝑓2,−𝑠m−imn𝑠in(𝑓𝑠(2,𝑓𝑠)2,𝑠) ar1at0ut5er,e(1𝛼0is6=,k1e00p.7t8,,ca0on.9nd,s0t1a.90n58t,);0(.𝐿n9u8m,=0b.9e9r5),;o1pf0r,iot5eb0ra,ab1tii0loi0tny,s2o0tfh0)em;tacenomeouplivenergr- +0.1mamx𝑠a(x𝑁𝑠(𝑠𝑁)−𝑠)𝑁𝑠 +0.3max𝑅𝑠𝑠(𝑅𝑠) (33) (c𝑞han=ge0.(0𝜑5,0=.11,0/2.2,,10/.33,,01/.44,);1/a5n,danmdax1/i6m)u;mmangunmitubdeer ochfaitnegre- +0.45max𝑠(𝐻𝑠)−𝐻𝑠, ations without accepting solutions (convergence criterion) max𝑠(𝐻𝑠) {10,20,50,100,200}. Thefollowingindicatorswereusedtomeasurethealgo- where 𝑓1,𝑠 and 𝑓2,𝑠 are the values in 𝑓1 and 𝑓2 for the 𝑠th rithmperformance: execution, 𝑁𝑠 is the number of nondominated solutions derivedinthe𝑠thexecution,𝑅𝑠 istheredundancymeasure, (1)Minimumvalueforeachofobjectivefunctions𝑓1and and𝐻𝑠isthemeanhypervolumedominated,respectively. 𝑓2:althoughtheaimofamultiobjectiveoptimization Asecondmetricrepresentsthecomputationtime: algorithm is to find compromise solutions equidis- 𝑡𝑠 tantly spread around the Pareto front, this indicator 𝑒𝑗 𝑚 = . can give an idea of the exploitation capacity of the 2,𝑠 max𝑠(𝑡𝑒𝑠𝑗) (34) algorithm. (2)Numberofnondominatedsolutionsfound:asthepre- 4.ParameterSettingandPerformance vious indicator, it is not a good measure of the Analysiswith5Aircrafts algorithm performance, since many solutions could be found but which are not well spread around the Theresultsofseveralexecutionswithdifferentcombination Paretofront.However,itmayprovidearoughideaof ofparametersregardingbothmetricsareshowninFigure7, theoperatingcapabilityofthealgorithm. whereasFigure8showsnondominatedpoints(executions)in Figure7. (3)Dispersionofsolutions:thismeasureprovidesinfor- Twelve nondominated points were found. The combi- mation about the good or bad distribution of solu- nation of parameters we chose (pointed out in red) due to tionsaroundtheParetofront.However,itcouldoccur its very good performance and very low computation time that the algorithm derives a set of solutions which areverywellspreadaroundtheParetofrontbutthat (0.1156seconds)is𝑇0 =100,𝐿=5,𝛼=0.99,𝑞=0.4,𝜑=1/4, andconvergenceiterations=20. set includes very few solutions. Thus, this measure The valuations in the indicators for the selected com- shouldbeinterpretedtogetherwiththepreviousone. Similarly, solutions could be well distributed but far bination are 𝑓1,𝑠 = 0.1947, 𝑓2,𝑠 = −960.7386 (0.0012 normalized), number of nondominated solutions = 115.2 fromtheParetofront. Aredundancymeasure(𝑅)isused,whichconsistsof (0.0938 normalized), 𝑅𝑠 = 78.7474 (0.0629 normalized), takingthedifferencebetweentheaveragedistanceof 𝐻𝑠 = 60.9578(0.1343normalized),and𝑡𝑒𝑠𝑗 = 0.1156seconds eachsolutionanditsnearestone: (𝑚2,𝑠 = 0.00073),with3583.8iterationscarriedout(0.0064 normalized). 𝑅= 𝑝1∑𝑝 󵄨󵄨󵄨󵄨𝑚−𝑑𝑖,min󵄨󵄨󵄨󵄨, (32) tionWofeswomillenionwterpersotvinidgeedleimffeernetnstafingdutrheesdshiffoewreinngtitnhdeievvidoulua-l 𝑖=1 objectives along the algorithm execution with the selected