Dual Quaternion Based Close Proximity Operation for In-Orbit Assembly via Model Predictive Control

This paper studies the problem of guidance and control for autonomous in-orbit assembly. A six-degree-of-freedom (6-DOF) motion control for in-orbit assembly close proximity operation between a service satellite and a target satellite is addressed in detail. The dynamics based on dual quaternion are introduced to dispose the coupling effect between translation and rotation in a succinct frame, in which relevant perturbation and disturbance are involved. With the consideration of economical principle for fuel consume, a generic control system based on model predictive control (MPC) is then designed to generate a suboptimal control sequence for rendezvous trajectory considering actuator output saturation. The stability and robustness issues of the MPC-based control system are analyzed and proved. Numerical simulations are presented to demonstrate the effectiveness and robustness of the proposed control scheme, while additional comparisons for diverse horizons of the MPC are further conducted.


Introduction
As the renewal and progress of astronautic science and technology, regular scaled spacecraft could hardly satisfy the increasing demand to explore the universe. Confronting this challenge, on-orbit service was proposed and has been improved to target on constructing large space systems [1]. Lots of space missions such as Intelligent Building Blocks for On-Orbit Satellite Servicing (iBOSS) were investigated consecutively to develop a novel spacecraft structure which can be easily implemented by means of cubic modules in orbit [2], which is illustrated in Figure 1. Hence, the aid of in-orbit assembly, particularly relevant technologies on autonomous rendezvous and docking are crucial to ensure the success of such operations, namely, the guidance and control for proximity operations.
To implement a smooth proximity operation, one of the necessities is to take full consideration of attitude and orbit motions simultaneously, particularly the coupling characters in the 6-DOF dynamics [3]. Besides, orbit perturbation and other disturbances merit weighty attention which could bring barrier to maneuver accuracy. With all these factors taken into consideration, a compact model for controller design is necessary.
Another core of such missions is to design a reliable, economically viable control algorithm. Nowadays, control theories and technologies have been developed with higher potential to incorporate cutting edge algorithms. Engineers in many industrial fields are more inclined to apply robust and transparent control process to target plant with simple controllers, and this is also the main reason that the PID controller can remain the wide popularity compared with the extensive development of other advanced control methodologies [4]. Yet, it is generally acknowledged that PID controller has limitation in performance adjustment and constraint accommodation [5,6]. From control performance perspective, the method of MPC is superior to the PID method when targeting on optimality, especially when control objective constraints are involved. Besides, practically control instructions may not response to the control instructions immediately. This effect of delay will cause overshoot problem for PID controller design. However, MPC is capable of handling this problem via prediction process. The method of MPC was proposed in 1960s and has been popularized in various industrial process and products, particularly in autopilot vehicles [7][8][9][10]. Inspired by such motivation, this research attempts to design a MPC for spacecraft proximity operation using a compact model, while considering disturbances and maneuver capability during the inorbit assembly operation.
As for proximity operation, previous works examined the relative motion by isolating the orbit motion and the attitude motion [11,12]. This technique of "divide and conquer" neglects the coupling influence between the translation and the rotation. Some works established simultaneous equations of translational and rotational motion to fix this deficiency [13]. For instance, Sun et al. designed adaptive robust controllers concerning the coupling effects and model uncertainties based on such nonlinear dynamics [14]. However, this type of dynamics indicates a mathematical representation of the superimposed motion lacking of an explicit physical meaning.
This paper takes advantages of dual quaternion for modeling based on the screw theory [15] associating with a geometric entity by involving combined position and attitude intrinsically to describe a 6-DOF motion in a compact form. Recent evidence has suggested that dual quaternion applied to robotics and satellite sheds new light on 6-DOF motion modeling for state control and estimation [16,17]. Xza et al. designed an adaptive iterative learning control for flexible spacecraft rendezvous for a noncooperative tumbling target [18]. Sun et al. examined a dual quaternionbased model for electromagnetic collocated satellites for diffraction imaging missions [19]. Furthermore, another similar mathematical representation, Special Euclidean group SE(3), can be also applied to formulate the coupled dynamics. Related works including Zhang et al. investigated a robust adaptive control method with the aid of exponential coordinates on Lie group SE(3) [20]. Peng et al. studied a predictor-based pose stabilization control for unmanned vehicles on SE(3) considering actuator delay and saturation [21]. Wang et al. designed a consensus extended Kalman filter to estimate the motion of a rigid body in a communicate network utilizing SE(3) dynamics [22]. This type of descrip-tion is essentially the homogeneous matrix of unit dual quaternion [23].
From control theory perspective, massive literature has made achievements for spacecraft pose maneuver. Many works aiming at optimizing pose maneuver trajectory were mostly covered in the frame of optimal control theory. However, most relevant research wanted to find a global optimal solution which is unnecessary in the entire process for practical use [24]. The common solutions for optimal control theory can be approximately inducted by variation method [25], maximum principle [26], and dynamic programming [27]. Under normal circumstances, the variation method and maximum principle are easier to handle linear model without complex constraints. While dynamic programming based on the Bellman optimality principle has become a fundamental core of many intelligent algorithms such as neural network, reinforced learning, and the quadratic programming in the MPC optimization [28,29]. It should be noted that another optimal control method, linear quadratic regulator (LQR), has the similar structure with MPC. The main difference between LQR and MPC is that the optimized control output is solved in a fixed time window, while MPC calculates the feasible control in a receding time horizon. In other words, LQR is an offline control method whereas MPC is an online control method. The characteristic of MPC has a better performance for a system with external disturbance and constraints. In literature [30], the advantage and disadvantage of MPC have been discussed compared with LQR for a pursuit maneuvering situation. As a stateof-art control method, MPC is designed to solve complex multivariable optimal control problem with constraints. Compared with typical robust control and adaptive control, MPC does not excessively rely on the precision of system model, and MPC can be utilized to dispose various control schemes according to the model dynamics. The proposed work outlines MPC as a general prospective method, which settles for optimal/suboptimal solution in the prediction horizon to balance the optimality and practicability.
Contributions in this paper are threefold. First, the dual quaternion is involved to model the translation and rotation with θ being the Euler angle and n being the unit Euler axis. This definition denotes a unit quaternion as well, which has the following conjugation and norm properties: Specially, the identity quaternion is addressed as q I = ½1, 0, 0, 0 T , which is a significant denotation for error attitude.
The basic quaternion operations are implemented as follows: Addition: Multiplication: It can also be represented in matrix form as: where ½⋅ ½⊗ denotes the cross product operator defined as: Quaternion cross products: Similar, can also be represented in matrix form as: 2.2. Dual Quaternion. Before defining dual quaternion, it is necessary to introduce dual number. Similar to complex numbers, dual numbers consist of real and dual parts denoted as: In terms of Clifford algebra, dual quaternion fills the real and dual part with two quaternions denoted as: Dual quaternion has a physical interpretation according to Chasles' theorem [32], and a rigid body displacement can be seen as a screw motion containing both translation and rotation in the hereunder description: in whichn = n + εðp × nÞ is the screw axis, p is a vector 3 International Journal of Aerospace Engineering vertical to n. b θ = θ + εðdÞ is the dual angle, where d is the pitch representing translation, θ representing rotation.
Referring back to quaternion, dual quaternion has similar properties to quaternion. First, we briefly address the dual quaternion basic operation. Addition: Multiplication: It can be represented in matrix multiplication in conformance with quaternion multiplication: Cross product: Similar, The conjugation and norm of a dual quaternion denote as: when kqk = 1, it is defined as unit dual quaternion which usually denotes asq

