EURASIP Journal on Applied Signal Processing 2003:8, 766–779 c ○ 2003 Hindawi Publishing Corporation Application of Evolution Strategies to the Design of Tracking Filters with a Large Number of Specifications

This


INTRODUCTION
A tracking filter has the double goal of reducing measurement noise and consistently predicting future values of signal. This kind of problems has efficient solutions in the case of stationary signals, but solutions for nonstationary problems are not so consolidated yet. This is the case in the field we are dealing with in this paper, tracking aircraft trajectories from radar measurements in air traffic control (ATC) applications.
The design of tracking filters for the ATC problem demands complex algorithms, like the modern interacting multiple model (IMM) filter [1]. These algorithms depend on a high number of parameters (seven in the IMM design presented here) which must be adjusted in order to achieve, as much as possible, the desired tracking filter performance. IMM has proven certainly satisfactory performance for tracking maneuvering targets, in relation to previous approaches. However, the relation between its input parameters and final performance is far from clear due to strongly nonlinear interactions among all parameters. Therefore, no direct design methodology has been proposed to generate the best solution for a specific application to date, apart from manual parameterization and evaluation with simulation.
Besides, real-world applications of tracking filters for ATC usually address performance specifications defined over an exhaustive set of realistic operational scenarios and covering a number of conflicting figures of merit. These two characteristics, large table of specifications and application of complex algorithms, make the design of modern tracking filter a very complex problem.
In this paper, the authors expose a new methodology to design and adjust tracking filters for ATC applications based on the use of evolution strategies (ES) as an optimization problem over a customized cost function (fitness function). The method has been demonstrated by the design of a realworld engineering application: a modern ATC system promoted by EUROCONTROL for Europe, the ARTAS system. Due to the high dimensionality of parameters' space and the large number of defined constrains (the operational scenarios and performance figures sum up to 264 specifications for ARTAS), an automatic procedure to search and tune the final solution is mandatory. Classical techniques, such as those based on gradient descent, were discarded due to the high number of local minima presented by the fitness function. ES have been selected for this problem due to their high robustness and immunity to local extremes/discontinuities in the fitness function.
However, the selection of a fitness function taking account of all specifications is not so direct since all of them should be simultaneously considered to guide the search. The performance of ES has been analyzed in previous works for sets of test functions, but its application to a real engineering problem with hundreds of specifications, where the fitness landscape's properties are not well known, is a harder task. A procedure has been proposed to build this function, exploiting specific knowledge about the domain. Objectives with similar behavior in the search are grouped first to select the worst cases for each group, and then combine all of them in the final cost function. Results show that this procedure is able to find acceptable solutions lowering the excess over some specifications as much as possible.
The paper starts by presenting the design performance constrains for ATC problems in Section 2 (particularized for an industrial application, the ARTAS system) and a description of the IMM algorithm in Section 3. In Section 4, we explain the proposed optimization method based on ES. Finally, Sections 5 and 6 are aimed at discussing optimization results and characteristics of solutions minimizing the fitness function, and summarizing the main conclusions.

