P Scala
email: p.m.scala@hva.nl

M Mujica

J Ma
email: ma@recherche.enac.fr

D Delahaye
email: delahaye@recherche.enac.fr

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 1

Keywords: Airport, Airspace, Optimization, Simulation, Uncertainty. I

Airport capacity has become a constraint in the air transportation networks, due to the growth of air traffic demand and the lack of resources able to accommodate this demand. This paper presents the algorithmic implementations of a decision support system for making a more efficient use of the airspace and ground capacity. The system would be able to provide support for air traffic controllers in handling large amount of flights while reducing to a minimum the potential conflicts. In this framework, airspace together with ground airport operations are considered. Conflicts are defined as separation minima violation between aircraft for what concerns airspace and runways, and as capacity overloads for taxiway network and terminals. The methodology proposed in this work consists of an iterative approach that couples optimization and simulation to find solutions that are resilient to perturbations due to the uncertainty present in different phases of the arrival and departure process. An optimization model was employed to find a (sub)optimal solution while a discrete event-based simulation model evaluated the objective function. By coupling simulation with optimization, we generate more robust solutions resilient to variability in the operations, this is supported by a case study of Paris Charles de Gaulle Airport.

are the main constraints for airport expansion. Air traffic growth rate for the next 20 years is forecasted as 4.4% [START_REF] Airbus | Global Market Forecast Growing Horizon 2017-2036[END_REF], this calls out for extra capacity in airports to accommodate the future traffic. Moreover, this growth requires solutions in the short-term horizon rather than the long-term one. Short-term solutions focus normally in improving airport operational efficiency. This improvement will result in an increase of the airport throughput and will reduce delays which is one of the most common indicator of inefficiencies.

The focus of this study is on air traffic management operations in the vicinity of airports, where airside and airspace operations are tackled in an integrated fashion. Airside operations are defined as operations carried out by aircraft that involve ground components like runways, taxiways and gates. Airspace operations are defined as those carried out by aircraft in the terminal manoeuvring area (TMA), a portion of the airspace which surrounds airports. In the literature regarding this research area, it is possible to find several studies about the improvement of these operations. In some cases researchers focused on individual operations [START_REF] Sama | Rolling horizon approach for aircraft scheduling in the terminal contro, area of busy airports[END_REF], [START_REF] Farhadi | Runway capacity management-an empirical study with application to Doha International Airport[END_REF], while other focused on the integration of them [START_REF] Guepet | Integration of aircraft ground movements and runway operations[END_REF]. In this context, many analytical models have been developed and applied to real case studies.

The problem studied in this work aims at providing an efficient decision support system (DSS) for managing the traffic and capacity of an airport. The objective is to avoid any type of conflict (e.g. separation minima violations, runway, taxiway, terminal capacity violations), and harmonize the sequence of arrivals and departures. Such a problem is defined by the authors as extended conflict detection and resolution (ECD&R), which is an extension of the air traffic conflict detection and resolution (CD&R) problem [START_REF] Kuchar | A review of conflict detection and resolution modeling methods[END_REF] considering also ground operations such as runways, taxiway and terminals operations. The main issue of most of the models about airport operations present in literature is that they do not consider the variability of real the system. Therefore, when solutions are applied to real conditions, they are very rigid and can become unfeasible due to disturbances of the real system. With the objective of making the solutions more resilient and applicable, the uncertainty coming from real-world conditions is considered in this approach. To this end, simulation techniques have been incorporated in the optimization algorithm to include the inherent variability of the real traffic. The purpose of the problem is to generate solutions which are Tackling Uncertainty for the Development of Efficient Decision Support System in Air Traffic Management Paolo Scala, Miguel Mujica, Ji Ma, and Daniel Delahaye M more resilient to perturbations due to uncertainty. The remaining of the paper is organized as follows: in section Ⅱ a literature review about optimization and simulation techniques applied to aviation-related problems is given, in section Ⅲ the problem tackled by this study is described. Section Ⅳ describes the methodology proposed in this study explaining the architecture of optimization and simulation and how it is applied for the specific problem. In section Ⅴ the implementation of the methodology on a real case study is presented. Conclusions and future work are presented in section Ⅵ.

II. LITERATURE REVIEW

This section presents a literature review about various techniques that include the use of optimization, simulation and a combination of the two, for problems that involve air traffic management operations at airports.

A. Optimization techniques

In the literature, one can find approaches that integrate arrival, departure and surface management as presented in [START_REF] Kjenstad | Optimizing AMAN-SMAN-DMAN at Hamburg and Arlanda Airport[END_REF]. In this work, it was implemented an optimization algorithm for solving in a sequential way the aircraft ground routing, arrival and departure sequencing, and conflict resolution problems. The ground aircraft routing problem was solved applying a shortest route algorithm, the arrival and departure schedule were solved with a time indexed formulation and the conflict resolution with a linear program model. The algorithm was tested with a real case study, achieving general improvements in surface performance indicators such as taxi times and punctuality. In [START_REF] Bertsimas | Unified optimization of traffic flows through airports[END_REF], near-terminal airspace and airside operations were unified in a two phase approach where, at the beginning, runway configuration selection and runway sequence problems were solved, and then, gate holding and ground routing decisions were made for airside operations. A binary optimization model was developed and results showed reductions in delays, fuel use and emissions. In [START_REF] Ma | Integrated optimization of terminal maneuvering area and airport[END_REF], airspace aircraft sequencing operations were integrated with airport ground operations at a macroscopic level. In this work, taxiways and terminals were modelled in low detail, considering them as nodes of a network characterized only by their capacity. The problem was solved by a heuristic employing a sliding window approach. The model was tested with a real case study and results showed the capability of reaching a conflict-free scenario for the sequencing of aircraft while mitigating ground congestion. In [START_REF] Guepet | Integration of aircraft ground movements and runway operations[END_REF] departure process management integrated with ground routing problem was proposed, with the aim of increasing runway efficiency and reducing taxi times. A mixed integer problem was formulated for both runway sequencing problem and ground routing problem. A sequential heuristic to solve both problems was employed. The methodology was able to provide high quality results in terms of completion times and taxi time while achieving short computational time. In all the work mentioned in this section, the solutions generated were based on optimization models using deterministic values, this represents a drawback since they will not be able to perform well in a real world scenario. In this work this drawback has been overcome by coupling simulation techniques to an optimization model. Furthermore, a limitation found in [START_REF] Bertsimas | Unified optimization of traffic flows through airports[END_REF], is linked to the fact that the operations were modeled in a deterministic way, by not addressing the uncertainty related to these kinds of operations. In this work, in order to overcome this limitation, we included different sources of uncertainty by coupling optimization techniques with simulation techniques. In the present work the optimization part is based on [START_REF] Ma | Integrated optimization of terminal maneuvering area and airport[END_REF], while and the simulation is based on discrete event simulation (DES) model where stochasticity is included to represent uncertainty from the real world.

B. Simulation techniques

Simulation is a technique widely used for many different types of problems in different fields. Its use is supported by the fact that it allows to include the variability and causality of real processes replicating with high detail the evolution of real systems. Simulation has been traditionally used for what-if analysis and the modeling of facilities during their planning phase or for infrastructures that are difficult or too expensive to experiment with in real life, like manufacturing plants, sea ports, airports or supply chain systems. In recent years, due to its capabilities for including variability and different abstraction levels, scientific community has done extensive research in techniques that couple simulation techniques with optimization ones with good results [START_REF] Amaran | Simulation optimization: a review of algorithms and applications[END_REF]. Within the air transportation framework, simulation has been widely used for many different types of operations involving airspace, airport airside and terminals. In [START_REF] Odoni | Existing and required modeling capabilities for evaluating ATM systems and concepts[END_REF] an extensive survey on the existing airport simulators at the time is presented. Although it is quite outdated, the main software that are described still exist, and they are still utilized to date. In this survey, it was made a classification between the airport models in macroscopic, mesoscopic and microscopic based on the level of detail and their purpose (strategical, tactical and operational). The most relevant airport software presented in the survey are SIMMOD [START_REF] Simmod | [END_REF] and RAMS [12], developed by the FAA and EUROCONTROL, respectively. SIMMOD is able to simulate airspace and airfield operations at airports, it also considers disruptive events, human activities as well as advanced operating concepts among other features. RAMS is a fast-time simulation software that focuses on airspace operations with the objective to quantify performance benefit for air traffic management (ATM) decision support. The majority of these programs are based on discrete event simulation paradigm [START_REF] Banks | Discrete-event System Simulation[END_REF], in which the state of the system changes based on specific event that are not time dependent. In [START_REF] Delcaire | Development of an on-site ground operations model for Logan International Airport[END_REF] was made an evaluation of SIMMOD simulation software for developing a model of the ground operations at Boston Logan Airport. The purpose of their study was to develop a model to be employed as a support decision tool for short term planning, where short-term planning considers operations that needs to be planned after alterations made by gate management or small disruptions. In [START_REF] Wei | Simulation study on closed spaced parallel runway analysis using SIMMOD plus[END_REF], SIMMOD was used as a tool to make an analysis of a closely spaced parallel runway, several scenarios were built based on the runway configurations assessing the overall performance according to delay and runway departure queue values. Beside simulation software used for specific purposes like SIMMOD and RAMS, there are also other types of simulation programs that for general purpose, that can be used for different application in different fields. These simulation software allow analysts to be more flexible when building the logic of the system, contrary to specific simulation software, where analysts are restricted to use the predefined objects and logic. In [START_REF] Scala | A divide and conquer approach for simulating an airport system[END_REF] a DES model for the evaluation of the integrated operations of both airspace and airside was proposed, by combining four different modules, each of them representing a specific operation (airspace aircraft sequencing, turn around, taxiway routing, runway movements). In this way an analysis of the airport performance from a holistic view was realized, identifying potential problems that might affect the capacity of the airport. This approach allows to run each module as a stand-alone model for conducting a more specific analysis, narrowing the area of interest to a specific area of the airport. This simulation model was developed with a general purpose simulation software called SIMIO [17]. SIMIO, like other general purpose simulation software present on the market, provides an extendable library of objects useful for modeling a wide range of dynamic systems. This program is used in the current study of the paper. In this work, a simulation model was developed, similarly as [START_REF] Scala | A divide and conquer approach for simulating an airport system[END_REF], however, in [START_REF] Scala | A divide and conquer approach for simulating an airport system[END_REF], the conflict detection and resolution (ECD&R) was not implemented, as it represents an innovative feature implemented in this work. An additional feature in which this work goes beyond the state of the art, is the inclusion of uncertainty in order to represent the variability of the operations both in the airspace and on the ground.

C. Simulation combined with optimization

Simulation together with optimization has proven to be a good combination of techniques for tackling problems were uncertainty plays a crucial role. Simulation on the one hand, allows representing the uncertainty that affects the system under study, which is not possible to model by using analytical models. On the other hand, optimization gives an optimal or near-optimal solution even for complex systems. By coupling the two techniques, a solution that is more feasible for implementation in real life can be found. In [18] examples of applications of simulation and optimization techniques for problems that include aeronautical, logistics and industrial cases study are presented. Within airport operations, uncertainty is a crucial factor and therefore it should always be considered when making an analysis. Reference [START_REF] Odoni | The importance of probability theory in the airport and air traffic control sectors[END_REF] highlights the importance of considering uncertainty in airport and air traffic control operations and the importance of including probabilistic aspects. It also mentions the use of computer models as a tool for the evaluation of such systems. In the literature we can find many applications that employ optimization combined with simulation techniques to solve specific airport operations: for instance [START_REF] Kleinman | Simulation optimization of air traffic delay cost[END_REF] deals with the problem of reducing the cost of delay in air traffic; [START_REF] Arias | A methodology combining optimization and simulation for real applications of the stochastic aircraft recovery problem[END_REF] focuses on stochastic aircraft recovery problem; [START_REF] Lee | Airport surface traffic optimization and simulation in the presence of uncertainties[END_REF] considers airport surface traffic; and [START_REF] Mujica | Check-in allocation improvements through the use of a simulation-optimization approach[END_REF] copes with the check-in allocation at airport terminals. Regarding similar problems as the one faced in this work, in [START_REF] Mujica | Improving airport performance through a model-based analysis and optimization approach[END_REF] a DES model integrating TMA and airside operations has been developed. Furthermore, they integrated the optimization tool OptQuest [25] in order to find the best combination for the number of ground vehicles with the objective of minimizing the turnaround time. To the best of our knowledge there are no studies that tackled the aforementioned ECD&R problem employing simulation together with optimization, therefore the purpose of this work is to fill this gap in the literature. In most of the previous references, simulation is used in order to include uncertainty and evaluating the objective function of an optimization model, our work differs from them in the sense that a better solution is found by applying an iterative process where some input parameter of the optimization are tuned after the optimized solution has been evaluated by using simulation. The paper presents a novel methodology which introduces for the first time simulation and optimization in a coordinated fashion. We present throughout the paper all the details on how to tackle real problems using this approach and we finally support it by presenting a case study based on real data from Paris Charles da Gaulle Airport.

