Trajectory tracking algorithms have many parameters that in uence their behavior. Adjustment of these parameters is complex because of their heterogeneous nature and the fact that the tracking quality depends on the parameters in a complicated and mutual manner. In the air tra c control area, there are many tracking quality criteria. An optimization program for this problem is elaborated on the basis of the genetic approach. With its use, the tracking quality of the real trajectory tracker has been signi cantly improved.
y(ti) = xG(ti) + (ti): The measurement y(ti) includes the random error (ti). A trajectory tracking algorithm maps the sequence fy(ti)gin=1 of the measurements up to the current instant tn to an estimate x^(tn) of the state vector x(tn).
The classical algorithm for the trajectory tracking is the Kalman lter [
The authors collaborate with the NITA company (Russia, Saint-Petersburg, http://nita.ru) which has the leader position in Russia in the ATM solutions. This company elaborated its own trajectory tracking program on the basis of the IMM algorithm. So, NITA is interested in optimizing the program work.
In the present paper, we consider the trajectory tracking program as a black
box and optimize its parameters without using any speci c features of the IMM
algorithm (in contrast to, for example, [
Unfortunately, we cannot fully describe the IMM implementation in the
NITA program and the IMM parameter list since we have the non-disclosure
agreement. But the structure of the IMM variables' dependence on parameters
is as follows. There are two main parts of the parameters: the rst part governs
the process noise covariance matrices Qk (for all k = 1; m, i.e., for all the
dynamics in the IMM method); the second one tunes the transition probabilities
matrix (TPM), whose elements k` are the probabilities of switches from kth
to the lth dynamics. Due to heterogeneity of the parameters and peculiarities of
the parametrization, the in uence of the parameters on the estimate sequence
fx^(ti)g is complex. For this reason, we decide to use the genetic approach [
In the optimization, simulated tracks are used. The model dynamics of an aircraft is as follows: 8x_ N = v cos '; > > ><x_ E = v sin '; >'_ = uv ; > >:v_ = w: (1) Here, xN , xE are the north and east coordinates of the aircraft in the plane, they make up together the \geometric" coordinate vector xG(ti) = [xN (ti) xE (ti)]T. The symbols v and ' are the velocity magnitude and the path angle, they make up together the vector xV (ti) = [v(ti) '(ti)]T of the velocity elements. The u and w are the lateral and longitudinal accelerations.
This dynamics is compliant with the common concepts of the aircraft
motion [
Performance Criteria in Target Tracking
In the ATM, there are conventional standards of the trajectory tracking quality.
All of them are based on the analysis of the di erences fx^(ti) x(ti)g between
the estimated x^(ti) and the true x(ti) states, which are collected along the
trajectory. Following the standard [
Unfortunately, there are some questions about the usage of the standards in real and model circumstances. 1. The accuracy of sensors is very di erent and di ers from the one described in the standards. 2. The maneuverability of aircraft is di erent too. This in uences the estimate error level and duration of the transitional processes. 3. Transitions between some motion segments are not described in the standards.
An attempt to answer these questions leads the authors to the use of the Cramer{ Rao lower bound (CRLB) instead of the standards. 2.2
Cramer{Rao Lower Bound as Universal Standard
The CRLB is the lower bound for the mean squared error (in a scalar case) and
the covariance matrix of the error (in a vector case) for the unbiased estimates.
The following matrix inequality holds [
E (x^(ti)
x(ti))(x^(ti)
x(ti))T
< J (ti):
The sign < denotes a partial order for symmetric matrices. Let A; B 2 Sm m,
x 2 Rm, then A < B , xTAx xTBx. The matrix J (ti) re ects the information
in the preceding measurements fyj gij=1 and can be calculated in a recurrent way
(see [
Though the estimates of the real algorithms including the IMM are biased, the CRLB shows the adequate behavior of the estimation error according to maneuvers and time. In the segments of a constant control, the CRLB speci es the error level that is hard to outperform. In Figure 1, the CRLB (red) and RMS deviation of the IMM algorithm (magenta) for the longitudinal channel are shown for the trajectory with one switch from the CV motion (u = w = 0) to the CT motion (u 6= 0, w = 0) at t = 50 s. As one can see in the gure, the decrease of the RMS error of the IMM lter after the maneuver is in accordance with the CRLB behavior.
A recurrent procedure is elaborated that calculates the matrices fJ (ti)g for
model tracks in the simulations. The analytical integrals from [
As the criteria for the optimization program, we decide to use the RMS deviation for values of this type: where
Xl(ti) =
xl(ti) qelT(ti)J (ti)el(ti)
; xl(ti) = elT(ti)(x^(ti) x(ti)): Here, Xl is the relative longitudinal deviation, that is, the ratio of the \absolute" (in meters) longitudinal deviation xl(ti) and the corresponding \projection" of the CRLB matrix onto the longitudinal channel de ned by the longitudinal direction vector el(ti) 2 R4. The vector el(ti) consists of the tangential vector to the curve xG( ) in the plane xN ; xE in the rst two coordinates and two zeros in the second two coordinates (the xV elements). The analogous formulas are introduced for the lateral (Xn), the velocity magnitude (Xv), and the path angle (X') channels using the corresponding directional vectors (en, ev, e').
Besides these \one-dimensional" criteria, we introduce additional two ones. The RMS deviation of the Mahalanobis distance in the \geometric" coordinates and the RMS deviation of the Mahalanobis distance in the total state vector space is
X4d(ti) = q (x^i x(ti))T(J (ti)) 1(x^i x(ti)): Here, JG(ti) is the upper-left block (two rst rows and columns) of the matrix J (ti), x^G(ti) is the \geometric" coordinates estimate (two rst elements of x^).
The total number of the criteria (cl, cn, cv, c', c2d, c4d) is nc = 6. In the case of iterating the criteria, they will be indexed as ci; i 2 1; nc. 2.3
Con dence Interval of Criterion In this section, we consider the longitudinal channel (deviation Xl, criterion cl) as an example. Calculations for other channels (criteria cn, cv, c', c2d, c4d) can be made in the same way.
The \true" criterion value cl is the expectation cl = q
E fXl2(t)g; with averaging over the distribution of random measurement errors, measurement instants, and probable trajectories.
For practical usage, we replace the expectation by its plug-in estimator ^cl
(see [
PN k=1 Nk k=1 j=1 1
N Nk X X Xl2(tj ; xk( )) =
n 1 X n i=1
Xl2;i ; c^l;n = qc^l2;n: (2) Here, N is the number of the trajectories xk( ) that have been processed up to the instant of c^l calculation, Nk is the number of the measurements of xk( ), n = PN
k=1 Nk is the total number of the measurements. The Xl values are indexed in
two ways: rstly, by the trajectory and the time instant as Xl(tj ; xk( )), secondly,
by the total index i over all the trajectories as Xl;i. The c^l;n satis es (see [
P c^l;n n!1! cl ;
E fc^l;ng = cl:
Since c^l;n 6= cl, we need a con dence interval for c^l;n for precise calculations.
An interval [c^ll;n; c^lu;n] 3 c^l;n is a 1 con dence interval if
lim inf P cl 2 [c^ll;n; c^lu;n]
n!1
1
:
In our work, we use the normal-based con dence interval [
s^en(c^l;n)z =2; c^l;n + s^en(c^l;n)z =2 ; where z =2 is the (1 (1 )=2)-quantile of the standard normal distribution N (0; 1), and s^en(c^l;n) is the estimated standard error of c^l;n, which can be calculated as follows:
n m^ l(;4n) = 1 X n i=1
Xl4;i ; s^en(c^l;n) =
1 2c^l;npn q m^ l(;4n) (c^l2;n)2 : 2.4
Optimization Problem Denote the parameter vector by a. In our case, there are 16 parameters of the IMM implementation by NITA; therefore, a 2 R16. All the parameters are subject of box constraints: amin a amax. The multicriteria optimization problem can be formulated as follows: ((cl(a); cn(a); cv(a); c'(a); c2d(a); c4d(a)) ! min; amin a amax: (3) 3
The genetic optimization algorithm is based on the evolution process simulation. In these terms, there is a population P = fag consisting of individuals. Every individual a is connected with some point in the space of parameters to be optimized (of the trajectory tracker in our case).
In Figure 2, the principal scheme of our genetic optimization program is
shown. The scheme is standard enough [
The \true" criterion values in (3) cannot be calculated; therefore, the estimates c^i;n(a)(a) are used instead: ( c^l;n(a)(a); c^n;n(a)(a); c^v;n(a)(a); c^';n(a)(a); c^2d;n(a)(a); c^4d;n(a)(a) amin a amax: ! min; Note that, here, the number of samples n for the estimates c^i;n(a)(a) depends on a: n = n(a). This re ects the fact that the estimates c^i of the criteria are calculated using di erent subsets of the whole training trajectory set X = fxk( )gkN=1.
In the program, we refuse the idea to calculate the estimates ^ci for a newly formed individual a using the whole set X . The set X is expected to be very big; therefore, if we did so, the computational load would be too large. Instead of that, at every generation b (this is the epoch of the evolution process), we calculate the \preliminary" estimates c^i;Nbatch (a) on a small batch Xb = fxk( )gkN=baNtcbhat(cbh+b1+)1 and then upgrade them at the succeeding generations using formulas (2).
The estimates c^i;n(a)(a) of the criteria are in uenced by random errors in
data; therefore, they have uncertainties, which are especially big if the
number n(a) is small. For this reason, the lower and upper con dence bounds ^cli;n(a)(a),
u
c^i;n(a)(a) are used in the selection procedure in order to prevent deletion of the
individuals that might have better values of c^i in the future. Other operations of
the genetic algorithm (the directed breeding, the crossover, and the mutation)
do not di er signi cantly from those in [
a2P a2P
Only non-immortal individuals of age greater than 3 participate in the second step of the selection. An individual as is removed from the population if its lower con dence bound of some criterion is greater than the population minimum of upper con dence bound of this criterion 9i 2 1; nc : c^i;n(as) > min c^iu;n(a): l a2P (4) In some cases (at the initial steps), condition (4) is not checked for all the criteria, but only for the criterion c4d. 4
There are two sets of simulated trajectories. The rst one is the training set for
optimization, the second one is the test set for validation of the optimization
results (see [
For every trajectory x( ) in the training set, a number of measurement realizations are made. The measurement instants are simulated as follows: initially, they are assigned in a grid with a certain step (can be from 3 to 11 s), then they vary randomly with the magnitude 0.45 of the grid step. To simulate misses of measurement, randomly chosen 5% of measurements are excluded from the collection. The measurement random errors i are drawn from the centered normal distribution with the variance matrix that is randomly chosen from 3 versions: 1) the major semi-axis is 200 m, the minor one is 100 m, the angle of the major semi-axis is 60 , 2) 300 m, 80 m, 30 , 3) 500 m, 150 m, 0 .
The genetic optimization program is adapted for parallel computation. The simulations were performed using the \Uran" supercomputer of N. N. Krasovskii Institute of Mathematics and Mechanics, Ural Branch of the Russian Academy of Sciences. Using 16 CPU cores, the total computation time were 58 hours. The computations was stopped at the 2320th generation of the evolution.
The test set consists of 400 trajectories which are created independently in the same way as the training ones (but without random variations of measurement instants). In Table 1, the criteria values on the test trajectories are shown.
An example of the trajectory tracker output in the plane xN ; xE for one test trajectory is shown in Fig. 4. The blue graph corresponds to the initial parameters, the green one corresponds to the optimized parameters at the generation 2320 (the best of the c4d criterion), the black dots are the measurements, and the true (ideal) trajectory is the magenta solid line. The estimated track for the optimized parameters is evidently closer to the true trajectory than that for the initial parameters. In Figure 5, the graphs of the RMS deviation of X4d is shown as a function of time along the trajectory shown in Fig. 4. The averaging of X4d is performed over 100 realizations of the random errors with the same grid of measurement instants. The blue solid line is for the initial parameters of the tracker and the red one is for the optimized parameters at the generation 2320. The dashed lines are the 0.95-con dence bands. The sub gure below the main gure shows the controls u(t) (green) and w(t) (blue). Using this sub gure, one can see that all the places where the behavior of X4d(t) changes are connected with the switches of the controls.
In Figure 6, the analogous RMS error graph for the absolute (in meters, without division by the value of the CRLB) lateral deviations xn(t) are shown (the trajectory is the same). As in the case of Fig. 5, the initial parameters correspond to the blue solid line and the optimized ones correspond to the red line. The dashed lines show the 0.95-con dence bands. The black line is for the RMS error of the measurements. 5
As a result of computer simulations, we can say that the trajectory tracking program with the optimized parameters shows a smaller RMS error (both absolute and relative) in each channel and has a faster transition process.
In further work on this project, we plane to use real ATC data for the training and test trajectory sets. Another point important for our study is the selection procedure, since its current version claims too high requirements on individuals.