1 Introduction

With the rapid development of the science and technology of spacecraft, the control task and type of the orbited spacecraft are diversifying increasingly, and the system function and structure become more complicated, and lead to higher technical requirements. The model of a simple rigid body or gyrostat will no longer be used to reflect the actual mechanical model. Thus, a great number of complex models and corresponding research topics have emerged, such as multi-rigid-body systems, rigid systems with flexible attachment bodies, liquid fuels. As a kind of classic models of multi-rigid-body spacecraft, the coupled-rigid-body spacecraft can maneuver global reorientation along with the movement of the ball-in-socket joint. The coupled-rigid-body spacecraft system belongs to a type of the underactuated systems, the attitude dynamics and control of which have attracted extensive attention from researchers and scholars around the world. Krishnaprasad ^[1] was the first pioneer who began to study the optimal control problem (OCP) of multi-rigid-body system in early years. Chen and Sreenath ^[2] considered the controllability and control synthesis of two coupled rigid bodies with motions in space. Kolmanovsky et al. ^[3] studied a system of rigid bodies interconnected by one degree of freedom (DOF) rotational joints in the planar. Reyhanoglu and McClamroch ^[4] assumed that the total angular momentum of a plannar rigid body system in space was zero and developed a reorientation maneuvering strategy. Fernandes et al. ^[5] addressed the motion planning of coupled rigid bodies based on the Ritz approximation theory. Xia et al. ^[6] investigated the problem of attitude control for a rigid spacecraft which was nonlinear in dynamics with inertia uncertainty and external disturbance by designing a sliding-mode controller. Ge and Guo ^[7] used the particle swarm algorithm and spline approximation method to research the motion planning of the dual-rigid-body spacecraft. Ge and Chen ^[8] researched the nonholonomic motion planning of a free-falling cat with the Gauss-Newton method based on the Ritz approximation theory.

Recently, the pseudospectral method (PM) has become a kind of important approaches to treat the OCP and nonlinear equations ^[9-19]. The basic idea of the PM is that it utilizes the orthogonal polynomial to approximate state and control variables at the same time, then, an OCP can be converted to a discrete nonlinear programming (NLP) problem. Benson^[9] and Benson et al.^[10] proposed the Gauss pseudospectral method (GPM) for solving the OCP and derived strictly the costate estimate. Fahroo and Ross ^[11] and Gong et al.^[12] presented the Chebyshev pseudospectral method (CPM) for directly solving a generic OCP with state and control constraints, and they gave the rigorous derivation for the costate estimate. Darby et al. ^[13-14] presented an hp-adaptive PM for numerically solving OCPs of varying complexity. Tang et al. ^[15] presented the Chebyshev-Gauss pseudospectral method (CGPM) which collocated at Chebyshev-Gauss (CG) points for direct trajectory optimization of OCPs, also they rigorously derived the costate and constraint multiplier estimates by comparing it with the Karush-Kuhn-Tucker conditions of the NLP problem. Shang and Guo ^[16] adopted the CPM to study the homogeneous initial boundary value problem of a class of multi-dimensional generalized symmetric regularized long wave equations. Afify and Elgazery ^[17] used the CPM to solve a set of coupled nonlinear ordinary differential equations transformed from a non-Darcy mixed convective heat and mass transfer flow. Zhuang and Huang ^[18] proposed a hybrid algorithm combining the particle swarm optimization algorithm with the Legendre pseudospectral method (LPM) to solve the time-optimal trajectory planning problem of the underactuated rigid spacecraft. Huang et al. ^[19] transformed the propellantless rendezvous problem of a charged spacecraft subjected to the Lorentz force about an arbitrary reference orbit into an OCP, and utilized the GPM to solve this NLP problem.