System Kinematics and Dynamics
In the first place, it is necessary to introduce the operation to be implemented and the necessary coordinate frames. The mission is assumed that the service spacecraft will maneuver to the target spacecraft in a set pose. Throughout this paper, 2 coordinate frames are shown in Figure 2 defined as: Earth Centered Inertial (ECI) frame O − X I Y I Z I : this frame is a nonaccelerating frame originated at the center of the Earth. Its X I axis points to the vernal, and the Z I axis points to the north parallel to the Earth rotation axis. Y I is given by the right-hand rule.
Spacecraft Body Frame for service and target spacecraft O − X S Y S Z S and O − X T Y T Z T : the center is located at the center of the mass. The other 3 axes of the frame are fixed to the body and aligned to the principle axes of the spacecraft.
Suppose the service spacecraft implements a pose maneuver to the target spacecraft between frame O − X S Y S Z S and frame O − X T Y T Z T . The aftermentioned subscripts indicate the frames, respectively. It can be represented by a dual quaternion containing a rotation q and a translation r as follows:q The translation vector r represented in two frames has the coordinate transformation relation r T = q * ∘ r S ∘ q, where the 3D translation vector r S extends to a quaternion form as r S = ½0, r S x , r S y , r S z T . Particularly, this transformation is applied in the latter definition of dual velocity twist and dual force for dual quaternion dynamic modeling.
The kinematic equation is given as: in which b ω = ½ω T , v T T is the dual velocity twist defined by angular velocity and translational velocity in the target frame.
As for the dynamic equation, it is given as: in which the dual forceF = ½ f , τ T represents the external forces f and torques τ. The dual inertia matrix is denoted as:Ĵ In this dynamic equation, the total dual force includes dual control force, dual force J 2 by spherical harmonic perturbation, two-body gravity, gravity gradient torque, and disturbance described as: where the dual forces due to spherical harmonic perturbation, two-body gravitation, and gravity gradient are where the r S S/I is the position vector represented in the ECI. R e is the radius of the Earth. J 2 is the index of harmonic terms.

