Trans. Japan Soc. Aero. Space Sci. Vol. 46, No. 154, pp. 271–274, 2004 Technical Note A Low-Cost Orbit Determination Method for Mobile Communication Satellites By Sangwoo CHO, Jinho B AE and Joohwan CHUN Scientific Computing Laboratory, Korea Advanced Institute of Science and Technology, South Korea (Received January 14th, 2003) Key Words: 1. Orbit Determination, Simulated Annealing Introduction The problem of satellite orbit determination using Earth’s magnetic field measurements can be approached from various directions. Shorshi and Bar-Itzhack proposed an autonomous navigation method using the extended Kalman filter (EKF).1) Psiaki suggested a batch filter scheme for autonomous orbit and magnetic field determination,2) where he used the Gauss-Newton technique for the minimization of a certain least-squares cost function. The EKF provides a suboptimal solution of the nonlinear estimation problem. However, an initial state far from the true one may cause a divergence problem due to the linearization of a highly nonlinear system equation. With the optimization techniques using (local) gradient information, convergence to the global minimum is dependent upon an educated guess of the initial state. With an inappropriate initial state, the iteration may converge to a local minimum.3) On the other hand, (global) optimization methods such as the simulated annealing technique have the ability to escape from local minima and to converge to the global minimum.4,5) This paper presents the generalized simulated annealing technique to find the initial state for the autonomous orbit determination using an onboard magnetometer. It uses the magnetometer measurement only, and does not require additional information, such as the attitude of the satellite, ground station intervention, or the Global Positioning System (GPS) signal. 2. Algorithm Description Our choice of state-vector elements is the classical orbital elements: x ¼ ½a e i ! M 0
T ; ð1Þ where a is the semi-major axis, e is the eccentricity, i is the inclination, is the longitude of the ascending node, ! is the argument of periapsis, and M0 is the mean anomaly at epoch. Because the mean anomaly varies uniformly in time, the mean anomaly at epoch is chosen as the sixth element instead of the true anomaly , the eccentric anomaly, or the time of periapsis passage. This simple linear property reÓ 2004 The Japan Society for Aeronautical and Space Sciences duces the complexity of our equation for orbital dynamics.6) In this paper, the J2 effect and drag influence are ignored for the sake of simple formulation. Therefore, the backward state-vector propagation from the current time t to the epoch time t0 depends only on the mean anomaly M, and is formulated by the continuous time equation: M ¼ M0 þ nðt  t0 Þ, where n is the mean motion. The true anomaly  can be obtained by a function of the eccentric anomaly which is determined by solving Kepler’s equation numerically.7) The Cartesian position vector p of a satellite in the Earthcentered Earth-fixed (ECEF) frame is obtained by p ¼ T2 ðt  t0 ÞT1 ði; ; !Þpw ; ð2Þ where pw is the Cartesian position vector in the perifocal frame, T1 is the transformation matrix from the perifocal frame to the Earth-centered inertial (ECI) frame, and T2 is the transformation matrix from ECI frame to ECEF frame. The spherical position ps ¼ ½r  
T of the vector p is the nonlinear function of state elements. Here, r is the geocentric distance,  is the colatitude, and  is the east longitude from the Greenwich meridian. If we consider the discrete measurements, then the dynamics of the spherical position become ½rk k k
T ¼ f ðx; tk Þ; ð3Þ where f is a nonlinear vector function and tk is the measurement time. The elements of the spherical position vector are used to determine the geomagnetic field vector bigrf ðr; ; Þ, obtained from the International Geomagnetic Reference Field (IGRF) model. The system measurement is the magnitude of the geomagnetic field vector bm;k obtained from the onboard magnetometer at time tk . We choose the magnitude instead of the individual components since the magnitude of the geomagnetic field is independent of the attitude of the satellite. The cost function JðxÞ we have chosen is the sum of the squares of the difference between the 2-norm of the measured geomagnetic field vector and the IGRF model: JðxÞ ¼ N X 1 ðkbm;k k2  kbigrf ðx; tk Þk2 Þ2 ;
2 N maxk kbm;k k2 k¼1 ð4Þ 272 Trans. Japan Soc. Aero. Space Sci. where N is the number of the measurements and the denominator Nðmaxk kbm;k k2 Þ2 is the normalizing factor that reduces the influence of the altitude on the cost function magnitude. To minimize the cost function, we use the generalized simulated annealing technique, which finds the global minimum along iterations with a fixed control parameter and a random walk. After each iteration, the control parameter increases and the size of the random walk decreases. The new state caused by the random walk can be accepted or rejected by the Metropolis criteria. If the value of cost function decreases at the new state, the state is accepted. The new state with an incremental value of cost function is accepted with a probability p ¼ expð J01
JÞ, which is able to escape from local minima. The initial state x0 can be chosen arbitrarily anywhere in the bounded search space  R6 . Since we consider LEO satellites with direct and polar orbital motions, the search space is al < a < au , 0  e  1, 0  i  =2, and 0  , !, M0 < 2. Here, the lower and upper bounds of the semi-major axis are al ¼ 6900 km and au ¼ 7900 km, respectively, because the altitude of the LEO satellite is between 500 km and 1500 km. The simulated annealing technique is summarized below: Vol. 46, No. 154 Algorithm (Simulated Annealing) Choose x0 Set t0 Compute J0 = J (x0 ) using eq. (4) If |J0 | < ǫ, then break End If For k1 = 1, 2, . . . , κ1 k2 = 1 While (|J0 | > ǫ) and (k2 < κ2 ) I: Generate r ⊙ u Set x1 = x0 + r ⊙ u If x1 ∈ Ŵ, then set J1 = J (x1 ) and J = J1 − J0 Else, goto step I End If k2 = k2 + 1 If J1 ≤ J0 , then set x0 = x1 and J0 = J1 If |J0 | < ǫ, then break Else, goto step I End If Else, set p = exp(−β J0−1 J ) If V ≥ P, then goto step I Else, set x0 = x1 and J0 = J1 , then goto step I End If End If End While r = γ r β = β/γ End For Table 1. Summary of the true state and estimates for each cases. Case A1 A2 A3 I1 I2 I3 O1 O2 O3 a (km) e i ( ) ( ) ! ( ) M0 ( ) J ð105 Þ 1.515 true state 7000.000 0.05000 60.000 120.000 70.000 50.000 estimates 6999.379 0.05005 59.301 119.256 69.571 50.521 1.129 true state 7400.000 0.05000 60.000 120.000 70.000 50.000 1.767 estimates 7400.624 0.05048 60.247 120.029 69.684 50.403 1.299 true state 7800.000 0.05000 60.000 120.000 70.000 50.000 2.177 estimates 7801.165 0.05027 60.115 120.902 69.564 50.037 0.495 true state 7400.000 0.05000 15.000 120.000 70.000 50.000 1.012 estimates 7400.329 0.04974 14.949 120.032 70.186 49.713 0.943 true state 7400.000 0.05000 45.000 120.000 70.000 50.000 0.902 estimates 7398.131 0.04991 45.419 118.563 70.762 49.759 0.527 true state 7400.000 0.05000 85.000 120.000 70.000 50.000 0.763 estimates 7398.094 0.04995 85.585 120.555 70.353 49.417 0.595 0.743 true state 7000.000 0.05000 60.000 100.000 70.000 50.000 estimates 7001.061 0.05026 60.227 99.643 70.704 49.534 0.563 true state 7000.000 0.05000 60.000 220.000 70.000 50.000 1.327 estimates 6999.373 0.05080 60.035 221.093 68.874 50.961 1.084 true state 7000.000 0.05000 60.000 340.000 70.000 50.000 1.536 estimates 6998.957 0.04982 60.208 341.187 70.523 49.390 1.280 In the algorithm,  is a threshold,
r is a step size, u is a unit random direction whose element is drawn from uniform distribution on [1, 1],  denotes the elementwise product, is a control parameter, V is a uniform variate on ½0; 1
, 1 is the number of control steps, 2 is the number of iterations per a fixed control parameter, and is a reduction factor of less than one. 3. Simulation Results To verify that the algorithm can be applied for any low Earth orbits, we consider nine low Earth orbits whose state elements are listed in Table 1. We classify the nine orbits into three cases according to element variation. At case A, we vary the semi-major axis for three subcases. Cases I and O are considered by varying the inclination and the longitude of the ascending node, respectively. For each subcase, the Feb. 2004 S. CHO et al.: A Low-Cost Orbit Determination Method for Mobile Communication Satellites 273 Fig. 1. The two-dimensional cost function for case A3.  indicates the true state. 7.9 0.8 e a (× 103 km) 0.6 7.4 0.4 0.2 6.9 0 0 5 10 15 20 25 30 90 0 5 10 15 20 25 30 360 Ω (°) i (°) 45 180 0 0 0 5 10 15 20 25 30 360 0 5 10 15 20 25 30 360 ω (°) M0 (°) 180 180 0 0 0 5 10 15 20 Control Step (κ1) 25 30 0 5 10 15 20 Control Step (κ1) 25 30 Fig. 2. The convergence behavior of all state elements. first and second rows are the true and estimated state values, respectively. In our simulations, it is assumed that the Earth is rotating eastward at a rate of 7:29211593  105 rad/sec. The time interval between measurements is 10 minutes, and the number of measurements N is 30. The 10th degree, 10th order IGRF model is used to determine the geomagnetic field. White Gaussian measurement noises with zero-mean and  ¼ 100, 150, and 200 nT standard deviations are added to each component of the measured geomagnetic field vector according to the semi-major axes (a ¼ 6900, 7400, and 7900 km) of the orbit, respectively. Figure 1 shows the selected two-dimensional cost functions for case A3. The  mark denotes the true values. The two-dimensional cost functions depicted in Fig. 1 ascertain that several local minima occur in the search space. In Fig. 1(a), the variation of the cost function along the semimajor axis (x axis) is about twice as much as other state elements because the magnitude of the geomagnetic field rapidly decreases as the altitude increases. The annealing schedule used for all cases is composed of thepcontrol paffiffiffi rameter ¼ 0:45 with reduction factor ¼ 1= 2, the step size
r ¼ ½300 km 0:3 27 108 108 108
T , the control step 1 ¼ 28 and the number of iterations 2 ¼ 500. The threshold  is 5  106 . We start the iteration with the initial state x0 ¼ ½7400 km 0:5 45 180 180 180
T . It is the center of the search space. The major computational load comes from the calculation of the cost function. For each iteration, 117,840 flops (floating point multiplications and additions) and 238 trigonometric calculations are needed. The estimates of each case are given in Table 1. The last column of Table 1 means the discrepancy between the true and the estimated states. Although function value at the final stage is smaller than the true state, it does not mean that estimate is more accurate because of measurement noise and uncertainty of the IGRF model. Figure 2 depicts the convergence behavior of each state element for case I3. The estimate converges near the global minimum prior to the 5th control step. Note that all elements have the same convergence rate, since the relative element magnitudes for the step size to the search space are equal. For case A3, the convergence behavior of the proposed method is displayed in Trans. Japan Soc. Aero. Space Sci. 360 270 270 Vol. 46, No. 154 M0 (°) 360 0 M (°) 274 180 180 90 90 0 0 0 90 180 270 360 0 90 Ω (°) 180 270 360 ω (°) Fig. 3. The visual convergence behavior of the state. Fig. 3. The dots denote the trajectory states of the accepted steps in the iteration stages. Note that the estimate near local minima does not get trapped. 4. Conclusions We have proposed a low-cost, robust and autonomous initial orbit determination method for an LEO-satellite using the simulated annealing technique. Our method has the following advantageous features: First, it only uses an inexpensive magnetometer. Second, it does not require attitude information, ground intervention, or a GPS signal. Finally, our method does not need a-priori information on the orbital elements, because the simulated annealing technique finds the global minimum. Acknowledgments This work was supported, in part, under grants from the Ministry of Science and Technology managed by SaTReC, MICROS and by KOSEF (R01-2003-000-10829-0). References 1) Shorshi, G. and Bar-Itzhack, I.: Satellite Autonomous Navigation Based on Magnetic Field Measurements, J. Guid. Control Dynam., 18 (1995), pp. 843–850. 2) Psiaki, L.: Autonomous Orbit and Magnetic Field Determination Using Magnetometer and Star Sensor Data, J. Guid. Control Dynam., 18 (1995), pp. 584–592. 3) Press, H., Teukolsky, A., Vetterling, T. and Flannery, P.: Numerical Recipes in C, Cambridge University Press, New York, 1992. 4) Bohachevsky, O., Johnson, E. and Stein, L.: Generalized Simulated Annealing for Function Optimization, Technometrics, 28 (1986), pp. 209–217. 5) Corana, A., Marchesi, M., Martini, C. and Ridella, S.: Minimizing Multimodal Functions of Continuous Variables with the Simulated Annealing Algorithm, ACM Trans. Math. Software, 13 (1987), pp. 262–280. 6) Chobotv, A.: Orbital Mechanics, AIAA Inc., Reston, Virginia, 1996. 7) Wertz, R.: Spacecraft Attitude Determination and Control, Reidel, Boston, MA, 1978.