Using the Decomposition-Based Multi-Objective Evolutionary Algorithm with Adaptive Neighborhood Sizes and Dynamic Constraint Strategies to Retrieve Atmospheric Ducts.

The traditional method of retrieving atmospheric ducts is to use the special sensor of weather balloons or rocket soundings to obtain information intelligently, and it is very expensive. Today, with the development of technology, it is very convenient to retrieve the atmospheric ducts from Global Navigation Satellite System (GNSS) phase delay and propagation loss observation data, and then the GNSS receiver on the ground forms an automatic receiving sensor. This paper proposes a hybrid decomposition-based multi-objective evolutionary algorithm with adaptive neighborhood sizes (EN-MOEA/ACD-NS), which dynamically imposes some constraints on the objectives. The decomposition-based multi-objective evolutionary algorithm (MOEA/D) updates the solutions through neighboring objectives, the number of which affects the quality of the optimal solution. Properly constraining the optimization objectives can effectively balance the diversity and convergence of the population. The experimental results from the Congress on Evolutionary Computation (CEC) 2009 on test instances with hypervolume (HV), inverted generational distance (IGD), and average Hausdorff distance ∆2 metrics show that the new method performs similarly to the evolutionary algorithm MOEA/ACD-NS, which considers only the dynamic change of the neighborhood sizes. The improved algorithm is applied to the practical problem of jointly retrieving atmospheric ducts with GNSS signals, and its performance further demonstrates its feasibility and practicability.


Introduction
In both technical practice and theoretical research, there are a large number of multi-objective optimization (MOP) problems that need to be optimized synchronously for multiple objectives [1][2][3]. When solving multi-objective optimization problems, the traditional method is to convert this problem into a single-objective problem based on prior knowledge and to use single-objective optimization methods to obtain a satisfactory Pareto optimal solution [4]. The Pareto solution contains a lot of information, which can be used to analyze the correlation between multiple objectives; furthermore, perform the search along the Pareto set, even for high-dimensional models but are of local nature [21]. However, there are still some defects in their algorithm. If the evolution regions are too large, a single new solution may replace several old solutions and further diminish the population diversity. Based on the ENS-MOEA/D framework, and in order to balance the population diversity and convergence, this paper combines the constrained decomposition approach with the adaptive strategy proposed by Wang et al. to reduce the improvement regions and control the evolution direction [22] and apply it to solve the problem of the joint inversion of atmospheric ducts based on global navigation satellite system (GNSS) signals.
The propagation of an electromagnetic wave in the atmospheric boundary layer, especially in the surface layer, is affected by atmospheric refraction, and the propagation path will bend to the ground. Under certain meteorological conditions, the curvature of the path may exceed the curvature of the Earth's surface, and the electromagnetic wave will be partially trapped in a thin layer of atmosphere with a certain thickness. This phenomenon is called the atmospheric ducts [23]. The atmospheric ducts can capture an electromagnetic wave and change its propagation path, so that the electromagnetic wave can be propagated beyond the horizontal sight distance of the electromagnetic system, realizing the over sight distance detection and over sight distance reception. Therefore, it is very important to study the atmospheric ducts for electromagnetic communication, and their influence on the propagation of electromagnetic wave has a broad application prospect. So far, with the development of computers, more and more intelligent algorithms are being applied in the field of meteorology [24][25][26][27][28][29]. The traditional method of retrieving atmospheric ducts is to use the special sensor of weather balloons or sounding rockets to obtain information intelligently; however, weather balloons and sounding rockets are expensive in terms of expendables and manpower. Today, with the development of technology, it is very convenient to obtain atmospheric ducts from GNSS phase delay and propagation loss observation data [30], so we tried to retrieve the atmospheric ducts by placing a receiver on the ground to receive the signal from the GNSS and use a new algorithm to analyze it. This method will would save a lot of costs in the retrieval of atmospheric ducts.
In this paper, we propose a hybrid decomposition-based multi-objective evolutionary algorithm with adaptive neighborhood sizes (EN-MOEA/ACD-NS), which dynamically imposes some constraints on the objectives, and compare the improved algorithm with the MOEA/D, MOEA/ACD, MOEA/ACD-NS, and EN-MOEA/D on the MOPs and unconstrained function (UF) problems in the Congress on Evolutionary Computation (CEC) 2009 standard test instances [31]. To demonstrate the practical application of our algorithm and solve the practical problem of jointly retrieving atmospheric ducts, we applied the algorithm to solve the problem of the joint retrieval of atmospheric ducts based on GNSS signals [32]. This new algorithm will provide a new method and save a lot of costs for the retrieval of atmospheric ducts, which is of great significance for the research of atmospheric ducts.