III. PROBLEM DESCRIPTION: INTEGRATING AIRSPACE AND AIRSIDE OPERATIONS

The problem tackled in this work, consists in the integration of airspace and airside operations. Airspace operations involve landing sequencing while airside operations involve runway, taxiway and terminal operations. The objective, regarding the landing sequencing, is to resolve airspace conflicts and have a smooth arrival flow of aircraft. Regarding airside operations, the objective is to regulate the departure rate aiming at mitigating the congestion at the airport surface. The overall objective is to produce a solution that balances both systems in terms of conflicts. In Fig. 1, a schematic representation of the main components modeled is given, showing how these components are linked to each other. In this section the airspace and airside operations which have been modeled are explained, giving the necessary details to the reader for understanding the level of abstraction that has been applied. For further details about the problem definition please refer to [START_REF] Ma | Integrated optimization of terminal maneuvering area and airport[END_REF].

A. Airspace operations

The airspace operations considered in this work concern the sequencing of landing aircraft inside the TMA airport airspace. The TMA is a portion of airspace that surrounds airports, its size is usually between 30 and 40 nautical miles radius, and its center is located at the airport location as it is represented in Fig. 2. According to [START_REF] Oag | OAG punctuality league[END_REF], airports are classified by the number of departing seats transported per year. Especially for the ones that transport more than 30 million of seats, classified as mega airports, TMA can be a very congested zone due to the converging traffic flow from and to different origin/destination airports. For that reason, the management of the aircraft flow in the TMA, ensured by the work of the air traffic controllers (ATCOs), can be very demanding. In literature, there are many studies aiming at supporting ATCOs work by the development of optimal arrival schedules [START_REF] Sama | Rolling horizon approach for aircraft scheduling in the terminal contro, area of busy airports[END_REF], [START_REF] Beasley | Scheduling aircraft landings: the static case[END_REF], [START_REF] Michelin | Merging of air traffic flows[END_REF], [START_REF] Zuniga | Integrating and sequencing flows in terminal maneuvering area by evolutionary algorithms[END_REF], [START_REF] Furini | Improved rolling approaches to the aircraft sequencing problem[END_REF], [START_REF] Toratani | Simultaneous optimization method for trajectory and sequence for receding horizon guidance in terminal area[END_REF].

In this work, the optimal scheduling of the aircraft during their landing phase is ensured by making decisions regarding: entry time in the TMA, initial speed in the TMA and runway choice. The objective is to avoid conflicts derived from the violation of separation minima between aircraft.

1) Airspace conflicts detection

Aircraft incur conflicts every time a loss of separation minima between two consecutive aircraft is detected. ICAO [START_REF]Air Traffic Management[END_REF] provides values regarding separation minima based on the aircraft wake turbulence category. The wake turbulence category is an indicator related to the aircraft maximum certificated take-off mass, and it classifies aircraft in three categories: light, medium, and heavy. In Table Ⅰ the values about separation minima are shown. Formalizing: let be A a given set of landing aircraft, i and j consecutive pair of aircraft, 𝑆𝑆 𝑖𝑖𝑖𝑖 the separation minima between the leading aircraft type i and the trailing aircraft type j and 𝐷𝐷 𝑖𝑖𝑖𝑖 the distance between aircraft type i and aircraft type j. In this way conflicts are detected and calculated as:

𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴 𝐴𝐴𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝐴𝐴𝐴𝐴𝑐𝑐 = � 1, 𝐷𝐷 𝑖𝑖𝑖𝑖 < 𝑆𝑆 𝑖𝑖𝑖𝑖 , ∀𝐴𝐴, 𝑗𝑗 ∈ 𝐴𝐴 0, 𝑐𝑐𝑐𝑐ℎ𝐴𝐴𝐴𝐴𝑒𝑒𝐴𝐴𝐴𝐴𝐴𝐴 (1)
Another type of conflict is detected for aircraft that do not respect the order of the sequence along the landing route.

In this context, this type of conflict is detected by checking at the entry and exit of each link which belong to a route, that the order of the entry sequence is equal to the order of the exit sequence. Formalizing: let be 𝑐𝑐(𝑢𝑢, 𝑣𝑣), 𝑐𝑐 = 1, . . , 𝑐𝑐, a given set of links 𝑐𝑐𝐴𝐴𝑜𝑜𝐴𝐴𝐴𝐴 𝑓𝑓𝑓𝑓 𝑙𝑙 𝐴𝐴𝑐𝑐𝑜𝑜 𝑐𝑐𝐴𝐴𝑜𝑜𝐴𝐴𝐴𝐴 𝑓𝑓𝑓𝑓 𝑙𝑙 the positions of the aircraft 𝐴𝐴 on the link entry 𝑢𝑢 and on the link exit 𝑣𝑣, respectively. The order of sequence conflicts are calculated as:

𝑂𝑂𝐴𝐴𝑜𝑜𝐴𝐴𝐴𝐴 𝑐𝑐𝑐𝑐 𝑐𝑐ℎ𝐴𝐴 𝐴𝐴𝐴𝐴𝑠𝑠𝑢𝑢𝐴𝐴𝑐𝑐𝐴𝐴𝐴𝐴 𝐴𝐴𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝐴𝐴𝐴𝐴𝑐𝑐 = ∑ (∑ 𝑐𝑐𝐴𝐴𝑜𝑜𝐴𝐴𝐴𝐴 𝑖𝑖𝑓𝑓 𝑙𝑙 -𝑐𝑐𝐴𝐴𝑜𝑜𝐴𝐴𝐴𝐴 𝑖𝑖𝑓𝑓 𝑙𝑙 𝑖𝑖∈𝐴𝐴 ) 𝑛𝑛 𝑙𝑙=1
(2)

B. Airside operations

According to [START_REF] Ma | Integrated optimization of terminal maneuvering area and airport[END_REF] decisions about landing operations can be made in a time frame of two hours in advance, while keeping the system with a certain level of detail and preserving the system from any loss of feasibility. In such a time frame, due to the uncertainty present in the airspace it is not relevant to give instructions for airside operations such as taxiway routing logics, and that is why, in this study, airside operations were modeled following a so-called macroscopic approach. In this context, a macroscopic approach for airside operations is achieved by modeling airside components in a high level, which means, by identifying the main airside components such as taxiways and terminals and considering them as nodes of a network characterized by their own capacity. In this study a capacity for the taxiway network is defined as the maximum number of aircraft that can simultaneously cross it without incurring surface congestion; capacity of the terminals is defined by the number of available parking gates. A factor which is considered as a decisions variable for airside operations is the choice of the push back time, that is, the time when an aircraft receives clearance to leave from its gate. Following this approach, the objective for airside operations is to keep the level of congestion under a certain threshold, which is the taxiway and terminal declared capacity.

1) Airside conflicts detection

Airside conflicts are detected, for taxiway and terminal components, when their declared capacity are exceeded (capacity overload). Conflicts values are derived from the weighted sum of two factors: the maximum overload (max OL) and the average overload (avg OL). Formalizing: let be C the declared capacity, 𝑂𝑂 𝑡𝑡 the aircraft occupancy for each discrete time increment 𝑐𝑐 ∈ 𝑇𝑇; and T the entire time frame considered, the max OL is given by:

𝑀𝑀𝐴𝐴𝑀𝑀 𝑂𝑂𝑂𝑂 = � 𝑚𝑚𝐴𝐴𝑀𝑀 𝑡𝑡∈𝑇𝑇 (𝑂𝑂𝑣𝑣𝐴𝐴𝐴𝐴𝑐𝑐𝑐𝑐𝐴𝐴𝑜𝑜 𝑡𝑡 ), 𝐴𝐴𝑐𝑐 𝑂𝑂 𝑡𝑡 > 𝐶𝐶 0, 𝑐𝑐𝑐𝑐ℎ𝐴𝐴𝐴𝐴𝑒𝑒𝐴𝐴𝐴𝐴𝐴𝐴 (3)
and avg OL is given by:

𝐴𝐴𝑣𝑣𝐴𝐴 𝑂𝑂𝑂𝑂 = ∑ 𝑜𝑜𝑓𝑓𝑜𝑜𝑜𝑜𝑙𝑙𝑜𝑜𝑜𝑜𝑜𝑜 𝑡𝑡 𝑡𝑡∈𝑇𝑇 𝑇𝑇 (4)
where,

𝑐𝑐𝑣𝑣𝐴𝐴𝐴𝐴𝑐𝑐𝑐𝑐𝐴𝐴𝑜𝑜 𝑡𝑡 = � 𝑂𝑂 𝑡𝑡 -𝐶𝐶, 𝐴𝐴𝑐𝑐 𝑂𝑂 𝑡𝑡 > 𝐶𝐶 ∀ 𝑐𝑐 ∈ 𝑇𝑇 0, 𝑐𝑐𝑐𝑐ℎ𝐴𝐴𝐴𝐴𝑒𝑒𝐴𝐴𝐴𝐴𝐴𝐴 (5)
the total conflicts value for each airside component is:

𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝑜𝑜𝐴𝐴 𝐴𝐴𝑐𝑐𝑚𝑚𝐴𝐴𝑐𝑐𝑐𝑐𝐴𝐴𝑐𝑐𝑐𝑐 𝐴𝐴𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝐴𝐴𝐴𝐴𝑐𝑐𝐴𝐴 = 𝛼𝛼𝑀𝑀𝐴𝐴𝑀𝑀 𝑂𝑂𝑂𝑂 + 𝛽𝛽𝐴𝐴𝑣𝑣𝐴𝐴 𝑂𝑂𝑂𝑂(6)
The weights 𝛼𝛼 and 𝛽𝛽 can guide the optimization minimization objective function towards different directions since they can prioritize the reduction of the max OL or of the avg OL. In this work these two weights have the same value, but in the future they can be exploited for further investigations.

C. Runway operations

Runway operations act as the connections between airspace and airside operations. The runway system is considered as the bottleneck of the airport system [START_REF] Idris | Observations of departure processes at Logan Airport to support the development of departure planning tools[END_REF], so a crucial aspect to improve is the performance of the runway system in terms of rate (arrivals and departures). In this work, landing and departure operations at runway are considered, so an additional objective of this study is to find a balance in the rate of arrivals and departures.

1) Runway conflicts detection

Runway conflicts are defined in a similar way as conflicts for the airspace. The only difference is in the values of separation minima which depend on the type of operation on the runway and also on the aircraft type. In Table Ⅱ, values of separation minima for the different runway operations (A-arrivals, D-departures) are presented [START_REF] Frankovich | Air traffic flow management at airports: a unified optimization approach[END_REF]. Formalizing: let be A a given set of aircraft, i and j a consecutive pair of aircraft, 𝑆𝑆 𝑖𝑖𝑖𝑖 the separation minima between aircraft type i and aircraft type j and 𝐷𝐷 𝑖𝑖𝑖𝑖 the detected time difference between aircraft type i and aircraft type j. In this way conflicts are formalized with the following function:

