Comparison of Global Optimization Methods for Insertion ... - ijass

Paper Int’l J. of Aeronautical & Space Sci. 15(3), 267–280 (2014) DOI:10.5139/IJASS.2014.15.3.267

Comparison of Global Optimization Methods for Insertion Maneuver into Earth-Moon L2 Quasi-Halo Orbit Considering Collision Avoidance Sang-Cherl Lee*, Hae-Dong Kim*, Do-Chul Yang*, Dong-Hyun Cho* and Jeong-Heum Im* Korea Aerospace Research Institute, Daejeon, Republic of Korea

Tae-Soo No** Department of Aerospace Engineering, Chonbuk National University, Jeonju, Republic of Korea

Seungkeun Kim*** and Jinyoung Suk**** Department of Aerospace Engineering, Chungnam National University, Daejeon, Republic of Korea

Abstract A spacecraft placed in an Earth-Moon L2 quasi-halo orbit can maintain constant communication between the Earth and the far side of the Moon. This quasi-halo orbit could be used to establish a lunar space station and serve as a gateway to explore the solar system. For a mission in an Earth-Moon L2 quasi-halo orbit, a spacecraft would have to be transferred from the Earth to the vicinity of the Earth-Moon L2 point, then inserted into the Earth-Moon L2 quasi-halo orbit. Unlike the near Earth case, this orbit is essentially very unstable due to mutually perturbing gravitational attractions by the Earth, the Moon and the Sun. In this paper, an insertion maneuver of a spacecraft into an Earth-Moon L2 quasi-halo orbit was investigated using the global optimization algorithm, including simulated annealing, genetic algorithm and pattern search method with collision avoidance taken into consideration. The result shows that the spacecraft can maintain its own position in the Earth-Moon L2 quasi-halo orbit and avoid collisions with threatening objects. Key words: Circular Restricted Three Body Problem, Quasi-Halo, Collision Avoidance, Space Objects, Simulated Annealing, Genetic Algorithm, Pattern Search, Maneuver Planning, Orbit Insertion Maneuver

1. Introduction In 1760, Leonard Euler formulated the restricted three body problem in a rotating frame handmade prediction on the premise that three collinear equilibrium points (L1, L2, L3) exist. In his 1772 essay “essaisur le probleme des Trois Corps”, Lagrange confirmed Euler’s prediction and solved for two additional equilibrium points (L4, L5) which form an equilateral triangle with the primaries. Henri Poincare later demonstrated in his “Methodes Nouvelles de la Mechanique Celeste” that there could be no such analytic quantities in positions, velocities and mass ratio. He concluded that the

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/bync/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received: April 2, 2014 Revised : August 2, 2014 Accepted: September 15, 2014 Copyright ⓒ The Korean Society for Aeronautical & Space Sciences

three body problem could not be solved using algebraic formulas and integrals[1]. The Sun-Earth and Earth-Moon Lagrange points as well as their Halo orbits are shown in Fig.1. The triangular equilibrium points, L4 and L5, could be used as parking regions because no station keeping is necessary on those points. On the other hand, the collinear points L1, L2 and L3 can be selected for space observation and planetary exploration missions. They can provide good observation orbits toward the Sun in the Sun-Earth system during Sun observation missions. A spacecraft for the space observation mission orbiting around the Sun-Earth L2 point can provide constant observation geometry with half of the entire celestial

* Senior Researcher Professor ** *** Professor **** Professor, Corresponding author : [email protected]

267

http://ijass.org pISSN: 2093-274x eISSN: 2093-2480

Int’l J. of Aeronautical & Space Sci. 15(3), 267–280 (2014)

They provided the first ever three-dimensional information on how energetic particle acceleration takes place near the Moon’s orbit in the distant magnetosphere and in the solar wind[4]. Roughly 50 man-made objects had been placed in the lunar orbit since the dawn of space exploration up to mid2005. At least 34 of them have crashed onto the Moon or been maneuvered out of the lunar orbit. Several of the remaining spacecrafts and propulsion units are also believed to have reached the lunar surface[5]. Although the number of manmade debris or asteroids in the lunar orbit is still small, a collision avoidance maneuver system would be useful to reduce the probability of colliding with uncontrolled threatening objects from the outer space. The scope of this paper focuses on an insertion maneuver to the quasi-halo orbit using global optimization methods including simulated annealing, genetic algorithm and pattern search. Prior works related to this research are

sphere available at any given time. Additionally, an orbit around the Sun-Earth L2 point is useful for non-cryogenic missions that require great thermal stability, suitable for highly precise visible light telescopes. Periodic orbits around the Earth-Moon L2 point can be used to establish a permanent communication link between the Earth and the far side of the Moon as suggested by A.C. Clark in 1950 and R. Farquhar in 1966. An orbit around a Lagrange point provides the ability for Earth transfer, return trajectory and interplanetary transport. Recent studies have also demonstrated that even formation flying is possible using the Lagrange points[2]. Table 1 lists past, present and future spacecraft missions on the Lagrange points[3]. The ARTEMIS mission was the first to navigate and perform station keeping operations on the Earth-Moon L1 and L2 Lagrange points. ARTEMIS P1 and P2 spacecrafts were used for simultaneous measurements of particles and electromagnetic fields from the two locations.

Fig.1. Lagrange Point Orbits and Halo orbits Fig. 1. Lagrange Point Orbits and Halo orbits

Fig.1. Lagrange Point Orbitson andthe HaloLagrange orbits Table 1. Spacecraft missions Points

Table 1. Spacecraft Missions on the Lagrange Points

Table 1. Spacecraft missions on the Lagrange Points

Point

L1

Status Mission Status Point Completed Mission Completed On Operation

L4 L5

Earth-Moon

Earth-Sun ISEE-3, Genesis

Earth-Moon

ISEE-3, Genesis

Not reported

Not reported

WIND, SOHO, ACE Not reported WIND, SOHO, ACE Not reported LISA Pathfinder, KUAFU, Deep Space Climate Planned Not reported LISA Pathfinder, KUAFU, Solar-C Deep Space Climate Observatory, Not reported Planned Observatory, Solar-C Mission WIND, Chang'e2 ARTEMIS P1/P2 Completed Mission WIND, Chang'e2 ARTEMIS P1/P2 Completed On Operation WMAP, Planck Surveyor, Herschel Space Observatory Not reported WMAP, Planck Surveyor, Herschel Space On Operation L2 JWST,Gaia, Terrestrial Observatory Planet Finder, Darwin, Wide Field Not reported Planned Exploration Gateway Platform Infrared Survey Telescope, Euclid JWST,Gaia, Terrestrial Planet Finder, Darwin, Wide Exploration Gateway Platform Planned Field InfraredSTEREO Survey Telescope, Euclid On Operation A Not reported L4 On Operation STEREO A Not reported On Operation STEREO B Not reported L5 On Operation STEREO B Not reported L1

L2

Earth-Sun

On Operation

50 man-made objects had been placed in the lunar orbitlunar since the dawnfrom of space exploration up to mid- exploration Roughly 50 Roughly man-made objects had been placed in the orbit the dawn of space

DOI:10.5139/IJASS.2014.15.3.267 268 2005. At least 34 of them have crashed onto the Moon or been maneuvered out of the lunar orbit. Several of the

to mid of 2005. At least 34 of them have crashed onto the Moon or been maneuvered out of the lunar remaining spacecrafts and propulsion units are also believed to have reached the lunar surface[5]. Although the

orbit, and several of the remaining spacecrafts and propulsion units are also believed to have reached 4

the lunar surface[5]. Although the number of man-made debris or asteroids in the lunar orbit is small,

Sang-Cherl Lee Comparison of Global Optimization Methods for Insertion Maneuver into Earth-Moon L2 ...

paper compares three global optimization algorithms. The simulated annealing method uses stochastic iteration or meta-heuristics which is derivative-free and can handle problems without constraints or boundary conditions. This method can find a global solution given a slowenough cooling schedule. The genetic algorithm is based on evolutionary computation and also uses stochastic iteration without gradients. The pattern search method uses deterministic iteration, which is known to be more robust and efficient in computation time. The simulated annealing method is a generic probabilistic meta-heuristic method proposed by Kirkpatrick, Gelett and Vecchiinin in 1983 to find a global minimum of a cost function. The method originated from thermodynamics and models the cooling process of a metal, with the initial state fed into a neighbor-generating function to create a random neighboring state. Any downhill step is accepted for minimizing functions, and the process is repeated from a new point afterward. If an uphill step is accepted, the simulated annealing can escape from a local minimum. This uphill decision is made using the Metropolis criteria. As optimization progresses, the length of the steps decreases and the solution gets closer to a global optimum[9, 12]. The Genetic algorithm introduced in 1975 by John Holland is a meta-heuristic global search method for solving optimization problems. The genetic algorithm creates a population and applies genetic operators such as mutation and crossover for finding the best solution. The three most important requirements for using the genetic algorithms are the definition of an objective function, the implementation of a genetic representation, and the implementation of genetic operators. The genetic algorithm is an iterative scheme where the population is modified using the best features of the ‘genes’ from previous generations. The selection, crossover and mutation operators are applied to identify the best solution[6, 12, 13]. Pattern search is an optimization process also known as the direct-search, derivative-free and black-box method. It can be used on problems that are not continuous or differentiable. It was named by Hooke and Jeeves in 1961 and does not require an optimized function gradient. Fermi and Metropolis at the Los Alamos nuclear laboratory later developed a more advanced method with an adjustable mesh size. Pattern search evaluates the objective function at the lattice point centered on the current iteration to find the best objective function value then defines a new iteration at the resulting lattice point. The algorithm also expands or contracts the lattice depending on whether a more appropriate solution has been found. The generation point

