AIAA-99-4144 A PROBABILISTIC FRAMEWORK FOR AIRCRAFT CONFLICT DETECTION Maria Prandini, John Lygeros, Arnab Nilim and Shankar Sastry Department of Electrical Engineering and Computer Sciences University of California at Berkeley - Berkeley CA 94720 fprandini, lygeros, nilim,
[email protected] ABSTRACT We describe a general con ict detection/resolution scheme, focusing on the con ict detection component for a pair of aircraft ying at the same altitude. The proposed approach is formulated in a probabilistic framework, thus allowing uncertainty in the aircraft positions to be explicitly taken into account when detecting a potential con ict. The computational issues involved in the application of the proposed con ict alerting system are addressed by resorting to randomized algorithms. Finally, the validation of the proposed detection scheme is performed by Monte Carlo simulation on a stochastic ODE model of the aircraft motion. Our simulations show very promising results. The detection scheme will be used as the basis for a con ict resolution scheme in forthcoming work.
INTRODUCTION Safety, of which con ict prediction and resolution form an integral part, is the primary concern of all advanced air trac management systems. Con ict prediction and resolution are considered at three di erent levels in the air trac management process: Long range: Some form of con ict prediction and resolution is carried out at the level of the entire National Airspace System (NAS), over a horizon of several hours. It involves composing ight plans and airline schedules (on a daily basis, for example) to ensure that airport and sector capacities are not exceeded. This is typically accomplished by large scale integer and linear programming techniques1;2 . Research supported by DARPA under grant F33615-98-C-3614, by NASA under grant NAG 2-1039 and by ARO under grant MURI DAAH 04-96-1-0341. The authors would like to thank Dr. George J. Pappas for the stimulating discussions. Copyright c 1999 by the Regents of the University of California. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.
Conflict probability information
Prediction
Resolution Updated flight plan
Sampled estimates of aircraft positions
Radar
ATC Aircraft
Figure 1: Con ict prediction/resolution functions
Mid-range: Con ict prediction and resolution is car-
ried out by air trac control (ATC), over horizons of the order of tens of minutes. It involves modifying the predetermined ight plan on-line, to ensure adequate aircraft separation. Semi-automated tools have been developed to assist air trac controllers with these decisions3 . The algorithms and techniques described in References 4 and 5 can also operate at this level. Short range: Finally, con ict prediction and resolution is also carried out on board the aircraft, over horizons of seconds to minutes. This is typically considered as a last resort solution. The Trac Alert and Collision Avoidance System (TCAS) currently operating on all commercial aircraft is such a prediction/resolution algorithm6;7. The algorithms of References 4 and 5 also belong to this category. The work presented here concentrates primarily on midrange con ict prediction. It is assumed that the proposed algorithms will eventually be implemented as computational tools to provide assistance to ATC. A con ict prediction and resolution framework consists of a number of components (Figure 1). If one thinks of these components as arranged in a control feedback loop,
the ATC, aircraft, and radar would correspond to the
plant (the ATC and radar playing the role of actuators and sensors respectively), while the prediction and resolution components would correspond to the controller. Our
goal is to design the two \controller" modules and verify that the closed loop system possesses certain desirable properties (safety, ATC workload reduction, etc.). In the present paper, we restrict our attention to speci c suggestions for the prediction and validation schemes. However, in ongoing work, we are investigating a number of design alternatives for prediction, resolution, and validation. The algorithms proposed here (highlighted in the remainder of this section and detailed in the section on con ict detection) make extensive use of randomized algorithms for estimating integrals and carrying out optimizations. The advantage of randomized techniques is that they tend to be computationally more ecient. They also provide analytical bounds on the accuracy of the approximation involved, provided one makes appropriate design choices (for the cost functions for con ict resolution for example). CONFLICT DETECTION
We consider two aircraft moving on the same horizontal plane, each following its individual ight plan. The ight plan is assumed to consist of a sequence of way points on the plane and a sequence of speeds for moving between them. Based on this information and estimates of the current positions of the aircraft, the models developed in References 3 and 8 can be adapted to provide empirically motivated estimates of the probability distribution for the projected position of the two aircraft in the near future. One can then de ne the probability of con ict (PC) at a future time, t, as the probability that at time t the two aircraft will be within an unsafe distance from one another. Con ict detection consists of estimating PC and initiating some action when PC is high. The design choices that enter into con ict detection are the di erent ways of computing estimates of PC and the di erent ways of weighting the value of PC at various times. The parameters that one needs to set are the prediction horizon and the threshold for declaring a con ict. Here we use randomized integration to estimate the value of PC at each time instant and the maximum value of PC over the prediction horizon. Con icts are then declared whenever the latter estimate exceeds a certain threshold. Other approaches for obtaining a measure of the likelihood of a con ict over a nite horizon are proposed in References 3 and 4, based respectively on on-line, analytical approximation of the maximum value of PC over the prediction horizon and on lookup tables created by intensive, o -line Monte Carlo simulations. Systematic
guidelines for setting threshold values can be found in Reference 9. CONFLICT RESOLUTION
Once a con ict has been declared, con ict resolution consists of modifying the ight plans of the two aircraft to ensure that the probability of con ict falls below a speci ed threshold. We are currently studying a model predictive control approach to achieve this goal. The idea is to use the current radar measurements, the prediction module and an optimization algorithm to minimize an appropriate resolution cost function over all possible ight plans. ATC is noti ed if changes in the upcoming way points are imminent. The process is repeated every time a new radar measurement becomes available. The design choices that enter into con ict resolution are the di erent resolution cost functions and the di erent optimization techniques. The parameters one needs to set are how far in advance should ATC be noti ed of ight plan changes (partly determined by human capabilities) and the bounds on the allowable ight plans (partly determined by aircraft capability and ATC protocol). The optimization can again be carried out either by randomized algorithms or combinatorially, aided by appropriate heuristics. VALIDATION
The performance of a con ict prediction/resolution scheme is measured in terms of safety, impact on ATC workload, and eciency (impact on fuel consumption, deviations from schedule, etc.). Here we restrict our attention only to the issue of safety . Validation of the prediction scheme requires one to model the radar and the aircraft. Here we use a simple model for the radar (additive white noise) and we introduce a stochastic di erence equation to model the aircraft movement. Validation of the detection scheme is then carried out by Monte Carlo simulation. At a later stage we hope to be able to provide theoretical bounds on the safety of the proposed design - including the resolution component - and estimate its impact on ATC workload by means of human-in-the-loop simulations. The rest of the paper is organized in four sections. We rst introduce the probabilistic framework in which we develop our prediction scheme and perform the validation. We then introduce our con ict detection scheme and describe how we deal with the involved computational issues. Finally, Monte Carlo simulation results obtained through the proposed validation model are reported in the nal section. In fact safety is closely coupled to ATC workload, as the controllers may choose to disregard disruptive advisories.
PROBABILISTIC MODELS A situation of con ict occurs when two aircraft come within a certain minimum allowable distance from one another. The minimumhorizontal separation for en-route airspace is 5 nautical miles (nmi) and the minimum vertical separation is 2000 ft above the altitude of 29000 ft, 1000 ft below 29000 ft. Inside the Terminal Radar Approach Control (TRACON), the horizontal separation may be reduced to 3 nmi. Here we restrict our attention to the 2D case, i.e. the case when the two aircraft y at the same altitude. This is for ease of explanation, since the generalization to the 3D case is straightforward. The following probabilistic models for the aircraft motion are inspired by the work of Reference 3. PREDICTION MODEL
Consider two aircraft, denoted by 1 and 2, moving on the same horizontal plane. We assume that the ight plan of each aircraft i is known and it is given by a sequence of way points, fPjigj =0;:::;ni , Pji 2 R2, i = 1; 2 (where the rst way point P0i represents the current position of aircraft i) and a sequence of speeds fvji gj =1;:::;ni , vji 2 R+, i = 1; 2, where vji is the speed of aircraft i between Pji?1 and Pji. Let fTjigj =0;:::;ni denote the nominal time required for aircraft i to arrive at way point Pji. Then, fTji g can be recursively computed by Tji =
kPji ? Pji? k + T i ; j = 1; : : : ; n i j? vi 1
j
1
(1)
with initial condition T0i = 0, since P0i corresponds to the current position of aircraft i. According to References 3, 4, and 10, the predicted position of aircraft i is a ected by uncertainty due to wind, prediction errors, and tracking, navigation, and control errors. Consequently, we assume that the position, xi (t) := [xi1(t) xi2(t)]T , of aircraft i at time t can be modeled as a bivariate Gaussian random variable, with mean value given by the nominal position of the aircraft along its ight plan. To re ect the fact that our uncertainty about the position of the aircraft increases the further we try to project into the future, we choose the covariance matrix such that the along-track standard dei (t), and the cross-track standard deviation, viation, AT i CT (t), increase with time. More speci cally, following i (t) grows linReferences 3, 4 and 10, we assume that AT early with time: i (t) = c1 + c2 t; AT (2) i (t) grows linearly with the traveled distance, whereas CT
si (t), until it reaches a saturation value according to equation i (t) = minfc4 ; c1 + c3 si (t)g; CT
(3)
where si (t) can be computed by si (t) = vji (t ? Tji?1) + si (Tji?1), 8t 2 (Tji?1; Tji], j = 1; : : : ; ni, initialized with si (T0i ) = 0. The values of the initial tracking error standard deviations c1, their growth rates c2 and c3, and the nal cross-track error standard deviation c4 are set equal to c1 = 50=1850 nmi, c2 = 0:25 nmi=min, c3 = 1=57 and c4 = 1 nmi, which have been obtained empirically3 based on air trac data with a time horizon of length T = 20 min. In particular, the error covariance associated with the current position xi (0) represents the e ect of the noise a ecting the radar measurements. Overall, the distribution of xi(t) for i = 1; 2 is given by xi (t) N (pi (t); Qi (t)); where the mean value pi(t) is the nominal aircraft position on its ight plan, i.e., Pi ? Pi (4) pi (t) = Pji?1 + vji (t ? Tji?1 ) kPji ? Pji?1k ; j j ?1 t 2 [Tji?1; Tji ), for j = 1; : : : ; ni, and the covariance matrix Qi (t) is obtained through the expression
2;i Qi (t) = R(ji ) AT0 (t) 2;i0 (t) R(ji )T ; (5) CT t 2 [Tji?1; Tji), where the rotational matrix R(ji ) 2 SO(2) associated with the heading ji (ji being the angle that vector Pji ? Pji?1 makes with the x1 axis of the global coordinate frame in which the Pji are given) is
cos(ji ) ?sin(ji ) : R(ji ) = sin( (6) i i j ) cos(j ) Figure 2 gives an overall idea of the probabilistic model for the motion of each aircraft. The positions of the two aircraft, x1(t) and x2(t), are assumed to be uncorrelated. It is important to observe that even though this assumption is commonly made in the literature8 , it may be rather unrealistic. The tracking noise is in fact primarily due to wind, which may be strongly correlated between the two aircraft, especially near the con ict point where the aircraft are close one to the other. Modeling the tracking errors correlation properties is not dealt with in this paper, though it obviously needs to be evaluated in future work.
P1 1
Uncertainty at first way point (before the turn)
Uncertainty at first way point (after the turn)
θ v
1
1 2
Uncertainty at second way point
2
P
θ
P
1 0 v
1 2
where v is the aircraft linear velocity. The error vector := [ 1 2 ]T , on the other hand, is obtained through the stochastic di erential equation: _ = A + ; (0) N (0; V ) (9)
1 1
1 1
Initial uncertainty
Figure 2: Probabilistic model for the motion of aircraft 1 χ1
x2 v
χ2
θ
p x1
Figure 3: The body and inertial frames VALIDATION MODEL
As before, consider an aircraft moving on a horizontal plane. Let (x1; x2) denote its position with respect to an inertial coordinate frame. Assume its velocity has magnitude v and makes an angle with respect to the x1 axis (Figure 3). Consider a body coordinate frame ( 1; 2) with 1 aligned with the aircraft velocity (along-track) and 2 perpendicular to it (cross-track). The two frames are related through a coordinate transformation: x1 = R() 1 + p; (7) 2 x2 where R() is the rotation matrix given in equation (6) and p denotes the position of the body reference frame origin with respect to the inertial reference frame. In light of the above discussion, equation (7) can be reinterpreted as follows
p represents the nominal position of the aircraft; and represent the variation in the aircraft po1
according to the deterministic di erential equation p_ = R() 10 v; (8)
2
sition with respect to the nominal one in the alongtrack and cross-track directions respectively; x := [x1 x2 ]T represents the actual position of the aircraft.
The nominal position p of the aircraft evolves in time
where := [1 2]T is a white Gaussian noise process with (t) N (0; V (t)), independent of (0). We shall show below that by appropriately selecting matrix A and the covariance matrices V and V (t), we are able to obtain uncorrelated Gaussian along track and cross-track errors whose standard deviation respectively keeps growing with time and saturates according to equations (2) and (3). Using the above set of equations, we can derive a kinematical model of the aircraft trajectory in terms of the inertial frame coordinates. The aircraft position x evolves in time according to the following equation _ + R _ + p_ x_ = R by (7) _ = (R + RA) + R + p_ by (9) 1 T = (R_ + RA)R (x ? p) + R 0 v + R by (7) & (8) where we set R := R() to simplify the notation. By _ observing that R_ = dR d , we nally get the equations governing the aircraft motion 8 T ! + RART ]x ? [ dR RT ! + RART ]p > x_ = [ dR R > > < d d +R [0 1]T v + R > > p_ = R [0 1]T v > : _ = ! where (v; !) are the linear and angular velocities. Note that the equations describing the aircraft motion are nonlinear. The situation can be simpli ed somewhat by assuming that ! 0, and modeling changes in heading as discrete events, taking place at the way points of the
ight plan. This makes the aircraft dynamics piecewise linear, with switching times determined by the occurrence of the turn events Turn(), the e ect of Turn() being the addition of to the current value of (Figure 4). Our validation model is then represented by the stochastic di erential equation x_ = A(j )(x ? p) + B(j )vj + C(j ) ; (10) p_ = D(j )vj
Proposition 1 Consider the stochastic di erential equation:
where > 0? and () is a white Gaussian noise V ) , independent of with (t) N [0 0]T ; diag(V (t); ? (0) N [0 0]T ; diag(V ; V ) . Then, () satis es the following properties
and x(0) is independent of the white Gaussian noise determining the tracking errors. We assume that estimates of the position of the aircraft, y := [y1 y2 ]T , are made available to the ATC every seconds (typically 12 seconds) through radar measurements. These measurements are a ected by noise (12)
where the measurement noise f(k)gk0 is described as a sequence of independent identically distributed (i.i.d.) Gaussian random variables with (k) N (0; V ) and it is assumed to be independent of all the other random variables involved in the validation model. For our stochastic validation model to resemble the statics derived from air trac data, we need to appropriately tune its parameters, more speci cally, the covariance matrices V and V(t) , and the matrix A entering the matrix A() of the aircraft trajectory model (10), and the covariance matrix V for the radar model (12). In the following, we shall use diag() as shorthand for the diagonal matrix with diagonal entries listed in the brackets. For the radar model parameters, we set V = diag(V ; V ), where V = V = c21, to reproduce the statistical characteristics estimated in Reference 3 for the uncertainty in the current aircraft position. The tuning of the aircraft trajectory model is performed in light of the following proposition. 2
1
2
2
1
for t 2 [Tj ?1; Tj ), j 1, where A() = R()AR()T , B() = D() = R() [0 1]T , C() = R(), and the switching times fTj g are the nominal arrival times at the corresponding way points fPj g, computed according to equation (1). The initial conditions for (10) are: x(0) N (P0 ; R(1)V R(1 )T ) ; (11) p(0) = P0
TUNING OF THE VALIDATION MODEL
2
1
Figure 4: Piecewise linear model for the aircraft motion
y(k) = x(k) + (k);
_ = 00 ?0 + ; t 2 [0; T]
. x = A(θ)(x-p)+B( θ)u+C( θ) η . p = D(θ )u
1
_ θ := θ + θ
_ Turn(θ)
1. 1() and 2 () are independent Gaussian processes; 2. 1() has zero meanR and variance var[ 1(t)] = V + ot V (u)du; 3. 2() has zero mean and variance var[ 2(t)] = V2 [1 ? e?2 t] + V e?2 t. 1
1
2
2
Proof: The proof is based on standard results for the
solution of stochastic ordinary di erential equations and is omitted due to space limitations. For the noise () to reproduce the empirically observed characteristics of the along-track and cross-track errors, we set V = V = c21 ; V = 2c1c2 + 2c22 t; V2 = c24; (13) A = 00 ?0 ; with = (c c?c ) v1; which leads to var[ 1 (t)] = (c1 + c2 t)2 c (14) var[ 2 (t)] = c24 + [c21 ? c24]e? c ?c v t : 1
2
1
2
3
4
1
2 3 ( 4 1)
1
Thus, var[ 1(t)] is equal to the along-track variance 2 AT (t) in equation (2), whereas var[ 2(t)] is monotonically increasing from the initial value c21 to the sat2 uration value c24 as the cross-track variance CT (t) in equation (3). Moreover, the time constant of the decaying exponential in (14) is equal to 0.5 of the time 2 (t) in (3) to reach the satut = cc ?vc required for CT 2 ration value c4 in the case when the aircraft trajectory is a straight line traveled with speed v1 . For typical speeds (v1 = 500 nmi=h) such a time is 6 minutes, therefore by setting = (c c?c ) v1 we get that var[ 2 (t)] reaches 86% of the saturation value in about 6 minutes. As for the prediction model, for the validation model we again assume that the positions and radar measurement noise of di erent aircraft are uncorrelated. As discussed above this assumption may be unrealistic in practice. 4 1 3 1
3
4
1
CONFLICT DETECTION GENERAL CONFLICT DETECTION
Consider again two aircraft, 1 and 2, ying at the same altitude, and let x1 and x2 denote their positions in an inertial reference frame and d := x1 ? x2 denote their horizontal distance. Given that a con ict occurs when the aircraft are closer than 5 nmi to one another, the probability of con ict at time instant t is given by PC(t) = Pr(kd(t)k 5) =
Z
kuk5
f(u; t)du;
(15)
where f is the probability density function of the distance d at time t. Following the prediction model proposed in the previous section, x1 and x2 are independent, with xi (t) N (pi (t); Qi(t)), i = 1; 2. Therefore, d(t) is a Gaussian random variable with mean (t) := p1 (t) ? p2 (t) and variance Q(t) := Q1 (t) + Q2(t), hence T ? e? (d?(t)) Q(t) (d?(t)) (16) f(d; t) = p 1 2 det(Q(t)) In the con ict prediction algorithm proposed here, the detection of a situation of con ict in the near future is based on the maximum value assumed by the probability of con ict over a certain time horizon T. An \alert" is issued when the maximum probability of con ict exceeds a certain threshold, P. The process should is repeated every time the con guration, , of the two aircraft (consisting of their ight plans) changes, i.e. every time a new radar measurement becomes available (hence the P0i change), every time the ATC changes a ight plan, etc. The idea is that only when an alert is issued should a con ict resolution algorithm be executed to compute a modi cation of the two aircraft ight plans such that the probability of con ict is appropriately reduced. The value for the time horizon T is chosen to be equal to 20 minutes, following References 3, 8. As for the threshold on the probability of con ict, it can be set either by heuristic arguments such as those used in Reference 8, or by referring to the procedure suggested in Reference 9, where the main objective is to strike a compromise between the number of false alarms and missed detections. Other aspects that should be taken into account for the tuning of the threshold are the detection of con ict a certain amount of time before it occurs and the prediction of its occurrence time (as to avoid \last minute" maneuvers). These aspects highly in uence the e ectiveness of every prediction/resolution scheme involving the humanin-the-loop component. The detection scheme can be encoded in pseudo-code: 1 2
1
Algorithm 1 (Con ict Detection) when changes do begin end
Compute C( ) := tmax PC(t) 2[0;T ] if C( ) P declare a con ict
The major obstacle in the implementation of this prediction scheme is the computation of C( ), since PC() cannot be expressed in analytical form, and even the evaluation of PC(t) for a given t is time consuming. This obstacle turns out to be overwhelming when it comes to on-line application of the prediction scheme, where computation is subject to time constraints. A way around these complicated implementation issues is to abandon the ambitious objective of nding an exact solution, and compute instead an approximate solution, together with a methodology for providing quantitative information on the level of approximation introduced. Here this is accomplished using results from the theory of empirical processes, whose main aim is to study how to estimate unknown quantities through experimentation12. RANDOMIZED ESTIMATION OF C AND P C
Suppose for the time being that given a t 2 [0; T], we are able to compute PC(t) with no error. Let us introduce the uniform probability distribution Q on [0; T]. We can then resort to the following optimization algorithm to estimate the maximum of PC(t) over [0; T].
Algorithm 2 (Randomized Optimization) Initialization: Choose an integer N C 0( ) = 0 for i = 1; : : : ; N begin Extract at random ti 2 [0; T] according to Q if C 0( ) < PC(ti) then C 0( ) = PC(ti) end
Since we are testing just N values of t, we cannot expect that C 0 ( ) is the global maximum of PC(t). In addition, the quality of the estimate is random due to the stochastic selection of tj . Nevertheless, it can be shown that C 0( ) is a \good estimate" in some probabilistic sense.
Theorem 1 (Randomized Optimization) Let f : Z ! R be a measurable function on the probability space (Z; F; Q). Select N independent z N := (z1 ; : : : ; zN ) 2 Z according to Q and set f = maxz2fzj gNj f(z). Fix an arbitrary real number > 0. Then, Qfz 2 Z : f(z) > fg ; (17) =1
with probability greater than 1 ? (1 ? )N .
P^ (A) := M1
Proof: See Lemma 11.1 in Reference 12.
In the previous statement, f is a random variable de ned on the product space Z N with the product probability measure QN hosting the random extraction (z1 ; : : : ; zN ). Thus, Qfz 2 Z : f(z) > fg may or may not be satis ed depending on the random extraction and, so, Qfz 2 Z : f(z) > fg de nes a probabilistic event in the space Z N . Theorem 1 says that the probability QN of this event is at least 1 ? (1 ? )N . A possible interpretation of equation (17) is the following Suppose that an opponent would like to determine a z such that f(z ) > f by randomly selecting a z in Z according to Q. Then, her probability of \beating" f is not greater than . In Theorem 1, parameter is arbitrary and, therefore, the probability of success left to the opponent can be reduced at will. On the other hand, (17) is not a deterministic statement and it only holds true with probability at least 1 ? (1 ? )N . As approaches 0, 1 ? (1 ? )N tends to 0 and so the statement becomes evanescent. In addition, statement (17) only quanti es the probability that one can improve the estimate f by randomly selecting a new parameter z, but says nothing about how large such an improvement can be. Here, we are interested in maximizing the probability of con ict PC(t) over the time interval [0; T]. Observe rst that, given a con guration, , the corresponding probability of con ict PC() is a piecewise smooth function; it only loses smoothness on points in the nite set fTj 2 T : Tj T g, where T = fTj1gj =0;:::;n [ fTj2gj =0;:::;n [ fT1; T2g, Ti are the time instants when the cross-track error variances saturate (c1 + c3 si (Ti ) = c4 ) and fTji gj =0;:::;ni are the nominal arrival times at the way points. This implies that if Z is taken to be the interval [0; T], F the Borel -algebra on [0; T] and Q the uniform probability distribution, Theorem 1 can be applied with function f(z) given by the probability of con ict PC(z). Next, we introduce a method to approximately compute function PC(t), which allows us to circumvent the dif culty of numerically integrating the probability density function f(d; t) in (16) over the circle C = fz 2 R2 : kz k 5g. Given the measurable space (Z; F) and a probability measure P on (Z; F), suppose that we want toRcompute the measure of a set A 2 F, given by P (A) := A P (dz): An estimate, P^ (A), of this quantity can be constructed by picking some points z M := (z1 ; : : : ; zM ) at random from Z according to the distribution P and setting 1
2
M X
i=1
IA (zi );
(18)
where IA (z) is the indicator function of A, that is IA (z) = 1 if z 2 A, 0 otherwise. This makes the problem computationally tractable, and also allows us to resort to results from the theory of empirical processes to quantify the level of approximation introduced when we substitute the integral by a nite sum.
Theorem 2 (Estimation of Probability Measures) Fix an accuracy parameter 2 (0; 1). Suppose that (Z; F) is a measurable space and P is a probability measure on (Z; F). Fix a set A 2 F . Extract z M i.i.d. samples drawn from Z in accordance with P and compute P^ (A) according to (18). Then,
P M fz M 2 Z M : jP^ (A) ? P (A)j > g 2e? M : Proof: By the Cherno bound. 2
2
Suppose now that we are given, not just a single set A, but a collection of sets A F. By drawing M i.i.d. samples z M from Z in accordance with P , we can form empirical estimates for each P (A), A 2 A as before, i.e. P^ (A) := M1 PMi=1 IA (zi ); 8A 2 A: If the set A has nite cardinality jAj < 1, from Theorem 2 we get that P Mfz M 2 Z M : sup jP^ (A) ? P (A)j > g 2jAje?2M (19) 2
A2A
which means that the collection of sets A has the property of uniform convergence of empirical probabilities (UCEP), since for each xed the empirical probabilities converge uniformly to their true values as the number of samples M goes to in nity. To be able to apply (19) (and also to improve the conditioning of the problem) we can transform the problem of computing PC(t) to that of computing the probability measure of a time dependent set according to a xed probability distribution. Using a change of variables, the measure of the xed set C = fz 2 R2 : kz k 5g according to the probability distribution N ((t); Q(t)), can be viewed as the measure of an appropriate time dependent set Et according to the xed (standard) probability distribution N (0; I), By computing the Cholesky factorizationy of the covariance matrix Q(t) = L(t)L(t)T and setting v = L(t)?1 [u ? (t)], we get Z 1 e? vT v dv; PC(t) = 2 Et 1 2
yA
similar procedure is followed in Reference 3.
where Et = fv 2 R2 : L(t)v + (t) 2 Cg. The estimate PC 0(t) can therefore be calculated by extracting M i.i.d. samples vM from P N (0; I) and computing M X IEt (vi ); PC 0(t) := M1 i=1
(20)
Theorem 3 (Estimation of C ) Given 2 (0; 1), 2 (0; 1) and 2 (0; 1), C 0( ) is an approximate estimate of C( ) to accuracy 2 and level with con dence 1 ? in the sense that
RANDOMIZED CONFLICT DETECTION ALGORITHM
Qft 2 [0 T] : PC(t) > C 0 ( ) + 2g ; with probability greater than 1 ? , where Q denotes the
Algorithm 3 (Randomized Con ict Detection) Initialization: Fix P 2 [0; 1] Fix 2 (0;l 1), 2m(0; 1), 2 (0; 1)
[0; T]N R2M . Since the extractions of all the ti and all the vj are independent of one another the probability measure in [0; T]N R2M is just the product probability measure QN P M , where Q and P respectively denote the uniform probability distribution on [0; T] and the standard Gaussian probability distribution on R2. Now, according to Theorem 1 relation
A joint use of the randomized optimization algorithm with the randomized estimate of PC(t) leads to a fully randomized algorithm for con ict detection. Let dz e denote the smallest integer greater than a real number z.
ln( ) 1 4N Set N = ln(1 ? ) , M = 2 ln when changes do 2
2
uniform probability distribution on [0; T]. Proof: C 0 is a random variable on the probability space
Qft 2 [0 T] : PC(t) > t2fmax PC(t)g ; t gN
begin
Extract at random vj 's 2 R2, j = 1; : : : ; M according to N (0; I) for i = 1; : : : ; N do
begin
Extract at random ti 2 [0; T] according to the uniform density function Compute (ti ) = p1 (ti) ? p2 (ti ) by (4) Compute Q(ti) = Q1(ti ) + Q2 (ti) by (5) Compute the Cholesky decomposition Q(ti) = L(ti )L(ti )T PC 0(ti ) = 0 for j = 1; : : : ; M do
begin if L(ti )vj + (ti ) 2 C then PC 0(ti) = PC 0(ti ) + 1 end 0
end
end
PC 0(ti ) = PCM(ti ) if PC 0(ti ) P declare a con ict
i i=1
holds true with probability QN greater than or equal to 1 ? (1 ? )N . Set q := P M fvM 2 R2M: supN jPC 0(t) ? PC(t)j > g: (22) t2ftigi=1
Putting together these two results we conclude that the following holds with a product probability QN P M not less than 1 ? [(1 ? )N + q]
Qft 2 [0; T] : PC(t) > PC(arg t2fmax PC 0(t)) + 2g tigNi Qft 2 [0; T] : PC(t) > PC(arg t2fmax PC(t))g : t gN =1
i i=1
where the rst inequality follows from (22) and the second from (21). Thus, maximizing PC 0(t) over fti gNi=1 leads to an approximate maximum of PC(t) to accuracy 2 and level with con dence 1 ? where = (1 ? )N + q. If we now use estimate q 2Ne?2M in equation (19), we ln( ) e and M= can easily conclude that setting N = d ln(1 ? ) 4N 1 d 2 ln e suces to approximately maximize PC(t) to accuracy 2 and level with con dence 1 ? . It is interesting to observe that the number of samples needed to achieve a certain approximation level in terms of accuracy , con dence and level , is independent of the nature of the sample space and of the probability distribution. Thus, by extending the proposed approach to the 3D case, the computational load does not significantly increase. This will not be the case if one would resort to numerical methods based on gridding. 2
2
Note that the algorithm declares a con ict if and only if the estimate of the measure of criticality M X C 0 ( ) := maxN M1 IEti (vj ); t2ftigi j =1 exceeds the threshold, P. It is assumed that all the random extractions are made independently of one another. The following theorem provides an estimate of the accuracy of our approximation (see also Reference 13). =1
(21)
2
VALIDATION
1. Given the ight plans of the two aircraft, generate pairs of aircraft trajectories over a 20 minutes time horizon according to a discretized version of the stochastic ODE (10) initialized with (11) (sample time = 1 second), and the corresponding radar measurements according to equation (12) ( = 12); 2. For each pair of simulated trajectories, execute the detection algorithm at every radar measurement time, each time using the updated ight plans and time horizon. These are respectively obtained by removing the way points which have been surpassed, setting the rst way point equal to the current radar measurement, and subtracting the elapsed time from the 20 minutes initial horizon. A con ict is declared as soon as the estimated value of the maximum probability of con ict exceeds the prescribed threshold; 3. Compute the probability of false alarm P(FA) and the probability of successful alert P(SA), i.e., the ratio between the number of alerts issued when there was no con ict and the total number of cases when there was no con ict (P(FA)), and the ratio between the number of alerts issued at least 60 seconds before a situation of con ict e ectively happens and the total number of con icts (P(SA)). It is evident that the estimates of P(FA) and P(SA) should both be decreasing functions of the threshold. The System Operating Characteristic (SOC) curve9 is a plot of P(SA) versus P(FA) parameterized by the threshold (see for example Figures 6 and 9). An ideal con ict detection scheme should operate near the point (0; 1), where there are no false alarms and all the situation of con icts are detected. On the other hand, a real con ict detection scheme cannot operate at this point - irrespective of the threshold chosen - due to the uncertainty a ecting the aircraft positions. However, the more the SOC curve approaches the point (0; 1), the better the performance of the system is likely to be. The threshold can, therefore, be chosen to ensure operation as close as possible to the ideal point, in an attempt to strike an \optimal" compromise between the number of false alarms and the number of successful alerts.
100
80
60
40 x2
The proposed con ict detection scheme was validated using Monte Carlo simulations of the stochastic ODE model of the aircraft trajectories and the model of the radar noise. The protocol for evaluating the detection scheme performance consists of the following steps:
120
20
0
−20
t=0
−40
−60
−80 −20
t=0
0
20
40
60
80 x1
100
120
140
160
180
Figure 5: Sample pair of simulated trajectories for Ex. 1. The stand for way points.
In the remainder of this section we describe the validation results for two cases: straight ight plans crossing at 90 and zig-zag ight plans. For both cases, we draw the SOC curve and compare the performance obtained by the detection algorithm proposed here and the detection algorithm proposed by Erzberger et al.8. The latter algorithm is based on the same description of the along-track and cross-track errors as our model for con ict prediction. Unlike our approach, however, the estimate of the probability of con ict on which detection is based, is the one corresponding to the minimum deterministic separation (computed based on the nominal ight plans). This estimate of the probability of con ict is computed by an analytical approximation of the integral of the distance probability density function obtained by extending the circular domain of integration to a whole strip (see Reference 3 for more details). Example 1 (90 path crossing angle con guration) Consider the case where the two aircraft are ying at the same altitude, along straight paths crossing at an angle of 90, and at speeds of 480 nmi/h and 500 nmi/h. For the sake of clarity, in Figure 5 we have drawn a realization of the trajectories obtained by performing step 1 of the described veri cation protocol. By running our detection algorithm and the one proposed in Reference 8 on 1000 Monte Carlo samples, we estimate the probability of successful alert and false alarm and draw the SOC curves. The results are shown in Figure 6 (the solid curve corresponds to the detection algorithm proposed in this paper while the dotted curve to the algorithm of Reference 8. We can compute the \optimal" threshold value P as the one corresponding to the point of the SOC curve nearest to ideal operating point (0; 1). This leads to P = 0:78 (P(FA)=0.235, P(SA)=0.737) and P = 0:82 (P(FA)=0.716, P(SA)=0.765) respectively for our detection procedure and for the one proposed in Reference 8. Thus, our algorithm seems to have a slightly lower proba-
1
1
0.95 0.9 0.9
0.85 0.8
P(SA)
P(SA)
0.8
0.75
0.7
0.7 0.6 0.65
0.6 0.5 0.55
0
0.1
0.2
0.3
0.4
0.5 P(FA)
0.6
0.7
0.8
0.9
0.4
1
Figure 6: SOC curves for our algorithm (solid) & Erzberger
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.1
0.2
0.3
0.4
0.5 P(FA)
0.6
0.7
0.8
0.9
1
Figure 9: SOC curves for our algorithm (solid) & Erzberger algorithm (dotted). The stand for the \optimal" threshold points (Ex. 2: = 0:05; = 0:1; = 0:05; 1000 simulations).
P(FA)
P(SA)
P(FA)
algorithm (dotted). The stand for the \optimal" threshold points (Ex. 1: = 0:05; = 0:1; = 0:05; 1000 simulations).
0
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6 P(SA)
0.5
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.2
0.1
0.1
0
0
0.1 0
0.1
0.2
0.3
0.4
0.5 threshold
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5 threshold
0.6
0.7
0.8
0.9
1
0
0.1
0
0.1
0.2
0.3
0.4
0.5 threshold
0.6
0.7
0.8
0.9
1
0
0
0.1
0.2
0.3
0.4
0.5 threshold
0.6
0.7
0.8
0.9
1
Figure 7: Plots of P(FA) and P(SA) vs threshold (Ex. 1).
Figure 10: Plots of P(FA) and P(SA) vs threshold (Ex. 2).
bility of successful detection and a false alarm probability of 23:5%, which is 1/3 that of Reference 8. 2 Example 2 (zig-zag ight paths con guration) Finally, we consider the case when the sequence of way points in the ight plans describes a zig-zag con guration at a xed altitude with the speeds being v11 = 505 nmi/h and v21 = 455 nmi/h, v12 = 480 nmi/h and v22 = 470 nmi/h. Figure 8 illustrates a realization of the trajectories. The estimates of the probability of successful alert and false alarm obtained by 1000 Monte Carlo samples were used to plot the SOC curves in Figure 9 (again the
solid line SOC curve corresponds to the detection algorithm proposed in this paper). In this case, the \optimal" threshold values for detection detection procedure proposed here and the one proposed in Reference 8 are respectively P = 0:92, corresponding to P(FA)=0.319 and P(SA)=0.686, and P = 0:97, corresponding to P(FA)=0.494 and P(SA)=0.543. Thus, our algorithm has a higher probability of successful detection and a false alarm probability which is about 3/5 that of Reference 8. 2 In these two examples, the SOC curve obtained by the algorithm of Reference 8 is lower than the one obtained by the algorithm proposed in this paper. In our opinion, the reason for this rests on the fact that the parameter used for detecting a situation of con ict in Reference 8 is an over approximation of the probability of con ict at the time of minimum deterministic distance (which is assumed to be the most critical time instant over the considered time horizon). This is re ected in the higher value assumed by the probability of false alarm in the algorithm of Reference 8 with respect to our approach, irrespective of the chosen threshold. On the other hand, this increase of P(FA) is not compensated by an adequate increase of the probability of successful alert as it is shown in Figure 7 and 10 representing P(FA) and P(SA) as a function of
140 120 100 80 60 40 20 0 −20 −40 −60 −80 −15
−10
−5
0
5
10
15
20
25
30
Figure 8: Sample pair of simulated trajectories for Ex. 2. The stand for way points.
the threshold for the two examples, with the solid line corresponding to our detection scheme and the dotted line corresponding to the detection scheme in Reference 8. It should however be noted that the detection algorithm of Reference 8 is computationally much simpler, and therefore may be better for real time implementation. The computational performance of our algorithm (but also its accuracy) can be a ected by choosing the tolerance parameters ; and . It is important to observe that, as expected, di erent
ight plans lead to di erent SOC curves. A sensitivity analysis of the dependence of the threshold on the ight plans should be performed by parameterizing them as a function for example of the crossing angles, minimum deterministic distance, and time of minimum distance, in order to set a value for the threshold which is appropriate for the typically encountered path con gurations. This issue is currently under investigation.
CONCLUSIONS We outlined a framework for con ict prediction and resolution, and speci cally dealt with the prediction component for pairs of aircraft moving on the same horizontal plane. We are currently working on formulating algorithms for resolution and investigating the di erent design alternatives listed in the introduction. Further objectives are extending our approach to the three dimensional case and multiple aircraft. The extension to three dimensions is straightforward, while the correct way of extending this approach to multiple aircraft is less clear. It may turn out to be adequate to resolve con icts pairwise (as proposed in Reference 3), but this approach has the potential for combinatorial explosion and it is unclear whether it will always lead to an overall resolution. Once the prediction/resolution algorithm has stabilized, we hope to be able to test it in human-in-the-loop simulations. This will allow us to tune the various parameters, and assess its impact on ATC workload. In parallel we are working towards a methodology for formally evaluating the safety properties of the proposed algorithm. This will hopefully lead to a more general probabilistic veri cation methodology for hybrid systems.
REFERENCES 1. D. Bertsimas and S. Stock Patterson. The air trac
ow management problem with enroute capacities. Operations Research, 46:406{422, 1998. 2. P.B.M Vranas, D. Bertsimas, and A.R. Odoni. Dynamic ground-holding policies for a network of airports. Transportation Science, 28:275{291, 1994.
3. R.A. Paielli and H. Erzberger. Con ict probability and estimation for free ight. In Proc. 35th Meeting of the American Institute of Aeronautics and Astronautics, AIAA-97-0001, 1997.
4. L. Yang and J. Kuchar. Prototype con ict alerting logic for free ight. In Proc. 35th AIAA Airspace Science Meeting & Exhibit, AIAA-97-0220, January 1997. 5. C. Tomlin, G.J. Pappas, and S. Sastry. Con ict resolution for air trac management: a study in multiagent hybrid systems. IEEE Trans. on Automatic Control, 43:509{521, 1998. 6. Radio Technical Commission for Aeronautics. Minimum operational performance standards for traf c alert and collision avoidance system (TCAS) airborn equipment. Technical Report RTCA/DO-185, RTCA, September 1990. Consolidated Edition. 7. John Lygeros and Nancy Lynch. On the formal veri cation of the TCAS con ict resolution algorithms. In Proc. 36th IEEE Conf. on Decision and Control, San Diego, pages 1829{1834, 1997. 8. H. Erzberger, R.A. Paielli, D.R. Isaacson, and M.M. Eshow. Con ict detection and resolution in the presence of prediction error. In 1st USA/Europe Air Trac Management R & D Seminar, 1997. 9. James K. Kuchar. A Uni ed Methodology for the Evaluation of Hazard Alerting Systems. PhD thesis, Massachusetts Institute of Technology, 1995. 10. R.A. Paielli and H. Erzberger. Con ict probability estimation for free ight. Journal of Guidance, Control and Dynamics, 20:588{596, 1997. 11. G. Chen, G. Chen, and S. Hsu. Linear stochastic control systems. CRC Press, Inc., 1995. 12. M. Vidyasagar. A theory of learning and generalization: with applications to neural networks and control systems. Springer-Verlag, London, 1997.
13. M. Vidyasagar. Statistical learning theory and its applications to randomized algorithms for robust controller synthesis. In Plenary Lectures and MiniCourses at the European Control Conf., pages 161{ 189, 1997.