SPECIFICATIONS FOR TRACKER MODULE OF ARTAS SYSTEM
ARTAS [2] is the concept of a Europe-wide distributed surveillance system developed by EUROCONTROL, relying on the implementation of interoperable units coordinated together. Each ARTAS unit will be in charge of processing all surveillance data reports (i.e., primary and secondary radar reports, ADS reports, etc.) to form a good estimate of the current air traffic situation in its responsibility volume.
Each of the ARTAS units should fulfill a set of well defined interoperability requirements to ensure a very high quality of the assessed air situation that will be delivered to the rest of the units. ARTAS defines, with a highly detailed level, the required performance for all components, and especially for the tracker systems which process radar data. To do this, it considers that the worst case of track performance will be expected in the case that a tracker receives only monoradar data, while other cases of fusion with extra data situations lead to relatively better performance. Therefore, the main emphasis is given to this monoradar case, leaving the definition of performance for other cases as a matter of specifying improvement factors. The most important aspect considered for tracker quality definition is the specification of track output quality in a set of well-defined representative input conditions. These conditions are classified with respect to radar and aircraft characteristics because of the very different behavior of any tracker for varying input conditions. Radar parameters represent the accuracy and quality of available data, while target conditions are the distance and orientation of the flight with respect to radar, motion state of aircraft (uniform velocity, turning, accelerating), and specific values of speed and acceleration.
Since it would not be possible to specify the performance for all possible input situations, which would require an enormous amount of figures, an area is defined in which the performance is described by a limited amount of parameters and some simple relations. Besides, since ARTAS will provide radar data processing basically for the control of civil aircraft, the specifications consider the most representative situations and the upper and lower limits of speed and accelerations in these conditions. ARTAS differentiates scenarios for two basic types of controlled areas in ATC terminal maneuvering area (TMA), covered by sensors with shorter refresh period (4 seconds), moderate range (up to 80 nautical miles or NM), and enroute area, and by sensors with longer period (12 seconds) and larger coverage (up to 230 NM). We have considered in this study the enroute area since the difficulty is higher to achieve the performance figures specified in this situation, being the design process for other situations completely similar.
Out of all possible combinations, ARTAS has carried out a choice containing the most important and realistically worst cases. It comprises a number of simple input scenarios on which the nominal track quality requirements are defined. The methodology specified for this evaluation is based on Monte Carlo simulation with the input parameters (radar and trajectory parameters) particularized for each scenario. The trajectories in different scenarios vary in the following features: (i) orientation with respect to the radar (radial or tangential starting courses, starting at a short, medium, or maximum range); (ii) sequence of different modes of flight (uniform, turns, and longitudinal accelerations); (iii) values of accelerations (upper and lower limits); (iv) values of speeds (upper and lower limits).
There are eight specified simple scenarios with uniform motion, and twelve complex scenarios including initialization with uniform motion, transition to transversal maneuver, and a second transition to come back to uniform motion. When the target is far enough from the radar, a pure radial approach to the radar leads to the worst case for transversal and heading errors during maneuver transitions, since azimuth error (much higher than radial error) is projected over these components. With a similar reasoning, a pure tangential approach is the worst case for longitudinal and groundspeed errors during maneuvers. So, the scenarios basically contain these two types of situations, varying in distance, velocities, and acceleration magnitudes. The authors have considered a couple of scenarios with longitudinal maneuvers although ARTAS does not specify performance for that type of situations. The reason for this is that these operations appear in civil operations (especially in the TMAs) and the filter is conceived to operate in real conditions. Otherwise, the resulting tracking filter could be overfitted to transversal maneuvers, but developing undesirable systematic errors with longitudinal maneuvers. The specifications for longitudinal scenarios were obtained extrapolating the ARTAS relations for the new input conditions. The resulting 22 scenarios, to be taken into account in the design of tracking filter are shown in Figure 1 (a circle represents radar position and a square the initial position of target trajectory). Since the specifications depend tightly on the input conditions, there is no a priori worst case scenario whose attainment would guarantee all cases, but all of them have to be considered simultaneously in the design process. It must be taken into account that the design of tracker will be done considering that all requirements will be met without intermediate adaptation of the tracker parameters once the tracker has been tuned for the typical radar characteristics and controlled volume (in this case, enroute area). The design will provide a single set of parameters that would allow the filter to accomplish all the specifications in all the scenarios considered.
For each of these scenarios, the performance of the tracker should approach listed performance goal values un-der the defined conditions. The accuracy requirements are expressed as a function of several input parameters depending on each specific-tested scenario: groundspeed, range, orientation of the trajectory with respect to the radar (radial and tangential projection of velocity heading), magnitude of the transversal acceleration, and magnitude of the groundspeed change. There are four quality parameters in which the requirements are defined: two for position (errors measured along and across trajectory direction, resp., longitudinal and transversal errors) and velocity (errors expressed in the groundspeed and heading components). All of them are expressed with the root mean square errors (RMSE), estimated by means of Monte Carlo simulation. Similarly, accuracy requirements are also defined for vertical coordinates, but this work will address only the 2D (horizontal) filtering, although similar ideas could be used for the design of a vertical tracker.
There are three basic parameters characterizing the desired shape of the RMS functions: peak value (RMSpv), convergence value (RMScv), and time period of RMS convergence to a certain level close to the final convergence value (RMSpv + c * RMSpv). These values are specified for different situations: initialization, transition from uniform motion to turn, and transition to come back from turn to uniform motion. Therefore, for each type of situation, the specifications are particularized according to the target evolution, defining a bounding mask for each magnitude and scenario. An example is indicated in Figure 2, with the transversal error obtained through simulation and the ARTAS bounding mask for the scenario 10. Instead of measuring performance along the whole trajectory in each scenario, only some interest points in the aircraft trajectory will be assessed to guarantee that the measured performance attains the bounding mask: convergence RMSE in rectilinear motion before and after maneuver segments (CV1 and CV2), and maximum RMSE during maneuver (PV).
The design of a tracking filter aims at attaining a satisfactory trade-off among all specifications. The quality of the design will be evaluated by means of simulation over 22 test scenarios, producing several types of trade-offs to be considered. First, the different transitions in modes of flight (uniform and maneuvers) impose a trade-off between steady-state smoothing and peak error during maneuvers, which always lead to conflicting requirements (the higher the smoothing factor the higher the filter error during transitions and vice versa). This is considered with the three representative values for each scenario and magnitude: CV1, CV2, and PV. Secondly, each one of the magnitudes evaluated (transversal, longitudinal, heading and groundspeed RMS errors) could individually shift the design towards different solutions, and so all magnitudes must be considered at the same time to arrive to a certain compromise. Finally, different design scenarios impose harder conditions for different magnitudes (radial trajectories for transversal and heading errors, etc.) so that all scenarios should be taken into account. In Table 1, we indicate the arrangement of specifications as they will be considered in the design. Specifications s(·) are particularized for the three evaluation

