Review of risk assessment status for air traffic

Review of risk assessment status for air traffic Henk Blom, Jaroslav Krystul, Pascal Lezaud, G.J. (bert) Bakker, Margriet B. Klompstra, Bart Klein Obbing, Sybert H. Stroeve, Hans H. de Jong To cite this version: Henk Blom, Jaroslav Krystul, Pascal Lezaud, G.J. (bert) Bakker, Margriet B. Klompstra, et al.. Review of risk assessment status for air traffic. 2007. ฀hal-00940849฀ HAL Id: hal-00940849 https://hal-enac.archives-ouvertes.fr/hal-00940849 Submitted on 5 May 2014 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à 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. Project no. TREN/07/FP6AE/S07.71574/037180 IFLY iFly Safety, Complexity and Responsibility based design and validation of highly automated Air Traffic Management Specific Targeted Research Projects (STREP) Thematic Priority 1.3.1.4.g Aeronautics and Space iFly Deliverable D7.2.a Review of risk assessment status for air traffic Due date of deliverable: 22 November 2007 Actual submission date: 8 November 2007 Start date of project: 22 May 2007 Duration: 39 months Version: Draft 0.3 Project co-funded by the European Commission within the Sixth Framework Programme (2002-2006) PU PP RE CO Dissemination Level Public Restricted to other programme participants (including the Commission Services) Restricted to a group specified by the consortium (including the Commission Services) Confidential, only for members of the consortium (including the Commission Services) DOCUMENT CONTROL SHEET Title of document: Review of risk assessment status for air traffic Authors of document: Deliverable number: Project acronym: Project title: D7.2.a iFly Project no.: Safety, Complexity and Responsibility based design and validation of highly automated Air Traffic Management TREN/07/FP6AE/S07.71574/037180 IFLY Instrument: Specific Targeted Research Projects (STREP) Thematic Priority: 1.3.1.4.g Aeronautics and Space DOCUMENT CHANGE LOG Version # Issue Date Sections affected Relevant information 0.1 14 June 2007 All 1st draft 0.2 18 October 2007 All 2nd draft 0.3 8 Nov. 2007 All 3rd draft … 1.0 Version 0.3 Organisation Authors Henk A.P. Blom NLR Jaroslav Krystul TWEN Pascal Lezaud DTI/SDER G.J. (Bert) Bakker NLR Margriet B. Klompstra NLR Bart Klein Obbing NLR Sybert H. Stroeve NLR Hans H. de Jong NLR Internal reviewers External reviewers Signature/Date Contents 1 Introduction 1.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Accident risk assessment in ATM . . . . . . . . . . . . . . . . . . . . . . . 1.3 Organisation of this report . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 3 6 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Safety Risk Assessment Steps . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Monte Carlo Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Active Runway Crossing Example . . . . . . . . . . . . . . . . . . 2.3.2 Safety Relevant Scenarios . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Multi-Agent Situation Awareness in the Simulation Model . . . . 2.3.4 Dynamic Stochastic Modelling . . . . . . . . . . . . . . . . . . . . 2.4 Use of Simulation Model in Risk Assessment . . . . . . . . . . . . . . . . 2.4.1 Does the Simulation Model Cover the Identified Hazards? . . . . 2.4.2 Parametrisation of the Simulation Model . . . . . . . . . . . . . . 2.4.3 Speeding up Monte Carlo Simulations . . . . . . . . . . . . . . . . 2.4.4 Validation of the assessed risk level . . . . . . . . . . . . . . . . . 2.5 Monte Carlo Simulation Results . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Assessed risk levels . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Who contributes to safety risk reduction? . . . . . . . . . . . . . 2.5.3 Comparison against expert based results . . . . . . . . . . . . . . 2.6 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 7 7 13 13 14 15 16 18 18 19 19 20 21 21 22 24 24 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Safety verification of free flight air traffic . . . . . . . 3.1.2 Probabilistic reachability analysis . . . . . . . . . . . . 3.1.3 Sequential Monte Carlo simulation . . . . . . . . . . . 3.1.4 Development of Monte Carlo simulation model . . . . 3.2 Sequential Monte Carlo estimation of collision risk . . . . . . 3.2.1 Stochastic hybrid process considered . . . . . . . . . . 3.2.2 Risk factorization using multiple conflict levels . . . . 3.2.3 Characterization of the risk factors . . . . . . . . . . . 3.2.4 Interacting Particle System based risk estimation . . . 26 26 26 27 28 29 30 30 31 33 33 1 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Contents 3.3 3.4 3.5 3.2.5 Modification of IPS resampling step 4 . . . . . . . . . . . . Development of a Petri net model of free flight . . . . . . . . . . . 3.3.1 Specification of Petri net model . . . . . . . . . . . . . . . . 3.3.2 High level interconnection arcs . . . . . . . . . . . . . . . . 3.3.3 Agents and LPNs to represent AMFF operations . . . . . . 3.3.4 Interconnected LPNs of ASAS . . . . . . . . . . . . . . . . 3.3.5 Interconnected LPNs of “Pilot Flying” . . . . . . . . . . . . 3.3.6 Verification, parameterization, and validation of simulation model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.7 Dimensions of Monte Carlo simulation model . . . . . . . . Simulated scenarios and collision risk estimates . . . . . . . . . . . 3.4.1 Parameterization of the IPS simulations . . . . . . . . . . . 3.4.2 Eight aircraft on collision course . . . . . . . . . . . . . . . 3.4.3 Free flight through an artificially constructed airspace . . . 3.4.4 Reduction of the aircraft density by a factor four . . . . . . 3.4.5 Discussion of IPS simulation results . . . . . . . . . . . . . Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . A Mathematical Background A.1 Stationary Point Processes, Palm measure and applications to A.1.1 Point process . . . . . . . . . . . . . . . . . . . . . . . A.1.2 Palm measure . . . . . . . . . . . . . . . . . . . . . . . A.1.3 Residual waiting time and spent waiting time . . . . . A.1.4 Application to the collision risk . . . . . . . . . . . . . A.2 Stochastic approach of the collision risk . . . . . . . . . . . . A.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . A.2.2 Errors Modelling . . . . . . . . . . . . . . . . . . . . . A.2.3 The Kramers equation . . . . . . . . . . . . . . . . . . A.2.4 Solution for a linear force . . . . . . . . . . . . . . . . A.2.5 Randomization approach . . . . . . . . . . . . . . . . A.2.6 Extreme cases . . . . . . . . . . . . . . . . . . . . . . A.2.7 Low-Friction Limit . . . . . . . . . . . . . . . . . . . . A.2.8 High Friction . . . . . . . . . . . . . . . . . . . . . . . A.3 Rare Event Simulations . . . . . . . . . . . . . . . . . . . . . A.3.1 Multilevel Feynman-Kac distributions . . . . . . . . . A.3.2 Interacting particle system approximations . . . . . . 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 36 36 37 38 40 43 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 46 47 47 48 49 51 52 52 risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 54 54 55 55 56 60 60 62 65 66 67 68 69 74 76 77 79 collision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Introduction 1.1 Objectives The objective of work package WP7 is to assess the A3 operations developed by work packages WP1 and WP2, through hazard identification and Monte Carlo simulation on accident risk as a function of traffic demand, to assess what traffic demand can safely be accommodated by this advanced operational concept, and to assess the efficiency of the flights. The accident risk levels assessed should be in the form of an expected value, a 95% uncertainty area, and a decomposition of the risk level over the main risk contributing sources. The latter verifies which of these sources should have been mitigated during the 2nd design cycle by the A4 concept. In order to accomplish this assessment through Monte Carlo simulation, the complementary aim of this WP is to further develop the innovative HYBRIDGE speed up approaches in rare event Monte Carlo simulation. The aim of the current report is to give a thorough review of risk assessment status for advanced air traffic operations. 1.2 Accident risk assessment in ATM Among the class of complex and safety critical industries, air traffic is an interesting example that poses exceptional challenges to advanced design. By its very nature, each aircraft has its own crew, and each crew is communicating with several human operators in different air traffic management (ATM) and airline operational control (AOC) centres on the ground in order to timely receive instructions critical to a safe flight. In addition, from an organisational perspective, air traffic involves interactions between many stake holders: pilots, air traffic controllers, airline operation centres, airport authorities, government regulators and the public travelling. Figure 1.2 highlights this characteristic feature of interplay between distributed decision making and safety criticality both for air traffic and for other complex or safety-critical industries, such as finance and nuclear and chemical plants. Among the safety critical industries, air traffic stands out regarding the many distributed levels of interactions in control and decision making. The implication is that safety of air traffic is the result of interactions between multiple human operators, procedures (including spacing and separation criteria), and technical systems (hardware and software) all of which are highly distributed. Since safety depends crucially on the interactions between the various elements of a system, providing safety is more than making sure that each of these elements function properly. It is imperative 3 1 Introduction Potential number of fatalities per accident Thousands Hundreds Tens Few Localised interactions Distributed interactions Highly distributed interactions Figure 1.1: Air traffic compared with other complex and/or safety-critical industries in terms of potential number of fatalities per accident and the level of distributed interactions to understand the safety impact of these interactions, particularly in relation to nonnominal situations. Traditional ATM design approaches tend to be bottom-up, that is starting from developing concept elements aimed at increasing capacity, and next to extend the design with safety features. The advantage of the traditional approach is that advanced design developments can be organised around the clusters of individual elements, i.e., the communication cluster, the navigation cluster, the surveillance cluster, the automation tools cluster, the controllers/pilots and their human machine interfaces (HMIs), the advanced procedures, etcetera. The disadvantage of this traditional approach is that it fails to fully address the impact of interactions between controllers, pilots and ATM systems on safety. A goal oriented approach would be to design ATM such that safety has been built in at the capacity-level required. From this perspective, safety assessment forms a primary source of feedback (Figure 1.2) in the development of advanced ATM designs. An early guidance of ATM design development on safety grounds can potentially avoid a costly redevelopment program, or an implementation program that turns out to be ineffective. Although understanding this idea is principally not very difficult, it can be brought into practice only when an ATM safety assessment approach is available that provides appropriate feedback to the ATM designers from an early stage of the concept development 4 1 Introduction (Figure 1.2). This feedback should provide information on which safety-capacity issues are the main contributor to unsafety. Safety/Capacity Assessment ATM design Figure 1.2: Safety feedback based ATM design Collision risk modelling for air transportation was initially developed in the 1960s to address the safety of proposed separation standard reductions in the North Atlantic Organized Track Structure. The civil aviation community has developed a mathematical model to estimate mid-air collision risk levels as a function of spacing values in route structures [56]. This model is known as the Reich collision risk model; it assumes that the physical shape of each aircraft is a box, having a fixed orientation, and the collision risk between two aircraft is approximated by an expression that has proven to be of practical use in designing route structures [54]. Apart from the approximation, the most severe shortcoming is that the Reich model does not adequately cover situations where interaction between pilots and controllers play a crucial role, e.g. when controllers monitor the air traffic through surveillance systems and provide tactical instructions to the aircraft crews. These shortcomings of the methodology have led to the development of new approaches and relevant tools for the safe separation assessment of advanced procedures in air traffic, e.g Analytic Blunder Risk Model (ABRM), Airspace Simulation and Analysis for Terminal instrument procedures (ASAT), ICAO’s Collision Risk Model (CRM), Reduced Aircraft Separation Risk Assessment Model (RASRAM) [89] and Traffic Organization and Perturbation AnalyZer (TOPAZ) [9]. An extensive overview of these approaches can be found in [26]. TOPAZ appeared to be the most advanced in going beyond established approaches. It was developed by National Aerospace Laboratory NLR as a safety risk assessment methodology which provides safety risk feedback to advanced air traffic operation design. Within TOPAZ, Petri net modelling and Monte Carlo simulation has proven to deserve a key role in modelling and assessment of advanced air traffic operations on safety risk [5, 10, 13, 14, 15, 16, 38, 39, 41, 90]. In this respect it is relevant to notice that the use of Petri nets has been shown to work well in modelling safety critical operations in nuclear and chemical industries (e.g. [72]). In this report we will explain in details the TOPAZ methodology and present recent extensions which were developed during the HYBRIDGE project. 5 1 Introduction 1.3 Organisation of this report This report is organized as follows. The aim of Chapter 2 is to explain how the TOPAZ methodology effectively uses Monte Carlo simulation in safety risk assessment of an advanced air traffic operation. Emphasis is on how Monte Carlo simulation of safety risk works and how this is embedded within a complete safety risk assessment cycle. First, Section 2.2 provides an overview of the steps of the TOPAZ safety risk assessment cycle and for which step Monte Carlo simulation is of direct use. Next, Section 2.3 provides an overview of how to develop a Monte Carlo simulation model of a given operation. In order to keep the explanation concrete, a particular example is introduced first. Subsequently, Section 2.4 provides an overview of key issues that have to be taken into account when using a Monte Carlo simulation supported safety risk assessment. Section 2.5 presents Monte Carlo simulation results for the particular example identified in Section 2.3. Finally, conclusions are drawn in Section 2.6. Chapter 3 presents novel sequential Monte Carlo simulation methods which were developed within HYBRIDGE project for a much efficient estimation of collision risk in advanced ATM scenarios. The approach is demonstrated on application of these novel Monte Carlo techniques to a free flight air traffic concept of operations. Section 3.2 develops the sequential Monte Carlo simulation approach toward probabilistic reachability analysis of a Generalised Stochastic Hybrid System (GSHS) model of free flight air traffic. Section 3.3 explains how an initial GSHS model has been developed for a specific free flight air traffic concept of operation. Section 3.4 applies the sequential Monte Carlo simulation approach of Section 3.2 to the GSHS model of Section 3.3. Section 3.5 draws conclusions. 6 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations 2.1 Introduction This chapter gives an overview of performing safety risk assessment of a safety critical operation with support of Monte Carlo simulation. The approach is outlined for an air traffic example involving aircraft departing from a runway, which is occasionally crossed by taxiing aircraft. At the airport considered, a Runway Incursion Alert System (RIAS) is installed to warn the air traffic controller in case of impending runway incursions. The chapter explains the key issues to be mastered in performing a Monte Carlo simulation supported safety risk assessment of this kind of operation. To begin with, one has to develop an appropriate simulation model, and a sound way to speed up the Monte Carlo simulation based on this model. Complementary, one has to validate the simulation model versus the real operation, and the simulation supported approach has to be embedded within the safety risk assessment of the total operation. For this application example Monte Carlo simulation results are given and the way of feedback to the design of the operation is outlined. 2.2 Safety Risk Assessment Steps An overview of the steps in a TOPAZ safety risk assessment cycle is given in Figure 2.1. Although the cycle itself is very much in line with the established risk assessment steps (e.g. [70]), some of these steps differ significantly. In step 0, the objective of the assessment is determined, as well as the safety context, the scope and the level of detail of the assessment. The actual safety assessment starts by determining the operation that is assessed (step 1). Next, hazards associated with the operation are identified (step 2), and aggregated into safety relevant scenarios (step 3). Using severity and frequency assessments (steps 4 and 5), the safety risk associated with each safety relevant scenario is classified (step 6). For each safety relevant scenario with a (possibly) unacceptable safety risk, the main sources contributing to unsafety (safety bottlenecks) are identified (step 7), which help operational concept developers to learn for which safety issues they should develop improvements in the ATM design. If the 7 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations 2 Identify hazards 3 Construct scenarios 1 Determine operation 0 Identify 4 severities Identify objective Iterate (option) Assess 5 frequency Operational development Assess risk Decision making 6 tolerability Identify safety 7 bottlenecks Figure 2.1: Steps in TOPAZ safety risk assessment cycle ATM design is changed, a new safety risk assessment cycle of the operation should be performed in order to investigate how much the risk posed by previous safety issues has been decreased, but also to assess any new safety issues that may have been introduced by the enhancements themselves. The following subsections present the risk assessment steps of a TOPAZ cycle in more detail. Then it also becomes clear that Monte Carlo simulation plays a key role in step 5: assess frequency. Step 0: Identify objective Before starting the actual safety assessment, the objective and scope of the assessment are set. This should be done in close co-operation with the decision makers and designers of the advanced operation. Also, the safety context must be made clear, such that the assessment is performed in line with the appropriate safety regulatory framework. An important issue for setting the safety context is the choice of safety criteria with respect to which the assessment is performed. Depending of the application, such criteria are defined for particular flight condition categories (e.g. flight phases or sub-phases) and for particular severity categories (e.g. accident, serious incident). Typically, within the chosen context, these criteria define which flight condition/severity categories have to be evaluated and which frequency level forms the Target Level of Safety (TLS) threshold per flight condition/severity category. 8 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations Step 1: Determine operation Step 1 serves for the safety assessors to obtain a complete and concise overview of the operation, and to freeze this description during each safety assessment cycle. Main input to step 1 is a description of the operational concept from the designers, while its output is a sufficiently complete, structured, consistent and concise description of the operation considered. The operation should be described in generic terms, the description should provide any particular operational assumption to be used in the safety assessment, and the description has to be verified by the operational concept designers. Typically during this step, holes and inconsistencies in the concept as developed are also identified and immediately fed back to the design team. Step 2: Identify hazards The term hazard is used in the wide sense; i.e. an event or situation with possibly negative effects on safety. Such a non-nominal event or situation may evolve into danger, or may hamper the resolution of the danger, possibly in combination with other hazards or under certain conditions. The goal of step 2 is to identify as many and diverse hazards as possible. Hazard identification brainstorming sessions are used as primary means to identify (novel) hazards. In system engineering, the functional approach to hazard identification is well-known. In this approach it is attempted to determine all possible failure conditions and their effects, for each function that plays a role in the operation, including the human operator tasks. Unfortunately, the approach cannot identify all hazards related to an operation that involves human operators. An important reason for this is that the performance of air traffic controllers and pilots depend on their (subjective) situational awareness. From a human cognition perspective a particular act by an air traffic controller or pilot can be logical, while from a function allocation perspective the particular act may be incorrect. Such incidents are often called “errors of commission” [91]. An example of error of commission in the crossing operation is that because of the complicated taxiway structure, the pilot thinks to be taxiing far from the runway, while in reality, he starts crossing the runway without noticing any of the runway signs. Another well-known technique of hazard identification is the HAZOP (HAZard and OPerability) method. With this method, hazards are identified and analyzed using sessions with operational experts. At the same time, the experts come up with potential solutions and measures to cope with the identified hazards [63]. The advantage of HAZOP with respect to the functional approach is that also non-functional hazards are identified during the brainstorm with operational experts. However, in applying HAZOP, one needs to take care that hazard analysis and solution activities do not disturb the hazard identification process, which could leave certain hazards unidentified or inappropriately “solved”. Leaving such latent hazards in a design typically is known to be very costly in safety critical operation. 9 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations Based on the experience gained in using the hazard identification part of HAZOP in a large number of safety analyses and on scientific studies of brainstorming, NLR has developed a method of hazard identification for air traffic operations – by means of pure brainstorming sessions [60]. In such a session no analysis is done and solutions are explicitly not considered. An important complementary source is formed by hazards identified in previous studies on related operations. For this purpose, hazards identified in earlier studies are collected in a TOPAZ database. Step 3: Construct scenarios When the list of hazards is as complete as reasonably practicable, it is processed to deal with duplicate, overlapping, similar and ambiguously described hazards. First, per flight condition selected in Step 0, the relevant scenarios which may result from the hazards are to be identified using a full list of potentially relevant scenarios, such as for instance ‘conflict between two aircraft merging onto one route’ or ‘aircraft encounters wake vortex of parallel departure’. Per flight condition, each potentially relevant scenario is subsequently used as crystallisation point upon which all applicable hazards and their combined effects are fitted. If hazards are not appropriately addressed by the crystals developed so far, then additional crystallisation points are defined. The output of such crystallisation process is a bundle of event/condition sequences and effects per crystallisation point, and these are referred to as a safety relevant scenario. This way of constructing scenarios aims to bring into account all relevant ways in which a hazard can play a role in each flight condition/severity category. In order to cope with the complexity of the various possible causes and results, clusters of similarly crystallised hazards are identified. A cluster of hazards could for instance be the set of ‘events causing a missed approach to deviate from the normal path’. An example is given in Figure 2.2. It should also be noticed that the same cluster of hazards may play a role in one or more safety relevant scenarios. Each of the identified hazards can be of the following types: • a root hazard (cluster), which may cause a safety relevant scenario; or • a resolution hazard (cluster), which may complicate the resolution of a safety relevant scenario. For an appropriate safety risk assessment, all combinations of root and resolution hazards have to be evaluated in the next steps. Step 4: Identify severities For each of the safety relevant scenarios identified in step 3, it is determined which of the severity categories selected in step 0 are applicable to its possible effects. Safety experts should assess which of the severities are applicable for each safety relevant scenario, by 10 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations Figure 2.2: Example of a safety relevant scenario diagram consultation of and review by operational experts. For each safety relevant scenario the effects and their severities depend on many factors, such as the conditions under which the scenario starts and evolves, the geometry of the scenario, and the possibilities of (timely) resolution of the conflict. Therefore, a range of severities may apply to a safety relevant scenario. If necessary, the structuring of the events in the safety relevant scenarios of step 3 are updated such that each applicable severity category is linked to the occurrence of specific event sequences. Step 5: Assess frequency Next, for each possible severity outcome of each safety relevant scenario, the occurrence frequency is evaluated by making use of an appropriate tree per safety relevant scenario. The probability of the top event in the tree is expressed as a sum of a product of probabilities of applicable conditional events. For assessing the factors in these trees, primary sources of data are operational experts and databases. Examples of databases are aviation safety databases, local controller reporting system(s), et cetera. For appro- 11 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations priate use of such data dedicated operational expertise is taken into account. Hence, important input for the frequency assessments is always formed by interviews with operational experts (including experienced pilots and controllers) who are as much as is possible familiar with the operation under consideration. Qualitative expressions are to be translated in quantitative numbers when the selected safety criteria of step 0 also are expressed in numbers. Complicating factors in assessing at once the frequency of a conflict ending in a given severity can be that there is often little or no experience with the new operation, and that the situation may involve several variables. This holds especially for the more severe outcomes of a safety relevant scenario, since these situations occur rarely, and consequently little information is at hands about the behaviour of air traffic controllers and pilots in such situations. For these difficult safety relevant scenarios it is logical to make use of Monte Carlo simulation of safety risk. This approach has three clear advantages: 1) the quality of the risk estimate improves; 2) it is possible to estimate a 95% confidence interval; and 3) once a Monte Carlo simulation tool for a particular application has been developed it can be re-used for assessing safety risk for similar applications. The next sections explain for an example safety risk assessment by Monte Carlo simulation. Step 6: Assess risk tolerability The aim of this step is to assess the tolerability of the risk for each of the flight condition/severity categories selected in step 0. First the total risk per flight condition/severity category is determined by summing over the assessed risk contributions per safety relevant scenario for that flight condition/severity category. This summation takes into account both the expected value and the 95% confidence interval of the risk summation. Next for each flight condition/severity category the total risk expected value and the 95% confidence interval are compared against the TLS selected in step 0. Step 7: Identify safety bottlenecks From the risk tolerability assessment, it follows which safety relevant scenario(s) contribute(s) most to the expected value and the 95% confidence interval of the risks that has been qualified as being not below the TLS. For each safety relevant scenario the hazards or conditions that contribute most to the high risk level or confidence interval are identified and localised during step 7. These are referred to as the safety bottlenecks. If desired, this may also be done for assessed risk levels that are just below the TLS. The identification and localisation of safety bottlenecks is important as it gives operational concept designers directions for searching potential risk mitigating measures of the operation, and it gives the safety assessment experts the hazards and conditions for which the reduction of uncertainty has priority. 12 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations 2.3 Monte Carlo Simulation Model 2.3.1 Active Runway Crossing Example The Monte Carlo simulation-based risk assessment approach will be illustrated for an active runway crossing operation. This example accounts for a number of interacting human agents (pilots and controllers). The runway configuration of the active runway crossing operation considered is shown in Figure 2.3. The configuration takes into account one runway, named runway A, with holdings for using the runway from two sides (A1 and A2) and with crossings (C1, C2, D1 and D2) and exits (E1, E2, E3 and E4). The crossings enable traffic between the aprons and a second runway, named runway B. Each crossing has remotely controlled stopbars on both sides of the runway. Also the holdings have remotely controlled stopbars and each exit has a fixed stopbar. Holding Aprons A1 E1 C1 Runway B Aprons C2 E2 E3 D1 Aprons Runway B D2 E4 A2 Holding Aprons Figure 2.3: Runway configuration of active runway crossing procedure 13 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations The involved human operators include the start-up controller, the ground controller, the runway A controller, the runway B controller, the departure controller, and the pilots flying (PF’s) and pilots not flying (PNF’s) of aircraft taking off and aircraft crossing. Communication between controllers and aircraft crews is via standard VHF R/T (Very High Frequency Receiver/Transmitter). Monitoring by the controllers can be by direct visual observation under sufficiently good visibility conditions; it is supported by ground radar surveillance. The runway A controller is supported by a runway incursion alert system and a stopbar violation alert system. The runway A controller manages the remotely controlled stopbars and the runway lighting. Monitoring by the aircraft crews is by visual observation, supported by the VHF R/T party-line effect. In the runway crossing operation considered, the control over the crossing aircraft is transferred from the ground controller or the runway B controller (depending on the direction of the runway crossing operation) to the runway A controller. If the runway A controller is aware that the runway is not used for a take-off, the crew of an aircraft intending to cross is cleared to do so and subsequently the appropriate remotely controlled stopbar is switched off. The PNF of the crossing aircraft acknowledges the clearance and the PF subsequently initiates the runway crossing. When the crossing aircraft has vacated the runway, then the PNF reports this to the runway A Controller. Finally, the control over the aircraft is transferred from the runway A controller to either the runway B controller or the ground controller. 2.3.2 Safety Relevant Scenarios Prior to the development of a quantitative accident risk model for the active runway crossing operation considered, all risk assessment steps had been performed using an expert-based approach. In this study the following safety relevant scenarios were found: • Scenario I: Aircraft erroneously in take-off and crossing aircraft on runway; • Scenario II: Aircraft erroneously crossing and other aircraft in take-off; • Scenario III: Aircraft taking off and runway unexpectedly occupied; • Scenario IV: Aircraft crossing and runway unexpectedly occupied by aircraft; • Scenario V: Aircraft crossing and vehicle on runway; • Scenario VI: Collision between aircraft sliding off runway and aircraft near crossing; • Scenario VII: Aircraft taking off and vehicle crossing; • Scenario VIII: Jet-blast from one aircraft to another; • Scenario IX: Conflict between aircraft overrunning/climbing out low and aircraft using a nearby taxiway. From this expert-based study it followed that of all identified safety relevant scenarios, for scenarios I, II and III it was difficult to assess the risk sufficiently accurate using an 14 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations expert based approach. For these three scenarios it is therefore useful to assess the risk through performing Monte Carlo simulations. In this chapter, we focus on the details of a Monte Carlo simulation accident risk model for scenario II. In this scenario there is one aircraft that takes off and has been allowed to do so and there is one aircraft that crosses the runway while it should not. Taxiing along a straight line over one of the standard runway crossings (i.e., via C1, C2, D1 or D2 in Figure 2.3) is considered. 2.3.3 Multi-Agent Situation Awareness in the Simulation Model The safe organisation of co-operation between pilots and controllers in air traffic depends to a large extent on the “picture” or situation awareness (SA) maintained by each of the pilots and controllers. When a difference, even a small one, sneaks into the individual pictures and remains unrecognised, this may create unnoticed miscommunication and a subsequent propagation and increase in differences between the individual pictures. Eventually the situation may spiral out of control, with potentially catastrophic results. Hence any mismatch between individual pictures forms a serious hazardous condition in maintaining a safe organisation. Many hazards identified for the runway crossing operation were of this type. Endsley [36] has defined human SA as the perception of elements in the environment within a volume of time and space, the comprehension of their meaning, and the projection of their status in the near future. Stroeve et al. [90] and Blom and Stroeve [15] have captured these perception, comprehension and projection notions of SA mathematically in terms of three components: State SA, Mode SA and Intent SA. They also extended this single (human) agent SA concept to a multi-agent SA concept for operations involving multiple humans and systems, inclusive the basic updating mechanisms of such multi-agent SA. As depicted in Figure 2.4, for the active runway crossing operation we identified a need to model 10 agent types (7 humans and 3 systems) and their interactions: • Pilots flying; • Pilots not flying; • (Each) aircraft; • Aircraft’s flight management systems (FMS); • Runway A controller; • Runway B controller; • Ground controller; • Departure controller; • Start-up controller; 15 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations Aircraft iS Pilot Flying iS Pilot Flying iT Aircraft iT FMS iS A Pilot Not Flying iS Pilot Not Flying iT FMS iTA Ground ATCo RWY A ATCo Start Up ATCo RWY B ATCo Departure ATCo ATC System Figure 2.4: Relations between agents identified for the active runway crossing operation • ATC system, which we broadly define to include airport manoeuvre control systems, air traffic communication and surveillance systems, airport configuration and environmental conditions. 2.3.4 Dynamic Stochastic Modelling The Monte Carlo simulations are based on dynamic stochastic models of all relevant agents. These simulation models are mathematically specified using the Dynamically Coloured Petri Net (DCPN) formalism [39, 41]. A high-level overview of the agents modelled is provided next. Taking-off Aircraft The model of the taking-off aircraft represents the ground run, airborne transition and airborne climb-out phases and includes the possibility of a rejected take-off. The takingoff aircraft initiates its take-off from a position near the runway threshold and may have a small initial velocity. The aircraft may have diminished acceleration or deceleration power. Two types of aircraft are included in the model: medium-weight aircraft and heavy-weight aircraft. Taxiing Aircraft The model of the taxiing aircraft represents aircraft movements (hold, acceleration, constant speed, deceleration) during taxiing. The taxiing aircraft enters the taxiway 16 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations leading to a runway crossing at a position close to the remotely controlled stopbar, with a normal taxiing speed or initiates taxiing from stance. The entrance time of the crossing aircraft is uniformly distributed around the take-off start time. The taxiing aircraft may have diminished deceleration power. Two types of aircraft are included in the model: medium-weight aircraft and heavy-weight aircraft. Pilot Flying of Taking-off Aircraft Initially, the pilot flying (PF) of a taking-off aircraft has the SA that taking-off is allowed and initiates a take-off. During the take-off the PF monitors the traffic situation on the runway visually and via the VHF communication channel. The PF starts a collision avoiding braking action if a crossing aircraft is observed within a critical distance from the runway centre-line or in reaction to a call of the controller, and if it is decided that braking will stop the aircraft in front of the crossing aircraft. Further details of taking-off aircraft PF model are given by [90]. Pilot Flying of Taxiing Aircraft Initially, the PF expects that the next airport way-point is either a regular taxiway or a runway crossing. In the former case the PF proceeds taxiing and in the latter case the PF may have the SA that crossing is allowed. The characteristics of the visual monitoring process of the PF depend on the intent SA. In case of awareness of a conflict, either due to own visual observation or due to a controller call, the PF stops the aircraft, unless it is already within a critical distance from the runway centre-line. Further details of taxiing aircraft PF model are given by [90]. Runway Controller The runway A controller visually monitors the traffic and has support from a stopbar violation alert and a runway incursion alert. If the controller is aware that a taxiing aircraft has passed the stopbar, a hold clearance is given to both taxiing and taking off aircraft. Further details of the runway controller model are given by [90]. Radar Surveillance System The model of the radar surveillance system represents position and velocity estimates for both aircraft. There is a probability that radar surveillance is not available, resulting in track loss. Radar surveillance data is used as basis for ATC stopbar violation alerting and ATC runway incursion alerting. ATC Alerts Two types of ATC alerts are included in the model: a stopbar violation alert and a runway incursion alert. A stopbar violation alert is presented to the controller if surveillance data indicates that an aircraft has passed an active stopbar. There is a probability that the stopbar violation alert system does not function, implying that there will be no alert. A runway incursion alert is presented to the controller if radar surveillance data indicates that the taxiing aircraft is within a critical distance of the runway centre-line and the taking-off aircraft has exceeded a velocity threshold in front of the runway crossing. There is a probability that the runway incursion alert system does not function, implying that there will be no alert. 17 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations VHF Communication Systems The model for the VHF communication system between the runway controller and the aircraft crews accounts for the communication system of the aircraft, the communication system of the controller, the tower communication system, the frequency selection of aircraft communication system and the VHF communication medium. The nominal status of these communication systems accounts for direct non-delaying communication. The model accounts for a probability of delay in or failure of the communication systems. 2.4 Use of Simulation Model in Risk Assessment Once the simulation model has been specified, there are several important aspects that have to be taken into account during the preparation, execution and interpretation of the Monte Carlo simulations. This section explains these aspects. 2.4.1 Does the Simulation Model Cover the Identified Hazards? During step 2 of the safety assessment cycle, a lengthy list of hazards, including nonnominal situations, has been identified. These hazards contribute individually and possibly in combination with other hazards to the safety risk of the operation considered. Hence it is quite important to verify prior to performing the simulations that the hazards identified in step 2 of the assessment cycle are covered by the model. The verification process consists of specifying per hazard how it is captured by the simulation model. A special class of hazards is formed by the situation awareness related hazards. Table 2.1 shows three of such situation awareness related hazards and includes a short explanation how these hazards are covered by the simulation model. Table 2.1: Examples of situation awareness related hazards and their simulation model. Hazard Pilots become confused about their location at the airport because of complexity of the airport layout. Crew of taxiing aircraft is lost and therefore not aware of starting to cross a runway. RIAS is switched off by maintenance and controllers are not informed. Simulation model State SA of the PF of a taxiing aircraft is that its aircraft is at a location that differs from the actual location. Intent SA of PF is and stays taxiing while PF starts crossing the runway. RIAS working or not is not connected to Mode SA of controllers. Inevitably this verification of each hazard against the model will lead to the identification of hazards that are not (yet) covered by the simulation model. For non-covered hazards the simulation model developers should consider to further extend the simulation model prior to performing Monte Carlo simulations. 18 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations 2.4.2 Parametrisation of the Simulation Model During the mathematical specification of the simulation model there is no need to bother about the correct parameter values to be used during the Monte Carlo simulation. Of course, this is addressed prior to running the simulations. In principle there are three kinds of sources for parameter values. The ideal source would consist of sufficient statistical data that has been gathered under the various contextual conditions for which the risk assessment has to be performed. In practice such ideal sources almost never exist. Instead one typically has to work with limited statistical data that has been gathered under different conditions. Fortunately there often are two complementary sources: domain expertise and scientific expertise (on safety and human factors). In the context of Monte Carlo simulation this means one fuses statistical and expertise sources into a probability density function for the possible values of each parameter. Typically the mean or mode of such a density function is then used as the best estimate of the parameter value to be used when running the Monte Carlo simulation. 2.4.3 Speeding up Monte Carlo Simulations Air traffic is a very safe means of transport. Consequently, the risk of collision between two aircraft is extremely low. The assessment of such low collision risk values through straightforward Monte Carlo simulation would need extremely lengthy computer simulation periods. In order to reduce this to practicable periods, five to six orders of magnitude in speeding up the Monte Carlo simulation are needed. The basis for realizing such speed-up factors in Monte Carlo simulation consists of decomposing accident risk simulations in a sequence of conditional Monte Carlo simulations, and then to combine the results of these conditional simulations into the assessed collision risk value. For the evaluation of logical systems good decomposition methods can often be obtained by Fault and Event Tree Analysis. Because air traffic operations involve all kinds of dependent, dynamic and concurrent feedback loops, these logic-based risk decomposition methods cannot be applied without adopting severe approximations, typically by assuming that events/entities happen/act independent of each other. The stochastic analysis framework, that has shown its value in many applications (e.g. in financial mathematics [50]), is exploited by the TOPAZ methodology to develop Monte Carlo simulation models and appropriate speed-up factors by risk decomposition. The power of these stochastic analysis tools lies in their capability to model and analyse in a proper way the arbitrary stochastic event sequences (including dependent events) and the conditional probabilities of such event sequences in stochastic dynamic processes [10, 65]. By using these tools from stochastic analysis, a Monte Carlo simulation based risk assessment can mathematically be decomposed into a well-defined sequence of conditional Monte Carlo simulations together with a subsequent composition of the total risk out of these conditional simulation results. The latter composition typically consists of a tree with conditional probabilities to be assessed at the leaves, and nodes which either 19 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations add or multiply the probabilities coming from the subbranches of that node. Within TOPAZ such a tree is referred to as a collision risk tree [9, 10]. For the active runway crossing example, the particular conditions taken into account for this risk decomposition are: • The type of each aircraft (either a medium-weight or a heavy-weight); • The intent SA of the PF of a crossing aircraft concerning the next way-point (Taxiway/Crossing) and concerning allowance of runway crossing (Allowed/Not Allowed); • The alert systems (functioning well or not); • The remotely controlled stopbar (functioning well or not); and • The communication systems (functioning well or not). Based on the simulation model and the accident risk decomposition, Monte Carlo simulation software is developed to evaluate the event probabilities and the conditional collision risks, and to compose this with the help of the collision risk tree into the collision risk value assessed for the simulation model. 2.4.4 Validation of the assessed risk level For operations as complex as the active runway example considered, a simulation model will always differ from reality. Hence, validation of the Monte Carlo simulation results does not mean that one should try to show that the model is perfect. Rather one should identify the differences between the simulation model and reality, and subsequently analyse what the effects of these differences are in terms of bias and uncertainty at the assessed risk level of the model. If the bias and uncertainty fall within acceptable bounds, then the assesed risk levels are valid for the specified application. Otherwise one should improve the Monte Carlo simulation model on those aspects causing the largest bias and uncertainty influence on the assessed risk level. Five types of differences between simulation model and the real operation can be distinguished [38]: • Numerical approximations; • Parameter values; • Assumptions on the model structure; • Non-covered hazards; • Differences between the real operational concept and the operational concept modelled. Thinking in terms of these differences makes it possible to consider the validation problem as a problem of making the differences specific, assessing each difference and its effect 20 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations on the collision risk, and subsequently decide if this is accurate enough (valid) or not (invalid) for the purpose aimed at. The effects of differences on the collision risk can mathematically be expressed in terms of bias and uncertainty that has to be taken into account when using the simulation model assessed risk value for decisions about reality: • Bias. The accident risk as defined by the simulation model is systematically higher or lower than it is for the real operation. • Uncertainty. In addition to a systematic bias, the differences between simulation model and reality may induce uncertainty in the difference between the safety risk of the real operation and the safety risk resulting from the simulation model. With this, the validation of a simulation based accident risk assessment has largely become a bias and uncertainty assessment process. Within TOPAZ, a bias and uncertainty assessment method has been developed which consists of the following steps: • Identify all differences between the simulation model and reality; • Assess how large these differences are, or how often they happen; • Assess the sensitivity (or elasticity) of the risk outcome of the simulation model to changes in parameter values; • Assess the effect of each difference on the risk outcome, using model sensitivity knowledge and complementary statistical and/or expert knowledge; • Combine the joint effects of all differences in bias and uncertainty factors, and compensate the risk value of the model with these bias and uncertainty factors. The result is an expected value of risk for the real operation, including a 95% confidence interval of other possible risk values. If the bias or the 95% confidence interval of the combined effects, or the bias and uncertainty of individual differences is too large, then these differences have to be taken into account in the decision making process regarding the acceptability and/or further design of the operation considered. 2.5 Monte Carlo Simulation Results This section presents collision risk results obtained by Monte Carlo simulation with a computer implementation of the mathematical model of the active runway example of Section 2.3. In order to relate these results to an actual operation, a bias and uncertainty assessment remains to be performed; however, this falls outside the scope of this chapter. 2.5.1 Assessed risk levels Figure 2.5 shows the accident risk as function of the position of the runway crossing with respect to the runway threshold. The probability of a collision decreases for positions 21 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations -7 Collision risk [occurrence / take-off] 10 Total SA = Proceed taxiway SA = Cross runway -8 10 -9 10 -10 10 500 1000 1500 Crossing distance [m] 2000 Figure 2.5: Contributions to the total collision risk by the simulation model for the cases that the SA of the PF of the taxiing aircraft is to proceed on a taxiway, or to cross the runway of the crossing distances further from the threshold. Figure 2.5 also shows the decomposition of the total risk for the cases that the pilot flying of the taxiing aircraft either thinks to be proceeding on a normal taxiway (without being aware to be heading to a runway crossing) or where the pilot intends to cross the runway (without being aware that crossing is currently not allowed). The largest contribution to the risk is from the situation that the pilot thinks to be proceeding on a normal taxiway. The relative size of this contribution depends on the crossing distance and varies from 64% for crossing at 500 m to about 83% for crossing at 1000 or 2000 m. A more complete overview of the contributions to the collision risk is provided by a projected version of the collision risk tree in Figure 2.6. It shows the contributions of events related to the situation awareness of the pilot of the taxiing aircraft (Cross runway/Proceed runway) and the functioning of ATC alert and communication systems (Up/Down). The collision risk results in the leaves of the tree are the product of the probability of the event combination indicated and the Monte Carlo simulation based collision risk given the event combination. The results in Figure 2.6 show that the risk is dominated by situations with a pilot flying of a taxiing aircraft having an erroneous situation awareness and the ATC alert and communication systems working nominally. 2.5.2 Who contributes to safety risk reduction? Based on results of the accident risk model, it is possible to attain insight in the accident risk reducing performance of involved human operators and technical systems. Table 2.2 shows conditional collision risks for the situation that an aircraft taxies towards a runway crossing at a distance of 1000 m from the runway threshold while the pilot has 22 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations Total collision risk Initial SA difference No <<7e -9 SA pilot taxiing a/c ATC alert systems Communication systems Cross runway Up 1.1e -9 Up 1.1e - 9 Down 2.7e -13 7.1e -9 Yes 7.1e -9 1.1e -9 Proceed taxiway 6.0e -9 Down 1.8e -13 Up Down 1.8e -13 1.4e -17 Up 6.0e -9 Up 6.0e -9 Down 4.1e -13 Down 1.4e -13 Up 1.4e -13 Down 2.2e -17 Figure 2.6: Projected version of the collision risk tree for the active runway crossing example, showing the contributions to the collision risk for various combinations of events related to pilot situation awareness and functioning of ATC alert and communication systems. The values are for a crossing distance of 1000 m the situation awareness to taxi on a normal taxiway. The conditional collision risks in Table 2.2 refer to cases where the model either does (‘yes’) or does not (‘no’) involve the indicated human operators actively monitoring for traffic conflicts. A risk reduction percentage is determined by comparing the conditional collision risk with the situation in which none of the human operators is actively monitoring. In this case, a collision is only avoided by the lucky circumstances that the taxiing aircraft just passes in front of or behind the taking-off aircraft (case 0 in Table 2.2). From the results in Table 2.2 a number of model-based insights into the operation can be attained: • It follows from case 1 that 99.8% of the accidents can be prevented by the combined effort of all human operators and alert systems. • It follows from a comparison of cases 1 and 5 that in the normal situation that all human operators are actively monitoring, ATC alert systems (runway incursion or stopbar violation) have a modest effect on the achieved risk. • It follows from a comparison of cases 1 and 4, and cases 5 and 8, that the risk reduction that can be achieved by the tower controller in addition to the risk reduction of both pilots is very small. • It follows from comparison of cases 1 and 3, and cases 5 and 7 that the pilot of the taxiing aircraft has the largest capability to prevent a collision in this context. Thus, resolution of the conflict is most likely to be by the human operator whose wrong situation awareness initiated the conflict. 23 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations Table 2.2: Risk reduction achieved in the simulation model by various combinations of involved human operators when the PF of a taxiing aircraft intends to proceed on a normal taxiway under good visibility (crossing is at 1000 m from runway threshold.) Case PF taxiing PF taking-off Runway Conditional Risk reaircraft aircraft controller collision duction risk 0 no no no 8.9 10-2 ATC alert systems on 1 yes yes yes 1.7 10-4 99.8% 2 yes no yes 4.0 10-4 99.6% 3 no yes yes 9.4 10-3 89.4% 4 yes yes no 2.3 10-4 99.7% ATC alert systems down 5 yes yes yes 2.2 10-4 99.8% 6 yes no yes 1.7 10-3 98.1% 7 no yes yes 1.1 10-2 87.9% 8 yes yes no 2.3 10-4 99.7% 2.5.3 Comparison against expert based results In the earlier conducted expert based safety risk assessment of the active runway crossing operation, it was concluded that both the pilots and the runway controller make large contributions to the prevention of a collision in the scenario aircraft erroneously crossing and other aircraft in take-off. In hindsight, it can be concluded that in the expert based safety risk assessment, the total effect of the pilots and the runway controller in preventing a collision turns out to be overestimated under good visibility condition. It is the simulation based approach that makes clear that although the runway controller identifies a good share of the conflicts, its contribution to timely conflict resolution is relatively small. One significant part of the instruction issued by the runway controller appears to concern conflicts that are already solved by the pilots. And another significant part of the instructions issued by the runway controller appear to arrive too late for the pilots to successfully avoid a collision. Because of this, the effective contribution by the runway controller towards reducing collision risk is relatively small. 2.6 Concluding remarks This chapter has given an overview of performing safety risk assessment and providing feedback to the design of advanced air traffic operations with support of Monte Carlo simulation. The motivation for developing such a Monte Carlo simulation approach towards safety risk assessment was the identified need for modelling stochastic dynamic 24 2 Safety Risk Assessment by Monte Carlo Simulation of Complex Safety Critical Operations events and interactions between multiple agents (humans and systems) in advanced air traffic operations. The distributed and dynamical interactions pose even greater challenges than those seen in, for instance, nuclear and chemical industries (e.g. [72]). The chapter has explained the key issues to be mastered in performing a Monte Carlo simulation supported safety risk assessment of air traffic operations, and how this fits within a full safety risk assessment cycle. The steps to be followed in developing an appropriate Monte Carlo simulation model has been outlined, including a short overview of multiagent situation awareness modelling, which plays a key role in the safe organization of cooperation between many pilots and controllers in air traffic. The chapter also has explained the need for using stochastic analysis tools in order to develop the necessary speed-up of the Monte Carlo simulations, and has shown a feasible way to validate the simulation model versus the real operation. This assessment approach has been applied to an air traffic example involving aircraft departing from a runway that is occasionally crossed by taxiing aircraft. The results obtained demonstrate the feasibility and value of performing Monte Carlo simulation in accident risk assessment for safety relevant scenarios that are difficult to assess expert based, because of many interacting agents. 25 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation 3.1 Introduction 3.1.1 Safety verification of free flight air traffic Technology allows aircraft to broadcast information about its own-ship position and velocity to surrounding aircraft, and to receive similar information from surrounding aircraft. This development has stimulated the rethinking of the overall concept for today’s Air Traffic Management (ATM), e.g., to transfer responsibility for conflict prevention from ground to air. As the aircrews thus obtain the freedom to select their trajectory, this conceptual idea is called free flight [87]. It changes ATM in such a fundamental way, that one could speak of a paradigm shift: centralized control becomes distributed, responsibilities transfer from ground to air, fixed air traffic routes are removed, and appropriate new technologies are brought in. Each aircrew has the responsibility to timely detect and solve conflicts, thereby assisted by navigation means, surveillance processing, and conflict resolution systems. Due to the potentially many aircraft involved, the system is highly distributed. This free flight concept idea has motivated the study of multiple operational concepts and implementation choices [47, 52, 58, 64, 78]. One of the key outstanding issues is the safety verification of free flight design, and in particular when air traffic demand is high. For en-route traffic, the International Civil Aviation Organization (ICAO) has established thresholds on the acceptable probability of a mid-air collision. Hence, the en-route free flight safety verification problem consists of estimating the collision probability of free flight operations, and subsequently comparing this estimated level with the ICAO established thresholds [48]. The civil aviation community also has established some approximate models to estimate (an upperbound of) the risk of collision between aircraft flying within a given parallel route structure [46, 54, 57]. Additional methods have been exploited to develop some valuable extensions of this approach, e.g., using fault trees see [34] and using stochastic analysis and Monte Carlo simulation [8, 10, 43]. Andrews et al. [2] have shown how statistical data in combination with a fault tree of the functionalities of the advanced operation can serve to predict how reliability of free flight supported systems impact contributions to collision risk of an advanced operation [37], neglecting other contributions to collision risk. The challenge is to analyze the risk of collision between aircraft in free flight without the limitation of a fixed route structure. We aim 26 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation to improve this situation by developing a novel approach toward collision risk assessment for advanced air traffic designs. An initial shorter paper on this development is [11]. 3.1.2 Probabilistic reachability analysis In air traffic, a mid-air collision event happens at the moment in time that the physical shapes of two airborne aircraft hit each other. Such event can be represented as a moment in time that the joint state of aircraft involved hit a certain subset of their joint state space. With this, the problem to estimate the probability of collision between two aircraft within a finite time period is to analyze the probability that this collision subset is reached by their joint aircraft state within that time period. In systems theory, the estimation of the probability of reaching a given subset of the state space within a given time period is known as a problem of probabilistic reachability analysis, e.g., see [71]. Hu et al. [55] and Prandini and Hu [82, 83] apply probabilistic reachability analysis for the development of a grid based computation to evaluate the probability that two aircraft come closer to each other than some established minimum separation criteria. The numerical challenge of this problem, however, differs from the free flight collision risk estimation problem on the following aspects: • The collision subset is more than three orders smaller in volume than the conflict subset is. • A safety directed model of a free flight operation includes per aircraft also the states of the technical systems and the pilot models, which increases the size of the state space by many orders in magnitude. • In free flight there are multiple aircraft, not just two, that can trigger domino effects. If we would follow the numerical approach of Prandini and Hu [83] to estimate collision risk in free flight operations, then these aspects would imply a blow up of the number of grid points to a practically unmanageably large number. In most safety critical industries, e.g., nuclear, chemical, etc., reachability analysis is addressed by methods that are known as dynamical approaches towards Probabilistic Risk Analysis (PRA). For an overview of these dynamical methods in PRA, see [72]. These dynamical PRA methods make explicitly use of the fact that between two discrete events the dynamical evolution satisfies an ordinary differential equation. Essentially this means that these dynamical PRA methods apply to the class of stochastic hybrid system models that do not involve Brownian motion. In the hybrid systems control community these are known as piecewise deterministic Markov process [18, 29]. For proper safety modeling of air traffic operations, however, it is needed to incorporate Brownian motion in the piecewise deterministic Markov process models, e.g., to represent the effect of random wind disturbances on aircraft trajectories [81]. The class of systems which incorporates Brownian motion within piecewise deterministic Markov processes, 27 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation has been defined as a Generalized Stochastic Hybrid System (GSHS) [17]. GSHS is the class of non-linear stochastic continuous-time hybrid dynamical systems, having a hybrid state consisting of two components: a continuous valued state component and a discrete valued state component. The continuous state evolves according to an SDE whose vector field and drift factor depend on both hybrid state components. Switching from one discrete state to another discrete state is governed by a probability law or occurs when the continuous state hits a pre-specified boundary. Whenever a switching occurs, the hybrid state is reset instantly to a new state according to a probability measure which depends itself on the past hybrid state. GSHS contain, as a subclass, the switching diffusion process, the probabilistic reachability of which is studied in [83]. Important complementary dynamics is induced by the interaction between the hybrid state components. 3.1.3 Sequential Monte Carlo simulation Shah et al. [88] explain very well that the advantage of using Monte Carlo simulation in evaluating advanced operations is its capability to identify and evaluate emergent behavior, i.e., novel behavior which is exhibited by complex safety critical systems and emerges from the combined dynamical actions and reactions by individual systems and humans within the system. This emergent behavior typically cannot be foreseen and evaluated by examining the individuals’ behavior alone. Shah et al.[88] explain that agent based Monte Carlo simulation is able to predict the impact of revolutionary changes in air transportation; it integrates cognitive models of technology behavior and description of their operating environment. Simulation of these individual models acting together can predict the results of completely new transformations in procedures and technology. Their Monte Carlo simulations reach up to the level of novel emerging hazardous events. For safety risk assessment however, it is required to go further with the Monte Carlo simulations up to the level of emerging catastrophic events. In en-route air traffic these catastrophic events are mid-air collisions. A seemingly simple approach toward the estimation of mid-air collision probability is to run many Monte Carlo simulations with a free flight stochastic hybrid model and count the fraction of runs for which a collision occurs. The advantage of a Monte Carlo simulation approach is that this does not require specific assumptions or limitations regarding the behavior of the system under consideration. A key problem is that in order to obtain accurate estimates of rare event probabilities, say about 10−9 per flying hour, it is required to simulate 1011 flying hours or more. Taking into account that an appropriate free flight model is large, this would require an impractically huge simulation time. Del Moral and co-workers [23, 25, 31] developed a sequential Monte Carlo simulation approach for estimating small reachability probabilities, including a characterization of convergence behaviour. The idea behind this approach is to express the small probability to be estimated as the product of a certain number of larger probabilities, which can be 28 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation efficiently estimated by the Monte Carlo approach. This can be achieved by introducing sets of intermediate states that are visited one set after the other, in an ordered sequence, before reaching the final set of states of interest. The reachability probability of interest is then given by the product of the conditional probabilities of reaching a set of intermediate states given that the previous set of intermediate states has been reached. Each conditional probability is estimated by simulating in parallel several copies of the system, i.e., each copy is considered as a particle following the trajectory generated through the system dynamics. To ensure unbiased estimation, the simulated process must have the strong Markov property. Hence, we extend the approach of [23, 25] for application to free flight, and illustrate its application to free flight scenarios. 3.1.4 Development of Monte Carlo simulation model For the modeling of accident risk of safety-critical operations in nuclear and chemical industries, the most advanced approaches use Petri nets as model specification formalism, and stochastic analysis and Monte Carlo simulation to evaluate the specified model, e.g., see [72]. Since their introduction as a systematic way to specify large discrete event systems that one meets in computer science, Petri nets have shown their usefulness for many practical applications in different industries, e.g., see [28]. Various Petri net extensions and generalisations and numerous supporting computer tools have been developed, which further increased their modeling opportunities. Nevertheless, literature on Petri nets appeared to fall short for modeling the class of GSHS [17] that was needed to model air traffic safety aspects well [81]. Cassandras and Lafortune [22] provide a control systems introduction to Petri nets and a comparison with other discrete event modeling formalisms like automata. Both Petri nets and automata have their specific advantages. Petri net is more powerful in the development of a model of a complex system, whereas automata are more powerful in supporting analysis. In order to combine the advantages offered by both approaches, there is need for a systematic way of transforming a Petri net model into an automata model. Such a transformation would allow using Petri nets for the specification and automata for the analysis. For a timed or stochastic Petri net with a bounded number of tokens and deterministic or Poisson process firing, such a transformation exists [22]. In order to make the Petri net formalism useful in modeling air traffic operations, we need an extension of the Petri net formalism including a one-to-one transformation to and from GSHS. Everdij and Blom [40, 42] have developed such extension in the form of (Stochastically and) Dynamically Coloured Petri Net, or for short (S)DCPN. Jensen [59] introduced the idea of attaching to each token in a basic Petri net (i.e., with logic transitions only), a color which assumes values from a finite set. Tokens and the attached colors determine which transitions are enabled. Upon firing by a transition, new tokens and attached colors are produced as a function of the removed tokens and colors. Haas [51] extended this color idea to (stochastically) timed Petri nets where the time period between enabling and firing depends of the input tokens and their attached colors. 29 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation In [51, 59] a color does not change as long as the token to which it is attached remains at its place. Everdij and Blom [40, 41] defined a Dyamically Coloured Petri Net (DCPN) by incorporating the following extensions: (1) a color assumes values from a Euclidean state space, its value evolves as solution of a differential equation and influences the time period between enabling and firing; (2) the new tokens and attached colors are produced as random functions of the removed tokens and colors. An SDCPN extends an DCPN in the sense that colors evolve as solutions of a stochastic differential equation [42]. 3.2 Sequential Monte Carlo estimation of collision risk 3.2.1 Stochastic hybrid process considered Throughout this and the following sections, all stochastic processes are defined on a complete stochastic basis (Ω, F, F, P, T) with (Ω, F, P) a complete probability space, and F an increasing sequence of sub-σ-algebra’s on the positive time line T = R+ , i.e., ∆ F = {J , (Ft , t ∈ T), F}, J containing all P-null sets of F and J ⊂ Fs ⊂ Ft ⊂ F for every s < t. We assume that air traffic operations are represented by a stochastic hybrid process {xt , θt } which satisfies the strong Markov property. In [17, 19, 68, 69], this property has been shown to hold true for the processes generated as execution of a GSHS. For an N -aircraft free flight traffic scenario the stochastic hybrid process {xt , θt } consists of ∆ Euclidean valued components xt = Col{x0t , x1t , . . . , xN t } and discrete valued components ∆ θt = Col{θt0 , θt1 , . . . , θtN }, where xit assumes values from Rni , and θti assumes values from a finite set (M i ). Physically, {xit , θti }, i = 1, . . . , N , is the hybrid state process related to the i-th aircraft, and {x0t , θt0 } is a hybrid state process !N of all non-aircraft "N components. n The process {xt , θt } is R × M -valued with n = i=0 ni and M = i=0 Mi . In order to model collisions between aircraft, we introduce mappings from the Euclidean valued process {xt } into the relative position and velocity between a pair of two aircraft (i, j). The relative horizontal position is obtained through the mapping y ij (xt ), the relative horizontal velocity is obtained through the mapping v ij (xt ). The relative vertical position is obtained through the mapping z ij (xt ), and relative vertical rate of climb is obtained through the mapping rij (xt ). The relation between the position and velocity mappings satisfies: dy ij (xt ) = v ij (xt ) dt ij ij dz (xt ) = r (xt ) dt. (3.2.1) (3.2.2) A collision between aircraft (i, j) means that the process {y ij (xt ), z ij (xt )} hits the boundary of an area where the distance between aircraft i and j is smaller than their physical size. Under the assumption that the length of an aircraft equals the width of an aircraft, and that the volume of an aircraft is represented by a cylinder the orientation of which 30 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation does not change in time, then aircraft (i, j) have zero separation if xt ∈ Dij with: # $ Dij = x ∈ Rn ; |y ij (x)| ≤ (li + lj )/2 AND |z ij (x)| ≤ (si + sj )/2 , i %= j (3.2.3) where lj and sj are length and height of aircraft j. For simplicity we assume that all aircraft have the same size, by which (3.2.3) becomes: # $ Dij = x ∈ Rn ; |y ij (x)| ≤ l AND |z ij (x)| ≤ s , i %= j (3.2.4) Although all aircraft have the same size, notice that in (3.2.4), Dij still depends of (i, j). If xt hits Dij at time τ ij , then we say a collision event between aircraft (i, j) occurs at τ ij , i.e., τ ij = inf{t > 0 ; xt ∈ Dij }, i %= j (3.2.5) Next we define the first moment τ i of collision with any other aircraft, i.e., τ i = inf {τ ij } = inf {t > 0 ; xt ∈ Dij } = inf{t > 0 ; xt ∈ Di }, j"=i j"=i (3.2.6) ∆ with Di = ∪j"=i Dij . From this moment τ i on, we assume that the differential equations for {xit , θti } stop evolving. An unbiased estimation procedure of the risk would be to simulate many times aircraft i amidst other aircraft over a period of length T and count all cases in which the realization of the moment τ i is smaller than T . An estimator for the collision risk of aircraft i per unit T of time then is the fraction of simulations for which τ i < T . 3.2.2 Risk factorization using multiple conflict levels Cérou et al. [23, 25] developed a novel way of speeding up Monte Carlo simulation to estimate the probability that an Rn -valued strong Markov process xt hits a given “small” subset D ∈ Rn within a given time period (0, T ). This method essentially consists of taking advantage of an appropriately nested sequence of closed subsets of Rn : D = Dm ⊂ Dm−1 ⊂ . . . ⊂ D1 , and then start simulation from outside D1 , and subsequently simulate from D1 to D2 , from D2 to D3 , . . ., and finally from Dm−1 to Dm . Krystul and Blom [65, 66] extended this Interacting Particle System (IPS) approach to hybrid state strong Markov processes, and illustrated the effectiveness for a simple example. In order to apply this IPS approach to air traffic we first have to develop an appropriate sequence of nested subsets. Prior to a collision of aircraft i with aircraft j, a sequence of conflicts ranging from long term to short term always occurs. In order to incorporate this explicitly in the Monte Carlo simulation, we formalize this sequence of conflict levels through a sequence ij ij ⊂ Dm−1 ⊂ . . . ⊂ D1ij with for k = 1, . . . , m: of closed subsets of Rn : Dij = Dm # Dkij = x ∈ Rn ; |y ij (x) + ∆v ij (x)| ≤ dk AND $ |z ij (x) + ∆rij (x)| ≤ hk , for some ∆ ∈ [0, ∆k ] , (3.2.7) 31 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation for i %= j, with dk , hk and ∆k the parameters of the conflict definition at level k, and with dm = l, hm = s and ∆m = 0, and with dk+1 ≤ dk , hk+1 ≤ hk and ∆k+1 ≤ ∆k . If xt hits Dkij at time τkij , then we say the first level k conflict event between aircraft (i, j) occurs at moment τkij , i.e., τkij = inf{t > 0 ; xt ∈ Dkij }. (3.2.8) Similarly as we did for reaching the collision level by aircraft i, we consider the first moment τki that aircraft i reaches conflict level k with any of the other aircraft, i.e., τki = inf {τkij } = inf {t > 0 ; xt ∈ Dkij } = inf{t > 0 ; xt ∈ Dki }, j"=i j"=i (3.2.9) ∆ with Dki = ∪j"=i Dkij . Next, we define {0, 1}-valued random variables {χik , k = 0, 1, . . . , m} as follows: χik = 1, if τki < T or k = 0 = 0, else. By using this χik definition we can write the probability of collision of aircraft i with any of the other aircraft as a product of conditional probabilities of reaching the next conflict level given the current conflict level has been reached: )m + m * * ( % i & ' i ( ' i P τm < T = E χm = E χk = E χik | χik−1 = 1 k=1 = m * k=1 % & ∆ i <T . with γki = P τki < T | τk−1 % k=1 & i P τki < T | τk−1 <T = m * γki , (3.2.10) k=1 With this, the problem can be seen as one to estimate the conditional probabilities γki in such a way that the product of these estimators is unbiased. Because of the multiplication of the various individual γki estimators, which depend on each other, in general such a product may be heavily biased. The key novelty in [25] was to show that such a product may be evaluated in an unbiased way when {xt } makes part of a larger stochastic process that satisfies the strong Markov property. This approach is explained next. 3.2.3 Characterization of the risk factors Let us denote E ′ = Rn+1 × M , and let E ′ be the Borel σ-algebra of E ′ . For any B ∈ E ′ , ∆ πki (B) denotes the conditional probability of ξk = (τk , xτk , θτk ) ∈ B given χil = 1 for 1 ≤ l ≤ k. 32 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Define Qik = (0, T ) × Dki × M , k = 1, . . . , m. Then the estimation of the probability for ξk to arrive at the k-th nested Borel set Qik is characterized through the following recursive set of transformations prediction i πk−1 (·) conditioning pik (·) −→ −→ πki (·) ↓ γki where pik (B) is the conditional probability of ξk ∈ B given χil = 1 for 0 ≤ l ≤ k − 1. Because {xt , θt } is a strong Markov process, {ξk } is a Markov sequence. Hence the one step prediction of ξk satisfies a Chapman-Kolmogorov equation: , i i pk (B) = pξk | ξk−1 (B|ξ) πk−1 (dξ) for all B ∈ E ′ . (3.2.11) E′ Next we characterize the conditional probability of reaching the next level % & i γki = P τki < T | τk−1 <T ' ( = E χik | χik−1 = 1 , = 1{ξ∈Qi } pik (dξ), (3.2.12) k E′ and the conditioning satisfies: πki (B) - =- B 1{ξ∈Qi } pik (dξ) k i ′ E ′ 1{ξ ′ ∈Qi } pk (dξ ) k for all B ∈ E ′ . (3.2.13) With this, each of the m terms γki in (3.2.10) is characterized as a solution of a sequence of “filtering” kind of equations (3.2.11)–(3.2.13). However, an important difference with “filtering” equations is that (3.2.11)–(3.2.13) are ordinary integral equations, i.e., they have no stochastic term entering them. 3.2.4 Interacting Particle System based risk estimation Based on this theory, an Interacting Particle System (IPS) simulation algorithm is explained next for an arbitrary hybrid state strong Markov process model of air traffic. The i < T ), transformations (3.2.11)–(3.2.13) lead to the IPS algorithm of [25] to estimate P(τm i i i i i i where γ̄k , p̄k and π̄k denote the numerical approximations of γk , pk and πk respectively. i When simulating from Dk−1 to Dki , a fraction γ̄ki of the Monte Carlo simulated trajeci tories only will reach Dk within the time period (0, T ). IPS Step 0. Initial sampling for k = 0. 33 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation • For l = 1, . . . , Np generate initial state value outside Qi1 by independent drawings (xl0 , θ0l ) from px0 ,θ0 (·) and set ξ0l = (0, xl0 , θ0l ). • For l = 1, . . . , Np , set the initial weights: ω0l = 1/Np . !Np l • Then π̄0i = l=1 ω0 δ{ξl } . 0 IPS Iteration cycle: For k = 1, . . . , m perform step 1 (prediction), step 2 (assess fraction), step 3 (conditioning), and step 4 (resampling). i IPS Step 1. Prediction of πk−1 −→ pik , based on (3.2.11); • For l = 1, . . . , Np simulate a new path of the hybrid state Markov process, starting l at ξk−1 until the k-th set Qik is hit or t = T (the first component of ξkl counts time). Np l • This yields new particles {ξˆkl , ωk−1 . }l=1 • p̄ik is the empirical distribution associated with the new cloud of particles: p̄ik = !Np l l=1 ωk−1 δξ̂ l . k IPS Step 2. Assess fraction γki , based on (3.2.12); • The particles that do not reach the set Qik are killed, i.e., we set ω̂kl = 0 if ξˆkl ∈ / Qik l and ω̂kl = ωk−1 if ξˆkl ∈ Qik . !Np l • Approximation: γki ≈ γ̄ki = l=1 ω̂k . If all particles are killed, i.e., γ̄ki = 0, then i the algorithm stops without P(τ < T ) estimate. IPS Step 3. Conditioning of pik −→ πki , based on (3.2.13); The non-killed particles form a set Ski , i.e., iff ξˆkl ∈ Qik , then particle {ξˆkl , ω̂kl } is stored in Ski . NS Renumbering the particles in Ski yields a set of particles {ξ˜kl , ω̃kl }l=1k with NSk the number !NS of particles in Ski . Hence, we also have γ̄ki = l=1k ω̃kl . IPS Step 4. Resampling of πki Draw Np particles ξkl independently from the empirical measure π̄ki = of which gets weight ωkl = 1 Np . !NSk l=1 ω̃kl δ{ξ̃l } each k N p After step 4, the new set of particles is {ξkl , ωkl }l=1 .. If k < m then repeat steps 1, 2, 3, i i 4 for k := k + 1. Otherwise, stop with P(τ < T ) ≈ m k=1 γ̄k . Remark 3.2.1 Cérou et al. [23, 25] have proven, under certain conditions, how to i manage the simulations from Dk−1 to Dki , such that the product of these fractions γ̄ki forms an unbiased estimate of the probability of xt to hit the set Di within the time period (0, T ), i.e., )m + m * * i E γ̄k = γki = P(τ i < T ), k=1 k=1 34 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation and there also is some bound on the expected estimation error, i.e., / E( m * γ̄ki k=1 − m * 01 q γki )q k=1 aq bq ≤1 , Np for some finite constants aq and bq , which depend on the simulated scenario and the sequence of intermediate levels adopted. 3.2.5 Modification of IPS resampling step 4 A well known problem with particle systems is the possibility of particle depletion or impoverishment. In order to reduce the sensitivity of the above algorithm on these points, we modify step 4 in two ways: (1) we reduce the chance of impoverishment by not throwing away any particle; and (2) we make copies of particles, but avoid that these copies take away too much weight from the original particles. Modified IPS Step 4 Resample Np particles from Ski as follows: NS If 21 Np ≤ NSk ≤ Np , then copy the NSk particles, i.e., ξkl = ξ˜kl and set ωkl = ω̃kl i k γ̄k Np for l = 1, . . . , NSk ; the total weight of these particles is NSk Np . Subsequently, draw Np − NSk ξkl NSk particles ωkl = !NSk independently from the empirical measure π̄ki = 2 l=1 l l=1 ω̃k γ̄ki Np = 1 ; the total weight of this is 1 − Np ω̃kl δ{ξ̃l } and set k NSk Np . Else, i.e., if NSk < 12 Np , then copy the NSk particles, i.e., ξkl = ξ˜kl , and set ωkl = 21 ω̃kl /γ̄ki for l = 1, . . . , NSk ; the total weight of these particles is 12 . For the remaining weight of 12 , we NSk independently draw Np − NSk particles and set ωkl = 1 !NSk l l=1 ω̃k 2 i γ̄k (Np − NSk ) = ξkl 1 2 N p − N Sk from the empirical measure π̄ki = 2 l=1 ω̃kl δ{ξl } k . 3.3 Development of a Petri net model of free flight In order to apply the IPS algorithm toward the assessment of collision risk of free flight, we need to develop a Monte Carlo simulator of these operations such that the simulated trajectories constitute realizations of a hybrid state strong Markov process. Everdij 35 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation and Blom [40, 42] have developed a Stochastically and Dynamically Coloured Petri Net (SDCPN) formalism that ensures the specification of a free flight Monte Carlo simulation model which is of the appropriate class. This section explains how the SDCPN formalism has been used to develop a Monte Carlo simulation model of a particular free flight design. The specific free flight design for which we wish to estimate the collision risk by sequential Monte Carlo simulation is the Autonomous Mediterranean Free Flight (AMFF) operation [79]. AMFF has been developed to study the introduction of autonomous free flight operation in Mediterranean airspace. In parallel to the current study, the safety of the AMFF operation has been addressed in [77], following a fault tree analysis approach. These results show that application of AMFF seems feasible to accommodate low en-route traffic conditions over the Mediterranean. However, this study also concludes that the fault tree approach has limited analysis capabilities in showing whether AMFF can safely accommodate a higher traffic density. For this, a need was identified to use a more advanced safety risk assessment approach that considers complex situations involving dynamic interactions between multiple human and systems. The current study addresses this for AMFF under relatively high traffic densities. For the development of a Petri net model of the AMFF operation, two key challenges have to be addressed: a syntactical challenge of developing a model that is consistent, complete, and unambiguous; and a semantics challenge of representing the AMFF operation sufficiently well. This section shows how the SDCPN formalism has been used to address the syntactical challenge. Addressing the semantics challenge falls outside the scope of this study. 3.3.1 Specification of Petri net model In using the (S)DCPN formalism [40, 41, 42] in modeling more and more complex multiagent hybrid systems, it was found that the compositional specification power of Petri nets reaches its limitations. More specifically, the following problems were identified: 1. For the modeling of a complete Petri net for complex systems, a hierarchical approach is necessary in order to be able to separate local modeling issues from global or interaction modeling issues. 2. Often the addition of an interconnection between two low-level Petri nets leads to a duplication of transitions and arcs in the receiving Petri net. 3. The number of interconnections between the different low level Petri nets tends to grow quadratically with the size of the Petri net. Everdij et al. [45] explained which Petri net model specification approaches from literature solve problem 1, and developed novel approaches to solve problems 2 and 3. Together, these approaches are integrated into a compositional specification approach for SDCPN, which is explained below. In order to avoid problem 1, the compositional specification of an SDCPN for a complex process or operation starts with developing a Local Petri Net (LPN) for each agent 36 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation that exists in the process or operation (e.g., air traffic controller, pilot, navigation and surveillance equipment). Essential is that these LPNs are allowed to be connected with other Petri net parts in such a way that the number of tokens residing in an LPN is not influenced by these interconnections. We use two types of interconnections between nodes and arcs in different LPNs: • Enabling arc (or inhibitor arc) from one place in one LPN to one transition in another LPN. These types of arcs have been used widely in Petri net literature. • Interaction Petri Net (IPN) from one (or more) transition(s) in one LPN to one (or more) transition(s) in another LPN. In order to avoid problems 2 and 3, high level interconnection arcs have been introduced that allow, with well-defined meanings, arcs to initiate and/or to end on the edge of the box surrounding an LPN [45]. The meaning of these interconnections from or to an edge of a box allows several arcs or transitions to be represented by only one arc or transition. 3.3.2 High level interconnection arcs As an illustration of how high level interconnection arcs avoid duplication of arcs and transitions within an LPN and duplication of arcs between LPNs, we give three examples of these high level interconnection arcs. See [45] for a complete overview of these high level interconnection arcs. In the first example, Figure 3.1, an enabling arc starts on the edge of an LPN box and ends on a transition in another LPN box, means that enabling arcs initiate from all places in the first LPN and end on duplications of this transition in the second LPN. The duplicated transitions should have the same guard or delay function and the same firing function and their input places should have the same color type. This high level interconnection arc is not defined for inhibitor or ordinary arcs instead of enabling arcs. Figure 3.1: High level enabling arc starts at the edge of an LPN box. In the second example, Figure 3.2, an enabling arc ends on the edge of an LPN box.1 This means that for each transition in the receiving LPN a copy of this enabling arc should be in place. Figure 3.2 shows an example of this high level interconnection arc. This type of high level arc can also be used with inhibitor arcs instead of enabling arcs. 1 Figures 3.2, 3.3, 3.4, and 3.5 are from [7], with can kind permission of Springer Science and Business Media. 37 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation It cannot be used with ordinary arcs, due to the restriction that the number of tokens in an LPN should remain the same. In the third example, Figure 3.3, an ordinary arc starts on the edge of an LPN box and ends on a transition inside the same box. This means that ordinary arcs start from all places in the LPN box to duplications of this transition. The duplicated transitions should have the same guard or delay function and the same firing function and their set of input places should have the same set of color types. Figure 3.3 illustrates how this avoids both the duplication of transitions and arcs within an LPN, and the duplication of arcs between LPNs. Figure 3.2: High level enabling arc ends at the edge of an LPN box. Figure 3.3: High level ordinary arc starts on the edge of an LPN box and ends on a transition inside the same LPN box. 3.3.3 Agents and LPNs to represent AMFF operations In the Petri net modeling of AMFF operations for the purpose of an initial collision risk assessment, the following agents are taken into account: • Aircraft • Pilot-Flying (PF) • Pilot-Not-Flying (PNF) • Airborne Guidance, Navigation and Control (AGNC) • Airborne Separation Assistance System (ASAS) 38 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation • Communication, Navigation and Surveillance (CNS) It should be noticed that our initial model representing AMFF operations, does not yet incorporate other relevant agents such as Airborne Collision Avoidance System (ACAS), Airline Operations Centre (AOC), Air Traffic Control (ATC), or an environmental model. This should be taken into account when interpreting the simulation results obtained with this initial model. Per agent, particular LPNs and IPNs have been developed and subsequently the interactions between these LPNs and IPNs have been specified. A listing of LPNs per agent reads as follows: • Aircraft LPNs: – – – – Type Evolution mode Systems mode Emergency mode • Pilot-Flying (PF) LPNs: – – – – – – State situation awareness Intent situation awareness Goal memory Current goal Task performance Cognitive mode • Pilot-Not-Flying (PNF) LPNs: – Current goal – Task performance • Airborne Separation Assistance System (ASAS) LPNs: – – – – – – – – Processing Alerting Audio alerting Surveillance System mode Priority switch mode Anti-priority switch mode Predictive alerting (of other aircraft) • Airborne Guidance, Navigation and Control (AGNC) LPNs: – Indicators failure mode for PF 39 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation – – – – – – – – – – – – – – – – – Engine failure mode for PF Navigation failure indicator for PF ASAS failure indicator for PF Automatic Dependent Surveillance-Broadcast (ADS-B) receiver failure indicator for PF ADS-B transmitter failure indicator for PF Indicator failure mode for PNF Guidance mode Horizontal guidance configuration mode Vertical guidance configuration mode Flight Management System (FMS) flightplan Airborne Global Positioning System (GPS) receiver Airborne Inertial Reference System (IRS) Altimeter Horizontal position processing Vertical position processing ADS-B transmission ADS-B receiver • Communication, Navigation and Surveillance (CNS) LPNs: – Global GPS / satellites – Global ADS-B ether frequency – Secondary Surveillance Radar (SSR) mode-S frequency The actual number of LPNs in the whole model then equals 38N + 3, where N is the number of aircraft. In addition there are 35 IPNs per aircraft; hence the number of IPNs equals 35N . Brownian motion enters in each of the aircraft evolution mode LNPs. In this initial model these Brownian motions are assumed to be independent. 3.3.4 Interconnected LPNs of ASAS The approach taken in developing the AMFF concept of operation [79] is to avoid much information exchange between aircraft and to avoid dedicated decision-making by artificial intelligent machines. Although the conflict detection and resolution approach developed for AMFF has its roots in the modified potential field approach [52], it has some significant deviations from this. The main deviation is that conflict resolution in AMFF is intentionally designed not to take the potential field of all aircraft into account. The resulting AMFF design can be summarized as follows: • All aircraft are supposed to be equipped with Automatic Dependent SurveillanceBroadcast (ADS-B), which is a system that periodically broadcasts own aircraft state information, and continuously receives the state information messages broadcasted by aircraft that fly within broadcasting range (∼ 100 Nm). 40 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation • To comply with pilot preferences, conflict resolution algorithms are designed to solve multiple conflicts one by one rather than according to a full concurrent way, e.g., see [52]. • Conflict detection and resolution are state-based, i.e., intent information, such as information at which point surrounding aircraft will change course or height, is supposed to be unknown. • The vertical separation minimum is 1000 ft and the horizontal separation minimum is 5 Nm. A conflict is detected if these separation minima will be violated within 6 minutes. • The conflict resolution process consists of two phases. During the first phase, one of the aircraft crews should make a resolution maneuver. If this does not work, then during the second phase, both crews should make a resolution maneuver. • Prior to the first phase, the crew is warned when an ASAS alert is expected to occur if no preventive action would be timely implemented; this prediction is done by a system referred to as P-ASAS (Predictive ASAS). • Conflict co-ordination does not take place explicitly, i.e., there is no communication on when and how a resolution maneuver will be executed. • All aircraft are supposed to use the same resolution algorithm, and all crew are assumed to use ASAS and to collaborate in line with the procedures. • Two conflict resolution maneuver options are presented: one in vertical and one in horizontal direction. The pilot decides which option to execute. • ASAS related information is presented to the crew through a Cockpit Display of Traffic Information (CDTI). ASAS is modeled through the SDCPN depicted in Figure 3.4. The ADS-B information received from other aircraft is processed by the LPN ASAS surveillance. Together with the information about its own aircraft state information (from AGNC), the LPN ASAS processing uses this information to perform conflict detection and resolution functionalities. Subsequently, the LPN ASAS alerting and the LPN P-ASAS processing informs the PF and PNF through ASAS audio alerting about any aircraft that is in potential ASAS conflict with its own aircraft, and suggests resolution options including a prioritization. Three complementary LPNs represent non-nominal behavior modes, each combination of which has a specific influence on the ASAS alerting LPN: • ASAS system mode may be working, failed, or corrupted (failed or corrupted mode also influences the ASAS processing LPN). • ASAS priority switching mode; under emergency, the PF switches this from “off” to “on.” • ASAS anti-priority switch; this is switched from “off” to “on” when own ADS-B is not working. 41 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Figure 3.4: The agent ASAS in AMFF is modelled by eight LPNs, a number of ordinary and enabling arcs, and two IPNs (with one place each). 3.3.5 Interconnected LPNs of “Pilot Flying” This subsection illustrates the specific Petri net model developed for the Pilot Flying. For the semantical basis behind this type of model, see [12, 13, 16, 27, 90]. A graphical representation of all LPNs the Pilot-Flying consists of, is given in Figure 3.5. 42 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Figure 3.5: The agent Pilot-Flying in AMFF is modelled by six LPNs, and a number of ordinary and enabling arcs and some IPNs, consisting of one place and input and output arcs. The interconnections with other agents are not shown. The Human-Machine-Interface where sound or visual clues might indicate that attention should be paid to a particular issue, is represented by a LPN that does not belong to the Pilot-Flying as agent and is therefore not depicted in the figure. Similarly, the arcs to or from any other agent are not shown in Figure 3.5. Because of the very nature of Petri nets, these arcs can easily be added during the follow-up specification cycle. To get an understanding of the different LPNs, a good starting point might be the LPN “Current Goal” (at the bottom of the figure) as it represents the objective the Pilot-Flying is currently working on. Examples of such goals are “Collision Avoidance,” “Conflict Resolution,” and “Horizontal Navigation.” For each of these goals, the pilot executes a number of tasks in a prescribed or conditional order, represented in the LPN “Task 43 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Performance.” Examples of such tasks are “Monitoring and Decision,”“Execution,” and “Execution Monitoring.” If all relevant tasks for the current goal are considered executed, the pilot chooses another goal, thereby using his memory (where goals deserving attention might be stored, represented by the LPN “Goal Memory”) and the Human-MachineInterface. His memory where goals deserving attention might be stored is represented as the LPN “Goal Memory” in Figure 3.5. So, the LPNs “Current Goal,” “Task Performance,” and “Goal Memory” are important in the modeling of which task the Pilot-Flying is executing. The other three LPNs are important in the modeling on how the Pilot-Flying is executing the tasks. The LPN “State SA”, where SA stands for Situation Awareness, represents the relevant perception of the pilot about the states of elements in his environment, e.g., whether he is aware of an engine failure. The LPN “Intent SA” represents the intent, e.g., whether he intends to leave the free flight airspace. The LPN “Cognitive mode” represents whether the pilot is in an opportunistic mode, leading to a high but error-prone throughput, or in a tactical mode, leading to a moderate throughput with a low error probability. 3.3.6 Verification, parameterization, and validation of simulation model The compositionally specified SDCPN model enables a systematic implementation, verification and validation of the resulting Monte Carlo simulator. This is done through the following systematic steps: • Software code testing. This is done through conducting the following sequence of testing: random number generation, statistical distributions, common functions, each LPN implementation, each agent implementation, interactions between all agents, full Monte Carlo simulation. • Numerical approximation testing. This is needed to learn choosing an appropriate numerical integration step size and an appropriate number of particular Monte Carlo simulations. • Graphical user interface testing. This is to verify that the input and output of data works well. • Parameterization. This is done through a search for literature and statistical sources, and complemented by expert interviews. The fusion of these different pieces of information is accomplished following a Bayesian approach. • Initial model validation through studying Monte Carlo simulator behavior and sensitivities to parameter changes under dedicated scenarios. • Overall validation, which is directed to the evaluation of differences between model and reality and what effect these differences have at the assessed risk level. The last validation step typically is done at a later stage in the risk assessment process, with the help of active participation of operational experts [38, 44]. 44 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Table 3.1: Dimensional analysis of agent PF. Pilot-Flying (PF) Number of places Maximum color LPNs and IPNs state space Pilot Flying (PF) LPNs: State Situation Awareness 1 R3 Intent Situation Awareness 1 R5 Goal memory 1 R Current goal 7 R Task performance 7 R2 Cognitive mode 2 R Pilot Flying (PF) internal IPNs: Int-PF-GM1 1 R Int-PF-GM2 1 R Int-PF-GM3 1 R Int-PF-GM4 1 R Int-PF-GM5 1 R Int-PF-TP1 1 R2 Int-PF-TP2 1 R Int-PF-ISA 1 R Pilot Flying (PF) external IPNs: Ext-PF-Audio-PF 5 Ext-PF-PNF 1 R Ext-PF-PASAS 1 ∅ Ext-PF-SSA-1 1 ∅ Ext-PF-SSA-2 1 R Ext-PF-SSA-3 1 R Ext-PF-SSA-4 1 R Ext-PF-SSA-5 1 R PRODUCT 490 R28 3.3.7 Dimensions of Monte Carlo simulation model Now, we analyze the dimensions of the joint state space of the resulting Monte Carlo simulation model. In Table 3.1 and Table 3.2, this is done for the agents PF and ASAS respectively, including all LPNs and all IPNs that end on one of these LPNs. The second column gives the number of places in the LPN or IPN. The third column gives the maximum state space of the color used within an LPN or IPN. We also perform this analysis to the LPNs and IPNs of the other agents. The resulting number of product places and product state spaces is given in Table 3.3. This table brings into account that of each agent, except global CNS, there is one per aircraft. 45 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Table 3.2: Dimensional analysis of agent ASAS. ASAS Number of places Maximum color LPNs and IPNs state space ASAS LPNs: Processing 1 R13+12N Alerting 2 R7 Audio alerting 2 ∅ 11+9N Surveillance 1 R System mode 3 ∅ Priority switch mode 2 ∅ Anti-priority switch mode 2 ∅ Predictive alerting 1 R3 ASAS internal IPNs: Int-ASAS-Resolution 1 ∅ Int-ASAS-Audio 1 ∅ ASAS external IPNs: Ext-ASAS-PF 1 R3 Ext-PASAS-PNF 1 ∅ Ext-ASASProc-PNF 1 ∅ Ext-ASASurv-ADSB-Global 1 R Ext-ASASprio-PNF PRODUCT 1 ∅ 38+21N PRODUCT 48 R The product places of the global CNS agent form the θt0 state space M 0 . The correspond0 ing state space is empty, which means places of the "N thati there is no xt . The product other agents form the state space i=1 M of the proces components θti , i = 1, . . . , N . Per aircraft, the number of product places is |M i | ≈ 0.777 × 1012 . The colors attached to the places in the other agents form the process components xit , i = 1, . . . , N , each of which assumes values in R126+21N . Each of the scenarios considered in the next subsection has eight aircraft, so N = 8. This means that the number of product places equals ≈ 16 × (0.777 × 1012 )8 ≈ 2.13 × 1096 , and that the product of the color state spaces equals R2352 . 3.4 Simulated scenarios and collision risk estimates The IPS algorithm of Section 3.2 is now applied to three hypothetical AMFF air traffic scenarios. The first scenario has eight aircraft that fly at the same flight level and their flight plans cause them to fly through the same point in airspace at the same moment in time. The second scenario has one aircraft flying through a virtual infinite airspace 46 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Table 3.3: Dimensional analysis of complete SDCPN. Agent Number of product Maximum color places product state space Aircraft 24N R13N Pilot Flying (PF) 490N R28N N Pilot-not-Flying (PNF) 7 R3N 16 N AGNC (15 × 2 ) R45N 2 ASAS 48N R37N +21N Global CNS 16 ∅ 12 N 126N +21N 2 PRODUCT ≈ 16 × (0.777 × 10 ) R of randomly distributed aircraft, with a density 3 times as high as in a current high capacity en route area. The third scenario is the same as the second, except that the aircraft density is four times lower. Prior to describing these scenarios and simulation results, we explain the parametrization of the IPS algorithm used. Table 3.4: Parameter values of free flight enabling technical systems. Model Parameter Probability Global GPS down 1.0 × 10−5 Global ADS-B down 1.0 × 10−6 Aircraft ADS-B receiver down 5.0 × 10−5 Aircraft ADS-B transmitter down 5.0 × 10−5 Aircraft ASAS system mode corrupted 5.0 × 10−5 Aircraft ASAS system mode failure 5.0 × 10−5 3.4.1 Parameterization of the IPS simulations The main safety critical parameter settings of the free flight enabling technical systems (GPS, ADS-B and ASAS) are given in Table 3.4; in the table, global ADS-B down refers to frequency congestion/overload of the data transfer technology used for ADS-B. The IPS conflict levels k are defined by parameter values for lateral conflict distance dk , conflict height hk , and time to conflict ∆k . These values have been determined through two steps. The first was to let an operational expert make a best guess of proper parameter values. Next, during initial simulations with the IPS some fine tuning of the number of levels and of parameter values per level has been done. The resulting values are given in Table 3.5. 47 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Table 3.5: k 1 dk (Nm) 4.5 hk (ft) 900 ∆k (min) 8 IPS conflict level parameter values. 2 3 4 5 6 7 8 4.5 4.5 4.5 2.5 1.25 0.5 0.054 900 900 900 900 500 250 131 2.5 1.5 0 0 0 0 0 3.4.2 Eight aircraft on collision course In this simulation eight aircraft start at the same flight level, some 135 Nm (250 km) out of each other, and fly in eight 45 degrees differing directions with a ground speed of 466 knots (= 240 m/s), all aiming to pass through the same point in airspace. By running ten times the IPS algorithm the collision risk is estimated ten times. The number of particles per IPS simulation run is 12,000. The total simulation time took about 20 hours on two machines, and the load of computer memory per machine was about 2.0 gigabyte. For the first eight IPS runs, the estimated fractions γ̄ki are given in Table 3.6 for each of the conflict levels, k = 1, . . . , 8, and aircraft i = 1. It can be seen that the first and sixth IPS runs have zero particles that reach the last (8th level). Hence the first and sixth IPS runs yield γ̄8i = 0. This is a clear example of particle depletion. The IPS estimated mean probability for one aircraft to collide with any of the other seven aircraft equals 2.2 × 10−5 . The minimum and maximum values now are respectively a factor 250 lower and a factor 4 higher than the mean value. We also verified that this risk value was not sensitive at all to the failure rates of the ASAS related technical systems. In [52] a similar eight aircraft encounter scenario has been simulated some thousand times, without experiencing any collision event. At a collision probability value of 2.2 × 10−5 , the chance to count at least one collision would be in the order of 5 %. As such the current results agree quite well with the fact that in these earlier simulations for the eight aircraft scenario no collision has been observed. We also verified that the novel simulation results for the eight aircraft scenario agreed quite well with the expectation of the designers of the AMFF operational concept. 3.4.3 Free flight through an artificially constructed airspace In this simulation the complete airspace is divided into packed containers. Within each container a fixed number of seven aircraft (i = 2, . . . , 8) fly at arbitrary position and in arbitrary direction at a ground speed of 466 Nm/hr. One additional aircraft (i = 1) aims to fly straight through a sequence of connected containers, at the same speed, and the aim is to estimate its probability of collision with any of the other aircraft per unit time of flying. Per container, the aircraft within it behave the same. This means that we have to simulate each aircraft in one container only, as long as we apply the ASAS conflict prediction 48 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Table 3.6: Fractions counted during eight IPS runs of scenario 1. Level 1 2 3 4 5 6 7 8 Product of fractions 1st IPS 1.000 0.528 0.426 0.033 0.175 0.267 0.150 0.000 0.0 2nd IPS 1.000 0.529 0.429 0.036 0.180 0.158 0.268 0.009 5.58× 10−7 3rd IPS 1.000 0.539 0.424 0.035 0.183 0.177 0.281 0.233 1.67× 10−5 4th IPS 1.000 0.533 0.431 0.037 0.181 0.144 0.427 0.043 4.01× 10−6 5th IPS 1.000 0.537 0.421 0.039 0.142 0.255 0.645 0.455 9.33× 10−5 6th IPS 1.000 0.538 0.428 0.031 0.157 0.138 0.208 0.000 0.0 7th IPS 1.000 0.536 0.426 0.044 0.181 0.295 0.253 0.006 8.00× 10−7 8th IPS 1.000 0.539 0.418 0.039 0.147 0.146 0.295 0.815 4.48× 10−5 and resolution also to aircraft copies in the neighboring containers. In principle this can mean that an aircraft experiences a conflict with its own copy in a neighboring container. This also means that the size of a container should not go below a certain minimum size. By changing container size we can vary traffic density. To choose the appropriate traffic density, our reference point is the highest number (17) of aircraft counted at 23rd July 1999 in an en-route area near Frankfurt of size 1 degree × 1 degree × FL290-FL420. This comes down to 0.0032 a/c per Nm3 . For our simulation we assume a 3 times higher traffic density, i.e., 0.01 a/c per Nm3 . This resulted in choosing containers having a length of 40 Nm, a width of 40 Nm, and a height of 3000 feet and with 8 aircraft flying in such a container. By running the IPS algorithm ten times (+ one extra later on) over 25 minutes (5 minutes to allow convergence and 20 minutes to estimate collision probability) the collision probability per unit time of flying has been estimated. The number of particles per IPS simulation run is 10,000. The total simulation time took about 300 hours on two machines, and the load of computer memory per machine was about 2.0 gigabyte. For the first eight IPS runs, the estimated fractions γ̄ki are given in Table 3.7 for each of the conflict levels, k = 1, . . . , 8, for aircraft i = 1. The estimated mean probability of collisions per 20 minutes of aircraft flight equals 5.22 × 10−5 , which is equal to a probability of collisions per aircraft flight hour of 1.6 × 10−4 , with minimum and maximum values respectively a factor four lower and higher. We also verified that this risk value was not sensitive at all to the failure rates of the ASAS related technical systems. One should be aware that this value has been estimated for the simulation model of the intended AMFF operation. Hence the question is what this means for the intended AMFF operation itself? By definition a simulation model of AMFF differs from the intended AMFF operation. If it can be shown that the combined effect of these differences on the risk level is small, then the results obtained for the simulation model may be 49 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation Table 3.7: Fractions counted during eight IPS runs of scenario 2. Level 1 2 3 4 5 6 7 8 Product of fractions 1st IPS 0.922 0.567 0.665 0.319 0.370 0.181 0.130 0.067 6.42× 10−5 2nd IPS 0.917 0.551 0.666 0.331 0.367 0.158 0.209 0.005 6.76× 10−6 3rd IPS 0.929 0.560 0.674 0.323 0.371 0.162 0.174 0.094 1.11× 10−4 4th IPS 0.926 0.559 0.676 0.321 0.379 0.171 0.145 0.066 6.99× 10−5 5th IPS 0.925 0.554 0.672 0.328 0.363 0.164 0.162 0.002 2.57× 10−6 6th IPS 0.925 0.551 0.673 0.321 0.345 0.181 0.170 0.150 1.75× 10−4 7th IPS 0.925 0.561 0.664 0.334 0.366 0.148 0.214 0.015 1.99× 10−5 8th IPS 0.921 0.556 0.670 0.331 0.343 0.191 0.215 0.019 2.98× 10−5 considered as a good representation of the accident risk of the intended operation. In order to assess the combined effect of these differences there is need to perform a bias and uncertainty assessment [38]. In order to better learn understanding of what causes the collision risk of the simulation model to be relatively high, we performed an extra IPS run, and memorized in static memory for each particle the ancestor history at each of the eight levels. This allowed us to trace back what happened for the particles that hit the last level set (i.e., collision). There appeared to be five different collision events. Evaluation of these five collision events showed that all five happened under nominal safety critical conditions. Four of the five collisions were due to a growing number of multiple conflicts that could not be solved in time under the operational concept adopted. The fifth collision was of another type: at quite a late moment finally a conflict between two aircraft was solved with a maneuver by one of the two aircraft. However because of this maneuver there was a sudden collision with a third nearby aircraft. These detailed evaluations of the five collision events of the 11th IPS run also showed that a significant increase of collision risk is caused by the relatively small height (4000 ft) of a container. Because of this small height it happened that an aircraft in one container came in conflict with a copy of its own in a neighboring container, and in such a situation there was an undesired limitation in conflict resolution options, and thus an undesired artificial increase in collision risk. The results in this section seem to indicate that the key factor in the increased risk of collision for encounters with homogeneous traffic in the background — as opposed to the eight encountering aircraft only scenario — are the multiple conflicts. Under the far higher traffic densities than what the AMFF operational concept was designed for, it is not always possible to timely solve a sufficiently high fraction of those multiple conflicts. On the basis of this finding one would expect that the collision risk would decrease faster than linear with a decrease in traffic density. The validity of this expectation is verified by the next scenario. 50 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation 3.4.4 Reduction of the aircraft density by a factor four Now we enlarge the length and width of each container by a factor two. This means that the traffic density has gone down by a factor four. Hence the density is now 43 of the density counted on 23rd July 1999 in the en-route area near Frankfurt. This still is a factor 2.5 higher than current average density above Europe. At the same time simulated flying time has been increased to 60 minutes (with 10 minutes prior flying to guarantee convergence). By running four times the IPS algorithm the collision risk is estimated four times. The number of particles per IPS simulation run is 10,000. The total simulation time took about 280 hours on two machines, and the load of computer memory per machine was about 2.0 gigabyte. For these IPS runs, the estimated fractions γ̄ki are given in Table 3.8 for each of the conflict levels, k = 1, . . . , 8, for aircraft i = 1. Table 3.8: Fractions counted during four IPS runs of scenario 3. Level 1 2 3 4 5 6 7 8 Product of fractions 1st IPS 0.755 0.295 0.476 0.263 0.321 0.068 0.156 0.011 1.07 × 10−6 2nd IPS 0.750 0.292 0.475 0.258 0.315 0.088 0.367 0.059 1.61 × 10−5 3rd IPS 0.752 0.286 0.497 0.266 0.300 0.082 0.290 0.021 4.31 × 10−6 4th IPS 0.749 0.285 0.487 0.267 0.328 0.096 0.254 0.005 1.07 × 10−6 The estimated mean probability of collision per aircraft flight hour equals 5.64 × 10−6 , with minimum and maximum values respectively a factor five lower and higher. This is about a factor 30 lower than the previous scenario with a four times higher aircraft density. Thus, for the model there is a steep decrease of collision probability with decrease of traffic density, and this agrees well with the expectation at the end of the previous section. 3.4.5 Discussion of IPS simulation results Because of the IPS simulation approach we were able to estimate collision risk for complex multiple aircraft scenarios. The large increase in handling complexity of multiple aircraft encounter situations is a major improvement over what was feasible before for two aircraft flying in a parallel route structure [10, 43]. Inherent to the IPS way of simulation, the dynamic memory of the computers used appeared to pose the main limitation on the full exploitation of the novel sequential Monte Carlo simulation approach. This also prevented performing a bias and uncertainty assessment for the differences between 51 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation the simulation model and the AMFF operation. As long as such a bias and uncertainty assessment has not been performed, any conclusion drawn from the simulation apply to the simulation model only, and need not apply to the intended AMFF operation. The simulations performed for a model of AMFF allow free flight operational concept developers to learn characteristics of the simulation model. Because of the sequential Monte Carlo simulation based speed up, these simulations can show events that have not been observed before in Monte Carlo simulations of an AMFF model. Under far higher traffic densities than what the AMFF operational concept has been designed for, the simulations of the model shows it is not always possible to timely solve multiple conflicts. As a result of this, at high traffic levels there is a significant chance that multiple conflicts are clogging together, and this eventually may cause a non-negligible chance of collision between aircraft in the simulation model. It has also been shown that by lowering traffic density, the chance of collision for the model rapidly goes down. 3.5 Concluding remarks We studied collision risk estimation of a free flight operation through a sequential Monte Carlo simulation. First a Monte Carlo simulation model of this free flight operational concept has been specified in a compositional way using the Stochastically and Dynamically Coloured Petri Net (SDCPN). Subsequently a novel sequential Monte Carlo simulation method [23, 25] has been extended for application to collision risk estimation in air traffic, and has subsequently been applied to an SDCPN model of free flight. The results obtained show that the novel simulation model specification and collision risk estimation method allow to speed up the Monte Carlo simulations for much more complex air traffic encounter situations than what was possible before, e.g., [10, 43]. Moreover, for the simulation model of the free flight operational concept considered, behavior has been made visible that was expected by free flight concept designers, but could not be observed in straightforward Monte Carlo simulations of free flight concepts (e.g., [52]): the rare chance of clogging multiple conflicts at far higher traffic density levels than where the particular concept has been designed for. Hence, further attention has to be drawn toward the development and incorporation in the particular operational concept design of advanced methods in handling multiple conflicts. For example, Hoekstra [52] studied a conflict resolution approach that performs better than the one adopted in the AMFF concept. In addition, there are some complementary developments that aim to develop complex conflict resolution solvers with some guaranteed level of performance [33, 76] under nominal conditions, and ways to incorporate situation awareness views by human operators (pilots and/or controllers) in these combinatorial conflict resolution problems [30]. The initial collision risk estimation results obtained with our sequential Monte Carlo simulation of free flight provides valuable feedback to the design team and allows them to learn from Monte Carlo simulation results they have never seen before. This allows 52 3 Free Flight Collision Risk Estimation by Sequential Monte Carlo Simulation them to significantly improve their understanding of when and why multiple conflicts are not solved in time anymore in the simulation model. Subsequently the operational concept designers can use their better understanding for adapting the free flight design such that it can better bring into account future high traffic levels. In its current form the sequential Monte Carlo simulation approach works well, but at the same time poses very high requirements on the availability of dynamic computer memory and simulation time. The good message is that in literature on sequential Monte Carlo simulation, e.g., see [31, 32, 35, 50, 67, 74], complementary directions have been developed which remain to be explored for application to free flight collision risk estimation. These potential improvements of sequential Monte Carlo simulation will be studied in the work package WP7. Acknowledgement. The authors thank Mariken Everdij (NLR) for valuable discussions and a thorough review of a draft version of this chapter. 53 A Mathematical Background A.1 Stationary Point Processes, Palm measure and applications to collision risk A.1.1 Point process In many applications, we observe some discrete events occurring at times T0 , T1 ,· · · , Tn ; more formally these observations can be encompassed in the sequence (Tn , n ∈ Z), where the Tn ≤ Tn+1 are random variables. We call this process a point process on R+ , and we say that the point process is simple if Tn < Tn+1 for each n. Another point of view consists in counting the!number N ((a, b]) of events observed during the time interval (a, b]; then N ((a, b]) = n∈Z 1(a,b] (Tn ) and Tn = inf{t : N (−∞, t]) = n}. We say that the point process is stationary if the joint probability of the number of events in m disjoint intervals I1 , · · · , Im is invariant by translation, ie for all m ∈ N and all t ∈ R P(N (I1 ) = k1 , · · · , N (Im ) = km ) = P(N (I1 + t) = k1 , · · · , N (Im + t) = km ). Let introduce for each t, the shift mapping θt defined on the probability space by N (θt ω, I) = N (ω, I + t). Therefore, the point process should be stationary if and only if P ◦ θt = P. These mappings satisfy the following properties: θt ◦ θs = θt+s , θt−1 = θ−t and for each Tn ≤ t < Tn+1 , θt Tn = T0 and θTn+1 = T1 with the convention that T0 ≤ 0 (see Figure (A.1)). 0 T0 (ω) T−1(θt ω) T0 (θt ω) T1(ω) t T2 (ω) T1(θt ω) 0 Figure A.1: The shift. Note the numbering of the points Ti In the following, we will consider only stationary point processes. The positive number λ = E(N ((0, 1])) is called the intensity of the point process; this number can be infinite, so we assume now that λ > 0 and is finite. Since the stationarity, for all interval I, the measure λ(I) := E(N (I)) is invariant by translation, so it is proportional to the Lebesgue measure l on R and we deduce that λ(I) = λl(I), with λ > 0. 54 A Mathematical Background A.1.2 Palm measure An important area in which point processes are used is queueing theory. For example, the number of clients at a bank counter is a point process and the study of the counting process N is useful to decide to open another counter or not. Nevertheless, the point of view of the client is different, since the more pertinent for him is the waiting time at the instant he arrives in the queue or more generally all information related to the counting process given its arrival time. This suggests us to introduce a new probability PN0 , called Palm measure, defined only on the event T0 = 0, i.e. PN0 (T0 = 0) = 1. Due to the stationarity of the point process, the random variables Sn = Tn − T0 are identically distributed, therefore we could define the empirical function F0 (t) of T1 = S0 and we will have F0 (t) = PN0 (Sn ≤ t) for all n. Thus, the Palm measure need to be a measure PN0 such that PN0 (A) = 0 for all measurable sets A such A ∩ {T0 = 0} = ∅. For that, we need to count the number of events in A, not by observing all the process but ! by shifting successively all the Tn at the origin of the time, and so we consider the sum n∈Z 1A ◦ θTn . Nevertheless, the Palm measure being a probability we have to normalize the sum, from which the probability 3, 4 1 1 0 PN (A) = E [(1A ◦ θTn )1C (Tn )] := E (1A ◦ θs )N (ds) . λl(C) λl(C) C It can be proved that this probability is independent of C and is supported by the event {T0 = 0}, so we adopt the following definition PN0 (A) = ( 1 ' E (1A ◦ θTn )1(0,1] (Tn ) . λ Knowing the Palm probability, it is also possible to obtain the law of the Point process by the following relation , ∞ P(A) = λ PN0 (T1 > t, θt ∈ A)dt. (A.1.1) 0 A.1.3 Residual waiting time and spent waiting time Let us come back to the empirical function F0 (t) of the Sn under the probability PN0 , i.e. F0 (t) = PN0 (Sn ≤ t) and introduce the residual waiting time at time t ≥ 0 defined by W (t) = TN (t)+1 − t and the spent waiting time at time t A(t) = t − TN (t) . 55 A Mathematical Background For a stationary point process, you can take t = 0, in that case W (0) = T1 and A(0) = −T0 . Using (A.1.1) we get for A = {T1 > v, −T0 > u} with u, v ≥ 0, , ∞ (1 − F0 (s))ds, P(T1 > v, −T0 > u) = λ v+u from which we derive P(T1 > v) = λ , ∞ v (1 − F0 (s))ds, P(−T0 > u) = λ , ∞ u (1 − F0 (s))ds. We conclude that −T0 and T1 are identically distributed for the probability P . Moreover, the law of S0 under P is given by , P(−T0 + T1 ≤ x) = λyF0 (dy), [0,x] which is in general different from F0 (x). For instance, 5 6 , ∞ Var0N (S0 ) 2 0 λy F0 (dy) = EN (S0 ) 1 + E(S0 ) = 0 (S ))2 , (EN 0 0 0 (S ) = 1/λ, so E(S ) = E 0 (S ) iff under the Palm probability, the variance of since EN 0 0 0 N S0 is zero, therefore iff S0 is a constant and thus also all the Sn . A.1.4 Application to the collision risk We consider two aircraft A1 and A2 , each flying on straight-line paths at constant velocity v1 and v2 respectively. If the aircraft A1 is considered fixed, the aircraft A2 can be considered having a velocity vr = v2 − v1 relative to the aircraft A2 and a relative position r(t) at time t given by, r(t) = r(0) + vr t. A conflict is declared when the predicted position of the two aircraft are such that both horizontal and vertical separation parameters are infringed. Therefore, we associate a conflict volume to the aircraft A1 which is an airspace formed around aircraft A1 into which, if aircraft A2 enters, a conflict is signalled. When this volume is a sphere a radius R, a conflict occurs when the predictive distance at the Closest Point of Approach (CPA) d is smaller than or equal to R. The CPA is determined by the condition dt (r(t) · r(t)) = 0 or r(t) · vr = 0, that means the relative distance r is orthogonal to the relative velocity vr . The time tCPA of the CPA and the distance dCPA at CPA are given respectively by tCPA = − r(0) · vr , vr · vr d2CPA = r(0) · r(0) − 56 (r(0) · vr )2 . vr · vr A Mathematical Background Let θr ∈ [0, 2π[ be the angle between vr and r(0), i.e. r(0) · vr = r(0)vr cos θr , with r(0) = 1r(0)1 and vr = 1vr 1. If θr = 0 then the aircraft A2 is passing from the aircraft A1 an this situation is not relevant, so by now we assume that θr %= 0. Introducing the distance a = R/(| sin θr |), a conflict will occur if r(0) ≤ a. We consider now an aircraft A2 whose the straight-line path crosses a flow of aircraft each of which are flying on straight-line path too. Moreover, we assume that aircraft A1 of the flow are distributed according to a stationary point process with intensity λ. Let choose the time zero as the moment the aircraft A2 crosses the flow and let T0 and T1 the time separations between the two closest aircraft A1 of the flow to the aircraft A2 . That means, at time t = 0, an aircraft A1 will arrive at the crossing point at time −T0 and the other crossed this same point T1 time units ago. If F0 is the distributed function of T1 under the Palm measure of the point process, the probability Pc to have a conflict between the aircraft A2 and one of aircraft A1 is given by , ∞ Pc = 1 − P(−T0 > a/v1 , T1 > a/v1 ) = 1 − λ (1 − F0 (u))du =λ , 0 since -∞ 0 2a/v1 2a/v1 (1 − F0 (u))du, (1 − F0 (u))du = E0N (T1 ) = 1/λ. The quantity 2a/v1 being very small, we obtain the first order approximation Pc ≈ (2λa)/v1 and if in addition, F0 has a density f0 , we obtain a second order approximation 5 6 2λa a Pc ≈ 1 − f0 (0) . v1 v1 Poisson Process The Poisson process corresponds to the case F0 (t) = 1 − e−λt ; so 5 6 λa 2λa −2λa/v1 1− . Pc = 1 − e ≈ v1 v1 As vr sin θr = v2 sin θ with θ the crossing angle between the two trajectories, we obtain the well-known formula 2λa . Pc ≈ v1 Gamma Law Now, we assume that F0 is a Gamma Law with parameters p > 0 and θ > 0 whose the density is given by θp −θu p−1 f0 (u) = e u u ≥ 0. Γ(p) 57 A Mathematical Background The expectation of this law being p/θ, we have to set λ = θ/p, therefore Pc = = = = = , 2a/v1 , t θp −θu p−1 e u du 0 Γ(p) 0 , 2a/v1 p , 2a/v1 θ λ2a −λ e−θu up−1 du dt v1 Γ(p) u 0 3 5 6 4 , 1 2aθ λ2aθ p 1 λ2a p−1 − v1 v 1− (1 − v)v e dv v1 v1 Γ(p) 0 + ) 5 6 6n 5 ∞ 1 λ2a λ2aθ p 1 2 n λ2a B(p + n, 2) 1− (−1) v1 v1 Γ(p) v1 n! n=0 + ) 6n+p 5 ∞ 2 λ2a 1 λ2aθ . 1− (−1)n v1 v1 n!(p + n + 1)(p + n)Γ(p) λ2a −λ v1 dt n=0 where B(x, y) = is given by -1 0 (1 − v)y−1 v x−1 dv = Γ(x)Γ(y) Γ(x+y) . Thus the second order approximation 3 5 6 4 2aθ 2aθ p 1 Pc ≈ . 1− pv1 pv1 Γ(p + 2) Regular events with random translatories Here, we consider events which are regularly spaced out but with a possibly random translatory. So, the i-th event occurs at time i = ·, −1, 0, 1, · · · , ti = a0 + ic + bi , where c is spacing between the events in absence of random perturbations and the bi are independent and identically distributed random variable whose the empirical function, under the Palm measure, will be denoted by FB . Then the i-th event is planned at time a0 + ic, but its real time is moved by bi . This type of process has been studied by T. Lewis and Govier (1960,1963) about the arrivals of tankers at a terminal. We choose the event with indix i = 0 as initial time ; this event being planned at time −b0 with regard to the origin of time, we deduce that, given b0 , the i-th event occurs at time ic + bi − b0 , so P0N (ti ≤ t|b0 ) = FB (t + b0 − ic) − FB (b0 − ic). It is important to note that the times ti are not the times Ti of the point process induced by the events considered, since the j-th event can occur before the i-th event, even if i < j. 58 A Mathematical Background To obtain F0 (t), it is enough to observe that PN0 (T1 > t) = PN0 (N (0, t] = 0), that means no event, except the 0-th event, has occurred during the time interval (0, t]. Thus P0N (T1 > t|b0 ) = ∞ * ′ [1 − (FB (t + b0 − ic) − FB (b0 − ic))] , i=−∞ F0 (t) = 1 − where *′ ∞ i=−∞ , ∞ * ′ ∞ −∞ i=−∞ [1 − (FB (t + b0 − ic) − FB (b0 − ic))] FB (db0 ), denotes the product over all the indices i except i = 0. Let assume now, that FB has a density fB with a finite value in zero. Then, for t small enough we get the following approximation F0 (t) ≈t=0 1 − ≈t=0 1 − =t =t , ∞ −∞ i=−∞ , ∞ , 2 ∞ −∞ ∞ i=−∞ −∞ ∞ 2 i=−∞ ∞ * ′ ) [1 − tfB (b0 − ic)] fB (b0 )db0 1−t ∞ 2 ′ i=−∞ + fB (b0 − i) fB (b0 )db0 fB (b0 − i)fB (b0 )db0 − t , ∞ −∞ fB2 (b0 )db0 fZ (−ic) − tfZ (0), where fZ is the density of a random variable Z with same law that the difference of two random variables i.i.d. with density fB . We deduce that 4 3 , ∞ −1 fZ (u)du − fZ (0) F0 (t) ≈ t c −∞ 6 5 1 − fZ (0) . =t c Therefore, the probability Pc can be approximated by Pc ≈ λ , 0 2a/v1 64 3 5 ' ( 2λa a 1 −1 1 − u(c − fZ (0)) du = − fZ (0) . 1− v1 v1 c For instance, if the B is distributed as a centered Gaussian variable with variance σb2 , then 5 6 x2 1 exp − 2 , fZ (x) = 4σb 2σb π 1/2 59 A Mathematical Background and Pc ≈ 64 3 5 2λa a 1 1 . 1− − v1 v1 c 2σb π 1/2 It remains to estimate the parameter λ, nevertheless, the process being stationary the mean number of events, during a period of t unit of time, is t/c, so λ is approximatively equal to 1/c. A.2 Stochastic approach of the collision risk In the previous section we have taken into account only random perturbations along the trajectory; now we generalize our approach by considering three-dimensional perturbations. A.2.1 Introduction In safety requirement, risk assessment and conflict prediction studies, the need to model and estimate the deviation of an aircraft from its assigned trajectory, leads to implement a probabilistic approach often based on Gaussian deviations from the predicted aircraft trajectories. We would like to mention two approaches: the geometric approach developed for conflict prediction [80] and the National Aerospace Laboratory (NLR) approach developed for collision risk evaluation [4]. In [80], the authors try to estimate the conflict probability from a modelling of lateral, longitudinal and vertical deviations. They approximate the prediction errors by a Gaussian distribution and use a geometric approach, based on uncertainty ellipsoid. However, this model does not take into account the time evolution of the uncertainty in the lateral and vertical direction while it assumes in the longitudinal direction a linear increasing standard deviation with time (i.e. 0.25NM/min in cruise situations). The other approach is the NLR approach which is used in their simulation tool TOPAZ [6] but also in probabilistic conflict probe applications [4]. The aim there is to evaluate collision risk. Therefore the deviations are modelled by a certain type of stochastic differential equations (SDE’s), the so-called Kramers’ equation in Physics. Initially linear SDE’s were used for modeling Gaussian deviations [4]. With the further development of TOPAZ the extension to non-linear SDE’s was made and implemented in the TOPAZ tool. These non-linear SDE’s result in non-Gaussian deviation distributions (like the Double Exponential stationary distribution, the Double-Double Exponential stationary distribution and the Gaussian Double Exponential stationary distribution, etc). These types of distributions are often used in fitting observed data for modelling the tails of data. More precisely, the two common stationary densities adopted for modeling errors (see [20, 21, 56, 86]) are a stationary Gaussian Double Exponential density 7 1 7 1 % d 2 &8 % |d| &8 +α exp − , f (d) = (1 − α) √ exp − 2 2σ 2λ λ σ 2π 60 A Mathematical Background or a stationary Double-Double Exponential density 7 1 71 % |d| &8 % |d| &8 exp − exp − +α , f (d) = (1 − α) 2κ κ 2λ λ where the parameter α determines the relative importance of each of the functions, α is usually small, so in case of the Gaussian-Double Exponential density, the Gaussian density models the core of the data, while the Double Exponential density models the tails of the data. The context of the NLR’s approach being the collision risk evaluation, all the 4 mentioned densities can be taken into account depending on which scenario is evaluated within TOPAZ. Hence, it is possible to use the NLR’s approach to obtain a Gaussian (stationary) density or to obtain a Gaussian Double Exponential stationary density. One of the advantages of this approach is its time evolution formulation, by using stochastic differential equations. A SDE is a differential equation with a random term and hence whose the solution is, in some sense, a random function. Thus, each solution of a SDE represents a different random trajectory. When the relative position xt between two aircraft is modelled by a stochastic process, or more exactly a diffusion, the probability of conflict during a time interval [0, T ] is given by Pc (T ) = P(∃t ∈ [0, T ] : xt ∈ D), where D is the volume of conflict. This probability is very difficult to evaluate in general; but let us introduce the first hitting time of D, defined by TD = inf{t > 0 : xt ∈ D}, then Pc (T ) = P(TD ≤ T ). Thus, it is enough to known the law of TD ; of course this law is also difficult to estimate, but it is easier to study, mainly by the so-called stochastic potential theory which allows to obtain lower and upper bounds of Pc in terms of a particular measure of D. In [5], the authors introduces two another approaches : the absorbing or the transient boundary situation. In the absorbing boundary approach, instead of considering the process xt , the authors study the absorbed process x̃t defined by x̃t = xt for all t ≤ TD and x̃t = xTD for all t ≥ TD . Thus Px (TD ≤ t) = Px (x̃t ∈ D), where x = x0 . We assume that the Kolmogorov backward equation is satisfied by Px (x̃t ∈ D), that means ∂Px (x̃t ∈ D) = LPx (x̃t ∈ D), ∂t where L is a second order differential operator ; the two boundary conditions being lim Px (x̃t ∈ D) = 1 x→∂D for t > 0, x ∈ Dc , lim Px (x̃t ∈ D) = 0 t↓0 for x ∈ Dc . If the initial position of xt is not deterministic, but distributed according to the initial law π0 (x), then the probability of collision in time period (t1 , t2 ] is given by , t2 , ∂Px (x̃t ∈ D) P(x̃t1 ∈ Dc , x̃t2 ∈ D) = π0 (x)dxdt. ∂t c t1 D The coefficient µ(t) = , Dc ∂Px (x̃t ∈ D) π0 (x)dx, ∂t 61 A Mathematical Background is called in [5] the conditional collision rate at time t. Let us notice that under good hypothesis, we get , Px (x̃t ∈ D) − Px (x̃t−∆ ∈ D) π0 (x)dx µ(t) = lim ∆↓0 ∆ c D , Px (x̃t ∈ D, xt−∆ ∈ Dc ) lim = π0 (x)dx, ∆ Dc ∆↓0 P(x̃t ∈ D, x̃t−∆ ∈ Dc ) = lim . ∆↓0 ∆ In the transient boundary situation, the authors introduce the in-crossing rate ϕ(t) at time t defined by P(xt ∈ D, xt−∆ ∈ Dc ) ϕ(t) = lim . ∆↓0 ∆ We notice that the conditional collision rate at time t of the process xt is the in-crossing rate at time t of the absorbed process x̃t . In fact, the existence de ϕ(t) needs the differentiability of xt , that means we must consider the diffusion (xt , vt ) with vt = dxt /dt. Now, let x ∈ ∂D, we introduce the in-crossing rate ϕ(x, t) at x, given by ϕ(t, x)|dB|dt = P( the process xt enters in D during [t, t + dt] through dB), where dB is some domain in ∂D whose area is |dB|. If ∂D is smooth, we obtain that , ϕ(t) = ϕ(t, x)dx. ∂D Under some assumptions [5] about the stochastic process (xt , vt ), the Rice formula [75] can be applied and gives ' ( ϕ(t, x) = E 4nx , vt 5− |xt = x ftX (x), where nx is the outwards normal in x at ∂D, ftX is the density of the xt and 4nx , vt 5− = 0 if 4nx , vt 5 ≥ 0 and 4nx , vt 5− = 4nx , vt 5 otherwise. In [5], an explicit derivation of ϕ(t) is obtained for a domain D = [−s1 , s1 ] × [−s2 , s2 ] × [−s3 , s3 ] in terms of ftX (x). A.2.2 Errors Modelling In this section, we describe a model of perturbation introduced by Bakker and Blom. Let consider an aircraft following a straight-line path and denote by x(t) and v(t) its position and velocity; the couple (x(t), v(t)) is a 6-components vector called the state of the aircraft. For simplicity, we use a referential oriented according its trajectory: so the current state of the aircraft is (x// , v// , x⊢ , v⊢ , x⊥ , v⊥ ), where (x// (t), v// (t)) is the longitudinal component, (x⊢ (t), v⊢ (t)) the lateral one and (x⊥ (t), v⊥ (t)) the vertical one. Denote by χ̄ the reference course and by v̄G the reference ground speed. The deviations 62 A Mathematical Background ~ V ~ Vg + Vg T ~ V// V ~ χ T T of the aircraft from its trajectory is modelled by random perturbations of the course and the ground speed. Nevertheless, modelling the deviations of the course seems difficult; but we assume that we know the stationary distribution of the lateral deviations from the assigned trajectory. This distribution is a mixture between a Gaussian distribution and a Double Exponential one. So, instead of modelling the deviations of the course, it seems easier to model the lateral deviations, since by using a certain SDE, called Kramers’ equation, we can obtain any stationary distribution for the lateral deviation (for instance a stationary Double Exponential distribution or a Gaussian distribution for the lateral deviation). Vg χ V// // Figure A.2: Description of deviations from the trajectory in the plane Let us introduce now the state vector (χ̃, ṽG , x̃⊢ , ṽ⊢ , x̃⊥ , ṽ⊥ ) that describes the deviations from the trajectory. The lateral deviations (x̃⊢ , ṽ⊢ ) are taken relatively to the assigned trajectory (see Figure A.2). We have the following identities:   v// = cos(χ̄ + χ̃)(v̄G + ṽG )     v⊢ = sin(χ̄ + χ̃)(v̄G + ṽG ) v⊥ = v̄⊥ + ṽ=  ⊥  >   ṽ⊢  √ χ̃ = arctan  2 2 (v̄G +ṽG ) −ṽ⊢ We do not get into the aircraft performance modelling, since for simplicity, we consider only a horizontal straight-line path, with a constant reference course (dχ̄(t) = 0), a constant reference ground speed (dv̄G (t) = 0) and a vertical reference speed zero (v̄⊥ (t) = 0). Hence, χ̄ may be fixed to zero, so we see that the deterministic motion (i.e. without random perturbations) is described by the following equations: ? dχ̄(t) = 0, dv̄G (t) = 0, v⊢ (t) = 0, v⊥ (t) = 0, v// (t) = v̄G , where the origin of the referential is assumed to be the location of the aircraft, at time t = 0. For more details on aircraft performance modelling, see [53] for example. To 63 A Mathematical Background obtain the actual position of the aircraft, we must add the random effects, whence the use of the stochastic differential equations. Considering the NLR’s approach, we have the following SDE for the lateral deviations (which result in the stationary Double Exponential distribution): ? dx̃⊢ = ṽ⊢ dt dṽ⊢ = −a1 sign(x̃⊢ )dt − a2 ṽ⊢ dt + noise To understand these equations, we can think about the motion of a particle suspended in a fluid, which is influenced by two principal forces. First, a nonrandom (deterministic) motion which is engendered by an external force derived of the potential function −a1 sign(x̃⊢ (t)) and by the nature of the underlying fluid flow with the friction coefficient a2 . Second, collisions and more general interaction relationships with the fluid cause generally random perturbations of the velocity. Over a short time, these perturbations are often well described by Brownian fluctuations. Thus for a small duration from time t to t + ∆t, the variation of the velocity is approximated by ṽ⊢ (t + ∆t) − ṽ⊢ (t) ≈ −a1 sign(x̃⊢ (t))∆t − a2 ṽ⊢ (t)∆t + b1 Normal(0, ∆t), where Normal(0, ∆t) is a Gaussian distributed noise with expectation 0 and variance ∆t. Following [62], we state that x̃⊢ (t) is a continuous-path Markov process such that the velocity ṽ⊢ (t) exists and is continuous. Moreover, (x̃⊢ (t), ṽ⊢ (t)) has a joint density p(x, v, t), so the probability of finding x̃⊢ (t) in ]x, x + ∆X] and ṽ⊢ (t) in ]v, v + ∆v] is p(x, v, t)∆x∆v as ∆x → 0 and ∆v → 0. It is remarkable that this density converges as t → ∞ to a simple stationary distribution π(x, v) = K exp(− 2a2 a1 |x| a2 v 2 ) exp(− ), b21 b21 in which position and velocity are independent; velocity is Normally distributed, and position deviation is distributed according as a Double Exponential density with variance b41 /(2a21 a22 ). Denote respectively σx̃2⊢ and σṽ2⊢ the variance of x̃⊢ and ṽ⊢ ; these quantities can be estimated using statistical data or expert judgement. So, we obtain √ σṽ2⊢ a1 = 2 , σx̃⊢ a2 = 1 = b1 > 2 . 2 σṽ⊢ Nevertheless, we note that the model depends on three parameters a1 , a2 , b1 linked together by two equations; thus one of these parameters remains undetermined by the only one observation of the stationary measure. In fact, the friction coefficient a2 influences only on the convergence speed to equilibrium and not on the stationary measure. A drawback of this model is the singular character of the function sign which is not derivable at the origin. The choice of this function has been guided by the form of the stationary measure you desired, despite its non physical aspect. Nevertheless, it 64 A Mathematical Background seems possible to overcome this drawback using the remark, already mentioned by J. B. Parker, that the Laplace distribution with parameter a can be obtain from a centered Gaussian distribution whose the variance is distributed according to an exponential law with parameter a2 /2, i.e. P(u < σ 2 ≤ u + du) = This follows from the identity , ∞ 0 a2 − a2 u e 2 du. 2 (A.2.1) a2 − a2 u − x2 du a = e−a|x| . e 2 e 2u √ 2 2 2πu This situation can be regarded as representing a population of errors made by different aircraft with varying degrees of proficiency (different σ) and perhaps under different operating conditions (again, different σ). A.2.3 The Kramers equation Many stochastic approaches of the collision risk are based on a brownian perturbation of the position or more generally on diffusions under a potential function (see for example [84, 85]). Nevertheless diffusions under a potential function are physically unrealistic because, recalling Newton’s laws of motion, a potential really acts to cause a change in velocity rather than in position. We now describe a more realistic model; the Kramers equation, which describes the motion of a particle under the influence of the external force potential H, of viscosity, and of random perturbations of velocity. Let consider only the case of one dimensional motion and denote by Xt the position, Vt the velocity of the particle and H ′ (x) the derivative of the potential H. The Kramers equation is a stochastic differential equation of the form ? dXt = Vt dt (A.2.2) √ dVt = −H ′ (Xt )dt − γVt dt + 2γvth dBt , where Bt is a standard Brownian motion. Thus, the pair (Xt , Vt ) is a 2-dimensional diffusion and the position is a smoothed stochastic process, since the Brownian component, which is not differentiable, acts upon the velocity of the particle only. In other terms, the particle possesses a velocity though in general no acceleration. This 2-dimensional diffusion process (Xt , Vt ) has a simple stationary distribution π(x, v) = K exp(− H(x) v2 ) exp(− 2 2 ), vth 2vth in which position and velocity are independent. 65 A Mathematical Background Example 1 [Laplace stationary distribution] Let H(x) = a|x|, where a is a positive real constant. This potential function satisfies the conditions that ensure existence and uniqueness of a solution. Here, the position has a 2 ). This example is bilateral exponential stationary distribution, π(x) = K exp(−a|x|/vth interesting, since it leads to a stationary Laplace distribution, which is used in modelling the tails of the deviation errors in the air traffic data. A.2.4 Solution for a linear force For a linear force H ′ (x) = ω 2 x, the Kramers equation may be written in the form 5 6 5 6 5 6 1 x x 0 dt + 2γvth = −A d dBt , v v 1 where A= 5 0 −1 ω2 γ 6 The general solution, given X0 , V0 is 5 6 5 6 5 6 , t 1 Xt −A(t−s) 0 −At X0 dBs . =e e + 2γvth 1 Vt V0 0 Furthermore the exponential matrix e−At is given by 6 5 1 −λ2 e−λ1 t + λ1 e−λ2 t e−λ2 t − e−λ1 t −At e = . ω 2 (e−λ1 t − e−λ2 t ) λ1 e−λ1 t − λ2 e−λ2 t λ1 − λ 2 Thus, the transition probability P (x, v, t|x′ , v ′ , 0) is given by the two-variable gaussian distribution with means values m(t) = (x(t), v(t)) and covariance matrix 6 5 Σxx (t) Σxv (t) , Σ(t) = Σxv (t) Σvv (t) where x′ v′ (λ1 e−λ2 t − λ2 e−λ1 t ) + (e−λ2 t − e−λ1 t ), λ1 − λ2 λ1 − λ2 x′ ω 2 v′ (e−λ1 t − e−λ2 t ) + (λ1 e−λ1 t − λ2 e−λ2 t ), v(t) = λ1 − λ2 λ1 − λ2 > @λ + λ 2 γvth 4 = −(λ1 +λ2 )t 1 −2λ2 t A 1 −2λ1 t 1 2 , e − 1 − + e − e Σxx (t) = (λ1 − λ2 )2 λ1 λ2 λ1 + λ2 λ1 λ2 (A.2.3) = >2 2 γvth −λ1 t −λ2 t e − e , Σxv (t) = (λ1 − λ2 )2 > @ A 2 γvth 4λ1 λ2 = −(λ1 +λ2 )t −2λ1 t −2λ2 t e − 1 − λ e − λ e λ + λ + , Σvv (t) = 1 2 1 2 (λ1 − λ2 )2 λ 1 + λ2 x(t) = 66 A Mathematical Background with λ1,2 = 21 (γ ± 1 γ 2 − 4ω 2 ). For ω 2 > γ 2 /4 the real parts of the eigenvalues λ1,2 , and for ω 2 ≤ γ 2 /4 the eigenvalues λ1,2 are larger than zero. Therefore, the covariance matrix converges for t → ∞ to 0 / 2 5 2 6 vth 0 def σx 0 2 ω = . Σ(∞) = 2 0 σv2 0 vth For t ∼ 0, we obtain with η = √ 2γvth Σxx (t) ∼ η 2 t3 /3, Σxv (t) ∼ η 2 t2 /2, Σvv (t) ∼ η 2 t. Initial Gaussian velocity We assume now that X0 = 0 and V0 is a Gaussian centered random variable with variance 2 , independent of the Brownian motion. The solution is now the sum of the previous vth solution with x′ = v ′ = 0 and the centered Gaussian vector 5 6 5 6 1 0 (e−λ2 t − e−λ1 t )V0 −At = e . V0 λ1 − λ2 (λ1 e−λ1 t − λ2 e−λ2 t )V0 So the solution is always a centered Gaussian vector whose the covariance matrix is now given by   where 2 vth (1 ω  2  − y 2 (t)) 2 y(t)z(t) vth 2 y(t)z(t) vth   2 (1 − ω 2 z 2 (t)) vth y(t) = (λ2 e−λ1 t − λ1 e−λ2 t )/(λ1 − λ2 ) z(t) = (e−λ1 t − e−λ2 t )/(λ1 − λ2 ). A.2.5 Randomization approach As mentioned previously, to obtain an asymptotic Laplace distribution, instead of considering the Kramers’ equation with the potential function H(x) = a|x| order to obtain an asymptotic Laplace distribution, we can use a parabolic potential H(x) = ω 2 x2 /2 and randomize the parameter ω 2 . Since the asymptotic distribution of the position is 2 /ω 2 , we consider a random potential such that the law of Gaussian with variance vth 2 2 ) (A.2.1). Of course, we 1/ω is an Exponential distribution with parameter a2 /(2vth does not obtain necessarily the same transition density than in the previous model, and so the same transient behavior. The transition density p̃(x, v, t) of this model, called the “annealed” transition density may be obtained by integrating the Gaussian density (so called the “quenched” transition density) with the Exponential density, , ∞ 2 2 a − 2va2 u −2 th p(x, v, t; u)du, where u = ω p̃(x, v, t) = . 2 e 2v 0 th 67 A Mathematical Background A.2.6 Extreme cases We consider the two extreme cases corresponding to large and low friction constant γ (see [1][§I14]. For large γ, let Zt = Xt + γ −1 Vt . Then equations (A.2.2) yield dZt = −γ −1 H ′ (Xt )dt + 1 2/γvth dBt . By using the scaling property of Brownian motion, which states that cBt/c2 is a Brownian motion for every c > 0, we obtain that the time-changed process Ztγ of Zt is solution of √ dZu = −H ′ (Xu )du + 2vth dBu . Thus, as γ → ∞, the speeded-up process Xγt converges to the diffusion with drift 2 . µ(x) = −H ′ (x) and σ 2 (x) = 2vth As γ → 0, the motion approximates that of deterministic frictionless motion under a potential. For instance, assume that H(x) attains its minimum at a unique point x0 . If γ = 0, starting from x1 > x0 with velocity 0, the particle moves to x̂1 < x0 such that H(x̂1 ) = H(x1 ) and then returns to x1 : call this a x1 -cycle. As the energy 1 Et = H(Xt ) + Vt2 , 2 (A.2.4) is conserved, we can calculate the duration of this x1 -cycle (i.e. the trajectory on the level set Et ≡ H(x̂1 )) , x1 , x1 −1 D(x1 ) = 2 ( velocity at x) dx = 2 (2(H(x1 ) − H(x))−1/2 dx. x̂1 x̂1 In the stochastic case (γ > 0), by using It calculus, we obtain the following SDE satisfied by the energy 1 2 (A.2.5) dEt = γ(vth − Vt2 )dt + 2γvth Vt dBt . Let X̂1 , X̂2 , · · · , the right-most extremes of successive cycles, then integrating (A.2.5) over a cycle gives 2 2 H(X̂n+1 ) ≈ H(X̂n ) + γ(vth D(X̂n ) − I(X̂n )) + Normal (0, 4γvth I(X̂n )), where I(x1 ) = , V (t)dt over a x1 − cycle = 2 , x1 x̂1 (2(H(x1 ) − H(x))1/2 dx. This approximation is the starting point of the heuristic developed in [1] to estimate the hitting time by Xt of a high level. 68 A Mathematical Background A.2.7 Low-Friction Limit For the limit γ = 0, the motion of the particle is described by the trajectories of the dynamical system in the phase space (x, v), corresponding to the energy level curves. For very small friction, the energy varies slowly in the course of time (cf. (A.2.4)). So the random motion consists of fast rotations along the unperturbed trajectories of the deterministic system and slow motion across these trajectories. The nature of this motion thus suggests a set of new coordinates which splits the two components of motion: the so-called action-angle coordinates. The action part is defined by the area enclosed by the level curves of the energy, hence, it captures the slow component of the motion. Whereas the angle describes the uniform motion along the level curves, and is therefore related with the fast motion [3]. With these new coordinates (I, φ), the new energy or Hamiltonian is a constant h(I) and the angle coordinate φ increases by 2π after each complete period T (x, v) = T (I) of the motion on the level curve corresponding to the energy h(I). Let assume that H(0) = 0 and H(x) = H(−x), so, we obtain the following: , 1 x0 (I) 1 I= 2(h(I) − H(ξ))dξ π −x0 (I) , x 1 ∂ φ= 2(h(I) − H(ξ))dξ ∂I −x0 (I) , x 1 ∂ ′ = h (I) 2(h − H(ξ))dξ ∂h −x0 (I) , x 1 1 = ω(I) dξ 2(h − H(ξ)) −x0 (I) The unperturbed Hamiltonian equations give that for I > 0 dφ = ω(I); dt dI = 0, dt thus I is constant and φ(t) = ω(I)t + φ(0). For I > 0 and φ ∈] − π, π[, the Jacobian of this transformation is given by [3] 0 / 0 / v ∂x ∂x vβ(φ, I) (∂x, v) ω(I) ∂I ′ (x) := ∂φ = ω(I) ∂v ∂v ′ (∂φ, I) − Hω(I) ∂φ ∂I v − H (x)β(φ, I) where β(φ, I) = ?- φ ' π 2 1 v 2 (ξ,I) − ω ′ (I) ω(I) ]dξ, −β(−φ, I) for I > 0, φ ∈ 0, π[ for I > 0, φ ∈ 0, π[ The inverse of the Jacobian may be immediately obtained by using the identity ∂x ∂v ∂x ∂v − = 1. ∂φ ∂I ∂I ∂φ 69 (A.2.6) (A.2.7) A Mathematical Background We obtain [3] (∂φ, I) := (∂x, v) 5 ∂φ ∂x ∂I ∂x ∂phi 6 ∂v ∂I ∂v = / ω v − H ′ (x)β H ′ (x) ω −vβ v ωv 0 Finally, the Kramers’equation in action-angle variables is given by 5 2 6 Vt ∂ It H ′ (Xt ) ′ 2 Vt dt + [−γVt − H (Xt )]dt + γvth dt dIt = ω(It ) ω(It ) ∂Vt2 1 Vt + 2γvth dBt ω(It ) 1 def = γfI (φt , It )dt + 2γvth gI (φt , It )dBt where dφt = ω(It )dt − H ′ (Xt )β(φt , It )Vt dt − Vt β(φt , It )[−γVt − H ′ (Xt )]dt 5 2 6 1 2 ∂ φt dt − 2γvth Vt β(φt , It )dBt + γvth 2 ∂Vt 1 def = ω(It )dt + γfφ (φt , It )dt + 2γvth gφ (φt , It )dBt 5 2 6 5 2 6 Vt2 (φ, I) ∂ It 2 2 ∂ φt 2 , fφ (φ, I) = Vt (φ, I)β(φ, I) + vth fI (φ, I) = − + vth ω(I) ∂Vt2 ∂Vt2 Vt (φ, I) gI (φ, I) = , gφ (φ, I) = −Vt β(φt , It ). ω(It ) Accordingly, as γ ↓ 0, we see that the first motion is slow since its speed is of order γ whereas the second motion is relatively fast. Thus before the diffusing particle makes a small displacement transversally to the level sets direction it makes many rotation along the periodic trajectories. We now apply the “averaging principle” in [49][§7.9] consisting in averaging the slow motion over a period of the fast motion, giving an averaged system 1 dI¯t = γ f¯I (I¯t )dt + 2γvth ḡI (I¯t )dBt , where , > v 2 = Iω ′ (I) 1 π fI (φ, I)dφ = −I − th −1 f¯I (I) = π 0 ω(I) ω(I) , π 1 I ḡI2 (I) = gI2 (φ, I)dφ = π 0 ω(I) We used the following identity , , 1 π1 1 π 2 2(h(I) − H(ξ))dξ = Iω(I). V (φ, I)dφ = ω(I) π 0 π 0 The main result proved in [49] is that with probability 1 the process It converges to the process I¯t when γ → 0. So, if γ is so small that the effect of the perturbative term in 70 A Mathematical Background γ has a negligible effect during a period, the trajectory It can be approximated by the averaged system I¯t . Finally, we get the new Fokker-Planck equation for the probability density p(I, t) (for convenience, we use the notation It ) ∂p(I, t) I ∂p ∂ 2 ∂ = γ (Ip) + γvth . ∂t ∂I ∂I ω(I) ∂I 2 − Iω( ¯ I)) ¯ and variance Let remark that h(I¯t ) is a one-dimensional diffusion with drift γ(vth 2 ¯ I), ¯ which corresponds to the averaged equation (A.2.5). 2γvth Iω( Example 2 When the potential is parabolic, i.e. H(x) = ω 2 x2 /2, the level curve h(I) = √ √ E is an ellipse with axis 2E and ( 2E)/ω. Accordingly, the area of this ellipse being (2πE)/ω, we get I = E/ω and ω(I) = ω. So, for small γ we obtain the averaged system ∂p(I, t) ∂ 2 −1 ∂ ∂p = γ (Ip) + γvth ω I . ∂t ∂I ∂I ∂I (A.2.8) 2 )I in (A.2.8), we get Introducing the new variable I˜ = (2ω/vth ˜ t) ∂p(I, ∂ ˜ ∂ ∂p = γ (Ip) + 2γ I˜ . ∂t ∂ I˜ ∂ I˜ ∂ I˜ So, a straightforward computation gives that the stationary distribution of I˜ is the Exponential law with parameter 1/2. This is the exact distribution, since following Section 2 is distributed as the sum of two A.2.4, the stationary state of the scaled energy 2E/vth squares of independent standard Gaussian variables, whose the law is precisely the Exponential law of parameter 1/2. The transition density p(x, y, t) of the averaged diffusion I¯t satisfies the Kolmogorov backward equation v2 v2 ∂ 2p ∂p ∂p = γ( th − x) + γ th x 2 . ∂t ω ∂x ω ∂x By the method of separation of variables, we attempt a solution of the form p(x, y, t) = e−λγt ϕ(x). We obtain the following equation (we omit the fixed variable y) xϕ′′ (x) + (1 − ω ω x)ϕ′ (x) + λ 2 ϕ(x) = 0. 2 vth vth As the Laguerre polynomials Ln (x) satisfy the differential equation xL′′n (x) + (1 − x)L′n (x) + nL(x) = 0, 71 0 ≤ x < ∞, (A.2.9) A Mathematical Background ω 2 x shows that vth Ln ( vω2 x). So th the transformation x → by λ = n and ϕn (x) = p(t, x, y) = 5 the solutions of the equation (A.2.9) are given 6 ∞ @ 5 ω 6 A2 = ω > = ω > ω −nγt exp − y e L n 2 2 2 x Ln v 2 y . vth vth vth th n=0 It turns out that the density of the averaged squared velocity Vt2 is ωp(t, v0 , ωv), hence 2 (1 − e−γt ). Let compare with the exact variance Σ (t) (cf. the expectation E0 (Vt2 ) = vth vv Section A.2.4). Let assume that γ ≪ 4ω 2 , so the first order approximation of Σvv (t) gives γv 2 2 Σvv (t) ∼ vth (1 − e−γt ) + th e−γt sin(2ωt), 2ω 2 /ω 2 )(1 − e−γt ). hence, the same value by averaging. Finally, we get E0 (Xt2 ) = (vth Example 3 Now, let consider the potential function H(x) = a|x|. The action-angle coordinates may be obtained directly by (A.2.6, A.2.7) I= 5 6 2 (2E)3/2 , 3aπ ω(I) = 5 6 πa 2/3 (3I)−1/3 , 2 hence the following Fokker-Planck equation 5 √ 62/3 ∂ ∂ 4/3 ∂p ∂p(I, t) 2 2 3 = γ (Ip) + γvth I . ∂t ∂I aπ ∂I ∂I The stationary distribution is given by p(I) = > = (3β)2/3 β 2/3 I , exp − √ 2 3 2vth 23/2 πvth with β = aπ , 2 2 has a chi-squared which means that the distribution of the random variable (3βI)2/3 /vth density with 3 degrees of freedom. −2 2/3 I , the function ut (x, y) = Accordingly, let consider the transformation x = (3β)2/3 vth 1/2 x p(t, x, y) satisfies the following equation 5 6 3 ∂u ∂ 1 ∂2 = − (3 − x)u + (4xu), 2γ ∂t ∂x 2 ∂x2 whose the adjoint form is 5 6 3 ∂u ∂u ∂2u = (3 − x) + 2x 2 . 2γ ∂t ∂x ∂x 72 A Mathematical Background As previously, we attempt a solution of the form u(x, y, t) = e− 2γλ t 3 ϕ(x, y). We get the following equation xϕ′′ (x) + (3/2 − z/2)ϕ′ (x) + λ/2ϕ(x) = 0. (A.2.10) By substituting x/2 for x, the equation (A.2.10) becomes an equation of the form zϕ′′ (z) + (3/2 − z)ϕ′ (z) + λϕ(z) = 0. A corresponding solution is for λ = n ϕ(x) = c1 L1/2 n (z) + c2 F (−n, 3/2, z), where Lαn (z) are the generalized Laguerre polynomials n 5 6 Γ(α + n + 1) 2 n (−1)k α Ln (z) = zk , Γ(n + 1) k Γ(λ + k + 1) k=0 and F (α, γ, z) is the degenerated hypergeometric function F (α, γ, z) = 1 + α(α + 1) z 2 αz + + ··· . γ 1! γ(γ + 1) 2! As F (−n, b, z) ∝ z n as z → ∞, we have c2 = 0. The Laguerre system Lαn (z) comprises the unique orthogonal polynomials with respect H the weight function w(z) = z α e−z Γ(α + 1), pour z > 0; w(0) = 0 for z < 0, since the identity 5 6 , ∞ n+α α α Ln (z)Lm (z)w(z)dz = δnm . n 0 For example, we have Lαn (0) = 5 6 n+α , α Lα0 (z) = 1, Lα1 (z) = −z + α + 1. H 2 2 is given Finally, the transition density of the diffusion Zt = (3βIt )2/3 vth = (2E t )/vth by ∞ 2 2γnt 1 1/2 p(t, x, y) = 3/2 y 1/2 e−y/2 e− 3 L1/2 n (x/2)Ln (y/2)wn , 2 Γ(3/2) n=0 % & with wn−1 = n+1/2 . n A straightforward but arduous separation of variables [62][§15.13] leads to 7 (x + y)e−bt 8= 1 = 1xye−bt > >−1/4 y 1/2 e−y/2 exp − I1/2 e−bt xy , p(t, x, y) = 3/2 4 1 − e−bt 1 − e−bt 2 Γ(3/2)(1 − e−bt ) 73 A Mathematical Background where b = 2γ/3 and I1/2 is the modified Bessel function of order 1/2. It follows that given E 0 = e Ee (E t ) = 2 % 2γt 3vth 2e & 1 − e− 3 (1 − 2 ) . 2 3vth Furthermore, let notice that V 2 t = ω(It )It = (2/3)E t , so a|Xt | = (2/3)E t = 2V 2 t . We deduce that ' 2γt 2e ( 2 1 − e− 3 (1 − 2 ) , Ee (V 2 t ) = vth 3vth Ee (|X|t ) = 2 % 2γt vth 2e & 1 − e− 3 (1 − 2 ) . a 3vth A.2.8 High Friction In Section A.2.6, we have seen that for large γ the diffusion Zt = Xt + γ −1 Vt ≈ Xt was a solution of 1 dZt = −γ −1 H ′ (Xt )dt + 2γ −1 dBt . By using scaling property of Brownian motion, we obtain that the time-changed process Ztγ/2v2 is solution of th a (A.2.11) dZu = − 2 sgn(Zu )du + dBu , 2vth and accordingly, for large γ, we can approximate Xu by Zu . To obtain the transition density, we follow [61][§6.5]. Firstly, we consider the new function sgn defined by ? 1; x>0 , sgn(x) = −1 ; x ≤ 0 which is left continuous (this is nothing else that the left-derivative of |x|). Let intro2 ). By Girsanov transformation, duce the notation u(x) = −θ sgn(x), with θ = a/(2vth equation (A.2.11) has a weak solution (Z, B), (Ω, Ft , P̃x ) which is unique in the sense of probability law, with finite-dimensional distribution given by P̃x [(Xt1 , · · · , Xtn ) ∈ Γ] , @ 7, t 8A 1 t 2 = Ex 1{(Bt1 ,··· ,Btn )∈Γ} exp u(Bs )dBs − u (Bs )ds , 2 0 0 where Ex is the expectation with respect the Brownian probability Px . The transition density p̃(x, z, t) takes the form , @ 7, t 8A 1 t 2 u(Bs )dBs − p̃(x, z, t) = Ex 1{Bt ∈dz} exp u (Bs )ds . (A.2.12) 2 0 0 To eliminate the stochastic integral in (A.2.12), let set ? , z θz; z<0 u(y)dy = f (z) = −θz; z ≥ 0 0 74 A Mathematical Background The It formula generalized for convex functions [61][Th 3.6.22] gives , t a u(Bs )dBs + 2 L(t), f (Bt ) = f (B0 ) + vth 0 where L(t) is the local time of B at the origin Lt (0) = lim ǫ↓0 1 meas{0 ≤ s ≤ t; |Bt | ≤ ǫ}, 4ǫ where meas(A) denotes the Lebesgue measure of the set A. On the other hand, , t a2 u2 (Bs )ds = 4 t, 4vth 0 and (A.2.12) becomes , @ t a2 ( ∞ a p̃(x, z, t) = exp f (z) − f (x) − exp{b 2 }Px [Bt ∈ dz; L(t) ∈ db] 4 2 4vth v b=0 th Finally, we get [61][Remark 6.5.2]  8 @ 7 (x−z−θt)2 1 −2θz ∞  √ +θe exp −   2t x+z 2πt    8 @ 7 p̃(x, z, t) = (x−z+θt)2 1 2θz ∞  √ +θe exp 2θx −  2t  x−z 2πt    7 exp − (v−θt)2 2t x ≥ 0, z ≤ 0 Accordingly a tedious computations gives for the function v(t, x) = Ẽx (Zt2 ) I 7 (|x| − θt)2 8 1 t 1 (|x| − θt − ) exp − v(t, x) = 2 + 2θ 2π θ 2t 7 8 = 1 |x| − θt > √ + (|x| − θt)2 + t − 2 Φ 2θ t = |x| >@ = |x| + θt >A 1 √ +t− 2 1−Φ + e2θ|x| θ 2θ t Φ(x) := , x e−y −∞ We can verify that lim v(t, x) = x2 , t↓0 2 /2 dy √ . 2π lim v(t, x) = t→∞ 2v 4 1 = 2th . 2 2θ a On the other hand v(t, x) satisfies the equation ∂v 1 ∂2v ∂v = − θsgn(x) . 2 ∂t 2 ∂x ∂x 75 8 A dv ; x7≥ 0, z > 0 8 A 2 exp − (v−θt) dv ; 2t and for x ≤ 0, we use the relation p̃(x, z, t) = p̃(−x, −z, t). where (A.2.13) A Mathematical Background A.3 Rare Event Simulations Let assume that the dynamic of the studied system is well described by a continuoustime strong Markov process X = (Xt )t≥0 with values in Rd for some d > 0. We suppose that the paths of X are cdlg (i.e. right-continuous with left-hand limits). Let B ⊂ Rd be some closed critical region, in which the system could enters but with a very small probability, say about 10−9 . Let TB = inf{t ≥ 0 : Xt ∈ B} the entrance time into B, our objective is to compute the probability of the critical event B and the probability distribution of the critical paths PB = P[TB ≤ T ], E[f (Xt ; 0 ≤ t ≤ TB )|TB ≤ T ], where T is a finite deterministic time, or an a.s. finite stopping time. In practice, none of the simulated trajectories will ever reach B, hence a crude Monte Carlo simulation fails. To speed-up the simulations, there exists two main methods: Importance Sampling and Importance Splitting. Importance Sampling works by changing the stochastic process X into a new process for which the set B has a high probability; this makes the Monte-Carlo estimator of the probability more efficient, mainly by decreasing its variance. Nevertheless, this estimator is biased, so to unbiased it we need to weight the sample values to get to account for the change in stochastic process we are simulating. The well-know problem to Importance Sampling method is that it is hard to derive the optimal change and in practice usually this optimal change cannot be found analytically and/or adaptively. In contrast to importance sampling type algorithms, in the trajectory splitting methodology, the step-by-step evolution of the system follows the original probability measure. Entering the intermediate states, which is usually characterized by crossing a threshold by a control parameter, triggers the splitting of the trajectory. The current system state is held, and a number of independent subtrajectories are simulated from that state. For example, let us consider (n + 1) sets Bi such that Bn ⊂ Bn−1 ⊂ · · · ⊂ B0 . When the rare event B coincide with Bn , we have the product formula P(B) = P(B|Bn−1 )P(Bn−1 |Bn−2 ) · · · P(B1 |B0 )P(B0 ). (A.3.1) On the right hand side of (A.3.1), each conditioning event is “not rare”. The branching splitting technique proceeds as follows. Make a {0, 1} Bernoulli trial to check whether or not the set event B1 has occurred. If B1 has occurred, then we split this trial in R1 Bernoulli additional trials, and for each of them we check again wether or not the event B2 has occurred. This procedure is repeated at each level. More precisely, each time 76 A Mathematical Background the event Bi has occurred, we sample Ri trials and we repeat this splitting technique each time Bi+1 has occurred. If an event level is not reached, neither is B , then we stop the current retrial. Using R0 independent replications of this procedure, we have then considered R0 R1 · · · Rm trials, taking into account for example, that if we have failed to reach a level Bi at the i-th step, the Ri · · · Rm possible retrials have failed. An unbiased estimator of P(B) is clearly given by the quantity PJ = R0 NB .m i=1 Ri , where NB is the total number of trajectories having reached the set B . It can be proven [73] that in some sense the optimal simulation is obtained if m = ⌊−0.6275 log P (B) − 1⌋, P(Bi |Bi−1 ) ≈ 1/5, Ri = 5. Nevertheless, in practice the trajectory splitting method may be difficult to apply. For example, the case of the estimation of the probability of a rare event in dynamical system is more complex, since the difficulty to find theoretically the optimal Bi , and Ri for each level i . Furthermore, the probability to reach Bi varies generally with the state of entrance in level Bi−1 . Finally, but not the least, the conditional probabilities P(Bi |Bi−1 ) are of course generally unknown! Hence, the random evolution of the process must be taken into account; and for that a well adapted tool is the so-called Feynman-Kac models. A.3.1 Multilevel Feynman-Kac distributions We consider a continuous-time strong Markov process X = {X(t); t ≥ 0} with value in Rd for some d > 0 and with cdlg trajectories. Let introduce an embedded sequence of closed regions [0, ∞) × B = An ⊂ · · · ⊂ A1 ⊂ A0 = [0, ∞) × Rd , which are defined by Ak = {(t, x) ∈ [0, ∞) × Rd : Fk (t, x) ≤ ck }, where Fk is some lower semi-continuous function. Let the corresponding hitting times Tk = inf{t ≥ 0 : (t, X(t)) ∈ Ak }, which satisfy 0 = T0 ≤ T1 ≤ · · · ≤ Tn = TB . To capture the precise behavior of the process X in each region, we consider the random excursions Xk of X between the successive random times Tk−1 and Tk . More precisely, 77 A Mathematical Background we introduce the discrete-time Markov chain X = {Xk , k = 1, · · · , n} with value in the excursion set E, defined by Xk = ((t, X(t)), Tk−1 ∧ T ≤ t ≤ Tk ∧ T ) , with t ∧ T = inf{T, t}. We observe that these excursions may have different random lengths. To check wether or not a given path e = ((t, X(t)), t′ ≤ t ≤ t′′ ) , starting at X(t′ ) ∈ Bk−1 at time t′ , has succeeded to reach the level Bk at time t′′ , it is convenient to introduce the terminal point π(e) = (t′′ , X(t′′ )) of the excursion and the indicator functions gk defined by gk (e) = 1{π(e)∈Ak } . With this notation, and for each k, we have Tk ≤ T iff gk (Xk ) = 1, hence 1{Tk ≤T } = gk (Xk ) = k * gp (Xp ). p=1 Now, let us introduce the so-called Feynman-Kac distributions γk on the path space E defined in such a way that the integral of all bounded measurable functions f relatively to this measure, denoted by γk (f ) is defined by the formula   k * ( ' gp (Xp ) = E f ((t, X(t)), Tk−1 ≤ t ≤ Tk )1{Tk ≤T } . γk (f ) = E f (Xk ) p=1 This measure is not a probability, since its normalizing constant is given by γk (1) = P[Tk ≤ T ]. Therefore, the corresponding probability µk is given by µk (f ) = γk (f ) = E [f ((t, X(t)), Tk−1 ≤ t ≤ Tk )|Tk ≤ T ] . γk (1) The entrance distribution with support included in Ak is the image of µk under π, i.e. the measure πk ◦ π −1 defined by µk ◦ π −1 (dt, dx) = P [Tk ∈ dt, X(Tk ) ∈ dx)] . Thus, for any bounded measurable function φ defined on Ak , we get µk ◦ π −1 (φ) = µk (φ ◦ π) = E [φ(Tk , X(Tk ))|Tk ≤ T ] . Similarly, introducing the measure γk− defined by   k−1 * ' ( γk− (f ) = E f (Xk ) gp (Xp ) = E f ((t, X(t)), Tk−1 ≤ t ≤ Tk ∧ T )1{Tk−1 ≤T } , p=1 78 A Mathematical Background leads us to define the so-called Feynman-Kac distributions ηk obtained by normalizing γk− , and so given by ηk (f ) = E [f ((t, X(t)), Tk−1 ≤ t ≤ Tk ∧ T )|Tk−1 ≤ T ] . In particular for f = gk , ηk (gk ) = P[Tk ≤ T |Tk−1 ≤ T ] := Pk , hence P[Tk ≤ T ] = γk (1) = k * ηp (gp ) = p=1 k * Pp . p=1 The interacting particles methods, we will discuss hereafter, provide numerical approximation for the rare event probability PB = P[TB ≤ T ] = γn (1), the transition probabilities from one level to the next Pk and the entrance distribution in one level µk ◦ π −1 . A.3.2 Interacting particle system approximations The idea of the Interacting Particle Systems (IPS) is derived from the following remarks. First at all, the so-called Feynman-Kac flow (ηk ; 0 ≤ k ≤ n) is the solution of a nonlinear measure-valued dynamical system ηk = Φk (ηk−1 ). (A.3.2) The mappings Φk from the set of measures Pk (E) = {η : η(gk ) > 0} into the set P(E) of measure on E are defined by Φk (η)(f ) = Ψk−1 (η) (Mk f ) , where the non-homogeneous Markov kernels Mk , which describe the Markovian transitions of the Markov chain X , are defined for all excursion e by Mk f (e) = E [f (Xk )|Xk−1 = e] , and the updating mappings Ψk from Pk (E) into Pk (E) are defined for any bounded function f by Ψk (η)(f ) = η(f gk )/η(gk ). 79 A Mathematical Background All these follow from  γk− (f ) = E gk−1 (Xk−1 ) k−2 * p=1  − (gk−1 Mk f ). gp (Xp )E [f (Xk )|Xk−1 ] = γk−1 Thus, we see that the recursion (A.3.2) involves two separate selection/mutation transitions selection mutation ηk ∈ Pk (E) −−−−−→ ηJk := Ψk (ηk ) ∈ P(E) −−−−−→ ηk+1 = ηJk Mk+1 ∈ P(E). Secondly, from a pure mathematical point of view, particles methods can be interpreted as a kind of stochastic linearization technique for solving nonlinear equations in measure space. The idea is to associate to (A.3.2) a sequence ξ = (ξ 1 , · · · , ξ N ) of N E-valued Markov chains such that the empirical measures ηN = N 1 2 δξ i N i=1 converge as N → ∞ to the desired distribution η. The state components of the E N valued Markov chains are called particles. Now, we describe in more details the evolution of the particles according to the diagram selection mutation 1 N ξk = (ξk1 , · · · , ξkN ) −−−−−→ ξJk = (ξJk1 , · · · , ξJkN ) −−−−−→ ξk+1 = (ξk+1 , · · · , ξk+1 ). At time k = 0, the initial configuration consists in N independent and identically distributed random variables ξ0i , i = 1, · · · , N with common law η0 . Since we have g0 (x) = 1 for η0 -almost every x, we may discard the selection at time k = 0 and set ξJ0i = ξ0i for i = T i = 0 and we introduce a cemetery point ∆ each i. We use the conventions T−1 0 which takes into account the possible stopping of the algorithm.. The mutation transition ξJk → ξk+1 at time k + 1 is defined as follows. If ξJk = ∆, we set ξJk+1 = ∆. Otherwise during mutation, independently of each other, each selected particle % & i ξJki = (t, X i (t)), Tk−1 ∧ T ≤ t ≤ Tki ∧ T evolves randomly according to the Markov transition Mk+1 , so that % & i i ξk+1 = (t, X i (t)), Tki ∧ T ≤ t ≤ Tk+1 ∧T , where i i Tk+1 = inf{t ≥ Tki ; ξk+1 ∈ Ak+1} . Thus, for each particle i, we start a trajectory from ξJki at time Tki , and let it evolve i . randomly as a copy of the process {Xs , s ≥ Tki }, until the stopping time Tk+1 80 A Mathematical Background The selection transition ξk → ξJk is defined as follows. From the N particles ξki only some of these have succeeded to reach the desired set Ak ; we denote by IkN the set of their labels. If IkN = ∅, then none of the particles have succeeded to reach the desired region; the algorithm is stopped and ξJk = ∆. Otherwise, we draw N particles ξJki uniformly among the |IkN | pieces of trajectories {ξki , i ∈ IkN }. Consequently, we obtain N independent random variables ξJki with common distribution N % & 2 gk (ξ i ) Ψk ηkN = !N k j δξki , j=1 gk (ξk ) i=1 1 2 = N δ i. |Ik | N ξk i∈Ik where ηkN = N 1 2 δξ i . k N i=1 An alternate scheme for the selection step, with less randomness in some sense, can be proposed. The key idea is to notice that the updating mapping Ψk can be rewritten in the following form , ′ ′ Ψk (η)(de ) = (ηSk (η))(dx ) := η(de)Sk (η)(e, de′ ), where Sk (η)(e, de′ ) = (1 − gk (e))Ψk (η)(de′ ) + gk (e)δe (de′ ). We easily verify that (ηSk (η))(f ) = (1 − η(gk ))Ψk (η)(f ) + η(f gk ) = Ψk (η)(f ). In this notation, (A.3.2) can be rewritten as ηk+1 = ηk Kk+1 (ηk ), with the composite Markov kernel Kk+1 (η) defined by , ′ ′ Kk+1 (η)(e, de ) = (Sk (η)Mk+1 )(e, de ) = Sk (η)(e, de′′ )M(e′′ , de′ ). E Note that the corresponding evolution equation is now decomposed into two separate transitions selection mutation ηk −−−−−→ ηJk = ηk Sk (ηk ) −−−−−−→ ηk+1 = ηJk Mk+1 , (A.3.3) In the alternative N -particle model associated with this new description each particle ξJki is sampled according to the selection distribution % & Sk ηkN (ξki , de) % N& = 1{X i (T i )∈Bk } δξi (de) + 1{X i (T i )∈B / k } Ψn ηk (de) . k k 81 k A Mathematical Background More precisely, if the i-th particle has reached the desired region, then we set ξJki = ξki . In the opposite case, the particle has not reached the k-th region. In this case, ξJki is chosen randomly and uniformly in the set {ξkj ; j ∈ IkN } of all excursions having entered into Bk . In other words, each particle that doesn’t enter into the k-th level is killed, and instantly a different particle in the Bk+1 level splits into two offsprings. We denote by τ N the lifetime of the N -particle algorithm, i.e. τ N = ∞ or τ N is the first time all the particles are killed. For each k < τ N , the N -particle approximating measure γkN defined by k−1 * γkN (f ) = ηkN (f ) ηpN (gp ), p=0 when applied to the constant function 1 gives N γk+1 (1) = γkN (gk ) k * |IpN | = . N p=1 Therefore, if τ N > n, we obtain the following approximation of the rare event probability PB = P[TB ≤ T ] = γn (1) ≈ γnN (1) n * |IpN | . = N p=1 This estimator is unbiased and the central limit theorem applies [31]: as N → ∞ 5 N 6 √ γn (1) N − 1 =⇒ N (0, σn2 ), PB where the asymptotic variance of the estimator is given by 6 2 n n 5 2 1 Varµk ◦π−1 (∆nk ) 1 2 −1 + , σn = Pk Pk E2µ ◦π−1 (∆nk ) k=1 k=1 k for the first algorithm, and 6 2 n n 5 2 1 Varµk ◦π−1 (∆nk ) 1 2 −1 + (1 − Pk2 ), σn = Pk Pk E2µ ◦π−1 (∆nk ) k=1 k=1 k for the alternate one, where ∆nk (t, x) = P [TB ≤ T |Tk = t, X(Tk ) = x] . In particular, if for each k the function ∆nk is constant, i.e. P [TB ≤ T |Tk , X(Tk )] does not depend on the hitting time and point of the level set Ak , given Tk ≤ T , then the second right-hand term is null and the asymptotic variance reduces to the expression 6 n 5 2 1 2 σn = −1 , Pk k=1 82 A Mathematical Background as given in [73]. Ideally, the level set Ak should be chosen such that the functions ∆nk are constant; even if this is clearly unrealistic for most practical problems, this observation gives an insight on how to choose them. To see that the asymptotic variance is better than the one obtained with a crude Monte Carlo simulation, let us consider the simplified framework of a one dimensional Markov process and of a N -particle model in which all the probabilities Pk are equal to the same P . In this case, the asymptotic variance of the estimator is nPB2 (1 − P )/P , whereas in naive Monte Carlo the variance of the estimator is PB (1 − PB ). It is also worth noting that the n factor does not mean that the variance increases with n, since for a given 1/n problem the rare event probability PB is fixed. So P ≈ PB , and we have n 1−P 1 1 1 ≈ n exp{− log PB − 1} ≈ − log PB + log2 PB + o( 2 ), P n 2n n which means that as n → ∞, the variance is decreasing to −PB2 log PB . Finally, let us mention the two adaptive algorithms [24, 74] developed to cope with the possible extinction of the particle system and the difficulty in determining the splitting levels at hand. 83 Bibliography [1] D Aldous. Probability Approximations via the Poisson Clumping Heuristic, volume 77 of Applied Mathematical Sciences. Springer Verlag, 1989. [2] J.W. Andrews, J.D. Welch, and H. Erzberger. Safety analysis for advanced separation concepts. In Proceedings of USA/Europe ATM R&D Seminar, Baltimore, MD, 27–30 June 2005. [3] L. Arnold, P. Imkeller, and N. Sri Namachchivaya. The Asymptotic Stability of weakly perturbed two dimensional Hamiltonian Systems. In Proceedings of the 18-th Biennal of Conference on Mechanical Vibration and Noise ASME Design Engineering Technical Conferences, Pittsburgh, Pennsylvania, September 9-12 2001. [4] B. G. J. Bakker, H. Kremer, and H. A. P Blom. Geometric and probabilistic approach towards conflict prediction in free flight, 1999. [5] G.J. Bakker and H.A.P Blom. Air traffic collision risk modelling. In Proceedings of the 32nd Conference on Decision and Control, volume 2, pages 1464–1469, San Antonio, Texas, December 1993. [6] H. A. P. Blom, G. J. Bakker, P. J. G. Blanker, J. Daams, M. H. C. Everdij, and M.B. Klompstra. Accident risk assessment for advanced atm. In Proceedings of the 2nd U.S.A/Europe Air Traffic Management R&D seminar, Orlando, December 1998. [7] H. A. P. Blom and J. Lygeros, editors. Stochastic Hybrid Systems: Theory and Safety Critical Applications. LNCIS series. Springer, Berlin, July 2006. [8] H.A.P. Blom and G.J. Bakker. Conflict probability and incrossing probability in air traffic management. In Proceedings IEEE Conference on Decision and Control, pages 2421–2426, Las Vegas, NV, December 2002. [9] H.A.P. Blom, G.J. Bakker, P.J.G. Blanker, J. Daams, M.H.C. Everdij, and M.B. Klompstra. Air Transport Systems Engineering, chapter Accident risk assessment for advanced air traffic management, pages 463–480. AIAA, 2001. [10] H.A.P. Blom, G.J. Bakker, M.H.C. Everdij, and M.N.J. van der Park. Collision risk modeling of air traffic. In Proc. European Control Conf., Cambridge, UK, 2003. [11] H.A.P. Blom, G.J. Bakker, B. Klein Obbink, and M.B. Klompstra. Free flight safety risk modelling and simulation. In Proceedings of International Conference on Research in Air Transportation (ICRAT), Belgrade, 26–28 June 2006. 84 Bibliography [12] H.A.P. Blom, K.M. Corker, and S.H. Stroeve. On the integration of human performance and collision risk simulation models of runway operation. In Proceedings of the 6th USA/Europe Air Traffic Management R&D Seminar, pages 27–30, Baltimore, MD, June 2005. [13] H.A.P. Blom, J. Daams, and H.B. Nijhuis. Air Transportation Systems Engineering, volume 193 of Progress in Astronautics and Aeronautics, chapter Human cognition modeling in Air Traffic Management safety assessment, pages 481–511. 2001. [14] H.A.P. Blom, M.B. Klompstra, and G.J. Bakker. Accident risk assessment of simultaneous converging instrument approaches. Air Traffic Control Quarterly, 11:123– 155, 2003. [15] H.A.P. Blom and S.H. Stroeve. Multi-agent situation awareness error evolution in air traffic. In Proc. 7th Conference on Probabilistic Safety Assessment & Management, Berlin, Germany, 2004. [16] H.A.P. Blom, S.H. Stroeve, M.H.C. Everdij, and M.N.J. Van der Park. Human cognition performance model to evaluate safe spacing in air traffic. Human Factors and Aerospace Safety, 2:59–82, 2003. [17] M. L. Bujorianu. Extended stochastic hybrid systems. In O. Mahler and A. Pnuelli, editors, Proceedings of Hybrid Systems Computation and Control, LNCIS number 2993, pages 234–249, Berlin, 2004. Springer. [18] M. L. Bujorianu and J. Lygeros. Reachability questions in piecewise deterministic markov processes. In O. Mahler and A. Pnuelli, editors, Proceedings of Hybrid Systems Computation and Control, LNCIS number 2623, pages 126–140, Berlin, 2003. Springer. [19] M. L. Bujorianu and J. Lygeros. Stochastic Hybrid Systems: Theory and Safety Critical Applications, chapter Towards a general theory of stochastic hybrid systems, pages 3–30. LNCIS series. Springer, Berlin, 2006. [20] L Burt. Results of route space calculations, progress report. Eurocontrol report NSSG/22-DP/5, Brussels, December, 1995. [21] CAA. Outline of a method for the determination of separation standards for future air traffic systems. CAA paper 91010, DORA report 9111, CAA London, June 1991. [22] C. G. Cassandras and S. Lafortune. Introduction to Discrete Event Systems. Kluwer Academic Publishers, Boston, 1999. [23] F. Cérou, P. Del Moral, F. Le Gland, and P. Lezaud. Limit theorems for the multilevel splitting algorithm in the simulation of rare events. In M.E. Kuhl, N.M. Steiger, F.B. Armstrong, and J.A. Joines, editors, Proceedings of the 2005 Winter Simulation Conference, Orlando, December 2005. 85 Bibliography [24] F. Cérou and A. Guyader. Adaptive multilevel splitting for rare event analysis. To appear in Stochastic Analysis and its Applications, 2007. [25] F. P. Cérou, P. Del Moral, F. Le Gland, and P. Lezaud. Genetic genealogical models in rare event analysis. Publications du Laboratoire de Statistiques et Probabilites, Toulouse III, 2002. [26] S. Cohen and S. Hockaday. A Concept Paper for Separation Safety Modeling, an FAA/EUROCONTROL Cooperative Effort on Air Traffic Modeling for Separation Standards. FAA and EUROCONTROL, Brussels, May 1998. [27] K. Corker. Cognitive Engineering in the Aviation Domain, chapter Cognitive Models & Control: Human & System Dynamics in Advanced Airspace Operations. Lawrence Earlbaum Associates, Hillsdale, New Jersey, 2000. [28] R. David and H. Alla. Petri nets for the modeling of dynamic systems - a survey. Automatica, 30(2):175–202, 1994. [29] M.H.A. Davis. Markov models and optimization. Chapman & Hall, London, 1993. [30] E. De Santis, M. D. Di Benedetto, S. Di Gennaro, A. D. ’Innocenzo, and G. Pola. Stochastic Hybrid Systems: Theory and Safety Critical Applications, chapter Critical observability of a class of hybrid systems and application to air traffic management, pages 141–170. LNCIS series. Springer, Berlin, 2006. [31] P. Del Moral. Feynman-Kac Formulae. Genealogical and Interacting Particle Systems with Applications. Probability and its Applications. Springer-Verlag, New York, 2004. [32] P. Del Moral and P. Lezaud. Branching and interacting particle interpretation of rare event probabilities. In H. A. P. Blom and J. Lygeros, editors, Stochastic Hybrid Systems: Theory and Safety Critical Applications, LNCIS series, pages 351–389. Springer-Verlag, Berlin, 2006. [33] D. V. Dimarogonas, S. G. Loizou, and K. J. Kyriapoulos. Stochastic Hybrid Systems: Theory and Safety Critical Applications, chapter Multirobot navigation functions II: towards decentralization, pages 209–256. LNCIS series. Springer, Berlin, 2006. [34] DNV. Safety assessment of P-RNAV route spacing and aircraft separation. Final report TRS 052/01, Eurocontrol, 2003. [35] A. Doucet, N. de Freitas, and N. Gordon. Sequential Monte Carlo Methods in Practice. Springer-Verlag, 2001. [36] M.R. Endsley. Towards a theory of situation awareness in dynamic systems. Human Factors, 37(1):32–64, 1995. 86 Bibliography [37] H. Erzberger. Transforming the NAS: The next generation air traffic control system. In Proceedings of the 24th Int. Congress of the Aeronautical Sciences (ICAS), Yokohoma, Japan, 2004. [38] M. H. C. Everdij and H. A. P. Blom. Bias and uncertainty in accident risk assessment. NLR report CR-2002-137, National Aerospace Laboratory NLR, 2002. [39] M. H. C. Everdij and H. A. P. Blom. Analysis and design of hybrid systems, chapter Petri-nets and hybrid-state Markov processes in a power-hierarchy of dependability models, pages 313–318. Elsevier, 2003. [40] M. H. C. Everdij and H. A. P. Blom. Petri nets and hybrid state Markov processes in a power-hierarchy of dependability models. In Proceedings of IFAC Conference on Analysis and Design of Hybrid Systems, pages 355–360, Saint-Malo Brittany, France, June 2003. [41] M. H. C. Everdij and H. A. P. Blom. Piecewise deterministic markov processes represented by dynamically coloured petri nets. Stochastics, 77:1–29, 2005. [42] M. H. C. Everdij and H. A. P. Blom. Stochastic Hybrid Systems: Theory and Safety Critical Applications, chapter Hybrid Petri nets with diffusion that have into-mappings with generalised stochastic hybrid processes, pages 31–64. LNCIS series. Springer, Berlin, 2006. [43] M. H. C. Everdij, H. A. P. Blom, and G. J. (Bert) Bakker. Modeling lateral spacing and separation for airborne separation assurance using petri nets. Simulation, Transactions of the Society for Modeling and Simulation International, 82, 2006. [44] M. H. C. Everdij, H. A. P. Blom, and S. H. Stroeve. Structured assessment of bias and uncertainty in Monte Carlo simulated accident risk. In Proceedings of the 8th Int. Conf. on Probabilistic Safety Assessment and Management (PSAM8), New Orleans, LA, May 2006. [45] M. H. C. Everdij, M. B. Klompstra, H. A. P. Blom, and B. Klein Obbink. Stochastic Hybrid Systems: Theory and Safety Critical Applications, chapter Compositional specification of a multi-agent system by stochastically and Dynamically Coloured Petri Nets, pages 325–350. LNCIS series. Springer, Berlin, 2006. [46] FAA/Eurocontrol. A concept paper for separation safety modeling. Cooperative R&D Action Plan 3 report. Available at http://www.faa.gov/asd/iaor/pdf/cpcomplete.pdf, 20 May 1998. [47] FAA/Eurocontrol. Principles of operations for the use of asas. Cooperative R&D Action Plan 1 report, Version 7.1, 2001. [48] FAA/Eurocontrol. Safety and asas applications. Co-operative R&D Action Plan 1 report, version 4.1, 2004. 87 Bibliography [49] M.I. Freidlin and A.D. Wentzell. Springer-Verlag, New-York, 1984. Random perturbations of dynamical systems. [50] P. Glasserman. Monte Carlo methods in financial engineering, volume 53 of Stochastic Modeling and Applied Probability. Springer, New York, NY, 2003. [51] P.J. Haas. Stochastic Petri Nets, Modeling, Stability, Simulation. Springer-Verlag, New York, Simulation. [52] J. Hoekstra. Designing for Safety, the Free Flight Air Traffic Management concept. PhD thesis, Delft University of Technology, November 2001. [53] E Hoffman. Contribution to Aircraft Performance Modeling for ATC. Air Traffic Control Quaterly, 2(2):103 –130, 1994. [54] D.A. Hsu. The evaluation of aircraft collision probabilities at intersecting air routes. J. of Navigation, 34:78–102, 1981. [55] J. Hu, M. Prandini, and S. Sastry. Probabilistic safety analysis in three-dimensional aircraft flight. In Proc. 42nd IEEE CDC, Maui, USA, December 2003. [56] ICAO. Review of the general concept separation panel. 6th meeting, 28 November15 December, Montreal. Doc 9536, RGCSP/6 Volume1, December 1988. [57] ICAO. Manual on airspace planning methodology for the determination of separation minima. ICAO Doc. 9689-AN/953, 1998. [58] ICAO. Airborne separation assistance system (asas) circular. Draft, version 3, SCRSP, WGW/1 WP/5.0, International Civil Aviation Organization, May 2003. [59] K. Jensen. Coloured Petri Nets: Basic Concepts, Analysis Methods and Practical Use, volume 1. Springer, London, UK, 1992. [60] H.H. De Jong. Guidelines for the identification of hazards; how to make unimaginable hazards imaginable? National Aerospace Laboratory NLR, Contract report for Eurocontrol, NLR-CR-2004-094, 2004. [61] I. Karatzas and S. Shreve. Brownian Motion and Stochastic Calculus. Springer, 2000. [62] S. Karlin and H. Taylor. A Second Course in Stochastic Processes. Academic Press, New York, 1982. [63] T. Kletz. Hazop and Hazan; identifying and assessing process industry hazards. The Institution of Chemical Engineers, 4th edition, 1999. [64] J. Krozel. Free flight research issues and literature search. Under NASA contract NAS2-98005, 2000. 88 Bibliography [65] J. Krystul and H. A. P. Blom. Monte Carlo simulation of rare events in hybrid system. Hybridge Report D8.3. Available at http://www.nlr.nl/public/hostedsites/hybridge, June 2004. [66] J. Krystul and H. A. P. Blom. Sequential Monte Carlo simulation of rare event probability in stochastic hybrid systems. In Proceedings of the 16th IFAC World Congress, Prague, Czech Republic, June 4-8 2005. [67] J. Krystul and H. A. P. Blom. Sequential Monte Carlo simulation for the estimation of small reachability probabilities for stochastic hybrid systems. In Proceedings of IEEE-EURASIP Int. Proceedings of IEEE-EURASIP Int. Symposium on Control, Communications and Signal Processing, Marrakech, Morocco, March 13-15 2006. [68] J. Krystul, H. A. P. Blom, and A. Bagchi. Stochastic hybrid processes as solutions to stochastic differential equations. In C. G. Cassandras and J. Lygeros, editors, Stochastic Hybrid Systems: Recent Developments and Research Trends. CRC Press, 2006. [69] J. Krystul and H.A.P. Blom. Generalized stochastic hybrid processes as strong solutions of stochastic differential equations. Hybridge Report D2.3. http://www.nlr.nl/public/hosted-sites/hybridge/, June 2005. [70] H. Kumamoto and E.J. Henley. Probabilistic Risk Assessment and management for engineers and scientists. IEEE Press, 1996. [71] A. B. Kurzhanski and P. Varaiya. On reachability under uncertainty. SIAM Journal on Control and Optimization, 41:181–216, 2002. [72] P. E. Labeau, C. Smidts, and S. Swaminathan. Dynamic reliability: Towards an integrated platform for probabilistic risk assessment. Reliability Engineering and System Safety, 68:219–254, 2000. [73] A. Lagnoux. Rare Events Simulation. PEIS, 20(01):45–66, January 2006. [74] F. Le Gland and N. Oudjane. A sequential algorithm that keeps the particle system alive. In H. A. P. Blom and J. Lygeros, editors, Stochastic Hybrid Systems : Theory and Safety Critical Applications, LNCIS series, pages 351–400. Springer, Berlin, 2006. [75] M. R. Leadbetter, G. Lindgren, and H. Rootzen. Extremes and Related Properties of Random Sequences and Processes. Springer-Verlag, New-York, 1983. [76] A. Lecchini, W. Glover, J. Lygeros, and J. Maciejowski. Stochastic Hybrid Systems: Theory and Safety Critical Applications, chapter Monte Carlo optimisation for conflcit resolution in air traffic control, pages 257–276. LNCIS series. Springer, Berlin, 2006. [77] MFF. MFF Final safety case. Report MFF D734, ed. 1.0. Available at http://www.medff.it/public/index.asp, November 2005. 89 Bibliography [78] NASA. Concept definition for distributed air-/ground traffic management (DAGTM). Version 1.0, Advanced Air Transportation Technologies project, Aviation System Capacity Program, National Aeronautics and Space Administration, NASA, 1999. [79] B. Klein Obbink. MFF airborne self separation assurance OSED. Report MFF R733D. Available at http://www.medff.it/public/index.asp, April 2005. [80] R. A. Paielli and H. Erzberger. Conflict Probability Estimation For Free Flight. AIAA Journal of Guidance, Control and Dynamics, 20(3):588–596, 1997. [81] G. Pola, M. Bujorianu, J. Lygeros, and M. D. Di Benedetto. Stochastic hybrid models: an overview with applications to air traffic management. In Proceedings of IFAC Conf. Analysis and Design of Hybrid Systems (ADHS), Saint Malo, Brittany, France, 2003. [82] M. Prandini and J. Hu. Stochastic Hybrid Systems: Theory and Safety Critical Applications, chapter A stochastic approxmation method for reachability computations, pages 107–139. LNCIS series. Springer, Berlin, 2006. [83] M. Prandini and J. Hu. Stochastic reachability: Theory and numerical approximation. In C.G. Cassandras and J. Lygeros, editors, Stochastic Hybrid Systems: Recent Developments and Research Trends. CRC Press, 2006. [84] M. Prandini, J. Hu, and S. Lygeros, J. amd Sastry. A probabilistic approach to aircarft conflict detection. IEEE Tr. on Intelligent Transportation Systems, 1(4):199– 220, 2000. [85] M. Prandini, J. Lygeros, A. Nilim, and S. Sastry. A Probabilistic Framework for Aircraft Conflict Detection. In AIAA Guidance Navigation and Control Conference, 1999. [86] H. J. Rome and V. Krishnan. Causal probabilistic model for evaluating future transoceanic airline separation. In IEEE Transactions on Aerospace and Electronics Systems, volume 26, pages 804–817, September 1990. [87] RTCA. Free Flight Implementation. Task Force 3 Final Technical Report, Washington DC, 1995. [88] A. P. Shah, A. R. Pritchett, K. M. Feigh, S. A. Kalarev, A. Jadhav, K. M. Corker, D. M. Holl, and R. C. Bea. Analyzing air traffic management systems using agentbased modeling and simulation. In Proceedings of the 6th USA/Europe Seminar on Air Traffic Management Research and Development, Baltimore, MD, 27-30 June 2005. [89] R. Sheperd, R. Cassell, R. Thava, and D. Lee. A reduced aircraft separation risk assessment model. In Proceedings of the AIAA Guidance, Navigation and Control Conference, AIAA, Reston, VA August 1997. ” 90 Bibliography [90] S. H. Stroeve, H. A. P. Blom, and M. N. van der Park. Multi-agent situation awareness error evolution in accident risk modeling. In Proceedings of the 5th USA/Europe Seminar on Air Traffic Management Research and Development, Budapest, Hungary, 23-27 June 2003. [91] O. Strter, V. Dang, B. Kaufer, and A. Daniels. On the way to assess errors of commission. Reliability Engineering and System Safety, 83:129–138, 2004. 91