Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Optimization of relative motion rendezvous via modified particle swarm and primer vector theory
(USC Thesis Other)
Optimization of relative motion rendezvous via modified particle swarm and primer vector theory
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
OptimizationofRelativeMotionRendezvousvia ModifiedParticleSwarmandPrimerVectorTheory by Colter Rylan Danekas _____________________________________________________ Thesis SubmittedinPartialFulfillmentofthe RequirementsfortheDegree MasterofScience AstronauticalEngineering ViterbiSchoolofEngineering UniversityofSouthernCalifornia August2018 Copyright2018 ColterRylanDanekas Contents ListofFigures ...................................................................................... iii ListofTables........................................................................................ iv Nomenclature....................................................................................... vi 1 Introduction .................................................................................... 2 1.1 PreviousRelevantWork .................................................................. 3 1.2 Objective .................................................................................. 4 1.3 OutlineofThesis.......................................................................... 5 2 OrbitalMechanics ............................................................................. 6 2.1 ClassicalOrbitElements.................................................................. 6 2.2 RelativeMotionofTwoSpacecraft ...................................................... 8 2.2.1 EquationsofRelativeMotion.................................................... 8 2.2.2 Clohessey-WiltshireEquations .................................................. 11 3 OptimizationMethodology.................................................................... 13 3.1 InitialStateSelection ..................................................................... 13 3.2 PrimerVectorTheory..................................................................... 14 3.2.1 CostateVectorEquations ........................................................ 15 3.2.2 AdjointEquation ................................................................. 17 3.2.3 CostFunction..................................................................... 18 i 3.2.4 AddingCoasts .................................................................... 21 3.2.5 CriteriaforanIntermediateImpulse ............................................ 21 3.2.6 NecessaryConditionsforOptimality............................................ 23 3.3 StochasticOptimization................................................................... 23 3.4 DifferentialEvolution..................................................................... 24 3.5 ParticleSwarmOptimization............................................................. 25 3.5.1 BasicParticleSwarmAlgorithm ................................................ 25 3.5.2 StandardParticleSwarmAlgorithm............................................. 26 3.5.3 UnifiedParticleSwarmAlgorithm .............................................. 27 3.6 ModifiedParticleSwarm ................................................................. 28 3.7 AlgorithmVerification.................................................................... 30 3.8 ExampleProblem ......................................................................... 30 4 Results........................................................................................... 32 4.1 Two-ImpulseTrajectory .................................................................. 32 4.1.1 ModifiedParticleSwarmSolution............................................... 33 4.1.2 PrimerVectorTheorySolution .................................................. 34 4.2 Three-ImpulseTrajectory................................................................. 37 4.2.1 ModifiedParticleSwarmSolution............................................... 37 4.2.2 ModifiedParticleSwarmSolutionWithCoast ................................. 38 4.2.3 PrimerVectorTheorySolution .................................................. 41 4.3 GeneralResults............................................................................ 43 5 Summary........................................................................................ 46 5.1 Conclusion................................................................................. 46 5.2 FutureWork ............................................................................... 47 Bibliography ........................................................................................ 49 ii ListofFigures Figure2.1 Angularclassicalorbitalelementsofanarbitraryorbit .................... 7 Figure2.2 Planarclassicalorbitalelementsofanarbitraryorbit ...................... 7 Figure2.3 Local-verticallocalhorizontal(LVLH)coordinatesystem ................ 9 Figure4.1 Optimaltwo-impulsetrajectoryviaMPSO ................................. 34 Figure4.2 Optimaltwo-impulsetrajectoryviaMPSOwithcoastaddedviaPVT .... 36 Figure4.3 Optimalthree-impulsetrajectoryviaMPSO ................................ 38 Figure4.3 Optimalthree-impulsetrajectoryviaMPSOwithcoastaddedviaPVT .. 40 Figure4.4 Optimalthree-impulsetrajectoryviaPVT .................................. 42 Figure4.4 Optimalthree-impulsetrajectoryviaPVTmagnifiedview ................ 43 iii ListofTables Table3.1 Exampleproblem1MPSOparametersfortrajectoryanalysis ................. 31 Table4.1 PositionandVelocityValuesforOptimalTwo-ImpulseMPSOSolution ..... 33 Table4.2 V ValuesandTimesforOptimalTwo-ImpulseMPSOSolution ............ 33 Table4.3 PositionandVelocityValuesforOptimalTwo-ImpulseMPSOSolutionwith CoastAddedviaPVT ............................................................ 35 Table4.4 V ValuesandTimesforOptimalTwo-ImpulseMPSOSolutionwithCoast AddedviaPVT ................................................................... 35 Table4.5 PositionandVelocityValuesforOptimalThree-ImpulseMPSOSolution ... 37 Table4.6 V ValuesandTimesforOptimalThree-ImpulseMPSOSolution ........... 37 Table 4.7 Position and Velocity Values for Optimal Three-Impulse MPSO Solution withCoastAddedviaPVT ...................................................... 39 Table 4.8 V Values and Times for Optimal Three-Impulse MPSO Solution with CoastAddedviaPVT............................................................. 39 Table4.9 PositionandVelocity ValuesforOptimalThree-ImpulsePVTSolution ...... 41 Table4.10 V ValuesandTimesforOptimalThree-ImpulsePVTSolution ............ 41 Table4.11 AverageDifferenceBetweenTwo-ImpulsePVTsolutionandParticleSwarm, swarmsize=150,numberofgenerations=1000 .............................. 44 Table4.12 Mean,Minimum,andStandardDeviationofTwo-ImpulseTrajectoryCosts Achieved by Each PSO Method, swarm size = 150, number of generations =1000 ............................................................................ 44 iv Table4.13 AverageDifferenceBetweenTwo-ImpulsePVTsolutionandParticleSwarm, swarmsize=200,numberofgenerations=3000 ............................. 45 Table4.14 MeanandMinimumTwo-ImpulseTrajectoryCostsAchievedbyEachPSO Method,swarmsize=200,numberofgenerations=3000 ................... 45 v Nomenclature x positionofthedeputyinthex-direction y positionofthedeputyinthey-direction z positionofthedeputyinthez-direction positionvectorofthedeputyrelativetothechief x velocityofthedeputyinthex-direction y velocityofthedeputyinthey-direction z velocityofthedeputyinthez-direction v velocityvectorofthedeputyrelativetothechief x statevectorofthedeputyrelativetothechief x accelerationofthedeputyinthex-direction y accelerationofthedeputyinthey-direction z accelerationofthedeputyinthez-direction v changeinvelocityvectorduetoimpulsivethrust v * velocityvectorofthedeputyrelativetothechiefjustpriortoinstantaneousthrust v + velocityvectorofthedeputyrelativetothechiefjustafterinstantaneousthrust a semi-majoraxisoftheorbit e eccentricityoftheorbit i inclinationoftheorbit ! argumentofperigeeoftheorbit rightascensionoftheascendingnodeoftheorbit f trueanomalyofspacecraftonorbit vi n meanmotion(angularvelocity)ofchief h specificangularmomentumvectorofspacecraft CW statetransitionmatrixforClohessy-Wiltshireequations costatevector p primervector q companionvector H Hamiltonianfunction t 1 timeoffirstimpulse t m timeofintermediateimpulse t 2 timeofsecondimpulse totaltimeofflightforrendezvousmaneuver C C cognitivecoefficientusedinparticleswarmoptimizationalgorithm C S socialcoefficientusedinparticleswarmoptimizationalgorithm p best bestvalueexperiencedbyagivenparticle g best bestvalueexperiencedwithinaneighborhoodoftheswarm G best bestvalueexperiencedbytheentireswarm constrictioncoefficient u unificationfactor CR crossoverprobability F weightingcoefficient randomnumbertakenfromauniformdistribution randomnumberfromuniformdistribution vii Abstract Rendezvousproximityoperationsareacrucialcomponentofspacecraftmissions,andwill become even more so for Earth-based missions as the number of satellites and orbital debris in- creases. By the optimization of these rendezvous operations, a mission’s life may be extended substantially. In addition, more missions may become possible due to this higher level of effi- ciency. Amodifiedparticleswarmalgorithmisproposedtoidentifyfuel-efficienttrajectoriesmore quickly, with a higher level of accuracy, and with a higher level of reliability or stability. The Clohessy-Wiltshire equations are used to describe the relative motion between the two spacecraft. The results of the modified particle swarm algorithm are then improved, if possible, by primer vectortheory(PVT)toensureoptimalityofthemaneuver. The proposed modified particle swarm optimization (MPSO) algorithm is proposed as a more sophisticated alternative to the basic particle swarm optimization (PSO) algorithm. The MPSO algorithm incorporates a component of mutation, based on a differential evolution (DE) algorithm, to help avoid stagnation and improve performance. The proposed MPSO algorithm successfullyachievednearoptimaltrajectorieswhencomparedtothoseachievedviaPVT.Thepri- mary benefits of the MPSO over its simplified PSO and DE counterparts is it’s improved ability inpreventstagnationinamultimodalsolutionspaceandenhancedlocalsearch,orexploitation,of the solution space. These benefits result in the MPSO algorithm locating optimal, or near optimal solutions with a higher level of efficiency and with a higher consistency. The higher level of effi- ciency is achieved by locating these near optimal solutions with a smaller swarm size and number of generations required when compared to both the PSO and DE algorithms, and the higher con- sistencyisachievedbysimplyanalyzingthemeanandstandarddeviationofsolutionsoveralarge numberofsimulations. 1 Chapter1 Introduction Satellite rendezvous proximity operations are a vital component of the future of the space industry. These proximity operations are relative motion maneuvers of one spacecraft relative to another. Optimizing these proximity operations, particularly if many are required, can extend a mission’slifesignificantly. Thereareseveraltechniquesusedtooptimizedthesetrajectories,which includebothdirectandindirectmethodologies. Theobjectiveistominimizethemagnitudeofthe sumoftheburns,i.e. ñv total ñ. Thiscanberelatedtotheamountoffuelrequiredforthemaneuver via the rocket equation. Minimizing the total burn cost, can ultimately lead to a longer and/or a cheapermission,bothhighlydesirablecharacteristicsofamission. The initial motivation for this work comes from proposed mission by Bento da Silva and Barnhartinwhichasinglespacecraftwouldrendezvouswithmultiplerocketbodies,attachadevice to each rocketbody that would induce orbital decay of the rocketbody and ultimately causes the rocketbody to burn up upon re-entry [1]. This is a novel approach, and given the rocketbodies selected, relative motion rendezvous proximity operations would be vital to the mission’s success. Themethodologyproposedinthisworkisrobustandcanbeappliedtomission’soutsideoforbital debrisremediationmissions. In satellite proximity operations the motion of an active satellite is typically described rel- ative to a stationary satellite. The active and stationary satellites are commonly referred to as the 2 “deputy”satelliteand“chief”satelliterespectively. Thisterminologywillbeusedintheremainder of this work. To describe the relative motion between the deputy and a chief in a circular orbit, a special set of linear equations are commonly used, which are known as the Clohessy-Wiltshire (CW) equations. These equations enable a quick and highly efficient means of analyzing relative motiontrajectories. The optimization is performed using a modified particle swarm algorithm (MPSO). This algorithm is a part of a larger collection of stochastic optimization techniques that directly deter- mine the best possible value obtainable. Stochastic optimization techniques do not guarantee that theglobalminimumwillbelocated,andthereforetheparticleswarmalgorithmisjuxtaposed,and improved, if possible, with primer vector theory, which is an indirect technique optimization tech- niqueconstructedusingcalculusofvariations. Ifthesolutionobtainedusingthemodifiedparticle swarm algorithm is sufficiently close to the solution determined via primer vector theory, then the particle swarm solution is said to be optimal or near optimal. The term of sufficiently close is subjecttoeachmission’srequirementsandconstraints. 1.1 PreviousRelevantWork Previousworkhasbeenperformedinoptimizingimpulse-thrusttrajectoriesviaprimervec- tor theory when utilizing the CW equations in the local-vertical local-horizontal frame. Jezewski providesathoroughanalysisofapplyingprimervectortheoryandprovidingconditionsforoptimal- ity of the transfer [2]. However, few details are included in the numerical examples leaving holes thatthereadermustfillin. Conwayalsoprovidesadetailedanalysisofdeterminingconditionsfor when initial and final coasts will lower the transfer cost. In addition, Conway also provides a way ofdeterminingtheperturbationorvariationofthepositionvectorwhenexaminingtheadditionof anintermediateimpulse[3]. Theparticleswarmoptimization(PSO)algorithmhasbeenmodifiedandimprovedsinceit wasfirstintroducedbyKennedyandEberharttomodelbirdsflockingandfishschooling[4]. Clerc 3 and Kennedy built upon the original work and further improved its stability and performance by introducinganewcoefficient. Finally,ParsopoulosandVrahatisfurtheredtheworkdonebyClerc and Kennedy by further exploiting the desired characteristics of the algorithm and fine-tuning the parameterstoimproveconvergenceaccuracyandconvergencerate[5]. Very little work has been done in the application of evolutionary algorithm in optimizing these rendezvous proximity operations. Even less work has been done in the means of applying particle swarm optimization algorithms to this type of problem. Aubin applied the basic particle swarm algorithm to an orbit-to-orbit transfer for work relating to satellite clusters. The work by Aubindidsuccessfullyshowthattheparticleswarmalgorithmcanlocatenear-optimaltrajectories determined by primer vector theory, however the analysis done required a large swarm size (500 particles) and number of generations (2000 generations) to be executed to achieve optimal results [6]. Additionally, this work does not investigate the fuel-minimizing trajectory for a rendezvous proximityoperation,whereafinalcoastingarcisnotpossible. 1.2 Objective Theobjectiveofthisworkistofurtherdevelopaneffective,reliable,andefficientprocedure for identifying fuel-minimizing rendezvous trajectories in the relative motion frame through the utilization of a MPSO technique. Primer vector theory, a well established means of checking for optimality,isthenusedtocomparetheresultsachievedviaMPSO,andfurtheroptimizetheMPSO results, if possible. If the optimal MPSO solution marginally differs from the true optimal value thentheMPSOsolutioncanbeconsideredasnearoptimal,andthusaeffectivetoolforidentifying optimalfuel-minimizingtrajectories. Thebasic(PSO)algorithmisanadvantageousoptimizationtechniqueduetoitseaseofim- plementationandstochasticdesign,whichenablesthealgorithmtoidentifyglobalextrema,ornear global extrema, when the solution set contains a plethora of local extrema. This results in a larger portionofthesolutionspacebeingsearchedopposedtonon-stochasticmethod. TheMPSOdiffers 4 fromthebasicPSOinitsconstruction,whichismoresophisticatedandexploitsthecharacteristics of the PSO algorithm. This in turn can lead to faster convergence, while maintaining a high level ofaccuracy,anadvantageouscharacteristicforanyoptimizationtechnique. While some previous work has been done in optimizing relative trajectories using evolu- tionary algorithms, very little work has been done utilizing PSO. Aubin optimized orbit transfers between periodic orbits via PSO, however large computational resources were used as previously discussed[6]. Unfortunately,nocomputationalresources,suchasprocessorused,numberofcores, machine, etc. were included in Aubin’s work. This work builds upon that work and enables near optimal trajectories to be achieved utilizing a smaller swarm and number of generations required. Achievingthiswouldequatetoafaster,morereliable,andmoreefficientmeansofcomputingfuel- optimaltrajectoriesviathePSOalgorithm. 1.3 OutlineofThesis Chapter2includesderivationsofthegoverningequationsusedinthiswork: theorbitalele- ments,relativemotionbetweentwosatellites,andtheClohessy-Wiltshire(CW)equations. Chapter 3discussesthegeneraloptimizationmethodologyusedinthiswork,detailsofprimervectortheory and how the basic PSO algorithm works, some of the variations of the PSO algorithm, and finally the construction and methodology of the MPSO algorithm. Chapter 4 discusses the results of an example problem and some general results. Finally, Chapter 5 draws conclusions on the work and givesrecommendationsforfuturework. 5 Chapter2 OrbitalMechanics 2.1 ClassicalOrbitElements To describe the position and velocity of a spacecraft in orbit a coordinate system must be selected in order to provide a reference. The Cartesian frame is a common selection in many ap- plications, but can be difficult to visualize the shape and orientation of an orbit. Therefore, it is often easier to work with a set of six parameters called the classical orbital elements. These six parameters replace the three positions and three velocity components required by the Cartesian frame, in order to fully define the spacecraft’s state, and define the orbit’s shape, orientation, and thepositionofthespacecraftontheorbit. Oneadvantageofusingtheclassicalorbitalelementsis thatonlyoneofthesixparametersvaries,makingiteasiertoanalyzethespacecraft’smotion. The angular orbital elements and the planar orbital elements can be seen in Figure 2.1 and Figure 2.2 respectively. Notethatthetermplanarisinreferencetobeingintheorbitalplane. 6 Figure2.1: Angularclassicalorbitalelementsofanarbitraryorbit[7] Figure2.2: Planarclassicalorbitalelementsofanarbitraryorbit 7 Thesixparametersaregivenby a: The semi-majoraxisdescribesthegeneralsizeoftheorbit. e: Theeccentricityoftheorbitdescribestheshapeoftheorbit,thelowerlimitof zeroindicatesacircle,whiletheupperlimitofoneindicatesaparabola. i: Theinclination,orangle,oftheorbitalplanefromtheequatorialplane. !: Theargumentofperigeeistheanglefromthenodallinetoperigee. : The rightascensionoftheascendingnode,istheanglefromthevernalequinox./ tothenodalline. f: Thetrueanomalyistheanglefromperigeetothespacecraft’scurrentposition intheorbit. 2.2 RelativeMotionofTwo Spacecraft 2.2.1 EquationsofRelativeMotion First, consider a chief in a circular orbit about the Earth. Then a local-vertical local hori- zontal (LVLH) coordinate system is attached to the chief and rotates with the chief’s radius vector atanangularrateequaltothemeanmotionofthechief. Thisreferenceframemayalsobereferred to as the Hill frame or the Clohessy-Wiltshire frame. In this coordinate system,x lies in the radial direction, i,awayfromtheEarth,yliesinthedirectionoftheon-trackvelocity, j,andzcompletes theright-handedorthogonalbasis, kasshowninFigure. 2.3. Notethatthe j directioninthefigure isintothepage,denotedby ¼ . 8 Figure2.3: Local-verticallocalhorizontal(LVLH)coordinatesystem IntheLVLHframethexandycomponentscorrespondtothein-planemotion,andthezcomponent referencestheout-of-planemotion. IntheLVLHframethepositionofthedeputyisgivenby r d =r c += ñr c ñ+x i+y j+z k ; (2.1) where ññ is taken as the Euclidean norm and the subscriptsc andd correspond to the chief and deputy respectively. For compactness, the subscript c will be omitted for the remainder of this work. From kinematics, the acceleration of the deputy in the chief’s frame is given by Equation (2.2). r d = r+ +./+2. /+ (2.2) The angular velocity of the chief’s frame, or the rate of rotation of the LVLH frame, is defined by Equation(2.3),wherehisthetheangularmomentumdefinedbyEquation(2.4). 9 = f k= ñhñ ñrñ 2 k (2.3) h=rv (2.4) Itfollowsthattheangularaccelerationoftheframeisgivenas = d dt ./= f k=* 2ñrñ f ñrñ 2 k : (2.5) SubstitutingEquation(2.2)andEquation(2.3)intoEquation(2.1)andsimplifyingyields r d = ñ rñ* f 2 x* fy*2 f f+ x i+ fx* f 2 y+2 f x+ y j+ z k : (2.6) TheaccelerationofthedeputyandchiefaboutthecentralbodyaregivenbyEquations(2.7). r d =* ñr d ñ 3 r d ñ rñ=* ñrñ 2 (2.7) The acceleration of the chief and the deputy can be rewritten by substituting Equations (2.7) into Equation(2.6),thusyielding * ñr d ñ 3 r d = 0 * ñrñ 2 * f 2 x* fy*2 f f+ x 1 i+ fx* f 2 y+2 f x+ y j+ z k : (2.8) This can be broken down in component form and then rearranged, thus yielding the full nonlinear equationsofrelativemotionasgiveninEquations(2.9). x*2 f y* fy* f 2 x* ñrñ 2 =* ñr d ñ 3 .ñrñ+x/ y+2 f x+ fx* f 2 y=* ñr d ñ 3 y (2.9) z=* ñr d ñ 3 z 10 The potential forces can be approximated such that a set of linear equations results as given in Equation(2.10)[7]. x*2 f y* 0 f 2 +2 ñrñ 3 1 x* fy=0 y+2 f x+ fx* 0 f 2 * ñrñ 3 1 y=0 (2.10) z+ ñrñ 3 z=0 This set of equations is called the linear equations of relative motion and provide the basis for the remainderofthiswork. 2.2.2 Clohessey-WiltshireEquations The previous section discussed the general nonlinear, and linear equations of relative mo- tion,whicharevalidforanygeneralKeplerianorbit. Thelinearequationscanbefurthersimplified when assuming the chief is in a circular orbit. The simplifications of assuming the chief is in a circularorbitaregivenbyEquations(2.11). f =0 f =n (2.11) n 2 = ñrñ 3 wheren isthemeanmotionofthechief. ThesesimplificationsyieldEquations(2.12). x*2n y*3n 2 x=0 y+2n x=0 (2.12) z+n 2 z=0 11 ThissetoflinearequationsofrelativemotionforacircularorbitarecalledtheClohessey-Wiltshire (CW) equations, or sometimes referred as the Hill-Clohessey-Wiltshire equations. In state-space form,theCWequationsaregivenbyEquations(2.13). x=Bx= b f f f f f f f f f f f d 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 3n 2 0 0 0 2n 0 0 0 0 *2n 0 0 0 0 *n 2 0 0 0 c g g g g g g g g g g g e x (2.13) Notethatthestatematrix,B,inEquation(2.13)istimeindependent. Thestate-spacematrixrelates the state vector to its derivative. However, by integrating Equation (2.13) and solving solving for constants of integration by applying the boundary conditions, it is possible to construct a state transition matrix, which can relate the state vector at one time to the state vector at another time. ThisrelationisgivenbyEquation(2.14). x.t/= CW t;t o x t o (2.14) Thisstatetransitionmatrixisgivenas CW t;t o = b f f f f f f f f f f f f f f f f f f d 4*3c 0 0 s n 2 n .1*c/ 0 6.s*nt < / 1 0 2 n .c*1/ 1 n .4s*3nt < / 0 0 0 c 0 0 s n 3ns 0 0 c 2s 0 6n.c*1/ 0 0 *2s 4c*3 0 0 0 *ns 0 0 c c g g g g g g g g g g g g g g g g g g e = b f f f f d ' rr 33 ' rv 33 ' vr 33 ' vv 33 c g g g g e (2.15) whereforcompactness,t < =t*t o ,c=cosnt < ,ands=sinnt < . 12 Chapter3 OptimizationMethodology 3.1 InitialStateSelection Inordertoeffectivelyutilizetheparticleswarmalgorithmtheinitialpositionandvelocityof the deputy were randomly constructed from a uniform distribution. The position was constructed such that the magnitude of the initial distance between the deputy and chief was d. A uniform distributionwasselectedoveranormaldistributiontodiscouragetheclusteringofinitialpositions, andallowforalargerdomainofthesolutionspacetobesearched. Thepositionandvelocityvectors ofthedeputyrelativetothechiefaregivenasdefinedbyEquations(3.1). =[x;y;z] andv=[ x; y; z] (3.1) ThecomponentsofthepositionvectoraredefinedbyEquation(3.2). x=d 2 1 *1 y=d 2 2 *1 (3.2) z=d 2 3 *1 Where, i = 1, 2, 3 are independent random numbers between zero and one taken from uniform 13 distribution. The velocity components must be scaled in such a way that the initial velocity is not so large to cause the deputy’s trajectory to leave the CW frame, but to give the maximum number of potential solutions. Aubin derived a scaling factor based on the maximum potential difference betweenthechiefanddeputy[6]. Thisscalingfactor,,isgivenbyEquation(3.3). = ó ó ó ó ó u a+ññ * u a ó ó ó ó ó (3.3) Thecomponentsofthevelocityvector arethenconstructed,andgivenbyEquations(3.4). x= 2 1 *1 y= 2 2 *1 (3.4) z= 2 3 *1 Again,,i=1,2,3areindependentrandomnumbersbetweenzeroandonetakenfromuniformdis- tribution. The random numbers used here are independent of those used in computing the random position. The index,i, simply indicates that the random number used in each velocity component isindependentoftheotherrandomvalues. 3.2 PrimerVectorTheory Firstly,thestatevectorisgivenasdefinedinEquation(3.5). x T = T ; v T (3.5) where andv are given by Equations (1.1). Note that here the superscriptT indicates the trans- pose of the vector. Recall from Section 2.2.2 that the state vector at some initial time,t o , can be propagatedthroughtimeviathestatetransitionmatrix. ThisisshowninEquation(2.14). Thestate vectorcanbedifferentiatedwithrespecttotime,whichyieldsEquation(3.6). 14 x= d dt .x/=[ x; y; z; x; y; z] T (3.6) SubstitutingintheCWequationsofmotionfromEquation(2.12)yieldsEquation(3.7). x= x; y; z; 2n y+3n 2 x; *2n x; *n 2 z T (3.7) 3.2.1 CostateVectorEquations ALagrangemultiplier,amathematicaltermincorporatedintotheproblemtohelpwiththe optimization process, can be incorporated into the state equations to create the costate equations, whichcreateasystemofequationsadjointtothestateequations. Thedifferentialformofthecostate equationsaregivenasdefinedinEquation(3.8). =* 0 ) x )x 1 T (3.8) AlgebraicallysolvingEquation(3.8)andrearrangingyieldsthesystemofequationsgivenbyEqua- tions(3.9). 1 =*3n 2 5 2 =0 3 =n 2 6 4 =* 1 +2n 5 (3.9) 5 =* 2 *2n 4 6 = 3 15 Eqns. (3.9)caneasilybeintegrated,thusyieldingEquations(3.10). 1 =a 1 c+a 2 s+2 a 3 +a 4 t 2 =* a 4 3n 3 =n a 5 s*a 6 c 4 =* 2 3n a 1 c+a 2 s * 1 n a 3 +a 4 t (3.10) 5 = 1 3n a 1 s*a 2 c* 2 n a 4 6 =a 5 c+a 6 s Wherea i ,.i=1;§;6/,areconstantsthataredeterminedbytheboundaryconditionsoftheprob- lem. Heres = sin.nt/ andc = cos.nt/. The costate vector,, can then be defined by Equation (3.11). T = q T ; p T (3.11) Thevector,p,iscalledtheprimervectorandisusedindeterminingoptimalityofthesolu- tion,andthevectorq,iscalledthecompanionvector. ThesevectorsaredefinedbyEquation(3.12). The primer vector,p, at the time of the impulse is defined to a unit vector in the direction of the impulse,ormathematicallythisisgivenbyEquation(3.13). q= 1 ; 2 ; 3 T and p= 4 ; 5 ; 6 T (3.12) p= v ñvñ (3.13) Consider that the primer vector is taken at two times, t 1 and t 2 , such that t 1 < t 2 . Now takingp T 1 =p T t 1 = l 1 l 2 l 3 andp T 2 =p T t 2 = m 1 m 2 m 3 , then the coefficients,a i can bevectorizedandcomputedusingEquation(3.14). 16 b f f f f f f f d a 1 a 2 a 3 a 4 c g g g g g g g e = *1 b f f f f f f f d l 1 l 2 m 1 m 2 c g g g g g g g e (3.14) where = b f f f f f f f d * 2 3n 0 * 1 n 0 0 * 1 3n 0 * 2 3n 2 * 2c 3n * 2s 3n * 1 n * n s 3n * c 3n 0 * 2 3n 2 c g g g g g g g e : (3.15) Theout-of-planecomponents,a 5 anda 6 ,canbesolvedindependentlyandaredeterminedbyEqua- tions(3.16). a 5 =l 3 a 6 = m 3 *l 3 c s (3.16) Ifthevectoraisknown,thentheprimervector,pandtherelatedvectorqcanbedetermined analytically for any timet, which can enable faster and less computationally intensive calculation oftheprimervectoropposedtonumericallycomputingit. Itshouldbenotedthata 6 ,andtherefore theresultingout-of-planesolutionisonlyvalidforn N; N =0;1;2;§. 3.2.2 AdjointEquation Inordertocontinuewiththedevelopmentofthecriteriaofthenecessaryconditionsforan optimal,N-impulsesolutionandforimprovingaCWtrajectoryanothermathematicalrelationship is required. This relationship relates the costate vector and the perturbations in the state. This 17 differential equationisgivenbyEquation(3.17). x= 0 ) x )x 1 x (3.17) Where is the contemporaneous variational parameter. Since x is linear with respect to the state vector, then the above differential equation represents the exact variation. Equation (3.17) and Equation (3.8) can be manipulated to yield a useful result. These manipulations include premul- tiplying Equation (3.17) by T and Equation (3.8) byx T . Performing these manipulations and summingtogetheryieldsEquation(3.18). T x+x T = d dt T x = T 0 ) x )x 1 x*x T 0 ) x )x 1 (3.18) Thevalueoftherighthandsidecanreadilybecomputedtobezero,andthereforetheabovecanbe writtenmoresuccinctlyasseeninEquation(3.19). d dt T x =0 (3.19) This implies that the dot, or scalar, product of the costate vector and the variational state vector is constantoveranycoastingarc,asseenbyEquation(3.20). T x=C; C ËR (3.20) TheaboveequationiscommonlyreferredtoastheadjointequationfortheCWequations. 3.2.3 CostFunction Thecostfunction,J,issimplythesumofthemagnitudesoftheimpulsevectors,whichcan beseenin Equation(3.21). 18 J = ñv 1 ñ+ k É i=1 ñv m i ñ+ñv 2 ñ : (3.21) Wherev 1 ;v m i .i=1;§;k/;v 2 aretheimpulsevectorsattimest 1 ;t m 1 ;§t m k ;t 2 suchthatt 1 < t m 1 <§<t m k <t 2 . Theimpulsevectorisdefinedasv=v + *v * . Forsimplicityathree-impulse solution will be assumed, however the solution can readily be reduced to a two-impulse solution or extended for multiple intermediate impulses. Therefore, Equation (3.21) reduces to Equation (3.22). J = ñv 1 ñ+ñv m ñ+ñv 2 ñ (3.22) TakingthevariationyieldsEquation(3.23). J = v T 1 ñv 1 ñ v 1 + v T m ñv m ñ v m + v T 2 ñv 2 ñ v 2 (3.23) Utilizingthedefinitionoftheprimer vector,theaboveiswrittenasseeninEquation(2.24). J =p T 1 v + 1 *v * 1 +p T m v + m *v * m +p T 2 v + 2 *v * 2 (3.24) Wherethesuperscriptminusandpluscorrespondtotheevaluationofthevectorimmediatelybefore andaftertheimpulse,respectively. Notethatp T 1 v + 1 andp T 2 v * 2 canbewritteninanalternativeform viatheadjointequation. ThisyieldsEquations(3.25). p T 1 v + 1 =q * m T r * m +p * m T v * m *q T 1 r + 1 p T 2 v * 2 =q + m T r + m +p + m T v + m *q T 2 r * 2 (3.25) UsingEquations(3.24)andsubstitutingthemintoEquation(3.23)yieldsEquations(3.26). 19 J =*q T 1 r + 1 *p T 1 v * 1 +q T 2 r * 2 +p T 2 v + 2 +q * m T r * m *q + m T r + m (3.26) Note that since the primer vector is continuous at an interior, or intermediate impulse, thenp * m = p m =p + m ,thereforethetermsinvolvingp m havebeensimplifiedandcanceledoutofEquation(3.23). Thevariationofthepositionandvelocityvectorscanbeapproximatedtothefirstorderas r=dr*vdt and v=dv* vdt (3.27) usingtheresultsofEquation(3.25)andevaluatingthevariationsofthepositionandvelocityvectors attimest + 1 ; t * m ; t + m ; andt * 2 yields J =* q T 1 dr 1 +p T 1 dv 1 + q T 2 dr 2 +p T 2 dv 2 + q T 1 v + 1 +p T 1 v * 1 dt 1 * q T 2 v * 2 +p T 2 v + 2 dt 2 (3.28) + H + m *H * m dt m * q + m *q * m T dr m : WhereHis theHamiltonianfunction definedbyEquation(3.29). H =p T v+q T r (3.29) Thegeneralcostfunctionforthetime-openorbit-to-orbittransferisgivenbyEquation(3.30). Time- open implies that the total flight time is not fixed and may be adjusted accordingly to reduce the cost. This ability to adjust the total flight time may or may not be possible given the constraints of theproblemathand. J =*H * 1 d o +H + 2 d f +H + 1 dt 1 *H * 2 dt 2 + H + m *H * m dt m * q + m *q * m T dr m (3.30) Where d o and d f are the time durations of the initial and final coasts. The rendezvous cost function,thecostfunctionrelevanttothiswork,canthenbeconstructedusingEquation(3.30),the definitions of the Hamiltonian function, and the fact thatv + i *v * i = ñv i ñP i . This cost function 20 isdefinedbyEquation(3.31). J =ñv 1 ñ p T 1 q 1 dt 1 +ñv 2 ñ p T 2 q 2 dt 2 (3.31) H + m *H * m dt m * q + m *q * m T dr m 3.2.4 AddingCoasts Inageneraltransfer,aninitialcoast,afinalcoast,orbothmaybeaddedtothethetrajectory toreducethecostofthetransfer. UsingEquation(3.31)theoptimaldepartureandarrivaltimesof thetransferareachievedwhenp T 1 andq 1 areorthogonal,andp T 2 andq 2 areorthogonalrespectively. Therefore, the cost of the transfer can be decreased by propagating the product ofp T q forward or backwardsintimeuntilitachievesavalueofzero. Conwayalsoshowedthatif p 1 >0aninitialcoast (given byd o >0) will lower the cost. Similarly, if p 2 <0 a final coast (given byd f <0) will lower thecost[3]. 3.2.5 CriteriaforanIntermediateImpulse In addition to adding coasts to a trajectory, an intermediate impulse may be added as well to decrease the cost of the transfer. Adding intermediate transfers is more complex than adding coasts,astheoptimaltime,location,impulsemagnitude,andimpulsedirectionmustbecomputed. In order to determine these parameters consider a two-impulse trajectory with no added coasts. ThenthecostfunctionisgivenbyEquation(3.32). J =v 1 +v 2 (3.32) An intermediate impulse is then added, therefore the variation in the cost is given by Equation (3.33). dJ =p T 1 dv 1 +v m +p T 2 dv 2 (3.33) 21 TheanalysiscompletedbyLionandHandelsmanyieldsEquation(3.34)[8]. dJ = ñv m ñ 0 1*p T m v m ñv m ñ 1 (3.34) If the dot product between the primer vector and the impulse vector at the intermediate impulseisgreaterthanonethanthantheperturbedtrajectorywillresultinalowercost,andtherefore a intermediate impulse will improve the two impulse trajectory. If this dot product is greater than one thenp m >1. Therefore an intermediate impulse added at any point whenp m >1 will lower the cost. Thelargestdecreasewilloccuriftheintermediateimpulseisappliedwhenp m isamaximum. Thetimeoftheintermediateimpulsehasbeendetermined,nowthelocation,directionand magnitudeoftheimpulsemustbedetermined. Thedirectionoftheimpulsemustbeparalleltop m , as this will ensure the maximum value of the dot product between the primer vector and impulse vector. TheanalysisdonebyHughes,Mailhe,andGuzmanyieldsEquation(3.35)[9]. r m =M *1 p m ñp m ñ (3.35) wherer m is the variation of the perturbed trajectory from the nominal, two-impulse trajectory, isasmall,positivevalue,andthematrixM isgivenbyEquation(3.36). M =' vv t 2 ;t m *1 ' rv t 2 ;t m *' vv t m ;t 1 ' rv t m ;t 1 *1 (3.36) Inorderto havev m paralleltop m ,itisnecessarythat v m = p m ñp m ñ ; (3.37) where is some small, positive number. The new position of the impulse is computed by adding the deviation to the nominal position. The velocity vectorsv * m andv + m can be computed via the state transition matrix, CW . Then, the perturbed trajectory is then computed and trajectory cost isanalyzed. IfdJ>0thenthevalueof mustbedecreaseandtheprocessisreiterateduntilthethe 22 valueofdJ isnegative,andtherefore resultinginalowertrajectorycost. 3.2.6 NecessaryConditionsforOptimality Using this methodology, Lawden derived four necessary conditions for optimal N-impulse trajectories[3]. Theyareasfollows 1. pand parecontinuousontheentireinterval. 2. ñpñf1,andñpñ =1atimpulses. 3. Ateachimpulse,theprimervectorisaunitvectorintheoptimalthrustdirection. 4. Asaresultoftheabove, dñpñ dt = ñ pñ= p T p=0ataninteriorimpulse. 3.3 StochasticOptimization Stochasticoptimizationtechniquesutilizerandomnesstoidentifyalocal,orglobalextrema. In this work, the a global minima is being determined. Stochastic optimization techniques are commonly used for their ease of implementation and enhanced ability of finding global extrema when many local extrema exist. Stochastic optimization techniques work to balance exploration and exploitation. The former is responsible for locating the most promising regions of the search spaceintermsoffindingtheglobalminima. Theexploitationcomponentpromotesconvergenceand relatestofinelysearchingthesearchspaceimmediatelysurroundingagivenpoint. Anunbalanced algorithmcanleadtoprematureconvergence,theglobalextremanotbeinglocated,andmayrequire a larger population size or number of iterations to counteract these potential issues. One negative aspectofstochasticoptimizationtechniquesisthattheydonotguaranteeconvergencetotheglobal minima, and therefore further analysis may be necessary to conclude, or verify if the converged valueisindeedtheglobalminima. 23 3.4 DifferentialEvolution DifferentialEvolution(DE)isastochasticoptimizationtechniquethatutilizesa"survivalof the fittest" approach to identify the global minima of a function in a given domain. The algorithm incorporates a mutation, which in turn is used to create a test vector. The said test vector will then be compared to the original vector and the vector yielding a lower cost value will remain in the population, and the other is discarded. The construction of the mutated vector is defined by the followingprocess. First, the crossover probability,CR, and the weighting coefficient,F, are prescribed such that CR Ë [0;1] and F Ë [0;2]. Typically, these variables are assigned values of 0.9 and 0.8 respectively, as these values promote a high mutation probability. Then each particle, or sample, x j i , is assigned a random number, j i , which is taken from a uniform distribution from zero to one. Thenewmutatedposition,z j i ,isgivenbyEquation(3.38)[12]. z j i =a j i +F b j i *c j i (3.38) Thetestvectory j i isdefinedbyEquation(3.39). y j i = h n l n j x j i ; j i >CR z j i ; j i fCR (3.39) Thevectorsa j i ,b j i ,andc j i arerandomlyselectedsolutionvectorstakenfromthejthiteration. Itis importantthatx j i ,a j i ,b j i ,andc j i mustnotbeequaltooneanotherinordertoensureamutatedsolu- tionsignificantlydifferentthanthepresentsolution,x j i . Ifthetransfercostofthemutatedposition islowerthanthereferencepositionthenthemutatedpositionreplacesthereference,otherwisethe mutatedpositionisdiscarded. Theresultofthisalgorithmisthattheparticlesconvergetothebest obtainedposition,orthepositionthatyieldsthelowestcost. Thisalgorithmcanbeappliedtoboth thepositionandvelocityvectorsofthedeputytodeterminetheoptimalinitialstate. 24 3.5 ParticleSwarmOptimization Particle Swarm Optimization (PSO) is a stochastic optimization technique first used by Kennedy and Eberhart to model the social behavior of birds flocking and fish schooling [4]. PSO isanattractiveoptionwithinthestochasticoptimizationcollectionduetoitssimplicitywhencom- paredtootherevolutionaryalgorithms,suchasageneticalgorithm. PSOtakesapopulation,called a swarm, of randomly constructed particles in the solution state space, evaluates their fitness, then modifiestheparticleswithinthesolutionspacetoimprovetheirfitnessforeachgeneration,oritera- tion. Thealgorithmthenmodifieseachparticleforthesubsequentgenerationbasedontheparticle’s bestachievedvalueandthebestvalueachievedbytheentireswarmduringanygeneration. 3.5.1 BasicParticleSwarmAlgorithm The initial PSO algorithm developed by Kennedy and Eberhart was rather simplistic, yet effectiveforalargerangeofproblems. Ifconsideringtheoptimizationofaparticle’sposition,then theupdate forthebasicPSOalgorithm(BPSO)isgivenasdefinedinEqns. (3.40-3.41). v j+1 i =v j i + C C 1 p best;i *x j i + C S 2 G best *x j i (3.40) x j+1 i =x j i +v j+1 i (3.41) Wherex andv are the particle’s position and velocity respectively, i andj are the particle’s in- dex in the swarm and the index of the generation respectively, x, p best;i andG best correspond to the i th particle’s current value, the best value that the i th particle has experienced, and the best value experienced by and particle in the swarm respectively,C C andC S are weighted coefficients corresponding top best;i andG best respectively, and referred to as the cognitive and social learning coefficients respectively. These coefficients are typically both taken to be approximately two. Fi- nally, 1 and 2 areindependentrandomnumbersfromzerotoone,whicharetakenfromauniform distribution. There are several stopping criteria, namely when a certain fitness level is reached, or 25 when a prescribed number of generations have been completed. Note that the PSO algorithm can be applied to both scalars and vectors. The only required change for the vector case is to ensure thatdimensionsareconsistentbetweenthedifferentparameters. Theparametersdeterminedbythe algorithm are initial state of deputy,x o and the time of flight,. In general, the final state of the deputyneedstobedetermined,butinaproximityoperationrendezvoustheseareprescribedbased onthemission’sdesign. 3.5.2 StandardParticleSwarmAlgorithm Several variations and modifications have been made to the BPSO algorithm to improve performance and convergence accuracy. Clerc and Kennedy established one such modification, which is known as the default contemporary PSO variant [10]. This modified algorithm will be referredtoasthestandardPSO(SPSO)algorithmfortheremainderofthiswork. Thismodification incorporates more chaos and an inertial coefficient to assist in preventing stagnation of the parti- cles, as well as preventing the velocities from becoming so large that the particles are no longer converging. TheSPSOvelocityisupdatedasasdefinedinEquation(3.42). v j+1 i = v j i + C C 1 p best;i *x j i + C S 2 G best *x j i (3.42) Where is called the constriction coefficient or constriction factor. The remainder of parameters are defined in the same manner as the BPSO algorithm. The position is then updated by Equation (3.41). TheconstrictioncoefficientisdefinedbyEquation(3.43). = 2 ó ó ó 2*'* ù '.'*4/ ó ó ó (3.43) Where'=C C +C S ,and'>4. 26 3.5.3 UnifiedParticleSwarmAlgorithm One note to make regarding the algorithms previously mentioned is that the best global valuegreatlyaffectsthestepsizeforeachparticlefromonegenerationtothenext. Thismaywork in a negative manner for the worst performing particles as they may take a step size sufficiently large that they may step over a even better value than what the best known global value is. To helplimitstepsizesandtoensureabalancebetweenexplorationandexploitationParsopoulosand Vrahatis proposed the unified particle swarm optimization (UPSO) scheme [5]. The first step of thealgorithmfollowsEquation(3.42),butisdenotedslightlydifferentasgivenbyEquation(3.44). G j+1 i = v j i + C C 1 p best;i *x j i + C S 2 G best *x j i (3.44) This equation incorporates the exploration component to the particle swarm algorithm. However, the UPSO algorithm differs from the SPSO algorithm in the fact that it incorporates a term that utilizes the exploitation characteristic of local search methods. This term is given by Equation(3.45). L j+1 i = v j i + C C 3 p best;i *x j i + C S 4 g best *x j i (3.45) Whereg best is the local best value achieved within a given neighborhood or subset of the swarm. Commonly a neighborhood is defined by a distance metric from one particle to the others around in a given radius. However, this can lead to clusters of particles and lower diversity in searching the solution space. Therefore the neighborhoods in this work are defined by their index, almost creating a set of smaller swarms within the total swarm. This method encourages that several of thebestperformingparticlescancommunicate,ratherthanjustthesinglebest. Thevelocityisthen updatedby Equation(3.46). v j+1 i =.1*u/L j+1 i +uG j+1 i (3.46) 27 Whereu is called the unification factor, which acts as a control for the balance between the global and local velocity update. The position is then updated using Equation (3.41). Parsopoulos and Vrahatisimprovedupontheirownworkaddanotherstochasticparameter,,toupdatethevelocity as defined in Equations (3.47). The parameter is taken from a Gaussian distribution. The pa- rameterismultipliedbythelocalorglobaltermbasedontheschemeprescribed. Inaddition,they implemented variance in the unification factor to further improve the UPSO algorithm. There are numerous schemes that let the unification vary from zero to unity. For example, one such scheme istolettheunificationfactorvarylinearlygivenbyEquation(3.48). [11]. v j+1 i =.1*u/L j+1 i +uG j+1 i v j+1 i =.1*u/L j+1 i +uG j+1 i (3.47) u j = j J max (3.48) WhereJ max isthemaximumnumberofprescribedgenerations,oriterations. Byvaryingtheunifica- tionfactor,thealgorithmpromotesexplorationinfindingthemostpromisingareasinthebeginning andthenswitchestopromotingexploitationtowardsthelatergenerationstofinelysiftthroughthe solutionspacetolocatetheoptimalvalue. 3.6 ModifiedParticleSwarm Both DE and PSO algorithms have their respective strengths and weaknesses, thus by cre- ating a hybrid of the two or a modified PSO (MPSO) algorithm the desired properties of each can beobtainedwhilediscouragingormitigatingthenegativeproperties. Firstaslightmodificationto theunifiedPSOalgorithmasgivenbyEquations(3.49-3.50) 28 G j+1 i = v j i + C C 1 ä p best;i *x j i + C S 2 ä G best *x j i (3.49) L j+1 i = v j i + C C 3 ä p best;i *x j i + C S 1* 3 ä g best *x j i (3.50) Now 1 , 2 and 3 arenowvectorsofthetheappropriatesize. Also,1isavectorofonesand has the same dimensions as 3 andä signifies element-wise vector multiplication. This increases the stochastic nature of the algorithm, a desired property when using a stochastic optimization technique. The velocity is then updated by Equation (3.51), the position update follows Equation (3.41),andtheunificationfactorisdefinedbyEquation(3.52). v j+1 i = 1*u j L j+1 i ä 5 +u j G j+1 i ä 1* 5 (3.51) u j =sin H 2 jmod j crit+1 j crit I (3.52) Wheremod is the modulo operation andj crit is a prescribed value for the number of times theunificationfactorvariesfromzerotoone. TheDEcomponentalsohasbeenmodifiedforfurther performanceimprovement. Thesamegeneralprocedureisused,howeverF andthemutatedvector arefoundinaslightlydifferentway. FirstF isgivenbyEquation(3.53). F j i =.0:8/1+0:05& (3.53) WhereF j i nowvariesforeachparticleinthepopulationandeachgenerationoriteration. In addition,& is vector with each element taken from a normal distribution from [-1,1]. The mutated vector,z j i isthencomputedusingEquation(3.54). z j i =G best +F j i ä p best;i *a j i + 1*F j i ä b j i *c j i (3.54) 29 Similar to the traditional DE algorithm, if the mutated vector produces a lower cost then it is used, and the initial vector is discarded. However, now if the mutated vector produces a higher cost,thenthenewvectoristheinitialvectorupdatedaccordingtotheMPSOprocedurepreviously outlined. Theaimofthisprocedureistoensuretheparticlemovestowardsabettersolutionateach iteration. This in turn results in a lower number of iterations required and a lower computational time. 3.7 AlgorithmVerification ToverifytheMPSOalgorithmisindeedanimprovementupontheotherPSOandDEalgo- rithms,allthealgorithmsweretestedonseveralcommonlyusedbenchmarkfunctions. Thediffer- entalgorithmscanbethenbejuxtaposedbasedontheirrelativeperformanceforeachtestfunction. Thisprocesscanidentifyasetoffunctionsinwhichoneormorealgorithmsaresuperiortotherest. The benchmark functions were the Ackley, Rosenbrock, Booth, and Drop Wave function. Each of theaforementionedbenchmarkfunctionsposesadifferentchallengetotheoptimizationalgorithm. 3.8 ExampleProblem To analyze the effectiveness of the MPSO algorithm, an example problem is proposed as- suming the chief has a semi-major axis value of 7000 km, and an eccentricity of zero. The initial stateisfreeandisdeterminedbytheMPSO,howeverthefinalstateisprescribedtobezerovector, i.e. x./ T = ./ T ; v./ T =[0; 0; 0; 0; 0; 0] T : (3.55) Notethatbyprescribingthisfinalstatethatafinalcoastisnotpossibleasthespacecraftisassumed tomaintainthisstateforafinitetimeratherthanforaninstant. Thiscanreadilybeseenbyapplying Newton’s laws of motion. The initial conditions determined by the two-impulse MPSO trajectory willbeusedforeachcasetoshowtheimprovementbetweeneachcase. 30 In this example problem the initial velocity is found as outlined earlier in the chapter. This in turn will enable the MPSO algorithm to determine the optimal rendezvous trajectory upon the finalapproachsuchthatthefinalapproachisaprescribeddistance,d,away. TheMPSOprescribed parametersaregiveninTable3.1 Table3.1: Exampleproblem1MPSOparametersfortrajectoryanalysis d 10km C C ,C S 2.05,2.05 N (swarmsize) 24 J max (iterationlimit) 2000 S (numberofneighborhoods) 3 CR 0.9 t [0;1*Period] 31 Chapter4 Results The two-impulse MPSO solution is treated as base case, then in the subsequent sections, further optimization is done via adding an initial coast and/or an intermediate burn The basis of adding the coast and/or intermediate burn is to determine if any significant decreases to the cost function through means of primer vector theory (PVT). The results of the example problem are summarized in the following sections. Note that magnitudes smaller than1:010 *15 are taken to be zero. The results were achieved using Matlab 2017a, Intel i7-67000 HQ processor at 2.60 GHz and16.0GBofRAM.Theauthordoesnotbelievethecodeusedwasscriptedinanoptimalfashion, however the code was consistent and therefore the correlations between the different methods are believedtobevalid. 4.1 Two-ImpulseTrajectory The transfer trajectory was first analyzed utilizing a two-impulse transfer. The optimal tra- jectory determined by the modified particle swarm optimization (MPSO) solution was first ex- amined, followed by the same MPSO optimal two-impulse trajectory with an initial coasting arc determinedbyPVT. 32 4.1.1 ModifiedParticleSwarmSolution Thetwo-impulsePSOtrajectorywasfoundusingtheprocedureoutlinedinsection3.3. This servesasthebaselinecasefortheremainderofthework. Theoptimalpositionandvelocityvalues forthetwo-impulseMPSOtrajectoryaregiveninTable4.1. Table4.1: PositionandVelocityValuesforOptimalTwo-ImpulseMPSOSolution v * i v + i 1 =(-7.1169,-9999.9,5.7303)m (-13.661, -56.101, 6.6207) cm/s (-13.661, -56.101, 6.6207) cm/s 2 =(0,0,0)m (-2.8215, -57.656, 6.6501) cm/s (0, 0 , 0) cm/s Theimpulsevectors,magnitudes,andtimesofimpulsefortheoptimalMPSOtwo-impulsetrajec- toryaregiveninTable. 4.2. Table4.2:V Valuesand TimesforOptimalTwo-ImpulseMPSOSolution t(s) v i (cm/s) ñv i ñ (cm/s) v 1 0 .6:3205;25:692;*2:9774/ E-11 2.6624E-10 v 2 5424 .2:8215;57:656;*6:6501/ 58.107 total 5424 N/A 58.107 33 Theoptimaltwo-impulsetrajectoryfoundbytheMPSOalgorithmisseeninFig. 4.1 Figure4.1: Optimaltwo-impulsetrajectoryfoundviaMPSO 4.1.2 PrimerVectorTheorySolution The optimal two-impulse MPSO trajectory was then further optimized by the means of an initial coasting arc determined via primer vector theory. The optimal position and velocity values forthetwo-impulseMPSOtrajectorywithacoastingarcaregiveninTable4.3. 34 Table 4.3: Position and Velocity Values for Optimal Two-Impulse MPSO Solution with Coast AddedviaPVT v * i v + i 1 =(-1986.1,-6743.1,30.125)m (-54.597, 376.18, -5.7793) cm/s (-54.597, 376.18, -5.7793) cm/s 2 =(0,0,0)m (-2.8214, -57.656, 6.6501) cm/s (0, 0, 0) cm/s Theimpulsevectors,magnitudes,andtimesofimpulsefortheoptimalMPSOtwo-impulsetrajec- torywitha coastingarcaregiveninTable. 4.4. Table4.4:V ValuesandTimesforOptimalTwo-ImpulseMPSOSolutionwithCoastAddedvia PVT t(s) v i (cm/s) ñv i ñ (cm/s) v 1 2316 .895:12;117:42;*5:6649/ E-12 9.0281E-10 v 2 5424 .*2:8214;57:656;*6:6501/ 58.107 total 5424 N/A 58.107 35 The optimal two-impulse trajectory found by the MPSO algorithm with the addition of a coastingarcisseeninFig. 4.2 Figure4.2: Optimaltwo-impulsetrajectoryviaMPSOwithcoastaddedviaPVT It can be seen that while the coasting arc is compared to the entire transfer duration, the magnitude ofV 1 of the MPSO solution is on the order of approximately1:010 *9 cm/s and thereforetheburnisalmostnegligible. 36 4.2 Three-ImpulseTrajectory Theoptimaltwo-impulseMPSOtrajectorywasthenoptimizedviaanimpulsiveburn. This newthree-impulseMPSOtrajectorywasthenfurtheroptimizedbymeansofaninitialcoastdeter- mined by primer vector theory. Finally, these two trajectories were compared to the optimal three impulsetrajectorydeterminedbyprimervectortheory. 4.2.1 ModifiedParticleSwarmSolution The original optimal two-impulse MPSO trajectory was then improved through the use of an intermediate impulse. The intermediate impulse location, direction and magnitude were all determinedviatheMPSOalgorithm. Theoptimalpositionanvelocityvaluesforthethree-impulse MPSOtrajectoryaregiveninTable4.5. Table4.5: PositionandVelocityValuesforOptimalThree-ImpulseMPSOSolution v * i v + i 1 =(-7.1169, -9999.9, 5.7303) m (-13.661,-56.101,6.6207)cm/s (-13.661,-56.101,6.6207)cm/s m =(-1435.7, -8838.4, 57.334) m (-107.63,255.94,-2.2394)cm/s (-107.63,255.94,-2.2394) cm/s 2 = (0, 0, 0) m (-2.8212,-57.655,6.6501)cm/s (0,0,0)cm/s The impulse vectors, magnitudes, and times of impulse for the optimal MPSO three-impulse tra- jectoryaregiveninTable. 4.6. Table4.6:V ValuesandTimesforOptimalThree-ImpulseMPSOSolution t(s) v i (cm/s) ñv i ñ(cm/s) v 1 0 .2:8126;25:801;*3:0364/ E-11 2.6131E-10 v m 1666.4 .157:65;*142:11;6:0021/ E-14 2.1233E-12 v 2 5424 .2:8212;57:655;*6:6501/ 58.106 total 5424 N/A 58.106 37 Theoptimalthree-impulsetrajectoryfoundbytheMPSOalgorithmisseeninFig. 4.3. Figure4.3: Optimalthree-impulsetrajectoryviaMPSO 4.2.2 ModifiedParticleSwarmSolutionWithCoast Theoptimalthree-impulseMPSOtrajectorywasthenfurtheroptimizedbythemeansofan initial coasting arc determined via primer vector theory. The optimal position an velocity values forthethree-impulseMPSOtrajectorywithacoastingarcaregiveninTable4.7. 38 Table 4.7: Position and Velocity Values for Optimal Three-Impulse MPSO Solution with Coast AddedviaPVT v * i v + i 1 =(-36.720, -10071, 14.541) m (-30.278,-49.635,6.4577)cm/s (-30.278,-49.635,6.4577) cm/s m =(-1435.7, -8838.4, 57.334) m (-107.627,255.942,-2.2394)cm/s (-107.627,255.942,-2.2394) cm/s 2 = (0, 0, 0) m (-2.8212,-57.655,6.6501)cm/s (0,0,0)cm/s Theimpulsevectors,magnitudes,andtimeofimpulsefortheoptimalMPSOthree-impulsetrajec- torywitha coastingarcaregiveninTable. 4.8. Table 4.8: V Values and Times for Optimal Three-Impulse MPSO Solution with Coast Added viaPVT t(s) v i (cm/s) ñv i ñ(cm/s) v 1 134.602 .*8:9789;*53:669;*2:9579/ E-11 5.4496E-10 v m 1666.4 .1750:9;618:31;4:5311/ E-10 1.8568E-9 v 2 5424 .2:8212;57:655;*6:6501/ 58.106 total 5424 N/A 58.106 39 Theoptimalthree-impulsetrajectoryfoundbytheMPSOalgorithmwiththeadditionofacoasting arcisseeninFig. 4.3. Figure4.4: Optimalthree-impulsetrajectoryviaMPSOwithcoastaddedviaPVT Note that the results between the three-impulse MPSO trajectory and the three-impulse MPSO trajectory with a coast result in nearly the same cost. The added coast added a negligible differenceofapproximately1E-10cm/s. Thisresultcanalsobeseenbycomparingthemagnitudes ofthethreeimpulses. 40 4.2.3 PrimerVectorTheorySolution The optimal three-impulse solution found using primer vector, with the addition of a coast was computed to juxtapose the trajectories found using the MPSO algorithm. The position and velocityvaluesfortheoptimalthree-impulsePVTsolutionaregiveninTable. 49. Table4.9: PositionandVelocityValuesforOptimalThree-ImpulsePVTSolution v * i v + i 1 =(-14.509, -10025, 8.6949) m (-19.259,-54.487,6.5820)cm/s (-19.259,-54.487,6.5820) cm/s m =(-1281.7, -752.81, -59.144) m (112.68,222.31,-1.5810)cm/s (112.68,222.31,-1.5810)cm/s 2 = (0, 0, 0) m (-2.8212,-57.655,6.6501)cm/s (0,0,0)cm/s Theimpulsevectors,magnitudes,andtimeofimpulsefortheoptimalPVTthree-impulsetrajectory withacoastingarcaregiveninTable. 4.10. Table4.10:V ValuesandTimesforOptimalThree-ImpulsePVTSolution t(s) v i (cm/s) ñv i ñ(cm/s) v 1 44.901 .1337:6;*2772:2;1:0060/ E-9 1.3660E-6 v m 4008.6 .27399;*9956:3;1:5078/ E-10 2.9152E-6 v 2 5424 .2:8212;57:655;*6:6501/ 58.106 total 5424 N/A 58.106 41 The optimal three-impulse trajectory found by PVT with the addition of a coasting arc is seen in Fig. 4.. Figure4.5: Optimalthree-impulsetrajectoryviaPVT 42 For clarity, a magnified view of the coasting arc of the optimal three-impulse trajectory found by PVTisseeninFig. 5 Figure4.6: Optimalthree-impulsetrajectoryviaPVTmagnifiedview 4.3 GeneralResults Theexampleproblemwasrunatotaloffiftytimesforeachmethodtodeterminetheconver- genceaccuracyoveralargenumberofcases. TheBPSO,SPSO,andMPSOmethodswerethenall comparedtotheoptimaltwo-impulsetrajectorywithanaddedcoastviaprimervectortheory. This 43 analysiswasdonefortwocases,asmallswarmandnumberofgenerationsandalargerswarmwith an increased number of generations. This was done to compare the accuracy and efficiency when fewer particles and generations are used, a likely case when computational resources are limited. Due to limitations of computational resources, the general results were not carried out for a large number of iterations or re-runs. The author believes that similar results would be achieved based on the two-impulse results and the example problem results. however the results achieved on the example problem and the general results for the two-impulse trajectory . The average difference between the two-impulse PVT solution and each PSO algorithm, and the percent difference are giveninTable. 4.11. Table 4.11: Average Difference Between Two-Impulse PVT solution and Particle Swarm, swarm size=150,numberofgenerations=1000 AbsoluteV difference(cm/s) PercentDifference AverageDifferenceBetween PVTandBPSO -0.94807 -1.6437% AverageDifferenceBetween PVTandSPSO -0.63745 -1.1052% AverageDifferenceBetween PVTandMPSO -0.31165 -0.54033% The mean, minimum, and standard deviation of the maneuver cost for each particle swarm algo- rithmaregiveninTable. 4.12. Table4.12: Mean,Minimum,andStandardDeviationofTwo-ImpulseTrajectoryCostsAchieved byEachPSOMethod,swarmsize=150,numberofgenerations=1000 MeanCost(cm/s) MinimumCost(cm/s) StandardDeviation(cm/s) BPSO 58.626 58.166 1.6867 SPSO 58.315 57.763 0.89474 MPSO 57.989 57.663 0.23566 The two-impulse case was then re-run using a larger swarm size and number of iterations. The average difference between the three-impulse PVT solution and each PSO algorithm, and the 44 percentdifferencearegiveninTable. 4.13. Table 4.13: Average Difference Between Two-Impulse PVT solution and Particle Swarm, swarm size=200,numberofgenerations=3000 AbsoluteV difference(cm/s) PercentDifference AverageDifferenceBetween PVTandBPSO -0.15768 -0.47437% AverageDifferenceBetween PVTandSPSO -0.08765 -0.12057% AverageDifferenceBetween PVTandMPSO -0.01165 -0.06571% The mean, minimum, and standard deviation of the maneuver cost for each particle swarm algo- rithmaregiveninTable. 4.14. Table 4.14: Mean and Minimum Two-Impulse Trajectory Costs Achieved by Each PSO Method, swarmsize=200,numberofgenerations=3000 MeanCost(cm/s) MinimumCost(cm/s) StandardDeviation(cm/s) BPSO 58.105 57.766 0.24692 SPSO 58.011 57.663 0.22474 MPSO 57.889 57.601 0.13496 Itshouldbenotedthattheinitialstatedeterminedbythetwo-impulseMPSOmethodology mayvaryslightlydependingonthenumberofnumberofparticlesintheswarmandthenumberof generations used. This in turn, is the cause of any discrepancy between minimum value achieved in the general results and the optimal value found in the example problem. It should be noted that theexampleproblemresultsdocorrelatewellwiththegeneralresultsfound. 45 Chapter5 Summary 5.1 Conclusion Themodifiedparticleswarmoptimization(MPSO)algorithmsuccessfullydemonstratesthe ability in finding near optimal, or even optimal, rendezvous trajectories while utilizing a smaller swarm and fewer generations compared to previous work using particle swarm optimization. The near optimality of these trajectories are evident in the fact that the the results obtained by the two- impulse MPSO solution were only marginally improved by adding an initial coast and/or an in- termediateimpulseusingprimervectortheory(PVT).Theincreasedcomputationalefficiencyand reliability is seen by the fact that the MPSO algorithm achieved a lower transfer cost, a smaller mean transfer cost, and a smaller standard deviation compared to the other particle swarm opti- mization(PSO)algorithmswhenusingasmallerswarmsizeandnumberofgenerations. Thiscan becomeparticularlyusefulwhenalargersolutionspaceissearched. However,allPSOalgorithms successfullylocatednearoptimaltrajectorieswhentheswarmsizeandnumberofgenerationswere increased. Oneobservationmadeisthatthemagnitudeofthein-trackpositionandvelocitycompo- nentsdominatetheradialandcross-trackvalues. Thisshouldbeexpectedifonelookstooptimize thefuelandutilizethedeputy’sinitial motionandmomentumtoitsadvantage. In the first example problem the MPSO algorithm essentially determines the approach tra- 46 jectorysuchthattheapproachstatevectorisoptimalsothatessentiallyalloftherequiredv total is required at the terminal location to create the zero state vector. The second example is somewhat more interesting as it adds a constraint upon the incoming velocity and therefore the deputy can no longer simply coast or apply a negligible burn at the beginning. This leads to more significant improvementbyaddinganadditionalintermediateimpulse. The computational resources required to identify an optimal four-impulse trajectory are exhaustive, costing magnitudes longer in computational time and offer a negligible improvement (less than 0.5%) in most problems. Therefore, an optimized three-impulse trajectory with a coast, determined via PVT), is considered sufficient in many cases. Moreover, in many cases utilizing theMPSOalgorithmresultsinatwo-impulsetrajectorythatisnearlyasfuel-efficientasthethree- impulsetrajectory. TherearesomelimitationsoftheMPSOalgorithm,asitrequiresasufficientlylargenumber of particles in order to locate a near optimal solution. Through empirical testing, this number is estimated to be approximately 24, however results may vary for different problems and scales of solutionspaces. Similarly,thenumberofgenerationsmustbesufficientlylargetoenableathorough localsearchoftheunifiedparticleswarmalgorithm. 5.2 FutureWork There are many potential avenues of future research in the field of trajectory optimization, specifically rendezvous proximity operations. While this work focuses on the use of the MPSO algorithm for rendezvous operations in circular orbits, one potential avenue of future work is to investigate the MPSO’s ability to identify fuel-minimizing trajectories for rendezvous operations in elliptical orbits. Additionally, the incorporation of theJ 2 perturbation effects provides a natural nextstep. Theseperturbationeffectswouldhaveaconsiderableimpactfortrajectorieswithlonger transfertimes. Therearealsoopportunitieswithinthealgorithmthatmaybefurtherdeveloped,forexample 47 tuningofthecoefficients,C C andC S ,swarmsize,numberofgenerations,etc. asthiswouldfurther improve the algorithm’s efficiency. Additionally, burns with a finite-duration may be incorporated tohelpcreateamorerealisticmodel. 48 Bibliography [1] G. D. Bento da Silva and D. Barnhart, “Optimized r/bs rendezvous study for active de- bris removal”, in Proceedings of the Center for Obital Debris Education and Research (CODER2016),ser.LectureSeriesonRemediationTechnologiesandArchitectures, 2016. [2] D.Jezewski,“Primervectortheoryandapplications”,NASALyndonB.JohnsonSpaceCen- ter,Houston,Texas,Tech.Rep.NASATRR-454,Nov.1975. [3] B.Conway, SpacecraftTrajectoryOptimization.CambridgeUniversityPress,2010. [4] J.KennedyandR.Eberhart,inProceedingsoftheIEEEInternationalConferenceonNeural Networks,1995,pp.1942–1948. [5] K.E.ParsopoulosandM.N.Vrahatis,“UPSO:Aunifiedparticleswarmoptimizationscheme”, in Proceedings of the International Conference of Computational Methods in Sciences and Engineering (ICCMSE ‘04), ser. Lecture Series on Computer and Computational Sciences, vol.1,Zeist,TheNetherlands: VSPInternationalSciencePublishers,2004,pp.868–873. [6] B.Aubin,“Optimizationofrelativeorbittransfersviaparticleswarmandprimervectorthe- ory”,MSThesis,UniversityofIllinoisatUrbana-Champaign,2011. [7] R.Sherill,“Dynamicsandcontrolofsatelliterelativemotioninellipticorbitsusinglyapunov- floquettheory”,PhDDissertation,AuburnUniversity,2012. [8] P. Lion and M. Handelsman, “Primer vector on fixed-time impulsive trajectories”, AIAA Journal,vol.6,no.1,pp.127–132,1968. 49 [9] S. P. Hughes, L. M. Mailhe, and J. J. Guzman, “A comparison of trajectory optimization methodsfortheimpulsiveminimumfuelrendezvousproblem”,NASAGoddardSpaceFlight Center,Greenbelt,Maryland,Tech.Rep.,Feb.2003. [10] M. Clerc and J. Kennedy, “The particle swarm-explosion, stability and convergence in a multidimensionalcomplexspace”,IEEETransactionsonEvolutionaryComputation,vol.6, no.1,pp.58–73,2002. [11] K.E.ParsopoulosandM.N.Vrahatis,“Parameterselectionandadaptationinunifiedparticle swarmoptimization”,MathematicalandComputerModelling,vol.46,no.1-2,pp.198–213, 2007. [12] C. Zhang, J. Ning, S. Lu, D. Ouyang, and T. Ding, “A novel hybrid differential evolution andparticleswarmoptimizationalgorithmforunconstrainedoptimization”,OperationsRe- searchLetters,vol.37,no.2,pp.117–122,2009. 50
Abstract (if available)
Abstract
Rendezvous proximity operations are a crucial component of spacecraft missions, and will become even more so for Earth-based missions as the number of satellites and orbital debris increases. By the optimization of these rendezvous operations, a mission’s life may be extended substantially. In addition, more missions may become possible due to this higher level of efficiency. A modified particle swarm algorithm is proposed to identify fuel-efficient trajectories more quickly, with a higher level of accuracy, and with a higher level of reliability or stability. The Clohessy-Wiltshire equations are used to describe the relative motion between the two spacecraft. The results of the modified particle swarm algorithm are then improved, if possible, by primer vector theory (PVT) to ensure optimality of the maneuver. ❧ The proposed modified particle swarm optimization (MPSO) algorithm is proposed as a more sophisticated alternative to the basic particle swarm optimization (PSO) algorithm. The MPSO algorithm incorporates a component of mutation, based on a differential evolution (DE) algorithm, to help avoid stagnation and improve performance. The proposed MPSO algorithm successfully achieved near optimal trajectories when compared to those achieved via PVT. The primary benefits of the MPSO over its simplified PSO and DE counterparts is it’s improved ability in prevent stagnation in a multimodal solution space and enhanced local search, or exploitation, of the solution space. These benefits result in the MPSO algorithm locating optimal, or near optimal solutions with a higher level of efficiency and with a higher consistency. The higher level of efficiency is achieved by locating these near optimal solutions with a smaller swarm size and number of generations required when compared to both the PSO and DE algorithms, and the higher consistency is achieved by simply analyzing the mean and standard deviation of solutions over a large number of simulations.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Relative-motion trajectory generation and maintenance for multi-spacecraft swarms
PDF
Spacecraft trajectory optimization: multiple-impulse to time-optimal finite-burn conversion
PDF
Techniques for analysis and design of temporary capture and resonant motion in astrodynamics
PDF
Discrete geometric motion control of autonomous vehicles
PDF
Solution of inverse scattering problems via hybrid global and local optimization
PDF
I. Asynchronous optimization over weakly coupled renewal systems
PDF
Trajectory mission design and navigation for a space weather forecast
PDF
The next generation of power-system operations: modeling and optimization innovations to mitigate renewable uncertainty
PDF
Efficient bounded-suboptimal multi-agent path finding and motion planning via improvements to focal search
PDF
Numerical and experimental investigations of ionic electrospray thruster plume
PDF
Kinetic studies of collisionless mesothermal plasma flow dynamics
PDF
Optimizing element & system compliance of robotic, gecko adhesion-based grippers to unprepared space debris targets
PDF
Trajectory planning for manipulators performing complex tasks
PDF
Speeding up trajectory planning for autonomous robots operating in complex environments
PDF
Mathematical characterizations of microbial communities: analysis and implications
PDF
Artificial Decision Intelligence: integrating deep learning and combinatorial optimization
PDF
Modern nonconvex optimization: parametric extensions and stochastic programming
PDF
Robot trajectory generation and placement under motion constraints
PDF
Hierarchical planning in security games: a game theoretic approach to strategic, tactical and operational decision making
PDF
Automatic conversion from flip-flop to 3-phase latch-based designs
Asset Metadata
Creator
Danekas, Colter Rylan
(author)
Core Title
Optimization of relative motion rendezvous via modified particle swarm and primer vector theory
School
Viterbi School of Engineering
Degree
Master of Science
Degree Program
Astronautical Engineering
Publication Date
08/14/2018
Defense Date
08/14/2018
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
OAI-PMH Harvest,optimization,orbital mechanics,particle swarm,primer vector theory,relative motion
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Barnhart, David (
committee chair
), Giuliano, Paul (
committee member
), Lauda, Aaron (
committee member
)
Creator Email
colterdanekas@gmail.com,danekas@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-70679
Unique identifier
UC11672297
Identifier
etd-DanekasCol-6754.pdf (filename),usctheses-c89-70679 (legacy record id)
Legacy Identifier
etd-DanekasCol-6754.pdf
Dmrecord
70679
Document Type
Thesis
Format
application/pdf (imt)
Rights
Danekas, Colter Rylan
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
optimization
orbital mechanics
particle swarm
primer vector theory
relative motion