𝑅𝑅𝑢𝑢𝑐𝑐𝑒𝑒𝐴𝐴𝑅𝑅 𝐴𝐴𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝐴𝐴𝐴𝐴𝑐𝑐 = � 1, 𝐴𝐴𝑐𝑐 𝐷𝐷 𝑖𝑖𝑖𝑖 < 𝑆𝑆 𝑖𝑖𝑖𝑖 , ∀𝐴𝐴, 𝑗𝑗 ∈ 𝐴𝐴 0, 𝑐𝑐𝑐𝑐ℎ𝐴𝐴𝐴𝐴𝑒𝑒𝐴𝐴𝐴𝐴𝐴𝐴 (7)

D. Model network

Airspace and airside operations are modeled by means of a network of nodes and links. Landing routes in the airspace are identified by the connection of nodes called waypoints forming a trajectory, Fig. 3 shows an example of landing routes or standard instrument arrival routes (STARs) taken from the case study considered in this work. Aircraft fly with a constant deceleration from one of the STAR entry points, in order to slow down until they reach a waypoint called final approach fix (FAF), from there, they continue descending with a constant speed until they approach the runway. In this work, it was possible to modify the speed at the entry point of the STAR. The routes considered in this work are assumed to be vertically well spaced, therefore, there are no possibilities of conflicts due to vertical separation minima violation. Another assumption regards to the departure routes, which in this work are omitted, therefore, they do not interfere with the landing routes. Airspace conflicts are detected on the nodes and on the links [START_REF] Ma | Integrated optimization of terminal maneuvering area and airport[END_REF], [START_REF] Zuniga | Integrating and sequencing flows in terminal maneuvering area by evolutionary algorithms[END_REF]. Airside components such as runway, taxiway and terminals are modeled as nodes inside the network. In the case of the runway, conflicts are detected in the node using formula [START_REF] Bertsimas | Unified optimization of traffic flows through airports[END_REF]. For the airside components, each node is characterized by a capacity value, and conflicts are detected according to the maximum number of aircraft on the nodes. and maximum values that can be assigned to 𝑐𝑐 𝑓𝑓 , and ∆𝑇𝑇 is a discretized time increment. In this case ∆𝑇𝑇 𝑚𝑚𝑖𝑖𝑛𝑛 is -5 minutes, ∆𝑇𝑇 𝑚𝑚𝑜𝑜𝑚𝑚 is +30 minutes and ∆𝑇𝑇 is 5 seconds. , 𝑗𝑗 ∈ ℤ }, ∆𝑃𝑃𝑃𝑃 𝑚𝑚𝑜𝑜𝑚𝑚 is the maximum value that can be assigned to 𝐴𝐴𝑝𝑝 𝑖𝑖 , and ∆𝑇𝑇 is a discretized time increment. In this case ∆𝑃𝑃𝑃𝑃 𝑚𝑚𝑜𝑜𝑚𝑚 is +15 minutes and ∆𝑇𝑇 is 5 seconds. The first decision can assign a delay or an advance to the aircraft which is about to enter inside the TMA. The second decision modifies the entry speed of the aircraft that is about to enter one of the STARs in the TMA, in this case aircraft can enter with a lower or higher speed. The third decision chooses the landing runway. The last decision is the push back time, in this case it is possible only to assign a delay, establishing the time when the aircraft can leave the gate.

The objective function evaluated in this problem is composed by the weighted sum of two factors: airspace conflicts and airside conflicts. It is important to mention that runway conflicts are included in airside conflicts computation. 𝑂𝑂𝑝𝑝𝑗𝑗𝐴𝐴𝐴𝐴𝑐𝑐𝐴𝐴𝑣𝑣𝐴𝐴 𝑐𝑐𝑢𝑢𝑐𝑐𝐴𝐴𝑐𝑐𝐴𝐴𝑐𝑐𝑐𝑐 = 𝛾𝛾 𝑜𝑜𝑖𝑖𝑜𝑜𝑎𝑎𝑎𝑎𝑜𝑜𝑎𝑎𝑜𝑜 * (𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴 𝐴𝐴𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝐴𝐴𝐴𝐴𝑐𝑐𝐴𝐴 + 𝑂𝑂𝐴𝐴𝑜𝑜𝐴𝐴𝐴𝐴 𝑐𝑐𝑐𝑐 𝑐𝑐ℎ𝐴𝐴 𝐴𝐴𝐴𝐴𝑠𝑠𝑢𝑢𝐴𝐴𝑐𝑐𝐴𝐴𝐴𝐴 𝐴𝐴𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝐴𝐴𝐴𝐴𝑐𝑐𝐴𝐴) + 𝛾𝛾 𝑜𝑜𝑖𝑖𝑜𝑜𝑎𝑎𝑖𝑖𝑜𝑜𝑜𝑜 * (𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝑜𝑜𝐴𝐴 𝐴𝐴𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝐴𝐴𝐴𝐴𝑐𝑐𝐴𝐴 + 𝑅𝑅𝑢𝑢𝑐𝑐𝑒𝑒𝐴𝐴𝑅𝑅 𝐴𝐴𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝐴𝐴𝐴𝐴𝑐𝑐𝐴𝐴) (8) The weight coefficients 𝛾𝛾 𝑜𝑜𝑖𝑖𝑜𝑜𝑎𝑎𝑎𝑎𝑜𝑜𝑜𝑜𝑎𝑎𝑜𝑜 and 𝛾𝛾 𝑜𝑜𝑖𝑖𝑜𝑜𝑎𝑎𝑖𝑖𝑜𝑜𝑜𝑜 , are assigned to the airspace and airside conflicts, respectively. Depending on which component the optimization process should focus on, these weights can be adjusted accordingly. In this work, for both 𝛾𝛾 𝑜𝑜𝑖𝑖𝑜𝑜𝑎𝑎𝑎𝑎𝑜𝑜𝑎𝑎𝑜𝑜 and 𝛾𝛾 𝑜𝑜𝑖𝑖𝑜𝑜𝑎𝑎𝑖𝑖𝑜𝑜𝑜𝑜 , we assigned 1 as their values.

F. Data input

The data input used in both optimization and simulation model is similar to an airport flight schedule. The data available is taken from real flight data and it covers an entire day of operations. The optimization model modifies this flight schedule according to the value of its decision variables provided by the (sub)optimal solution. The data input contains the following information: Ⅲ contains an example of data input from the flight schedule, if we consider for instance, the flight corresponding to the first row, we see that this flight is of type arrivaldeparture meaning that the same aircraft will first land, and then, after turnaround operations will depart. Looking at the data input it is known that this aircraft will enter the TMA from the TMA entry point number 3 at 10:30:20 as it is denoted by the entry point in the TMA and the entry time in the TMA values. It will fly at a speed of 430 Kt as the entry speed in the TMA value shows, and will land at the landing runway number 2. Once landed, the aircraft will reach terminal number 2 and then leave the terminal at 12:10:30 as it is set by the terminal number and the push back time values present in the flight schedule. Finally, it will take off from the departing runway number 3 as it is set in the flight schedule. Values for entry time in the TMA, entry speed in the TMA, landing runway and push back time are modified accordingly to the solution provided by the optimization model. Looking at the second row, it shows a flight of type arrival, meaning that the aircraft will land at the runway, taxi to the terminal and then occupying a gate parking position for the entire day. It can be noticed that for this flight the records corresponding to the push back time and departing runway are empty, because this aircraft will not depart.

IV. METHODOLOGY DESCRIPTION

Optimization and simulation techniques are combined in this study in order to develop a DSS for ATM operations efficiency. Simulation allows to include the uncertainty present in the ATM framework, and supported by an optimization algorithm, the combination of both techniques provides more robust and feasible solutions. Simulation is used for evaluating the optimized solution when uncertainty is also considered, the methodology aims at finding a more robust solution by tuning the input parameters of the optimization. This procedure is furthermore iterated until a better solution is found. Fig. 4 shows a schematic representation of this methodology, where firstly, an optimization model is run, providing an (sub)optimal solution evaluated by its objective function value Z_opt. Then a simulation model is employed in order to evaluate the same objective function but including uncertainty in the problem, in this case, the output is the objective function Z_sim value. At this stage, the robustness of the solution is measured by the value of Z_sim. In this study, we define robustness as the capability of the system of being resilient to perturbations. In this context, perturbations are defined as sources of variability that affect the airport operations, such as flight inefficiencies in the airspace and on the ground due to human factor miscommunication, weather conditions, and so on. These perturbations, which represent the uncertainty of the system, are implemented with the simulation in order to evaluate the objective function and measure the robustness of the model. Since the optimization always provides a conflictfree solution, the absolute value of the objective function evaluated by the simulation, Z_sim, is taken as indicator of the robustness. A small value corresponds to a robust solution, while a big value means that the solution is less robust since it is away from the optimal solution.

A. Optimization approach

The optimization approach is based on the implementation of the sliding window method. This method has been employed in other studies with good results [START_REF] Ma | Integrated optimization of terminal maneuvering area and airport[END_REF], [START_REF] Furini | Improved rolling approaches to the aircraft sequencing problem[END_REF], [START_REF] Toratani | Simultaneous optimization method for trajectory and sequence for receding horizon guidance in terminal area[END_REF], [START_REF] Hu | Receding horizon control for aircraft arrival sequencing and scheduling[END_REF], [START_REF] Hu | Multiairport capacity management: genetic algorithm with receding horizon[END_REF], [START_REF] Zhan | An efficient ant colony system based on receding horizon control for the aircraft arrival sequencing and scheduling problem[END_REF], [START_REF] Liang | A novel approach to automated merge 4D arrival trajectories for multi-parallel runways[END_REF]. The sliding window approach consists in considering time frame windows of small size and shifting them along the entire time horizon of the study. By using this approach, a small instance of the problem is solved considering only the elements within the time window. Fig. 4 illustrates the approach. The main parameters to set are the window length and the shift length. The main advantages coming from the use of this method are: less computational time required; the possibility of treating the problem in a dynamic way, by considering, as time passes, new information updates due to changes in the environment and to the interactions between entities with the surrounding environment. As highlighted by [START_REF] Ma | Aircraft merging and sequencing problems in TMA[END_REF] it is crucial to choose the right values for the window and shift length, since big values of window length might lead to worsen the overall performance, as decisions that are too far ahead in time are made. Moreover, there will also be an increase in the computational time. On the other hand, small values of window length might lead not to consider enough information about aircraft, affecting the quality of the solution. In the same way, a big window shift will reduce the effect of the dynamic behavior of the method, since in this case less decisions will be updated from one window to another. Small window shift will increase the computational time performance. As it can be appreciated from the Fig. 5, depending on the relative time position of the aircraft entities with respect to the current active window, there can be different categories of decisions: fixed, updated, to make or to be made. • Fixed decisions are the ones that fall behind the window and it is not anymore possible to take decisions for them. • Updated decisions are the ones that fall partly before and partly inside the window, these decisions were already made in the previous window and they will be updated in the current window. • To make decisions are the ones that are fully inside the widows and these decisions will be made for them in the current window. • To be made are the ones that are after the windows and these decisions are not yet made for them.

The capability of the method to update decisions ensures its dynamic behavior in solving the problem. This type of problem is proven to be NP-hard according to [START_REF] Oag | OAG punctuality league[END_REF], therefore a meta-heuristic has been applied to solve the problem and to find a (sub)optimal solution. The metaheuristic employed is the simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. The metaheuristic will solve any instance of the problem in the boundaries of each window. This approach has been successfully applied to similar problems like the one presented in [START_REF] Liang | A novel approach to automated merge 4D arrival trajectories for multi-parallel runways[END_REF]. For further details about the development of the optimization model mentioned in this work we refer the reader to [START_REF] Ma | Integrated optimization of terminal maneuvering area and airport[END_REF].

B. Simulation approach