Definition of Multi-Objective Optimization
Some problems that we encounter in life have only one objective function and are known as single-objective problems, such as the classic traveling salesman problem (TSP) problem. Some problems, called multi-objective problems, have more than two objective functions whose goals are contradictory, such as logic circuit design, and it is difficult for each function to achieve optimal solutions at the same time [33]. A problem with m optimization goals and n decision variables can be expressed as follows [32,34]: Here, → x ∈ R n represents the n-dimensional decision space and f i → x is the ith objective. The k inequality constraints and l equality constraints are the conditions constraining the above formula, which satisfy the following: The objectives cannot be optimized synchronously. In order to achieve an optimal solution of the function F → x , the objectives that conflict with each other are considered comprehensively to meet the constraint conditions and the objective functions are optimized. Generally, the solution that satisfies the constraints is defined as a feasible solution. In order to meet the optimal conditions, Pareto dominance is introduced to describe the relationship between two feasible solutions, and then the Pareto optimal solution is found. Suppose there are two feasible solutions, x a is called the Pareto optimal dominating solution. The set composed of such solutions is the Pareto optimal solution set (PS), and the edge of the target space mapped by PS is called the Pareto front (PF) [4].

Evaluation Metrics
It is a very important and difficult task to evaluate the performance of MOEAs. At present, evaluation focuses mainly on the following three aspects: (a) the approximation degree between the solution set and the Pareto optimal surface (convergence), (b) the distribution uniformity of the solution set in the target space (distribution uniformity), and (c) the extensive degree of the solution set (widespread distribution) [14,31,35]. In this paper, the performance of the improved algorithm in different test problems is evaluated by HV, inverted generational distance (IGD), and the average Hausdorff distance (∆ 2 ).

Hypervolume
HV was first proposed by Zitzler in 1998, and is called a set of quality test indicators of the "space coverage area" [36]. It is widely used for quantitative measurement of population quality because HV is still sensitive to Pareto dominance and population diversity in high dimensions and does not decrease with increasing dimensions. In general, HV is used to measure the hypervolume between non-dominated points and reference points, which can be expressed as follows [36]: Here, p * = p 1 , p 2 , . . . , p n , p * is the number of reference points, and vol is the Lebesgue metric; the larger the value of HV, the better the performance of the algorithm.

Inverted Generational Distance
IGD is a comprehensive performance evaluation index. It mainly evaluates the convergence and distribution performance of the algorithm by calculating the sum of the minimum distance between each point (individual) on the real Pareto front and the group of individuals obtained by the algorithm. It can be expressed as follows [37]: Here, p * = p 1 , p 2 , . . . , p n , p * is the number of reference points, and Q is the optimal Pareto solution set obtained by the algorithm. D (v, Q) is the minimum Euclidean distance between individual v and group Q in P * ; the smaller the value of IGD, the better the performance of the algorithm.

Average Hausdorff Distance
The classical Hausdorff distance is a measurement criterion of the maximum and minimum distance, and is used to describe the similarity between two sets of points. However, it is susceptible to external interference, resulting in large errors in the final results. As a supplementary measure, Schütze et al. revised the original generational distance (GD) and IGD and proposed an average Hausdorff distance based on these two indicators [12]. Assuming that there is a set of non-dominated points X = {x 1 , x 2 , . . . , x n } and a set of reference points P * = p 1 , p 2 , . . . , p m which is a finite non-empty set. Then, the Euclidean norm of the average Hausdorff distance is defined as follows [12]: The smaller the value of ∆ 2 , the better the performance of the algorithm.

