Abstract
The weighted sum and genetic algorithm-based hybrid method (WSGA-based HM), which has been applied to multiobjective orbit optimizations, is negatively influenced by human factors through the artificial choice of the weight coefficients in weighted sum method and the slow convergence of GA. To address these two problems, a cluster and principal component analysis-based optimization method (CPC-based OM) is proposed, in which many candidate orbits are gradually randomly generated until the optimal orbit is obtained using a data mining method, that is, cluster analysis based on principal components. Then, the second cluster analysis of the orbital elements is introduced into CPC-based OM to improve the convergence, developing a novel double cluster and principal component analysis-based optimization method (DCPC-based OM). In DCPC-based OM, the cluster analysis based on principal components has the advantage of reducing the human influences, and the cluster analysis based on six orbital elements can reduce the search space to effectively accelerate convergence. The test results from a multiobjective numerical benchmark function and the orbit design results of an Earth observation satellite show that DCPC-based OM converges more efficiently than WSGA-based HM. And DCPC-based OM, to some degree, reduces the influence of human factors presented in WSGA-based HM.
1. Introduction
Earth observation satellites provide essential information on ocean, land, and atmosphere, which are very important in the environment protection and resources management. The first step of satellite mission design is usually the determination of a suitable orbit. The objective of orbit design for Earth observation satellites is to ensure that all target sites are best visited, including observation sites and ground stations. The quality of an orbit can be measured with key orbit performance indices [1]. The key orbit performance indices of an Earth observation satellite include the total coverage time, the frequency of coverage, the average time per coverage, the maximum coverage gap, the minimum coverage gap, and the average coverage gap [1, 2]. Thus, orbit design is a typical multiobjective optimization problem. Numerical methods for multiobjective orbit design optimization can be classified into three primary groups: indirect methods, direct methods, and evolutionary algorithms [3]. The last group is currently receiving research attention because of the capability of achieving global optima in very large search spaces.
In the evolutionary optimization for multiobjective orbit design, multiobjective functions are usually transformed into a single-objective function using the weighted sum method, and then a mature single-objective optimization method, such as genetic algorithm (GA), is employed to optimize the single-objective function to obtain the optimal orbit [4–9]. Abdelkhalik and Mortari [4, 5] employed GA to optimize the weighted sum of squares of the distances between each target site and the satellite at the nearest ground track point, taking five orbital elements plus all visiting times as the design variables. Abdelkhalik and Gad [6] applied a weighted function of the total number of covered sites and the ground track repetition period as fitness function and adopted GA to optimize eccentricity, inclination, space-craft’s true anomaly above the first ground site, and the ground track repetition period, to design space orbits for Earth orbiting missions. Vtipil and Newman [7] and Vtipil [8] employed the sum of all time slot values of visiting as the cost function and adopted GA to conduct optimizations. The effect of population sizes was further researched. Zhang et al. [9] used a hybrid-encoding GA to optimize the sum of absolute value of velocity increment for long-duration rendezvous phasing missions. In weighted sum method [10–12], weight coefficients are utilized to transform the multiobjective function into the single-objective function. One disadvantage of the WSGA-based HM is that the artificially set values of the weight coefficients are unreasonable and subjective and depend significantly on human factors. In addition, the other disadvantage of the WSGA-based HM is the inefficient convergence of GA [13, 14].
To address these two disadvantages, this study proposes a population-based optimization method named CPC-based OM, in which candidate orbits are gradually randomly generated until the optimal orbit is obtained using a clustering via principal components based data mining method. A sufficient number of candidate orbits could ensure that the global optimal solution is obtained. In addition, the influence of the human factors from the weighted sum method is reduced in the optimization procedure because the candidate orbits are clustered based on the principal components rather than the weighted functions of the optimization objectives. Many methods have been investigated to reduce the influences of human factors of weighted sum method in multiobjective optimization [15–17]. Among them, the principal component analysis [18] is one of the most feasible methods, which transforms the original variables into a new set of variables, referred to as principal components, by using the eigenvalue-eigenvector method. The principal component analysis is thought to be the best way that explains the internal structure of the data [19, 20] and wildly applied by numerous researchers [21–26] to transform multiobjective functions for subsequent optimizations.
Methods must be introduced to accelerate convergence because the search procedure to obtain the optimal solution in CPC-based OM was a nearly exhaustive search with inefficient convergence. The methods of reducing feasible region are popular approaches [27–29]. In the methods of reducing feasible region, parts of the feasible region that do not include the optimum solution are deleted, and the subsequent optimization is accelerated because the remaining search space (feasible region) is smaller. Cluster analysis with the capability of dividing the feasible region into different regions has been used to reduce feasible regions [30–32]. Therefore, the second cluster analysis is introduced to CPC-based OM to accelerate convergence, and a novel population-based optimization method named DCPC-based OM is presented.
In this study, an orbit optimization model with constraints, six design variables, and eight optimization objectives is developed for Earth observation satellites. The process to obtain the optimal orbit using CPC-based OM is presented. To improve poor convergence of CPC-based OM, a more advanced DCPC-based OM is proposed by introducing cluster analysis based on six orbital elements. Finally, a test with numerical benchmark functions is conducted and the performances of DCPC-based OM, CPC-based OM, and WSGA-based HM on the orbit optimization of Earth observation satellites are compared.
2. Orbit Optimization Model for the Earth Observation Satellites
Abdelkhalik and Mortari [4, 5] explored the concept of developing an orbit based on target sites with no thrusters, in which design variables include five orbital elements and all the visiting times. The number of design variables increases with the increasing number of target sites. To avoid the increase of computational burden as the number of target sites increases, an optimization model with six orbital elements as design variables is employed. In addition, to deal with the increasing complexity of the observation mission, more orbit performance indices were taken into account than in prior studies [4–9].
2.1. Orbital Dynamics Model and Six Orbit Elements
For the orbit design of an Earth observation satellite without maneuvering, the relevant orbit dynamics equations in a geocentric equatorial inertial system (GEI) are as follows:where , , and are the components of the satellite position vector; , , and are the components of the velocity vector; and , , and are the components of external force, including Earth’s gravity (considering Earth nonspherical shape perturbation forces), atmospheric drag perturbation forces, and solar radiation pressure as well as lunar and solar perturbations forces [1]. The positions of the satellite at each moment can be calculated using (1) and six orbital elements at the initial moment. The six orbital elements include the semimajor axis , the eccentricity , the inclination , the argument of the perigee , the longitude of the ascending node , and the true anomaly . The key orbit performance indices of an Earth observation satellite are calculated using the positions of the satellite at each moment, the longitude and latitude data of the observation sites, the longitude and latitude data of the ground stations, and the right ascension of Greenwich at the initial moment [1].
2.2. Coverage and TT&C Performance Indices
Various key performance indices have been employed in the orbit design of the Earth observation satellites [4–9]. This paper adopts the key orbit performance indices which are systematically and comprehensively summarized by Wertz and Larson [1] and are applied by Wei et al. [2]. The key orbit performance indices of the Earth observation satellites can be separated into the coverage performance indices and the tracking telemetry and command (TT&C) performance indices. Among them, the coverage performance indices include the total coverage time (TCT), the frequency of coverage (FC), the average time per coverage (ATC), the maximum coverage gap (MCG), the minimum coverage gap (ICG), and the average coverage gap (ACG). And the TT&C performance indices include the average time interval of TT&C (ATI-TT&C) and the average time of each TT&C (AT-TT&C). The definitions of the eight orbit performance indices are as follows [1, 2].where is the total number of times of coverage in the simulation time , is the time of the th coverage, is the time of the th coverage gap, is the total number of coverage gaps, is the time of the th interval of TT&C, is the total number of the intervals of TT&C, is the time of the th TT&C, and is the total number of TT&C.
2.3. Orbit Optimization Model
The orbit optimization model of Earth observation satellites is shown in
The optimization objective is to make TCT, FC, ATC, and AT-TT&C be the maximum, the MCG, ICG, and ACG be the minimum, and the ATI-TT&C be within an expected range. The constraint is . The six orbital elements at the initial moment are the independent variables. The orbit optimization model in (3) is a typical multidimensional and multiobjective optimization problem. The purpose of this study is to provide an optimization method for orbit decision-making for an Earth observing satellite mission. The resulting optimal orbit may need to be refined in the case of additional system designs, including orbit stability, fuel consumption in orbital maneuvering, or launch site restrictions.
3. CPC-Based OM for Orbit Optimizations
To reduce the influences of human factors in orbit optimizations [4–9], CPC-based OM is presented, and the accelerating convergence approach will be introduced in Section 4 to develop DCPC-based OM. The process flow of orbit design optimization with CPC-based OM is shown in Figure 1, including five steps.
Step 1. Randomly generate candidate orbits according to the feasible regions of the six orbital elements, and then calculate the coverage and TT&C performance indices.
Step 2. Nondimensionalize the coverage and TT&C performance indices of candidate orbits.
Step 3. Calculate the principal components of nondimensionalized orbit performance indices of candidate orbits, divide candidate orbits into classes by performing cluster analysis based on the principal components, and evaluate all class centers using the weighted sum function of key orbit performance indices to obtain the optimal class.
Step 4. If the number of orbits in the optimal class is more than six, go to Step and take orbits in the optimal class as candidate orbits in Step .
Step 5. If the number of orbits in the optimal class is not more than six, determine the temporary optimal orbit from the optimal class by using the weighted sum method. Randomly generate additional candidate orbits and go to Step with candidate orbits, until the relative improvement of the temporary optimal orbit is less than a certain bound. The applied stopping criterion refers to that used in the GA [4–9].
Note that Steps 3 and 4 constitute “multilevel cluster analysis,” in which the optimal class with not more than six orbits is obtained. Steps 2~5 constitute “random exhaustive” process to obtain the optimal orbit. In addition, because the cluster analysis is based on the principal components of nondimensionalized key orbit performance indices rather than based on weighted function of the nondimensionalized key orbit performance indices in Step , the influence of human factors on the optimum result is reduced. The principal components, which explain the internal structure of the key orbit performance indices [19, 20], are suitable for cluster analysis. However, they do not have physical meaning and thus can not be used for evaluating the quality of orbits. Therefore, the optimal class is selected from all classes utilizing weighted sum method, rather than utilizing principal component analysis. Similarly, in Step , the temporary optimal orbit is obtained from the optimal class by using the weighted sum method, rather than by principal components analysis. Therefore, CPC-based OM, to some degree, can mitigate the influence of human factors and is advanced. Some details of CPC-based OM are briefly introduced below.
3.1. Index Nondimensionalization Method
When various candidate orbits are generated and orbit performance indices with different units are calculated, the performance indices are nondimensionalized using (4) [2]. If performance indices are to be maximized, the dimensionless coefficient is equal to the performance index value divided by the maximum value in all candidate orbits. If performance indices are to be minimized, the dimensionless coefficient is equal to the reciprocal of the performance index value divided by the minimum value in all candidate orbits. If performance indices are to be within an expected range, the dimensionless coefficient is equal to the performance index value divided by the median of the range when the performance index is less than the median of the range and dimensionless coefficient is equal to the reciprocal of the performance index value divided by the median of the range when the performance index is greater than the median of the range.where is the th performance index of the jth orbit, is the dimensionless coefficient of the performance index , is the maximum value of in all candidate orbits, is the minimum value of in all candidate orbits, and is the median of the range.
3.2. Principal Component Analysis
Supposing the dimensionless coefficients can be represented by vector , where is the number of performance indices of each orbit, then are the dimensionless coefficients of the performance indices of candidate orbits. In principal component analysis [2, 21], the elements of the covariance matrix are firstly calculated as follows:where
Next, the eigenvalues of the covariance matrix are calculated. The cumulative contribution ratio of the previous eigenvalues is as follows:
The number of principal components is the minimum , which makes greater than 88%. Then, by using the orthogonal normalized eigenvector of the covariance matrix , the principal component vector of the th candidate orbit is calculated as follows:
Because is generally smaller than , the total number of principal components to be clustered is less than the total number of performance indices.
3.3. Multilevel Cluster Analysis
The candidate orbits could be clustered according to the principal components to obtain the optimal class [33, 34]. The max–min distance method [35] is employed in this research. To avoid randomly selecting the first cluster center, the orbit most close to the zero point of principal components is selected as the first cluster center. The Euclidean distance is employed in cluster analysis, and the definition is as follows:where and are two arbitrary principal component vectors of candidate orbits, and are two arbitrary natural numbers, and is the th principal component in . A total class distance criterion is applied to determine the optimal number of classes in each clustering and the candidate values of include 4, 5, and 6.
The procedures of cluster analysis are as follows.
Step 1. Assume is set as one of the three candidate values.
Step 2. Select initial cluster centers. Calculate Euclidean distance from all to the zero point of principal components using corresponding to minimum is taken as the first cluster center . The second cluster center will be which has the maximum Euclidean distance to . The th cluster center will be corresponding to the maximum , where , and is the Euclidean distance between and . cluster centers are finally selected in this way.
Step 3. Calculate Euclidean distances from each of the remaining to all the cluster centers. Each will be assigned to the most close cluster center.
Step 4. Determine new cluster centers. Calculate the average value of principal components of classes using where is the average value of the th principal component in the kth class and is the total number of candidate orbits in the kth class. The new cluster center of each class will be most close to the average value of principal components. The distance from to average value of principal components is shown in where is the distance from to the average value of principal components of the th class.
Step 5. If one or more than one new cluster center is different from the corresponding last cluster center, reassign all remaining to the cluster centers, then repeat Step ; otherwise stop the clustering analysis corresponding to this candidate value of .
Step 6. Set to other candidate values, and relevant cluster analyses are carried out according to aforementioned Steps 2~5. Then, the total class distances in three cluster analyses corresponding to three candidates are calculated [2].where is the candidate value of in each cluster analysis (, 5 or 6), is the number of orbits in class , and is the Euclidean distance between the class center and the other orbits in this class. The final cluster result is the cluster result corresponding to the smallest .
The number of candidate orbits is very large, after one cluster analysis, the number of orbits in the optimal class is usually more than six. Therefore, cluster analyses are repeatedly carried out until the number of orbits in the optimal class is less than six. This is known as multilevel cluster analysis.
3.4. Weighted Sum Method
The weighted sum method [2, 10] is adopted to determine the optimal class from all the classes and the optimal orbit from the orbits in the last optimal class. The evaluation index is defined as follows:where () is the weight coefficient of each performance index, in the range . The representative coefficient is given by where is the expected range of and is the dimensionless coefficient of .
4. Novel DCPC-Based OM for Orbit Optimizations
4.1. The Advanced Characteristic of CPC-Based OM and Modification for Convergence Efficiency
In CPC-based OM, principal component analysis, rather than weighted sum method, is adopted to cluster orbits, and therefore the negative influence of the artificially set weight coefficients in weighted sum method is reduced. The detailed analyses are illustrated in contours resulting from principal component analysis and weighted sum method in Figure 2(a), assuming that two orbit elements are optimized and orbit performance indices are transformed into one principle component. Because artificial weight coefficients are usually different from the orthogonal normalized eigenvector (refer to Section 3.2), significantly different contours are shown in Figure 2(a). The artificial weight coefficients are unreasonable and easily influenced by human factors, while principal components are objective and therefore more reasonable. A comparison analysis of the advantage of principal component analysis is presented in Section 6.3.
(a)
(b)
(c)
(d)
(e)
In CPC-based OM, when enough candidate orbits are randomly generated, the global optimal orbit with high accuracy could be achieved. For a certain number of candidate orbits, as indicated by black stars in Figure 2(a), star A is the closest to the global optimal orbit indicated by a large red star. Therefore, star A is the optimal solution among these candidate orbits. Orbit A is not a local optimal solution but a global optimal solution with low accuracy. If the number of candidate orbits is smaller than that in Figure 2(a) and orbit A is not included in the candidate orbits, the optimal orbit obtained using CPC-based OM will be star B (as shown in Figure 2(b)); orbit B is a local optimal solution. Therefore, if there are a large number of candidate orbits, the global optimal solution can be achieved by CPC-based OM. To improve the accuracy of the optimal solution, generating more candidate orbits, is a feasible approach, as shown in Figure 2(c). When additional orbits (indicated with blue stars) are generated as candidate orbits, a better orbit marked with star C is the optimal solution. Star C is closer than star A to the global optimal orbit. Therefore, the accuracy of the global optimal solution is satisfied when sufficient candidate orbits are prepared in CPC-based OM. Additional candidate orbits are gradually generated in “random exhaustive” process in CPC-based OM as shown in Figure 1, which ensure global optimum and high accuracy in CPC-based OM. But “random exhaustive” leads to low computational efficiency for CPC-based OM.
The method of reducing the feasible region will be employed to improve the convergence. It could be known from Figure 2(c) that the additional orbits cover the entire feasible region of orbit elements. Supposing that all orbits in the green contour belong to the optimal class, the rectangular box R2 (shown in Figure 2(d)), which is envelop of green contour, is the feasible region of the optimal class, and it is far smaller than the entire feasible region R1 (shown in Figure 2(d)). If the same numbers of additional orbits as that in Figure 2(c) are generated and they cover the feasible region R2 of the optimal class, as shown in Figure 2(d), star D will be the optimal solution (better than star C). Unfortunately, some orbits (marked with a red dotted box in Figure 2(d)) in R2 do not belong to the optimal class. To avoid orbits generated in red dotted box, cluster analysis of orbits in optimal class based on the orbital elements could be introduced to divide optimal class to two parts, and the concentrated zones R3 and R4 will be the feasible regions of two parts of the optimal class, as shown in Figure 2(e). The areas of the new feasible regions R3 and R4 are smaller than the area of R2. If the same numbers of additional orbits in feasible regions R3 and R4 are generated, star E is the new optimal solution (shown in Figure 2(e)) and is better than star D (shown in Figure 2(d)).
4.2. Processing Flow
By introducing cluster analysis based on the six orbital elements, a novel population-based optimization method named DCPC-based OM is developed based on CPC-based OM. The processing flow is shown in Figure 3, where the gray zone indicates the cluster analysis based on the six orbital elements, including totally five steps.
Step 1. Randomly generate candidate orbits according to the initial feasible regions of the six orbital elements, and then calculate the coverage and TT&C performance indices.
Step 2. Nondimensionalize the coverage and TT&C performance indices of candidate orbits.
Step 3. Calculate the principal components of nondimensionalized orbit performance indices, divide candidate orbits into classes by performing cluster analysis based on the principal components, and evaluate class centers using the weighted sum method to obtain the optimal class.
Step 4. When the total number of candidate orbits in the optimal class is greater than six, the six orbital elements are nondimensionalized by the upper limit value of the relevant feasible region. Cluster analysis is conducted on the orbits belonging to the optimal class, by using the nondimensionalized values of the six orbital elements. The candidate number of classes here is not limited to 4, 5, or 6; that is, it could change from 1 to the total number of orbits in the optimal class. After orbits with similar six orbital elements are clustered into a class, the concentrated feasible region of a class is the envelope of six orbital elements of all orbits in this class. Then additional candidate orbits are generated randomly in the all concentrated feasible regions. The value of is as follows:where is the total number of candidate orbits in the ()th cluster analysis and is the total number of orbits in the optimal class of the th cluster analysis. An evolutionary rate is defined to ensure that the number of candidate orbits is reduced in the subsequent optimization, and the total number of candidate orbits in the th cluster analysis is .
Step 5. Choosing the orbits in the optimal class and additional orbits as candidate orbits (a total of orbits), repeat Step , until the number of orbits in the optimal class is no more than six and then determine the optimal orbit from the optimal class by using the weighted sum method.
The characteristics of DCPC-based OM can be concluded as follows.
() Cluster analysis based on principal components of orbit performance indices and six orbital elements are performed, respectively. Therefore, double cluster analyses are included in DCPC-based OM. By clustering of the orbits belonging to the optimal class using six orbital elements, the concentrated feasible regions of the optimal class are obtained. Additional candidate orbits are generated only in the concentrated feasible regions which is smaller than initial feasible region. An evolutionary rate () is defined to ensure that the candidate orbits is reduced in the subsequent optimization.
() According to the definition of evolutionary rate in Step and the flow chart in Figure 3, the total number of iterations N in DCPC-based OM is approximately as follows:where 5 is the average value of the number of classes (4, 5, and 6) and 3.5 is the average number of the orbits in the last optimal class (1, 2, 3, 4, 5, and 6). Based on the hypothesis that each class has the same number of orbits in cluster analysis, the number of additional orbits in N iterations can be calculated using
It can be seen from (18) that the numbers of additional orbits in all iterations are in a form of geometrical sequences. Therefore the total number of generated candidate orbits equals the number of initial orbits plus the sum of geometrical sequences in (18), as shown in
For orbit optimization, the calculation procedure of orbit performance indices is more time-consuming than optimization operation, because a large number of numerical computations are needed to calculate orbit performance indices. Therefore, the total optimization time approximately equals the total number of generated candidate orbits multiplied by the computation time required for the performance indices of one orbit, which means the time cost is nearly predictable. The proposed orbit design optimization method with predictable time cost is more convenient than other population-based optimization methods with unpredictable time cost for a scheduled project.
() In DCPC-based OM, the method of reducing the feasible region improved the computational efficiency while it might result in the reduction of capability of global optimization, because the feasible regions including the global optimum solution might be mistakenly deleted. The large value of and can improve the capability of global optimization, and the small value of and can improve computational efficiency but increase the risk of missing the global optimal solution.
5. Experiment and Analysis
A four-objective benchmark function (as shown in (20)) including two simple multimodal problems ( and ) and two unrotated multimodal problems ( and ) is adopted to testify the proposed method [36]. The Rosenbrock function is modified to to ensure that global optimal solution is [] and the minimum value is 0. The global optimal solution of the multiobjective benchmark function is and .
The optimizations with dimensions of 5 and 10 are conducted using DCPC-based OM, CPC-based OM, and WSGA-based HM [4–6]. In DCPC-based OM and CPC-based OM, the number of initial candidate elements was set as 10,000,000. With respect to CPC-based OM, the increasing number was 100,000. In DCPC-based OM, the evolutionary rate was 0.6. In WSGA-based HM, the individual number in a generation was 100,000, and the fitness values were calculated using weighted sum method, as shown in Section 3.4, with weight coefficients (). The initial range of was []. At each time the clustering in DCPC-based OM was finished, the minimum values of function calculated by the three methods are shown in Figure 4. Three methods started optimization process at the same time and terminated simultaneously when DCPC-based OM completed optimization. Because in WSGA-based HM was far less than in DCPC-based OM and CPC-based OM, many iterations in WSGA-based HM had occurred when the first clustering in DCPC-based OM was completed. In addition, with repeat clustering, the number of candidate elements in the optimal class decreases in DCPC-based OM and the number of completed iterations of WSGA-based HM in one clustering decreases. Therefore, in each clustering of DCPC-based OM, the number of iterations in WSGA-based HM varies. The number of times of completed iterations in WSGA-based HM at each time when clustering in DCPC-based OM was finished is shown in Figure 4.
(a)
(b)
Figure 4(a) shows that the function values in all three methods decreased as the number of times of clustering increased. A comparative analysis between the function values of DCPC-based OM and CPC-based OM shows that, in the first clustering, the function values of DCPC-based OM were similar to that of CPC-based OM. However, in the last clustering, the function values of DCPC-based OM were far lower than that of CPC-based OM. It validated the improved convergence by the method of reducing feasible region in DCPC-based OM. When the first clustering in DCPC-based OM was completed, iterations in WSGA-based HM had been finished 64 times. Thus, the function values in WSGA-based HM were less than that in DCPC-based OM. By the 8th clustering in DCPC-based OM, 104 iterations of WSGA-based HM had been finished. And the additional 2 iterations of WSGA-BASED HM were completed by the 22nd cluster. At the end of the last clustering, the function values of DCPC-based OM were less than that in WSGA-based HM. In WSGA-based HM, the feasible region is constant, the convergence is driven by evolutionary capacity of genetic operation, but genetic operation is actually random. In a large feasible region, the global optimum is difficult to be completely randomly generated. Hence, the convergence efficiency is very low in WSGA-based HM. The candidate elements in DCPC-based OM are randomly generated. However, because reducing the feasible region is an effective method to improve optimization efficiency [27–32], DCPC-based OM achieves more efficient convergence than WSGA-based HM.
6. Orbit Design Results and Analysis
6.1. Orbit Optimization Conditions and Convergence Analysis
For a certain Earth observation satellite, five observation targets are with latitude and longitude coordinates of (25°N, 120°E), (10°N, 110°E), (40°N, 130°E), (15°N, 90°W), and (20°S, 130°E) and the same vision field angle of 25°. The minimum elevation angle, latitude, and longitude coordinates of TT&C station are 5°, 40°N, and 120°E, respectively. The expected range of ATI-TT&C is 8000 s~40000 s. The range of the semimajor axis is 400 km~600 km and (). The simulation time for each candidate orbit is 3 days.
The orbit design optimization is conducted using DCPC-based OM, CPC-based OM, and WSGA-based HM [4–6]. In DCPC-based OM and CPC-based OM, the number of initial candidate orbits is set as 100,000. In DCPC-based OM, the evolutionary rate is 0.5. The number of additional orbits is set as 10,000 for CPC-based OM. The number of individuals in a generation is 10,000 for WSGA-based HM. The change of the evaluation indices for DCPC-based OM, CPC-based OM, and WSGA-based HM with the number of times of clustering in DCPC-based OM is shown in Figure 5.
The data in Figure 5 show that the evaluation indices for all three methods increase as the number of times of clustering increases. After the first clustering, DCPC-based OM and CPC-based OM have similar evaluation indices, which are less than that for WSGA-based HM. The optimization process for CPC-based OM and WSGA-based HM is stopped when DCPC-based OM finishes its optimization, and DCPC-based OM has the highest evaluation index of all three methods in the last clustering. The evaluation index of the orbit of DCPC-based OM is higher than that of CPC-based OM by 22.9% and higher than that of WSGA-based HM by 8.0%. In other words, DCPC-based OM has more efficient convergence than both CPC-based OM and WSGA-based HM.
6.2. Optimization Results with DCPC-Based OM
Details of orbit optimization results with DCPC-based OM are presented in this section. In the first clustering, principal component analysis was conducted after the performance indices of all candidate orbits are calculated and nondimensionalized. The contribution ratios of the first three principal components in the principal component analysis were 74.10%, 11.69%, and 8.38%, respectively, as shown in Figure 6.
The data in Figure 6 shows that the first three principal components contribute more than the others. The cumulative contribution ratio of the first three principal components is 94.17%. The number of principal components is three in all the clusters in this study, although the contribution ratios of the first three principal components vary slightly in various principal component analyses. Consistently, the multilevel cluster analysis is conducted.
According to total class distance criterion, the optimal number of classes in the first cluster was five. The principal components of the cluster centers in the five classes are shown in Table 1. The data in Table 1 illustrates the differences among the cluster centers of the five classes. Classes 4, 2, and 3 had the largest principal components 1, 2, and 3, respectively.
The data in Figure 7(a) show the distribution characteristics of the principal components in the five cluster centers. Because the principal components have no physical meaning, it is difficult to determine the optimal class according to the principal components. Therefore, the data in Figure 7(b) show the differences among the performance indices of the five cluster centers, and the weighted sum method is used to select the optimal class by evaluating the performance indices of the cluster centers. The eight dimensionless performance indices of the five cluster centers are listed in Table 2.
(a)
(b)
The data in Figure 7(b) indicate that the dimensionless orbit performance indices of Class 4 are relatively large, and the optimal class obtained by using the weighted sum method was also Class 4.
The orbits in Class 4 from the first clustering and additional candidate orbits are then selected to continue with the clustering process. The principal components and dimensionless performance indices in the cluster center after clustering fifteen times using DCPC-based OM are shown in Figures 8(a) and 8(b). The analysis results using the weighted sum method indicate that Class 3 is the optimal class.
(a)
(b)
The eight dimensionless performance indices for all orbits of Class 3 are shown in Figure 9. The optimal orbit obtained by the weighted sum method is Orbit 1.
6.3. Global Optimality and Principal Component Estimation in DCPC-Based OM
The evolution process of the maximum and minimum evaluation indices of the orbits in the optimal class is shown in Figure 10. The uptrend of the two types of indices indicates that the orbits are being optimized. The phenomenon of the two sets of indices approaching each other indicates that the differences among all orbits in the optimal class are decreasing. When the two sets of indices are close enough, the optimal orbit is achieved.
A decrease in the maximum evaluation index from O1 to O2 indicates that, in the 8th clustering, the maximum evaluation index is not contained in the optimal class and is filtered out. The coverage and TT&C performance indices of orbits O1 and O2 are listed in Table 3. Compared to orbit O2, orbit O1 exhibits better TCT and ATC, worse MCG, ICG, and AT-TT&C, and similar FC, ACG, and ATI-TT&C. Considering an obviously oversized MCG and ICG in orbit O1, orbit O2 is better than orbit O1, even though the evaluation index of orbit O2 determined by weighted sum method is less than that of orbit O1. The reason for the higher evaluation index of orbit O1 is that ATC is linearly correlated with TCT, and in the weighted sum method, both TCT and ATC are considered, that is, the coverage time is counted twice. It leads to the fact that effects of the oversized MCG and ICG are masked in weighted sum method. However, the effect of the linear correlation is eliminated in the principal component analysis, which leads to the fact that the coverage time is counted only once and the negative effects of the oversized MCG and ICG are reflected. Thus, the principal component analysis provides a more accurate estimation.
7. Conclusions
This paper proposes a population-based optimization method named DCPC-based OM, which consists of the index nondimensionalization method, principal component analysis, double cluster analysis, and the weighted sum method. Tests using numerical benchmark functions were conducted, and an example of orbit optimization for Earth observation satellites was analyzed. Both optimization results show that the proposed method, with characteristics of a predictable time cost, has the advantages of reducing the influence of human factors that commonly exist in the weighted sum method and bring more efficient convergence than genetic algorithm.
This paper describes the results of a preliminary research study of the developed optimization method. Further study will be performed, such as determining the quantity of initial orbits to ensure the capability of global optimization, understanding how evolutionary rate influences the accuracy of the optimal solution, and determining whether a dynamic evolutionary rate is necessary. In addition, the optimal capability when the optimal solution is on the boundary of the feasible region will be further investigated.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The research work is supported by the National Defense 973 Program (Grant no. 613237).