The simulation model for the airspace and airside operations has been developed with the use of DES. DES is a simulation paradigm which fits well for the modeling of systems in which the entities have already a specific set of processes to undergo. In DES the state of the system changes based on the occurrences of events and are not depending on time occurrences. In this specific problem, aircraft entities undergo a predefined set of processes in a concurrent fashion, for this reason, the DES paradigm was used to model the processes of arrival, landing, taxi and turnaround in an airport.

The main advantage of using simulation in this work, is that it includes sources of uncertainties in the system. Uncertainty can derive from different causes like external and internal factors such as the surrounding environment or the human factor. Often, the mere use of analytical models does not allow to represent these uncertainty, thus, to overcome this limit, simulation can be used as a tool for including uncertainty. In the specific problem tackled in this work, uncertainty refers to: • Entry time in the TMA, it can be affected by deviations due to speed, environmental variations and human factors.

• Taxiway time, it can vary based the speed of the aircraft on the taxiway, that in turn, it can be affected by the taxiway network configuration and by the level of congestion on the taxiway. • Pushback time for the aircraft, it can vary due to ground handler inefficiencies in performing the necessary maneuvers to the aircraft. Despite the fact that turnaround operations can be very fluctuant in terms of time consumed, uncertainty was not considered for them, and this represents the main assumption of the simulation model. Turnaround times have been derived from a real schedule, which is why uncertainty was not considered. The following subsection describe the conceptual modeling of the elements in the system under study.

1) Modeling airspace network

The airspace in the simulation model has been modeled by a network of nodes and links. Each node represents some of the main waypoints that aircraft follow while flying the landing route. Since uncertainty plays a crucial role in the system, in the airspace network is important to detect conflicts as accurately as possible. In the simulation model, unlike the optimization model, airspace conflicts are detected only on the nodes of the network. In this context, the number of nodes used for building the routes of the airspace network plays an important role because there is the potential risk of missing some conflicts on the links. Fig. 6 shows an example that compares two different ways of modeling the airspace routes by using different number of nodes. As it can be seen the comparison between the route with less nodes, placed at the bottom in Fig. 6, and the other one with more nodes, placed in the upper part of Fig. 6, illustrates that some conflicts might not be detected when the route has few nodes. To overcome this problem, we divided each route in links of 5 NM, which according to [START_REF] Kleinman | Simulation optimization of air traffic delay cost[END_REF] seems to be a good distance for achieving a good level of accuracy. An additional feature of the simulation model is that it is capable of evaluating the number of aircraft involved in the conflicts, which is an indicator that provides more insights about the level of congestion of the airspace.

V. IMPLEMENTATION

As a proof of concept of the approach, we applied it to a real case study. The case of study refers to one of the major airport with the biggest traffic in Europe, so to properly test the benefits of the methodology. The optimization and simulation models were run on a 2.5 GHz core i7 CPU, under Windows operating system PC. The optimization model was developed based on Java code, while the simulation model was developed by using a general purpose simulation software called SIMIO. The average computational time was similar for each of the different scenarios, and it was 20 minutes for the optimization model and 5 minutes for the simulation model. These computational times make the approach suitable for an ATC tactical decision making. Tactical decisions are the ones taken on the day of operations and include the update of daily plans and the continuous capacity optimization according to real traffic demand [START_REF]Strategic, pre-tactical, tactical and post-ops air traffic flow and capacity management[END_REF]. In the next sections the case of study is introduced and results for the different scenarios are presented.

A. Case study: Paris Charles de Gaulle Airport

In this work, Paris Charles de Gaulle Airport (PCDG) was taken as a case study. PCDG Airport is one of the major airports in Europe, for its size, passengers transported and air traffic movements. It is the hub of the French carrier Air France, and according to EUROCONTROL [START_REF]Industry monitor[END_REF] in the 2016 it transported 66M of passengers and carried 473K air traffic movements. The data used in the experiments refers to a flight schedule of an entire day of operations at the airport. In total there were 1116 movements. 1) Airport characteristics Airport TMA and runway system. PCDG Airport has 4 parallel runways. They operate as independent runways, in the sense that they can simultaneously accommodate air traffic movements. In real operations two of the runways are used only for landings (runway in) and other two only for departures (runway out). In the TMA four landing routes were identified, two coming from south and two coming from north. These routes were identified based on airspace surveillance records. Depending on which of the two landing runway is used, in any of the four available landing routes, it is possible to take two different route trajectories. These trajectories are spaced by a safe vertical distance, so to have in total 8 trajectories available. Airport airside. PCDG Airport is constituted by three terminals and a complex taxiway network. Following the macroscopic approach for modeling the airport airside, only capacity has been considered for terminals and taxiway network. One assumption was made for ground operations regarding taxiway routes potential aircraft conflicts due to physical position on the taxiway were not detected. Another assumption was also that aircraft could park at any gate of the terminal, regardless of the aircraft size. Tables Ⅴ andⅥ show the capacity of these components and the average times taxiway were used in the optimization model. In Fig. 7 it is possible to see a schematic representation of the model and the methodology applied to it.

B. Analysis of the original and the first optimized schedule

The conflict detection and resolution problem applied to airport operation such as landings and ground components has been solved by means of an optimization model. This model was based on a heuristic algorithm which provides a (sub)optimal solution. The objective of the model was to minimize an objective composed by the sum of the conflicts for the airspace and ground operations. Section Ⅲ gives a detailed description of the problem. Table Ⅶ shows the number of conflicts detected when the original flight schedule is run without implementing the heuristic algorithm, and the total value of objective function is 1972.76. Most of the conflicts were detected in the airspace, where we had 1114 conflicts on links and 461 conflicts on nodes. Runway was also affected, having 381 conflicts. Taxiway network had an overload peak of 8 and an average overload of 0.14, which shows that the system is always at the edge its capacity. Regarding the terminals, they had overload peaks of 2, 4 and 3 for terminal 1, 2 and 3, respectively. Average overload were 0.12, 0.3 and 0.2 for terminal 1, 2 and 3, respectively. From these results it can be seen that terminal 2 was the most congested in terms of both peak and average overload, although terminal 1 and 2 suffered from congestion as well. Overall, the results revealed that the airspace was more congested than the ground, meaning that the capacity of the ground was able to absorb the traffic better than the capacity of the airspace. Conflicts on the airspace routes (links and nodes) were transferred on the runway, which was also highly overloaded. Therefore, the optimization process needs to pay a lot of its effort in solving airspace conflicts, being a critical area to focus on. As first step in the methodology, the optimization model effectively finds a solution that reduces the conflict to zero, therefore it is considered as optimal. The solution obtained by the optimization model, is tested with the use of a simulation model where uncertainty is considered with their values showed in Table Ⅷ. Values have been chosen based on the experience of the authors, and they are continuous. Table Ⅸ shows the results obtained by the simulation model. The simulation model has been able to evaluate the conflicts for all the airspace and airside components, moreover, it has been able to evaluate further indicators like the number of aircraft involved in the airspace conflicts. Regarding the conflicts on the runway, it was possible to break them down in landing runway (runway in) and departing runway (runway out) conflicts. As can be seen from Table Ⅸ, when including uncertainty in the system, we detect some conflicts, especially for the airspace and the runway out components of the objective function. In Table Ⅸ are also shown minimum, average, maximum and standard deviation values since there were run 50 replications for each scenario. Looking at the average values, it can be noticed that terminals are not affected by uncertainty since they do not have any overload, while the taxiway network is slightly affected with max OL of 2.1 and an avg OL of 0.0057. This can be explained by the fact that the turnaround time, which ranges approximately between 30 minutes and 2 hours, absorbs the deviation of taxi times and pushback times, reducing at minimum the effect of the uncertainty. The component which is the most affected is the runway out with 90.44 conflicts, it is followed by the airspace with 19.18 conflicts and then runway in with 9.2 conflicts. These initial experiments pointed out that using optimization techniques alone to deal with such problems would not be enough for supporting decision makers in real operations. The deterministic nature of the optimization model represents a relevant limitation, as it was highlighted by the number of conflicts found after including uncertainty in the system by the simulation model.

C. Additional scenarios

From the previous analysis simulation raises the impact of uncertainty. With this in mind, the objective of the methodology is to improve the solution of the optimization and make it more resilient to uncertainty. One way of achieving this objective is by fine-tuning some of the optimization input parameters and obtaining a new solution that in turn can be tested with simulation for evaluating the performance. Following this line, a series of tests have been conducted based on the different values for the chosen input parameters to be tuned. The parameters chosen to be tuned are the separation minima between aircraft in the airspace and the one on the runways. By enlarging the separation minima in the optimization model we relax this constraint of the problem. With those new parameters, the optimization will achieve a more conservative solution, however, it might reduce the effect of the uncertainty on the number of conflicts because of the separation minima violation. Table Ⅹ shows the different value that have been tested for the input parameter. An indicator for measuring the goodness of a solution has been defined as the magnitude of the conflict detected. In this analysis, indicators like the average conflict size violation (%) is measured. Moreover, also the number of conflicts falling within a certain interval of separation violation percentage is measured and showed in a graph. With these two additional indicators, we conducted a more exhaustive analysis. In the next sections this analysis is presented individually for each objective function component: airspace, runway in, runway out, taxiway network and terminals conflicts. The charts depicted in Figs. 9(a) and 9(b), give information about the size of the conflicts. Fig. 9(a) in particular, shows the number of conflicts that fell in the different conflict size percentage intervals. Ideally a solution would be better if the conflicts size would be small, meaning to have high values on the left-hand side of the chart. In the chart it can be seen that the base case scenario is outperformed by all the other scenarios, and the best ones are the 20% and 30% separation minima increase scenarios. However, the new solutions obtained by the other scenarios produce more conflicts with high conflict size percentage than the base case scenario, and this is visible in Fig. 9 2) Runway in conflict analysis Fig. 10, shows that scenarios with 20% and 30% of separation minima increase outperform the base case and 10% separation minima increase scenarios, obtaining in average a value of 2 conflicts while the base scenario and 10% scenario have 9 and 6 conflicts, respectively. Concerning the conflict size we clearly see in both Figs.

11(a) and 11(b) that scenarios with 20% and 30% outperform the other two scenarios. Looking at Fig. 11(a), it is possible to find some conflicts with a high conflict percentage such as 25-30%, 30-40% and 75-100%. These circumstances would imply that the air traffic controllers would apply the same procedures as they would do for the airspace conflicts, performing missed approaches in the cases when aircraft would be very close to approach the runway [START_REF]Strategic, pre-tactical, tactical and post-ops air traffic flow and capacity management[END_REF]. Fig. 11(b) shows an improvement of around 1.5% from the base and 10% separation minima increase scenarios, with average conflict size of around 9.25% for the 20% and 30% separation minima increase scenarios, against and average conflict size of around 10.75% for the base and 10% separation minima increase scenarios.

3) Runway out conflict analysis

Fig. 12 shows the number of conflicts occurred at the departing runway (runway out). The graph shows that the best scenarios that obtains less conflicts is the 30% separation minima increase scenario, improving the performance of the base case scenario and achieving a reduction of almost 10 conflicts. Looking at Fig. 13(a) it can be noticed that most of the conflicts, regardless of the scenario simulated, are on the right-hand side of the chart, meaning to have high conflict size. The main reason can be found in the magnitude of the uncertainty that affects ground side, identified in the simulation model as deviation from off block time and taxi time out. The values of separation minima for the runway out range between 60 and 111 seconds, while values of deviation from off block time range between minus and plus 60 seconds. Considering also that the deviation in taxi out time ranges between minus and plus 10% of the average taxi out time, which in turn ranges between 835 and 1240 seconds, it can be clearly understood that this will affect severely the runway out conflicts performance. The consequence of these high number of conflicts and conflicts size for the runway out, are translated, in real operations, as number of aircraft queuing up and waiting times queuing at the departing runway entry. In Fig. 13(b) it can be seen that the average conflict size is quite high, between 51.7% and 52.5%, and also that the different scenarios obtain more or less similar values of conflict size. This results implies that the runway out is very sensitive to uncertainty and, on the contrary, it is not sensitive to the different scenarios, therefore the output does not change much between the different scenarios. Concerning taxiway network conflicts, we found the same average values for all the different scenarios, which is around 2. The uncertainty included in the taxi times makes it to be overloaded in all the scenarios, however, to a small extent. Regarding the terminals, all scenarios except scenario 10% of separation minima increase, show no conflicts. For these two components, uncertainty does not affect too much their performance. In real operations, having this overload would mean having aircraft waiting at the terminal apron area for a gate to be released, or aircraft waiting in some hot spot of the taxiway network. It can be also noticed that the avg OL for taxiway and terminals is almost zero, which means that the max OL has taken place for a small amount of time.

