Geostationary orbit (GEO) plays an important role in communication, navigation, and other areas   . With an orbital period equal to the Earth’s rotation period, the GEO satellites are always above the Earth’s equator, and are ﬁxed relative to an Earth-ﬁxed observer  . This principal characteristic of the GEO makes it suitable for communication and navigation satellites  . Meanwhile, it is essential of the frequency resources for construction of the communication and navigation satellites’ systems  . Hence, it is of signiﬁcance to monitor the usage of the frequency resources of GEO satellites. Based on this background, problem A of CTOC9 has considered the orbit design and optimization mission for the low earth orbit (LEO) satellites in formation to observe the GEO satellites’ beams  . The J2 perturbation must be considered when propagating the orbits and maintaining the flying formation  .
The problem is a multi-object optimization problem, which appears in many practical applications. For example, the multi-asteroid exploration, the multi-planet exploration and the multi-debris clearing mission   . At the end of the CTOC9, teams from NUAA  , NUDT  , and Tsinghua University  have obtained the first, second, and third prizes. And we have obtained the fourth. Most of them used the regressive orbit to solve the special problem, but our methods are different from theirs. We attempted to propose more general methods to solve not only this problem but also other multi-objective trajectory design problems.
This paper is organized as follows. Section 2 gives a brief description of the problem. Section 3 gives the general description of the methods. Application of the methods to the problem is presented in Section 4. Results and conclusions are finally presented in Sections 5 and 6.
2. Brief Description of the Problem
In the problem, the GEO constellation consists of 18 satellites, each of which is equipped with a transmitter which directs its beam to the ground station. The diagram of the projection of the GEO satellites’ beams is shown in Figure 1.
An LEO formation consisting of 3 monitoring satellites is inserted into circular orbits of altitude of 400 km and inclination of 42.8˚ at initial time. The task of the three monitoring satellites is to monitor all the beams of the GEO satellites. The “monitoring” operation of a GEO satellite beam is executed only when the three monitoring satellites are in the coverage of the beam. Monitoring proﬁt of a GEO satellite, which is deﬁned by the geometry between the three monitoring satellites and the GEO satellite, is accumulated in the process of monitoring execution. Assuming the three monitoring satellites are in the coverage of a GEO satellite’s beam at a specified time, Figure 2 shows the geometry of the three monitoring satellites and the GEO satellite. When the monitoring proﬁts of all 18 GEO satellites reach a certain threshold (In the problem, it is 1 × 107 km2), the whole monitoring task is completed. Please refer to Ref.  for more information of the monitoring profits’ computational method and the overall description of the problem.
There are two optimization objectives of diﬀerent priorities in the problem. The primary objective is to minimize the time duration of the whole monitoring
Figure 1. The diagram of the projection of the GEO satellites’ beams  . (a) Refers to the two-dimensional view and (b) Refers to the three-dimensional view.
Figure 2. Geometry of the three monitoring satellites and one GEO satellite  .
task. The secondary objective, in case two or more teams achieve the same primary objective, is to minimize the total fuel consumption of the three monitoring satellites.
Each monitoring satellite is subject to the J2 perturbation. And the motion equation in the J2000 earth centered inertial coordinate is as follows,
In which, is the gravitational constant of the earth. x, y, and z refer to the position of the satellites. vx, vy, vz refer to the velocity. refers to the range from the earth center to the satellite. Re is the radius of the earth, and J2 is the perturbation constant of the earth’s oblateness.
The initial masses of all satellites are 500 kg, and the allowed fuel consumption should not be more than 200 kg. The maneuvers of the monitoring satellites are instantaneous changes of the satellites velocities. After each maneuver, the fuel consumption Δm should be evaluated using the Tsiolkovsky equation.
where m is the mass before the maneuver, ΔV is the magnitude of the maneuver, ge is the acceleration of gravity at sea level of Earth, Isp is the specific impulse of the propulsion.
The GEO satellites are assumed to follow non-perturbed Keplerian orbits. The longitude and radius of each GEO satellite remain constant, and the latitude remains zero during the whole mission. For more information of the problem description, please refer to Ref.  .
3. General Description of the Method
3.1. Outline of the Method
In order to solve the problem effectively, we decomposed the methods used by us into following steps.
Step 1: The formation is designed according to the observation demands by analyzing the observation features.
Step 2: The flying sequence is firstly determined by a specified reference satellite using a proposed improved ephemeris matching method (IEMM). Then, the formation is maintained and transferred following the reference satellite by a multi-impulse control method (MICM). The total monitoring profit is computed by propagating the satellites’ orbits subject to J2 perturbation according to the sequence and transfer strategies.
Step 3: The result data is checked and submitted to the competition organizer.
The basic procedure of the methods employed in this paper for solving the problem A of CTOC9 is described in Figure 3.
3.2. Formation Designing
In order to design stable formations and reduce the maintaining cost, we analyzed the stability of the formations by propagating the differential equations of the elements difference between the satellites. In the presence of the J2 perturbation, some orbital elements change, including the drifts in perigee and mean anomaly, and the nodal regression. Employing orbit-averaged quantities   , differences in the orbital elements between the leader and the follower satellites are given as below:
where μ is the gravity constant of Earth. a, e, i, Ω, ω, Μ are the orbital elements of the satellites, and δa, δe, δi, δΩ, δω, δΜ are the elements’ differences between
Figure 3. The procedure of the methods employed for solving the problem A of CTOC9.
the leader and the follower satellites. The elements with subscript L refer to the leader satellites.
Equation (3) shows that the J2 perturbation does not change the differences of a, e, i but changes the differences of Ω, ω, Μ between the leader and the follower satellites. If the initial values of the elements’ differences with subscript “0” at time t0 are given, the values at any time t can be computed as follows,
Equations (3) and (4) obviously show that, when there are differences in the elements a, e or i, the formation is unstable, and changes quickly. As a result, in this paper, the three satellites’ initial elements a, e and i are chosen to be nearly the same.
The initial relations of the three satellites’ elements will be determined by the mission’s requirements, which will be shown in Section 4.1.
3.3. Flying Sequence Searching
For the flying sequence searching, the improved ephemeris matching method (IEMM) is proposed.
From the characteristics of the Lambert problem  , it is not hard to understand that, the closer distance from the satellite’s position (propagated to the target time without any control from the leaving beam bl) to the position of the target beam b (propagated to the target time) is, the less energy is needed for the satellite to transfer from bl to flyby b. Thereby we can measure the relative size of the transfer energy consumed by examining the closeness between the target beam’s ephemeris and the satellite’ ephemeris propagated from the leaving beam. In this paper, we define the reference satellite as the leader satellite. When propagating the reference satellite’s orbit, the satellite’s velocity when leaving the beam is not fixed. We add a small Δv along the velocity in range [−25 m/s, 25 m/s] at the departure epoch with a step size of 5m/s and propagate the satellite’s trajectory. We pick up the one closest to the matching target beam (if the monitoring profit for this beam is not fulfill the requirement), record the Δv value, the arrival time tf, and the position R and velocity V of the closest point from the central line of the beam to the reference satellite’s propagated orbit. This process is different from the ephemeris matching method (EMM) in Ref.  , in which the trajectory is propagated without Δv along the velocity, so we call it improved ephemeris matching method (IEMM).
For a segment of the transfer trajectories, the mathematical model of IEMM can be expressed as
where f (.) is a monotonic increasing function; R and V are respectively the position and velocity of the satellite’s propagated orbit at arrival time t; Rb and Vb are respectively the position and velocity of a point on the central line of a beam b. The point is chosen to be with the same height to the satellite’s propagated orbit at the arrival time. λr and λv are the weights of the positions’ difference and velocities’ difference. For this paper, the satellites need to flyby (not rendezvous) the target beam, so the velocities between them need not to be equal. Only the positions need to match. Therefore, λv is set to be 0, and λr could be 1.
For a transfer segment, when we have set a Δvi (which is determined by step 5m/s in the small range [−25, 25] m/s from lowest to highest), the target beam bi whose ephemeris matches the reference satellite could be found along the propagating satellite’s orbit till ti. Then we identify the target beam and record the corresponding parameters that meet the condition of top 5 values of the net accumulated monitoring profits when propagating the satellites orbits in the segment. Here net profit means when the accumulated monitoring profit of a GEO satellite’s beam has up to the threshold, it will not be added to the total accumulated value computed along this segment. When computing the monitoring profit, the formation is formed theoretically by the geometry relationship between the leader and the follower satellites, because in this sequence searching process, the formation maintaining is not considered.
The procedure of the IEMM algorithm is shown in Figure 4, and is described as follows:
Figure 4. The procedure of the IEMM algorithm. The solid arrows and circles refer to the recorded sequences, and the dotted arrows and circles refer to the abandoned sequences. The red arrows are connected to show an example sequence c.
Step 1: Propagate the reference satellite’s orbits from the last beams of the recorded sequences with the initial velocity offset Δvi from −25 m/s to 25 m/s with a step size of 5 m/s. There will be 11 orbits to be propagated.
Step 2: Find the 11 target beams whose central line is closest to the satellite’s orbits corresponding to Δvi at the end time ti.
Step 3: Propagate the reference satellite’s orbit and use the theoretical formation (by the geometry relationship between the leader and follower satellites) to compute the net accumulated monitoring profits.
Step 4: For a sequence, if the total monitoring profit meets the condition of all the beams’ monitoring profits up to the threshold, record the new sequence as one of the candidate sequences. If any one of the satellite’s fuel has all been consumed but not all the observation values of the beams have been up to the threshold, the sequence is abandoned. Otherwise, abandon the 6 sequences with the least monitoring profits, record the top 5 sequences for further searching, and turn to Step 1.
Step 5: If all recorded sequences have fulfilled the threshold, stop searching, and pick out the sequence with shortest mission time for further processing.
3.4. Formation Maintaining
For the formation maintaining and orbit transfer, a multi-impulse control method (MICM) is employed. The Δvs of all maneuvers are optimized by the differential evolution (DE) algorithm and the J2 perturbed Lambert problem solving algorithm. The algorithm procedure is shown in Figure 5.
Figure 5. The procedure of the MICM algorithm.
There are two layers of the loop. The inner layer is the solving of the J2 perturbed Lambert problem. The locally optimization shooting method is used to solve the boundary value problem. If the three satellites’ final positions Rf and velocities Vf equal the values predetermined by the theoretical relationship of the formation, the problem is solved. The predetermined values are given by the sequence searching and formation designing algorithms. When the position and velocity of the reference satellite are provided, based on the relationship between the three satellites in formation, the required positions and velocities of the follower satellites are determined, which are the objective values of the shooting method.
The outer layer of the loop is using the DE algorithm to optimize the velocity increments (Δvs) to be minimal. The optimizing problem is a time-fixed optimal transfer problem. DE is a heuristic searching algorithm introduced by Storn and Price  . Its remarkable performance as a global optimization algorithm on continuous numerical minimization problems has been extensively explored  . As with other evolutionary algorithms, DE solves optimization problems by evolving a population of candidate solutions using alteration and selection operators. For the problem in this paper, the objective function is as follows:
where Δvi,k refers to the sum of the velocity increments of the ith transfer of satellite k. The variables to be determined by the optimization are as follows,
in which Φ1,k indicates the variables relative to the three-impulse transfer for the first maneuver of satellite k; Δt1,k is the waiting time when the satellite is flying along the initial orbit; Ω1,k, ω1,k, and M1,k are respectively the right ascension node, the argument of perigee, and the mean anomaly of the initial orbit; Δt’1,k is the transfer time from the leaving point to a middle point on an initial J2 perturbed Lambert orbit; ΔR1,x, ΔR1,y, and ΔR1,z are the position offsets from the middle points. Then two Lambert orbits are computed again to fulfill the transfers from the leaving point to the middle point and from the middle point to the end point.
Φc,k indicates the variables relative to the three-impulse transfer for formation changing (Why the formation need to be changed? Please refer to section 4.1.) of satellite k; Δtc,k is the waiting time, at the end of which the satellite start to transfer; Δt’c,k is the transfer time from the leaving point to a middle point on an initial J2 perturbed Lambert orbit; ΔRc,x, ΔRc,y, and ΔRc,z are the position offsets from the middle points. Then two Lambert orbits are computed again to fulfill the transfers from the leaving point to the middle point and from the middle point to the end point.
Φi,k indicates the variables relative to the two-impulse transfers for the rest maneuvers of satellite k; Δωi,k is the difference of the argument of perigee between the leader satellite and the follower satellite. When propagating the orbits, only small difference of Δωi,k will be generated by the J2 perturbation. In order to save fuel, we do not require the argument of perigee for the 3 satellites exactly equal each other. But when solving the J2 perturbed Lambert problem, the offset of the target mean anomaly should be set to have an opposite value to Δωi,k. This will ensure a relatively stable shape to fulfill the observation. The bounds of the variables are shown in Table 1.
3.5. Monitoring Profit Computation
Incorporating all the maneuvers of the 3 monitoring satellites, the observation values Wj are computed by propagating the J2 perturbed orbits with a step time tm of 1 second using Equation (9)  :
where j refers to the label of the GEO satellite’s beam, Sj is the area of the triangle formed by 3 monitoring satellites, and gk is sine of the angle from the GEO satellite’s beam to the triangle plane. The results are organized to fit the submission format in order to pass the data checking process. M refers to the total duration (in second) when the triangle is in the beam j.
4. Application of the Method to the Problem
4.1. Formation Determination
In order to monitor the 11th GEO satellite whose beam is locating at the highest latitude (which is shown in Figure 6), we raise the apogees of the three orbits and make the formation be a triangle parallel with the orbit plane to reach the beam. Of course if we raise the inclinations of the orbits, the 11th beam could also be monitored, but it will consume too much fuel. So we choose the first way.
Table 1. Bounds of the variables.
Note: ΔT1 is the time range from initial to the first beam in the sequence. ΔTc is the time range from the leaving to the arriving point of the formation changing transfer.
Figure 6. The three-diminsional geometry of the 11th beam, which is the highest one in the figure. The yellow thin lines are the satellites’s initial orbits.
When the monitoring profit of the 11th GEO satellite’s beam is near 1 × 107 km2, the formation is changed from the first one to the second one which constitutes a triangle perpendicular to the orbit plane for the observing of the remaining GEO satellites’ beams locating at the low latitudes. The first and the second formations are shown in Figure 7.
When designing the formations, not only the monitoring requirements, but also the stability should be considered to reduce the formation maintaining cost. Finally, the relations of the 3 satellites are determined and shown in Table 2.
The leader satellites’ elements aL, eL, iL, ΩL, ωL, ΜL are determined by the flying sequence and transferring strategies. The follower satellites are transferred to fit the relations in Table 2. The values of differences between the elements of the leader and the followers are determined by optimization to find the highest mean monitoring profit for passes of the GEO satellite’s beams.
4.2. Flying Sequence Determination
For the flying sequence determination, the proposed IEMM described in Section 3.3 is employed. For every searched sequence, the labels of the beams, the epochs when the leader satellite arrive the beams and the positions & velocities of the orbits are recorded. Then the top 5 beams with high monitoring profits are chosen and added to the sequence to form new sequences. For a special sequence, if the total monitoring profit meets the condition of all the beams’ monitoring profits up to the threshold, record the new sequence as one of the candidate sequences. If any one of the satellite’s fuel has all been consumed but not all the monitoring profits of the beams have been up to the thresholds, the sequence is abandoned. Otherwise, abandon the last 6 values, and record the top 5 sequences for further searching. If all recorded sequences have fulfilled the threshold condition, stop searching, and pick out the sequence with shortest mission time for further processing.
The sequence searching result is shown in Figure 8, which shows that most of the durations between two target beams of the ephemeris matching processes are no longer than 1 × 105 s (27.8 hours). Most of the distances between the
Figure 7. The designed formations. The red star refers to the reference satellite which is also the leader satellite, and the other two are the follower satellites. (a) The first formation constitutes a triangle parallel to the orbit plane. (b) The second formation constitutes a triangle perpendicular to the orbit plane.
Figure 8. The sequence searching result. The first sub-figure refers to the beam numbers of the best sequence. The second sub-figure refers to the duration between two points of the sequence. The third sub-figure refers to the distance between the reference satellite’s orbit and the central line of the corresponding target beam. And the fourth sub-figure refers to the change of the reference satellite’s mass.
Table 2. The relations of the 3 satellites.
reference satellite’s orbits and the central line of the corresponding target beams are no more than 100 km, which means that most of the target beams in addition to the passed beams along the transfer orbit could all be observed by the formation. The consumed fuel by the reference satellite is lower than 160 kg, which means 40 kg left could be used to guarantee the fulfillment of the formation transfer and maintenance along the overall task.
4.3. Formation Transfer and Maintenance
For the formation transfer and maintenance, the MICM described in section 3.4 is employed. The time points, Δvs of all maneuvers are optimized by the DE algorithm and J2 perturbed Lambert problem solving algorithm. Therefore, the satellites can fly following the reference satellite and the formation can be maintained. The orbit of the reference satellite is shown in Figure 9, in which the red pluses refer to the maneuvers.
The changes of the reference satellite’s elements are shown in Figure 10, in which the semi-major and the orbital eccentricity are firstly raised by the initial transfer, and then oscillate in the following transfers. The orbital inclination nearly unchanged to save fuel. And the right ascension node changes from about 225 deg to 100 deg due to the J2 perturbation.
The consumed fuel for all the maneuvers are shown in Figure 11. For every transfer, the velocity increment and the consumed fuel are minimized by the DE algorithm. There are 3 transfers whose consumed fuels are much more than others. The first one is for the transfer from the initial circular orbit to the elliptical orbit, the second one is for the formation changing, and the third one is for the final flying by to monitor the last beam in order to make the total time to be the least while consuming nearly all the fuel of each satellite.
Figure 9. The orbit of the reference satellite (The red pluses refer to the maneuvers).
Figure 10. The changing of the reference satellite’s elements.
Figure 11. The consumed fuels for the satellites’ maneuvers.
4.4. Monitoring Profit Computation
Incorporating all the maneuvers of the 3 satellites, the monitoring profits (observation values) shown in Figure 12 are computed using Equation (9) by propagating the orbits. Figure 12 also shows that all the monitoring profits of the GEO satellites’ beams are larger than the threshold of 1 × 107 km, which fulfills the requirement of the overall monitoring task.
Figure 12. The monitoring profits of all the GEO satellite’s beams (The red line refers to the value of 1 × 107 km2). The distance between every two pluses for every beam refers to the incremented monitoring profit for every monitoring mission.
5. Competition Results
We have obtained the fourth prize of CTOC9. The total flying time is 26.9611 days, and the total consumed fuel is 549.3799 kg, as shown in Table 3. The reasons for the discrepancies of the results obtained in this paper and that obtained by the previous research teams at NUAA, NUDT and Tsinghua University may be the amount of calculation is much larger than them, thus missed the optimal solution in a limited duration of the competition.
The initial parameters, the number of maneuvers, and the final consumed fuels of the 3 satellites are shown in Table 4.
The final observation values are shown in Table 5.
This paper presented the methods and results for the problem A of CTOC9, the trajectory design and optimization problem for the low earth orbit (LEO) satellites in formation to observe the geostationary orbit (GEO) satellites’ beams. Two methods have been proposed to respectively solve the sequence searching and formation transfer & maintenance problems. The first method is the IEMM and the second one is MICM. Compared with other groups with better results, the methods in this paper may be more powerful in solving similar or more complex multi-objective trajectory design problems, because most of the methods used by other groups employed the regressive orbit to solve the special problem A of CTOC9. If the beams pointing to the surface of the earth to be
Table 3. The results of the objective values. (To make a clear comparison, we also presented the results of the other teams in the table.The contents in bold is our result.)
Table 4. Initial parameters and the final consumed fuels.
Table 5. The final monitoring profits.
monitored are not fixed, the regressive orbit may not be an effective candidate tool to solve the problem. But we have attempted to propose general methods which need not any special characteristics of the monitor or the target. The proposed methods may have a wider application area to solve more multi-objective trajectory design and optimization problems.
We are very grateful to the 9th China trajectory optimization competition organizers for the interesting problem, challenge, and fruitful discussions. This work was supported by the Youth Fund Project of the Space Systems Department (No. 2017SY26B0001), the Funds of Natural Science of China (No. 61673390, and No. 61503414 ), and the Fund of The State Key Laboratory of Astronautic Dynamics (No. 2016ADL-DW0202).
 Asgarimehr, M. and Hossainali, M.M. (2014) Optimization of Geosynchronous Satellite Constellation for Independent Regional Navigation and Positioning in Middle East Region. Acta Astronautica, 104, 147-158.
 Zhu, K., Jiang, F., Li, J., et al. (2009) Trajectory Optimization of Multi-Asteroids Exploration with Low Thrust. Transactions of the Japan Society for Aeronautical and Space Sciences, 175, 47-54.
 Storn, R. and Price, K. (1997) Differential Evolution—A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. Journal of Global Optimization, 11, 341-359.