Remark 1.
Other additional forces and torques such as atmospheric drag, solar drag disturbances are not included owing to the short operation period of time. Besides, according to the previous work in literature [33][34][35][36], the gravitation and J 2 perturbation have minimal impacts on orbit maneuver. Hence, the aforementioned disturbances can be approximatively taken as a sum of bounded disturbances exerted on the state in terms of a proper magnitude order, and this setting is also convenient for latter controller design.

Model Predictive Controller Design
In this section, based on the spacecraft dual quaternion dynamics indicated previously, we take advantage of the relevant matrix operator of dual quaternion and forward Euler method to represent a time-varying state space model for MPC design. The structure of model predictive controller for dual quaternion nonlinear system comprises discretization, prediction, and optimization. The schematic is depicted in Figure 3.
Due to the strong nonlinearity of dual quaternion dynamics, it is common to make an approximate linearization by first order Taylor expansion in a sampling period Δ t to obtain a discrete state space model as:   x (k+3 | k) x (k+N p | k) x (k+2 | k) Figure 4: Illustration of receding control approach.

International Journal of Aerospace Engineering
To guarantee the controllability, choose the state variable as x = ½ω∧, q∧ T and control input as uðtÞ =FðtÞ. Thus, we can get a concise discrete formula as: where Next, according to the initial state and control, the state information in p sampling instants for prediction can be recurred due to the discrete model as: Given the convenience for latter optimization, denote the state and control input as: We can obtain a new matrix form for prediction: where Q, R are positive constant diagonal matrices. The cost function can be deduced into a quadratic form for optimization as Hence, the close-loop system converts to a quadratic programming problem Subject to:Ĵ 4.1. Stability Analysis. Model predictive control is established using the principle of receding horizon control [37], in which the future control trajectory is optimized at each sample time by minimizing cost function subject to constraints. In this section, we elaborate the stability and robustness in two parts. First, we have to address necessary assumption for the proof.
Assumption 2. The terminal state of the receding horizon optimization can be seen as an additional constraint of xðk + N p | kÞ = 0 resulting from the control sequenceΔuðk + mÞ, m = 0, 1, 2, ⋯, N p . Assumption 3. There exists an optimal solution to the cost function for each sampling instant in which the cost function is minimized subject to the constraints.    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 Strictly based on the assumptions, the closed loop MPC system is asymptotically stable.
Proof. Consider the system without disturbance and perturbation. Unlike continuous system, the Lyapunov function for discrete system herein can be chosen as the minimum of the finite horizon cost function at each sampling instant. At time k, it can be denoted as: Namely, VðXðkÞ, kÞ = J min . It is obvious that VðxðkÞ, kÞ tends to infinity when xðkÞ tends to infinity. Similar, at time k + 1, the Lyapunov function becomes Recall the relationship between the state on time k and k + 1 which means the state of time k + 1 is driven by the previous recursive relation. Particularly, according to Assumption 3 and the Lyapunov function designed, the optimality is guaranteed by the optimization at every single instant. Hence, we have where the VðXðk + 1Þ, k + 1Þ uses the feasible control sequence ΔuðkÞ for state update for the sample time k + 1, k + 2, ⋯, k + N p − 1. It is noticeable the Vðxðk + 1Þ, k + 1Þ and VðxðkÞ, kÞ share the same state sequence as well. Thus, we have The right hand part of the unequal sign can be expanded as:

