Convex Optimization for Rendezvous and Proximity Operation via Birkhoff Pseudospectral Method

: Rapid and accurate rendezvous and proximity operations for spacecraft are crucial to the success of most space missions. In this paper, a sequential convex programming method, combined with the ﬁrst-order and second-order Birkhoff pseudospectral methods, is proposed for the autonomous rendezvous and proximity operations of spacecraft. The original nonlinear and nonconvex close-range rendezvous problem with thrust constraints and no-ﬂy zone constraints is converted into its convex version by using the sequential convexiﬁcation techniques; then, the Birkhoff pseudospectral method is used to transcribe the dynamic constraints into a series of linear algebraic equality constraints, in other words, a convex second-order conic programming problem with a relatively small condition number. Thus, the resulting problem can be accurately and efﬁciently solved by a convex solver. The simulation results indicate that the proposed methods, especially the second-order Birkhoff pseudospectral method, have obvious advantages over other methods in


Introduction
The rendezvous and proximity operation (RPO), a key technology in on-orbit space missions such as the assembly, supply, maintenance, or astronaut rotation, usually refers to the process in which the chaser spacecraft approaches the target spacecraft to complete docking or accompany a flight. In the last century, the early space programs carried out by the United States and Russia have accumulated rich RPO technological experiences, hence, some famous projects have naturally become the main research objects of scholars. Fehse [1] presented the background information, related concepts, and control strategies of early RPO tasks, such as the Gemini [2], Apollo [3], and Space Shuttle of the US, as well as the Soyuz/Progress [4] of the former Soviet Union, or Russia. Goodman [5] surveyed the overall development and flight history of the Space Shuttle from 1983 to 2005 and gave a detailed introduction to the procedural limitations and technical challenges encountered in the early space shuttle RPO missions. Based on the characteristics of different technical routes between the United States and Russia, Woffinden and Geller [6] comprehensively summarized the orbital rendezvous standards, tasks, and technologies, and concluded that autonomous technology would become the mainstream trend of RPO tasks in the future.
Generally, an RPO can be divided into five periods: the launch, phasing, close-range rendezvous, final approach, and docking [7,8]. The spacecraft's maneuver requires comprehensive navigation and control support from the ground during the phasing period. As a result, researchers are concentrating their efforts on the close-range rendezvous and final approach to prepare for final docking or berthing using autonomous rendezvous technology [9]. Autonomous rendezvous technology means that the entire rendezvous procedure relies solely on onboard equipment to create a full set of maneuver solutions, O(1), thanks to the advantage of the Birkhoff interpolation in dealing with higher-order differential equations [44][45][46].
In light of the progress of the Birkhoff interpolation [39,[43][44][45], we extend the paradigm of the Birkhoff pseudospectral method and supplement the research findings in computational efficiency and stability. Three main contributions can be summarized in our work. First, we provide a novel sequential convex programming method that combines convex techniques and the Birkhoff pseudospectral method, focusing on the rendezvous and proximity operation of the spacecraft. Second, we provide a modified pseudospectral discretization technique based on the second Birkhoff interpolation to deal with the dynamic constraints. According to the simulation results, the revised technique performs better than the first-order Birkhoff pseudospectral method. Finally, we discover that the Birkhoff pseudospectral methods have clear advantages in computing efficiency and sensitivity, in addition to the improvement in the condition number.
The content of this paper is arranged as follows. The close-range rendezvous problem is discussed in the second section, along with several convexification techniques. Section 3 presents the discretization of the optimal control problem by using the classic pseudospectral method. In Section 4, two well-conditioned pseudospectral discretization approaches, based on the Birkhoff interpolation, are provided for the rendezvous trajectory optimization. The computational efficiency and sensitivity of the proposed algorithm are described through simulation experiments and results analysis in Section 5. Finally, Section 6 provides a conclusion and evaluates the work of this paper.

Relative Equations of Motion
A close-range rendezvous maneuver scenario is discussed in this paper, in which the deputy spacecraft (called the Chaser) passes by the chief spacecraft (called the Target) without colliding. The motion of the Chaser relative to the Target is normally used to illustrate the above proximity. As shown in Figure 1, a local-vertical, local-horizontal (LVLH) coordinate system is defined at the center of the Target O T and moves with it. The x axis, or R-bar, is radially outward. The y axis, or V-bar, is perpendicular to the x axis and points in the direction of the Target , and, in some cases, it has even reached (1)  , thanks to the advantage of the Birkhoff interpolation in dealing with higher-order differential equations [44][45][46].
In light of the progress of the Birkhoff interpolation [39,[43][44][45], we extend the paradigm of the Birkhoff pseudospectral method and supplement the research findings in computational efficiency and stability. Three main contributions can be summarized in our work. First, we provide a novel sequential convex programming method that combines convex techniques and the Birkhoff pseudospectral method, focusing on the rendezvous and proximity operation of the spacecraft. Second, we provide a modified pseudospectral discretization technique based on the second Birkhoff interpolation to deal with the dynamic constraints. According to the simulation results, the revised technique performs better than the first-order Birkhoff pseudospectral method. Finally, we discover that the Birkhoff pseudospectral methods have clear advantages in computing efficiency and sensitivity, in addition to the improvement in the condition number.
The content of this paper is arranged as follows. The close-range rendezvous problem is discussed in the second section, along with several convexification techniques. Section 3 presents the discretization of the optimal control problem by using the classic pseudospectral method. In Section 4, two well-conditioned pseudospectral discretization approaches, based on the Birkhoff interpolation, are provided for the rendezvous trajectory optimization. The computational efficiency and sensitivity of the proposed algorithm are described through simulation experiments and results analysis in Section 5. Finally, Section 6 provides a conclusion and evaluates the work of this paper.

Relative Equations of Motion
A close-range rendezvous maneuver scenario is discussed in this paper, in which the deputy spacecraft (called the Chaser) passes by the chief spacecraft (called the Target) without colliding. The motion of the Chaser relative to the Target is normally used to illustrate the above proximity. As shown in Figure 1, a local-vertical, local-horizontal (LVLH) coordinate system is defined at the center of the Target      and completes the right-hand rule. The vector positions of the Chaser and Target relative to the center of the Earth O E are denoted by r T and r C , accordingly. The angular velocity of the Target is given by the vector ω, and the scalar ω indicates its magnitude. Therefore, the vector position of the Chaser satisfies Differentiating both sides of the Equation (1) twice with respect to time for the GEI coordinate system yields the Chaser's motion equation as ..
where the vectors, r represents the acceleration of the Chaser relative to the Target. The remaining three terms on the right side of Equation (2) are the Coriolis acceleration, the tangential acceleration, and the centripetal acceleration, respectively. The dot notation [·] on signals is defined as the derivative d/dt. The double dot [··] on signals represents the second-order derivative d 2 /dt 2 . The cross-product of the two vectors is indicated by the operator [×]. The following are the operation rules.
whereω andˆ. ω are the skew-symmetric matrices corresponding to ω and . ω, respectively. The following essential assumptions are made to simplify the model of the close-range rendezvous motion. Assume that 1.
the Target is in a circular orbit near the Earth, and the orbit radius of the Target r T . is much larger than the distance between the two spacecraft r; 3.
both spacecraft are regarded as mass points, and all the orbit perturbations are ignored.
The motion of the Chaser, relative to the Target, can be described by the Clohessy-Wiltshire(CW) equations, based on Equations (2) and (3) and the above assumptions. The CW equations are a set of linear, time-invariant differential equations [47]. The detailed derivation is shown in Chapter 7, reference [48], and the results are given here, directly.
where scalar T i (i = x, y, z) represents a thrust component of the Chaser's engine. The mass of the Chaser is given by m. Therefore, the relative equations of motion or dynamics equations can be written as a first-order system of variables (r, v) or a second-order system of r .. where The matrix 0 3×3 and I 3×3 denote the [3 × 3]-dimensional zero matrix and the identity matrix, respectively. The analysis of the relative motion of two bodies using the CW equations, which can offer analytical solutions, is widespread and effective [49]. However, the dynamical constraints in the close-range rendezvous process must be balanced against the events, path constraints, and other factors. The analytical solution of the CW equations becomes invalid in this instance. The RPO problem is treated as an optimal control problem (OCP), with various nonlinear and nonconvex constraints in the later subsection.

Rendezvous Trajectory Optimization
The Chaser completes the objective of approaching the Target with the shortest maneuvering time, or the least fuel consumption, under complicated constraints using the law of relative motion equations. This process is called rendezvous trajectory optimization (RTO). Let [r, v, m] ∈ R 7×1 be the state variables and T ∈ R 3×1 be the control variable. Therefore, the OCP version of the RTO in the LVLH coordinate system can be defined as follows: where scalar t f signifies the terminal time, which is set as a constant, i.e., the time of RTO is fixed. Scalar I sp is the engine-specific impulse, and g 0 is called the sea level gravitational acceleration. The operator · here denotes the modulus of a vector. The magnitude of the thrust vector, T , must not be greater than the maximum thrust amplitude T max . The objective function J is the integral of thrust T to time. The initial condition [r 0 ; v 0 ] and the terminal condition [r f ; v f ] for the Chaser are generally determined according to the mission requirements and serve as reference indexes for the accuracy of the subsequent algorithms. The radius of the no-fly zone, denoted by R s , is the minimum safe distance at which the Chaser will not collide with the Target when approaching it.
The existence of nonlinear or nonconvex terms, such as T , T/m and r , makes it challenging to obtain the optimal solution to Problem 1. Convexification techniques will be utilized in the following section to convert the original problem into a SOCP problem with a linear objective function and convex constraints.

Convexification Techniques Applied to RTO
Convexification techniques are used to handle the dynamic constraints and no-fly zone constraints in Problem 1, which was inspired by the study of Aciknese et al. [50,51]. In this study, the key convexification techniques are equivalent transformation, change of variables, and successive approximation, as stated in reference [52].

Equivalent Transformation of the Objective Function
The objective function in Problem 1 cannot become a linear function after discretization, because of the nonlinear factor T . Equivalent transformation is used to deal with this situation, as is done in [20]. A slack variable, ζ, is introduced to replace the norm, T . The new objective function, J 1 , is as follows: Aerospace 2022, 9, 505 The transformation from J to J 1 is equivalent and satisfies the following property at the optimal solution of the problem.
Proof of the equivalence and property is given in Proposition 1, reference [53].

Change of Variables for Dynamic Constraints
The nonlinear terms in dynamic constraints are replaced by linear terms, using a change of variables. The following variables are defined: The dynamics equations can be rewritten as .
. ρ = −µσ (15) where the constant 1/I sp g 0 is replaced by µ for brevity. Let x = [r, v, ρ] T ∈ R 7×1 be the state vector and u = [γ, σ] T ∈ R 4×1 be the control vector. Therefore, the new description form of the dynamic constraints is as follows: where The following is an update to Equation (10) Then, using the first-order Taylor expansion, the nonlinear term e −ρ in Equation (19) is approximated at the time, t, yielding: where ρ = ln m and the lower bound ρ 0 is ln(m 0 − µT max t). Now, we can replace Equation (9) with derived from Equation (15), Aerospace 2022, 9, 505 7 of 30 From Equation (22), minimizing J 2 is equivalent to maximizing m(t), which also implies the optimal fuel consumption. Therefore, replacing J 1 with J 2 as the objective function is identical.

Linearization of the No-Fly Zone
Setting a no-fly zone is a crucial step in ensuring the safety of the proximity process. The no-fly zone constraint in Problem 1 is defined in the LVLH coordinate system, as shown in Figure 2a. The magnitude of the position vector r also represents the distance between the two spacecraft, as the origin is set at the Target centroid. The distance must at least meet the minimum safe distance, R s , to avoid collisions. Now, we can replace Equation (9) with From Equation (22), minimizing 2 J is equivalent to maximizing ( ) m t , which also implies the optimal fuel consumption. Therefore, replacing 1 J with 2 J as the objective function is identical.

Linearization of the No-Fly Zone
Setting a no-fly zone is a crucial step in ensuring the safety of the proximity process. The no-fly zone constraint in Problem 1 is defined in the LVLH coordinate system, as shown in Figure 2a. The magnitude of the position vector r also represents the distance between the two spacecraft, as the origin is set at the Target centroid. The distance must at least meet the minimum safe distance, s R , to avoid collisions.
(a) (b) The no-fly zone constraint can also be expressed as follows:

.
T s R  r r (23) The affine function of the position vector r is denoted by ( )  The no-fly zone constraint can also be expressed as follows: The affine function of the position vector r is denoted by g(r) and approximated using the first-order Taylor expansion at a given reference point r re f . The left side of Equation (23) becomes a linear term.
The trust region constraint is added to guarantee a reasonable approximation of the original value.
r − r re f ≤ ε, (25) where the radius of trust-region ε is positive and small enough. The reference trajectory r ref commonly takes the solution of the last iteration using the SC algorithm, as mentioned in Section 4.4. Inspired by the separating hyperplane theorem [24], we define a separating hyperplane, as shown in Figure 2b, and the describing equation is where the vector r P is in the same direction as r re f and the endpoint P is on the boundary of the no-fly zone. The relationship can be denoted as Based on Equations (24)- (27), we can derive a more intuitive and linear expression of the no-fly zone constraint.
Therefore, the original Problem 1 is converted into a convex Problem 2 through the above convexification techniques.
For Problem 2, the linear time-varying system (LTV) in Equation (29) can be converted from a continuous infinite dimension to a discrete finite dimension by some discretization methods, including the zero-order hold (ZOH), first-order hold (FOH), the Runge-Kutta (RK) method, and the global pseudospectral (PS) method [26]. After discretization, Problem 2 with the convex constraints can obtain the optimal solution by solving a convex programming problem in each iteration until the convergence criterion is satisfied, via the sequential convex (SC) method. In this study, the PS method based on the Birkhoff interpolation is combined with the SC method to solve the RTO problem.

Time-Domain Transformation
In the standard PS method, the time domain is transformed from the physical time t 0 ∈[t 0 ,t f ] to the PS time τ ∈ [−1, 1]. The mapping function Γ(τ) is given as: where the time mapping factor is defined as κ = dΓ/dτ =(t f − t 0 )/2. The discretization method of dynamic constraints, or the transformation of differential equations into algebraic equations, is the focus of this section. For the convenience of explanation, we describe the convex RTO problem as a distilled optimal control problem [35], which simplifies the expression of the boundary constraints and path constraints in Problem 2. For multivariable differential systems, the dynamic constraints follow the vector representation. The time-domain transformation of dynamic constraints is defined as where the notation ' ' denotes the derivative d/dτ. On the one hand, the boundary constraints were simplified into unified equality constraints where the dimension of the zero vector 0 eq was determined by the number of equality constraints. On the other hand, the path constraints included the thrust constraints H(x(τ), u(τ)) ≤ 0 ineq (33) where the dimension of the zero vector 0 ineq depended on the number of inequality constraints. As a result, a distilled Problem 3 [39,54] is rephrased in the PS time domain as follows: Varibles :

Classical PS Method for RTO
Firstly, an arbitrary grid of distinct nodes is defined in the interval [−1, 1]. These nodes are chosen to be the roots of orthogonal polynomials, such as the Legendre polynomials and the Chebyshev polynomials [55] where the scalar N represents the order of the orthogonal polynomials and the number of grid points is N + 1. The well-conditioned Gauss-Lobatto (GL) points are chosen as the discrete grid points [35].
In the PS theory, the state function p(τ) can be approximated by the approximation function p N (τ) with the interpolation basis polynomial at these collocation points. For the GL points, the collocation points are the same as the discrete grid points [35]. Therefore, on the grid π N , we define where the Lagrange basis functions L k (τ) (k = 0, 1, 2, . . . , N) satisfy the Kronecker delta condition.
The symbol ' ∼ =' can be replaced by the symbol '=' when the order N tends toward positive infinity.
By introducing Equation (37) into Equation (36), we can get To simplify the representation, we use p k instead of p(τ k ), and Equation (38) can be written as p k = p N (τ k ).
One of the features of the PS methods, or any Galerkin method, is that the differentiation of the approximation function can be converted into the differentiation of the interpolation polynomial [56][57][58][59]. Consequently, taking the derivative of the time τ on both sides of Equation (36), we get Define D = D kj 0≤k,j≤N as the first-order PS differential matrix (PSDM), which is represented as where D in = D kj 1≤k,j≤N−1 is part of D. Matrix D is also referred to as the differential operator, and it possesses a key characteristic proposed in [14].
where D (λ) is called the λ-th order PSDM. In addition, p k (k = 0, 1, 2 . . . N) can also be denoted by the vector p.
Equation (43) is extended to a higher-order form by applying Equation (41) where p(λ) denotes the λ-th derivative of p with respect to the PS time τ. For example, for a second-order system, Equation (44) can be written as The discretized decision variables are expressed as follows: where Therefore, the dynamic constraints in Problem 3 can be discretized as Defining an overloaded function F (X, U), Equation (48) can be written as Based on Equation (43), a new matrix form of differential transformation is given by: Note that the transpose operation here is because the decision variable X is a combination of seven variables, and, after discretization, it is a matrix of 7 × (N + 1), as stated in Equation (46).
Substituting Equation (50) into Equation (49), we get Aerospace 2022, 9, 505 11 of 30 The dynamic constraints have been turned into a set of algebraic equations, using the differential transformation of the PS method. This means that the RTO problem has become an NLP problem that can be solved by specific solvers [27,28].
For the objective function Equation (21), we make an affine transformation on the time domain, according to Equation (30).
The integral term in the objective function is calculated using the Gauss quadrature formula, as described in [44].
where σ i = σ(τ i ) and The weight, w i (i = 0,1,2,· · · ,N), and the first-order differential matrix, D, as mentioned before, have many representations based on different grid π N values. Readers can refer to [39,44] for further information, owing to the article's length and major topics.
Due to the length of the paper, the specific derivation process is omitted and the results are directly used in Section 4.4. Readers can refer to Sagliano's research [58] for a detailed processing framework. Thus, the PS transformation of the boundary constraints and the path constraints on the grid π N remains as in Problem 3. As a result, Problem 3 on the discretized points τ ∈ π N can be rewritten, using the first-order differential PS method (FDPSM).

Basic Properties of Birkhoff Interpolation
As the basis of the Birkhoff PS method, the Birkhoff interpolation is regarded as a new extension of the Lagrange interpolation and Hermite interpolation [60,61]. Unlike the Lagrange interpolation in Equation (36), the Birkhoff interpolation function is a mix of various orders of the derivatives of a function [60]. Two general forms of the Birkhoff interpolation, based on the first-order or second-order boundary value problems (BVPs), were proposed in [44]. The dynamics system in the RTO can be a first-order system about (r, v), as indicated in Equation (5), or a second-order system about r, as shown in Equation (6).
Let φ(τ) be a first-order, continuously bounded function on the grid points τ ∈ π N . Therefore, a first-order Birkhoff interpolation polynomial φ N (τ) is defined as [39], where the Birkhoff interpolation basis polynomials, B 1 k (τ), k = 0, 1, . . . , N, can be regarded as the counterpart of the Lagrange polynomials [39] and must satisfy The interpolation polynomial φ N (τ) needs to meet the interpolation conditions, which can be regarded as an extension of Equation (38) 1≤j,k≤N (59) and satisfies that where the matrix D is generated by substituting the initial row of D with e 1 = (1, 0, · · · , 0). The calculation of B 1 and the property of B 1 and B 1 in in Equation (60) are omitted here, and the detailed derivation can be found in Section 4.1, reference [44]. Therefore, the Equation (56) can be rewritten in vector-matrix form: where Let ψ(τ) be a second-order, continuously bounded function on the grid points τ ∈ π N . Similarly, a second-order Birkhoff interpolation polynomial ψ N (τ) can be written as the interpolation conditions Equations (57) and (63) can be rewritten as and and satisfies that where the matrix D (2) can be produced by replacing the initial and last rows of the matrix D (2) with e 1 = (1, 0, · · · , 0) and e N = (0, 0, · · · , 1), respectively. The detailed derivation can be found in Sections 3.1 and 3.2, reference [44]. The vector-matrix form of Equation (63) is given as

The First-Order Birkhoff PS Method for RTO
The first-order Birkhoff PS method (FBPSM) is based on the first-order Birkhoff interpolation in Equation (56). For the first-order system, we approximate the state function p(τ) by applying Equation (56) to Problem 3. The approximation function p N (τ) on the grid π N can be rewritten as where the vector p = p 0 p 1 . . . p N T is regarded as the unknown optimization variable and B 1 k , k = 0, 1, . . . , N satisfies Equation (57). Taking the derivative of the time τ on both sides of Equation (70), we get As in Equation (61), the vector-matrix form of Equation (70) is given by Therefore, the vector-matrix form of Equation (71) can be given by where 1 jk ] 0≤j,k≤N represents the derivative of B 1 over τ. According to Equations (43) and (70), it can be deduced that where Applying the property of the FBPSM to Problem 3, the discretized unknown variables can be denoted by Combining the decision variables defined in Equation (46), the dynamic constraint in Equation (51) is rewritten as where Therefore, Problem 3 on the discretized points τ ∈ π N is given by using the first-order Birkhoff PS method (FBPSM).

The Second-Order Birkhoff PS Method for RTO
The second-order Birkhoff PS method (SBPSM) is based on the second-order Birkhoff interpolation in Equation (63). The Equations (13) and (14) in the dynamic constraint can be expressed in a second-order form. ..
Let set q(τ) be a component of r. Applying Equation (63) to Equation (81) on the grid π N , we get where the vector q (2) T is the unknown optimization variable, q N (τ) is the approximation function, and B 2 k (τ), k = 0, 1, . . . , N satisfies the condition in Equation (64). Differentiating twice on both sides of Equation (82), we get ..
The vector-matrix form of Equation (82) can be written as where The vector-matrix form of Equation (83) is given as where where Applying the SBPSM to Problem 3, the unknown optimization variables are defined as The new discretized representation of the decision variables is given as Thus, Equation (81) is rewritten as where . In addition, the FBPSM is used to deal with the remaining first-order equation, Equation (15). We define the unknown variable as The new discretized representation of the decision variables ρ is given as Thus, Equation (15) can be rewritten as where Therefore, Problem 3 on the discretized points τ ∈ π N is given by using the secondorder Birkhoff PS method (SBPSM).

Sequential Convex Algorithm for RTO
In this section, the sequential convex (SC) algorithm is used to solve the finite convex problem after the convexification and discretization. The main idea is to solve a series of convex subproblems to converge to the optimal solution {x * , u * }. The iteration process via the SC algorithm is described in Algorithm 1. Step 1: Calculate the initial trajectory r 0 without the no-fly constraint and set it as the reference trajectory r ref of the first iteration. Go to step 2.
Step 2: Start the loop until meeting the trust-region constraint.
1. For k = 1:Maxiter; 2. Calculate the current trajectory {r k ,u k }; 3. Calculate the maximum distance ∆d max between r k and r ref ; 4. Check the trust-region constraint: If ∆d max ≤ ε established, exit the loop and go to step 3; otherwise, set r ref = r k and continue. 5. End.
Step 3: Evaluate the accuracy of the solution e. Go to step 4.
Step 4: Output the result. The optimal trajectory and optimal control are {r*,u*} ; the optimal objective value is j*. Remark 1. The convergence of the SC algorithm is significantly influenced by a good initial trajectory. We propose to cancel the no-fly zone constraint to achieve this goal, which is also given by the subsequent simulation. Therefore, the first subproblem is given as Remark 2. In step 2, the subproblem required to be solved in the k-th iteration with all constraints is as follows: The parameter maxiter in step 2 represents the SC algorithm's maximum number of iterations. The SC algorithm does not converge to the optimal solution, fulfilling the trust-region condition, if k takes the value of maxiter. It's feasible to replace more relaxed convergence conditions or add the mesh points to improve the convergence of the SC algorithm.

Remark 3.
In step 3, we linearly interpolate the optimal control u k of the iterative result of the SC algorithm, integrate the original dynamic equation with fourth-order Runge-Kutta to obtain the terminal statex k , and then compare the terminal valuer k (τ f ) to the given accurate value r(τ f ) to calculate the terminal error e, which is used as the algorithm's accuracy index.
The stepsize of the fourth-order Runge-Kutta integral was 0.001 in the PS time domain.

Remark 4.
To simplify the problem, only one no-fly zone constraint is considered in step 2, and the geometric shape is a sphere or an ellipsoid. Readers can learn more about the excellent work on multiple no-fly zones and complex geometries through [62,63].
For the ease of comprehension, we provide the complete process framework of the SC algorithm, using the FBPSM discretization method as an illustration, in Figure 3.  ( ) , Step 1 Step 2 Step 3 Figure 3. The complete process framework of the SC algorithm takes the FBPSM discretization method as an example.

Numerical Simulation
In this section, two cases are given to illustrate the feasibility and advantages of the proposed SC algorithm, based on the first-order and second-order Birkhoff pseudospectral discretization methods (FBPSM and SBPSM). As comparison methods, other discretization techniques, such as the zero-order hold (ZOH) [26] and the classic pseudospectral method, as mentioned in Section 3 (FDPSM), are chosen.
All numerical simulations were performed on a desktop with an Intel (R) Core (TM) i5-10400 (CPU 2.90 GHz) and with MATLAB version 2020a. The modeling tool is Yalmip [64] and the solver tool is Mosek [21]. Additionally, we found through research that the preprocessing option of the MOSEK solver significantly affects the results of the solution. Therefore, we changed the value of the option "MSK_IPAR_PRESOLVE_USE" to OFF and the reason is given in the following section.
Nondimensionalization is an effective method to simplify the calculation and improve accuracy. Table A1 in Appendix A provides a list of the principal dimensionless units. The radius of the earth, e R , is 6378.14 km, and the gravitational acceleration, 0 g , is 9.807 m/s. The value of the geocentric gravitational constant, M G , is  Table A2, Appendix A. Note that the variables contained in this section tables and figures are defined in the LVLH coordinate system.

Case 1: Excellent Computational Efficiency and Solution Accuracy under Different Grid Points
We set up several numerical tests with different grid points to evaluate the accuracy of the solution and the computational efficiency of the proposed methods. In general, the grid points of the ZOH method are evenly distributed. The Runge phenomenon [54] becomes increasingly prevalent as the number of grid points rises, which reduces the

Numerical Simulation
In this section, two cases are given to illustrate the feasibility and advantages of the proposed SC algorithm, based on the first-order and second-order Birkhoff pseudospectral discretization methods (FBPSM and SBPSM). As comparison methods, other discretization techniques, such as the zero-order hold (ZOH) [26] and the classic pseudospectral method, as mentioned in Section 3 (FDPSM), are chosen.
All numerical simulations were performed on a desktop with an Intel (R) Core (TM) i5-10400 (CPU 2.90 GHz) and with MATLAB version 2020a. The modeling tool is Yalmip [64] and the solver tool is Mosek [21]. Additionally, we found through research that the preprocessing option of the MOSEK solver significantly affects the results of the solution. Therefore, we changed the value of the option "MSK_IPAR_PRESOLVE_USE" to OFF and the reason is given in the following section.
Nondimensionalization is an effective method to simplify the calculation and improve accuracy. Table A1 in Appendix A provides a list of the principal dimensionless units. The radius of the earth, R e , is 6378.14 km, and the gravitational acceleration, g 0 , is 9.807 m/s. The value of the geocentric gravitational constant, G M , is 3.986012 × 10 14 m 3 /s 2 . The satellite completes the maneuvering and rendezvous at an orbital altitude, H 0 , of 600 km. The initial mass of the Chaser, M 0 , is 1000 kg. The initial parameters of the simulations are dimensionless, as shown in Table A2, Appendix A. Note that the variables contained in this section tables and figures are defined in the LVLH coordinate system.

Case 1: Excellent Computational Efficiency and Solution Accuracy under Different Grid Points
We set up several numerical tests with different grid points to evaluate the accuracy of the solution and the computational efficiency of the proposed methods. In general, the grid points of the ZOH method are evenly distributed. The Runge phenomenon [54] becomes increasingly prevalent as the number of grid points rises, which reduces the precision of the solution. Although the FDPSM improves the situation to a certain extent, it destroys the sparsity of the system, which is manifested in the exponential growth of the condition number of the coefficient matrix and leads to a longer solution time. This situation can be improved by using the proposed methods, based on the Birkhoff interpolation. Given that reference [39] contains the outcomes of the first-order system, we just provide the change law of the condition number of the coefficient matrix under the second-order system, as depicted in Figure A1, Appendix A. The definitions of three types of coefficient matrices can be found in Section 5 of [39]. The condition number, based on the second-order Birkhoff PSIM B 2 , shows an obvious decrease from O N 4 to O √ N , even with an invariable constant O(1), which indicates that the Birkhoff PS method also has a significant effect on improving the condition number of the second-order system. First, we divided the tests into seven groups, each with a number of grid points ranging from 30 to 210. The total time, T total , required to arrive at the optimal solution is utilized as the evaluation criterion for computational efficiency, and it primarily involves the two stages of problem-modeling and solver-solving. The expression is as follows: where the term (T 0 model + T 0 solver ) represents the time of solving Subproblem 1 in Equation (95), the intermediate-term denotes the time of solving Subproblem 2 in Equation (96), and S iter is the number of iterations. The term T other is the time of the remaining parts in Figure 3 and is equal, in principle, for the above four methods. The terminal state error, e, which is stated in Equation (97) of remark 3, is used to determine the accuracy of the solution. The value of the performance index, J obj , is related to the optimality of the solution. The statistical results are shown in Table 1. The results in Table 1 demonstrate that the number of grid points, N, does have an impact on the computational efficiency and some new findings are given here. The indicators of the final three methods show a minimal difference, whereas the ZOH method's total time, T total , has nearly tripled and has a significant terminal error, according to the findings of the first four groups. However, when N is selected as 150 or more, the total solution time of the FDPSM rises sharply, which is roughly three times that of the Birkhoff PS methods.
According to our analysis, the Birkhoff PS methods decrease the condition number of the coefficient matrix, enhance the dynamics equations, and thus shorten the modeling time of the optimization problems. Additionally, we discover that the total solution time of the SBPSM is less than the FBPSM. In comparison to the FDPSM and FBPSM, the SBPSM reduces the number of unknown variables, which, in turn, reduces the dimension of the problem and ultimately speeds up the solving process. Figure 4 illustrates the model time and solver time at various grid points and supports the validity of the conclusion. For example, the model time of the FDPSM in Figure 4a is greater than that of the FBPSM and SBPSM, when N is set as 150 or more. In contrast, the two first-order approaches take longer to solve the problem than the SBPSM in Figure 4b. time of the optimization problems. Additionally, we discover that the total solution time of the SBPSM is less than the FBPSM. In comparison to the FDPSM and FBPSM, the SBPSM reduces the number of unknown variables, which, in turn, reduces the dimension of the problem and ultimately speeds up the solving process. Figure 4 illustrates the model time and solver time at various grid points and supports the validity of the conclusion. For example, the model time of the FDPSM in Figure 4a is greater than that of the FBPSM and SBPSM, when N is set as 150 or more. In contrast, the two first-order approaches take longer to solve the problem than the SBPSM in Figure 4b.  A perfect discretization method has a modest terminal state error, in addition to a brief total solution time. Figure 5 demonstrates that, although the terminal error of ZOH steadily lowers as the number of grid points increases, it is still higher than that of other discretization approaches. In terms of the solution accuracy, this result demonstrates that the ZOH is inferior to the pseudospectral discretization method. The outcomes also indicate that the new discrete approach, which is based on the Birkhoff interpolation, inherits the superior qualities of the conventional pseudospectral method. Figure 6 displays the position errors as measured by the difference between the optimal trajectory and the propagated trajectory in both the x-axis and y-axis directions. The solution accuracy of the SBPSM is higher because, during the entire maneuver, the error of the SBPSM is smaller than the result of the ZOH in both directions. A perfect discretization method has a modest terminal state error, in addition to a brief total solution time. Figure 5 demonstrates that, although the terminal error of ZOH steadily lowers as the number of grid points increases, it is still higher than that of other discretization approaches. In terms of the solution accuracy, this result demonstrates that the ZOH is inferior to the pseudospectral discretization method. The outcomes also indicate that the new discrete approach, which is based on the Birkhoff interpolation, inherits the superior qualities of the conventional pseudospectral method. Figure 6 displays the position errors as measured by the difference between the optimal trajectory and the propagated trajectory in both the x-axis and y-axis directions. The solution accuracy of the SBPSM is higher because, during the entire maneuver, the error of the SBPSM is smaller than the result of the ZOH in both directions. time of the optimization problems. Additionally, we discover that the total solution time of the SBPSM is less than the FBPSM. In comparison to the FDPSM and FBPSM, the SBPSM reduces the number of unknown variables, which, in turn, reduces the dimension of the problem and ultimately speeds up the solving process. Figure 4 illustrates the model time and solver time at various grid points and supports the validity of the conclusion. For example, the model time of the FDPSM in Figure 4a is greater than that of the FBPSM and SBPSM, when N is set as 150 or more. In contrast, the two first-order approaches take longer to solve the problem than the SBPSM in Figure 4b.  A perfect discretization method has a modest terminal state error, in addition to a brief total solution time. Figure 5 demonstrates that, although the terminal error of ZOH steadily lowers as the number of grid points increases, it is still higher than that of other discretization approaches. In terms of the solution accuracy, this result demonstrates that the ZOH is inferior to the pseudospectral discretization method. The outcomes also indicate that the new discrete approach, which is based on the Birkhoff interpolation, inherits the superior qualities of the conventional pseudospectral method. Figure 6 displays the position errors as measured by the difference between the optimal trajectory and the propagated trajectory in both the x-axis and y-axis directions. The solution accuracy of the SBPSM is higher because, during the entire maneuver, the error of the SBPSM is smaller than the result of the ZOH in both directions.  (a) the x-axis direction (b) the y-axis direction Furthermore, we investigated the impact of the grid points count on the no-fly zone and thrust constraints. The findings indicate that increasing the grid points was necessary to improving the practicality of the ideal solution and the maneuverability. Readers can refer to Appendix B for details.
In summary, to ensure that the constrained problem has an optimal solution, we need to set the number of grid points, N , to be large enough, which will reduce the computational efficiency and solution accuracy. The proposed methods based on the first-order and second-order Birkhoff interpolation can significantly improve the efficiency of the solution, while maintaining high accuracy.

Case 2: Sensitivity Analysis under Different Conditions
In this section, the performance of the proposed algorithm was verified by a sensitivity analysis under different test conditions. The main measurement indicators were the total solution time, the number of iterations, and the terminal state error. Firstly, two conditions were considered: the type of grid points, and the shape of the no-fly zone. By adding the perturbation y to the initial position in the y-axis direction, a Monte Carlo analysis was conducted. The grid points were chosen between LGL and CGL [26] for the fuel optimization problem over a fixed time interval. In addition to the sphere, the ellipsoid  (a) the x-axis direction (b) the y-axis direction Furthermore, we investigated the impact of the grid points count on the no-fly zone and thrust constraints. The findings indicate that increasing the grid points was necessary to improving the practicality of the ideal solution and the maneuverability. Readers can refer to Appendix B for details.
In summary, to ensure that the constrained problem has an optimal solution, we need to set the number of grid points, N , to be large enough, which will reduce the computational efficiency and solution accuracy. The proposed methods based on the first-order and second-order Birkhoff interpolation can significantly improve the efficiency of the solution, while maintaining high accuracy.

Case 2: Sensitivity Analysis under Different Conditions
In this section, the performance of the proposed algorithm was verified by a sensitivity analysis under different test conditions. The main measurement indicators were the total solution time, the number of iterations, and the terminal state error. Firstly, two conditions were considered: the type of grid points, and the shape of the no-fly zone. By adding the perturbation y to the initial position in the y-axis direction, a Monte Carlo analysis was conducted. The grid points were chosen between LGL and CGL [26] for the fuel optimization problem over a fixed time interval. In addition to the sphere, the ellipsoid Furthermore, we investigated the impact of the grid points count on the no-fly zone and thrust constraints. The findings indicate that increasing the grid points was necessary to improving the practicality of the ideal solution and the maneuverability. Readers can refer to Appendix B for details.
In summary, to ensure that the constrained problem has an optimal solution, we need to set the number of grid points, N, to be large enough, which will reduce the computational efficiency and solution accuracy. The proposed methods based on the firstorder and second-order Birkhoff interpolation can significantly improve the efficiency of the solution, while maintaining high accuracy.

Case 2: Sensitivity Analysis under Different Conditions
In this section, the performance of the proposed algorithm was verified by a sensitivity analysis under different test conditions. The main measurement indicators were the total solution time, the number of iterations, and the terminal state error. Firstly, two conditions were considered: the type of grid points, and the shape of the no-fly zone. By adding the perturbation ∆y to the initial position in the y-axis direction, a Monte Carlo analysis was conducted. The grid points were chosen between LGL and CGL [26] for the fuel optimization problem over a fixed time interval. In addition to the sphere, the ellipsoid was also selected for the shape of the no-fly zone. Table A2 in Appendix A displays the shape characteristics of no-fly zones. The value of the perturbation ∆y (m) was randomly generated within the range of [−1,1] m and the number of simulation groups was set as 300. The different simulation settings are denoted by the titles of the subgraphs in Figures 7  and 8, and three pseudospectral methods are described in the legend. The ZOH technique is not considered in this section, due to the uniform distribution of the grid points and the difficulty in meeting the accuracy requirements with the maximum number of iterations, due to the huge terminal position error.
Aerospace 2022, 9, x FOR PEER REVIEW 22 of 31 was also selected for the shape of the no-fly zone. Table A2 in Appendix A displays the shape characteristics of no-fly zones. The value of the perturbation (m) y was randomly generated within the range of [−1,1] m and the number of simulation groups was set as 300. The different simulation settings are denoted by the titles of the subgraphs in Figures  7 and 8, and three pseudospectral methods are described in the legend. The ZOH technique is not considered in this section, due to the uniform distribution of the grid points and the difficulty in meeting the accuracy requirements with the maximum number of iterations, due to the huge terminal position error.  It is clear from Figure 7 that the FBPSM and SBPSM have more stable total solution times than the FDPSM, as evidenced by the roughly linear growth and square growth was also selected for the shape of the no-fly zone. Table A2 in Appendix A displays the shape characteristics of no-fly zones. The value of the perturbation (m) y was randomly generated within the range of [−1,1] m and the number of simulation groups was set as 300. The different simulation settings are denoted by the titles of the subgraphs in Figures  7 and 8, and three pseudospectral methods are described in the legend. The ZOH technique is not considered in this section, due to the uniform distribution of the grid points and the difficulty in meeting the accuracy requirements with the maximum number of iterations, due to the huge terminal position error.  It is clear from Figure 7 that the FBPSM and SBPSM have more stable total solution times than the FDPSM, as evidenced by the roughly linear growth and square growth It is clear from Figure 7 that the FBPSM and SBPSM have more stable total solution times than the FDPSM, as evidenced by the roughly linear growth and square growth trends once the grid point exceeds 100, respectively. At the same time, the SBPSM performs marginally better than the FBPSM, which is in line with the analysis in the preceding section. The number of iterations of the three methods in Figure 8 displays a "steady" fluctuation, but the interpretation is different. The fact that the number of iterations fluctuates between five and eight for the FBPSM and SBPSM indicates that the two approaches can steadily converge to the optimal solution. The approach of the FDPSM cannot ensure a stable convergence to the optimal solution since the number of iterations varies between the maximum and the finite number of iterations (Maxiter = 20). The findings demonstrate that the SC algorithm, based on the FBPSM and SBPSM, is suitable for the four circumstances listed and can maintain a faster solution and more stable convergence when there are more grid points.
The purpose of the Monte Carlo simulation tests is to evaluate the sensitivity of these three methods in dealing with the initial state perturbations. The tests are divided into two scenarios: N = 90 and N = 150. Each test only changes the initial position in the y-axis direction, and the perturbations are randomly generated. The total solution time, T total , and the terminal position error, e y , are chosen as the evaluation indices. The results are presented in Appendix A.
A visual evaluation model is created in Figure A2 by using the total solution time as the vertical axis and the terminal position error as the horizontal axis. The overall performance of this method is better the closer the results are to the lower left corner. In scenario 1, as depicted in Figure A2a, the output from the three algorithms reveals two identical clustering zones, with an error line of 0.4 m serving as the dividing line. The increase in the terminal position accuracy is accompanied by a jump in the total solution time, which rises from the right side of the dividing line to the left side of the dividing line, with an increase of roughly two seconds. The clustering regions' findings indicate that, while the solution accuracy may increase in some ranges, the total solution time may remain unchanged. Figure A2b depicts the differential distribution between the FDPSM and the Birkhoff pseudospectral methods results in scenario 2. Although the majority of the FDPSM results are clustered around T total = 60 s, there are some scattered points in a wider range in both the horizontal and vertical directions. The results of the FBPSM and the SBPSM are more concentrated, and the advantage in the total solution time can be perceived. To some extent, the approach of the FDPSM is more sensitive to the initial position disturbance. Figure A3 illustrates the findings of our quantitative analysis of the Monte Carlo simulation data, which is an addition to the results of the qualitative study discussed above. The terminal position errors of the three methods are primarily spread in the range of [0,1.2] m depicted in Figure A3a, with a median value of roughly 0.5 m. The highest limit of the FBPSM is closer to 15 s in terms of total solution time, whereas the majority of the SBPSM outcomes (about 75%) are 1.5 s faster than the other two approaches. When the number of grid points is set to 150, in Figure A3b, the terminal position error of the FDPSM spreads to 1.5 m, the median value of the total solution time rises from 10 s to 60 s, and the properties of the centralized distribution deteriorate. The original error level is maintained by the FBPSM and SBPSM, and the median value of the total solution time is 21 s and 19 s, respectively. In terms of both the solution accuracy and computational efficiency, the Birkhoff pseudospectral approaches are more stable than the FDPSM.
Through sensitivity analysis, it is found that, compared with the FDPSM, the Birkhoff pseudospectral methods have the following advantages: (i) they can converge to the optimal solution steadily at a faster speed and are suited for two different types of grid points and no-fly zones. (ii) In the case of the initial position disturbance, they can maintain a more stable terminal position error distribution and higher computational efficiency. (iii) Under the same conditions, the second-order Birkhoff pseudospectral method (SBPSM) is better than the first-order Birkhoff pseudospectral method (FBPSM), both in solution time and the robustness against the initial position disturbances.

Conclusions
With the increase in the series of space missions, such as the space station construction, on-orbit maintenance, material supply, and astronaut rotation, in the future, the RPOs will play a more important role. The popularization of automatic rendezvous technology will be of great significance to save human and material resources. The method proposed in this paper combines convex optimization technology with the Birkhoff pseudospectral method. The advantages of convex optimization techniques, such as polynomial complexity and global optimal solution, are inherited. At the same time, the pseudospectral method, based on the Birkhoff interpolation, overcomes the shortcomings of the original, traditional differential operator pseudospectral method. With the increase of the polynomial order, N, the proposed methods can not only maintain the solution accuracy but also greatly improve the computational efficiency. Especially for the dynamic system, the second-order Birkhoff pseudospectral method will have a more obvious improvement effect. Under certain conditions, the improvement of the computational efficiency can be expanded by up to three times. This work discovery will provide a possible basis for realizing spacecraft autonomous online real-time trajectory planning. Further exploring the computational performance of the proposed method and striving for practical onboard applications will be the goal of our future work.   Table A2. Initialization parameters of the RTO problem for Section 5.1.

Initialization Parameters Value
Initial Mass m 0 = 1000/M u Initial and Terminal Time   Figure A4a illustrates that the no-fly zone constraint fails when N is 30 because the Chaser's optimal trajectory, as determined by the FDPSM, passes across the no-fly zone. The optimal trajectory, as seen in Figure A4b, flies around the boundary of the no-fly zone and finally reaches the desired position safely when there are 120 grid points. This demonstrates that the effectiveness of the no-fly zone constraint is correlated with the number of grid points.  Figure A4a illustrates that the no-fly zone constraint fails when N is 30 because the Chaser's optimal trajectory, as determined by the FDPSM, passes across the no-fly zone. The optimal trajectory, as seen in Figure A4b, flies around the boundary of the no-fly zone and finally reaches the desired position safely when there are 120 grid points. This demonstrates that the effectiveness of the no-fly zone constraint is correlated with the number of grid points.

Appendix B
Generally, the distribution characteristics of the Chebyshev-Gauss-Lobatto (CGL) grid points cluster at both ends and are sparse in the middle. As a result, when N is too small, there are few interpolation points near the no-fly zone, causing the trajectory to directly cross this area, as illustrated in Figure A5. Finally, the optimal trajectory fails to meet the no-fly zone constraint, and the two satellites may collide in the close-range rendezvous mission. There are enough interpolation nodes close to the no-fly zone's boundary when N is a suitable value. Through the successive iterative optimization process, the Chaser satellite avoids the no-fly zone along the optimal trajectory and reaches the endpoint, as shown in Figure A6. The results and conclusions obtained by using the FBPSM and SBPSM are consistent with the above findings using the FDPSM. Therefore, after assessing the computational efficiency and precision of the solution, the Birkhoff pseudospectral method is unquestionably a better option when the number of grid points must be raised to satisfy the requirements. Generally, the distribution characteristics of the Chebyshev-Gauss-Lobatto (CGL) grid points cluster at both ends and are sparse in the middle. As a result, when N is too small, there are few interpolation points near the no-fly zone, causing the trajectory to directly cross this area, as illustrated in Figure A5. Finally, the optimal trajectory fails to meet the no-fly zone constraint, and the two satellites may collide in the close-range rendezvous mission. There are enough interpolation nodes close to the no-fly zone's boundary when N is a suitable value. Through the successive iterative optimization process, the Chaser satellite avoids the no-fly zone along the optimal trajectory and reaches the endpoint, as shown in Figure A6. The results and conclusions obtained by using the FBPSM and SBPSM are consistent with the above findings using the FDPSM. Therefore, after assessing the computational efficiency and precision of the solution, the Birkhoff pseudospectral method is unquestionably a better option when the number of grid points must be raised to satisfy the requirements.   Generally, the distribution characteristics of the Chebyshev-Gauss-Lobatto (CGL) grid points cluster at both ends and are sparse in the middle. As a result, when N is too small, there are few interpolation points near the no-fly zone, causing the trajectory to directly cross this area, as illustrated in Figure A5. Finally, the optimal trajectory fails to meet the no-fly zone constraint, and the two satellites may collide in the close-range rendezvous mission. There are enough interpolation nodes close to the no-fly zone's boundary when N is a suitable value. Through the successive iterative optimization process, the Chaser satellite avoids the no-fly zone along the optimal trajectory and reaches the endpoint, as shown in Figure A6. The results and conclusions obtained by using the FBPSM and SBPSM are consistent with the above findings using the FDPSM. Therefore, after assessing the computational efficiency and precision of the solution, the Birkhoff pseudospectral method is unquestionably a better option when the number of grid points must be raised to satisfy the requirements.    Figure A7 depicts the variation curve of the component and the magnitude of the Chaser's thrust when the number, N , is taken as 30 and 120, respectively. First, from the two subgraphs on the right, it can be found that the slack variable,  , is completely consistent with the change curve of the thrust magnitude, T . This illustrates that the lossless convexity operation is effective and verifies that the solution of the relaxed Problem 2 is equivalent to the original Problem 1. From the subgraphs of the thrust component, it can be found that the thrust mainly acts on the x -axis and the y -axis, which is also in line with the characteristics of the maneuvering transfer in the xOy plane, with minimum fuel consumption. Moreover, when the number of grid points is small, the change  Figure A7 depicts the variation curve of the component and the magnitude of the Chaser's thrust when the number, N, is taken as 30 and 120, respectively. First, from the two subgraphs on the right, it can be found that the slack variable, ζ, is completely consistent with the change curve of the thrust magnitude, T . This illustrates that the lossless convexity operation is effective and verifies that the solution of the relaxed Problem 2 is equivalent to the original Problem 1. From the subgraphs of the thrust component, it can be found that the thrust mainly acts on the x-axis and the y-axis, which is also in line with the characteristics of the maneuvering transfer in the xOy plane, with minimum fuel consumption. Moreover, when the number of grid points is small, the change in thrust is relatively stable and slow, especially when passing through the no-fly zone, when the working time of the whole engine is about 80 s. In sharp contrast, when N is 120, the satellite makes a rapid maneuver response when approaching the no-fly zone, which lasts only 15 s. From the previous result analysis, it is known that, in the end, the former mission failed, while the latter perfectly avoided the no-fly zone and reached the destination.
Chaser's thrust when the number, N , is taken as 30 and 120, respectively. First, from the two subgraphs on the right, it can be found that the slack variable,  , is completely consistent with the change curve of the thrust magnitude, T . This illustrates that the lossless convexity operation is effective and verifies that the solution of the relaxed Problem 2 is equivalent to the original Problem 1. From the subgraphs of the thrust component, it can be found that the thrust mainly acts on the x -axis and the y -axis, which is also in line with the characteristics of the maneuvering transfer in the xOy plane, with minimum fuel consumption. Moreover, when the number of grid points is small, the change in thrust is relatively stable and slow, especially when passing through the no-fly zone, when the working time of the whole engine is about 80 s. In sharp contrast, when N is 120, the satellite makes a rapid maneuver response when approaching the no-fly zone, which lasts only 15 s. From the previous result analysis, it is known that, in the end, the former mission failed, while the latter perfectly avoided the no-fly zone and reached the destination.