VI. CONCLUSION

This paper presents a methodology that can be used for making a decision support system for ATM operations with the objective of making them more efficient. It proposes an innovative approach that couples optimization and simulation providing two main advantages; first it provides robust solutions for the addressed problem and second it provides more feasible solutions, which in turn will provide better support to ATC than the current technology. An analysis of the different components such as airspace, runways, taxiway network and terminals is conducted, pointing out that, some of the components are more sensitive to uncertainty like the departing runway, while other components are less sensitive to it like taxiway and terminals. In general, performance in terms of conflicts are considerably improved, confirming the benefits of the methodology. These results demonstrate that the methodology is able to improve the solution of the optimization when data input is tuned up. In this context, a crucial aspect of the methodology is the choice of the right input data to be tuned. The methodology presented is general enough and it can be adapted to any particular airport configuration, however, it is aimed at close-to-congestion airports with a big amount of traffic. As a future development of the methodology, the heuristic algorithm will be included into the simulation model to fully integrate the simulation and optimization in a loop that will progressively improve the solution at each time window of the sliding window approach. Further investigation will involve the choice of the uncertainty parameter values, and the choice of the parameter to fine tune at each iteration, so to represent uncertainty as much as possible in a realistic way and improve the solution robustness.

Figure 1 :

 1 Figure 1: Schematic representation of the airport components modeled, airspace and airside.

Fig. 2 .

 2 Fig. 2. TMA representation

Fig. 3 .,

 3 Fig. 3. Airspace network E. Decision variables and objective function of the problem The decision variables considered for this problem, for each flight 𝑐𝑐 ∈ 𝐹𝐹 are: • TMA entry time 𝑐𝑐 𝑓𝑓 ∈ 𝑇𝑇 𝑓𝑓 𝑒𝑒ℎ𝐴𝐴𝐴𝐴𝐴𝐴 𝑇𝑇 𝑓𝑓 = �𝑇𝑇 𝑓𝑓 0 + 𝑗𝑗∆𝑇𝑇 � ∆𝑇𝑇 𝑚𝑚𝑚𝑚𝑚𝑚 ∆𝑇𝑇 ≤

•

  Flight type (arrival, arrival-departure and departure), • Wake turbulence category of the aircraft • Terminal number • Entry point in the TMA (route) • Entry time in the TMA • Entry speed in the TMA • Landing runway • Push back time • Departing runway Table

Fig. 4 .

 4 Fig. 4. Methodology schema of the combined approach of optimization and simulation techniques.

Fig. 5 .

 5 Fig. 5. Sliding window representation

Fig. 6 .

 6 Fig. 6. Example of the detection of an airspace conflicts in the airspace network of the simulation model.

Fig. 7 .

 7 Fig. 7. Schematic representation of both model and methodology.

1 )Fig. 8 .

 18 Fig. 8. (a) Aircraft involved in airspace conflicts. (b) Number of airspace conflicts.

Fig. 9 .

 9 Fig. 9. (a) Number of airspace conflicts that are in a specific conflict size interval. (b) Average airspace conflict size percentage.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Runway in number of conflicts.

Fig. 13 .

 13 Fig. 13. (a) Number of runway out conflicts that are in a specific conflict size interval. (b) Average runway out conflict size percentage.

4 )

 4 Taxiway network and terminal conflict analysis Simulation results for taxiway network and terminals are depicted in the two charts of Figs. 14(a) and 14(b).

Fig. 14 .

 14 Fig. 14. (a) Number of taxiway network conflicts. (b) Number of terminals conflicts.

TABLE I ICAO

 I

				Leading Aircraft
		Wake turbulence category	Heavy	Medium	Light
	Trailing Aircraft	Heavy Medium Light	4 5 6	3 3 5	3 3 3

[START_REF] Toratani | Simultaneous optimization method for trajectory and sequence for receding horizon guidance in terminal area[END_REF]

SEPARATION MINIMA 𝑆𝑆 𝑖𝑖𝑖𝑖 IN NAUTICAL MILES (NM)

TABLE II RUNWAY

 II SEPARATION MINIMA S ij

  • TMA entry speed 𝐴𝐴 𝑓𝑓 ∈ 𝑆𝑆 𝑓𝑓 𝑒𝑒ℎ𝐴𝐴𝐴𝐴𝐴𝐴 𝑆𝑆 𝑓𝑓 = �𝑆𝑆 𝑓𝑓 0 + 𝑗𝑗∆𝑆𝑆 � |𝑗𝑗| ≤ (∆𝑆𝑆 𝑚𝑚𝑜𝑜𝑚𝑚 -∆𝑆𝑆 𝑚𝑚𝑖𝑖𝑛𝑛 )/∆𝑆𝑆, 𝑗𝑗 ∈ ℤ}, ∆𝑆𝑆 𝑚𝑚𝑖𝑖𝑛𝑛 and ∆𝑆𝑆 𝑚𝑚𝑜𝑜𝑚𝑚 are the minimum and maximum values that can be assigned to 𝐴𝐴 𝑓𝑓 and ∆𝑆𝑆 is a discretized increment. In this case ∆𝑆𝑆 𝑚𝑚𝑖𝑖𝑛𝑛 is 0.9𝑆𝑆 𝑓𝑓 0 , ∆𝑆𝑆 𝑚𝑚𝑜𝑜𝑚𝑚 is 1.1𝑆𝑆 𝑓𝑓 0 and ∆𝑆𝑆 is 0.01𝑆𝑆 𝑓𝑓 0 .• Landing runway choice 𝐴𝐴 𝑓𝑓 ∈ 𝑅𝑅 𝑓𝑓 𝑒𝑒ℎ𝐴𝐴𝐴𝐴𝐴𝐴 𝑅𝑅 𝑓𝑓 is the set of available landing runways. • Push back time 𝐴𝐴𝑝𝑝 𝑓𝑓 ∈ 𝑃𝑃𝑃𝑃 𝑓𝑓 𝑒𝑒ℎ𝐴𝐴𝐴𝐴𝐴𝐴 𝑃𝑃𝑃𝑃 𝑓𝑓 = �𝑃𝑃𝑃𝑃 𝑓𝑓 0 +

	𝑗𝑗∆𝑇𝑇 �0 ≤ 𝑗𝑗 ≤	∆𝑃𝑃𝑃𝑃 𝑚𝑚𝑚𝑚𝑚𝑚 ∆𝑇𝑇

TABLE III EXAMPLE

 III

OF THE DATA INPUT

  Table Ⅳ presents the type of traffic involved.

TABLE VII VALUE

 VII OF THE OBJECTIVE FUNCTION BROKEN DOWN INTO THE DIFFERENT COMPONENTS USING THE ORIGINAL SCHEDULE.

	Objective function component	Objective function subcomponent	Individual value	Aggregated value
	Airspace conflicts	Link conflicts Node conflicts	1114 461	1575
	Runway conflicts			381
			Max OL	8
	Taxiway network conflicts	Avg OL	0.14	8.14
		T1	Max OL Avg OL	2 0.12
	Terminal conflicts	T2	Max OL Avg OL	4 0.3	8.62
		T3	Max OL Avg OL	2 0.2
	Total conflicts			1972.76

TABLE VIII

 VIII

	UNCERTAINTY FACTORS AND THEIR RESPECTIVE VALUES
	Uncertainty factor	Value
	Entry time in the TMA deviation	[-30s, +30s]
	Taxiway time deviation	[-10%, +10%]
	Pushback time deviation	[-60s, +60s]

ACKNOWLEDGMENT

"The authors would like to thank the Aviation Academy of the Amsterdam University of Applied Sciences for supporting this study and the Dutch Benelux Simulation Society (www.dutchbss.org) and EUROSIM for the dissemination of the findings of this study."


Jijju Thomas

Laurentiu Hetel

Christophe Fiter

Nathan Van De Wouw

Jean-Pierre Richard

L2- Stability

L2-stability criterion for systems with decentralized asynchronous controllers

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

L 2 -Stability Criterion for Systems with Decentralized Asynchronous Controllers Jijju Thomas 1 , Laurentiu Hetel 2 , Christophe Fiter, 3 Nathan van de Wouw 4 , and Jean-Pierre Richard 5 Abstract-This paper deals with the stability analysis of decentralized sampled-data Linear Time Invariant (LTI) control systems with asynchronous sensors and actuators. We consider the case where each controller in the decentralized setting has its own sampling and actuation frequency which translates to asynchrony between sensors and actuators. The errors induced due to sampling and asynchronicity are modelled using two different operator approaches, leading to simple L2stability criteria for the overall decentralized control system. The simplicity of the obtained criteria is illustrated by an example and simulation results exhibit the effectiveness of the approach.

I. INTRODUCTION

Decentralized control systems can be found aplenty in technological, environmental or societal environments [START_REF] Bakule | Decentralized Design of Feedback Control for Large-scale Systems[END_REF][START_REF] Bakule | Decentralized control: An overview[END_REF]. In such systems, controllers are assigned to individual subsystems, using only local plant information (see Fig. 1 for a typical example). Since the feedback scheme involved in decentralized control is local, there are a few advantages of decentralization. Firstly, a substantial amount of wiring can be avoided. Secondly, owing to the decoupled nature of the controllers, the diagnosis and maintenance is easier. The aforementioned two points also translate to overall lower running costs [START_REF] Bakule | Decentralized Design of Feedback Control for Large-scale Systems[END_REF], [START_REF] Bakule | Decentralized control: An overview[END_REF] and [START_REF] Lunze | Feedback control of large-scale systems[END_REF]. However, the design of such a decentralized control scheme may be quite complex since the local design has to be done from a global perspective. In this paper, a particular problem within the sampled-data implementation [START_REF] Chen | Optimal Sampled-Data Control Systems[END_REF][START_REF] Åström | Computer-Controlled Systems: Theory and Design[END_REF] of decentralized controls is considered. More precisely, we will analyze the effect of asynchronism between the local sampled data controllers on the overall stability of the system. At an implementation level, con-This work was supported by project UCoCoS, funded by the European Union's EU Framework Programme for Research and Innovation, Horizon H2020, Grant Agreement No: 675080. The project was also funded by ANR Project ROCC-SYS (agreement ANR-14-CE27-0008).

1 Jijju Thomas(jijju.thomas@inria.fr) is with Univ. Lille, CNRS, Inria, Centrale Lille, UMR 9189 -CRIStAL, Inria Lille Nord-Europe, F-59000 Lille, France and Eindhoven University of Technology, The Netherlands 2 Laurentiu Hetel (laurentiu.hetel@ec-lille.fr) is with Univ. Lille, CNRS, Centrale Lille, UMR 9189 -CRIStAL -Centre de Recherche en Informatique Signal et Automatique de Lille, F-59000 Lille, France.

3 Christophe Fiter (christophe.fiter@univ-lille.fr) is with Univ. Lille, CNRS, Centrale Lille, UMR 9189 -CRIStAL -Centre de Recherche en Informatique Signal et Automatique de Lille, F-59000 Lille, France. 4 Nathan van de Wouw (n.v.d.wouw@tue.nl) is with Eindhoven University of Technology, Eindhoven 5600 MB, The Netherlands, with the Department of Civil, Environmental and Geo-Engineering, University of Minnesota, Minneapolis, MN 55455 USA, and also with the Delft Center for Systems and Control, Delft University of Technology, Delft 2628 CD, The Netherlands. 5 Jean-Pierre Richard (jean-pierre.richard@ec-lille.fr) is with with Univ. Lille, CNRS, Inria, Centrale Lille, UMR 9189 CRIStAL, Inria Lille Nord-Europe, F-59000 Lille, France.