IMM TRACKING FILTER FOR AIR TRAFFIC CONTROL
Since the specifications for ARTAS units require a very high quality of output, the tracker in the core will have to apply advanced filtering techniques (IMM filtering, joint probabilistic data association, etc.). In this section we briefly describe the basic principles of IMM trackers, the proposed structure for these application, and the basic aspects for the design process.

General considerations
The IMM tracking methodology maintains a set of different dynamic models, each one is matched to a specific type of motion pattern, and represents the target trajectory as a series of states, with the sequence of transitions modelled as a Markov chain. In our case, the states considered will be uniform motion, transversal maneuvers (both towards right and left), and longitudinal maneuvers. To estimate the target state (location, velocity, etc.), there is a bank of Kalman filters corresponding to the different motion models in the set, complemented with an estimation of the probabilities that the target is in each one of the possible states. So, the elementary module in the tracking structure is a Kalman filter [3] which sequentially processes the measurements z[k], combining them with predictions computed according to the target dynamic model, to update the estimation of target state and associated covariance matrixx[k], P[k], respectively (see Figure 3).
The IMM maintains tracks conditioned to each jth motion state, with different Kalman filters,x j [k], P j [k], and estimation of the probability that the target is in each of them, µ j [k]. One of the basic elements in this methodology is the interacting process, which keeps all of them engaged to the most probable one. The structure considered in this work is shown in Figure 4, with four Kalman filters corresponding to the four motion states considered. It takes as input the Kalman filter , and provides the estimation of target position and kinematic state, together with estimated covariance matrix of errors, The IMM algorithm develops the following four steps to process the measures received from the available sensors to estimate the target state: intermode interaction/mixing, prediction, updating, and combination for output. . So, the input to each Kalman filter is not directly the last update but a weighted combination of all modes taking into account the mode probabilities. This step is oriented to assure that the most probable mode dominates the rest. (ii) Then, the prediction and updating phases are performed with the Kalman filter equations according to the available models for target motion contained in each mode. (iii) The estimated probabilities of modes µ j [k] are updated, based on two types of variables: a priori transition probabilities of Markov chain p i j , and mode likelihoods computed with the residuals between each plot and mode predictions Λ j [k]. (iv) Finally, mode probabilities are employed as weights to combine partial tracks for final output. Besides, each individual output and probability is internally stored to process plots coming in the future.