International Journal of Aerospace Engineering
In terms of Assumption 2, it can be easily obtained: Hence, the asymptotic stability of the MPC is proved.☐ Remark 4. The attitude part in the set point of reference state denotes as 1 0 0 0 ½ T , which is not a strictly mathematical defined "zero." From physical perspective, this definition has the physical meaning of zero error Euler angle, which will not disobey Assumption 2 mentioned above.
Remark 5. Assumption 3, this precondition holds water if and only if the optimal problem of solving the cost function can be described as a convex one including all the constraints. Otherwise, the stability may not be guaranteed on account of the local convergence or the overconstrainted issue [38].  Assumption 6. All the system states can be observed.
Assumption 7. The 2-norm of system disturbance is bounded, denoting as kd k k 2 ≤ D.
Consider the discrete system with disturbance as: According to equation (22) and Remark 1, the state xðk + 1Þ can be specified as a domain which is centered in the undisturbed state with a varying radius of the norm of disturbances, defined by: where ρ k is a positive definite diagonal matrix representing the multidimensional radiuses. Take a single component of the input as an example illustrated in Figure 4, in the receding process the state space function can be rewritten based on equation (27) as: Then, this approach is similar to the stability proved. However, the system is not able to be stable on origin but converged into a neighborhood of origin, which satisfying the condition of robust control invariant set (RCIS): ∀x ∈ Ω, ∃u ∈ U, if it makes f ðx, uÞ + d ∈ Ω when ∀d ∈ D, both X and U are compact and convex.
Remark 8. The robustness analysis extremely relies on the Assumption 3 in the stability analysis. Except for the constraints, the condition for an optimal solution likewise depends on the prediction horizon when the origin system may not find a solution if the prediction horizon keeps maintaining, in which case the optimization has to reduce the predicting intervals [39].
Additionally, the robustness explanation points at the MPC in this work have the capability to resist disturbances to some extent, which is not targeting on a robust MPC for a higher order of magnitude of system interference.

Simulation
In this section, we present a simulation example to demonstrate the effectiveness of the proposed approach for the 6-

11
International Journal of Aerospace Engineering convergence enable the proposed approach a viable candidate for guidance and control strategy for proximity maneuver. In the simulation, it is assumed the target body frame is always aligned with the orbit frame. The orbit parameters of chaser satellite are depicted in Table 1 in which we can obtain the target state information. The physical property and other initial state parameters are illustrated in Table 2. The desired state for attitude and position for chaser satellite is defined as the identity unit dual quaternion. The MPC system set is shown in Table 3. The disturbances exerted on chaser satellite are considered as In Figure 5, control input for the close proximity maneuver is depicted. Note that the control forces and torques never vanish for compensating the external disturbances. Figure 6 represents the time histories for the attitude and position in unit dual quaternion of the target satellite observed, whereas Figure 7 depicts the 3-D overall trajectory for the maneuver via MPC. Figure 8 represents the traces of velocities of translation and attitude. The results prove the effectiveness of our proposed approach as well as the stability and robustness. The controller could satisfy the actuator amplitude in the maneuver process.
Additionally, when dealing with the adjustment of controller parameters, we also find an important phenomenon related to the control horizon and prediction horizon. Generally, the increase of control period will reduce the controller performances such as response and stabilization time. Relative literatures have given detailed theory demonstrations, which we also verified during our research [40,41]. Thus, we choose the minimum control horizon which is also an ideal situation to apply to our simulation. However, for the sake of fairness, considering the scenario of different prediction, we maintain other system parameters but change

12
International Journal of Aerospace Engineering the prediction horizon from 6 to 15 sampling instances for comparisons. Given the control output for translational and rotational motion, the 1-norm of the control optimized is used to demonstrate the fuel costume. As is shown in Table 4, when prediction horizon is set at 7 * Ts, the required fuel is better optimized than other control groups. When the prediction horizon increases, the maneuver trends depict that the process requires more time and fuel costume.
To provide an obvious contrast, the prediction with 15 sampling instances is given to indicate the impact of bigger prediction, which is illustrated in Figures 9-12. The optimized trajectory is planned requiring more distance to maneuver. In other words, the system requires more fuel cost and more time to drive the system to the set point for stabilization. This phenomenon indicates the significance of a proper prediction horizon to balance computing efficiency and cost.

Conclusion
In this paper, aiming at on-orbit service technology, we present an autonomous guidance and control strategy for 6-DOF close proximity operation. Dual quaternion is used to parameterize translational and rotational motion of rigid spacecraft. Based on the dual quaternion dynamics, a MPC-based pose control scheme is designed with considering actuator output constraints and external disturbances. The resulting approach is verified through numerical simulations. It is shown that it can provide a generic method for autonomous spacecraft operations. In addition, we discuss influences of different prediction horizons on MPC performance in the simulation via control variable method. The result depicts that prediction horizon could directly influence the optimization. Broader prediction horizon may cause more fuel consumption and more sensitive to the external disturbances, which may introduce new lead for further studies.

Data Availability
The numerical simulation data used to support the findings of this study are included within the article.