Basic Framework
The main idea of MOEA/D is to use an aggregate function to decompose a MOP problem into multiple single-objective optimization problems. The individual uses the information of adjacent solutions to evolve in order to obtain the optimal solution. At present, the improvement of MOEA/D is mainly concerned with improvements in the weight vector, decomposition method, matching selection, evolutionary operator, replacement and selection method, and computing resource allocation [38][39][40][41][42][43][44]. Various improvement methods have been carried out according to the basic framework of MOEA/D. The algorithm in this paper improves the decomposition method and matching selection. In the basic framework, some information needs to be saved each time, such as the solution set x 1 , x 2 , . . . , x p in which population size is p, the corresponding fitness set FV 1 , FV 2 , . . . , FV p , the current optimal solution set b 1 , b 2 , . . . , b p , and the external population EP. Here, x i (I = 1, 2, . . . , p) is the current solution to the ith objective. The basic framework can be expressed as follows [31,39]: Step 1: Initiation.
Step 1.1: Let EP = Ø Step 1.2: Calculate the Euclidean distance between two weight vectors according to which we select T individuals as a neighbor population. Set B(i) = {i 1 , i 2 , . . . , i T } for i = 1,2,...,n, and λ j as the neighbor vector to λ i ∀j ∈ B(i).
Step 1.3: Generate an initial population randomly.
Step 1.4: Step 1.5: according to the different problems.
Step 2.1: After generating new generation offspring, select two individuals x l and x k from the adjacent set B(i), and generate new children x new by a genetic algorithm or differential evolution [45].
Step 2.2: Adjust the x new by using a modified heuristic method for the particular problem.
Step 2.3: Stop loop if the stop condition is satisfied and output EP; otherwise, jump to Step 2.

Adaptive Neighborhood Sizes
The new solution generated by MOEA/D involves a number of neighboring objectives. The larger the number, the more likely the optimal individual is to be globally optimal. Otherwise, the smaller the number, the more likely it is to be locally optimal. Therefore, the determination of the number of neighboring problems (T) is related to the performance of the algorithm. Zhao et al. proposed a method to dynamically adjust the number of neighboring objectives [15]. This paper uses their methods to improve the original MOEA/D. First, we select several number values as candidates for the neighborhood sizes, such as {2, 4, 8, 12, 16}. The initial value is set to 2, and, when the number of iterations reaches a certain limit, such as 50, it is determined by probability p k,G whether the T is redefined. p k,G is expressed as follows: where S k,G represents the probability of successfully improving progeny when the neighboring number is k in the learning period (LP). TNS k,g is the number of new individuals generated by the children determined by the kth neighborhood size (NS). TNS_success k,g is the number of individuals that have been successfully improved from TNS k,g . Here, δ = 0.05 is a small constant used to avoid S k,G becoming zero.