Design of an IMM filter
The two basic aspects involved in the design of an IMM tracking system which determine its performance are the following: the number and type of models used in the set, and transition parameters. The first aspect is dependent on each tracking problem, and we have selected, as seen in Section 3.1, a particular structure composed of four tracking modes reflecting the most representative situations in civil air traffic: constant velocity, turns to right or left, and longitudinal accelerations. They correspond to target states θ = 1, 2, 3, 4 in Figure 4. All modes interact within the IMM structure to achieve the most proper response for each situation. Mode 1, θ = 1, is a simple constant velocity model  with zero plant variance noise. Modes for tracking transversal maneuvers (turns), θ = 2, 3, are filters with circular extrapolation dynamics [4,5], one for each possible direction. They provide a highly adaptive response to transversal transitions, being one of the parameters to fix, in this filter, the typical acceleration of target when performing turns. Finally, mode θ = 4 is a linear-extrapolation motion model with a plant noise component projected along longitudinal direction. Since the target deviations along transversal direction are covered by circular modes, this last model will quickly detect and adapt to variations in longitudinal velocity during accelerations and decelerations. Each mode in the structure has its own parameters to tune, and must be adjusted in the design process. Besides, the transition probabilities between all possible pairs of modes, modelled as a Markov chain, are directly related with the rate of change from any mode to the rest. They have a very deep impact in the tracker behaviour during transitions and the "purity" of output during each type of motion, so the design must also decide the most proper values for these parameters. Since there are four modes, the transition probability matrix p i j , being defined each term as probability of the target arriving to state j at time k, given that the state at time k − 1 was i, is The number of parameters have been simplified by considering only as possible transitions between uniform motion and the rest of modes. The parameters p UT , p UL are the probabilities of starting transversal and longitudinal maneuvers, given an aircraft at uniform motion, while the parameters p TU , p LU are the probabilities of transitions to uniform motion, given that the aircraft is performing, respectively, transversal and longitudinal maneuvers.
It is important to notice that all parameters, those in each particular model plus transition probabilities in Markov chain, are completely coupled through the IMM algorithm since partial outputs from each mode are combined and feedback all modes. So, there is a strongly nonlinear interaction between them, making the adjusting process certainly difficult. The whole set of parameters in the tracking structure is summarized in Table 2.

DESIGN OF FILTER PARAMETERS
The design of the particular IMM tracking structure addressed in this work, stated as adjusting the seven numeric input parameters to fit filter performance within ARTAS specifications, can be generally considered as a numerical optimization problem. We are searching for the proper combination of real input parameters that minimizes a real function assessing the quality of solutions as a cost f : The subspace V stands for the region of feasible solutions, defined as those vectors representing a valid IMM filter: parameters for probabilities must fall in the interval [0, 1] and parameters for variances must be positive. These are the only constraints to be accomplished by solutions during the search. Performance specifications are not considered as constraints here, but they will be used as penalty terms in the objective cost function. The cost would achieve a minimum value of zero only in the ideal case of a solution accomplishing all specifications, grading the rest of possible cases with a positive global cost function that will be detailed later.