Σ 1 S 1 s 1 k H 1 a 1 k F 1 Σ 2 S 2 s 2 k H 2 a 2 k F 2 Σ n S n s n k H n a n k F n x 1 u 1 x 2 u 2 x n u n
Fig. 1: A decentralized controller setup.

trollers are usually algorithms programmed on embedded processors which work at different frequencies. Moreover, sensors and actuators are distributed over different communication channels which function aperiodically. This renders the synchronization of the different elements in control loops quite challenging [START_REF] Freris | Fundamental limits on synchronizing clocks over networks[END_REF]. This may in turn affect the overall performance of the system and even its stability as illustrated in the following example. Consider the decentralized LTI system defined by

Σ 1 : ẋ1 (t) = -2x 1 (t) -x 2 (t) + u 1 (t) Σ 2 : ẋ2 (t) = 4x 2 (t) -2.8x 1 (t) + u 2 (t) (1)
where u 1 (t) = -x 1 (t), u 2 (t) = -4.6x 2 (t) are the decentralized control inputs to systems Σ 1 and Σ 2 respectively, and x1 (t), x2 (t) are the state values obtained through sampling and hold. In the event both systems Σ 1 and Σ 2 are sampled periodically as well as synchronously with a sampling period T = 0.59 (i.e., xi (t) = x i (kT ), ∀t ∈ [kT, (k + 1)T ), i = {1, 2}), the overall system is stable as illustrated in Figure 2a. However, as can be seen from Figure 2b, the stability is affected when the sampling is periodic but control loops are asynchronous. Figure 2b presents the case when a shift δ = 0.2 is introduced in the sampling of the second state, i.e., when x2 (t) = x 2 (kT + δ), ∀t ∈ [kT + δ, (k + 1)T + δ). The stability problem can become even more complex when both the sensors and actuators involved within individual control loops are asynchronous.

In this paper, we will address the problem of stability analysis for the case of LTI systems with decentralized sampled-data linear controllers subject to asynchronicity. More precisely, we consider that each sampled-data controller has its own sampling and actuation frequencies. This particular problem setting gives rise to complexities induced by sampling, asynchronicity, network effects, etc. The decentralized control problem that we introduce in this paper is unique to the best of our knowledge.

Although the problem considered in this paper is novel, stability analysis methods have been proposed for centralized controllers subjected to sampling and asynchronism between sensors and actuators [START_REF] Wakaiki | Stabilization of systems with asynchronous sensors and controllers[END_REF][START_REF] Wakaiki | Stabilization of networked control systems with clock offsets[END_REF][START_REF] Wakaiki | Control under clock offsets and actuator saturation[END_REF]. However, the period between sampling and actuation instants was treated to be constant and all the system state data was considered to be sampled at same time instants. In the scope of sampled-data and networked control systems [START_REF] Heemels | Networked control systems with communication constraints: Tradeoffs between transmission intervals, delays and performance[END_REF][START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF][START_REF] Fiter | Stability of LTI systems with distributed sensors and aperiodic sampling[END_REF], there are a few similarities with the problem we consider. For example, decentralized event-triggered control and delay introduced due to network effects can also be seen as a form of asynchronism [START_REF] Dolk | Output-based and decentralized dynamic event-triggered control with guaranteed Lpgain performance and zeno-freeness[END_REF]. In comparison to the very few existing results addressing problems similar to the one considered in this paper, we propose a novel and simple approach that guarantees stability.

The main contribution of this paper is to provide approaches for L 2 -stability analysis of decentralized sampleddata controllers. For the sake of generality, we consider the sampling and actuation intervals to be time-varying and possibly unknown (but bounded). We take into account the asynchronicity between individual controllers as well as the asynchronicity between sensors and actuators within a local control loop. By using tools based on input-output methods [START_REF] Megretski | System analysis via integral quadratic constraints[END_REF][START_REF] Zames | On the input-output stability of time-varying nonlinear feedback systems part one: Conditions derived using concepts of loop gain, conicity, and positivity[END_REF], related to the ones previously used for systems with delays [START_REF] Jun | Iqc robustness analysis for time-delay systems[END_REF][START_REF] Kao | Simple stability criteria for systems with time-varying delays[END_REF][START_REF] Kao | Stability analysis of systems with uncertain time-varying delays[END_REF][START_REF] Fridman | Input-output approach to stability and l2gain analysis of systems with time-varying delays[END_REF], we provide two novel and different stability analysis methods based on easy-to-check frequencydomain criteria.

The remainder of this paper has been structured as follows. In Section II, we introduce the problem formulation, followed by technical preliminaries. Sections III and IV deal with the transformation of the closed-loop sampled-data dynamics into a feedback interconnection model and provide a stability criterion. Section V provides a numerical example corroborating the presented results.

Notations: R is the set of all real numbers, implying R n is the set of all n-dimensional real vectors. Diag(M 1 , M 2 , ..., M n ) is the block diagonal matrix with elements M i of appropriate dimensions. L 2e (a, b) is the extended L 2 -space of all square integrable and Lebesgue measurable functions defined on the interval [a, b], with the L 2 -norm defined as q 2 L2 = q, q , and the inner product p, q = b a p(s) T q(s)ds.

II. PROBLEM FORMULATION AND TECHNICAL PRELIMINARIES

A. Motivating Problem 1) System Model: The system under consideration consists of a set of Linear Time Invariant (LTI) systems, wherein each individual system is influenced by its corresponding control input and other system states. Figure 1 depicts this decentralized control configuration. Consider that the dynamics of each LTI system (denoted Σ i ) is given by

ẋi (t) = A i x i (t) + B i u i (t) + M j=1,i =j A ij x j (t), ∀t ∈ R, (2) with i ∈ {1, 2, .., M }, x i (t) ∈ R ni and u i (t) ∈ R mi .
The matrices A i , B i and A ij are of appropriate dimensions. The term A ij x j (t) denotes the influence of the states of the j th plant Σ j on the dynamics of system Σ i . Here, we consider the case where the control of the global system is linear. Furthermore, we assume that it is decentralized in the sense that the control input u i (t) only depends on the local state variables x i (t). Furthermore, we consider that the control inputs are asynchronous. The system states x i (t) are sampled according to a sampling sequence {s i k } k∈Z defined by

{s i k : s i k+1 -s i k = h i k , k ∈ Z, i ∈ {1, 2, .., M }}. (3)
The sequence of sampling intervals

{h i k } k∈Z satisfying h i k ∈ [h i , hi ]
takes into account imperfection in sampling caused by e.g. jitter, data packet dropouts, etc. Note that the sampling instants of different systems are not necessarily synchronous. The control input u i (t) based on x i (s i k ) will be implemented at a time instant a i k at the level of the actuator of system Σ i . We consider that the sequence of actuation times {a i k } k∈Z satisfies

{a i k : a i k = s i k + η i k , η i k ≤ h i k , k ∈ Z, i ∈ {1, 2, . . . , M }}, (4)
where η k i ∈ [η i , ηi ] denotes the asynchrony between sensors and actuators. Such an asynchrony may be due to network delays, control computation delay, etc. Based on this consideration, the control input to the system Σ i is given by the sampled-data decentralized static state-feedback law

u i (t) = F i x i (s i k ), ∀t ∈ [a i k , a i k+1 ). (5)
The goal of this paper is to analyse the stability of the system defined by ( 2), (3), ( 4) and (5)1 .

B. Preliminaries

We introduce some basic concepts of linear operator theory that are used in this paper. An operator G : L 2e (a, b) → L 2e (c, d) receives an input p ∈ L 2e (a, b) and produces an output q ∈ L 2e (c, d). 2) Feedback interconnection: The standard feedback interconnection of two operators G 1 and G 2 , is given by

G 1 + f G 2 + g p 1 q 1 p 2 q 2
Σ G1G2 : p 2 = G 1 p 1 + f p 1 = G 2 p 2 + g. (6)
Figure 3 shows the graphical representation of the standard feedback interconnection.

3) Well-posed system: A feedback system is said to be well-posed if all the closed-loop transfer matrices are welldefined and proper [START_REF] Zhou | Essentials of Robust Control[END_REF]. The implication for the standard feedback interconnection Σ G1G2 given by ( 6) is that the well-posedness is guaranteed only if

I -G 1 G 2 is invertible.
4) L 2 -stability: An operator G is said to be L 2 -stable if it has a finite L 2 -gain [START_REF] Schaft | L2-Gain and Passivity Techniques in Nonlinear Control[END_REF].

5) Small-Gain Theorem: A feedback interconnection of two operators G 1 and G 2 given by ( 6), has a finite L 2 -gain for the mapping

f g → p 1 p 2 (7) if γ(G 1 )γ(G 2 ) < 1, (8)
where γ(G 1 ) and γ(G 2 ) are the L 2 -gain of the operators G 1 and G 2 respectively [START_REF] Khalil | Nonlinear systems[END_REF].

In this paper, we will use an operator approach to take into account the asynchrony in decentralized control loops.

Stability will be analysed in the L 2 -sense by modelling the system and the effects of sampling and asynchrony using operators. Two methods will be presented. The first method models the overall effect of sampling and asynchrony between sensors and actuators in a global manner using one operator. The second method takes the difference between the effects of sampling and asynchrony using two separate operators.

III. STABILITY ANALYSIS USING A SINGLE OPERATOR FOR SAMPLING AND ASYNCHRONY

The configuration shown in Figure 1, defined by (2) can also be expressed by the standard state-space equation

Ẋ(t) = AX(t) + BU (t), ∀ t ∈ R, (9)
wherein the system state X(t) ∈ R n and the control input U (t) ∈ R m can be decomposed as

X(t) = x T 1 (t) x T 2 (t) ... x T M (t) T , U (t) = u T 1 (t) u T 2 (t) ... u T M (t) T . (10)
with

x i (t) ∈ R ni , u i (t) ∈ R mi , M i=1 n i = n and M i=1 m i = m.
Similarly, the system matrix A ∈ R n×n and the input matrix B ∈ R n×m are given by

A =      A 1 A 12 . . . A 1M A 21 A 2 . . . A 2M . . . . . . . . . . . . A M 1 A M 2 . . . A M      , B = diag(B 1 , B 2 , . . . , B M ). (11)

A. System Model Reformulation

Let xi (t) represent the information used in computing the control input u i (t) under the influence of asynchrony and sampling:

xi (t) = x i (s i k ), ∀t ∈ [a i k , a i k+1 ), X(t) = xT 1 (t) xT 2 (t) ... xT M (t) T (12)
with k ∈ Z and i ∈ {1, 2, . . . , M }. Consequently, the decentralized control law can be defined as

U (t) = F X(t) (13)
with F = diag(F 1 , F 2 , ..., F M ). The closed-loop system model defined by ( 2), ( 3), ( 4), ( 5) can, therefore, be formulated as follows:

Ẋ(t) = AX(t) + BF X(t) = (A + BF )X(t) + BF ( X(t) -X(t)) = A cl X(t) + B cl E(t), (14)
where A cl := A + BF and B cl := BF . The vector E(t) in ( 14) represents the error induced in the system (9) by sampling and asynchrony, i.e,

E(t) := X(t) -X(t) = e T 1 (t) e T 2 (t) ... e T M (t) T , (15) where
e i (t) = xi (t) -x i (t), = x i (s i k ) -x i (t). (16)
Choosing an auxiliary output

Y (t) = Ẋ(t) = ẋ1 (t) T ẋ2 (t) T ... ẋM (t) T T , = y T 1 (t) y T 2 (t) ... y T M (t) T ( 17
)
and by using an integral operator ∆ i :