Adaptive Constraints Approach
In the previous section, dynamically adjusting NS only affects the number of fathers that may produce new individuals. However, how to select right individuals from parents to mate more quickly and efficiently is still a question that needs to be solved. Using the method of Reference [22], we could dynamically adjust the direction of paternal mating to further improve MOEA/D. After the calculations of the evolutionary algorithm, many new individuals will be generated. In order to reduce the distribution space of new individuals, the optimization direction should be limited. For example, in the ith objective, it must satisfy the following: where the constraint condition is is the aggregated objective function of each optimization objective determined by the weight vector λ i ; z * is a utopian point; a i , F(x) − z * is the acute angle between a i and F(x) − z * , which represents the basic optimized direction; and a i is the optimized direction vector of the objective. Additionally, let θ i be the control parameter that determines the size and direction of the optimizable region. Consequently, dynamically adjusting θ i means controlling the optimized direction and optimizable population.
The specific dynamic strategy is as follows. Assuming that x i is the solution set of the current objective, and defining the divergence function α( At this point, the number of neighbors will increase. When all α(x i ) are small, the solution set x is close to the real value and the population has good diversity. The method of Reference [22] also provides another measurement: the relative divergence degree α r (i), which is defined as α r (x i ) = [n · α(x i )]/ n k=1 α(x k ). Under such constraints, α(x i ) can quickly reach zero when θ i is small. However, the objective functions may not be optimal at this time. Therefore, we could adjust θ i according to the divergence degree. If α(x i ) is very large, the optimizable region can be narrowed; namely, by reducing θ i . This can be expressed as follows: where θ i (t) is the control parameter of the t generation in the ith objective. For convenience, the initial value of θ i is set to θa i , a ji,min , where j is the adjacent problem of the ith objective. θ i,max is the double maximum angle between a i and m coordinates.

Test Results
In order to test the performance of the new hybrid algorithm, two types of test function were selected to compare the performances of the original MOEA/D, MOEA/ACD, MOEA/ ACD-NS, EN-MOEA/D, and EN-MOEA/ACD-NS. One type of test is the MOP test function set in CEC 2009 [31]. The other type of test is the unconstrained function (UF) problem, whose Pareto set is easy to be obtained [46]. A set of unconstrained functions proposed by Van Veldhuizen in 1999 can comprehensively test a multi-objective optimization algorithm, and is in accordance with the design criteria of the MOP test function proposed by Whitley et al. [46,47]. The evolutionary algorithm based on MOEA/D involves three kinds of decomposition methods, namely: the weighted sum approach (WS), the weighted Tchebycheff approach (WT), and the penalty-based boundary intersection approach (PBI). In this paper, we selected the PBI approach and generated new individuals by differential evolution. The mutation operator is a polynomial variation, and the weight vector is generated by the old method.

Results and Discussion of the Classical Test Functions
In this paper, we used test sets from Reference [48,49], which are UF1-UF9 and MOP1-MOP7, respectively. The search space dimensions of MOP1-MOP6 and UF1-UF6 are 2, and the search space dimensions of MOP7 and UF7-UF9 are 3. For the purpose of fairness, the recombination operators and decomposition methods of all of the algorithms in the comparative experiment are the same. In the programming implementation, we used a differential evolution operator [50] and polynomial mutation operator [48] as the recombination operators; furthermore, the method of the Chebyshev decomposition was used to decompose the multi-objective function. Considering the size of the population and the parameter of the neighbors' number in ENS-MOEAD, we set the parameters NS in ENS-MOEAD, MOEA/ACD-NS, and EN-MOEA/ACD-NS to {10, 20, 25}.Through the experiment of independent parameters, we found that the range of θ i within Reference [2,10] is better for the average experiment of the test sets UF1-UF9 and MOP1-MOP7. In this experiment, we set the value of θ i as   Figures 1 and 2 show box plots of the metrics. The box plots contain the following of five values: maximum, minimum, the first quartile, the third quartile, and median. They show the shape, center, and dispersion of the data distribution. The "whiskers" above and below the rectangle show the positions of the minima and maxima, respectively, and the "+" in the figures indicate the outliers. In the MOPs, with HV and IGD as the evaluation indicators, the outliers of EN-MOEA/ACD-NS are significantly less than the other algorithms, indicating that our improved algorithm is more robust. Observing the quartiles and medians of the box plots, it is easy to see that EN-MOEA/ACD-NS is generally consistent with MOEA/ACD-NS except for MOP2 and MOP3 with IGD as the metric and the IGD of MOEA/ACD-NS is less than other algorithms, indicating that the two algorithms perform similarly in the calculation process and have a better accuracy. In the UFs, although the quartile and median of ENS-MOEA/ACD-NS differ from MOEA/ACD-NS, there are relatively few outliers in the various metrics of the hybrid algorithm, which further confirms the robustness of our algorithm. If IGD is taken as the metric, the median of EN-MOEA/ACD-NS is relatively small compared with the other algorithms, except for the problem UF10. Similarly, with HV as an indicator, the median of EN-MOEA/ACD-NS is larger than the other algorithms except for UF10. This result indicates that the EN-MOEA/ACD-NS has a better convergence and distribution. So, we can conclude that the hybrid algorithm performs better, as mentioned in the analysis of Tables 1-3.

Results and Discussion of the Joint Inversion of Atmospheric Ducts Problem
Today, with the development of technology, we can obtain the meteorological data by GNSS. In this paper, we placed a receiver on the ground and adjusted its antenna to a proper angle in order to receive the signal from the GNSS, and then recorded the GNSS phase delay and propagation loss observation in turn, which formed an automatic receiving sensor [30]. Considering the complexity of the atmospheric boundary layer, increasing the number of GNSS antenna sensors can effectively increase the information capacity of the atmospheric ducts, and then enhance the constraints of the retrieval of the atmospheric ducts' parameters. Theoretically, it can improve the retrieval accuracy of the atmospheric ducts; therefore, we proposed the idea of obtaining GNSS signals from multiple sensors on the ground in order to jointly retrieve the atmospheric ducts. After we collected the data, we used the new algorithm to solve the practical problem of jointly retrieving atmospheric ducts. The specific framework of the joint retrieval of atmospheric ducts has been elaborated in our previous work [32,51]. This paper takes the problem as a series of practical test problems to verify whether the improved algorithm can be applied in practical conditions. According to the classification of this problem in Reference [32], fifteen types of problems in three situations were used. The maximum number of iterations was set to 50, and the population number was 300. The learning period (LP) was defined as 20. We ran the algorithm 20 times independently and analyzed the results statistically. Table 4 shows the simulated values of the retrieved parameters, as well as the corresponding mean values and standard deviations. In order to show the results for diversity, we calculated the evaluation metric ∆ 2 . Comparing the metric values, it is easy to find that the gap between the simulated and retrieved parameters gradually became larger when the number of retrieved parameters increased under the same known conditions. When the number of inversion parameters is the same, the deviation of individual parameters will increase when the known conditions increase. This is also shown by the metric ∆ 2 . For example, when considering only one antenna to receive data, all ∆ 2 were less than 2, whereas, when three or four antennas were considered, the ∆ 2 in D = 5 and D = 6 were greater than 3. We can conclude that the convergence of the approximate solution set becomes weaker when the discretization degree of the decision space increases. Furthermore, the optimal solution obtained does not necessarily converge to the real value. In fact, similar conclusions can be reached from the results shown in Reference [51]. Therefore, we can conclude that increasing the known information may not effectively improve the retrieved effect in the actual operation process when the number of inversion parameters is certain.  Figure 3 shows the profiles retrieved and the corresponding absolute deviations versus heights in various problems. GPS1-5 represents the type of problem, and, in Table 4, M and D are the dimensions of the solution space and decision space, respectively. In GPS1, the retrieved profile obtained converges closely to the simulated profile, and the absolute deviation is less than 0.6N units. It could be claimed that the improved algorithm achieved better results. However, the absolute deviation is around 15N units or even 20N units in other cases, suggesting that the retrieved effect does not seem to be improved when the known information increases. After a careful comparison of the absolute deviation under the same D and different M conditions, we find that the retrieved parameters approximate to the real value. When the number of retrieved parameters D increased, the absolute deviation gradually increases between 100 and 400 m. The deviation between 400 and 600 m was slightly less than the deviation between 100 and 400 m in the same problem.  In order to better describe whether an increase in information content is effective for improving the retrieved results, we compared results for the retrieval of four parameters. As Figure 4 shows, the retrieved result is the best when the number of antennas is 1, and the worst when the number is 2.
Moreover, the inversion effect does not seem to be improved as the number of antennas increases, which is consistent with the analysis of Figure 3.

Conclusions
To solve multi-objective optimization problems, we improved the original MOEA/D algorithm by dynamically adjusting the neighborhood sizes and improvement regions and proposed EN-MOEA/ACD-NS. The larger NS can enhance the global search ability. In order to properly adjust the scope of NS, we chose p k,G as the parameter to judge whether to change the NS. At the same time, the expansion of NS is bound to affect the convergence of the population. We can control parameter θ i and dynamically adjust the convergence direction to improve the convergence accuracy and speed of our algorithm. To test its effectiveness and feasibility, the CEC 2009 standard test instances and the practical problem of jointly retrieving atmospheric ducts were used for numerical simulation experiments. The experimental results demonstrate that EN-MOEA/ACD-NS has a better computational performance, which is close to that of MOEA/ACD-NS, and it has a better convergence and distribution. When retrieving the refractivity parameters only, the improved algorithm can obtain a better retrieved value with the absolute deviation within 3N units. When the number of retrieved parameters increases, the value of ∆ 2 becomes larger and the diversity of the solution sets increases, which means that the difficulty of optimal solutions converging to the real values also increases. This new algorithm will provide a new method for the retrieval of atmospheric ducts, which is of great significance for the research of atmospheric ducts.
The performance of the proposed EN-MOEA/ACD-NS is better than other algorithms in the application of retrieving the atmospheric ducts. However, the NS and θ i need to be adjusted by experiments according to practical problems. When the parameter evaluation for a real multi-objective problem is complicated, this algorithm will consume a lot of time. In the future, we will make a more detailed analysis on the parameter setting of EN-MOEA/ACD-NS and provide the reference parameter setting intervals under different decision space dimensions and search space dimensions to solve this problem.