Evolution strategies
In numeric optimization problems, when f is a smooth, low-dimensional function, there are an available number of classic optimization methods. The best case is for lowdimensional analytical functions, where solutions can be analytically determined or found with simple sampling methods. If partial derivatives of function with respect to input parameters are available, gradient-descent methods could be used to find the directions leading to a minimum. However, these gradient-descent methods quickly converge and stop at local minima, so additional steps must be added to find the global minimum. For instance, with a moderated number of global minima, we could run several gradient-descent solvers to find the best solution. The problem is that the number of similar local minima increases exponentially with dimensionality, making these types of solvers unfeasible. In our particular case, besides a high-dimensional input space causing multimodal dependence, we do not have an analytical function to optimize. It is the result of a complex and exhaustive evaluation process implying the simulation and performance assessment of tracking structure on the whole set of 22 scenarios defined. The evaluation of a single point in the input space requires several minutes of CPU time (Pentium III, 700 MHz). Besides, the evaluation of quality after all simulations is not direct but it should take into account system performance in all scenarios and magnitudes in comparison with the whole table of specifications. As we will see later, multiple specifications (or objectives) will increase the number of solutions with similar performance, increasing therefore the complexity of the search.
For complex domains, evolutionary algorithms have proven to be robust and efficient stochastic optimization methods, combining properties of volume and path-oriented searching techniques. ES [6] are the evolutionary algorithms specifically conceived for numerical optimization, and have been successfully applied to engineering optimization problems with real-valued vector representations [7]. They combine a search process which randomly scans the feasible region (exploration) and local optimization along certain paths (exploitation), achieving very acceptable rates of robustness and efficiency. Each solution to the problem is defined as an individual in a population, codifying each individual with a couple of real-valued vectors: the searched parameters and a standard deviation of each parameter used in the search process. In this specific problem, one individual will represent the set of dynamic parameters in the IMM structure, as indicated in Table 2, (x 1 , . . . , x 7 ), and their corresponding standard deviations (σ 1 , . . . , σ 7 ).
The optimization search basically consists in evolving a population of individuals in order to find better solutions. The computational procedure of ES can be summarized in the following steps, according to the named "µ + λ" strategy defined by Bäck and Schwefel [8], and particularized for our problem: (1) generate an initial population with µ individuals uniformly distributed on the search space V ; (2) evaluate the objective value for each individual in population f ( − → x i ), i = 1, . . . , µ; (3) Select the best parents in population to generate a set of λ new individuals, by means of genetic operators of recombination and mutation. In this case, recombination follows a canonical discrete recombination [6], and mutation is carried out as follows: where x i and σ i are the mutated values and N(0, σ) stands for a normal distribution with zero mean and variance σ 2 ; (4) calculate the objective value of the generated offspring f ( − → x i ), i = 1, . . . , λ, and select the best µ individuals of this new set containing parents and children to form the next generation; (5) Stop if the halting criterion is satisfied. Otherwise, go to step (3).
We have implemented ES for this problem with a size of 50 + 30 individuals and mutation factor ∆σ = 0.9. The fitness function will directly depend on the differences between RMS values of errors, evaluated through Monte Carlo simulation, and ARTAS specifications for all scenarios and magnitudes, as will be detailed next. It is important to notice that simulations are carried out using common random numbers to evaluate all individuals in all generations, enhancing system comparison within the optimization loop. In other words, the noise samples used to simulate all scenarios in the RMS evaluation are the same for each individual in order to exploit the advantages coming from the use of a deterministic fitness function. Besides, the number of iterations was selected to guarantee that confidence intervals of estimated figures were short in relation to the estimated values.
A basic aspect to achieve successful optimization in any evolutionary algorithm is the control of diversity, but this appropriateness will depend on the problem landscape. If a population converges to a particular point in a search space too fast in relation to the roughness of its landscape, it is very probable that it will end in a local minimum. On the contrary, a too slow convergence will require a large computational effort to find the solution. ES give the higher importance to the mutation operator, achieving the interesting property of being "self-adaptive" in the sizes of steps carried out during mutation, as indicated in step (3) of the algorithm above. Before selecting an algorithm for optimization, it is interesting to consider the point of view of the "no free lunch" (NFL) theorem [9], which asserts that no optimization procedure is better than a random search if the performance measurement consists in averaging arbitrary fitness functions. The performance of ES has been widely analyzed under a set of well-known test functions [8,10]. They are artificial analytical functions used as benchmarks for comparison of representative properties of optimization techniques, such as convergence velocity under unimodal landscapes, robustness with multimodality, nonlinearity, constraints, presence of flat plateaus at different heights, and so forth. However, the performance on these test functions cannot be directly extrapolated to real engineering applications. The application of ES to a new problem, such as our complex IMM design against multiple specifications where the landscape properties are not known (it is not known even if there is a global minimum or not), is a challenge open to research.

