- Email: [email protected]
The optimal control of collision avoidance trajectories in air traffic management
TRANSPORTATION RESEARCH PART B
Transportation Research Part B 33 (1999) 265±280
The optimal control of collision avoidance trajectories in air trac management John C. Clements * Department of Mathematics, Statistics, and Computing Science, Dalhousie University, Halifax, Nova Scotia, B3H 3J5, Canada Received 8 August 1996; received in revised form 21 June 1998
Abstract A collision avoidance/con¯ict resolution model for aircraft traversing planar intersecting trajectories is developed based on the actual ¯ight dynamics of each aircraft and the criteria that a ¯ight regime correction is optimal if it represents a minimum-time safe deviation from a preassigned ¯ight program. The formulation of the model leads to a singular optimal control problem where the control variable constraint is the maximum acceptable aircraft turn rate for passenger comfort and safety and the state variable constraint is de®ned in terms of the radius of the protected zone about a potentially con¯icting aircraft. A robust, accurate real-time solution procedure is derived for computing the steering program for the execution of the optimal safe avoidance maneuver. Several representative examples of pairwise trajectory con¯icts are analysed in detail. A brief discussion of how the model might be generalized to incorporate changes in airspeed and altitude as well as more comprehensive dynamics is included. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Air trac management; Optimal control
1. Introduction Although current systems are already highly sophisticated, air trac control in the future will become even more complex and will he based on ``user-preferred trajectories'' and the continuous automatic monitoring and control of ¯ight trajectories both from within the aircraft and via the relevant air trac control centers (Bowers, 1995; Cotton, 1995; Friedman, 1989; Massoglia et al., 1989; Niedringhaus, 1991; NASPG, 1995; Lanzer and Jenny, 1995; Rigotti and Dean, 1992; Tomlin *
Tel.: +1-902-494-3381; fax: +1-902-494-5130; e-mail: [email protected].
0191-2615/98/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0191-2615(9 8 ) 0 0 0 3 1 - 9
266
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
et al., 1997, 1998; Williamson and Spencer, 1989). Identifying potential con¯icts and determining adjustments to the aected ¯ight trajectories which are ``optimal'' in some sense will be an important component of such systems. A number of collision avoidance models have been examined, notably the linear heading (Friedman, 1991) and speed control (Friedman, (1988) avoidance strategies developed by Friedman, the linear programing approach taken by Niedringhaus (1991) for the resolution of multiple interrelated separation problems and, more recently, the presentation of a con¯ict resolution architecture for multi-agent hybrid systems by Tomlin et al. (1998). In this work, we present an optimal control con¯ict resolution model for aircraft traversing planar intersecting trajectories which diers from the studies listed above in that it incorporates the ¯ight dynamics of each aircraft and the resulting avoidance maneuvers are smooth trajectories de®ned by the turn rate constraints of the avoiding aircraft. The model is based on the criteria that a ¯ight regime correction is optimal if it represents a minimum~time safe deviation from a preassigned ¯ight program. The formulation of the model leads to a singular optimal control problem where the control variable constraint (3) is the maximum acceptable aircraft turn rate for passenger comfort and safety under normal operating conditions and the non-autonomous state variable constraint (4) is de®ned in terms of the radius rp of the protective zone about a potentially con¯icting aircraft. In Section 2, we con®ne our attention to the simplest model which involves two aircraft operating at constant altitude and constant speci®ed airspeeds and for which only heading changes of one aircraft are permitted. Without loss of generality, we can assume that aircraft #1 is travelling at airspeed 1, in the positive x-direction while aircraft #2 must execute an avoidance maneuver to a speci®ed designated objective (xf , yf ) at airspeed 2 (cf. Fig. 1). Here (xf , yf ) could be the original ¯ight plan destination of aircraft #2 or merely a point chosen on the intended ¯ight path of that aircraft which is beyond the region of potential con¯ict. The positions of aircraft #1 and #2 in the xy coordinate plane at time t 0 are taken to be (0, 0) and (x0 ; y0), respectively, and 0 is the initial heading (and ground track) angle of aircraft #2. The hatched circle of radius rp in Fig. 1 is the protective zone about aircraft #1. The optimal con¯ict avoidance trajectory is the minimum-time path C : f x t; y t j 04t4tf ; x ; y 2 C1 0; tf g from (x0 ; y0 ) to (xf ; yf ) which does not violate the (moving) protective zone x ÿ 1 t2 y2 r2p about aircraft #1. A robust, accurate real-time solution procedure based on an analysis of the related Hamiltonian and costate equations for the problem is developed in Sections 3 and 4. The calculation of the interval of singular control and the determination of the precise boundary arc control law (i.e. when the state variable constraint is strictly satis®ed) are not at all straightforward and are the main focus of this work. In particular, it is shown that determination of the boundary arc control requires the solution of the coupled dierential±algebraic system of equations given by (12), (14) andÿ (18). of the optimal trajectories concerns the point An ÿ interesting and important ÿ property t2 at which they intercept the boundary of the moving (x t2 , y t2 ) and heading angle protective zone. This is best understood by considering the problem transformed to a coordinate system centred on aircraft #1 using x ÿ 1 t and y and is described in detail in Section 4. The numerical integration of system (12) and (18) is straightforward and is not a serious computational consideration. The model studied here has the advantage that solutions can he continuously recomputed to allow for corrections in cross-track error and positioning accuracy throughout the con¯ict resolution maneuver. Several representative examples are considered in Section 5. Examples 1 and 2 are concerned with the determination of the optimal avoidance
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
267
trajectory when 2 > 1 . Example 3 examines a con¯ict resolution con®guration where 2 < 1 . In the ®rst two examples, the boundary arc represents the envelope of feasible minimum-time trajectories and is plotted beyond the boundary arc turn to clearly indicate the general trajectory de®ned by the state variable constraint (cf. the hatched line in Figs. 3 and 4). The generalization of the model to incorporate possible changes in airspeed and altitude as well as more comprehensive dynamics is discussed brie¯y in Section 6. 2. Formulation of the simplest model Let Ck 0; tf and Cp 0; tf denote the linear spaces of k times continuously dierential functions and piecewise continuous functions, on [0; tf ], respectively. The ¯ight trajectory of aircraft #2
Fig. 1. Potential con¯ict con®guration in the xy coordinate plane.
Fig. 2. Relationship between t, t and t on a boundary arc.
268
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
Fig. 3. Feasible extremal trajectories for example 1.
Fig. 4. Feasible extremal trajectories for example 2.
from (x0 ; y0 ) to (xf ; yf ) is given by C : x t; y t j 04t4tf where the aircraft position coordinates x t and y t are elements of C1 0; tf and tf is the time required to reach (xf , yf ). Neglecting the inertial acceleration terms associated with changes in aircraft heading, the dynamics of the problem can be formulated in terms of the system of equations x_ t 2 cos t; with end conditions
y_ t 2 sin t;
_ t u t;
04t4tf
1
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
x 0 x0 ; y 0 y0 ;
0
0;
ÿ ÿ x tf xf ; y tf yf
269
2
where the heading or steering angle t 2 C 0; tf is given in radians measured positive counterclockwise from the x-axis. u t 2 Cp 0; tf , 04t4tf , is identi®ed as the control variable and the control constraint j u t j 4 M ;
04t4tf
3
is preset as the maximum acceptable turn rate under normal operating conditions. u t M corresponds to a maximum allowable rate of turn to the left and t ÿ M to a maximum turn rate to the right. The non-autonomous state variable constraint is given by S t; x t; y t r2p ÿ x t ÿ 1 t2 y2 t40;
04t4tf
4
and corresponds to the protective zone of radius rp about aircraft #1. Here rp might be assigned a dierent value for aircraft #1 in each pairwise con¯ict depending on the speed and maneuverability of that aircraft or it may be set to a ®xed value for all aircraft. Thus, the simplest optimal control problem is to determine the steering program t, 04t4tf , which guides aircraft #2 along that trajectory C which reaches the (xf , yf ) in minimum time tf without violating the protected zone of aircraft #1. That is, the problem is to minimize the ®nal time functional tf min J u dt
u2Cp 0;tf
5
0
subject to the state Eq. (1), the end conditions (2) and the control and state variable constraints (3) and (4). For computational simplicity, it is assumed here that the initial heading 0 of aircraft #2 represents the shortest ¯ight path to (xf , yf ). Thus, if the state variable constraint (4) is not violated, the optimal control is given by the default value u t 0; 04t4tf . If the initial heading 0 of aircraft #2 is maintained for a given con®guration, then the minimum separation distance dm between the two aircraft occurs at time tm x0 1 ÿ 2 cos 0 ÿ y0 2 sin 0 = 21 ÿ 2 2 1 cos 0 22 6 and is given by q dm x0 tm 2 cos 0 ÿ 1 tm 2 y0 tm u2 sin 0 2 :
7
Hence, the defaultÿ control u1 t 0; 04t4tf is optimal if and only if dm 5rp , in which case 2 tf xf ÿ x0 2 yf ÿ y0 2 =u2 . If dm < rp , then a feasible extremal is a minimum-time trajectory which satis®es the constraints of the problem. Regarding the initial aircraft con®guration, it is further assumed that j yf j> rp and that the initial aircraft separation is at least d0 rp 1 uM . As
270
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
shown in Section 6, this ensures (except in the case u t 0; 04t4tf ) there always exists at least one and possibly two feasible extremal trajectories for the optimal control problem. The ®rst of these extremal trajectories passes aft of the advancing protective zone x ÿ 1 t2 y2 4r2p about aircraft #1 while the second, if it exists, will pass in front of this advancing zone. Both extremal trajectories must be geodesics. The optimal trajectory is that feasible extremal requiring the least time. If 2 4 1 , then only the trajectory passing aft of the protective zone about aircraft #1 is feasible and hence optimal. If 2 > 1 , then both feasible extremals exist. 3. The optimal control u t; 04t4tf The formulation and analysis of the Hamiltonian for system (1)-(5) is similar to that found in except for the derivation of the boundary control law required to accommodate the time-dependent state variable constraint (4). We show here that this law is given by the coupled dierentialalgebraic system (12), (14) and (18). The Hamiltonian for (1)-(5) is given by H t 1 l1 tv2 cos t l2 tv2 sin t l3 tu t tS t; x t; y t
8
for Lagrange multipliers li t 2 C1 0; tf ; i 1; :::; 3 satisfying the adjoint equations l_ 1 t ÿ@H=@x 0 l_ 2 t ÿ@H=@y 0
9
l_ 3 t ÿ@H=@ l1 t 2 sin t ÿ l2 t 2 cos t for all t in [0; tf ] where the multiplier t50 satis®es ÿ ÿ tS t; x t; y t 0 on [0; tf ]. From the 0 and l t 0 , l transversality conditions, l 3 f i i tf , i 1; 2 and l3 0 are free parameters. A arc subarc of : x t; y t; t j 04t4tf for which S t; x t; y t < 0 is called an interior and a boundary arc is a subarc of II for which S t; x t; y t 0. t in (8) is in C 0; tf and is continuously dierentiable on the interior of any boundary arc (Clements, 1993; Maurer, 1977). This is a singular control problem with switching function t @H=@u l3 t. Since l_ 1 t l_ 2 t 0 in (9), l1 and l2 are constant for all t in [0; tf ]. Consequently, on any interior arc of , the optimal control law is u t ÿ M u t 0 u t M
l3 t > 0 l3 t 0 l3 t < 0:
10
and a feasible extremal must begin with a maximum rate turn (u t uM ; 1; 04t4t1 ) for radial (u t 0; t1 ? 4t4t2 ) which intercepts the moving zone at a time t1 toÿ that ÿ ÿ protective ÿ point (x t2 ; y t2 ). The details regarding the calculation of t2 and (x t2 , y t2 ) are given in Section 4 and show in particular that the radial (u t 0; t1 ? 4t4t2 ) is not tangent to the advancing circular zone at the ÿ point ÿ of interception. Here (t1 ?; t2 ) is an interval of singular control. The maneuver to (x t2 ; y t2 ) is followed by an opposite turn along the boundary of
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
271
the moving circle x ÿ 1 2 y2 t r2p to the ®nal radial to the ¯y-to-point (xf , yf ) [cf. Fig. 3(a)] unless the maximum allowable turn rate is exceeded at the point of interception, in which case a maximum rate turn (u t ÿuM ; 1; t2 4t4t3 ) to the ®nal radial to (xf ; yf ) is executed [cf. Fig. 4(a)]. If an interval (t2 ; t3 ) of boundary arc control does exist, then > 0 on this interval and the optimal control law can no longer be determined from the form of the Hamiltonian (8) and the adjoint equations (9). Determination of the precise form of the control law on any boundary arc of is fairly involved but essential for safety considerations, particularly in the case of a trajectory passing in front of the protective zone. Let (x t; y t) be the trajectory components of aircraft #2 on a boundary arc and let t denote the heading and T t the component of the speed of aircraft #2 in the direction tangent to the moving zone at (x t; y t) as in Fig. 2. Then the geometry gives 2 sin t T t sin t 2 cos t T t cos t 1
11
2 sin t ÿ t 1 sin t and it follows that t t ÿ sinÿ1 1 sin t= 2
12
and 2T t 2 cos t ÿ 1 2 2 sin t2 :
13
Dierentiation of (12) gives the form of the required control on the boundary arc 2
3
1 cos t= 2 6 7 u t _ t 41 ÿ q5 _ t 1 ÿ 1 sin t= 2 2
14
where it remains to solve for t. Let t denote the angle which (x t; y t) makes with the origin in Fig. 2. Since (x t; y t) must remain on the boundary of the moving circle, ÿ x t rp cos t 1 t; rp cos t =2 1 t; y t rp sin t rp sin t =2: 15 where the corresponds to a right (ÿ) or left turn (+) on the boundary trajectory. It then follows from (13) and (15) that 2T t r2p _ t2 22 21 ÿ 2 2 1 cos t which yields after some simpli®cation
16
272
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
ÿ r2p _ t2 22 21 ÿ 2 2 1 cos t ÿ sinÿ1 1 sin t= 2 q2 1 cos t ÿ 22 ÿ 21 21 cos2 t :
17
Thus, the boundary arc control is given by (14) where t is computed as the solution of the initial value problem q
_ t ÿ 1 cos t ÿ 22 ÿ 21 21 cos2 t =rp ; t > t2 18 ÿ
t2 2 ÿ where the ÿ corresponds to the turn direction and t2 2 is computed directly from (11) and ÿ ÿ (13) 2 sin t2
2 qÿ 19 ÿ ÿ 2 ÿ ÿ ÿ 2 : 2 cos t2 ÿ 1 2 sin t2 Hence, the minimum-time solution (x t, y t, t), 04t4tf is de®ned in terms of the optimal switching times ti , i 1; 2; 3 and will normally consist of: i. a maximum rate turn (u t M , 1, 04t4t1 ) for timeÿ t1 to that radial (u t 0, ÿ t1 4t4t2 ) which intercepts the advancing safety ÿ zone at (x t2 , y t2 ), followed by either _ ii. a boundary arc trajectory (provided j t2 j< uM ) de®ned by the dierential±algebraic system (12), (14) and (18) for t2 4t4t3 where t3 is the time of exit from the boundary arc onto the ®nal radial (uÿ t 0; t3 4t4tf ) to the ¯y-to-point (xf , yf ) [cf. Fig. 3(a)], or iii. (if j t2 j 5uM ) a maximum rate turn (u t ÿuM , t2 4t4t3 to the ®nal radial (u t 0, t3 4t4tf ) to the ¯y-to-point (xf , yf ) [Fig. 4(a)]. It should be emphasized here that (i), (ii) and (iii) above represent the optimal control regime for most but not all con¯ict avoidance con®gurations. Considerably more complex control histories are possible. For example, since _ ÿ tis either strictly increasing or strictly decreasing on a boundary arc, it may happen that j _ tB j uM for some tB > t2 on an arc where t is increasing, in which case the boundary arc turn would switch at tB to a maximum rate turn [i.e. add the additional switching time tB . An even more (iii)] to the ®nal radial to (xf , yf ). This ÿ would involved scenario would result if j _ t2 j 5uM [i.e. (iii)] and the radial from the second turn to (xf , yf ) intersects the boundary arc (which is displayed in Figs. 3 and 4 as a hatched line). In this case, the second turn would be to the radial to the boundary arc followed by a boundary arc turn [i.e. (ii)] to the ®nal radial to (xf , yf ) and would introduce two additional switching times. However, since all of these possibilities just result in rearrangements of (ii) and (iii) above, the con®gurations analysed in Section 4 will be restricted to the normal control history given by (i)±(iii). 4. Calculation of the extremal trajectories In what follows, we describe the calculations required to compute the feasible extremal solutions (x t, y t, t), 04t4tf for the normal control histories given in Section 3. Integration
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
273
of (1) gives x t x0 2 = uM sin
0
y t y0 2 = uM ÿ cos
t1 uM ÿ sin 0
0
t1 uM cos
2 cos 0
0
2 sin
t1 uM t ÿ t1 0
t1 uM t ÿ t1
20 21
for (x t, y t) on t1 4t4t2 where 1 de®nes the direction of the initial maximum rate ®rst turn and t1 is the time in the ®rst turn to the exit point (x t1 , y t1 ). t2 is the total time to the point (x t2 , y t2 ) of interception with the advancing protective zone (x t2 ÿ 1 t2 2 y2 t2 r2p [cf. Fig. 3(a)]. The calculation of t2 , (x t2 , y t2 ) and the corresponding angle of interception with the zone is also not straightforward. This is because the interception point must ®rst be calculated for the problem transformed to a coordinate system centered on aircraft #1 (i.e. using x ÿ 1 t and y). For the transformed problem, the optimal trajectory will consist of a turn to (1 ,1)=(x(t1ÿu1t1),y(t1)) followed by a radial which intercepts the now stationary zone 2 2 r2p at a point of tangency (2 ; 2 ) and then turns along a boundary arc of that circle. For any ®xed t1 , the tangent point (2 ,2 ) in the ±plane can be computed as the solution of the system 2 2 r2p ;
ÿ 1 ÿ 1 0
22
and is given by 2
r2p 1 r2p 21 21 ÿ r2p 12 1=2 12 21
;
2 r2p ÿ 2 1 =1
23
where 1 ; 1 x t1 ÿ 1 t1 ; y t1 ) and the sign in (23) is determined by whether the trajectory is to be executed fore or aft of the protective zone. Thus, t1 and t2 for each feasible extremal can be computed uniquely as the solution of the nonlinear system of algebraic equations ab ÿ1 f1 t1 ; t2 t2 ÿ t1 ÿ h1 t1 ; t2 0; g1 t1 ; t2 cos 0 24 k a kk b k where a cos
0
t1 M ; sin
0
t1 M ;
b x t2 ÿ x t1 ; y t2 ÿ y t1 ; ÿ
x t2 ÿ x t1 2 y t2 ÿ y t1 2 h1 t1 ; t2 2
1=2
and x t2 2 1 t2 and y y t2 2 are given by (23) for each t t2 . h1 t1 ; t2 is the time on the radial from the (x t1 , y t1 ) to x t2 , y t2 ). g1 t1 ; t2 is the angle between the heading vector a at (x t1 , y t1 ) and the heading vector b along the radial from (x t1 , y t1 ) to (x t2 , y t2 ). Here denotes the usual vector dot product. For the t1 t1 and t2 t2 for which g1 0 and f1 0, the heading vector out of the ®rst turn a lines up with the radial to the interception point b
274
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
and the time in the ®rst turn plus the time on the radial equals the time the safety zone has been moving. The feasible extremal corresponding to 1 are now completely deterÿ trajectories ÿ mined up to the point (x t2 , y t2 ) of interception with the protective zone by (20) and (21) for t1 t1 and t1 4t4t2 . While there are many formulations which could be used to solve for t1 and t2 , the particular form of the functions f1 and g1 in (24) is important for computational considerations. f1 t1 ; t2 is monotone increasing for t2 > t1 , t1 ®xed, and g1 t1 ; t2 is strictly convex in t1 for ®xed t2 . This means that a unique solution (t1 , t2 ) always exists and that very fast iterative solvers can be used with completeÿreliability ÿ to compute these values. Once the interception point (x t2 , y t2 ) has been determined, ÿ there are normally two possibilities for the remaining portion of a feasible extremal. If j t2 j 4uM , then the remaining trajectory consists of a boundary arc turn to the ®nal radial to (xf ; yf ). This ®rst case requires the numerical integration of the initial value problem (12), (18) until the ®nal radial condition cd ÿ1 g2 t3 cos 0 25 k c kk d k is satis®ed at t3 t3 where ÿ ÿ c cos t3 ; sin t3 ; and
t3 ÿ x t3 x t2 2 cos
ÿ d xf ÿ x t3 ; yf ÿ y t3
tdt;
t3 ÿ y t3 y t2 2 sin
t2
26
tdt
27
t2
are computed numerically using t from (12) and (18). g2 t3 is also strictly convex here and is just the condition that the heading vector c out of the second turn lines up with the radial vector d to (xf , yf ). Note again that for examples 1 and 2 considered in Section 5, the boundary arc is plotted as a hatched line for t > t3 . That is, the numerical solution of (12) and (18) is extended beyond t3 to indicate the envelope generated by constraint (4). The total time tf including the time on the ®nal radial is then given by tf t3
ÿ 2 ÿ 2 1=2 xf ÿ x t3 yf ÿ y t3 2
:
28
Thus, the optimal control law for 04t4tf in the case involving an interval t2 4t4t3 of singular boundary arc control is u t uM u t 0 u t 12; 14 and 18 u t 0
04t4t1
t1 4t4t2
t2 4t4t3
t3 4t4tf :
29
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
275
Theÿ second possibility for the remaining portion t2 4t4tf of a feasible extremal occurs when j _ t2 j> uM , that is, when the rate of turn required to follow the boundary arc of the advancing zone exceeds uM . Clearly, this will always hold when 1 > 2 . In this second case, the ®rst radial is followed by a maximum rate turn to the ®nal radial to (xf , yf ) and t3 is computed as the solution of the ®xed point Eq. (25) where now c cos
1
ÿ t3 uM ; sin
1
ÿ t3 uM
ÿ d xf ÿ x t3 ; yf ÿ y t3 ÿ x t3 x t2 2 = ÿuM sin
1
ÿ y t3 y t2 2 = ÿuM ÿ cos 1
0
ÿ t3 uM ÿ sin 1
1
ÿ t3 uM cos
1
t1 uM
and again g2 is a strictly convex function of t3 . tf is computed as in (28) and the optimal control law for this second case has the simple form u t M u t 0 u t ÿ M u t 0
04t4t1
t1 4t4t2 t2 4t4t3
30
t3 4t4tf :
5. Example applications In order to illustrate the solution procedures developed here, three example con¯ict avoidance con®gurations are considered. A Fortran 77 program was generated to carry out the calculations de®ned by Eqs. (8)±(30). Table 1 contains the data for each example con®guration ÿ x0 ; y0 ; xf ; yf ;
0 ; 1 ; 2 ; uM ; rp ; rM ; tm ; dm
which includes the turn radius rM u M2 of aircraft #2 and the time tm and distance dm of minimum separation in the event an avoidance maneuver is not executed. Table 2 contains the output information for each of the three examples. Distances and positions are given in nautical miles (nm) from a ®xed origin in the xy-plane where true North is de®ned by the direction of the positive y-axis. Time is measured in hours and speed is in knots. Headings are measured in radians. The output required for the optimal control and display (Figs. 3±5) of the con¯ict resolution maneuvers are the switching times ti , i 1; 2; 3 and tf which are computed to within an
276
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
Table 1 Aircraft con®guration data for examples 1±3 Data (x0 , y0 ) (xf , yf ) 0
v1 v2 uM rp rM v2 =uM tm dm
Example 1
Example 2
Example 3
Units
(15,ÿ15) (15,15) /2 250 350 80 9 4.37 175.13 3.48
(15,ÿ15) (10,15) /2 300 450 70 9 6.42 131.09 1.81
(15,ÿ15) (15,15) /2 300 200 80 9 2.50 207.70 4.16
nm nm radians knots knots radians/h nm nm s nm
Table 2 Switching times t (in s) for the feasible extremals in examples 1±3 Example 1
Example 2
Example 3
t
Aft
Front*
Aft
Front*
Aft*
Front
t1 t2 t3 tf
55.10 117.81 208.81 400.44
18.45 199.97 267.17 355.29
60.41 91.34 185.54 318.21
28.43 151.07 208.47 296.48
34.20 151.74 197.96 592.87
± ± ± ±
Fig. 5. Feasible extremal trajectories for example 3.
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
277
error tolerance of 3.010ÿ6 h for each feasible extremal and are given in Table 2 in seconds. In order to be consistent with the examples considered in (Friedman, 1991), the radius of the safety zone about aircraft #1 was take to be 9 nautical miles. The CPU times on a Macintosh IIci microcomputer required to compute the solutions for the examples considered here ranged between 450 and 1350 ms which is well within the time scale required for the ``real-time'' automatic calculation and control of the optimal con¯ict avoidance maneuvers. The results for each of the examples 1±3 are plotted in Figs. 3±5, respectively with the extremal trajectory passing aft of the zone labeled (a) and that passing in front of the zone labelled (b) in Figs. 3 and 4. Since 2 < 1 , in example 3, there is no feasible extremal which can pass on front of the moving zone. In each ®gure, the solid dot indicates the initial position (0, 0) of aircraft #1 which is moving to the right along the x-axis with airspeed 1 , the lower aircraft symbol indicates the initial position (x0, y0)=(15, ÿ15) of aircraft #2 and the x-box gives the position of the ¯y-topoint (xf , yf ). The safety zone at that instant of time (t2 ) when aircraft #2 intercepts the protective zone about aircraft #1 is plotted as a dashed-line circle and the position ( 1 t2 , 0) of aircraft #1 at symbol. The hatched-line trajectory which extends from time t2 is denoted by the ÿupper aircraft ÿ the interception point (x t2 ; y t2 ) is the locus of the boundary arc trajectory [i.e. the numerical solution of (12), and (18)] which in each case has been extended beyond the exit point from the boundary arc turn to indicate the full extent of the envelope determined by the state variable constraint. The solid line denotes the extremal trajectory C : f x t, y t j 04t4tf g and the arrow the direction of ¯ight. Tick marks associated with the switching times ÿ ti ?, ÿ i 1; 2; 3 for each of the feasible extremals are placed at the corresponding points (x ti , y ti ) along the trajectory C . In examples 1 and 3, uM is taken to be uM 80 radians/h which corresponds to a gentle (comfortable) turn radius for aircraft #2. In example 2, uM is taken to be uM 70 radians/ h which results in a slightly extended turn radius rM for aircraft #2. For each example, the optimal (minimum-time) solution is indicated with an asterisk in Table 2 and the minimum ®nal times are indicated in bold print. Example 1. Example 1 illustrates the control sequence (i)±(ii) in Section 3 where both feasible extremals [Fig. 3(a) and (b)] exist since 2 350 > 1 250 and each has an interval (t2 , t3 ) of boundary arc control. As can be seen from Table 2, the optimal trajectory is that extremal which passes in front of the zone [Fig. 3(b)] and requires 355.29 s to execute with 85 s of turning time. The second feasible extremal passes aft of the zone [Fig. 3(a)] and requires 400.44 s with 146 s of turning time. It should be noted again that in all of the ®gures the ground track of aircraft #1 is plotted for time [0, t2 ] while the protective zone is only plotted at the time t2 of interception with the ®rst radial. Consequently, while it may at ®rst appear thatÿthe zoneÿ about aircraft #1 has been violated, this is not the case and the ®rst point of contact (x t2 ; y t2 ) occurs at time t2 . Example 2. Example 2 illustrates the control sequence (i)±(iii) and (i)±(ii) in Section 3ÿ where the _ t2 j> uM ) ®rst feasible extremal [Fig. 4(a)] consists of two maximum rate turns (because j while the second extremal [Fig. 4(b)] involves an interval (t2 , t3 ) of singular boundary arc control. Here the optimal trajectory is again the extremal which passes in front of the zone and requires 296.48 s to execute with 86 s of turning time. The second feasible extremal passes aft of the zone [Fig. 2(a)] and requires 318.21 s with 154 s of turning time. Example 3. Since 2 < , here, example 3 illustrates the control sequence (i)±(iii) where the only feasible extremal is the trajectory which passes aft of the advancing protective zone (Fig. 3) and consists of two maximum rate turns. This extremal is the optimal trajectory and requires 592.87 s
278
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
to execute with only 46 s of turning time. Note that the hatched line corresponding to the boundary arc trajectory is not plotted in Fig. 3 because such a trajectory is not possible when 2 < 1 . 6. Generalizations and extensions The ÿassumptions that jyfj>rp and that the initial separation is at least d0 rp 2 =uM rp = 2 1 in the formulation of the simplest model in Section 2 are not restrictive in any real sense but are included only to exclude the case where aircraft #2 is initially so close (and directly in front of) the advancing protective zone about aircraft #1 that no avoidance maneuver would be possible within the control constraint set (3). Here 2 =uM is the turn radius of aircraft #2 and rp = 2 is the time required for aircraft #2 to ¯y a distance rp . This does not exclude the case, for example, where aircraft #2 is approaching aircraft #1 on a direct collision course along the x-axis, it merely guarantees that if this is the case, then aircraft #2 will have sucient maneuverability to avoid the safety zone about aircraft #1. Con®gurations which violate these assumptions would be unrealistic in an ATC setting and could be easily handled as special cases using the theory developed here. For any given pairwise con¯ict, the minimum-time safe avoidance maneuver would normally be computed for each aircraft under the assumption that the true heading of the other aircraft remains constant for the time required for the maneuver. In the literature this is referred to as ``cooperative'' con¯ict resolution. The obvious choice would then be to select that aircraft and maneuver corresponding to shortest of the ®nal times tf . Given that the precise ``optimal'' safe control of any potential pairwise con¯ict can be computed in milliseconds using only Eqs. (8)±(30), algorithms similar to those found in Bowers (1995) and Niedringhaus, (1991) can be generated to analyse and resolve a very large number of interrelated potential pairwise separation problems in a controlled airspace. There is the further possibility here of exploiting the inherent parallel structure of the problem. Optimal criteria other than minimum time and more general geometrical con®gurations for the protected zones are certainly possible but there are computational advantages if the generalized zones are de®ned in terms of smooth manifolds. The dynamics of the con¯ict model can readily be generalized to include both heading and airspeed changes as well as local wind speed Wi and direction i as in, for example, x_ i t i t cos
i t
ÿ Wi cos i
y_i t i t sin
i t
ÿ Wi sin i ; 04t4tf
_ i t ui t;
j ui t j 4ui;M ;
i t ai t;
i;m 4 i t4 i;M
i 1; :::; n ai;m 4ai t4ai;M ;
where the bounds i;M , i;m , i;M , ai;m and ai;m , i 1; :::; n de®ne the control and state variable constraint sets. While it is straightforward to also incorporate roll rates as was done in Clements, (1993) as well as optional constrained altitude changes, these generalizations add a signi®cant
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
279
level of complexity to the calculations required to solve the resulting dynamical system and their consideration has been left for a later study. 7. Discussion This work addresses the complex problem of resolving potential pairwise con¯icts in a multiple aircraft environment. The main objective of this study was to explore the application of some standard techniques from singular optimal control theory to the development of fast, accurate and robust calculation procedures for determining optimal safe collision avoidance/con¯ict resolution maneuvers for aircraft traversing planar intersecting trajectories. Here the avoidance maneuver was termed ``optimal'' if it represented the minimum-time safe deviation from a preassigned ¯ight program (see e.g. Friedman, 1991). Eqs. (6) and (7) yield a necessary and sucient condition for the existence of a pairwise con¯ict between two aircraft. If such a potential con¯ict does exist, Eqs. (8)±(30) can be used to determine the minimum-time safe steering program t, 04t4tf for the avoidance maneuver. In particular, the precise form of the boundary arc control u t _ t, t2 4t4t3 was found to be given by the solution of the dierential±algebraic system (12), (14) and (18). The corresponding boundary arc turns represent the minimum envelope that a feasible trajectory can execute and will be an essential calculation in any con¯ict avoidance program, particularly regarding trajectories which may pass in front of the moving protective zone. An interesting and somewhat counterintuitive result regarding the optimal radial which intercepts the moving zone is that the extremal trajectory is not tangent to the zone at the point of interception. The computed optimal steering program u t _ t; 04t4tf incorporates the speci®c ¯ight capabilities and constraints of each maneuvering aircraft and, as is shown for the three example con¯ict con®gurations considered in Section 5, requires only limited computing time and resources. Consequently, these solutions could be continuously recomputed on aircraft navigation computers or at ATC centres to adjust for any ¯ight program errors encountered during the avoidance maneuver. This precise control of the minimum-time avoidance trajectory ensures increased safety and could result in reduced operating expenses. While it has been shown here that optimal control theory and methods can be applied eectively to optimally resolve planar con¯icts, much more work is required before an informed assessment can be made regarding the feasibility of incorporating such an approach into a comprehensive air trac management program. Acknowledgements This research was partially supported by the Natural Sciences and Engineering Research Council of Canada.
References Bowers, K.L., 1995. Mathematical collision avoidance models for air trac in a free ¯ight environment. M.Sc. thesis, Dalhousie University, Halifax, NS, Canada.
280
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
Clements, J.C., 1993. Real-time optimal control of aircraft turn trajectories. J. Guidance 16(2), 385±387. Cotton, C.W.B., 1995. Free ¯ight in domestic ATM. The Journal of Air Trac Control 10±13. Friedman, M.F., 1988. Decision analysis and optimality in air trac control con¯ict resolution: I. Optimal timing of speed control in a linear planar con®guration. Transportation Research 22B(3), 207±216. Friedman, M.F., 1989. On the exact distribution of the number of perceived con¯icts with prescribed duration in air trac control. Transportation Research 23B(4), 245±256. Friedman, M.F., 1991. Decision analysis and optimality in air trac control con¯ict resolution: II. Optimal heading (vectoring) control in a linear planar con®guration. Transportation Research 25B(1), 39±53. Lanzer, N., Jenny, M.T., 1995. Managing the evolution to free ¯ight. The Journal of Air Trac Control, 18±20. Massoglia P.L., Pozesky M.T., Germana, G.T., 1989. The use of satellite technology for oceanic air trac control. Proceedings of the IEEE 77(11), 1695±1708. Maurer, H., 1977. On optimal control problems with bounded state variables and control appearing linearly. SIAM J. Control Optim. 5(3), 345±362. NASPG (North Atlantic Systems Planning Group), 1995. Progress in the Development of a Collision Risk Mode for the Future NAT CNS/ATM Environment. Technical report NAT SPG/31 WP/4. International Civil Aviation Organization, Paris. Niedringhaus, W.P., 1991. Maneuver option manager: automated simpli®cation of complex air trac control problems. IEEE Transactions on Systems, Man, and Cybernetics 22(5), 1047±1057. Rigotti K.D., Dean G.C., 1992. Algorithms for Con¯ict Prediction. Technical report RAD/902/15/4.1, AD4 Division, Dra Malvern. Computer Assistance for Enroute, London. Tomlin, C., Pappas, G., Sastry, S., 1997. Noncooperative con¯ict resolution. Proceedings of the 36th WEE Conference on Decision and Control, San Diego, CA. Tomlin, C., Pappas, G., Sastry, S., 1998. Con¯ict resolution in air trac management: a study in multi-agent hybrid systems. IEEE Transactions on Automatic Control 43(4), 509±521. Williamson, T., Spencer, N., 1989. Development and operation of the trac alert and collision avoidance system (TCAS). Proceedings of the IEEE 77(11), 1735±1744.
Transportation Research Part B 33 (1999) 265±280
The optimal control of collision avoidance trajectories in air trac management John C. Clements * Department of Mathematics, Statistics, and Computing Science, Dalhousie University, Halifax, Nova Scotia, B3H 3J5, Canada Received 8 August 1996; received in revised form 21 June 1998
Abstract A collision avoidance/con¯ict resolution model for aircraft traversing planar intersecting trajectories is developed based on the actual ¯ight dynamics of each aircraft and the criteria that a ¯ight regime correction is optimal if it represents a minimum-time safe deviation from a preassigned ¯ight program. The formulation of the model leads to a singular optimal control problem where the control variable constraint is the maximum acceptable aircraft turn rate for passenger comfort and safety and the state variable constraint is de®ned in terms of the radius of the protected zone about a potentially con¯icting aircraft. A robust, accurate real-time solution procedure is derived for computing the steering program for the execution of the optimal safe avoidance maneuver. Several representative examples of pairwise trajectory con¯icts are analysed in detail. A brief discussion of how the model might be generalized to incorporate changes in airspeed and altitude as well as more comprehensive dynamics is included. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Air trac management; Optimal control
1. Introduction Although current systems are already highly sophisticated, air trac control in the future will become even more complex and will he based on ``user-preferred trajectories'' and the continuous automatic monitoring and control of ¯ight trajectories both from within the aircraft and via the relevant air trac control centers (Bowers, 1995; Cotton, 1995; Friedman, 1989; Massoglia et al., 1989; Niedringhaus, 1991; NASPG, 1995; Lanzer and Jenny, 1995; Rigotti and Dean, 1992; Tomlin *
Tel.: +1-902-494-3381; fax: +1-902-494-5130; e-mail: [email protected].
0191-2615/98/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0191-2615(9 8 ) 0 0 0 3 1 - 9
266
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
et al., 1997, 1998; Williamson and Spencer, 1989). Identifying potential con¯icts and determining adjustments to the aected ¯ight trajectories which are ``optimal'' in some sense will be an important component of such systems. A number of collision avoidance models have been examined, notably the linear heading (Friedman, 1991) and speed control (Friedman, (1988) avoidance strategies developed by Friedman, the linear programing approach taken by Niedringhaus (1991) for the resolution of multiple interrelated separation problems and, more recently, the presentation of a con¯ict resolution architecture for multi-agent hybrid systems by Tomlin et al. (1998). In this work, we present an optimal control con¯ict resolution model for aircraft traversing planar intersecting trajectories which diers from the studies listed above in that it incorporates the ¯ight dynamics of each aircraft and the resulting avoidance maneuvers are smooth trajectories de®ned by the turn rate constraints of the avoiding aircraft. The model is based on the criteria that a ¯ight regime correction is optimal if it represents a minimum~time safe deviation from a preassigned ¯ight program. The formulation of the model leads to a singular optimal control problem where the control variable constraint (3) is the maximum acceptable aircraft turn rate for passenger comfort and safety under normal operating conditions and the non-autonomous state variable constraint (4) is de®ned in terms of the radius rp of the protective zone about a potentially con¯icting aircraft. In Section 2, we con®ne our attention to the simplest model which involves two aircraft operating at constant altitude and constant speci®ed airspeeds and for which only heading changes of one aircraft are permitted. Without loss of generality, we can assume that aircraft #1 is travelling at airspeed 1, in the positive x-direction while aircraft #2 must execute an avoidance maneuver to a speci®ed designated objective (xf , yf ) at airspeed 2 (cf. Fig. 1). Here (xf , yf ) could be the original ¯ight plan destination of aircraft #2 or merely a point chosen on the intended ¯ight path of that aircraft which is beyond the region of potential con¯ict. The positions of aircraft #1 and #2 in the xy coordinate plane at time t 0 are taken to be (0, 0) and (x0 ; y0), respectively, and 0 is the initial heading (and ground track) angle of aircraft #2. The hatched circle of radius rp in Fig. 1 is the protective zone about aircraft #1. The optimal con¯ict avoidance trajectory is the minimum-time path C : f x t; y t j 04t4tf ; x ; y 2 C1 0; tf g from (x0 ; y0 ) to (xf ; yf ) which does not violate the (moving) protective zone x ÿ 1 t2 y2 r2p about aircraft #1. A robust, accurate real-time solution procedure based on an analysis of the related Hamiltonian and costate equations for the problem is developed in Sections 3 and 4. The calculation of the interval of singular control and the determination of the precise boundary arc control law (i.e. when the state variable constraint is strictly satis®ed) are not at all straightforward and are the main focus of this work. In particular, it is shown that determination of the boundary arc control requires the solution of the coupled dierential±algebraic system of equations given by (12), (14) andÿ (18). of the optimal trajectories concerns the point An ÿ interesting and important ÿ property t2 at which they intercept the boundary of the moving (x t2 , y t2 ) and heading angle protective zone. This is best understood by considering the problem transformed to a coordinate system centred on aircraft #1 using x ÿ 1 t and y and is described in detail in Section 4. The numerical integration of system (12) and (18) is straightforward and is not a serious computational consideration. The model studied here has the advantage that solutions can he continuously recomputed to allow for corrections in cross-track error and positioning accuracy throughout the con¯ict resolution maneuver. Several representative examples are considered in Section 5. Examples 1 and 2 are concerned with the determination of the optimal avoidance
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
267
trajectory when 2 > 1 . Example 3 examines a con¯ict resolution con®guration where 2 < 1 . In the ®rst two examples, the boundary arc represents the envelope of feasible minimum-time trajectories and is plotted beyond the boundary arc turn to clearly indicate the general trajectory de®ned by the state variable constraint (cf. the hatched line in Figs. 3 and 4). The generalization of the model to incorporate possible changes in airspeed and altitude as well as more comprehensive dynamics is discussed brie¯y in Section 6. 2. Formulation of the simplest model Let Ck 0; tf and Cp 0; tf denote the linear spaces of k times continuously dierential functions and piecewise continuous functions, on [0; tf ], respectively. The ¯ight trajectory of aircraft #2
Fig. 1. Potential con¯ict con®guration in the xy coordinate plane.
Fig. 2. Relationship between t, t and t on a boundary arc.
268
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
Fig. 3. Feasible extremal trajectories for example 1.
Fig. 4. Feasible extremal trajectories for example 2.
from (x0 ; y0 ) to (xf ; yf ) is given by C : x t; y t j 04t4tf where the aircraft position coordinates x t and y t are elements of C1 0; tf and tf is the time required to reach (xf , yf ). Neglecting the inertial acceleration terms associated with changes in aircraft heading, the dynamics of the problem can be formulated in terms of the system of equations x_ t 2 cos t; with end conditions
y_ t 2 sin t;
_ t u t;
04t4tf
1
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
x 0 x0 ; y 0 y0 ;
0
0;
ÿ ÿ x tf xf ; y tf yf
269
2
where the heading or steering angle t 2 C 0; tf is given in radians measured positive counterclockwise from the x-axis. u t 2 Cp 0; tf , 04t4tf , is identi®ed as the control variable and the control constraint j u t j 4 M ;
04t4tf
3
is preset as the maximum acceptable turn rate under normal operating conditions. u t M corresponds to a maximum allowable rate of turn to the left and t ÿ M to a maximum turn rate to the right. The non-autonomous state variable constraint is given by S t; x t; y t r2p ÿ x t ÿ 1 t2 y2 t40;
04t4tf
4
and corresponds to the protective zone of radius rp about aircraft #1. Here rp might be assigned a dierent value for aircraft #1 in each pairwise con¯ict depending on the speed and maneuverability of that aircraft or it may be set to a ®xed value for all aircraft. Thus, the simplest optimal control problem is to determine the steering program t, 04t4tf , which guides aircraft #2 along that trajectory C which reaches the (xf , yf ) in minimum time tf without violating the protected zone of aircraft #1. That is, the problem is to minimize the ®nal time functional tf min J u dt
u2Cp 0;tf
5
0
subject to the state Eq. (1), the end conditions (2) and the control and state variable constraints (3) and (4). For computational simplicity, it is assumed here that the initial heading 0 of aircraft #2 represents the shortest ¯ight path to (xf , yf ). Thus, if the state variable constraint (4) is not violated, the optimal control is given by the default value u t 0; 04t4tf . If the initial heading 0 of aircraft #2 is maintained for a given con®guration, then the minimum separation distance dm between the two aircraft occurs at time tm x0 1 ÿ 2 cos 0 ÿ y0 2 sin 0 = 21 ÿ 2 2 1 cos 0 22 6 and is given by q dm x0 tm 2 cos 0 ÿ 1 tm 2 y0 tm u2 sin 0 2 :
7
Hence, the defaultÿ control u1 t 0; 04t4tf is optimal if and only if dm 5rp , in which case 2 tf xf ÿ x0 2 yf ÿ y0 2 =u2 . If dm < rp , then a feasible extremal is a minimum-time trajectory which satis®es the constraints of the problem. Regarding the initial aircraft con®guration, it is further assumed that j yf j> rp and that the initial aircraft separation is at least d0 rp 1 uM . As
270
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
shown in Section 6, this ensures (except in the case u t 0; 04t4tf ) there always exists at least one and possibly two feasible extremal trajectories for the optimal control problem. The ®rst of these extremal trajectories passes aft of the advancing protective zone x ÿ 1 t2 y2 4r2p about aircraft #1 while the second, if it exists, will pass in front of this advancing zone. Both extremal trajectories must be geodesics. The optimal trajectory is that feasible extremal requiring the least time. If 2 4 1 , then only the trajectory passing aft of the protective zone about aircraft #1 is feasible and hence optimal. If 2 > 1 , then both feasible extremals exist. 3. The optimal control u t; 04t4tf The formulation and analysis of the Hamiltonian for system (1)-(5) is similar to that found in except for the derivation of the boundary control law required to accommodate the time-dependent state variable constraint (4). We show here that this law is given by the coupled dierentialalgebraic system (12), (14) and (18). The Hamiltonian for (1)-(5) is given by H t 1 l1 tv2 cos t l2 tv2 sin t l3 tu t tS t; x t; y t
8
for Lagrange multipliers li t 2 C1 0; tf ; i 1; :::; 3 satisfying the adjoint equations l_ 1 t ÿ@H=@x 0 l_ 2 t ÿ@H=@y 0
9
l_ 3 t ÿ@H=@ l1 t 2 sin t ÿ l2 t 2 cos t for all t in [0; tf ] where the multiplier t50 satis®es ÿ ÿ tS t; x t; y t 0 on [0; tf ]. From the 0 and l t 0 , l transversality conditions, l 3 f i i tf , i 1; 2 and l3 0 are free parameters. A arc subarc of : x t; y t; t j 04t4tf for which S t; x t; y t < 0 is called an interior and a boundary arc is a subarc of II for which S t; x t; y t 0. t in (8) is in C 0; tf and is continuously dierentiable on the interior of any boundary arc (Clements, 1993; Maurer, 1977). This is a singular control problem with switching function t @H=@u l3 t. Since l_ 1 t l_ 2 t 0 in (9), l1 and l2 are constant for all t in [0; tf ]. Consequently, on any interior arc of , the optimal control law is u t ÿ M u t 0 u t M
l3 t > 0 l3 t 0 l3 t < 0:
10
and a feasible extremal must begin with a maximum rate turn (u t uM ; 1; 04t4t1 ) for radial (u t 0; t1 ? 4t4t2 ) which intercepts the moving zone at a time t1 toÿ that ÿ ÿ protective ÿ point (x t2 ; y t2 ). The details regarding the calculation of t2 and (x t2 , y t2 ) are given in Section 4 and show in particular that the radial (u t 0; t1 ? 4t4t2 ) is not tangent to the advancing circular zone at the ÿ point ÿ of interception. Here (t1 ?; t2 ) is an interval of singular control. The maneuver to (x t2 ; y t2 ) is followed by an opposite turn along the boundary of
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
271
the moving circle x ÿ 1 2 y2 t r2p to the ®nal radial to the ¯y-to-point (xf , yf ) [cf. Fig. 3(a)] unless the maximum allowable turn rate is exceeded at the point of interception, in which case a maximum rate turn (u t ÿuM ; 1; t2 4t4t3 ) to the ®nal radial to (xf ; yf ) is executed [cf. Fig. 4(a)]. If an interval (t2 ; t3 ) of boundary arc control does exist, then > 0 on this interval and the optimal control law can no longer be determined from the form of the Hamiltonian (8) and the adjoint equations (9). Determination of the precise form of the control law on any boundary arc of is fairly involved but essential for safety considerations, particularly in the case of a trajectory passing in front of the protective zone. Let (x t; y t) be the trajectory components of aircraft #2 on a boundary arc and let t denote the heading and T t the component of the speed of aircraft #2 in the direction tangent to the moving zone at (x t; y t) as in Fig. 2. Then the geometry gives 2 sin t T t sin t 2 cos t T t cos t 1
11
2 sin t ÿ t 1 sin t and it follows that t t ÿ sinÿ1 1 sin t= 2
12
and 2T t 2 cos t ÿ 1 2 2 sin t2 :
13
Dierentiation of (12) gives the form of the required control on the boundary arc 2
3
1 cos t= 2 6 7 u t _ t 41 ÿ q5 _ t 1 ÿ 1 sin t= 2 2
14
where it remains to solve for t. Let t denote the angle which (x t; y t) makes with the origin in Fig. 2. Since (x t; y t) must remain on the boundary of the moving circle, ÿ x t rp cos t 1 t; rp cos t =2 1 t; y t rp sin t rp sin t =2: 15 where the corresponds to a right (ÿ) or left turn (+) on the boundary trajectory. It then follows from (13) and (15) that 2T t r2p _ t2 22 21 ÿ 2 2 1 cos t which yields after some simpli®cation
16
272
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
ÿ r2p _ t2 22 21 ÿ 2 2 1 cos t ÿ sinÿ1 1 sin t= 2 q2 1 cos t ÿ 22 ÿ 21 21 cos2 t :
17
Thus, the boundary arc control is given by (14) where t is computed as the solution of the initial value problem q
_ t ÿ 1 cos t ÿ 22 ÿ 21 21 cos2 t =rp ; t > t2 18 ÿ
t2 2 ÿ where the ÿ corresponds to the turn direction and t2 2 is computed directly from (11) and ÿ ÿ (13) 2 sin t2
2 qÿ 19 ÿ ÿ 2 ÿ ÿ ÿ 2 : 2 cos t2 ÿ 1 2 sin t2 Hence, the minimum-time solution (x t, y t, t), 04t4tf is de®ned in terms of the optimal switching times ti , i 1; 2; 3 and will normally consist of: i. a maximum rate turn (u t M , 1, 04t4t1 ) for timeÿ t1 to that radial (u t 0, ÿ t1 4t4t2 ) which intercepts the advancing safety ÿ zone at (x t2 , y t2 ), followed by either _ ii. a boundary arc trajectory (provided j t2 j< uM ) de®ned by the dierential±algebraic system (12), (14) and (18) for t2 4t4t3 where t3 is the time of exit from the boundary arc onto the ®nal radial (uÿ t 0; t3 4t4tf ) to the ¯y-to-point (xf , yf ) [cf. Fig. 3(a)], or iii. (if j t2 j 5uM ) a maximum rate turn (u t ÿuM , t2 4t4t3 to the ®nal radial (u t 0, t3 4t4tf ) to the ¯y-to-point (xf , yf ) [Fig. 4(a)]. It should be emphasized here that (i), (ii) and (iii) above represent the optimal control regime for most but not all con¯ict avoidance con®gurations. Considerably more complex control histories are possible. For example, since _ ÿ tis either strictly increasing or strictly decreasing on a boundary arc, it may happen that j _ tB j uM for some tB > t2 on an arc where t is increasing, in which case the boundary arc turn would switch at tB to a maximum rate turn [i.e. add the additional switching time tB . An even more (iii)] to the ®nal radial to (xf , yf ). This ÿ would involved scenario would result if j _ t2 j 5uM [i.e. (iii)] and the radial from the second turn to (xf , yf ) intersects the boundary arc (which is displayed in Figs. 3 and 4 as a hatched line). In this case, the second turn would be to the radial to the boundary arc followed by a boundary arc turn [i.e. (ii)] to the ®nal radial to (xf , yf ) and would introduce two additional switching times. However, since all of these possibilities just result in rearrangements of (ii) and (iii) above, the con®gurations analysed in Section 4 will be restricted to the normal control history given by (i)±(iii). 4. Calculation of the extremal trajectories In what follows, we describe the calculations required to compute the feasible extremal solutions (x t, y t, t), 04t4tf for the normal control histories given in Section 3. Integration
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
273
of (1) gives x t x0 2 = uM sin
0
y t y0 2 = uM ÿ cos
t1 uM ÿ sin 0
0
t1 uM cos
2 cos 0
0
2 sin
t1 uM t ÿ t1 0
t1 uM t ÿ t1
20 21
for (x t, y t) on t1 4t4t2 where 1 de®nes the direction of the initial maximum rate ®rst turn and t1 is the time in the ®rst turn to the exit point (x t1 , y t1 ). t2 is the total time to the point (x t2 , y t2 ) of interception with the advancing protective zone (x t2 ÿ 1 t2 2 y2 t2 r2p [cf. Fig. 3(a)]. The calculation of t2 , (x t2 , y t2 ) and the corresponding angle of interception with the zone is also not straightforward. This is because the interception point must ®rst be calculated for the problem transformed to a coordinate system centered on aircraft #1 (i.e. using x ÿ 1 t and y). For the transformed problem, the optimal trajectory will consist of a turn to (1 ,1)=(x(t1ÿu1t1),y(t1)) followed by a radial which intercepts the now stationary zone 2 2 r2p at a point of tangency (2 ; 2 ) and then turns along a boundary arc of that circle. For any ®xed t1 , the tangent point (2 ,2 ) in the ±plane can be computed as the solution of the system 2 2 r2p ;
ÿ 1 ÿ 1 0
22
and is given by 2
r2p 1 r2p 21 21 ÿ r2p 12 1=2 12 21
;
2 r2p ÿ 2 1 =1
23
where 1 ; 1 x t1 ÿ 1 t1 ; y t1 ) and the sign in (23) is determined by whether the trajectory is to be executed fore or aft of the protective zone. Thus, t1 and t2 for each feasible extremal can be computed uniquely as the solution of the nonlinear system of algebraic equations ab ÿ1 f1 t1 ; t2 t2 ÿ t1 ÿ h1 t1 ; t2 0; g1 t1 ; t2 cos 0 24 k a kk b k where a cos
0
t1 M ; sin
0
t1 M ;
b x t2 ÿ x t1 ; y t2 ÿ y t1 ; ÿ
x t2 ÿ x t1 2 y t2 ÿ y t1 2 h1 t1 ; t2 2
1=2
and x t2 2 1 t2 and y y t2 2 are given by (23) for each t t2 . h1 t1 ; t2 is the time on the radial from the (x t1 , y t1 ) to x t2 , y t2 ). g1 t1 ; t2 is the angle between the heading vector a at (x t1 , y t1 ) and the heading vector b along the radial from (x t1 , y t1 ) to (x t2 , y t2 ). Here denotes the usual vector dot product. For the t1 t1 and t2 t2 for which g1 0 and f1 0, the heading vector out of the ®rst turn a lines up with the radial to the interception point b
274
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
and the time in the ®rst turn plus the time on the radial equals the time the safety zone has been moving. The feasible extremal corresponding to 1 are now completely deterÿ trajectories ÿ mined up to the point (x t2 , y t2 ) of interception with the protective zone by (20) and (21) for t1 t1 and t1 4t4t2 . While there are many formulations which could be used to solve for t1 and t2 , the particular form of the functions f1 and g1 in (24) is important for computational considerations. f1 t1 ; t2 is monotone increasing for t2 > t1 , t1 ®xed, and g1 t1 ; t2 is strictly convex in t1 for ®xed t2 . This means that a unique solution (t1 , t2 ) always exists and that very fast iterative solvers can be used with completeÿreliability ÿ to compute these values. Once the interception point (x t2 , y t2 ) has been determined, ÿ there are normally two possibilities for the remaining portion of a feasible extremal. If j t2 j 4uM , then the remaining trajectory consists of a boundary arc turn to the ®nal radial to (xf ; yf ). This ®rst case requires the numerical integration of the initial value problem (12), (18) until the ®nal radial condition cd ÿ1 g2 t3 cos 0 25 k c kk d k is satis®ed at t3 t3 where ÿ ÿ c cos t3 ; sin t3 ; and
t3 ÿ x t3 x t2 2 cos
ÿ d xf ÿ x t3 ; yf ÿ y t3
tdt;
t3 ÿ y t3 y t2 2 sin
t2
26
tdt
27
t2
are computed numerically using t from (12) and (18). g2 t3 is also strictly convex here and is just the condition that the heading vector c out of the second turn lines up with the radial vector d to (xf , yf ). Note again that for examples 1 and 2 considered in Section 5, the boundary arc is plotted as a hatched line for t > t3 . That is, the numerical solution of (12) and (18) is extended beyond t3 to indicate the envelope generated by constraint (4). The total time tf including the time on the ®nal radial is then given by tf t3
ÿ 2 ÿ 2 1=2 xf ÿ x t3 yf ÿ y t3 2
:
28
Thus, the optimal control law for 04t4tf in the case involving an interval t2 4t4t3 of singular boundary arc control is u t uM u t 0 u t 12; 14 and 18 u t 0
04t4t1
t1 4t4t2
t2 4t4t3
t3 4t4tf :
29
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
275
Theÿ second possibility for the remaining portion t2 4t4tf of a feasible extremal occurs when j _ t2 j> uM , that is, when the rate of turn required to follow the boundary arc of the advancing zone exceeds uM . Clearly, this will always hold when 1 > 2 . In this second case, the ®rst radial is followed by a maximum rate turn to the ®nal radial to (xf , yf ) and t3 is computed as the solution of the ®xed point Eq. (25) where now c cos
1
ÿ t3 uM ; sin
1
ÿ t3 uM
ÿ d xf ÿ x t3 ; yf ÿ y t3 ÿ x t3 x t2 2 = ÿuM sin
1
ÿ y t3 y t2 2 = ÿuM ÿ cos 1
0
ÿ t3 uM ÿ sin 1
1
ÿ t3 uM cos
1
t1 uM
and again g2 is a strictly convex function of t3 . tf is computed as in (28) and the optimal control law for this second case has the simple form u t M u t 0 u t ÿ M u t 0
04t4t1
t1 4t4t2 t2 4t4t3
30
t3 4t4tf :
5. Example applications In order to illustrate the solution procedures developed here, three example con¯ict avoidance con®gurations are considered. A Fortran 77 program was generated to carry out the calculations de®ned by Eqs. (8)±(30). Table 1 contains the data for each example con®guration ÿ x0 ; y0 ; xf ; yf ;
0 ; 1 ; 2 ; uM ; rp ; rM ; tm ; dm
which includes the turn radius rM u M2 of aircraft #2 and the time tm and distance dm of minimum separation in the event an avoidance maneuver is not executed. Table 2 contains the output information for each of the three examples. Distances and positions are given in nautical miles (nm) from a ®xed origin in the xy-plane where true North is de®ned by the direction of the positive y-axis. Time is measured in hours and speed is in knots. Headings are measured in radians. The output required for the optimal control and display (Figs. 3±5) of the con¯ict resolution maneuvers are the switching times ti , i 1; 2; 3 and tf which are computed to within an
276
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
Table 1 Aircraft con®guration data for examples 1±3 Data (x0 , y0 ) (xf , yf ) 0
v1 v2 uM rp rM v2 =uM tm dm
Example 1
Example 2
Example 3
Units
(15,ÿ15) (15,15) /2 250 350 80 9 4.37 175.13 3.48
(15,ÿ15) (10,15) /2 300 450 70 9 6.42 131.09 1.81
(15,ÿ15) (15,15) /2 300 200 80 9 2.50 207.70 4.16
nm nm radians knots knots radians/h nm nm s nm
Table 2 Switching times t (in s) for the feasible extremals in examples 1±3 Example 1
Example 2
Example 3
t
Aft
Front*
Aft
Front*
Aft*
Front
t1 t2 t3 tf
55.10 117.81 208.81 400.44
18.45 199.97 267.17 355.29
60.41 91.34 185.54 318.21
28.43 151.07 208.47 296.48
34.20 151.74 197.96 592.87
± ± ± ±
Fig. 5. Feasible extremal trajectories for example 3.
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
277
error tolerance of 3.010ÿ6 h for each feasible extremal and are given in Table 2 in seconds. In order to be consistent with the examples considered in (Friedman, 1991), the radius of the safety zone about aircraft #1 was take to be 9 nautical miles. The CPU times on a Macintosh IIci microcomputer required to compute the solutions for the examples considered here ranged between 450 and 1350 ms which is well within the time scale required for the ``real-time'' automatic calculation and control of the optimal con¯ict avoidance maneuvers. The results for each of the examples 1±3 are plotted in Figs. 3±5, respectively with the extremal trajectory passing aft of the zone labeled (a) and that passing in front of the zone labelled (b) in Figs. 3 and 4. Since 2 < 1 , in example 3, there is no feasible extremal which can pass on front of the moving zone. In each ®gure, the solid dot indicates the initial position (0, 0) of aircraft #1 which is moving to the right along the x-axis with airspeed 1 , the lower aircraft symbol indicates the initial position (x0, y0)=(15, ÿ15) of aircraft #2 and the x-box gives the position of the ¯y-topoint (xf , yf ). The safety zone at that instant of time (t2 ) when aircraft #2 intercepts the protective zone about aircraft #1 is plotted as a dashed-line circle and the position ( 1 t2 , 0) of aircraft #1 at symbol. The hatched-line trajectory which extends from time t2 is denoted by the ÿupper aircraft ÿ the interception point (x t2 ; y t2 ) is the locus of the boundary arc trajectory [i.e. the numerical solution of (12), and (18)] which in each case has been extended beyond the exit point from the boundary arc turn to indicate the full extent of the envelope determined by the state variable constraint. The solid line denotes the extremal trajectory C : f x t, y t j 04t4tf g and the arrow the direction of ¯ight. Tick marks associated with the switching times ÿ ti ?, ÿ i 1; 2; 3 for each of the feasible extremals are placed at the corresponding points (x ti , y ti ) along the trajectory C . In examples 1 and 3, uM is taken to be uM 80 radians/h which corresponds to a gentle (comfortable) turn radius for aircraft #2. In example 2, uM is taken to be uM 70 radians/ h which results in a slightly extended turn radius rM for aircraft #2. For each example, the optimal (minimum-time) solution is indicated with an asterisk in Table 2 and the minimum ®nal times are indicated in bold print. Example 1. Example 1 illustrates the control sequence (i)±(ii) in Section 3 where both feasible extremals [Fig. 3(a) and (b)] exist since 2 350 > 1 250 and each has an interval (t2 , t3 ) of boundary arc control. As can be seen from Table 2, the optimal trajectory is that extremal which passes in front of the zone [Fig. 3(b)] and requires 355.29 s to execute with 85 s of turning time. The second feasible extremal passes aft of the zone [Fig. 3(a)] and requires 400.44 s with 146 s of turning time. It should be noted again that in all of the ®gures the ground track of aircraft #1 is plotted for time [0, t2 ] while the protective zone is only plotted at the time t2 of interception with the ®rst radial. Consequently, while it may at ®rst appear thatÿthe zoneÿ about aircraft #1 has been violated, this is not the case and the ®rst point of contact (x t2 ; y t2 ) occurs at time t2 . Example 2. Example 2 illustrates the control sequence (i)±(iii) and (i)±(ii) in Section 3ÿ where the _ t2 j> uM ) ®rst feasible extremal [Fig. 4(a)] consists of two maximum rate turns (because j while the second extremal [Fig. 4(b)] involves an interval (t2 , t3 ) of singular boundary arc control. Here the optimal trajectory is again the extremal which passes in front of the zone and requires 296.48 s to execute with 86 s of turning time. The second feasible extremal passes aft of the zone [Fig. 2(a)] and requires 318.21 s with 154 s of turning time. Example 3. Since 2 < , here, example 3 illustrates the control sequence (i)±(iii) where the only feasible extremal is the trajectory which passes aft of the advancing protective zone (Fig. 3) and consists of two maximum rate turns. This extremal is the optimal trajectory and requires 592.87 s
278
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
to execute with only 46 s of turning time. Note that the hatched line corresponding to the boundary arc trajectory is not plotted in Fig. 3 because such a trajectory is not possible when 2 < 1 . 6. Generalizations and extensions The ÿassumptions that jyfj>rp and that the initial separation is at least d0 rp 2 =uM rp = 2 1 in the formulation of the simplest model in Section 2 are not restrictive in any real sense but are included only to exclude the case where aircraft #2 is initially so close (and directly in front of) the advancing protective zone about aircraft #1 that no avoidance maneuver would be possible within the control constraint set (3). Here 2 =uM is the turn radius of aircraft #2 and rp = 2 is the time required for aircraft #2 to ¯y a distance rp . This does not exclude the case, for example, where aircraft #2 is approaching aircraft #1 on a direct collision course along the x-axis, it merely guarantees that if this is the case, then aircraft #2 will have sucient maneuverability to avoid the safety zone about aircraft #1. Con®gurations which violate these assumptions would be unrealistic in an ATC setting and could be easily handled as special cases using the theory developed here. For any given pairwise con¯ict, the minimum-time safe avoidance maneuver would normally be computed for each aircraft under the assumption that the true heading of the other aircraft remains constant for the time required for the maneuver. In the literature this is referred to as ``cooperative'' con¯ict resolution. The obvious choice would then be to select that aircraft and maneuver corresponding to shortest of the ®nal times tf . Given that the precise ``optimal'' safe control of any potential pairwise con¯ict can be computed in milliseconds using only Eqs. (8)±(30), algorithms similar to those found in Bowers (1995) and Niedringhaus, (1991) can be generated to analyse and resolve a very large number of interrelated potential pairwise separation problems in a controlled airspace. There is the further possibility here of exploiting the inherent parallel structure of the problem. Optimal criteria other than minimum time and more general geometrical con®gurations for the protected zones are certainly possible but there are computational advantages if the generalized zones are de®ned in terms of smooth manifolds. The dynamics of the con¯ict model can readily be generalized to include both heading and airspeed changes as well as local wind speed Wi and direction i as in, for example, x_ i t i t cos
i t
ÿ Wi cos i
y_i t i t sin
i t
ÿ Wi sin i ; 04t4tf
_ i t ui t;
j ui t j 4ui;M ;
i t ai t;
i;m 4 i t4 i;M
i 1; :::; n ai;m 4ai t4ai;M ;
where the bounds i;M , i;m , i;M , ai;m and ai;m , i 1; :::; n de®ne the control and state variable constraint sets. While it is straightforward to also incorporate roll rates as was done in Clements, (1993) as well as optional constrained altitude changes, these generalizations add a signi®cant
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
279
level of complexity to the calculations required to solve the resulting dynamical system and their consideration has been left for a later study. 7. Discussion This work addresses the complex problem of resolving potential pairwise con¯icts in a multiple aircraft environment. The main objective of this study was to explore the application of some standard techniques from singular optimal control theory to the development of fast, accurate and robust calculation procedures for determining optimal safe collision avoidance/con¯ict resolution maneuvers for aircraft traversing planar intersecting trajectories. Here the avoidance maneuver was termed ``optimal'' if it represented the minimum-time safe deviation from a preassigned ¯ight program (see e.g. Friedman, 1991). Eqs. (6) and (7) yield a necessary and sucient condition for the existence of a pairwise con¯ict between two aircraft. If such a potential con¯ict does exist, Eqs. (8)±(30) can be used to determine the minimum-time safe steering program t, 04t4tf for the avoidance maneuver. In particular, the precise form of the boundary arc control u t _ t, t2 4t4t3 was found to be given by the solution of the dierential±algebraic system (12), (14) and (18). The corresponding boundary arc turns represent the minimum envelope that a feasible trajectory can execute and will be an essential calculation in any con¯ict avoidance program, particularly regarding trajectories which may pass in front of the moving protective zone. An interesting and somewhat counterintuitive result regarding the optimal radial which intercepts the moving zone is that the extremal trajectory is not tangent to the zone at the point of interception. The computed optimal steering program u t _ t; 04t4tf incorporates the speci®c ¯ight capabilities and constraints of each maneuvering aircraft and, as is shown for the three example con¯ict con®gurations considered in Section 5, requires only limited computing time and resources. Consequently, these solutions could be continuously recomputed on aircraft navigation computers or at ATC centres to adjust for any ¯ight program errors encountered during the avoidance maneuver. This precise control of the minimum-time avoidance trajectory ensures increased safety and could result in reduced operating expenses. While it has been shown here that optimal control theory and methods can be applied eectively to optimally resolve planar con¯icts, much more work is required before an informed assessment can be made regarding the feasibility of incorporating such an approach into a comprehensive air trac management program. Acknowledgements This research was partially supported by the Natural Sciences and Engineering Research Council of Canada.
References Bowers, K.L., 1995. Mathematical collision avoidance models for air trac in a free ¯ight environment. M.Sc. thesis, Dalhousie University, Halifax, NS, Canada.
280
J.C. Clements / Transportation Research Part B 33 (1999) 265±280
Clements, J.C., 1993. Real-time optimal control of aircraft turn trajectories. J. Guidance 16(2), 385±387. Cotton, C.W.B., 1995. Free ¯ight in domestic ATM. The Journal of Air Trac Control 10±13. Friedman, M.F., 1988. Decision analysis and optimality in air trac control con¯ict resolution: I. Optimal timing of speed control in a linear planar con®guration. Transportation Research 22B(3), 207±216. Friedman, M.F., 1989. On the exact distribution of the number of perceived con¯icts with prescribed duration in air trac control. Transportation Research 23B(4), 245±256. Friedman, M.F., 1991. Decision analysis and optimality in air trac control con¯ict resolution: II. Optimal heading (vectoring) control in a linear planar con®guration. Transportation Research 25B(1), 39±53. Lanzer, N., Jenny, M.T., 1995. Managing the evolution to free ¯ight. The Journal of Air Trac Control, 18±20. Massoglia P.L., Pozesky M.T., Germana, G.T., 1989. The use of satellite technology for oceanic air trac control. Proceedings of the IEEE 77(11), 1695±1708. Maurer, H., 1977. On optimal control problems with bounded state variables and control appearing linearly. SIAM J. Control Optim. 5(3), 345±362. NASPG (North Atlantic Systems Planning Group), 1995. Progress in the Development of a Collision Risk Mode for the Future NAT CNS/ATM Environment. Technical report NAT SPG/31 WP/4. International Civil Aviation Organization, Paris. Niedringhaus, W.P., 1991. Maneuver option manager: automated simpli®cation of complex air trac control problems. IEEE Transactions on Systems, Man, and Cybernetics 22(5), 1047±1057. Rigotti K.D., Dean G.C., 1992. Algorithms for Con¯ict Prediction. Technical report RAD/902/15/4.1, AD4 Division, Dra Malvern. Computer Assistance for Enroute, London. Tomlin, C., Pappas, G., Sastry, S., 1997. Noncooperative con¯ict resolution. Proceedings of the 36th WEE Conference on Decision and Control, San Diego, CA. Tomlin, C., Pappas, G., Sastry, S., 1998. Con¯ict resolution in air trac management: a study in multi-agent hybrid systems. IEEE Transactions on Automatic Control 43(4), 509±521. Williamson, T., Spencer, N., 1989. Development and operation of the trac alert and collision avoidance system (TCAS). Proceedings of the IEEE 77(11), 1735±1744.