L ni 2e (-∞, ∞) → L ni 2e (-∞, ∞)
, we can rewrite [START_REF] Jun | Iqc robustness analysis for time-delay systems[END_REF] as follows:

e i (t) = (∆ i y i )(t) := - t s i k y i (θ)dθ, ∀t ∈ [a i k , a i k+1 ), k ∈ Z. ( 18
)
The operator ∆ i accounts for the error induced in the system Σ i (in closed-loop with its local controller) by sampling and asynchrony.

B. Stability Analysis

Motivated by the problem under consideration in Section II-A, we study the stability of the feedback interconnection of G and ∆ in the standard form (see Figure 4), defined by

Σ G∆ : Y = GE + f E = ∆Y + g, (19)
where f, g ∈ L n 2e (-∞, ∞). The operator G is defined by the transfer function

G(s) = C cl (sI -A cl ) -1 B cl + D cl , (20)
where

C cl = A cl = A + BF, D cl = B cl = BF, (21)
and the operator ∆ = diag(∆ 1 , ∆ 2 , ..., ∆ M ), where ∆ i is defined by [START_REF] Kao | Stability analysis of systems with uncertain time-varying delays[END_REF]. The feedback interconnection Σ G∆ is equivalent to the decentralized control system given by ( 2), ( 3), ( 4) and ( 5), affected by perturbations in the measured state value. That is,

ẋi (t) = A i x i (t) + B i F i x i (s i k ) + w i (t) + M j=1,i =j A ij x j (t), w i (t) = g i (t) + t s i k f i (s)ds, x i (0) = 0.
(22) f i and g i are components of f and g, with appropriate dimensions. Before providing the stability criterion, we introduce the following technical lemma which is an adaptation of the result in [START_REF] Kao | Simple stability criteria for systems with time-varying delays[END_REF] for continuous-time systems with time delay.

Lemma 1: The L 2 induced norm of the operator ∆ is upper-bounded by γ, where

γ = M max i=1 { hi + ηi }.
(23) Proof: The proof is available in the technical report [START_REF] Thomas | L 2 -stability criterion for systems with decentralized asynchronous controllers[END_REF]. Using Lemma 1, the following stability result can be obtained.

Theorem 2: The feedback interconnection of operators G and ∆, denoted by

Σ G∆ in (19) is L 2 -stable if sup ω∈R σ G(jω) < M max i=1 { hi + ηi } -1 , ( 24
)
where G(jω) = C cl (jωI -A cl ) -1 B cl + D cl is the frequency response function matrix of the system defined by ( 14) and ( 17), with the matrices A cl , B cl , C cl and D cl defined in [START_REF] Schaft | L2-Gain and Passivity Techniques in Nonlinear Control[END_REF], and σ G(jω) is the largest singular value of the G(jω)

Proof: The proof is available in the technical report [START_REF] Thomas | L 2 -stability criterion for systems with decentralized asynchronous controllers[END_REF].

IV. STABILITY ANALYSIS USING TWO SEPARATE

OPERATORS We have seen that Theorem 2 provides an easy-to-check criterion for stability analysis of the closed-loop LTI system [START_REF] Wakaiki | Control under clock offsets and actuator saturation[END_REF], since it only requires a frequency-domain check of an LTI system. However, this result may be conservative since both the effects of sampling and asynchrony are modelled using a global operator. Below, we propose an alternative approach in which the error induced by the sampling and asynchrony are considered separately, in terms of an operator that represents the effect of sampling and hold, and an operator that represents the delay induced by asynchrony. This alleviates conservatism by providing flexibility in employing more accurate function bounding inequalities.

A. Feedback interconnection system representation

In this section, we show that the operator ∆ can be decomposed into two separate operators. The operator ∆ sam represents the error induced by sampling whereas the operator ∆ asy denotes the error induced by asynchrony between the sensors and actuators. Let us recall the definition of E(t) defined as in [START_REF] Zames | On the input-output stability of time-varying nonlinear feedback systems part one: Conditions derived using concepts of loop gain, conicity, and positivity[END_REF]. Let xi (t) denote the sampled version of x i (t), along the sampling sequence {s i k } k∈Z and be given by

xi (t) = x i (s i k ) ∀t ∈ [s i k , s i k+1 ). (25)
We care to stress the difference between xi (t) in (25) and xi (t) in ( 12) in terms of their domain of definition. Let us define

e sam i (t) := xi (t) -x i (t), ∀t ∈ [s i k , s i k+1 ). (26)
Note that e sam i (t) corresponds to the error between the signal x(t) and it's sampled version (see Figure 5 for a graphical illustration). Given the signal y i (t) in [START_REF] Kao | Simple stability criteria for systems with time-varying delays[END_REF], the sampling induced error can be characterized by

e sam i (t) = - t s i k y i (θ)dθ =: (∆ sam i y i )(t). (27)
Considering the rectangular signal xi (t), representing the sampled version of system state x i (t), and the signal xi (t) as given in [START_REF] Fiter | Stability of LTI systems with distributed sensors and aperiodic sampling[END_REF] representing the signal actually used at the level of actuators, the effect of asynchrony can be captured by introducing an error e asy i (t) = xi (t) -xi (t), ∀t ≥ 0 (28) as illustrated in Figure 6. Let us remark that e asy i (t) :=

x i (s i k-1 ) -x i (s i k ), ∀t ∈ [s i k , a i k ) 0, ∀t ∈ [a i k , s i k+1 ) (29)
Considering y i (t) as given in [START_REF] Kao | Simple stability criteria for systems with time-varying delays[END_REF], we define e asy i (t) = (∆ asy i y)(t)

:=    - s i k s i k-1 y i (θ)dθ, ∀t ∈ [s i k , a i k ) 0, ∀t ∈ [a i k , s i k+1 ). (30)
we have xi (t) = x i (t) + e sam i (t) + e asy i (t), which leads to the decomposition of ∆ i in [START_REF] Kao | Stability analysis of systems with uncertain time-varying delays[END_REF] given by

(∆ i y i )(t) = (∆ sam i y i )(t) + (∆ asy i y i )(t), (32)
as shown in Figure 7. Then, we have for the reformulated system ( 14) and ( 17),

E(t) = e T 1 (t) e T 2 (t) ... e T M (t) T (33)
with

e i (t) = (∆ i y)(t) = (∆ sam i y)(t) + (∆ asy i y)(t), ∀t ∈ [s i k , s i k+1 ). (34)

B. Stability Analysis

In the following lemma, we compute the L 2 -norm of the operator ∆ i by upper-bounding each of the operators introduced by the decomposition shown in (34), thereby providing a bound on the operator ∆. Lemma 3: The L 2 -induced norm of the operator ∆ is upper-bounded by γ 1 , where

γ 1 = M max i=1 2 hi π + hi ηi . (35)
Proof: The proof is available in the technical report [START_REF] Thomas | L 2 -stability criterion for systems with decentralized asynchronous controllers[END_REF]. Based on Lemma 3, we provide in Theorem 4 a less conservative and also easy-to-check stability criterion for the L 2 -stability of the feedback interconnection Σ G∆ .

Theorem 4: The feedback interconnection Σ G∆ of operators G and ∆ as defined in [START_REF] Fridman | Input-output approach to stability and l2gain analysis of systems with time-varying delays[END_REF], where ∆ satisfies the decomposition (34), is L 2 -stable if

sup ω∈R σ G(jω) < M max i=1 2 hi π + hi ηi -1 , (36)
where G(jω) = C cl (jωI -A cl ) -1 B cl + D cl is the frequency response function matrix of the closed-loop system defined by ( 14) and ( 17), with the matrices A cl , B cl , C cl and D cl defined in [START_REF] Schaft | L2-Gain and Passivity Techniques in Nonlinear Control[END_REF]. Proof: The proof is available in the technical report [START_REF] Thomas | L 2 -stability criterion for systems with decentralized asynchronous controllers[END_REF].

V. SIMULATION RESULTS

In this section, we apply the stability criteria provided in Theorems 2 and 4 to the decentralized LTI system previously considered in Section I, equation [START_REF] Bakule | Decentralized Design of Feedback Control for Large-scale Systems[END_REF]. Expressing the decentralized system (1) in the standard state-space model given by (9), we have

A = -2 -1 2.8 4 , B = 1 0 0 1 , F = -1 0 0 -4.6 , (37)
which provides

A cl = C cl = A + BF = -3 -1 2.8 -0.6 , B cl = D cl = BF = -1 0 0 -4.6 . (38)
The L 2 -norm of the operator G can easily be obtained from the transfer function in [START_REF] Zhou | Essentials of Robust Control[END_REF] Fig. 8: hi and ηi satisfying the stability criteria obtained using Theorem 2 (39) and Theorem 4 (40). Fig. 9: The evolution of (a) system states and (b) control input for (37), with parameters h1 , h2 , η1 and η2 chosen according to the feasibility region shown in Figure 8.

The feasible values of hi and ηi satisfying (39) and (40) are shown in Figure 8. It is quite clear that the criterion obtained using Theorem 4, given by (40), provides less conservative results in comparison to the criterion obtained using Theorem 2, given by (39). This corroborates our theoretical result that by encompassing the error induced due to sampling and asynchrony within two separate operators, we obtain a less conservative result. Additionally, the criterion is simple to employ since the only computation involved is in obtaining the L 2 -norm of the system operator, which can be computed easily using the H ∞ norm of it's transfer function [START_REF] Zhou | Robust and Optimal Control[END_REF]. Choosing the parameters h1 = 0.117, h2 = 0.1035, η1 = 0.0945 and η2 = 0.0405, we simulate the system by introducing a rectangular wave perturbation w i (t) = 0.5, ∀t ≤ a i 10 , i ∈ {1, 2}. We can see in Figure 9a that the system (37) is indeed stable.

VI. CONCLUSION

In this paper, the stability analysis problem for LTI systems with decentralized sampled-data linear controllers subjected to asynchrony has been studied. Two different stability analysis methods based on easy-to-check frequency-domain criteria have been provided. The method primarily included modelling the error induced by sampling and asynchrony using operators, and obtaining the L 2 -gain bounds on these operators. The effectiveness of the method was illustrated using numerical simulations.

Fig. 3 : 2 . 1 )

 321 Fig. 3: Standard feedback interconnection of operators G 1 and G 2 . 1) Bounded operators: The operator G : L 2e (a, b) → L 2e (c, d) is said to be bounded if there exists a constant γ ∈ R so that G(p) L2 ≤ γ p L2 for all p ∈ L 2e (a, b). The minimal constant γ satisfying the aforementioned inequality is called the induced L 2 -gain of the operator G, and is denoted by G L2 or γ(G).2) Feedback interconnection: The standard feedback interconnection of two operators G 1 and G 2 , is given by

Fig. 4 :

 4 Fig. 4: The system (14) as a feedback interconnection

Fig. 7 :

 7 Fig. 7: The decomposition of operator ∆ i into operators ∆ sam i

  for the given x i (t) and xi (t).

				1
			xi(t), xi(t)	0	xi(t) xi(t)		e sam i	(t)
				-1
				0	2	4	6	8	10
						t
				2		e sam i	(t) = xi(t) -xi(t)
			(t)			Zero
			i e sam	0
				-2
				0	2	4	6	8	10
						t
	Fig. 5: The top figure gives a sample signal x i (t) (in black)
	and its sampled version xi (t) (in red). The bottom figure
	shows the resetting nature (to zero) of the sampling-induced
	error e sam i xi(t), xi(t) (t) t 0 1 2 -1 1 0	3 xi(t) xi(t)	4	e asy i (t)	5
				2
			e asy i (t)	0 1		e asy i (t) = xi(t) -xi(t)
				-1
				0	1	2	3	4	5
						t
	Fig. 6: The top figure gives the sampled signal xi (t) (in red)
	and the version xi (t) (in blue) that is used at the actuator
	level. The bottom figure shows the piece-wise constant error
	e	asy(t) i	introduced due to asynchrony for the given x i (t),
	xi (t) and xi (t).
	Since			xi (t) = x i (t) + e sam i	(t),
					xi (t) = xi (t) + e asy i (t),

  and is 5.2143. By employing Theorem 2, we can state that the system remains stable for all hi and ηi satisfying

		0.12
							hi and ηi satisfying (39)
		0.1
							hi and ηi satisfying (40)
		0.08
	ηi	0.06
		0.04
		0.02
		0
		0	0.05	0.1	hi	0.15	0.2	0.25	0.3
									2 max i=1 { hi + ηi } <	1 5.2143	.	(39)
									Similarly, using Theorem 4 L 2 -stability is guaranteed if
									2 max i=1	2 hi π	+ hi ηi <	1 5.2143	.	(40)