Multiobjective optimization
The selection of the proper fitness function for this application is the problem-dependent feature with the highest impact on the algorithm (higher than the ES parameters such as population size or mutation factor). Really, we should regard this design as a multiobjective optimization problem, where each individual objective is the minimization of difference between desired specification and assessed performance in each specific figure of merit. When a problem involves simultaneous optimization of multiple, usually conflicting objectives (or criteria), the goal is not so clear as in the case of single-objective optimization. The presence of different objectives generates a set of alternative solutions, defined as Pareto-optimal solutions [11]. The presence of conflicting multiple objectives leads to the fact that different solutions cannot be directly compared and ranked to determine the best one, but the concept of domination appears for com- x 2 is better than − → x 1 simultaneously in all objective functions considered. In any other case, they could not be strictly compared. Taking into account this concept of domination, a Pareto-optimal set P is defined as the set of solutions such that there exists no solution in the search space dominating any member in P.
Some multiobjective optimization techniques have the double goal of guiding the search towards the global Paretooptimal set and at the same time covering as many solutions as possible. There are several proposed evolutionary methods [12] that address this goal by maintaining a population diversity to cover the whole Pareto front. This fact implies first the enlargement of population size and then specific procedures to guarantee guiding the search to the desired optimal set with a well-distributed sample of the front. Among these procedures, we can mention methods, such as selection by aggregation and so forth, switching the objectives during the selection phase to decide which individuals will appear in the mating pool. Zitzler et al. [12] analyze and compare, over some standard test analytical functions, some of the most outstanding multiobjective evolutionary algorithms.
From the authors point of view, the peculiarities of the problem dealt with, namely, the complexity and computational cost of evaluation function together with the considerable number of specifications, preclude the application of techniques to derive the whole Pareto set. We have considered a weighting sum on partial goals to build a global fitness function: As indicated by Deb [11], this type of approaches with weighted sums converge to particular solutions of Pareto front, corresponding to the tangential point in the direction defined by the vector of weights. The general idea is illustrated in Figure 5  However, a large number of specifications will make the weighted summation cumbersome, being difficult that all objectives are simultaneously considered to guide the search. In our specific problem, we should fix a weighting vector with 264 components. A variation is proposed to reduce the number of objectives in the sum by exploiting knowledge about the problem. Basically, objectives with similar behavior are grouped to select a "representative" per group, the one with the worst value, so that it guarantees that all objectives in the group are represented in the final function. If we consider Table 1 with the whole set of specifications, we are going to select the worst case for each column, leaving only 12 terms in the summation. It is important to notice that this maximum operation will break the linearity of function with respect to objectives and will make the landscape depend on each specific input vector. A trajectory of solutions in the search process may jump along different goal functions if the scenarios with the worst case change. The justification comes from the fact that each magnitude has certain dependence with the input parameters similar in all scenarios, so a single representative is enough to be considered in the optimization. Besides, the selection of the worst case assures that if the method can satisfy that term, all the scenarios will be simultaneously accomplished.
Taking into account this consideration, the fitness function, which assesses the quality of a solution as the degree of attainment performance figures with respect to specifications, is presented next. The following details have also been considered.
(i) It assesses the excess over the specification for each performance figure, penalizing a solution as the error increases, but once the error is below the specification, the cost is zero. This is so because there is no additional advantage if the RMSE decreases more after the required values are attained. This is implemented for each magnitude by means of the expression R(p i − s(p i )), where p i is the ith performance figure (RMSE), s(p i ) the specification, and R(·) the ramp function: and so are normalized with the specification value, defining a partial cost for ith figure, (iii) In order to add some flexibility in the trade-off between maneuver and uniform motion performances, weighting factors α t are included. They allow us to vary the priority of these performance figures, in the case where all of them cannot be attained at the same time, defining therefore a cost per jth scenario, where the subindex i represents each interest magnitude (longitudinal, transversal, groundspeed, and heading) and j the scenario index. instant (PV, CV1, and CV2). Therefore, the final goal function to be minimized is as follows: So, this function considers the relative excesses over specifications for all performance figures, each one assessed in the worst case scenario.