described as follows. Santos et al. studied the transfer trajectories between two coplanar orbits using several impulses, using genetic algorithms to find solutions that minimize fuel consumption[6]. Paluszek et al. used a downhill simplex, genetic algorithm and simulated annealing to investigate the problem of finding a minimum time trajectory from the Earth to the Mars orbit using a fixed thrust electric propulsion system[7]. Wang et al. implemented the pattern search method to reduce deviations between calculations performed by the onboard calculation and the ground elaboration orbit models with the goal of improving the accuracy of on-board attitude determination for satellites[8]. Vries developed an efficient method of satellite maneuver optimization based on a Monte Carlo approach in combination with simulated annealing[9]. Lastly, Kim et al. studied a genetic algorithm optimization of the lunar transfer trajectory using a weak stability boundary[10]. This paper investigates an insertion maneuver of a spacecraft into the Earth-Moon L2 quasi-halo orbit by using simulated annealing, genetic algorithm and pattern search method to decrease the probability of collision with threatening uncontrolled objects. The fitness function is used to minimize the distance between the reference L2 quasi-halo orbit and a maneuvering spacecraft. The STK/ Astrogator, STK/AdvCat and Matlab global optimization toolbox are utilized as simulation tools[12, 22]. Unlike the aforementioned works, this study aims to compare global optimization methods for an insertion maneuver into the Earth-Moon L2 quasi-halo orbit with collision avoidance taken into consideration. The most valuable benefit of this approach is that it can get a periodic orbit trajectory near the unstable Lagrange points where the gravitational attractions of the Earth and Moon equalize. The rest of this paper is organized as follows. Section 2 contains a brief introduction of the optimization algorithm used in this paper. The circular restricted three body problem is described in Section 3. Section 4 presents the maneuver plan for an insertion into the Earth-Moon L2 quasi-halo orbit using global optimization. Sections 5 and 6 describe simulation results and the conclusion.

2. Global Optimization Algorithm The non-linear solver based on the gradient method quickly converges to a local minimum near the initial guess, but may not find a global solution. In contrast, the global search method is able to search more effectively within an entire solution space and find a global minimum[12]. This

269

http://ijass.org

distance from the Moon to the spacecraft, respectively. μand (1-μ) represent the non-dimensional mass of the distance from the Moon to the spacecraft, respectively. μand to (1-μ) therespectively. non-dimensional the distance from the Moon the represent spacecraft, μand mass (1-μ) of represent the non-dimensional

Moon and the non-dimensional mass of the Earth, respectively.F� , F� andF� denote external perturbation forces Moon and the non-dimensional mass Moon of the and Earth, andFof� denote external perturbation therespectively.F non-dimensional the Earth, respectively.F � , F�mass � , Fforces � andF� denote external perturb

on the x and y axes and the thrust of the spacecraft respectively. on the x and y axes and the thrust of the on spacecraft the x and yrespectively. axes and the thrust of the spacecraft respectively.

