Collision Avoidance Maneuver Optimization with a Multiple-Impulse Convex Formulation
CC OLLISION A VOIDANCE M ANEUVER O PTIMIZATION WITH A M ULTIPLE -I MPULSE C ONVEX F ORMULATION
A P
REPRINT
Roberto Armellin ∗ Te P¯unaha ¯Atea – Auckland Space Institute,The University of Auckland20 Symonds Street, 1010 Auckland, New ZealandMarch 1, 2021 A BSTRACT
A method to compute optimal collision avoidance maneuvers for short-term encounters is presented.The maneuver is modeled as a multiple-impulse one to handle impulsive cases or to approximatefinite burn arcs associated either with short alert times or the use of low-thrust propulsion. Themaneuver design is formulated as a sequence of convex optimization problems solved in polynomialtime by state-of-the-art primal-dual interior-point algorithms. The proposed approach calculatesoptimal solutions without assumptions on thrust arc structure and thrust direction, with executiontime compatible with onboard use.
The number of catalogued Resident Space Object (RSO) is growing due to miniaturization of spacecraft (e.g., Cube-Sats) and the launch of mega-constellations (e.g., Starlink [1]). In parallel, the number of tracked space debris isexpected to increase due to the improvement in tracking systems (e.g., the space fence system [2]), resulting in a largernumber of conjunctions to be processed and more Collision Avoidance Maneuver (CAM)s designed and executed.Conjunction analysis and collision avoidance are currently performed by agencies and operators on the ground usingseveral tools and processes that were developed over the last twenty years [3]. These tools support operators’ activities;however, the decision process still requires human intervention. This approach will not be practical in the future whenconjunctions screening, collision avoidance decision process, and CAM design and execution need to be automated.Additionally, performing as many of these tasks onboard represents an enabling technology for safer and cheaperspace missions. Furthermore, more accurate conjunction services (i.e., limiting unnecessary maneuvers) and optimalCAMs will be required to reduce the propellant budget allocated for conjunction management. Within this context,this work aims to propose a method for CAM optimization that can be suitable for the use onboard of spacecraft witheither high- or low-thrust capabilities.A CAM is performed when, at the time of closest approach, a threshold on the miss distance, or the collision prob-ability, is exceeded [4]. In [5] a method to optimize an impulsive CAM based on the assumption of small maneuverwas introduced. This simplification allowed decoupling the problem into maneuver direction and magnitude determi-nation. The direction was determined by the gradient of the collision probability while the magnitude with an iterativeprocess. In [6] Alfano describes a tool for CAM analysis that can perform parametric studies of single-axis and dual-axis maneuvers. Collision probability contours for single-axis maneuvering are calculated based on an upper bound onthe impulse magnitude and a range of permissible maneuver times. By selecting a specific time, contours are producedfor dual-axis maneuvering. Agencies and operators use similar tools [3]. For example, a simple approach employingtangential maneuvers, thus sacrificing optimality, is adopted by the German Aerospace Center [7]. European SpaceAgency (ESA)’s tool CAMOS instead can deal with single and multiple impulses with arbitrary direction, differentobjective functions (e.g., collision probability, miss distance, total ∆ v ), bounds on the maneuvers, and constraints onthe post-maneuver orbital elements [8]. The drawback of this flexibility is that the obtained solutions are only locally ∗ email: [email protected] a r X i v : . [ m a t h . O C ] F e b PREPRINT - M
ARCH
1, 2021optimal, and therefore, the analyst must critically analyze the results. In [9] a multi-objective approach for CAMdesign was presented that enabled an exhaustive analysis of the problem building Pareto optimal solutions accordingto multiple criteria. This single-impulse approach also allows for merging station-keeping with CAM and assessingcollision risk for a window of one week after the maneuver. However, this approach is numerically intensive, as itrequires multiple evaluations of complex objective functions. Bombardelli and Hernando-Ayuso [10] developed ananalytical and semi-analytical method to find the fuel-optimal impulsive CAM. The proposed methods have provenconvergence as the problem is reduced either to an eigenvalue problem or a convex optimization one. However, as asingle impulse was considered, a single linearization was used, and the dynamics were Keplerian.The research on low-thrust debris avoidance is not very developed. A similar problem is studied in a great deal ofdetail for formation flying collision-free reconfiguration. Both direct and indirect optimal control approaches were alsoproposed and compatible with onboard use (see the introduction in [11] and the references therein for an exhaustiveoverview). Restricting the analysis to low-thrust CAM design, Reiter and Spencer [12] noted that the impulsiveapproximation could be inaccurate for short alert times; thus, finite burn analysis is needed. However, to cope with thisissue, a radial thrust assumption was made. In [13] an approach based on averaged dynamics and Gauss variationalequation is proposed with the underlying assumption of continuous tangential thrust. Four methods based on anindirect formulation of an optimal control problem are presented in [14]. It is concluded that a semi-analytical methodbased on the linearization of the dynamics offers the best compromise between accuracy and computational time.However, formulating an energy optimal control problem without bounds on thrust magnitude is a significant limitationfor this approach. In this work, we present a methodology for optimal CAM design with convergence properties andexecution time compatible with autonomous onboard use. The approach is suitable for both impulsive and low-thrust CAM design and can handle short-term encounters with warning times from few minutes to multiple orbitalrevolutions. Furthermore, no a priori assumption on the direction of the impulses is made, and an arbitrary dynamicalmodel can be used. Constraints either or miss distance, maximum collision probability, or an approximated value ofthe collision probability can be enforced while minimizing the total ∆ v . The approach is based framing a multiple-impulse CAM optimization problem as a convex optimization one [15]. Thanks to the proven existence and uniquenessof the solution and computational advantages ensured by polynomial complexity, convex optimization has found manyapplications in aerospace engineering in the last 15 years [16], including long-duration low-thrust transfers [17], [18].The multiple-impulse CAM optimization is a Nonlinear Programming Problem (NLP) problem, and three steps arerequired to formulate it as a convex optimization problem, in our case as a Second-Order Cone Programming (SOCP)problem. Three main issues are encountered:1. the minimum fuel problem is non-convex;2. the dynamics and the transformation of coordinates are nonlinear;3. the constraints are non-convex.The introduction of slack variables and a lossless convexification [19],[18] solve the first issue. The dynamics and thetransformations required to perform the b-plane analysis are linearized, using Differential Algebra (DA) implementedin DACE as a first-order automatic differentiation technique [20]. Due to the small deviations introduced by CAMs,linearized dynamics accurately describe the conjunction dynamics [10, 14]. As a result, the problems associated withlinearization (i.e., artificial infeasibility and unboundedness [21]) are not relevant here, and the successive convexi-fication is introduced for refinement purposes, with no need for virtual controls and complex trust-region strategies[22].The third issue is tackled by working on the squared Mahalanobis distance, whose contour lines describe an ellipse onthe b-plane. The squared Mahalanobis distance allows us to set a constraint either on minimum distance, maximumcollision probability, or an approximation of the collision probability [23]. However, the constraint on this quantityresults in a non-convex problem, as the admissible region on the b-plane is non-convex (as in the keep-out zone inrendezvous dynamics [24]). The projection and linearization technique proposed in [25] is adopted to tackle this issue,resulting in a second iterative procedure. Depending on the chosen initial guess, the iterations can converge to differentoptima. However, two suitably selected guesses allows for the identification of the global minimum.The final CAM optimization approach is solved with the state-of-the-art primal-dual interior-point algorithm imple-mented in the software MOSEK [26]. Solutions with hundreds of impulses are obtained robustly and efficiently,representing both high-thrust and low-thrust maneuvers. We test our algorithm on 2,170 real conjunctions derivedfrom the ESA Collision Avoidance Challenge https://kelvins.esa.int/collision-avoidance-challenge The paper is organized as follows. In Sec. 2 a brief overview of the short-term conjunction dynamics is providedtogether with a nonlinear approach to study the effect of a CAM on the conjunction geometry. Section 3 contains a The open-source software package is available at https://github.com/dacelib/dace PREPRINT - M
ARCH
1, 2021description of the methodology developed in this work. We first state the CAM design as NLP problem, followedby the description of the steps required for its convexification. The algorithm’s application to a large set of cases isdescribed in 4, leading to the conclusive section.
A brief introduction of the key quantities of a short-term encounter is provided. Only relevant concepts for the designof CAMs are summarized. Afterward, we introduce the DA-based approach to study the effect of maneuvers on theconjunction geometry used to define the optimization problem of Sec. 3.
We consider the conjunction between a primary (subscript p ) and a secondary (subscript s ). The primary is thespacecraft we control, whereas the secondary is assumed to be passive. We indicate the relative position and velocityvectors at the Closest Approach (CA) as ∆ r ∗ CA = r ∗ p,CA − r ∗ s,CA , (1) ∆ v ∗ CA = v ∗ p,CA − v ∗ s,CA (2)in which r and v are the absolute position and velocity vectors, and the asterisk indicates the ballistic motion. As atthe CA the distance between the two objects is minimum, it follows ∆ r ∗ CA · ∆ v ∗ CA = 0 . (3)To compute the collision probability it is useful to introduce a coordinate system referred to as the b-plane. The originof the axes of this frame lies at the centre of the secondary object at the time of CA; the η -axis is defined along thedirection of the relative velocity of the primary with respect to the secondary object; the ξζ plane is perpendicular tothat η -axis ˆ u ξ = v ∗ s,CA × v ∗ p,CA (cid:107) v ∗ s,CA × v ∗ p,CA (cid:107) , (4) ˆ u η = v ∗ p,CA − v ∗ s,CA (cid:107) v ∗ p,CA − v ∗ s,CA (cid:107) , (5) ˆ u ζ = ˆ u ξ × ˆ u η . (6)All the introduced quantities are shown in Fig. 1.Figure 1: Conjunction geometry and b-plane definition3 PREPRINT - M
ARCH
1, 2021The unit vectors define the rotation matrix from the inertial reference frame to the b-plane R D = [ˆ u ξ ˆ u η ˆ u ζ ] T , (7)while the projection in the η -axis is achieved by R D = [ˆ u ξ ˆ u ζ ] T . (8)The nominal position of the primary on the b-plane at the time of closest approach t ∗ CA is ∆ r ∗ CA = ( ξ ∗ , ζ ∗ ) . Theparticular case in which ∆ r ∗ CA = 0 is referred to as a direct impact.Under the short encounter approximation and Gaussian distribution of objects state vectors, the collision probability is P ∗ C = 12 π (det C ∗ CA ) / (cid:90) (cid:90) A e − ( ∆ r − ∆ r ∗ CA ) T ( C ∗ CA ) − ( ∆ r − ∆ r ∗ CA ) d ξ d ζ (9)in which C ∗ CA is the sum of the positional covariance of the two objects referred to a common reference frame andprojected onto the b-plane via (8), and A collision cross sectional area, a circle of radius R = R p + R s where R p/s is the radius of the sphere enclosing the primary/secondary respectively. Several methods have been developed overthe years to calculate P C . Among them, [27], [28], and [29] are worth of a particular merit as they provided analyticalsolutions. As explained in [23], the simplest approach consists in approximating the value of P C by assuming theprobability density is constant over the collision circle. The approximate value is P ∗ C = R C ∗ CA ) / e − ( d ∗ CA ) , (10)in which d ∗ CA = (cid:113) ∆ r ∗ CAT ( C ∗ CA ) − ∆ r ∗ CA is the Mahalanobis distance. Additionally, the maximum collision prob-ability can be obtained by optimally scaling the combined covariance ([23]) resulting in P ∗ C, max = R (cid:0) d ∗ CA (cid:1) (det C ∗ CA ) / e . (11)These approximations allow us to use the squared Mahalanobis distance to set constraints on an approximate valueof the collision probability through Eq. (10), the maximum collision probability through Eq. (11), or the miss dis-tance (by setting C ∗ CA = I , the identity matrix). Moreover, for constant covariance matrix, the contour lines of thesquared Mahalanobis distance describe ellipses on the b-plane, a type of constraint that can be dealt with efficientlyby successive convexifications, as it will be shown later. When the primary is maneuvered, the conjunction geometry changes. This aspect is here analyzed with the use ofarbitrary order Taylor expansions enabled by DA (the notation adopted in [30] is used). These effects are generallysmall for short-term encounters and small maneuvers. Nevertheless, a general treatment is provided as the proposedapproach could be potentially applied to longer encounters with minimal changes (provided that the distance functionis convex in time).The first step is to use DA to introduce perturbations to the closest encounter time t ∗ CA + δt , the primary position δ r ∗ p,CA + δ r p , and the primary velocity δ v ∗ p,CA + δ v p . Taylor expansions of the states of both the primary andsecondary are obtained by DA-based numerical integrations (expansion in time and state, see [31] for details) of theorbital dynamics, delivering r p = T r p ( δt, δ r p , δ v p ) v p = T v p ( δt, δ r p , δ v p ) r s = T r s ( δt ) v s = T v s ( δt ) . (12)Note that in Eq. (12), due to the effect of the introduced perturbations, we have dropped the CA subscript. Additionally,the state of the secondary is only affected by the time perturbation.From Eq. (12) the relative quantities can be calculated ∆ r = T ∆ r ( δt, δ r p , δ v p ) ∆ v = T ∆ v ( δt, δ r p , δ v p ) . (13)4 PREPRINT - M
ARCH
1, 2021The closest encounter conditions, Eq. (3), in DA formalism, reads ∆ r · ∆ v = T ∆ r · ∆ v ( δt, δ r p , δ v p ) = 0 . (14)This constraint is a parametric implicit equation that can be solved for δt by polynomial partial inversion techniques(see [31] and [32] for more details). The polynomial partial inversion provides δt = T δt ( ∆ r · ∆ v , δ r p , δ v p ) (15)and, by substitution of the closest encounter constraint ∆ r · ∆ v = 0 , we obtain δt CA = T δt C A ( δ r p , δ v p ) (16)This map gives a Taylor approximation of how the the closest encounter time changes due to a variation of the state ofthe primary (due to a maneuver), an aspect that is commonly ignored in CAM design. This polynomial can be insertedback in Eq. (12), obtaining r p,CA = T r p,CA ( δ r p , δ v p ) v p,CA = T v p,CA ( δ r p , δ v p ) r s,CA = T r s,CA ( δ r p , δ v p ) v s,CA = T v s,CA ( δ r p , δ v p ) . (17)These polynomial maps approximate the states of both objects at the different times of CA as a result of a CAM thatproduces a variation of position and velocity of the primary object at t ∗ CA . Note that in Eq. (17) we re-introduce theCA subscript, but we remove the asterisk, as now each perturbed solution has a different time of closest approach,determined by Eq. (16). Similarly, the projection matrix in Eq. (8) is expanded by using Eq. (17) for the calculationof the unit vectors, resulting in R D = T R D ( δ r p , δ v p ) (18)Equation (18) is used both to project the relative state and the combined covariance matrix on the b-plane, thus allowingfor the calculation of the expansion of all the relevant conjunction quantities: ∆ r CA = T ∆ r CA ( δ r p , δ v p ) , (19) d CA = T d CA ( δ r p , δ v p ) , (20) P C = T P C ( δ r p , δ v p ) , (21)and P C, max = T P C, max ( δ r p , δ v p ) . (22)Note that in Eqs. (19)-(22) the dependency on δ r p and δ v p is due to the covariance C CA = T C CA ( δ r p , δ v p ) as aresult of the dependency of the projection matrix (18) on the perturbation of the states at the conjunctions. On theother hand, it is always assumed that the state covariances provided in the Conjunction Data Message (CDM) are notdirectly altered by the implementation of the maneuver. The details of CAM design algorithm for multiple-impulse maneuvers is presented. Before framing it as a successiveconvexification problem, a general NLP formulation is described to provide the general setting on which the successiveconvexification is introduced.
A uniform N -point time grid is constructed by selecting a discretization time step ∆ t and starting from the earliestmaneuvering time t . Note that more refined discretization schemes are possible without adding complexity to the al-gorithm (e.g., by using a suitable angular variable to better deal with eccentric orbits), but these were not needed here asthe orbits in our test cases are almost circular. The number N is selected as N = min (cid:18) floor (cid:18) t CA ∗ − t ∆ t (cid:19) , N max (cid:19) ,where N max accounts for the maximum time span in which a maneuver can be implemented. At every discretizationpoint a maneuver can be added in the form of an instantaneous change in velocity ∆ v i with i = 0 , . . . , N − . Thenominal trajectory at discretization points is given by r ∗ i , v ∗ i with i = 0 , · · · , N , with r ∗ N = r ∗ CA , v ∗ N = v ∗ CA . For5 PREPRINT - M
ARCH
1, 2021each i = 0 , . . . , N − we calculate a o -th order Taylor approximation of the mapping between deviations in the initialstate and deviations of the final state δr − i +1 = T δr − i +1 ( δr + i , δv + i ) δv − i +1 = T δr − i +1 ( δr + i , δv + i ) . (23)In Eq. (23) the + superscript indicates chosen quantities at the beginning of an interval, whereas the − propagatedquantities from the previous interval. These maps are built with N DA integrations of a dynamical model of choiceusing the unperturbed trajectory as reference (e.g., obtained by backward propagation from the CDM with the samedynamical model). As the CAM ∆ V s are small, a low order (in most cases order 2) allows for sufficiently accurateapproximations without the need of iterations.The optimization variables are the set of N ∆ v i applied at nodes i = 0 , . . . , N − . The effect of these maneuvers onthe trajectory of the primary are obtained by the use of maps of Eq. (23). In particular, for i = 0 we can set δr +0 = 0 and δv +0 = ∆ v , and, by applying (23), we obtain the mapped perturbations δr − and δv − . For i = 1 , . . . , N − we proceed by defining the new perturbations δr + i = δr − i − and δv + i = δv − i − + ∆ v i and use the i -th set ofmaps of Eq. (23) to map these perturbations to the end of the segment. At i = N − we achieve how the set x = [ ∆ v ; . . . ; ∆ v N − ] is mapped into the final perturbations δr p and δv p , which then in turn allows us to computethe relevant quantities by Eqs. (19)-(22). The multiple-impulse CAM optimization problem, can be stated as follows:Minimize the sum of magnitudes of the impulses min x N − (cid:88) i =0 ∆ v i (24)subject to nonlinear inequality constraints ∆ v i ≤ ∆¯ v for i = 0 , . . . , N − ,P C ≤ ¯ P C or P C, max ≤ ¯ P C, max or ∆ r CA ≥ ¯ d min , (25)where the overline indicates assigned values, and ∆ v i the impulses magnitude. This optimization problem is a NLPthat can be solved with dedicated solvers. In this work we use the Sequential Quadratic Programming (SQP) algorithmimplemented in the MATLAB fmincon function providing analytical gradient of the objective function and Jacobianof the constraints. For the latter the derivatives included in the polynomials maps are used together with the chainrule to calculate the sensitivity of the constraints with respect to the optimisation vector. The summary of the NLPformulation is provided in Algorithm 1. As NLP problems, this formulation is non-deterministic polynomial-time hard(NP-hard), meaning that the computation time may be very long if the problem is solved at all [24]. Algorithm 1
Nonlinear programming formulation Get inputs from CDM: R , r ∗ p/s,CA v ∗ p/s,CA , t ∗ CA , C ∗ p/s,CA ; Assign t , ∆ t , ∆¯ v , N , ¯ P C or ¯ P C, max or ¯ d min , and x ; Back propagate the trajectories from t ∗ CA to t and save r p/s,t , and v p/s,t ; Define the time grid ( t : ∆ t : t + N ∆ t ) ; Build maps Eq. (23) by
N o -th order DA forward propagations; Solve the NLP problem defined by (24) and (25);
As described in the introduction, three main steps are required to formulate the CAM design problem as a convexoptimization one. Firstly, in Sec. 3.2.1 the objective function and the constraints on the velocity magnitude are re-formulated by introducing slack variables and lossless convexification. Afterward, the dynamics of the problem arelinearized in Sec. 3.2.2, followed by constraints linearization in Sec. 3.2.3. Only the constraint on the squared Maha-lanobis distance is taken into account as this type of constraint is effectively handled by a projection and linearizationapproach [25]. The details of the resolution algorithm are then presented in Sec. 3.2.4.
The ∆ v magnitudes are introduced as slack variables in the optimization problem. As a result, each impulseis described by four independent variables ∆˜ v i = [ ∆ v i ; ∆ v i ] , and the optimization vector becomes x = PREPRINT - M
ARCH
1, 2021 [ ∆ v ; . . . ; ∆ v N − ; ∆ v ; . . . ; ∆ v N − ] . The introduction of the slack variables renders objective function linear min x N − (cid:88) i =0 ∆ v i (26)and transforms the N constraints on the impulses magnitude in Second Order Cone (SOC) constraints (cid:113) ∆ v i,x + ∆ v i,y + ∆ v i,z ≤ ∆ v i for i = 0 , . . . , N − . (27)Lastly, the bounds on the slack variables ≤ ∆ v i ≤ ∆¯ v for i = 0 , . . . , N − (28)are added to the problem. This convexification step is referred to lossless as it can be proved that the optimal solutionof the covexified problem is also the optimal solution of the original one [24]. The introduction of the slack variables is not sufficient to make the problem convex due to the nonlinearities in Eq.(23). These are dealt with by successive linearizations, requiring an iterative process. We will refer to these iterationsas major iterations associated with index j in the remainder of this section.Assume a solution x j − is available providing a reference trajectory (only at the first major iteration the referencetrajectory is ballistic) about which the dynamics are linearized. At the j -th iteration, the linear part of Eq. (23) can beextracted resulting in (cid:20) δr − i +1 δv − i +1 (cid:21) j = [ A δr ,i A δv ,i ] j (cid:20) δr + i δv + i (cid:21) j = A ji (cid:20) δr + i δv + i (cid:21) j (29)The composition of all these linear maps results in (cid:20) δr N δv N (cid:21) j = (cid:20) δr p δv p (cid:21) j = [ A N − A N − . . . A δv , , A N − A N − . . . A δv , , . . . , A δv ,N − , × N ] j (cid:0) x j − x j − (cid:1) == A j (cid:0) x j − x j − (cid:1) , (30)where ( x j − x j − ) is due to the fact that the optimization vector x j contains absolute impulses, whereas the lin-earizations are about the optimal impulses of the ( j − -th iteration. The perturbed time of closest approach, t jCA ,is t jCA = t j − CA + B j A j (cid:0) x j − x j − (cid:1) = B j A j x j + (cid:16) t jCA − B j A j x j − (cid:17) , (31)where B is the the × linear part of Eq. (16). Similarly, the perturbed relative position vector on the b-plane, ∆ r jCA ,is given by ∆ r jCA = ∆ r j − CA + C j A j (cid:0) x j − x j − (cid:1) = C j A j x j + (cid:16) ∆ r j − CA − C j A j x j − (cid:17) , (32)in which C j is the × linear part of Eq. (19) at the j -th iteration. When j = 1 , all the ( j − quantities inthe equations above are relative to the un-maneuvered case, i.e. ∆ r CA = ∆ r ∗ CA and x = . Note that as theCAMs determine a variation of the position vector of only few kilometers with respect to the ballistic trajectory (i.e.,a relative variation of . ), the convexification of the dynamics does not come with issues like artificial infeasibilityand unboundedness. For this reason we did not introduce any artificial control or sophisticated trust-region algorithms,but only problem-driven bounds on the impulses magnitude and a maximum final state deviation on the b-plane. The last step consists of dealing with the squared Mahalanobis distance constraint that defines an elliptically shapedavoidance region (or keep-out region using rendezvous terminology) on the b-plane. This non-convex constraint isdealt with by a second iterative process, nested in each major iteration, consisting of a projection and a linearization.We refer to these iterations as minor iterations , with index k .Assume the solution x j,k − of the j -th major iteration and ( k − -th minor iteration is available. This defines therelative position vector on the b-plane ∆ r j,k − CA . The k -th iteration starts by the projection algorithm that finds the7 PREPRINT - M
ARCH
1, 2021point z k on the ellipse (cid:0) d CA (cid:1) j = ¯ d CA closest to ∆ r j,k − CA . This is a convex optimization sub-problem with objectivefunction min z || ∆ r j,k − CA − z || (33)subject to the inequality constraint z T ( C jCA ) − z T ≤ ¯ d CA , (34)that can be solved efficiently using convex optimization algorithms (the interested reader can refer to [32] for a ex-haustive analysis of this sub-problem). C jCA does not depend on k as this quantity is assumed constant within eachminor iteration.Once z k is computed, the squared Mahalanobis distance constraint is linearized. The linear constraint ensures that ∆ r j,kCA belongs to the half-plane tangent to the constraint in z k , i.e. ∇ d M j ( z k )( ∆ r j,kCA − z k ) ≥ . (35)By substituting Eq. (32), the final expression is obtained − ∇ d M j ( z k ) C j A j x j,k ≤ ∇ d M j ( z k )( ∆ r j − CA − z k − C j A j x j − ) , (36)where the ( j − indicates known quantities at the end of the minor iteration of the ( j − -th major iteration.Additionally, ∆ r j, CA = ∆ r j − CA .The selection of the first starting point determines the convergence of the algorithm to a local optimum of the original,non-convex, problem. However, as it will be illustrated in Sec. 4.2, the problem appears to have only two two localminima on the opposite side of the elliptical boundary of the avoidance region. Thus, two initial starting pointsdetermined by ∆ r , CA = ± ∆ r ∗ CA are sufficient to automatically identify the global minimum. With the introduction of the slack variables, the linearization of the dynamics, and the linearization of the constraints,the multiple-impulse CAM design problem can be solved by primal-dual interior-point methods. We use MOSEK [26]through its MATLAB interface.The Algorithm 2 provides a complete overview of the solution process. As previously described, two iterations areneeded: a major j -th iteration for the linearization of the dynamics and a minor k -th iteration for linearization of theconstraint. The two iterations stops when two conditions are met: the minor one ends when || ∆ r j,kCA − ∆ r j,k − CA || ≤ tol m , while the major one when || x j − x j − || ∞ ≤ tol M . The algorithm minimizes the total ∆ v while constrainingeither the risk, the maximum risk, or the closest approach’s distance. The latter case is dealt by setting C ∗ CA = I inall iterations. The methodology described in the previous sections is applied to test cases derived from the ESA Collision AvoidanceChallenge. For this competition, ESA provided the teams with real conjunction data extracted from 162,634 CDM,corresponding to 13,154 unique events. However, ESA did not distribute the full orbital elements set and provided thepositional covariances in the primary Radial, Transverse, and Normal (RTN) reference frame only. With a procedureomitted here, we managed to reconstruct the full orbital data of the objects except for the absolute value of the rightascension of the ascending node of the primaries, which are arbitrarily set to . These data were filtered to considerconjunctions with ∆ r ∗ CA ≤ km, P ∗ C > − , and P ∗ C, max > − , resulting in a new data file with 2,170 conjunc-tions, available for download at github.com/arma1978/conjunction . The distribution of the minimum distanceat CA, the collision probability and maximum collision probability are shown in Fig.2-4. All the conjunctions arerelative to objects in Low Earth Orbit (LEO) with 90 % of cases with relative conjunction speed ∈ [1 . , . km/s.8 PREPRINT - M
ARCH
1, 2021
Algorithm 2
Successive convexification optimization algorithm Get inputs from CDM: R , r ∗ p/s,CA v ∗ p/s,CA , t ∗ CA , C ∗ p/s,CA ; Assign t , ∆ t , ∆¯ v i N , ¯ P C or ¯ P C, max or ¯ d min , tol M , tol m ; if ¯ P C, max is defined then Calculate C ∗ CA and ¯ d CA by Eq. (10); else if P C, max is defined then Calculate C ∗ CA and ¯ d CA by Eq. (11); else Set C ∗ CA ← I and ¯ d CA ← ¯ d end if Back propagate the trajectories from t ∗ CA to t and save r p/s,t , and v p/s,t ; Define the time grid ( t : ∆ t : t + N ∆ t ) ; j ← , x ← , t CA ← t ∗ CA , ∆ r , CA ← ± ∆ r ∗ CA , C CA ← C ∗ CA ; while ( j = 0) or || x j − x j − || ∞ ≥ tol M do j ← j + 1 ; Perform a -st order DA propagation of the trajectories at ( t : ∆ t : t + N ∆ t ) and t j − CA with impulses extractedfrom x j − ; if P C is defined then Calculate C jCA and then ¯ d CA from Eq. (10); else if P C, max is defined then Calculate C jCA and then ¯ d CA from Eq. (11); else Set C jCA ← I and ¯ d CA ← d ; end if Assemble the matrices A j , B j , C j as in Eq. (29)-(32); k ← and ∆ r j, CA ← ∆ r j − CA ; while ( k = 0) or || ∆ r j,kCA − ∆ r j,k − CA || ≥ tol m do k ← k + 1 ; Calculate z k by solving the convex optimization sub-problem (33)-(34); Calculate x j,k by solving the convex optimization problem defined by (26)-(28) and (36); t j,kCA ← B j A j x j,k + (cid:16) t j − CA − B j A j x j − (cid:17) ; ∆ r j,kCA ← C j A j x j,k + (cid:16) ∆ r j − CA − C j A j x j − (cid:17) ; end while x j ← x j,k , t jCA ← t j,kCA , and ∆ r jCA ← ∆ r j,kCA ; end while All the simulations presented in the next sections were obtained with a dynamical model including J − J zonalharmonics ˙ x = v x ˙ y = v y ˙ z = v z ˙ v x = − µxr + µJ R e r ( z r − x + µJ R e xz r ( z r −
3) + µJ R c x r (1 − z r + z r )˙ v y = − µyr + µJ R e r ( z r − y + µJ R e yz r ( z r −
3) + µJ R e y r (1 − z r + z r )˙ v z = − µzr + µJ R e r ( z r − z + µJ R e r ( − z r + z r ) + µJ R e z r (5 − z r + z r ) (37)in which r = [ x, y, z ] T and v = [ v x , v y , v z ] T are the spacecraft position and velocity vectors; µ , R e , and J i arethe gravitational parameter, the mean equatorial radius, and the i -th zonal harmonic coefficient of the Earth. Anydynamical model can be selected without affecting the algorithm complexity as the required state transition matricesare automatically obtained with DA without the need to derive and integrate the variational equations.In Sec. 4.2–4.5 we offer a detailed analysis for the reference scenario, which is the case with the highest collisionprobability in the dataset (details in Appendix). In Sec. 4.6 we provide a summary of the method performance whenapplied to all the datasets. Unless differently specified, all simulations assume ∆ t = 1 min, a maximum impulse of E − m/s (corresponding to a constant acceleration of E − m/s in each segment), tol M = 1 E − m/s, and9 PREPRINT - M
ARCH
1, 2021tol m = 1 m. The one-minute discretization is deemed appropriate to approximate finite burn arcs by a sequence ofimpulses. Note that a uniform discretization in an angular variable, such as the true anomaly, should be in generalpreferred unless orbits are almost circular as in our test cases. The simulations are run on a MacBook Pro with a 3,5GHz Dual-Core Intel Core i7 and 16 GB Memory.Figure 2: Distribution of minimum distance at closest approachFigure 3: Distribution of collision probabilities Figure 5 provides details of the convergence of the successive iteration algorithm summarised in Algorithm 2 for acase in which the CAM can be applied between 8 and 6 orbits before the CA with a maximum of 200 impulses and atarget ¯ P C, max = 1 E − . In Fig. 5(a) the contour line P C, max = ¯ P C, max is plotted as a black solid line, the dot closeto the origin is the unperturbed primary position with ∆ r CA = 43 m, while the circles represent the perturbed primaryposition during the algorithm iterations colored according to the maneuver δv . The optimization requires two majoriterations, with 5 and 1 minor iterations, respectively. In Fig. 5(b) the first minor iteration starts from the red circlelabeled representing z , which is the point on the constraint line closest to the unperturbed solution. The outputof the first run of the optimizer delivers the solution ∆ r , CA indicated as a yellow filled circle in the figure, whichbelongs to the ellipse tangent in z . From ∆ r , CA the new z is calculated allowing for the updated solution ∆ r , CA with significant reduction of maneuver ∆ v . Figure 5(c) shows similar information for the third and fourth iterationsand 5(d) for the fifth iteration of the first major loop and the first and only iteration of the second major iteration. Itis worth highlighting that the contour line of equal maximum collision probability is slightly changed (dashed line)because, between the first and second major iteration, the b-plane and the combined 2D covariance are updated.10 PREPRINT - M
ARCH
1, 2021Figure 4: Distribution of maximum collision probabilities (a) Convergence summary (b) First and second iterations(c) Third and fourth iterations (d) Fifth iteration and first iteration of second major iteration
Figure 5: B-plane convergence analysisFigure 6(a) highlights the thrust arcs on the primary trajectory (with the pentagram indicating the conjunction point)and Fig. 6(b) shows the details of the impulses in the RTN directions. The thrust arcs’ optimal location is spread along11
PREPRINT - M
ARCH
1, 2021the portion of the orbit opposite to the conjunction. The total ∆ v is . m/s, requiring 34 impulses. Figure 6(b)shows that the optimal maneuver has components in the radial and, to a lesser extent, in the out of plane directions. It isremarkable how each minor iteration of the problem is solved in approximately 12 ms and 14 iterations, allowing for acomplete solution in less than 0.1 s. This computational time is required to solve a large optimization problem entailing800 optimization variables, 200 second-order cone constraints, one inequality constraint on the squared Mahalanobisdistance, and simple bounds. By contrast, the NLP method described in Sec. 1, requires, when it reaches convergence,thousands of iterations and a few minutes to achieve a comparable solution starting from a random feasible guess. (a) Thrust arcs on primary trajectory (b) Impulses timing and components Figure 6: Example of trajectory and impulses profileAs previously discussed, it is not guaranteed that the minor iterations converge to the original non-convex problem’sglobal optimum. Convergence to a local minimum can occur depending on the selected starting point for the minoriteration. However, as the optimal solution lies on the admissible region’s elliptical boundary, we can investigate theobjective function’s behavior on this boundary. This analysis is done by sampling the ellipse with M points and solving M convex optimization problems in which a terminal equality constraint substitutes the keep-out zone one. Figure7 shows the results for M = 300 , where the pentagram indicates a local minimum (0.2139 m/s) and the hexagramthe global one (0.2042 m/s) analyzed in the previous figures. The examined cases have similar objective functionstructure, with the two minima located at the opposite side of the ellipse and corresponding to thrust mainly alignedwith either the tangential or the anti-tangential direction. These two minima can be identified by starting the minoriterations with two points on the ellipse with opposite coordinates.Figure 7: Objective function profile on the boundaries of the avoidance region12 PREPRINT - M
ARCH
1, 2021
The effect of the different constraints on the CAM design is analyzed, considering impulses on the last two orbitsbefore the conjunction. Figure 8(a) shows the contour lines corresponding to ¯ d min = 2 km, ¯ P C, max = 1 E − (the case of the previous section), and ¯ P C = 1 E − . The collision probability constraint is the less restrictive one,corresponding to a ∆ v = 0 . m/s and five impulses. The maximum collision probability requires ∆ v = 0 . m/s and 48 impulses, slightly higher than the previous section’s case due to the shorter time to the conjunction.However, these two constraints are characterized by the same elliptical shape on the b-plane with different sizes. Theminimum distance constraint is a circle, and the chosen value of km results in the most constraining one resultingin ∆ v = 0 . m/s and 88 impulses. Changing the constraint value makes it possible to identify a line of optimaltarget points on the b-plane. An example is provided in Fig. 8(b) where different values of ¯ P C, max are considered,resulting in ∆ v ∈ [0 . , . m/s. A trade-off between safety and propellant consumption is enabled by therapid calculation of CAM. (a) Different constraints and corresponding solutions indicatedwith different markers (b) Different ¯ P C, max contour levels and corresponding solu-tions indicated with different markers Figure 8: B-plane analysis according to different constraints
We analyze the effect of the alert time on the CAM. In Fig. 9 four graphs with leading time decreasing from 16 to 4orbits are reported. The ∆ v is 0.1281 m/s for the first case, 0.1534 m/s for the second, 0.2042 m/s for the third and0.2681 m/s for the last. The optimal maneuver in this case is performed at the first passage through the orbit oppositethe conjunction. When the leading time reduces and the ∆ v increases, the optimal maneuver is split into multiple arcs,two in this case. Additionally, the optimal maneuver always has a radial component. Furthermore, a maneuver appliedseveral orbits before the conjunctions has a larger impact on the conjunction geometry, affecting the time of CA formore than 2 seconds.Although the behaviour described in Fig. 9 is the most common one, not always the optimal maneuver occurs atearliest time. Figure 10 shows that, for test case number 10 in our dataset, the maneuver is performed at the latestopposition with respect to the conjunction. This counter-intuitive behaviour is probably due to the characteristics ofthe conjunction in terms of geometry and covariances, but its explanation requires further work. The effect of varying the single impulse magnitude is studied. In Fig.11 the maximum impulse available is reduce from0.6 m/s to E − m/s, corresponding to a constant acceleration from E − m/s to E − m/s . This change hasa significant impact on the number of impulses and the efficiency of the maneuver. In the first case, a single impulsewith a magnitude less than half of the maximum allowed is sufficient, resulting in a ∆ v = 0 . m/s. In the lattercase, the maneuver requires 115 impulses (with thrust on for more than half of the time) split into two thrusting arcsand a total ∆ v = 0 . m/s. Closer to the conjunction, the radial component of thrust becomes more relevant, asreported in [12]. 13 PREPRINT - M
ARCH
1, 2021 (a) 16-14 orbits to CA (b) 12-10 orbits to CA(c) 8-6 orbits to CA (d) 4-2 orbits to CA
Figure 9: Impulses profile as a function of time to CAFigure 10: Test case in which the optimal maneuver occurs at the latest opposition14
PREPRINT - M
ARCH
1, 2021 (a) ∆¯ v = 0 . m/s (b) ∆¯ v = 6 E − m/s(c) ∆¯ v = 6 E − m/s (d) ∆¯ v = 3 E − m/s Figure 11: Impulses profile as a function of maximum impulse magnitude
The proposed approach is run on 2,170 test cases, using either P C, max ≤ E − , P C ≤ E − , or ∆ r CA ≥ kmas constraints. The maneuvers can be implemented starting from 2 revolutions before the encounter using up to 170impulses. The maximum impulse magnitude is × E − m/s, corresponding to a constant acceleration of E − m/s in each segment.Figures 12–14 show the distribution of relevant quantities when the three different constraints are applied. In thegraphs, data in the range − th percentile are plotted for the sake of readability. The P C, max and P C cases producesimilar results in term of median ∆ v (0.0212 m/s versus 0.0178 m/s), relative closest encounter distance (0.7818 kmversus 0.8627 km), and number of impulses (4 in both cases). The miss distance constraint is more demanding interms of median ∆ v (0.0689 m/s) and number of impulses (12), without providing on average neither a lower collisionrisk nor a lower maximum collision risk (median values of . E − and . E − , respectively). Thisindicates that when orbital knowledge statistics are reliable, collision probability constraints should be preferred overthe miss distance one. In terms of iterations, the three constraints share a similar behavior with a median of 2 majoriterations (with a maximum of 6) and 3-4 minor iterations. In most cases, only one minor iteration is sufficient for thesecond major iteration, showing that the dynamics are almost linear due to both the short time to conjunction and themaneuvers’ size. The variation of the time of closest approach is in most cases less than half-second, but in a few cases,it reaches 8 seconds. Lastly, 7 and 6 failures occur in the risk and miss distance cases, respectively. However, it isworth mentioning that these failures are not related to convergence issues in the successive convexification iterations,but they occur when the control authority is insufficient to meet the constraints.15 PREPRINT - M
ARCH
1, 2021 (a) Distribution of P C (b) Distribution of ∆ r CA (c) Distribution of ∆ v (d) Distribution of impulses Figure 12: Maximum risk constraint
A method based on lossless and successive convexification was proposed for the optimal design of collision avoidancemaneuvers (CAMs). The maneuver is modeled as a set of impulses and it is thus suitable for handling high- and low-thrust propulsion systems. Maneuvers with minimum propellant consumption were computed meeting constraintseither on an estimate of collision risk, maximum collision risk, or miss distance without any prior knowledge onthe thrust arc structure and thrust direction. Alert times from minutes to several orbital periods were considered.The proposed method’s convergence properties and efficiency were proven by optimizing CAMs for 2,170 realisticconjunctions, showing that this methodology is promising for future autonomous onboard usage. Future efforts will bedirected towards handling long-term and multiple encounters and the direct use of accelerations as decision variables.
Acknowledgement
The author acknowledges the help of S´ebastien Henry for the reconstruction of the conjunction geometries and thefeedback and advice from peers. 16
PREPRINT - M
ARCH
1, 2021 (a) Distribution of P C, max (b) Distribution of ∆ r CA (c) Distribution of ∆ v (d) Distribution of impulses Figure 13: Risk constraint17
PREPRINT - M
ARCH
1, 2021 (a) Distribution of P C (b) Distribution of P C, max (c) Distribution of ∆ v (d) Distribution of impulses Figure 14: Miss distance constraint18
PREPRINT - M
ARCH
1, 2021
Appendix
The data used for the simulations presented in Sec. 4.2–4.5 are reported below, those for the extensive simulations areavailable at github.com/arma1978/conjunction . . E + 00 − . E + 00 − . E + 03 − . E − . E + 03 3 . E − ] . E − − . E −
04 2 . E − − . E −
04 1 . E − − . E − . E − − . E −
05 1 . E − . E + 00 7 . E + 00 − . E + 03 − . E + 007 . E + 03 − . E − ] . E − − . E −
03 7 . E − − . E −
03 8 . E −
01 1 . E − . E −
05 1 . E −
03 2 . E − R = 29 . m d CA = 8 . E − km P C = 1 . E − Eq. 5a in [33] P C = 1 . E − Eq. (10) P C, max = 1 . E − Eq. (11)
References [1] Jonathan C. McDowell. The low earth orbit satellite population and impacts of the SpaceX starlink constellation.
The Astrophysical Journal , 892(2):L36, apr 2020.[2] Joseph A Haimerl and Gregory P Fonder. Space fence system overview. In
Proceedings of the Advanced MauiOptical and Space Surveillance Technology Conference . Curran Associates, Inc. Redhook, NY, 2015.[3] Fabian Schiemenz, Jens Utzmann, and Hakan Kayal. Survey of the operational state of the art in conjunctionanalysis.
CEAS Space Journal , 11(3):255–268, 2019.[4] H. Klinkrad.
Space debris models and risk analysis . Wiley Online Library, 2006.[5] Russell P. Patera and Glenn E. Peterson. Space vehicle maneuver method to lower collision risk to an acceptablelevel.
Journal of Guidance, Control, and Dynamics , 26(2):233–237, 2003.[6] Salvatore Alfano. Collision Avoidance Maneuver Planning Tool. In , pages 1–15, 2005.[7] Saika Aida. Conjunction Risk Assessment and Avoidance Maneuver Planning Tools. In , 2016.[8] Juan Antonio Cobo-Pulido, Noelia S´anchez-Ortiz, Ignacio Grande-Olalla, and Klaus Merz. CORAM: ESA’scollision risk assessment and avoidance manoeuvres computation tool. In
Advances in the Astronautical Sciences ,volume 153, pages 453–471, 2015.[9] Alessandro Morselli, Roberto Armellin, Pierluigi Di Lizia, Franco Bernelli-Zazzera, P. Di Lizia, and FrancoBernelli-Zazzera. Collision avoidance maneuver design based on multi-objective optimization. In , volume 152, 2014.[10] Claudio Bombardelli and Javier Hernando-Ayuso. Optimal Impulsive Collision Avoidance in Low Earth Orbit.
Journal of Guidance Control and Dynamics , 38(2):217–225, 2015.19
PREPRINT - M
ARCH
1, 2021[11] G. Di Mauro, D. Spiller, S. F. Rafano Carn`a, and R. Bevilacqua. Minimum-fuel control strategy for spacecraftformation reconfiguration via finite-time maneuvers.
Journal of Guidance, Control, and Dynamics , 42(4):752–768, 2019.[12] Jason A. Reiter and David B. Spencer. Solutions to Rapid Collision-Avoidance Maneuvers Constrained byMission Performance Requirements.
Journal of Spacecraft and Rockets , pages 1–9, 2018.[13] Juan Luis Gonzalo, Camilla Colombo, and Pierluigi Di Lizia. A semi-analytical approach to low-thrust collisionavoidance design. In , number October, pages 21–25, 2019.[14] Giuseppina Salemme, Roberto Armellin, and Pierluigi Di Lizia. Continuous-thrust collision avoidance manoeu-vres optimization. In
SciTech 2020 , number January, pages 1–21, 2020.[15] Stephen Boyd and Lieven Vandenberghe.
Convex Optimization . Cambridge University Press, New York, 2004.[16] Xinfu Liu, Ping Lu, and Binfeng Pan. Survey of convex optimization for aerospace applications.
Astrodynamics ,1(1):23–40, 2017.[17] Gao Tang, Fanghua Jiang, and Junfeng Li. Fuel-Optimal Low-Thrust Trajectory Optimization Using Indi-rect Method and Successive Convex Programming.
IEEE Transactions on Aerospace and Electronic Systems ,54(4):2053–2066, 2018.[18] Zhenbo Wang, Michael J. Grant, Camille E. Bergin, Gillian S. McGlothin, Spencer T. McDonald, and ZhenboWang. Minimum-Fuel Low-Thrust Transfers for Spacecraft: A Convex Approach.
IEEE Transactions onAerospace and Electronic Systems , 54(5):2274–2290, 2018.[19] Behc¸et Ac¸ıkmes¸e and Lars Blackmore. Lossless convexification of a class of optimal control problems withnon-convex control constraints.
Automatica , 47(2):341–347, 2011.[20] Mirco Rasotto, Alessandro Morselli, Alexander Wittig, Mauro Massari, Pierluigi Di Lizia, Roberto Armellin,C Y Valles, and G Ortega. Differential Algebra Space Toolbox for Nonlinear Uncertainty Propagation in SpaceDynamics. In , Germany, 2016.[21] Yuanqi Mao, Michael Szmuk, Xiangru Xu, and Behc¸et Ac¸ikmese. Successive convexification: A superlinearlyconvergent algorithm for non-convex optimal control problems. arXiv preprint arXiv:1804.06539 , 2018.[22] Yuanqi Mao, Michael Szmuk, and Behcet Ac¸ikmes¸e. A Tutorial on Real-time Convex Optimization BasedGuidance and Control for Aerospace Applications.
Proceedings of the American Control Conference , 2018-June:2410–2416, 2018.[23] Kyle Alfriend, Maruthi Akella, Joseph Frisbee, James Foster, Deok-Jin Lee, and Matthew Wilkins. Probabilityof Collision Error Analysis.
Space Debris , 1(1):21–35, 1999.[24] Xinfu Liu and Ping Lu. Solving nonconvex optimal control problems by convex optimization.
Journal ofGuidance, Control, and Dynamics , 37(3):750–765, 2014.[25] Yuanqi Mao, Daniel Dueri, Michael Szmuk, and Behc¸et Ac¸ıkmes¸e. Successive Convexification of Non-ConvexOptimal Control Problems with State Constraints.
IFAC-PapersOnLine , 50(1):4063–4069, 2017.[26] MOSEK ApS.
The MOSEK optimization toolbox for MATLAB manual. Version 9.0. , 2019.[27] F. Kenneth Chan.
Spacecraft Collision Probability . American Institute of Aeronautics and Astronautics, 2008.[28] Romain Serra, Denis Arzelier, Mioara Joldes, Jean Bernard Lasserre, Aude Rondepierre, and Bruno Salvy. Fastand accurate computation of orbital collision probability for short-term encounters.
Journal of Guidance, Control,and Dynamics , 39(5):1009–1021, 2016.[29] Ricardo Garc´ıa-Pelayo and Javier Hernando-Ayuso. Series for collision probability in short-encounter model.
Journal of Guidance, Control, and Dynamics , 39(8):1904–1912, 2016.[30] Roberto Armellin and Pierluigi Di Lizia. Probabilistic Optical and Radar Initial Orbit Determination.
Journal ofGuidance, Control, and Dynamics , 41(1):101–118, apr 2018.[31] Roberto Armellin, P Di Lizia, Franco Bernelli-Zazzera, Martin Berz, and Pierluigi Di Lizia. Asteroid closeencounters characterization using differential algebra: the case of Apophis.
Celestial Mechanics and DynamicalAstronomy , 107(4):451–470, aug 2010.[32] Yanchao He, Roberto Armellin, and Ming Xu. Bounded Relative Orbits in the Zonal Problem via High-OrderPoincar´e Maps.
Journal of Guidance, Control, and Dynamics , 42(1):1–13, nov 2018.[33] Salvatore Alfano. Review of Conjunction Probability Methods for Short-term Encounters.