The CG points^[20] are used to substitute for Legendre-Gauss (LG) collocation points of GPM, which is the reason why the term CGPM is used. Firstly, the attitude motion planning problem of two-rigid-body spacecraft coupled by a ball-in-socket joint can be converted into a discrete NLP problem based on the CGPM theory, and then the sequential quadratic programming (SQP) algorithm is used to handle this NLP problem. Since the collation point is the CG point, the Clenshaw-Curtis quadrature rule is presented to calculate the integral within the cost function. The corresponding quadrature weights can be calculated efficiently by the fast Fourier transform (FFT). For the aim of improving numerical stability and computational efficiency, the barycentric Lagrange interpolation is presented to substitute for the classic Lagrange interpolation in the approximation of state and control. The trigonometric identity is used to alleviate the numerical float errors when calculating the state differential matrix and the barycentric weights. A layered optimization strategy from a feasible solution to an optimal solution is proposed to avoid sensitivity to the initial values for the SQP algorithm. At first, the SQP algorithm is used to search for a feasible solution which is satisfied by all the constraint conditions whether they are equalities or inequalities. Then, more nodes of design variables can be obtained in those feasible points with the cubic spline interpolation. At last, the SQP algorithm is used again to solve the optimal solutions as the initial values.

In this paper, a drift free nonlinear OCP of a two-rigid-body spacecraft coupled by a ball-in-socket joint is formulated based on the theory of angular momentum conservation and the principle of minimum energy control. The rationale of CGPM is presented, and the optimal control model is converted into a discrete NLP problem. A layered optimization strategy is also introduced. The results and analysis of numerical simulation are displayed, and some concluding remarks are given.

2 OCP of coupled-rigid-body spacecraft 2.1 Attitude motion model of coupled-rigid-body spacecraft

In general, the attitude of the spacecraft can be controlled by the external torques (e.g., gas jet) or the internal torques (e.g., momentum wheel). However, we consider a more interesting case of the spacecraft model that is composed of double rigid bodies by frictionless ball-in-socket joint connection in the paper (see Fig. 1). The kinematics and dynamics of each individual body are tightly coupled with the motion of the other body in the system. The spacecraft attitude control can be carried out by the hinge joint control input. Such control inputs of hinge joint preserve the system angular momentum.

Suppose that the system consists of rigid bodies B_i(i=1, 2). The coordinate frame Γ_C is called the inertial reference frame with the origin located at the point C in space, while the coordinate frames Γ_i are called the body fixed frames to the rigid body B_i with the origin located at the mass center C_i(i=1, 2), respectively.

The meaning of each symbol in Fig. 1 is as follows:

r₀, the vector from the origin C of coordinate frame Γ_C to the mass center C₀ of the system;

r₁, the vector from the origin C of coordinate frame Γ_C to the mass center C₁ of body B₁;

r₂, the vector from the origin C of coordinate frame Γ_C to the mass center C₂ of body B₂;

d₁, the vector from the ball-in-socket joint O based on Γ₁ to the mass center C₁ of body B₁;

d₂, the vector from the ball-in-socket joint O based on Γ₂ to the mass center C₂ of body B₂;

r₁⁰, the vector from the mass center C₀ of system based on Γ_C to the mass center C₁ of body B₁;

r₂⁰, the vector from the mass center C₀ of system based on Γ_C to the mass center C₂ of body B₂;

A₁, the coordinate transformation matrix of Γ₁ relative to Γ_C;

A₂, the coordinate transformation matrix of Γ₂ relative to Γ_C.

In the absence of the mass of links connecting the coupled rigid bodies, it is assumed that the mass of rigid body B_i is m_i(i=1, 2), and the vector of any point on the rigid body B₁ that is relative to the inertial reference frame Γ_C is q₁. Then, the kinetic energy of rigid body B₁ is given by

(1)

where ρ₁ is the vector of any point on the rigid body B₁ that is relative to the body fixed frame Γ₁. Substituting q₁=r₁+A₁ρ₁ into Eq. (1) yields

(2)

where "tr" represents trace of square matrix. Similarly, we have the kinetic energy of rigid body B₂,

(3)

In Eqs. (2) and (3), I_i=∫_Biρ_iρ_i^Tdm_i(i=1, 2). Based on Eqs. (2) and (3), the kinetic energy of system can be written as follows:

(4)