RESULTS
In this section, the results obtained along the optimization process to adjust the filter parameters according to ARTAS specifications are presented and analyzed. They have been obtained particularizing expression (6) to the case of a weight of 1 for all magnitudes α PV = α CV1 = α CV2 = 1. First, Figure 6 summarizes the evolution of best individual in the population (the one with the lowest value of fitness), indicating graphically the accomplishment of specifications along the generations. Each design objective is presented by a row in the diagram, while the best individual for each generation appears in each column. The grey level of position (i, j) in the image indicates the quality of the fitting to the ith specification of the best individual for the jth generation. The grey level represents linearly the relative excess over the restriction (no excess is presented as white, 100% or higher excess as black), which is the partial cost function related with this constraint. Therefore, a completely white column means that the optimization process has found a set of parameters able to fulfil all design restrictions, while a complete white row means that all best individuals in this optimization exercise are able to fulfil the specification for this magnitude and situation. Below, the fitness function computed from the whole set of partial costs as indicated in Section 4 is plotted. This kind of figure serves not only to see the convergence of the optimization process graphically, but also to see the most demanding performance criteria to be accomplished and to compare the suitability of a predefined tracking scheme (with some free design parameters) for a certain tracking problem. Applying exactly the same proposed methodology, we could have performed the optimization exercise with an alternative IMM structure, or even with a different tracking technique with open design parameters, and compared after designing the process its capabilities against specifications.
As it can be seen, the optimization process makes the overall figure lighter from the initial generations (left) to the end of the optimization (right), achieving a trade-off point to accomplish as many specifications as possible. The highest improvement is carried out in the first 80 generations, with very slight modifications from that point until the end. The rows with a darker profile indicate higher difficulty to attain that specification together with the rest. So, scenarios 12 and 13, corresponding to specifications 133-156, present the worst performance after the optimization. The specific performance values and ARTAS bounding masks for these scenarios, corresponding to transversal maneuvers at 215 NM, v = 300 m/2, a = 2.5 m/s 2 , (scenario 12) and at 65 NM, v = 150 m/2, a = 6.0 m/s 2 , (scenario 13), are indicated in Figures 7 and 8. The magnitudes with worst performance are the transversal and heading errors (peak values) during transversal maneuvers. The peak value of heading error is the globally worst figure in the set, more than 100% over specification. Besides, as it can be seen, the convergence error values for some of the magnitudes in these scenarios are practically tangent to specifications, indicating that the optimization process has effectively considered all of them to arrive to the final trade-off solution. So, this method selects the parameters adapting system behavior to the bounding mask. This is apparent not only for the 1 Figure 9: Evolution of fitness and performance in each specific objective.
presented scenarios with worst cases but for all design scenarios as well.
Different runs of the global optimization process (using different random seeds to generate individuals in the initial population) were carried out to analyze the consistency of the solutions obtained. The results of ten independent runs are indicated in Figure 9, presenting only the best individual in population after optimization (instead of the whole evolution process) and the final values of fitness achieved.
As it can be seen, different runs led to solutions quite consistent in terms of overall fitness and whose specifications are presenting problems to the filter (always those in scenarios 12 and 13). However, the specific vector solutions found after optimization in each run had significant differences, indicating that fitness function probably has a multimodal landscape, even after having selected a particular set of weighting factors among specifications, α = 1.
Since it is not possible to represent fitness landscape with seven dimensions, the following analysis was carried out. The three solutions with closest fitness values, resulting from runs 1, 2, and 5, were selected to be combined and to generate a grid of linear combinations (convex hull) as follows: The fitness landscape for a grid with α, β varying in the interval [−1.5, 0.5], in steps 0.1 units, is indicated in Figure 10. It (really, there are much more solutions with similar fitness values that are presented in Figure 10; they are only particular cases of linear combinations). The particular solutions 1, 2, and 5 correspond to the points (0, 0), (1, 0), and (0, 1) in plane. In Figure 11, the projection on αβ plane is presented with grey levels, where the feasible region in αβ plane can be clearly separated. All solutions within the convex combinations of solutions α, β in [0, 1] are feasible. Besides, the 2D graphs corresponding to the paths connecting all pairs of solutions are presented in Figures 11 and 12. As it can be seen, the solutions found by algorithm are effectively local minima of fitness function in spite of the fact that the function is almost flat in this region of convex linear combinations. This shows that the algorithm is capable of finding appropriate solutions, and confirms the fact that we have a multimodal function even after having combined the multiple restrictions in a scalar function. Different runs arrived to different local minima in a region where the relative difference between minima can be practically neglected, so all solutions can be taken as good design points for the adopted criterion. The algorithm was carried out with different criteria (for instance, the penalty of RMSE peak values being ten times higher than convergence RMSE), achieving results consistent with the preferences: all specifications with the highest priority were first accomplished, leading to higher errors in the other specifications.

CONCLUSION
In this paper, we have described a methodology based on ES for the design of IMM-tracker techniques to accomplish a considerably large set of predefined specifications. An exhaustive set of test scenarios with performance specifications for each and a specific IMM structure with open parameters are the input to solver. The procedure may be summarized as performing an optimization over the pa-rameters space, using ES, defining as the fitness function one combination of partial excesses over specifications that takes into account some knowledge about the problem in the form described in Section 4. This fitness function summarizes the attainment of all interest accuracy statistics for the different interest times (steady state, start and end of maneuvers, etc.) in all design scenarios. The evaluation involved the costly Monte Carlo simulation, as specified by ARTAS, to calculate accuracy statistics, although the methodology is open for the inclusion of other possible evaluation methods for IMM tracking filters, such as the one described in [9]. This method has been successfully used in a monoradar application, leading to a significant improvement over previous nonsystematic approaches for the same problem. Even more, the form of fitness function described serves as a method for relaxing constraints: those more important for us are provided a higher weight in (6), and those not so important a lower weight.