On one of the three collinear liberation points (x� � in the planar case, it is assumed that the external On one of the three collinear liberation (x�three � in the planarliberation case, it is assumed theplanar external On onepoints of the collinear points (x� � that in the case, it is assumed that

and thrust are zero. Then, Eqs.(1)~(3) can be written as: Int’l J. of Aeronautical & Space Sci. 15(3), perturbation 267–280 (2014) perturbation and thrust are zero. Then,perturbation Eqs.(1)~(3)and can thrust be written as: Then, Eqs.(1)~(3) can be written as: are zero.

in a pattern is determined by polling, with the maximum number of points generated at each step. The efficiency of the solution can be improved by controlling the polling order ��μ� μ of the points[3, 12, 14]. where � � |��μ�� |� μ � . where � � |��μ| � � |��μ| � . � �μ| � ��μ|

The eigenvalue of Eqs.(4)~(6) is The eigenvalue of Eqs.(4)~(6) is

3. Circular Restricted Three Body Problem

x� � �y� � �� x x� � �y� � �� x

x� � �y� � �� x

y� � ��x� � �� y y� � ��x� � �� y z� � ��z z� � ��z

��μ�

where � � |� where

� �μ|

�

� |�

y� � ��x� � �� y μ

� ��μ|

�

..

z� � ��z

(4) (4)

(5) (5)

(5)

(6) (6) (6)

The eigenvalue of Eqs.(4)~(6) is The eigenvalue of Eqs.(4)~(6) is 8 8 � � � � �� β� �� β � � �β � λ �� β �� β �β ��β �βββ��β� � λ�� λλ�� β � � �β λλ�� β � � �β � �� β � �β � � β��

8 (7) (7) (7) (7) λ The general three body problem posits that three bodies � β� � (7) � (7) �� β � �β�� λ�� β (7) �β λ�� β � � � with arbitrary masses attract one another according to the � � � �� β (7) λ�� β� � β� � �� β� � �β � β��β (7) � β�� β� � β� (7) λ�� λ�β � �� universal law of gravitation. The gravitational threeλ�� body � � �β (7) � � �� β � � � � � � � � �β � β (8) λ �� β �� β �βββ��β�� (8) λ�� βββ��β�� λλλ�� (8) (8) � λ � �β � �� (8) � �� β � � �� (8) problem requires 18 integrals of motion, but the number �β � �� β � � β λ �� (8) λ�� β� � � �β� � � β� � λ�� β� � �β � � ��β� of integrals could be reduced to 8 by using conserved � � �� (8) λ�� β �� β � β �β �� (9) β� � � �� β�� (8) λ�� λ�� √√ �� β β� � (8) 3. Circular Restricted Three Body Problem λ�� (9) �� λλλ� cc�cccβ�� β�� β�� β�� √ quantities; six of these scalar integrals are determined λ �� λ�� (9) (8) (9) (9) √ � �� (9) 3. Circular Restricted Three Body Problem √ λ�� √c√c �� λ � �� c (9) λ�� √ 3. Circular Restricted Three Body Problem �� from conservation of linear momentum, three neral three body problem posits thatthe three bodies with arbitrary masses attractββ�βone according to � �� another ��from and β the ��c ��c ��c �� where � � � �� and ��c where ��c � � �� where neral that with arbitrary attract another according to λ�� ...�√ (9) � and � and β� � ��c � �� and � ��c �� ... c� where β�β�� c ��c �� c �� where ��c λ��βββ��ββ�� c � (9) where β � and β � ��c � ��c � �� . . where � ��and where √ neral three three body body problem problem posits posits that three three bodies bodies with angular arbitrary masses masses attractβone one another according to �� λ �� c (9) β� � �(9) � � � and � � ��c � �� where where β √ ��c conservation of total momentum, and one integral � �� and �β� � ��c � ��c � �� . � � λ�� √ c� � � Circular Restricted Body Problem versal law of gravitation. 3. The gravitational three Three body problem requires 18general integrals of motion, but as: the � is written solution versal gravitational three problem requires 18 integrals motion, the conservation of energy[15]. The � general is written as: �The The general solution iswritten written as: � ��c general solution is as: � �of �� and βbut � the �� solution . solution where ��c The general solution is written as: versal law law of of gravitation. gravitation. The Thefrom gravitational three body body problem 18� integrals of motion, but the The general solution is written The general solution is written as: �requires � �βThe and β�β�� .as: where β�where The is is written as:as: ��c � and �� c � � β ��c . general ��c general solution is written as:� �� The general solution written �� The � � � � and β � � ��c � �� . where βone � ��c three body problem posits The that bodies with arbitrary masses attract another according to three bodyconserved problem can be simplified into a circular rneral of integrals could be reduced to 8three by using quantities; six of these scalar integrals are � � λ�λ� �λ� � �λ� �� α e � α e � A c��ω � � φ� (10) λ� �λ� of of these scalar are � � c��ω �� α�α�α��αee�eλ� eλ��ααα��αee�e�λ� e�λ��A A �� φ� φ� (10) (10) (10) �� of integrals integrals could could be be reduced reduced to to 88 by by using using conserved conserved quantities; quantities; six six ofgeneral these solution scalar integrals integrals are � c��ω �� c��ω�� (10) The is written (10) �φ� �� α��αe�λ�eλ� �� α�αe��λ� �� A�Ac��ω � �� φ�φ� �� Theproblem general solution issolution written as:writtenof �� restricted three body defined as a system two as: �� α�� eλ�� α�� e�λ�� AA A��c��ω �φ� φ� (10) e�λ� The general is as: � c��ω�� c��ω versalfrom lawthe of conservation gravitation. The gravitational threethree bodyfrom problem requiresThe 18 general integrals of motion, but the is written as: ined of linear momentum, the conservation of total solution angular momentum, λ�λ� �λ� ned from the conservation of linear momentum, three from the conservation of total angular momentum, �λ� α e � k α e � k A ��ω � � φ� (11) �� k λ� �λ� �� ω �� in a circular their barycenters, and momentum, a third λ� λ�� ned from the conservation of bodies linear momentum, threeorbit fromabout the conservation of total angular eλ� � kαα�α�� αee�eA e�λ� �kkk��kA A �� φ� φ�λ� λ� (11) �� k �k e� kk�λ� � �� (11) �� λ� � λ� �A �eα � � ��ω �λ�� α α� e� ��ω � �φ� (10) α�α kk ��ω �λ� � �� (11) ��λ� � c��ω αφ� k �kα(10) e�λ� � (11) k(11) ��ω � �� φ�φ� �� k �α �� α�� e � α �� e�� φ� ��α ��ee � �φ� �(10) �kA��A �� e�� e�λ� � A�� ω � φ� (11) �� αA�k ��k �� αe� e � � � �� k eλ�� ec��ω Aλ� c��ω φ�k��A�� λ� �λ�α� � �� k � ��α� e � ��ω � A� � of integrals be reduced 8 by negligible using conserved quantities; six of these scalar integrals �are � � � α� �� er integral from could the conservation oftoenergy[15]. � � � α e � α e c��ω � � φ� (10) body with mass. In the Earth-Moon barycentric � � � �� ee integral from the conservation of energy[15]. �� A c��ω � � ψ� (12) integral from the conservation of energy[15]. � � c��ω �λ� �� ψ� ψ�α e�λ� � k A ��ω � � φ� (12) �� A�AA �� (12) � c��ω �� kψ� (11) �� λ� �k �� c��ω (12) �� ψ� � αk � e��A � �� A��Ac��ω ψ�(11) ��c��ω eλ�� k � k� α�� α�� e�A��λ� � φ� �� (11) �the � � �k theabarycenter of the two celestial bodies becomes � � �� ψ� (12) �� ψ� (12) (12) � α�� ω e� � k� eλ��λ� kα��Ae��λ� �� φ� �� inedbody fromproblem the conservation offrame, linear momentum, three from the conservation of total angular momentum, � c��ω � c��ω � �A �� αα�� ω ree can be simplified into circular restricted three body problem defined as�� a system of e � k � k A ��ω � � φ� (11) �� k � � � � � � �� ee body problem can be simplified into a circular restricted three body problem defined as a system of � � ��ω ��ω � � �� ee body problem can be simplified into a circular restricted three body problem defined as a system of ��ω � ��ω origin. The x axis is parallel to the line between the Earth �� ,�� k � and ��ω�� where ��ω ��ω α�� A and A��A can bebe calculated from the initial condition, and φφ �� ω ��ω ��ω ��ω �� ω �.. .α ��ω �� AA c��ω ψ�be (12) �� ,k � wherekkkk��k�� α�,,,�,,ψ� A and can calculated fromthe theinitial initial condition, and �� where αααα��α� ,,� ,A and A can calculated from condition, and φ ��ω ��ω �ω�ω �ω �� ω ��ω �� where kcalculated � , k � the � . α� �� ,�� c��ω �αwhere (12) where ��ω α��α and be calculated from initial condition, and φφ from where kkkk�� AA can be the initial condition, φ�ω edies integral from theorbit conservation of barycenters, energy[15]. and a third body with negligible �and �ω � �ω �..A c��ω �A ψ� (12) � k . α. �α,from α, �α, �A A�A can becan cond �ω �ω �� and in a circular about their In �k�,�,,k � ,�,,Earth-Moon �� .�where �the �� ψ� �A�and �ω where k�mass. �� A�� A�� and A�,�� can becalculated calculated the initial condition, and , can k� , the and can calculated from theinitial initial con and �ω � �� ω �ω�� , α � �� c��ω � Earth-Moon �� where � � �from � A �� A (12) �ω�ω �ω z be Moon, and is directed from the Earth toward the Moon. The �the �� dies in a circular orbit about their barycenters, and a third body with negligible mass. In � �� ω �ω �ω � Earth-Moon �� dies in a circular orbit about their barycenters, and a third body with negligible mass. In the � � ��ω ��ω � �� ω �angular ��ω and ψψxare are �� phases. �problem iscelestial 90a degrees from the x axis in Earth ��ω ree body problem can be simplified circular body as system ofand�be ��ω where kare �angular � . αcalculated , α�be , Acalculated A� from can be from the initial condition, φ phases. ntric frame, the barycenter of ytheaxis twointo bodiesrestricted becomes the origin. The axis is and parallel to the line and are angular phases. and angular � �defined � �� and �,a, k where kthree � , the kψ � αMoon’s , A��ω A thecalculated initial condition, and φ and and angular the initial condition, and areand and phases. �ω ψ�ψ are angular phases. �the ��can where � ,�ω ��.phases. α��ω , �� α�� , and and Aangular canfrom bephases. calculated from the initial condition, and ψφ are ntric two celestial bodies becomes origin. The isk��ω �� toα� �parallel � angular �ωk��and and ψx angular phases. are phases. where k�are , to k��the �. line .Aand α� ,ψ αψare �ω ntric frame, frame, the the barycenter barycenter of of the the two of celestial bodies origin. The x�ωaxis axis parallel the line � �is�� , A� and A� can be calculated from the initial condition, and φ plane motion. Thebecomes y axis isthe positive in the general �ω�direction �ω�� angular phases. To analyze the periodic motion, assume andα zero. Then the linearized solution becomes[15] dies a circular orbit about barycenters, a third body negligible mass. Inperiodic the Earth-Moon � are n theinEarth and Moon, and is their directed from the and Earth toward the with Moon. The y axis isthe 90 degrees from the Toanalyze analyze the motion, assume andα arezero. zero.Then Thenthe thelinearized linearizedsolution solutionbecomes[15] becomes[15] To the periodic motion, assume αααα��α�andα ::: :: �are � are analyze the periodic motion, assume To analyze periodic motion, assume andα zero. Then the linearized becomes[15] To analyze the periodic assume zero. Then the linearized becomes[15] and are nn the from Earth Moon. The y is 90 degrees from the ToTo analyze the periodic motion, assume α�solution andα� are zero. To Then the linearized solution becomes[15 � andα �analyze ��are of the motion thetoward Moon relative to the Earth. Thephases. zmotion, axis and ψ the are angular phases. To analyze the motion, α andα are zero. Then the linearized solution : the periodic motion, assumeassume αsolution are Then linearized solution becomes[1 the Earth Earth and and Moon, Moon, and and is is directed directed from the the of Earth toward the Moon. The yψaxis axis isangular 90 periodic degrees from theassume � �the � andα and ψ are angular phases. To analyze the periodic motion, α�1becomes[15] andzero. α2 are zero. and ψ are angular phases. �� A c��ω � � φ� (13) ntric frame, the barycenter of the two celestial bodies becomes the origin. The x axis is parallel to the line � � c��ω �� in the Earth and Moon’s planeisofdefined motion. The y axis is positive in the generalcoordinate direction of the motion of to complete a right-handed system �� A �� φ� (13) �� A �� φ� (13) � c��ω �� A c��ω φ� (13) nn the is in the general direction of motion ofare To analyze the periodic motion, assume α��the andα are �� zero. Then linearized becomes[15] �� Athe φ� (13) Then linearized solution becomes[15] : � c��ω �� solution �� φ�φ�: �linearized �� analyze the periodic motion, assume zero. Then solution :A�Ac��ω �� A��c��ω φ� the becomes[15] (13) �� the Earth Earth and and Moon’s Moon’s plane plane of of motion. motion. The The y y axis axis To is positive positive the general direction of the theαassume motion � andα�of � c��ω �� Toin analyze periodic motion, andα Then�� the� � linearized solution becomes[15] : �� are Tothe analyze theby periodic motion,α�assume α� zero. andα and is normal to the plane of motion spanned the x and � are zero. Then the linearized solution becomes[15] : nonthe Earth to and from thetoEarth toward the Moon. The y axis issystem 90 degrees the �� k � φ�φ� (14) relative theMoon, Earth.and Theiszdirected axis is defined complete a right-handed coordinate and isfrom normal �A � ��ω �� k � �k A ��ω �� φ� (14) �� A ��ω on is to aa right-handed coordinate system and �A � ��ω �� A�� c��ω �� (13) �� k φ� �� φ� (14) (13)��(14) �� k ��ω �(14) �� φ�φ� ��Aφ� ��ω �� is A�� φ� (13) on relative relative to to the the Earth. Earth. The The zzyaxis axis is defined defined to complete complete right-handed coordinate system and is� �normal normal � A(13) axes[16]. From a combination of the kinematic expansion, �� k �k (14) �� k � c��ω � � �� A� c��ω �� φ� � A� ��ω�� φ� ��A� ��ω ��

�� A c��ω � � φ� (13) � �� in theofEarth andspanned Moon’s by plane motion. The y From axis isa positive theofgeneral direction of the motion of plane motion thea xofscalar, and y axes[16]. combination the kinematic expansion, a scalar, �� A� A ψ�ψ� (15) � c��ω �� second order and innonlinear set ofexpansion, differential �� A c��ω ��k �ψ� (15) �� (15) lane of motion spanned by the x and y axes[16]. From a combination of the kinematic a scalar, �c��ω � c��ω � �� A� ��ω�� φ� (14) �� A A � ψ� (15) c��ω � ψ� (15) � lane of motion spanned by the x and y axes[16]. From a combination of the kinematic expansion, a scalar, � � � A c��ω � � ψ� � � � �� k A ��ω � � φ� (14) � � �� A c��ω � � ψ� (15) �� A c��ω � � ψ� � � �k A�� ω � � φ� (14)(14) (14) � �� k�� aretoderived as:a right-handed � A� ��ω�� φ� on relative to the Earth. The zequations axis is defined complete coordinate system and is normal order and nonlinear set of differential equations are derived as: The different types ofofperiodic periodic orbits near collinear point are; Thedifferent differenttypes typesof periodicorbits orbitsnear near collinear point are; order The aaaaacollinear point are; The different(15) types of periodic orbits nea The different types of periodic orbits near point are; The orbits near collinear point are; �� collinear � � A� c��ω � ψ� order and and nonlinear nonlinear set set of of differential differential equations equations are are derived derived as: as: The different types of periodic orbits near a collinear point are; � � of Thedifferent differenttypes typesof ofperiodic periodic orbits near adifferent collinear point are; types periodic orbits near a collinear point are; �� A�� The �c��ω ψ� (15) � c��ω �A � � � � � ψ� (15)(15) (15) � � ��μ��μ� μ��μ� ofF�the kinematic expansion, a scalar, �� A c��ω � � ψ� plane of motion spanned by the x x�and y axes[16]. From a combination � � A � ��A � �: Planar Lyapunov Orbits (1) A(1) � �y� � x� � ��μ��μ� � μ��μ� � FF� � ��A � ��:�: A ��A Planar Lyapunov Orbits Planar Lyapunov ��μ��μ� μ��μ� � �� A� � ��A � �� x�x� � � A ��A Planar Lyapunov Orbits A(1) ��point ��A �:�: Lyapunov Orbits different typesnear of periodic near a�collinear point are; Orbits A�A� ��A� � �: Planar Lyapunov Orbits � � A�orbits ��A �:Planar Planar Lyapunov Orbits � �y� �y� � � x� x� � � �� The � different types of periodic a collinear are; � � ��The (1) �� a� � � � � ��A � � �: Planar Lyapunov Orbits The oforbits periodic near collinear point are; point � Thetypes �� different different types oforbits periodic orbits near a collinear are; order and nonlinear set of differential equations are derived as: The different types of periodic orbits near a collinear point A � �� : Vertical Lyapunov Orbits F� ��μ�� μ� � �� : A �� :Vertical Vertical Lyapunov Orbits A Lyapunov Orbits � �� A� � �� A �� : ��: Vertical Lyapunov Orbits y� � ��x� � y� � ��μ�� μ�� FF� (2) A(2) Vertical Lyapunov ��A � �: Planar Lyapunov Orbits AOrbits �� : Vertical Lyapunov Orbits �� Orbits A� � ��A �(2) �: Planar Lyapunov Orbits A�� : Vertical Lyapunov A� �� μ� � � � � y�y� � �� μ�� Lyapunov � � �� : Vertical Lyapunov Orbits A�� A �: Planar Orbits are; F� � ��x� ��x� � � y� y� � ��μ��μ� � �� μ��μ� �� (2) �� A � ��A � �: Planar Lyapunov Orbits � � � �� x� � �y� � x� � � � � � � A(1) ��ω ��ω ��ω ωω Lissajous Orbits( ororQuasi-Halo Quasi-Halo Orbits) � � ��A � �� :�Lissajous A ��A ��ω : Lissajous Lissajous Orbits(or Quasi-HaloOrbits) Orbits) A Orbits( � � � ��A � �≠0 �� A� Orbits) � ��A ��A ��ω Orbits( orQuasi-Halo Quasi-Halo Orbits) ��μ�� μ� � F �� : ��A ωω��:Lyapunov or AVertical &��ω ALyapunov =0 :�Planar Lyapunov Orbits A�Orbits � ��A ��ω �Orbits) ω : �Lissajous Orbits( oror Quasi-Halo �� : ::Lissajous �x �� or� �� ω zVertical A� � �� :AA Vertical Lyapunov Orbits A�� A Lissajous Orbits) AOrbits( � ��A � ��ω � : Lissajous Orbits( Quasi-Halo Orbits z� � ��μ�� μ�� FF�� (3) � � ��ω �� ω�Orbits �Orbits( � Quasi-Halo A�� : ��μ�� μ� �� z�z� � �� (3) A� � �� : Vertical Lyapunov Orbits F� μ� � �� μ�� (3) (3) A � ��A � ��ω � ω : Halo Orbits A =0 & A ≠0 : Vertical Lyapunov Orbits � � �� x z A � ��A � ��ω � ω : Halo Orbits � � � A ��A ��ω ω� ::��:Halo Orbits y� � ��x� � y� � � � � � � (2) �� A� � ��A Lissajous Orbits( or Quasi-Halo Orbits) ��A ��ω Halo Orbits AA �� Orbits �� A�or ��A ��ω �� ω : �Halo Orbits � ��A � ��A ��ω ω� :��Lissajous Orbits( Quasi-Halo Orbits) � � �� ω A�� A ��A ��ω �ωω ω� �or :Halo Halo A� � ��A � ��ω : Halo Orbits � � non-dimensional A�� ω ωA� : �� Lissajous Orbits( Orbits) � � �� & �Quasi-Halo d and r denote the non-dimensional distance from �the Earth to�the spacecraft andAthe Ax�� ≠0 ωz : LOrbits issajous Quasi-Halo A�� ω ��A ��ω ωxy Orbits(Orbits(or or� Quasi-Halo Orbits) z≠0 � �& �� ω � :≠Lissajous dd and r denote the non-dimensional distance from the Earth to the spacecraft and the non-dimensional and r denote the non-dimensional distance from to the spacecraft and theorbits non-dimensional The Lyapunov orbits are divided into horizontal Lyapunov orbit and vertical Lyapunov orbit. The horizontal ��μ�� μ� the F Earth TheLyapunov Lyapunov orbits aredivided divided into horizontal Lyapunov orbit andaaaaavertical verticalLyapunov Lyapunovorbit. orbit. The horizontal The are into aa ahorizontal Lyapunov orbit and The horizontal where d and r denote non-dimensional distance from the The Lyapunov orbits are dividedorbit. into aTh h A � ��A � ω�are : are Halo Orbits The Lyapunov orbits are divided into horizontal Lyapunov orbit and vertical Lyapunov orbit. The horizontal The Lyapunov orbits horizontal Lyapunov orbit and vertical Lyapunov orbit. horizontal Orbits) z� � respectively. �� μand � �the (3) The Lyapunov orbits divided a horizontal Lyapunov orbit and a vertical Lyapunov �into � � ��ω �� A� � ��A �divided ��ω ��ω ωaaaThe Halo Orbits Thenon-dimensional Lyapunov orbits are divided into Lyapunov orbit andinto ainto vertical Lyapunov orbit.The The horizontal orbits divided a horizontal Lyapunov orbit and a vertical Lyapunov orbit. T �� :horizontal Aare ��A ��ω ω Orbits e from the Moon to the spacecraft, mass of the � (1-μ) represent the �� Lyapunov � : Halo A � ��A � � ω : Halo Orbits ee from the Moon to the spacecraft, respectively. μand (1-μ) represent the non-dimensional mass of the � with no �out-of-plane �� motion � Earth respectively. to the spacecraft and the non-dimensional distance from the Moon to the spacecraft, μand (1-μ) represent the non-dimensional mass of the Lyapunov orbit is a planar orbit and lies entirely within the orbital plane of the Ano & Az≠0 &motion ωmotion :and Halo Lyapunov orbitis isa aplanar planar orbitwith withno no out-of-plane and liesaOrbits entirelywithin withinthe theorbital orbital plane ofthe theis a planar orbit with no x≠0 xy=ωzand Lyapunov orbit orbit out-of-plane lies entirely plane of Lyapunov orbit Lyapunov orbit is planar orbit with out-of-plane lies entirely within the orbital plane the Lyapunov orbit isnon-dimensional orbit with no out-of-plane motion and lies entirely within the plane ofof the The orbits are divided into aLyapunov horizontal Lyapunov orbit and vertical Lyapunov orbit. The horizontal Lyapunov orbit amotion planar orbit with nono out-of-plane motion and lies entirely Lyapunov orbit isaaaaplanar planar with no out-of-plane motion and lies entirely within theorbital orbital plane of the orbit is a planar orbit with out-of-plane motion and lies entirelywithin withinthetheorbital orbita TheEarth Lyapunov orbits areLyapunov divided into horizontal Lyapunov orbit and aisvertical Lyapunov orbit. The horizontal d and denote the non-dimensional distance from theLyapunov spacecraft and the and ther non-dimensional massfrom of thethe Earth, respectively.F , F�to andF denote external perturbation forces Moon to the the spacecraft, respectively. μand (1-μ) The orbits are divided into aorbit horizontal Lyapunov orbit and aorbit vertical Lyapunov orbit. The horizontal � �denote The Lyapunov orbits are divided into a horizontal and the non-dimensional mass of the Earth, respectively.F , F andF external perturbation forces The Lyapunov orbits are divided into a horizontal Lyapunov and a vertical Lyapunov orbit. The horizontal � � � primaries. The vertical Lyapunov orbit is almost a vertical orbit with a typical eight-shape doubly symmetric and the non-dimensional mass of the Earth, respectively.F� , F� andF� denote external perturbation forces primaries. The vertical Lyapunov orbit is almost a vertical orbit with a typical eight-shape doubly symmetric primaries. The vertical Lyapunov orbit is almost a vertical orbit with a typical eight-shape doubly symmetric primaries. The vertical Lyapunov orbit i represent the non-dimensional mass of the Moon and the primaries. The vertical Lyapunov orbit is almost a vertical orbit with a typical eight-shape doubly symmetric primaries. The vertical Lyapunov orbit is almost a vertical orbit with a typical eight-shape doubly symmetric orbit is ano planar with nois out-of-plane motion andwith lies within orbital plane ofa the primaries. The orbit is isalmost aofthe vertical orbit with typical Lyapunov orbit and aLyapunov vertical Lyapunov orbit. horizontal primaries. The vertical Lyapunov orbit almost a vertical orbit aentirely typical eight-shape doubly symmetric primaries. The vertical Lyapunov orbit almost aThe vertical with a typicaleight-shape eight-shapedoubly doub Lyapunov orbit is the aLyapunov planar with out-of-plane motion and lies entirely within the orbital the fromy axes the Moon to thrust the spacecraft, respectively. μand (1-μ) represent non-dimensional mass oforbit theout-of-plane Lyapunov orbit is orbit a orbit planar orbit with no motion and lies entirely within theplane orbital plane of orbit the xe and and the of the spacecraft respectively. Lyapunov is,axisymmetric a and planar orbit withthe nox-axis out-of-plane motion and lies entirely within the orbital plane of the xx and y axes and the thrust of the spacecraft respectively. geometry which is about and symmetric with respect to the xz plane. non-dimensional mass of the Earth, respectively. F F F and y axes and the thrust of the spacecraft respectively. xis y z geometrywhich which axisymmetric about the x-axis and symmetric with respect tono the xzplane. plane. doubly geometry is axisymmetric about the and with respect to the xz Lyapunov orbit is a planar orbit out-of-plane motion geometry which is to axisymmetric about th geometry which axisymmetric about the x-axis and symmetric with respect to the xz plane. geometry which isis axisymmetric about the x-axis and symmetric with respect tothe the xz plane. primaries. Theperturbation vertical Lyapunov orbit isx-axis almost a symmetric vertical orbit with awith typical eight-shape symmetric geometry which axisymmetric about the and symmetric with respect thethe xzxz plane. axisymmetric about the x-axis and symmetric with respect tox-axis the xz plane. geometry which is axisymmetric about x-axis and symmetric with respect to plane. The vertical Lyapunov orbit is and almost vertical with ais typical doubly symmetric and of the Earth, respectively.F , F�primaries. andF external forces The vertical Lyapunov orbit isa almost aorbit vertical orbit with aeight-shape typical eight-shape doubly symmetric � denote e ofthethenon-dimensional three collinearmass liberation points (x� �perturbation in primaries. the �planar case, itgeometry isthe assumed that the external denote external forces on x which and yisaxes primaries. The vertical Lyapunov orbit is almost a vertical orbit with a typical eight-shape doubly symmetric ee of the three collinear liberation points (x � in the planar case, it is assumed that the external and lies with entirely withinand theout-of-plane orbital plane of thebutprimaries. The Lissajous orbit is defined as a trajectory an in-plane oscillation, the oscillation of the three collinear liberation points (x�� in the planar case, The itThe isLissajous assumed that the external Lissajous orbit isdefined definedas asaabout trajectory with ansymmetric in-planeand and out-of-plane oscillation, but theLissajous oscillation orbit is trajectory with an in-plane out-of-plane oscillation, but the oscillation The orbit is defined as abut traje The Lissajous orbit is defined as trajectory with an in-plane and out-of-plane oscillation, but the oscillation Lissajous orbit isthe as aaaasymmetric trajectory with an in-plane and out-of-plane oscillation, but the oscillation geometry which is axisymmetric the x-axis and with respect to the xz plane. The Lissajous orbit is to as a trajectory ananin-plane thrust ofcan thebe spacecraft respectively. The Lissajous orbit isdefined defined asx-axis trajectory with an in-plane and oscillation, but theand oscillation The Lissajous orbit isdefined defined as trajectory with in-plane andout-of-plane out-of-planeoscillation, oscillation, butth geometry which isThe axisymmetric about x-axisthe and with respect the xz plane. xation and and y axes andare thezero. thrustThen, of the theEqs.(1)~(3) spacecraft respectively. The vertical Lyapunov orbit is awith vertical orbit with which is axisymmetric about and symmetric with respect toout-of-plane theaalmost xz plane. thrust written as: geometry geometry which is axisymmetric about the x-axis and symmetric with respect to the xz plane. ation and thrust are zero. Then, Eqs.(1)~(3) can be written as: frequencies ofofthe the in-plane and out-of-plane motion differ over time. The Lissajous orbits can bebeseen seen asasquasiquasiation and thrust are zero. Then, Eqs.(1)~(3) bethree writtencollinear as: On one ofcan the liberation points (xin-plane ) isindefined the frequencies the andout-of-plane out-of-plane motion differ over time. The Lissajous orbitscan canbe seenoscillation quasifrequencies of and motion differ over time. The Lissajous orbits as ein-plane frequencies of the in-plane and can out-of-pl frequencies the in-plane and out-of-plane motion differ over time. The Lissajous orbits can be seen quasitypical eight-shape doubly symmetric geometry frequencies ofof the in-plane and out-of-plane motion differ over time. The Lissajous orbits can be seen asas quasiThe Lissajous orbit as trajectory with an in-plane and out-of-plane oscillation, but the ofout-of-plane the in-plane and out-of-plane differ over time. The Lissajous orbits bebe se the in-plane and out-of-plane differ over time. The Lissajous orbits can bewhich seen asis quasifrequencies of the in-plane and out-of-plane motion over time. The Lissajous orbits can orbit isisdefined asdefined athat trajectory with anaafrequencies in-plane and oscillation, but themotion oscillation e of the three collinear liberation points (x� in The the Lissajous planar itfrequencies assumed the Thecase, Lissajous orbit isof asdefined aexternal trajectory with an motion in-plane and out-of-plane oscillation, but thediffer oscillation �� x� � �y� � �� x (4)as � �assumed The Lissajous orbit is a trajectory with an in-plane and out-of-plane oscillation, but the oscillation planar case, it is that the external perturbation and x�x� � �y� � � ��x (4) symmetric to the xy-plane and xz-plane. axisymmetric about the x-axis and symmetric with respect � �y� � �� x (4) symmetricto tothe the xy-plane and xz-plane. symmetric xy-plane and xz-plane. symmetric to the xy-plane and xz-plane. symmetric the xy-plane and xz-plane. symmetric toto the xy-plane and xz-plane. frequencies the in-plane and out-of-plane motion differ overand time. The orbits can be seen as quasisymmetric to to the xy-plane xz-plane. symmetric toof the xy-plane and xz-plane. the xy-plane and xz-plane. frequencies of the in-plane out-of-plane motion differ over time. The Lissajous orbits canLissajous be seen thrust Then, Eqs.(1)~(3) can be as: ation and thrust are zero. Then, Eqs.(1)~(3) be��y written as: frequencies of written the and in-plane out-of-plane motion differ over time. The time. Lissajous orbits canas bequasiseen asbequasi�� y� � are ��x� zero. �can (5) to symmetric the xz plane. frequencies of theand in-plane and out-of-plane motion differ over The Lissajous orbits can seen as quasi9 �� y�y� � ��x� � � ��y (5) 99 � ��x� � �� y (5) 9 9 with an symmetric to the xy-plane and xz-plane. The Lissajous999orbit is defined as a trajectory to the xy-plane and xz-plane. (4) to the xy-plane and xz-plane. �y� � �� x symmetricsymmetric (4) z�x� � ��z (6) symmetric to the xy-plane and(6) xz-plane. z�z� � ��z � ��z (6) 9 9 9 ��μ� μ y� � ��x� � �� y (5) 9 μ � � |��μ� � . � � |� ��μ| � . ��μ� μ �μ| �� . |� � |� �μ|� � � |�� μ| � |�� μ| ��μ| � � z� � ��z (6) DOI:10.5139/IJASS.2014.15.3.267 270 genvalue of Eqs.(4)~(6) is envalue of Eqs.(4)~(6) is envalue of Eqs.(4)~(6) is ��μ� μ � � |� � � |� � . � �μ| � ��μ| 8 8 8 genvalue of Eqs.(4)~(6) is

8


Moon L2 quasi-halo orbit. There are three types of transfers to enter the quasi-halo orbit; direct transfer, low thrust transfer and ballistic lunar trajectory transfer. Parker[20] calculated the fuel savings for these different types of transfers for a spacecraft travelling from the Earth parking orbit at 185km altitude to the Earth-Moon L2 point. Parker showed that the direct transfer from LEO to L2 requires a total delta-V of 3.9823 km/s over 5 days of transfer time. This paper adopted the direct transfer because the low thrust and ballistic lunar trajectory transfers to the Moon require relatively long transfer times. Such long durations to reach the quasi-halo orbit would decrease the reliability of the mission[21]. Although halo orbit is both periodic and time independent in the CR3BP, it does not consider all the gravitational and non-gravitational effects of the full solar system, the solar radiation pressure and the non-uniform gravity of the Earth and the Moon. A differential correction method is used to compute maneuver time and delta-V to inject a satellite into an Earth-Moon L2 quasi-halo orbit departing from a circular low Earth parking orbit. A quasi-halo orbit insertion maneuver is also investigated using the differential corrector. The STK/Astrogator contains full-force models and interplanetary propagation/targeting techniques. Table 2 lists the initial state of a spacecraft on the low Earth orbit. The mission starts from an initial low Earth circular parking orbit followed by a quasi-halo orbit injection maneuver near the Earth-Moon L2 Lagrange point. The insertion point of the quasi-halo orbit is selected at 20,000 km in the z axis, and 0 km in the x and y axes with the origin of the EarthMoon L2 point. To make a periodic orbit for station keeping, the end condition for this propagation segment after one period requires that the spacecraft cross the Earth-Moon L2 Z-X plane with Vx=0. The delta-v required for quasi-halo orbit insertion using the differential-corrector method in STK Astrogator is obtained as shown in Table 3. Table 3 lists the maneuver information and state of the spacecraft on transfer orbit for insertion into the Earth-Moon L2 quasihalo orbit[22]. After the insertion maneuver into the quasi-halo orbit, it is assumed that collision is expected from an uncontrolled object on the Earth-Moon L2 quasi-halo orbit. Table 4 shows the initial state of the uncontrolled object. If a potential collision risk is detected, a precise orbit should be determined by using all possible tracking data from global observation sites on Earth. Table 5 lists orbit determination covariance and the collision information of the assumed spacecraft with the uncontrolled object. Once the collision risk is greater than a maximum threshold, a plan for avoidance maneuver should be generated. An

in-plane and out-of-plane oscillation, but the oscillation frequencies of the in-plane and out-of-plane motion differ over time. The Lissajous orbits can be seen as quasisymmetric to the xy-plane and xz-plane. The halo orbit has an in-plane and out-of-plane motion. However, the frequency stays the same and the orbit is periodic. The symmetricity to the xy-plane of the primaries no longer exists. The halo orbit only exists for the minimum in-plane amplitude and provides an exclusion zone located around the line connecting the primaries. The quasi-halo orbit is a “mixture” of the Lissajous and halo orbits. The orbit circulates through the same path as the halo orbit, but its shape is similar to torus. The quasihalo orbit originates from the Lissajous orbits with a specific minimum boundary value of the out-of-plane amplitude. At this boundary amplitude, the Lissajous orbit loses their symmetry with respect to the xy-plane and starts to develop an exclusion zone about the line connecting the primaries. It is very useful for science and observation missions with formation flights[3, 17, 18].

4. Insertion Maneuver Planning for EarthMoon L2 Quasi-Halo Orbit In CR3BP (Circular Restricted Three Body Problem), the third small body is located at the Lagrange point where the gravitational pull of two large bodies is balanced. There exist five Lagrange points in the Earth-Moon rotating coordinate system. L1, L2 and L3 are saddle points which are stable along the axis between the two primary bodies but unstable in the other direction. Station keeping maneuvers are therefore required to remain at these points. The L4 and L5 points are stable and do not require station keeping maneuvers. The Quasi-halo orbit can perform periodic motion but is not actually an orbit around the Lagrange point. The Quasi-halo orbit can also supply trajectories to transfer into and out of the Lagrange orbits by providing low energy manifolds[18, 19]. A spacecraft placed in an Earth-Moon L2 quasi-halo orbit can communicate between the Earth and the far side of the Moon. This quasi-halo orbit could be used to establish a lunar space station that can serve as a gateway for solar system exploration or as a fueling spacecraft to help other spacecrafts depart the Earth Moon system with lower maneuver costs. This can reduce the overall size of the spacecraft as well as the cost of the mission. For this mission, a spacecraft requires an impulsive burn to transfer from the Earth parking orbit to the point of insertion into the quasi-halo orbit. Additionally, an insertion maneuver is required to approach an Earth-

271

http://ijass.org

β � ∑�� erf �

�� ·R�� ·��·� σ� √�

(18)

�

R �� denotes the total cross sectional radius of two objects assumed to be circular. σ� , σ� , x� and y� denote the standard deviation and distance from the center in the direction of each axis.

The control states are impulsive burn time and delta-v in the radial, in-track and cross-track directions. Int’l J. of Aeronautical & Space Sci. 15(3), 267–280 (2014) The fitness value can be written as

acceptable avoidance maneuver plan should produce an orbit (19) (16) � � ∑N �� ∆L� �t � � � ��P� change that reduces the risk of collision while still meetingth where ∆L represents the i position difference between the orbit difference and an orbit generated thereference ith position between theby the where ∆Li represents positioning constraints to retain the quasi-halo� orbit. reference and an orbit generated by and the aoptimization Figure 2. shows a flowchart to optimization calculate the optimal algorithm. P� and K representorbit a maximum collision probability weighting factor, algorithm. Pc and K represent a maximum collision probability maneuver and evaluate its performance for collision respectively. The maximum collision probability is set to 10�� and collision miss distance is set to 100km. and a weighting factor, respectively. The maximum collision avoidance while remaining in the quasi-halo orbit. STK/Pro probability is set to 10-4 and collision miss distance is set to to is used as the starting point to model, analyze, visualize Global Opt. Tman 100km. Algorithm and communicate with other components. The STK/ STK/Pro VRIC (Matlab) Astrogator propagates the orbit of a spacecraft with delta-V No

Yes at the maneuver start time generated from the Matlab global Converge? End 5. Simulation Results optimization algorithm. It then integrates an accurate full perturbation model that includes solar radiation pressure Calculate Fitness Value For comparison, ten simulation trials were performed to and planet gravitational perturbation, where “Tman” and determine the influence of the cost function on each global “ΔVRIC” represent the maneuver time in UT, the delta-V in STK/Astrogator STK/CAT radial, in-track and cross-track respectively. The collision Table 2. The Initial State of the Spacecraft on the Earth Orbit for Injecprobability is obtained from STK/AdvCAT, which adopts the Fig.2. A Flowchart of Optimal Maneuver Planning to TransferonOrbit Alfano formula[23, 24, 25, 26, 27]. Alfano developed the twoTable 2.The Initial State oftion the Spacecraft the Earth Orbit for Injection to Transfer Orbit dimensional collision probability equation in the encounter Item Value plane given in Cartesian space as a combination of error Coordinate System Earth Inertial functions and exponential terms as shown in Eqs.(16)~(18). Coordinate Type Target Vector Outgoing Asymptote Alfano used a rectangular coordinate system� to obtain the R�� ·��·�� Orbit Epoch 1 Jul 2025 12:00:00.000UTCG �� R�� ·� collision probability value, but�� did� not��assume that the · �α � β� · ex� � � (16) P� � � Radius of Periapsis 6778 km �·σ� √�·π·σ� ·� velocity remains constant for each tube segment.

�

R�� ·��·�� ·R�� ·��·� ��R�� ·��·�� ·�� erf � αRR� � � �� ∑ ��R�� ·��·�� ·� �� σ � · �α � β� · ex� � P √� R�� ·� · �α � β� · ex�� √�·π·σ ��·σ� P ·� �� P�� √�·π·σ��·� · �α � β� · ex� � �·σ�� ·σ� √�·π·σ� ·� � �� ·R�� ·��·� β � ∑�� erf � � σ� √� �� ·R�� ·��·� � � ��·R ��·��·� � ��·R ��·��·� ∑ α � �� erf �� ∑�� erf �∑ erf �� σ √� αα � �� σ �√� �� σ �√�

(16)

C3 Energy

-1.62335 km^2/Sec^2

RA of Outgoing Asymptote (17)

252.107 deg

(16) (16) (16) Declination of Outgoing

(17)

R �� denotes the total cross sectional radius of� two objects assumed to be circular. σ� , σ� , x� and y� denote �� ·R�� ·��·�

�� ·R�� ·��·� � ∑��erf β� �∑ erf �� ·��·� � � ��·R�� σ √� in the direction the standard deviation and distance from of each axis. �� erf ∑�� the center � ββ � σ �√� σ��√�

-24.848 deg

Asymptote Velocity Azimuth at (17) (17) Periapsis (17)

(18)

(18) (18) (18)

(18)

82.9741 deg

True Anomaly13

control states are impulsive time and objects delta-v assumed in the radial, andσσ�cross-track directions. denotes the total total cross cross sectional burn radius of two two objects assumed to be bein-track circular. , σ� , xx� and and denote �� The � denote RR�� denotes the sectional radius of to circular. yy� � , σ�,, x � R radius of objects assumed to beradius circular. �� denotes the total cross sectional � , σ�objects � and y� denote Robj denotes thetwo total cross sectional ofσtwo The fitness value and can distance be written as the center in the direction of each axis. the standard standard deviation from assumed to the be center circular. σx,direction σy, xm and ym denote the deviation and distance distance from from the center in in the the direction of each each axis. the standard Spacecraft Properties the standard deviation and of axis. N from the center in the direction of each deviation and distance ∑�� ∆L �t � radial, � �delta-v (16) � ��P�in-track and cross-track directions. The control states are impulsive burn time and in the � � The control control states states are are impulsive impulsive burn burn time and delta-v in the radial, in-track and cross-track directions. The axis. th time and delta-v in the radial, in-track and cross-track directions. wherevalue ∆L� represents the as i position difference between the reference orbit and an orbit generated by the The fitness fitness value can be be written written as The can The The fitness value can be written as control states are impulsive burn time and delta-v in

0.000786134 deg Total Mass: 800kg(Dry Mass: 500 kg, Fuel Mass : 300 kg) Area /Mass Ratio: 1.25e-9 km^2/kg Tank Pressure : 5000 Pa, Fuel Density: 1000 kg/m^3 Cr: 2.0, Cd: 2.0, Drag Area:1m^2, SPR Area: 1m^2 Solar Radiation Pressure Area : 20m^2 Radiation Pressure Coefficient : 1.0

N the radial, and directions. optimization algorithm. P� and�in-track K collision probability and a weighting ∑represent �a �maximum �∑ (16)factor, ��P� N � �t � cross-track �� ∆L � (16) �� ∆L� �t � � � ��P� ∑N �� (16) �� ∆L� �t � � � ��P� The fitness value can be written as �� The maximum collision probability is set to and collision miss is set to 100km. whererespectively. ∆L� represents represents the iiththth position position difference between the10 reference orbit and and an distance orbit generated generated by the the where ∆L the difference between the reference orbit an orbit by Table 3. ITable nsertion Maneuver by Differential CorrectorCorrector 3. Insertion Maneuver by Differential � where ∆L � represents the i position difference between the reference orbit and an orbit generated by the

optimization algorithm. algorithm. PP� and and K K represent represent aa maximum maximum collision probability and and a weighting weighting factor, factor, Global Opt. probability optimization collision � optimization algorithm. P and aa weighting factor, � and K represent a maximum collision probability Tman

Algorithm �� respectively. The The maximum maximum collision collision probability probability is is set set to to 10 10�� distance is set set to to 100km. 100km. RIC �� and collision missV respectively. and collision collision miss miss distance distance is (Matlab) respectively. The maximum collision probability is set to 10 and is set to 100km. STK/Pro No

Converge?

STK/Pro STK/Pro STK/Pro

Yes Global Opt. End Global Opt. Global Opt. Algorithm Algorithm Algorithm (Matlab) (Matlab) (Matlab)

Calculate No No Fitness Value No Yes Yes Converge? Yes Converge? Converge?

STK/CAT

Calculate Calculate Calculate Fitness Value Fitness Value Fitness Value

man TTman man T V RIC V RIC VRIC

End End End

STK/Astrogator

Maneuver Type

Impulsive

Engine Model

Constant Trust and Isp(500N, 300 sec)

CoordinateSystem

Earth Inertial

Maneuver Time

7 Jul 2025 23:20:41.629 UTCG

Delta-V(km/sec)

Vx:0.88709, Vy: -0.29449, Vz: -0.12307

Estimated Fuel Used

219.342 kg 1290.61 sec

Duration

STK/Astrogator STK/Astrogator STK/Astrogator

Fig.2. A A Flowchart Flowchart of of Optimal Optimal Maneuver Maneuver Planning Planning Fig.2. DOI:10.5139/IJASS.2014.15.3.267 Fig.2. A Flowchart of Optimal Maneuver Planning

Value

Estimated Burn

Fig.2. A Flowchart ofManeuver Optimal Maneuver Fig. 2. A Flowchart of Optimal Planning Planning

STK/CAT STK/CAT STK/CAT

Item

Table 4.The Initial State of an Uncontrolled Object

272

Central Body

Moon

Orbit Epoch

15 Jul 2025 17:19:14.361 UTCG

Coordinate Type

Cartesian

Coordinate System

J2000

14

CoordinateSystem CoordinateSystem

Earth Earth Inertial Inertial

Maneuver Maneuver Time Time

7 7 Jul Jul 2025 2025 23:20:41.629 23:20:41.629 UTCG UTCG

Delta-V(km/sec) Delta-V(km/sec)

Vx:0.88709, Vx:0.88709, Vy: Vy: -0.29449, -0.29449, Vz: Vz: -0.12307 -0.12307

Estimated Estimated Fuel Fuel Used Used

219.342 219.342 kg kg

Estimated Estimated Burn Burn 5. Simulation Results 1290.61 1290.61 sec sec Duration Sang-Cherl Lee Comparison of Global Optimization Methods for Insertion Maneuver into Earth-Moon L2 ... Duration For comparison, ten simulation trials were performed to determine the influence of the cost function on each global optimization case. In this simulation, the optimization parameters ‘MaxIter’, ‘MaxFunEvals’ and ‘TolFun’ Table 4.The Initial of Table 4. The Initial State of anState Uncontrolled Object Object Table 4.The Initial State of an an Uncontrolled Uncontrolled Object parameters ‘MaxIter’, ‘MaxFunEvals’ and ‘TolFun’ is set to is set to 4000, 4000 and 10e-6, respectively. If the optimality measure less 10e-6, than ‘TolFun’, the iterations 4000, 4000isand respectively. If theend, optimality measure Central Moon Central Body Body Moon

is less than ‘TolFun’, the iterations end, and ‘TolFun’ can also be a relative bound. ‘MaxIter’ is a bound on the number of is a bound on the number of function evaluations. Coordinate Type Cartesian solver iterations. ‘MaxFunEvals’ is a bound on the number of Coordinate ‘MaxFunEvals’ Type Cartesian evaluations. Fig.3.a~c show the 3D view for one of 10 cases in view of function the Moon fixed frame to get an insertion maneuver Coordinate J2000 Coordinate System System J2000 Figure 3.a~c show the 3D view for one of 10 cases in Propagator HPOP into a quasi-halo orbit by a simulated annealing, a genetic algorithm and a pattern search algorithm. The Propagator HPOP view of the Moon fixed frame to get an insertion maneuver Position(km) X: 50665.7, Y: 3519.95, Z: -14372.3 Position(km) X: object 50665.7, Y: 3519.95, Z: the -14372.3 into toward a quasi-halo orbit on by the a simulated annealing, a genetic uncontrolled approaches from direction of the Moon the spacecraft quasi-halo orbit. Velocity(km/sec) Velocity(km/sec) Vx: Vx: 0.028431,Vy: 0.028431,Vy: 0.34545, 0.34545, Vz: Vz: 0.75294 0.75294 algorithm and a pattern search algorithm. The uncontrolled Possible collision is declared by the conjunction analysis software and the global optimal algorithm attempts to object approaches from the direction of the Moon toward find the best solution with the boundary constraints provided. fitness values at eachorbit. iteration. theThe spacecraft on are thechecked quasi-halo Possible collision Table 5.Orbit Determination Covariance and Collision Information Table 5. O rbit Determination Table 5.Orbit DeterminationCovariance Covarianceand andCollision CollisionInformation Information is declared by the conjunction analysis software and the The white line indicates the Moon trajectory and the green line stands for the spacecraft trajectory before the T(km) C(km) N(km) Size(m) T(km) C(km) N(km) Size(m) global optimal algorithm attempts to find the best solution L2 55 55 red line represents the quasi-halo orbit after the insertion maneuver orbit 55insertion. The L2 Probe Probe maneuver55 for quasi-halo with the boundary constraints provided. The fitness values HALO 55 55 55 50 HALO Object Object 50 are checked at each iteration. The white converged using the global optimization algorithm. The quasi-halo orbit is periodic in the Y-Z axis as shown in line indicates Collision Collision Exp. Exp. 15 Jul 2025 17: 19: 14 the Moon trajectory and the green line stands for the 15 Jul 2025 17: 19: 14 Time Time the left side images of Figs. 3. The right side images of Figs. 3. in the X-Z axis illustrate that inclinations are Minimum spacecraft trajectory before the maneuver for quasi-halo Minimum 0.534 0.534 spacecraft trajectory Dist.(km) orbit difference insertion.between The redtheline represents the quasi-halo orbit Dist.(km)formed. Figs. 4 show radial, in-track and cross-track distance after the insertion maneuver converged Collision Prob. 1.02e-2 Collision Prob. 1.02e-2corrector and the spacecraft trajectory generated by global optimization methodsusing the global generated by the differential optimization algorithm. The quasi-halo orbit is periodic for the first case 1. in the Y-Z axis as shown in the left side images of Figs. 3. optimization case. In this simulation, the optimization and ‘TolFun’15can also be a relativeUTCG bound. ‘MaxIter’ is a bound on the number of solver iterations. Orbit Orbit Epoch Epoch 15 Jul Jul 2025 2025 17:19:14.361 17:19:14.361 UTCG

HALO Transfer Object Orbit

Orbit

Object Orbit

15 15 HALO Transfer Orbit Moon Orbit

EarthMoonQuasi-Halo

Moon Orbit

EarthMoonQuasi-Halo

Fig.3a Simulation Results (Simulated Annealing Case 1) Fig. 3a. Simulation Results (Simulated Annealing Case 1) 16

Fig.3b Simulation Results (Genetic Algorithm Case 1) Fig. 3b. Simulation Results (Genetic Algorithm Case 1)

273

http://ijass.org


spacecraft trajectory generated by global optimization The right side images of Figs. 3. in the X-Z axis illustrate methods for the first case 1. that inclinations are formed. Figs. 4 show radial, in-track Fig.3b Simulation Results Algorithm Tables 6a~cCase and1) and cross-track distance difference between spacecraft Fig.3bthe Simulation Results (Genetic (Genetic Algorithm Case 1)Fig. 5a~c show the simulation results for avoidance maneuver planning using global optimization. trajectory generated by the differential corrector and the

Fig.3c Fig.3c Simulation Simulation Results Results (Pattern (Pattern Search Search Case Case 1) 1) Fig. 3c. Simulation Results (Pattern Search Case 1)

Fig.4a RIC Distance Difference (Simulated Fig. 4a. RIC and Range Distance Difference Annealing Case 1) Fig.4a(Simulated RIC and and Range Range Distance Difference (Simulated Annealing Annealing Case Case 1) 1)

17 17

Fig.4b RIC and Range Distance (Genetic Algorithm Case 1) Fig. 4b. RIC and Range Distance Difference (Genetic Algorithm Case 1) Difference Fig.4b RIC and Range Distance Difference (Genetic Algorithm Case 1)

Fig.4c RIC and Range Distance Difference (Pattern Search Case 1) Fig.4c RICSearch and Range Difference (Pattern Search Case 1) Fig. 4c. RIC and Range Distance Difference (Pattern CaseDistance 1) Tables 6a~c and Fig. 5a~c Tables 6a~c and Fig. 5a~c DOI:10.5139/IJASS.2014.15.3.267 optimization. The maximum optimization. The maximum

show the simulation results for avoidance maneuver planning using global show the simulation results for avoidance maneuver planning using global

274 iteration for each optimization method was set to 4,000. For the simulated iteration for each optimization method was set to 4,000. For the simulated

annealing, the minimum best-cost-function was 1.6704e4, while the maximum best-cost-function was 4.8500e4. annealing, the minimum best-cost-function was 1.6704e4, while the maximum best-cost-function was 4.8500e4. The minimum best-cost-function was 1.353e4 using the genetic algorithm, and the maximum was 2.902e4. With The minimum best-cost-function was 1.353e4 using the genetic algorithm, and the maximum was 2.902e4. With the pattern search method, the minimum best-cost-function was 1.53249e4, and the maximum was 2.8472e4. All the pattern search method, the minimum best-cost-function was 1.53249e4, and the maximum was 2.8472e4. All

but at 1000 but ended ended at the the maximum maximum 1000 iterations. iterations. Sang-Cherl Lee Comparison of Global Optimization Methods for Insertion Maneuver into Earth-Moon L2 ...

Table 6a. Simulation Result after Avoidance Maneuver Using Simulated Annealing Table 6a. Result after Maneuver Table 6a. Simulation Simulation ResultPlanning after Avoidance Avoidance Maneuver Planning Planning Using Using Simulated Simulated Annealing Annealing Case Case

Man_Time Man_Time

DelV_x DelV_x

DelV_y DelV_y

DelV_z DelV_z

No No

(UTCG) (UTCG)

(km/s,ECI) (km/s,ECI)



11

77 Jul Jul 2025 2025 23:10:01.920 23:10:01.920

0.850505 0.850505

-0.291555 -0.291555

22

77 Jul Jul 2025 2025 23:10:01.920 23:10:01.920

0.901398 0.901398

33

77 Jul Jul 2025 2025 23:10:01.920 23:10:01.920

44 55

TotalDelV TotalDelV

Cost Cost Fun Fun

Max. Max. Func. Func. Eval. Eval.

Collision Collision

-0.130624 -0.130624

0.908530 0.908530

3.2527e4 3.2527e4

4000 4000

No No Warn. Warn.

-0.289530 -0.289530

-0.131955 -0.131955

0.955907 0.955907

1.7037e4 1.7037e4

4000 4000

No No Warn. Warn.

0.889694 0.889694

-0.303757 -0.303757

-0.097427 -0.097427

0.945155 0.945155

1.6704e4 1.6704e4

4000 4000

No No Warn. Warn.

77 Jul Jul 2025 2025 23:10:01.920 23:10:01.920

0.864243 0.864243

-0.291628 -0.291628

-0.130641 -0.130641

0.921428 0.921428

2.4393e4 2.4393e4

4000 4000

No No Warn. Warn.

77 Jul Jul 2025 2025 23:28:19.200 23:28:19.200

0.865691 0.865691

-0.296707 -0.296707

-0.119394 -0.119394

0.922882 0.922882

2.3530e4 2.3530e4

4000 4000

No No Warn. Warn.

66

77 Jul Jul 2025 2025 23:29:36.960 23:29:36.960

0.838590 0.838590

-0.287033 -0.287033

-0.140430 -0.140430

0.897409 0.897409

3.7957e4 3.7957e4

4000 4000

No No Warn. Warn.

77

77 Jul Jul 2025 2025 23:19:06.240 23:19:06.240

0.837496 0.837496

-0.297909 -0.297909

-0.114908 -0.114908

0.896300 0.896300

4.2240e4 4.2240e4

4000 4000

No No Warn. Warn.

88

77 Jul Jul 2025 2025 23:10:27.840 23:10:27.840

0.921033 0.921033

-0.301070 -0.301070

-0.097284 -0.097284

0.973863 0.973863

3.4012e4 3.4012e4

4000 4000

No No Warn. Warn.

99

77 Jul Jul 2025 2025 23:23:25.440 23:23:25.440

0.826213 0.826213

-0.279156 -0.279156

-0.154978 -0.154978

0.885762 0.885762

4.8500e4 4.8500e4

4000 4000

No No Warn. Warn.

10 10

77 Jul Jul 2025 2025 23:20:24.000 23:20:24.000

0.869968 0.869968

-0.267055 -0.267055

-0.179721 -0.179721

0.927611 0.927611

3.4286e4 3.4286e4

4000 4000

No No Warn. Warn.

Prob. Prob.

Table 6b. Simulation Result after Avoidance Maneuver Planning Using Genetic Algorithm Table 6b. Result after Maneuver Planning Table 6b. Simulation Simulation Result after Avoidance Avoidance Maneuver Planning Using Using Genetic Genetic Algorithm Algorithm Total Total DelV DelV

Best Best Cost Cost Fun Fun

Max. Max. Fun. Fun. Eval. Eval.

Collision Collision

-0.094120 -0.094120

0.952380 0.952380

2.902e4 2.902e4

4000 4000

No No Warn. Warn.

-0.287613 -0.287613

-0.139542 -0.139542

0.940179 0.940179

1.353e4 1.353e4

4000 4000

No No Warn. Warn.

0.879809 0.879809

-0.285391 -0.285391

-0.144266 -0.144266

0.936123 0.936123

1.621e4 1.621e4

4000 4000

No No Warn. Warn.

77 Jul Jul 2025 2025 23:10:27.840 23:10:27.840

0.909170 0.909170

-0.296791 -0.296791

-0.112597 -0.112597

0.962992 0.962992

2.455e4 2.455e4

4000 4000

No No Warn. Warn.

77 Jul Jul 2025 2025 23:08:26.880 23:08:26.880

0.920401 0.920401

-0.293585 -0.293585

-0.117621 -0.117621

0.973224 0.973224

2.705e4 2.705e4

4000 4000

No No Warn. Warn.

66

77 Jul Jul 2025 2025 23:25:09.120 23:25:09.120

0.888628 0.888628

-0.294739 -0.294739

-0.122487 -0.122487

0.944211 0.944211

1.373e4 1.373e4

4000 4000

No No Warn. Warn.

77

77 Jul Jul 2025 2025 23:25:00.480 23:25:00.480

0.879775 0.879775

-0.296022 -0.296022

-0.120297 -0.120297

0.936005 0.936005

1.855e4 1.855e4

4000 4000

No No Warn. Warn.

88

77 Jul Jul 2025 2025 23:28:36.480 23:28:36.480

0.869536 0.869536

-0.293840 -0.293840

-0.126260 -0.126260

0.926486 0.926486

2.019e4 2.019e4

4000 4000

No No Warn. Warn.

99

77 Jul Jul 2025 2025 23:13:55.200 23:13:55.200

0.917336 0.917336

-0.285886 -0.285886

-0.137184 -0.137184

0.970596 0.970596

2.292e4 2.292e4

4000 4000

No No Warn. Warn.

10 10

77 Jul Jul 2025 2025 23:17:48.480 23:17:48.480

0.909015 0.909015

-0.288382 -0.288382

-0.133660 -0.133660

0.962984 0.962984

1.891e4 1.891e4

4000 4000

No No Warn. Warn.

Case Case

Man_Time Man_Time

DelV_x DelV_x

DelV_y DelV_y

DelV_z DelV_z

No No

(UTCG) (UTCG)




11

77 Jul Jul 2025 2025 23:23:34.080 23:23:34.080

0.897378 0.897378

-0.304766 -0.304766

22

77 Jul Jul 2025 2025 23:28:36.480 23:28:36.480

0.884162 0.884162

33

77 Jul Jul 2025 2025 23:20:15.360 23:20:15.360

44 55

Prob. Prob.

Table 6c. Simulation Result afterTable Avoidance Maneuver Planning Using Pattern SearchPlanning Using Pattern Search 6c. Simulation Result after Avoidance Maneuver Case

Man_Time

DelV_x

DelV_y

DelV_z

Total DelV

Best Cost Fun

Max. Func Eval.

Collision

No

(UTCG)

(km/s,ECI)

(km/s,ECI)

(km/s,ECI)

1

7 Jul 2025 23:14:21.120

0.886258

-0.283714

-0.147013

0.942104

1.6527e4

279

No Warn.

2

7 Jul 2025 23:11:36.960

0.902937

3

7 Jul 2025 23:19:58.080

0.917523

-0.296465

-0.114899

0.957282

2.1010e4

211

No Warn.

-0.298591

-0.106289

0.970722

2.8470e4

125

No Warn.

4

7 Jul 2025 23:19:58.080

0.919001

-0.297857

-0.107845

0.972067

2.8472e4

414

No Warn.

5

7 Jul 2025 23:20:24.000

0.890843

-0.295652

-0.119651

0.946217

1.5324e4

896

No Warn.

6

7 Jul 2025 23:14:21.120

0.899524

-0.282099

-0.148789

0.954391

1.8371e4

112

No Warn.

7

7 Jul 2025 23:19:58.080

0.915705

-0.293800

-0.119325

0.969058

2.3589e4

508

No Warn.

8

7 Jul 2025 23:12:54.720

0.886717

-0.283660

-0.147006

0.942519

1.6693e4

278

No Warn.

9

7 Jul 2025 23:19:58.080

0.899802

-0.298733

-0.110234

0.954482

2.1358e4

119

No Warn.

10

7 Jul 2025 23:19:58.080

0.915109

-0.292987

-0.121500

0.968519

2.2820e4

152

No Warn.

19 19

275

Prob.

http://ijass.org


2.8472e4. All cases produce periodic orbits. However, the pattern search converged faster than the others. Before avoidance maneuver, the predicted distance reached approximately 0.533 km and the collision probability was 1.02e-2. After collision avoidance, all of the cases detected no collision because the collision warning is set to 100km as the miss distance.

The maximum iteration for each optimization method was set to 4,000. For the simulated annealing, the minimum bestcost-function was 1.6704e4, while the maximum best-costfunction was 4.8500e4. The minimum best-cost-function was 1.353e4 using the genetic algorithm, and the maximum was 2.902e4. With the pattern search method, the minimum best-cost-function was 1.53249e4, and the maximum was

Fig. 5a. Fitness Function by Simulated Annealing

DOI:10.5139/IJASS.2014.15.3.267

Fig. 5a. Fitness Function by Simulated Annealing

276


converged more rapidly than the other cases but ended at the maximum 1000 iterations.

Figure 5.a~c show the best fitness function values for all iterations. The fitness function of the simulated annealing method decreased faster than the other cases, but it ended at the maximum iteration value. The minimum and maximum fitness function values of the genetic algorithm and the pattern search are similar, though the pattern search values converged more rapidly. The pattern search method

Fig. 5b. Fitness Function by Genetic Algorithm

6. Conclusion This study investigated an insertion maneuver considering collision avoidance against an uncontrolled

Fig. 5b. Fitness Function by Genetic Algorithm

277

http://ijass.org


Fig. 5c. Fitness Function by Pattern Search

Fig. 5c. Fitness Function by Pattern Search

the STK/Astrogator and Matlab global optimization toolbox. All three algorithms are used to obtain the maneuver start time and delta-V so as to reduce the collision probability from uncontrolled objects and to retain the reference trajectory calculated by the differential corrector with the

object in the Earth-Moon L2 quasi-halo orbit. Global optimization algorithms including simulated annealing, genetic algorithm and pattern search were used to minimize the differences between a reference orbit and an orbit after avoidance maneuver. The simulations were performed using

DOI:10.5139/IJASS.2014.15.3.267

278


[6] dos Santos, D.P.S., and Prado, A.F.B.A., Minimum Fuel Multi-Impulsive Orbital Maneuvers Using Genetic Algorithms, IAA-AAS-DyCoSS1-11-01. [7] Michael A. Paluszek, and Stephanie J. Thomas, “Trajectory Optimization Using Global Methods”, The 29th International Electric Propulsion Conference Princeton University, October 31 - November 4, 2005. [8] Lu-sha Wang, Zhi-gang Wang, and Wei Li, Optimization of Parameters Injecting into Orbits Based on the Pattern Search Method, 2009 First International Conference on Information Science and Engineering(ICISE), 2009, pp.38233826. [9] Willem H. de Vries, Maneuver Optimization through Simulated Annealing, Lawrence Livermore National Laboratory, LLNL-CONF-497728. [10] Hae-Dong Kim, and Eun-Sum Sim, “Lunar Transfer Trajectory Design Using a WSB”, Proceeding of the 2009 KSAS Fall Conference, 2009. 21-1201. [11] Michael A. Paluszek, and Stephanie J. Thomas, “Trajectory Optimization Using Global Methods”, The 29th International Electric Propulsion Conference Princeton University, October 31 - November 4, 2005. [12] MATLAB, Global Optimization ToolBox, Mathworks. [13] Sang-Cherl Lee, Hae-Dong Kim, and Jinyoung Suk, “Collision Avoidance Maneuver Planning Using GA for LEO and GEO Satellite Maintained in Keeping Area”, Int’l J. of Aeronautical & Space Sci., Vol. 13, No.4, 2012, pp. 474-483. [14] Magnus Erik Hvass Pedersen, Tuning & Simplifying Heuristical Optimization, Ph.D. Thesis for the degree of Doctor of Philosophy, January 2010 [15] Rosário Laureano, Determinism Versus Predictability in the Context of Poincare’s Work on the Restricted 3-Body Problem, University Institute of Lisbon, Portugal, Vol. 2, No. 2, March 2012.ISSN-L: 2223-9553, ISSN: 2223-9944. [16] Pierpaolo Pergola, Chiara Finocchietti, “Circular Restricted Tree Body Model”, Space System Lecture, Pisa, November 2012. [17] F.Renk, M.Landgraf, “Sun-Earth Libration point transfer options with intermediate HEO”, ActaAstronautica, Vol. 74, May-June 2012, pp. 1-19. [18] Infeld, Samantha I., Optimization of mission design for constrained libration point space missions, ProQuest Dissertations And Theses; Thesis (Ph.D.), Stanford University, 2006, Publication Number: AAI3209042; ISBN: 9780542570926; Source: Dissertation Abstracts International, Volume: 67-02, Section: B, page: 1001.; 133 p. [19] Tae Soo No, JiMarn Lee, Gyeong Eon Jeon, Daero Lee, and Ghangho Kim, “A Study on Earth-Moon Transfer Orbit Design”, Int’l J. of Aeronautical & Space Sci., Volume 13, No. 1, 2012, pp. 106-116.

STK/Astrogator. The genetic algorithm method could create new individuals by combining different solutions, while simulated annealing created a new individual by modifying only a solution with a local move. Consequently, the genetic algorithm method performed better than the simulated annealing method. From the performance point of view, genetic algorithm is a slow starter with more accuracy, while the simulated annealing method starts quickly and can obtain a relatively good solution in a short time. The convergence of the genetic algorithm took longer than that of the simulated annealing method. The genetic algorithm method uses a large search space and requires a lot of iterations for optimization, while the pattern search method finds a certain direction and requires relatively fewer steps for solving the problem. Additionally, the genetic algorithm method optimizes its search space in complicated problems, while the pattern search method displays the best performance with a relatively smaller search space and less computing time. This approach can be applied to an actual quasi-halo orbit mission for inserting a spacecraft into a halo orbit with collision avoidance maneuver for the Earth-Moon and the Sun-Earth Lagrange point mission.

Acknowledgements The research was supported by the “Electro-optic Space Surveillance System” of the National Agenda Project funded by the Korea Research Council of Fundamental Science & Technology. Authors would like to thank KRCFS and KARI for their support.

References [1] Dawn Perry Gordon., Transfers to Earth-Moon L2 HALO Orbits Using Lunar Proximity and Invariant Manifolds, M.S. Thesis, Purdue University, August 2008. [2] Bernelli F, Topputo F. and Massari M., Assessment of Mission Design Including Utilization of Libration Points and Weak Stability Boundaries, Final Report, ESTEC Contract No. 18147/04/NL/MV, 2004. [3] http://en.wikipedia.org/wiki/Lagrangian_point [4] Woodard, M., Folta, D., and Woodfork, D., “ARTEMIS: The First Mission to the Lunar Libration Points”, presented at the 21st International Symposium on Space Flight Dynamics (ISSFD), Toulouse, France, September 28 - October 2, 2009. [5] http://www.answers.com/topic/space-exploration

279

http://ijass.org


Utah State University, 2012. [25] Alfano S, “Collision avoidance maneuver planning tool”, in Proceedings of the 15th AAS/AIAA Astrodynamics Specialist Conference, Lake Tahoe, CA, 7-11 Aug 2005. [26] Alfano S, “Review of conjunction probability methods for short-term encounters”, in Proceedings of the AAS/AIAA Space Flight Mechanics Meeting, Sedona, AZ, 28 Jan-1 Feb 2007. [27] Eun-Hyouek Kim, Hae-Dong Kim, Hak-Jung Kim, “Optimal Solution of Collision Avoidance Maneuver with Multiple Space Debrisc, Journal of Space Operations— Letters, Vol. 9, No. 3, July—September 2012.

[20] Parker, Jeff., Dynamical Systems Theory Applied to Spacecraft Mission Design, Dec 6, 2012. [21] Kyle Wolma, The Use of Lagrange Points for Lunar Exploration, Settlement, and Support, ASEN5050, Colorado University. [22] STK Tutorial, Trajectory Design with STK/Astorgator : Earth-Moon L1 Trajectory Tutorial. [23] Alfano, S., “A Numerical Implementation of Spherical Object Collision Probability”, The Journal of Astronautical Sciences, Vol. 53, No. 1, January-March 2005, pp. 103–109. [24] Michael R. Phillips, Spacecraft Collision Probability Estimation For Rendezvous and Proximity and Operation,

DOI:10.5139/IJASS.2014.15.3.267

280