By the theorem of motion of mass center, we have m₁r₁⁰+m₂r₂⁰=0. The geometrical relationships as presented in Fig. 1 are expressed as follows:

(5)

Additionally, , where and represent the third-order skew symmetric matrices. Finally, the kinetic energy of coupled-rigid-body spacecraft can be written as

(6)

where , m=m₁+m₂, and ∈=m₁m₂/m.

The attitude motion of the coupled-rigid-body spacecraft is considered, and the influence of the rotational motion on translational motion is ignored. It is assumed that the total angular momentum of the system is conserved without the external torque, where L=T is the Lagrange function ^[21]. Suppose that the initial value of the angular momentum is zero. Then, it would be zero all the time without the external torque. Therefore, the projection of angular momentum relative to the mass center C₀ of the system to the inertial reference frame Γ_C can be derived as

(7)

The direction cosine matrix of the body fixed frame Γ₂ to the rigid body B₂ relative to the body fixed frame Γ₁ to the rigid body B₁ is A.

The angular velocity of coordinate matrix can be expressed by the direction cosine matrix of the body fixed frame Γ₂ relative to the body fixed frame Γ₁,

(8)

The angular velocity ω₂=A^Tω₁+ω of rigid body B₂ can be gained by the relationship A₂=A₁A and Eq. (8). Then, substituting it into Eq. (7) yields

(9)

where J_t=A₁J₁A₁^T+A₂J₂A₂^T+A₂J₁₂^TA₁^T+A₁J₁₂A₂^Tis the generalized inertia tensor of the system. Similarly, we can have the following equation of rigid body B₂ :

(10)

The two-rigid-body spacecraft coupled by a ball-in-socket joint has six DOFs. Three DOFs of them are the attitude angles of rigid body B₁ relative to the inertial reference frame Γ_C, while the rest are the attitude angles of rigid body B₂ relative to the body fixed frame Γ₁ of rigid body B₁. In this paper, the Rodrigues parameter is used to denote the attitude angles of rigid body B₁ relative to the inertial reference frame Γ_C. Suppose α₀=1/(1+α²), where α²=α^Tα. Therefore, the attitude matrix A₁ can be denoted in the following form ^[22]:

(11)

The angular velocity of rigid body B₁ can be written as

(12)

The Rodrigues parameter is also denoted for the rigid body B₂. Then, we have

(13)

where β₀ has the similar definition pattern with α₀.

The angular velocity of motion of the ball-in-socket joint is treated as a pseudocontrol input. And the attitude angle of two rigid bodies is defined as the state variable which is denoted as . Using Eqs. (12) and (13), Eqs. (9), and (10) can be written in the matrix-vector form,

(14)

where

Since G₂(q) is nonsingular, Eq. (14) can also be written again in the form,

(15)

where defines a state matrix of the whole system, and Eq. (15) is the attitude nonholonomic motion equation of the rigid-body spacecraft coupled by a ball-in-socket joint.

2.2 Description of OCP

Without loss of generality, this paper takes the OCP of Bolza into account:

(16)

where Φ(·) represents the Mayer cost function, G(·) represents the Lagrange cost function, f(·) is the state equation constraint, φ(·) is the initial and terminal boundary constraint, and C(·) is the path constraint.

Equation (15) can be treated as the control system of the coupled rigid spacecraft. According to the controllable rank condition (Chow theorem) ^[23], there exists a full rank controllable Lie algebra in the control system, which can handle the attitude motion planning problem of the system by principle. We assume that there exists a solution u^*∈L² ([0, T]) in the control system, where L²([0, T]) is the Hilbert space of measurable vector valued functions of the form u(t), t∈[0, T]. According to the optimal control theory of minimum energy, we select the cost energy of the ball-in-socket joint as the objective function,

(17)

while the state constraint equation of Bolza in Eq. (16) can be denoted by Eq. (15).

Given out the initial and terminal configurations q₀ and (namely, the state boundary constraints), the objective function can search for the pseudo-control input , t∈[0, T]^T. Then, the attitude motion trajectory maneuvered from q₀ to q_f can be confirmed.