The exact mathematical concept of stability that we use in this paper will be formalized in Section II-B


Jijju Thomas
email: jijju.thomas@phd.centralelille.fr

Christophe Fiter
email: christophe.fiter@univ-lille1.fr

Laurentiu Hetel
email: laurentiu.hetel@ec-lille.fr

Nathan Van De Wouw

Jean-Pierre Richard
email: jean-pierre.richard@ec-lille.fr

Dissipativity based stability criterion for aperiodic sampled-data systems subject to time-delay

Keywords: Sampled-data system, dissipativity based stability, system with time-delay

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The last few decades has seen a tremendous increase in research on sampled-data systems as well as their applications [START_REF] Åström | Computer-Controlled Systems: Theory and Design[END_REF]; [START_REF] Chen | Optimal Sampled-Data Control Systems[END_REF]; [START_REF] Hespanha | A survey of recent results in networked control systems[END_REF]]. The stability analysis of such systems is an important research track that has been explored over the years and an overview of the different approaches that have been proposed can be found in [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF]]. The Time-delay approach, embeds the sampling phenomenon into a time-delay problem and generally consists of Lyapunov-Krasovskii Functional (LKF) based stability criteria [START_REF] Mikheev | Asymptotic analysis of digital control systems[END_REF]; [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF]; van de Wouw et al. (2010); [START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF]]. Since sampleddata systems exhibit both continuous and discrete dynamics, the second approach, namely the Hybrid systems approach was developed [START_REF] Naghshtabrizi | Exponential stability of impulsive systems with application to uncertain sampled-data systems[END_REF]; [START_REF] Nesić | Input-to-state stability of networked control systems[END_REF]]. In the Discrete-time approach, stability criteria are obtained using system integration over sampling intervals and by convex embedding of the state transition matrix between sampling instants [START_REF] Fujioka | Stability analysis of systems with aperiodic sample-and-hold devices[END_REF]; [START_REF] Cloosterman | Controller synthesis for networked control systems[END_REF][START_REF] Van De Wouw | Tracking control for sampleddata systems with uncertain timevarying sampling intervals and delays[END_REF]]. The third, Input-output stability approach treats the error induced by sampling as a perturbation to the continuoustime control system, which aids in employing classical This work was supported by project UCoCoS, funded by the European Union's EU Framework Programme for Research and Innovation, Horizon H2020, Grant Agreement No: 675080. The project was also funded by ANR Project ROCC-SYS (agreement ANR-14-CE27-0008).

robust control tools to analyse the stability of the system [START_REF] Mirkin | Some remarks on the use of timevarying delay to model sample-and-hold circuits[END_REF]; [START_REF] Fujioka | Stability analysis of systems with aperiodic sample-and-hold devices[END_REF]; [START_REF] Omran | On the stability of input-affine nonlinear systems with sampled-data control[END_REF]]. The approach draws similarities with the input-output stability analysis of time-delay systems as provided in [START_REF] Fridman | Input-output approach to stability and l2-gain analysis of systems with time-varying delays[END_REF]; [START_REF] Kao | Stability analysis of systems with uncertain time-varying delays[END_REF]; [START_REF] Kao | Simple stability criteria for systems with time-varying delays[END_REF]] for the case of continuous time-varying delay. In this article, we focus on the input-output approach for stability analysis of aperiodic sampled-data systems subject to time-delay. A primary advantage of this approach is that it is intuitively simple to develop. The problem trickles down to the classical robust control framework that could possibly take into account various perturbations and non-linearities. However, several problems remain to be solved in the input-output framework. For example, the existing results provide only asymptotic stability criteria for sample-data systems. However, in numerous practical scenarios, it is desirable to at least have a measure of the system decay-rate so as to provide a basic performance. Furthermore, while the stability analysis of LTI systems in the presence of delay and sampling has been studied separately using the input-output approach, a result for sampled-data systems with time delay is yet to be provided. The extension towards this direction is challenging since directly embedding sampling as an additional delay leads to a discontinuous delay and the existing approaches in the input-output framework do not cover this scenario. In this paper, we close the aforementioned gap for Linear Time Invariant (LTI) systems. The main contribution of this article is to provide an exponential stability criterion for systems with aperiodic sampling and time-delay. We provide a detailed analysis of sampling and delay, by taking into account the specific discontinuities in delay. The provided result will lay the foundation for a possible extension to the stability analysis of Non-linear sampleddata systems with delay.

PROBLEM STATEMENT

We consider the following LTI system:

ẋ(t) = Ax(t) + Bu(t), ∀t ≥ 0, x(0) = x 0 , (1)
where x(t) ∈ R n , u(t) ∈ R m , and A and B are matrices of appropriate dimensions. The time instant t s 0 specifies the first sampling instant, and belongs to the sampling sequence {t s k } k∈N defined by

t s k+1 -t s k = h k , ∀k ∈ N.
(2) The possibly time-varying sequence of sampling intervals {h k } k∈N satisfying h k ∈ [h, h], with h > 0, considers imperfections in sampling caused by jitter, data packet dropouts, etc. The actuation timing sequence {t a k } k∈N , implying the sequence of time-instants at which the control input u(t) based on x(t s k ) will be implemented at the level of the actuator, is given by

t a k = t s k + τ k , t a k ≤ t a k+1 , ∀k ∈ N (3) where τ k ∈ [τ , τ ]
represents the time-varying delay between sampling and actuations instants. Without loss of generality, we consider that the first actuation instance occurs at time t a 0 = τ + h, while the first sampling instant occurs at t s 0 = t a 0 -τ 0 . Based on the sequences ( 2) and (3), we have the control input u(t) in (1) satisfying the state-feedback and zero-order hold strategies given by

u(t) = 0, ∀t ∈ [0, t a 0 ), Kx(t s k ), ∀t ∈ [t a k , t a k+1 ), k ∈ N, (4)
where K ∈ R n×n is the controller gain. The objective of this paper is to analyse the stability of the system defined by (1), ( 2), ( 3) and (4) using input-output stability analysis approach.

MAIN RESULTS

System Reformulation

In this section, we reformulate the sampled-data system model to include the effects of sampling and delay as perturbations to the system defined by ( 1)-( 4). We have,

for all t ∈ [t a k , t a k+1 ), k ∈ N, ẋ(t) = Ax(t) + Bu(t), = Ax(t) + BKx(t) -BKx(t) + BKx(t s k ), = A cl x(t) + B cl e(t), (5)
where A cl = A + BK, B cl = BK, and e(t) = x(t s k ) -x(t), ∀t ∈ [t a k , t a k+1 ), k ∈ N (6) represents the error induced by sampling and delay. By choosing an auxiliary output y(t) = ẋ(t), we can reformulate system (1)-(4) as follows:

ẋ(t) = A cl x(t) + B cl e(t), ∀t ≥ t a 0 y(t) = ẋ(t) = A cl x(t) + B cl e(t), e(t) = x(t s k ) -x(t) = - t t s k y(s)ds, ∀t ∈ [t a k , t a k+1 ), k ∈ N. ( 7
)
For t ∈ [0, t a 0 ), the system ( 1) is in open-loop and therefore can be remodelled as ẋ(t) = Ax(t), y(t) = ẋ(t) = Ax(t), e(t) = 0.

(8)

Stability criterion

By exploiting the properties of the error e(t) and using arguments inspired by the dissipativity theory, we derive the following result. Theorem 1. The system defined by ( 1)-( 4) or ( 7) and ( 8) is exponentially stable with a decay rate α/2 if there exist symmetric positive definite matrices P , R 1 , and R 2 such that the following linear matrix inequalities are feasible,

Q + C T cl (γ 2 1 R 1 + γ 2 2 R 2 )C cl P B cl + C T cl (γ 2 1 R 1 + γ 2 2 R 2 )D cl D T cl (γ 2 1 R 1 + γ 2 2 R 2 )D cl -(R 1 + R 2 ) < 0, (9)
and

Q + γ 2 2 C T cl R 2 C cl P B cl + γ 2 2 C T cl R 2 D cl γ 2 2 D T cl R 2 D cl -R 2 < 0, (10)
where

Q = A T cl P + P A cl + αP, (11) and γ 2 1 = ( h + τ ) 2 , γ 2 2 = ( h + τ ) 2 e α( h+τ ) . (12)
Proof: We provide as follows, a sketch of the proof. The result is inspired by the frequency domain criteria provided by [START_REF] Fujioka | Stability analysis of systems with aperiodic sample-and-hold devices[END_REF]] for sampled-data systems without delay and further extended in [START_REF] Omran | Stability analysis of some classes of input-affine nonlinear systems with aperiodic sampled-data control[END_REF]]. The proof is based on the existence of a storage function V : R n → R + , V (x) = x T P x, and a supply function S t, y(t), e(t) : R + × R n × R n → R, such that the exponential dissipativity inequality V (x(t))+αV (x(t)) ≤ e -ατ (t) S t, y(t), e(t) , ∀t ≥ t a 0 , (13) with τ (t) = t -t a 0 , is satisfied. For the supply function given by S θ, y(θ), e(θ) = e T (θ) R 1 + e ατ (θ) R 2 e(θ)

-y T (θ) γ 2 1 R 1 + e ατ (θ) γ 2 2 R 2 y(θ), (14)
with

R 1 = R T 1 > 0, R 2 = R T 2 > 0,
and γ 1 , γ 2 given by ( 12), it can be shown that the inequalities

t 0 S θ, y(θ), e(θ) dθ ≤ 0, ∀t ≥ 0, ( 15
) - t a 0 0 S θ, y(θ), e(θ) dθ ≤ ηV (x(t a 0 )), (16)
hold with

η = h + τ eigmin(P ) max θ∈[0,t a 0 ]
eigmax (e Aτ (θ) θ) ) .

) T A T (γ 2 1 R 1 + e ατ (θ) γ 2 2 R 2 )A(e Aτ (
(17) Using simple manipulations, ( 13), ( 15) and ( 16) imply V (x(t)) ≤ Ce -αt V (x(t a 0 )), ∀t ≥ t a 0 , (18) for some C > 0. To conclude, the set of LMI conditions ( 9) and ( 10) are sufficient for the dissipation inequality (13) to

Numerical example

Consider the system (1)-( 4) characterized by the parameters [START_REF] Zhang | Stability Analysis of Networked Control Systems[END_REF]]

A = 1 3 2 1 , B = 1 0.6 , K = -[1 6] . (19)
By virtue of Theorem 1, we can compute the maximum allowable values of h + τ with respect to α. The trade-off region thus obtained is shown in Figure 1. For example, considering α = 0.001, we obtain the stability domain shown in Figure 2. Choosing the following parameters in accordance with the stability region, h = 0.1,

τ = 0.06, (20)
along with an initial condition x 0 = [0.8, 0.4] T , we can see from Figure 3 that the system is indeed exponentially stable. x 1 (t)

x 2 (t)

CONCLUSION

In this paper we have considered the stability analysis of aperiodic sampled-data systems subjected to time-varying delay. A novel stability criterion based on the dissipativity theory has been provided. The criterion is obtained by designing a suitable supply function that satisfies a general dissipativity condition for stability. The effectiveness of the obtained result has been corroborated through simulation results.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Trade-off between α and h + τ for system (19), satisfying Theorem 1

Fig. 3 .

 3 Fig. 3. Evolution of system (19) with parameters (20).