At the moment of start and shut of actuating motor controlling the motion of the ball-in-socket joint, the values of control should be zero. That is to say, u₀=u_f=0 (control boundary constraints). On the other hand, it is impossible that the output control values reach the infinite (it should have a maximum). Therefore, we are supposed to take ||u||_∞≤u_max into consideration.

3 CGPM 3.1 Rationale of CGPM

The CGPM converts the continuous OCP into a discrete parameter optimization problem contained algebraic constraints through a series of transformation, namely, the NLP problem, and it can be realized through the following steps.

(ⅰ) Time affine transformation

The CGPM should import a new time interval to convert the actual time interval [t₀, t_f] of system into the interval τ∈[-1, 1], and the affine transformation can be written as

(18)

(ⅱ) Variable approximation

The approximation of state and control variable CGPM takes the K CG points and τ₀=-1 as the collocation points, where the CG points are defined as follows:

(19)

The CG points are also called Fejér-2 or Filippi points ^[24]. And the advantage of distribution over the interval (-1, 1) of CG points can suppress effectively the Runge phenomenon occurred in the Lagrange interpolation, which is characterized by the dense distribution on both ends and sparse distribution in the middle.

Using a basis of (K+1) th-order Lagrange interpolating polynomials to approximate the state variables,

(20)

where the classic Lagrange interpolating polynomials are defined as

(21)

In order to improve the computational efficiency and numerical stability, the barycentric Lagrange interpolation is used to substitute for the classic Lagrange interpolation ^[25],

(22)

where the barycentric weights ξ_i are defined as

(23)

It can be found in the denominator of Eq. (23) that the calculation of (τ_i-τ_j) would lead to floating-point cancellation errors for a large K. For this reason, an effective method is presented. Firstly, the initial and terminal point of barycentric weights ξ_i^' for the CG point can be separated as follows ^[25]:

(24)

By comparing Eq. (23) with Eq. (24), it can be seen that ξ_i is equal to ξ_i^' when the denominator is i=K. Therefore, using the trigonometric identity and considering τ₀=-1 and τ_K+1=1, Eq. (23) can be expressed again in another form as

(25)

There is no derivative of control in the CGPM. Therefore, the approximation of control is simpler than the state variables. For consistency in the form, CG points are selected as collocation points for the CGPM, and the Kth-order barycentric Lagrange interpolating polynomials L_i(τ)(i=1, 2, …, K) are used as basis functions to approximate the control variables,

(26)

(ⅲ) Transformation of dynamic differential equation

In the PM, the state variables are approximated by global orthogonal polynomials. Therefore, the time derivative of them can be derived by Eq. (20),

(27)

where the classic Lagrange interpolating polynomials can be denoted by a differential matrix as

(28)

Since barycentric Lagrange interpolating is used in this paper, the differential matrix should be given by Eq. (29). It has a better numerical stability than Eq. (28) ^[26],

(29)

where k=1, 2, …, K, and i=0, 1, …, K. Substituting Eq. (29) into the dynamics equation yields the following discrete form:

(30)

It can also be found in the denominator of the first row in Eq. (29) that the calculation of (τ_k-τ_i) will lead to floating-point cancellation errors. On the other hand, τ_i (see Eq. (19)) is a triangle cosine function. For these reasons, the trigonometric identity is used to avoid phenomena and yield Eq. (31). While ξ_i and ξ_k of Eq. (29) can be calculated efficiently by Eq. (25),

(31)

(ⅳ) Terminal constraint of state

It is known in the Bolza OCP (see Eq. (16)) that the state equations contain terminal constraints. However, the time interval of Eq. (20) is [-1, 1). Thus, the approximated state does not include the terminal state X(τ_f). Therefore, the integration of the dynamics equation over the time interval [-1, 1] yields

(32)

The Clenshaw-Curtis quadrature rules are used to substitute for Gauss quadrature in the approximation of Eq. (32) ^[27],

(33)

The direct expansion of the left side of Eq. (32) can gain a linear result. The nonlinear terms in Eq. (33) are commonly used in the GPM ^{[9-10, 14]}. However, the former has a better computational efficiency than the latter, and the linear result can be written as

(34)

where τ_k presents the CG point. And μ_k^cc denotes Clenshaw-Curtis quadrature weights which can be efficiently and numerically stable by using the FFT ^[28],

(35)

where c_k is defined by

(36)

and μ_k^f₂ represents the Fejér-2 quadrature weight,

(37)

where F_K+1^-1 represents the inverse FFT (IFFT), and υ_k is defined as follows:

(38)

Finally, we assume that μ_K+1^f₂=μ₀^f₂=0 and μ_K+1^cc=μ₀^cc. And we take the real part after calculating Eq. (37).

(ⅴ) Objective function

Utilizing the Clenshaw-Curtis quadrature rules to approximate the integral part of Bolza OCP (see Eq. (16)) yields the following equation:

(39)

3.2 NLP problem

In this section, we utilize the CGPM to convert the OCP of the coupled-rigid-body spacecraft into a discrete NLP problem. Firstly, we should define the state variables of each collocation point as , the control variables as U_1K, U_2K, U_3K∈. Then, the discrete approximated cost function is given by

(40)

where μ^cc is a sequence of Clenshaw-Curtis quadrature weights.

Using the differential approximated matrix to denote the discrete state equation in the form of matrix expression,

(41)

where G_ij(i=1, 2, …, 6, and j=1, 2, 3) are the units of state matrix in Eq. (15).

Also, the terminal constraints of state are given by

(42)

Finally, the initial and terminal constraints of state and control, as well as the control constraints, are set as

(43)

where X₀ and X_f represent the first and last columns of the approximated state, respectively, and it is the same with U₀ and U_f. Here, U_maxdenotes the control maximum.

From Eqs. (40) -(43), a general NLP problem of the OCP can be expressed as

(44)

where y represents the design variables containing state variables and control variables.

3.3 Optimization strategy and process

Theoretically, we can obtain the attitude trajectories of the coupled-rigid-body spacecraft by employing the SQP algorithm.

However, the actual troubles are that the number of design variables is huge when large CG points are chosen. For example, when we take 60 CG points, there will be 558 design variables. For so many design variables, it is difficult for us to set the initial guess values for the SQP, and if the inappropriate initial values are used, we may not get a feasible solution for the SQP algorithm. In addition, the SQP algorithm does not possess the capability to search for a global optimal solution. It depends on the initial values. Generally, it only can converge to a local optimal solution that is close to the initial value. Given these reasons, a layered optimization strategy is proposed as follows.

(ⅰ) Calculating feasible solution

The essences of a feasible solution are that it does not search for the optimal solution which is satisfied with all constraint conditions (including equality and inequality constraints), and that it does not consider the actual cost function. Contrarily, it treats the converted equality constraints as the new pseudo-objective function. Therefore, the new NLP problem does not contain equality constraints, and it can be obtained from Eq. (44),

(45)

Since the PM method can obtain the higher solution accuracy at lower nodes, also, it has lower sensitivity to initial values and has fewer design variables at this time, we are supposed intensively to choose fewer nodes K₁, and utilize the SQP algorithm of the toolbox of MATLAB to obtain a feasible solution.

(ⅱ) Calculating optimal solution

Take the discrete solution gained from the first step as the interpolating nodes of cubic spline interpolation, and then we can get more CG points by applying these discrete scaled time CG points into the spline function. In other words, K₂ initial values are yielded by K₁ interpolating nodes. Then, the optimal solution can be achieved now by substituting the K₂ initial values into the SQP algorithm to search for the NLP problem (see Eq. (44)).

The layered optimization strategy is presented in Fig. 2.

4 Results of numerical simulation

In this section, the attitude motion planning problems of the two-rigid-body spacecraft coupled by a ball-in-socket joint are numerically simulated. The mass and geometric parameters of the system are given as follows: m₁=m₂=2 k, d₁=d₂=[0, 0, 1], and

The initial and terminal attitudes of the system, as well as the control constraints, are given as

The layered optimization strategy is adopted during the optimization. Here, 6 CG points are selected for the feasible solution and 60 CG points for the optimal solution. Set the maneuver time T=5 s. After simulating, the curves of the optimal solution are shown in Figs. 3, 4, and 5, where the circles denote the CG points of the optimal solution. In Fig. 3, the black solid line denotes the attitude trajectory of rigid body B₁. Similarly, the black solid line in Fig. 4 denotes the attitude trajectory of rigid body B₂. The start and end points of curves represent the initial and terminal attitude configurations, respectively. Moreover, Fig. 5 provides the curve of pseudo-control input. It is obvious that the start and end points of control curve are zero. Therefore, this solution can satisfy the requirements of zero boundary control at the moment of start and shut of actuating motor perfectly. Moreover, the control can switch between the minimal and maximal control values.

In order to validate the attitude curves, the control values in Fig. 5 are applied to Eq. (15) with the integral operation which is the fourth-order and fifth-order Runge-Kutta method. The results of integral operation are shown in Figs. 3 and 4 in the form of black dotted lines. As seen in the local enlarged figures, the integral curves from ODE45 almost overlap with the curves by the CGPM. Therefore, it is rational and valid.

For showing intuitively and comparing the feature of the CGPM, this paper also presents the calculation results of GPM. The collection points of GPM are LG points, the GPM utilizes the Gaussian quadrature to approximate the integration item in the index function, also, it uses the classic Lagrange interpolation to approximate state and control variables, and it makes use of the nonlinear terms in Eq. (33) to compute the terminal constraints of the states. The GPM also makes use of the layered optimization strategy, and the setting of related parameters is similar to the CGPM. The initial guess values of the feasible solution for both CGPM and GPM are created randomly by the "rand" function included in MATLAB. The data are recorded in Table 1, which are gained by calculating the mean of the 6th running results. In Table 1, the parameters t_{d, f}, J_f, and n_max denote the iterating time, the value of object function, and the maximal absolute value of inequality constraints after solving the feasible solutions, respectively. t_{d, o} and J_o are the iterating time and the objective function value of the optimal solution, respectively. Finally, the parameter A_e represents the maximum absolute error of terminal attitude, and T_CPU denotes the whole running time. It can be seen that the objective function value of CGPM is smaller than the one of GPM for both the feasible solution and the optimal solution. The running time of CGPM is half of the one with the GPM. The terminal attitude has already a very high accuracy 10^-14 of the maximal absolute error, and the CGPM can also reach this level easily. Thus, the comparisons in Table 1 can embody intuitively that the CGPM adopted in this paper has the higher computational efficiency and solving accuracy.

5 Conclusions

In this paper, the CG points are used to substitute for the LG collocation points of GPM. Thus, the method is called the CGPM. Firstly, the attitude motion planning problem of two-rigid-body-spacecraft coupled by a ball-in-socket joint can be converted into a discrete NLP problem based on the CGPM theory. Then, the SQP algorithm is used to handle this NLP problem. Moreover, a layered optimization strategy from a feasible solution to an optimal solution is proposed to avoid sensitivity to the initial values for the SQP algorithm. A set of simple yet efficient methods are used to improve the numerical stability and computational efficiency and avoid the floating-point cancellation errors. The results of numerical simulation and analysis in this paper show the following conclusions: (ⅰ) The nonlinear control problem of the coupled-rigid-body spacecraft can be converted into a motion programming problem of draft free system; (ⅱ) In absence of the external torque, the rotation motion of the joint interconnecting the two-rigid-body can maneuver the motion of the whole system; (ⅲ) The optimal control can meet the requirements of zero boundary control, and switch between the minimal and maximal control values; (ⅳ) The inverse integration demonstrates that the attitude curves are rational and effective; (ⅴ) The comparison between CGPM and GPM indicates that the CGPM has the higher computational efficiency. The method adopted in this paper provides another fast and efficient approach for solving problems of direct trajectory optimization and attitude optimization.