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

Authors: Chuqi Sun; Yan Xiao; Zhaowei Sun; Dong Ye
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

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.


---

## Body

## 1. 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.Figure 1
Solar array and platform extension of iBoss through in-orbit assembly.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–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 in-orbit 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 quaternion-based 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 description 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 state-of-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 simultaneously, which can pave the way for the controller. Second, considering actuator output constraints and external disturbance and perturbation, a model predictive controller is designed to track a set-point reference trajectory. Additionally, we discuss the influences of different horizons of the prediction and control for the controller performance through numerical simulation. The rest of the paper is organized as follows. Section2 devotes mathematical preliminaries of quaternion and dual quaternion. Section 3 establishes the dynamics in the framework of dual quaternion for the proximity operation. Section 4 gives detailed MPC design process and then address the stability proof and robustness analysis. Section 5 presents the simulation results to verify the effectiveness of the proposed control scheme, followed by the conclusion of this work in Section 6.


## 2. Mathematical Preliminaries

As a preliminary of this work, we concisely present relevant physical concepts and mathematical rules of quaternion and dual quaternion. Readers can also refer to literature [31] for more details.


### 2.1. Quaternion

Quaternion is a hypercomplex number to spaceℝ4.which is defined as the following:
(1)q=q0+q1i+q2j+q3k,where the imaginary units follow:
(2)ij=−ji=k,jk=−kj=i,ki=−ik=j,i2=j2=k2=−1.Such description also indicates as a combination of scalar and vector components, which can be expressed asq=q0q1q2q3T=qsqvT. From physical perspective, the quaternion is used to define a spatial rotation transformation of a coordinate frame in three dimensions constructed as:
(3)q=cosθ2nsinθ2T=cosθ2nxsinθ2nysinθ2nzsinθ2T,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:
(4)q∗=qs−qvT,q−1=q∗,q=q⋅q=q∘q∗=q∗∘q=1.Specially, the identity quaternion is addressed asqI=1,0,0,0T, which is a significant denotation for error attitude.The basic quaternion operations are implemented as follows:Addition:(5)q1+q2=qs1+qs2qv1+qv2T.Multiplication:(6)q1∘q2=qs1qs2−qv1⋅qv2,qs1qv2+qs2qv1+qv1×qv2T.It can also be represented in matrix form as:(7)q1∘q2=qs1−qv1qv1Tqv1⊗+qs1I3qs2qv2=q1∘q2,where ⋅⊗ denotes the cross product operator defined as:
(8)a⊗=0−a3a2a30−a1−a2a10.Quaternion cross products:(9)q1×q2=0,qs1qv2+qv1qs2+qv1×qv2T.Similar, can also be represented in matrix form as:(10)q1×q2=001×3qv1Tqv1⊗+qs1I3qs2qv2=q1×q2.


### 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:(11)z=a+εb,whereε2=0,ε≠0.In terms of Clifford algebra, dual quaternion fills the real and dual part with two quaternions denoted as:(12)q^=q1+εq2=qrTqdTT.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:
(13)q^=cosθ∧2sinθ∧2n∧T,in which n^=n+εp×n is the screw axis, p is a vector vertical to n. θ^=θ+ε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:(14)q^1+q^2=qr1T+qr2T,qd1T+qd2TT.Multiplication:(15)q^1⊙q^2=qr1∘qr2,qr1∘qd2+qd2∘qr2T.It can be represented in matrix multiplication in conformance with quaternion multiplication:(16)q^1⊙q^2=qr1∘04×4qr2∘qr1∘qr2qd2=q^1⊙q^2.Cross product:(17)q^1⊗q^2=qr1×qr2,qr1×qd2+qd2×qr2T.Similar,(18)q^1⊗q^2=qr1×04×4qr2×qr1×qr2qd2=q^1⊗q^2.The conjugation and norm of a dual quaternion denote as:(19)q∧∗=qr∗Tqd∗TT,q^=q^⊙q∧∗,when q^=1, it is defined as unit dual quaternion which usually denotes as q^I=qIT01×4T.


## 2.1. Quaternion

Quaternion is a hypercomplex number to spaceℝ4.which is defined as the following:
(1)q=q0+q1i+q2j+q3k,where the imaginary units follow:
(2)ij=−ji=k,jk=−kj=i,ki=−ik=j,i2=j2=k2=−1.Such description also indicates as a combination of scalar and vector components, which can be expressed asq=q0q1q2q3T=qsqvT. From physical perspective, the quaternion is used to define a spatial rotation transformation of a coordinate frame in three dimensions constructed as:
(3)q=cosθ2nsinθ2T=cosθ2nxsinθ2nysinθ2nzsinθ2T,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:
(4)q∗=qs−qvT,q−1=q∗,q=q⋅q=q∘q∗=q∗∘q=1.Specially, the identity quaternion is addressed asqI=1,0,0,0T, which is a significant denotation for error attitude.The basic quaternion operations are implemented as follows:Addition:(5)q1+q2=qs1+qs2qv1+qv2T.Multiplication:(6)q1∘q2=qs1qs2−qv1⋅qv2,qs1qv2+qs2qv1+qv1×qv2T.It can also be represented in matrix form as:(7)q1∘q2=qs1−qv1qv1Tqv1⊗+qs1I3qs2qv2=q1∘q2,where ⋅⊗ denotes the cross product operator defined as:
(8)a⊗=0−a3a2a30−a1−a2a10.Quaternion cross products:(9)q1×q2=0,qs1qv2+qv1qs2+qv1×qv2T.Similar, can also be represented in matrix form as:(10)q1×q2=001×3qv1Tqv1⊗+qs1I3qs2qv2=q1×q2.


## 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:(11)z=a+εb,whereε2=0,ε≠0.In terms of Clifford algebra, dual quaternion fills the real and dual part with two quaternions denoted as:(12)q^=q1+εq2=qrTqdTT.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:
(13)q^=cosθ∧2sinθ∧2n∧T,in which n^=n+εp×n is the screw axis, p is a vector vertical to n. θ^=θ+ε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:(14)q^1+q^2=qr1T+qr2T,qd1T+qd2TT.Multiplication:(15)q^1⊙q^2=qr1∘qr2,qr1∘qd2+qd2∘qr2T.It can be represented in matrix multiplication in conformance with quaternion multiplication:(16)q^1⊙q^2=qr1∘04×4qr2∘qr1∘qr2qd2=q^1⊙q^2.Cross product:(17)q^1⊗q^2=qr1×qr2,qr1×qd2+qd2×qr2T.Similar,(18)q^1⊗q^2=qr1×04×4qr2×qr1×qr2qd2=q^1⊗q^2.The conjugation and norm of a dual quaternion denote as:(19)q∧∗=qr∗Tqd∗TT,q^=q^⊙q∧∗,when q^=1, it is defined as unit dual quaternion which usually denotes as q^I=qIT01×4T.


## 3. 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 Figure2 defined as:Figure 2
Depiction of in-orbit assembly and coordinate frames.Earth Centered Inertial (ECI) frameO−XIYIZI: this frame is a nonaccelerating frame originated at the center of the Earth. Its XI axis points to the vernal, and the ZI axis points to the north parallel to the Earth rotation axis. YI is given by the right-hand rule.Spacecraft Body Frame for service and target spacecraftO−XSYSZS and O−XTYTZT: 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 frameO−XSYSZS and frame O−XTYTZT. 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:
(20)q^=q12q∘r.The translation vectorr represented in two frames has the coordinate transformation relation rT=q∗∘rS∘q, where the 3D translation vector rS extends to a quaternion form as rS=0,rSx,rSy,rSzT. 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:(21)q^̇=12q^⊙ω^,in whichω^=ωT,vTT 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:(22)J^⊙ω^̇+ω^⊗J^ω^=F^,in which the dual forceF^=f,τT represents the external forces f and torques τ. The dual inertia matrix is denoted as:
(23)J^=03×3010000mI30100J003×30.In this dynamic equation, the total dual force includes dual control force, dual forceJ2 by spherical harmonic perturbation, two-body gravity, gravity gradient torque, and disturbance described as:
(24)F^=f^u+f^J2+f^g+f^∇g+f^d,where the dual forces due to spherical harmonic perturbation, two-body gravitation, and gravity gradient are obtained by:(25)f^g=−μmrS/ISrS/IS301×4T,f^∇g=01×43μrS/IS×JrS/ISrS/IS5T,f^J2=3μJ2Re22rS/IS501−5rSz/ISrS/IS2rSx/IS1−5rSz/ISrS/IS2rSy/IS3−5rSz/ISrS/IS2rSz/IS01×4,where the rS/IS is the position vector represented in the ECI. Re is the radius of the Earth. J2 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–36], the gravitation and J2 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.


## 4. 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 Figure3.Figure 3
Basic framework of MPC.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:
(26)ω^t+1−ω^tΔt=J∧−1−ω^t⊗J^ω^t+F^t,q^t+1−q^tΔt=12q^t⊙ω^t+1.To guarantee the controllability, choose the state variable asx=ω∧,q∧T and control input as ut=F^t. Thus, we can get a concise discrete formula as:
(27)xk+1=Atxk+Btuk,where
(28)A=I8×8−J∧−1Δtω∧⊗J^08×812Δtq∧⊙I8×8,B=ΔtI8×808×8.Next, according to the initial state and control, the state information inp sampling instants for prediction can be recurred due to the discrete model as:
(29)xk+p∣k=Atpxk+∑i=0p−1Atp−1−iBtuk+i∣k.Given the convenience for latter optimization, denote the state and control input as:(30)Xk=xk+1∣kTxk+2∣kT⋯xk+p∣kTT,Uk=uk∣kTuk+1∣kT⋯uk+p−1∣kTT.We can obtain a new matrix form for prediction:(31)Xk=Ψxk+ΘUk,whereQ,R are positive constant diagonal matrices. The cost function can be deduced into a quadratic form for optimization as
(32)minJUk=Xk−RkTQXk−Rk+UkTRUk=Ψxk−Rk+ΘUkTQΨxk−Rk+ΘUk+UkTRUk=E+ΘUkTQE+ΘUk+UkTRUk=UkTΘTQΘ+RUk+2ETQΘUk+ETQE.Hence, the close-loop system converts to a quadratic programming problem(33)minUkJxt∣0,ut.Subject to:(34)J^⊙ω^̇+ω^⊗J^ω^=F^,q^̇=12q^⊙ω^,τi≤τmax,i=1,2,3,fi≤fmax,i=1,2,3,xt∣0=q∧0,ω∧0T=q∧0,ω∧0T.


### 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 ofxk+Np∣k=0 resulting from the control sequenceΔuk+m,m=0,1,2,⋯,Np.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.
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 timek, it can be denoted as:
(35)Vxk,k=∑m=1Npxk+m∣kTQxk+m∣k+∑m=0Np−1Δuk+mTRΔuk+m.
Namely,VXk,k=Jmin. It is obvious that Vxk,k tends to infinity when xk tends to infinity. Similar, at time k+1, the Lyapunov function becomes
(36)Vxk+1,k+1=∑m=1Npxk+m+1∣k+1TQxk+m+1∣k+1+∑m=0Np−1Δuk+m+1TRΔuk+m+1.
Recall the relationship between the state on timek and k+1(37)xk+1=Axk+BΔuk,
which means the state of timek+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
(38)Vxk+1,k+1≤V¯xk+1,k+1,
where theV¯Xk+1,k+1 uses the feasible control sequence Δuk for state update for the sample time k+1,k+2,⋯,k+Np−1. It is noticeable the V¯xk+1,k+1 and Vxk,k share the same state sequence as well. Thus, we have
(39)Vxk+1,k+1−Vxk,k≤V¯xk+1,k+1−Vxk,k.
The right hand part of the unequal sign can be expanded as:(40)V¯xk+1,k+1−Vxk,k=xk+Np∣kTQxk+Np∣k−xk+1TQxk+1−ΔukTRΔuk.
In terms of Assumption2, it can be easily obtained:
(41)V¯xk+1,k+1−Vxk,k=−xk+1TQxk+1−ΔukTRΔuk<0.
Hence, the asymptotic stability of the MPC is proved.☐Remark 4.
The attitude part in the set point of reference state denotes as1000T, 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.
Assumption3, 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].


### 4.2. Robustness Analysis

Robustness analysis requires the following assumptions:Assumption 6.
All the system states can be observed.Assumption 7.
The 2-norm of system disturbance is bounded, denoting asdk2≤D¯.Consider the discrete system with disturbance as:(42)xk+1=fxk,uk,dk.According to equation (22) and Remark 1, the state xk+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:
(43)ℵk=x~k∣x~kTρk−2x~k≤1,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:
(44)x~k+1=Ax~k+BΔvk.Figure 4
Illustration of receding control approach.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 fx,u+d∈Ω when ∀d∈D, both X and U are compact and convex.Remark 8.
The robustness analysis extremely relies on the Assumption3 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.


## 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 ofxk+Np∣k=0 resulting from the control sequenceΔuk+m,m=0,1,2,⋯,Np.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.
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 timek, it can be denoted as:
(35)Vxk,k=∑m=1Npxk+m∣kTQxk+m∣k+∑m=0Np−1Δuk+mTRΔuk+m.
Namely,VXk,k=Jmin. It is obvious that Vxk,k tends to infinity when xk tends to infinity. Similar, at time k+1, the Lyapunov function becomes
(36)Vxk+1,k+1=∑m=1Npxk+m+1∣k+1TQxk+m+1∣k+1+∑m=0Np−1Δuk+m+1TRΔuk+m+1.
Recall the relationship between the state on timek and k+1(37)xk+1=Axk+BΔuk,
which means the state of timek+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
(38)Vxk+1,k+1≤V¯xk+1,k+1,
where theV¯Xk+1,k+1 uses the feasible control sequence Δuk for state update for the sample time k+1,k+2,⋯,k+Np−1. It is noticeable the V¯xk+1,k+1 and Vxk,k share the same state sequence as well. Thus, we have
(39)Vxk+1,k+1−Vxk,k≤V¯xk+1,k+1−Vxk,k.
The right hand part of the unequal sign can be expanded as:(40)V¯xk+1,k+1−Vxk,k=xk+Np∣kTQxk+Np∣k−xk+1TQxk+1−ΔukTRΔuk.
In terms of Assumption2, it can be easily obtained:
(41)V¯xk+1,k+1−Vxk,k=−xk+1TQxk+1−ΔukTRΔuk<0.
Hence, the asymptotic stability of the MPC is proved.☐Remark 4.
The attitude part in the set point of reference state denotes as1000T, 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.
Assumption3, 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].


## 4.2. Robustness Analysis

Robustness analysis requires the following assumptions:Assumption 6.
All the system states can be observed.Assumption 7.
The 2-norm of system disturbance is bounded, denoting asdk2≤D¯.Consider the discrete system with disturbance as:(42)xk+1=fxk,uk,dk.According to equation (22) and Remark 1, the state xk+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:
(43)ℵk=x~k∣x~kTρk−2x~k≤1,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:
(44)x~k+1=Ax~k+BΔvk.Figure 4
Illustration of receding control approach.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 fx,u+d∈Ω when ∀d∈D, both X and U are compact and convex.Remark 8.
The robustness analysis extremely relies on the Assumption3 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.


## 5. Simulation

In this section, we present a simulation example to demonstrate the effectiveness of the proposed approach for the 6-DOF close proximity problem. The simulations are carried out with Matlab R2021a software and modeled by Simulink. Although the simulations were not carried out using a custom solver for optimization, the execution time and the 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 Table1 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
(45)fd=0.06+0.01sin0.5t0.05+0.03sin0.5t0.04+0.02sin0.5tN,τd=0.002+0.001sinπ10t0.003−0.001sinπ10t0.002+0.003sinπ10tN⋅m.Table 1
Orbit parameters of chaser satellite.
ParametersMajor semi-axisEccentricityInclinationRAANArgument of perigeeTrue anomalyValues6998455m0.0245 (deg)0030 (deg)Table 2
Initial values in simulation.
ParametersValuesMass100kgMoment of inertiadiag19.2,17,20kg⋅m2Initial error attitude0.3772,−0.4329,0.6645,0.4783TInitial error position0,−100,0TmInitial error velocity0,0,0Tm/sInitial error angular velocity0,0,0Trad/sTable 3
MPC system set.
ParametersValuesSampling period0.1sPredicting horizon6∗TsControl horizon1∗TsThrust boundFx,Fy,Fz≤10NTorque boundτx,τy,τz≤0.1N⋅mQdiag1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1Rdiag1,1,1,1,1,1In Figure5, 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.Figure 5
Time histories of control forces and torques expressed in the body-fixed frame.Figure 6
Time histories for the position and attitude of chaser satellite in target body fixed frame.Figure 7
3-D overall trajectory of close proximity maneuver via MPC.Figure 8
Time histories for the translational and rotational velocities in the body fixed frame.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 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.Table 4
Cost comparison according to the prediction horizons.
Prediction horizonsTranslational costRotational costConvergence time6∗Ts949.912.341297∗Ts886.811.71258∗Ts875.413.611379∗Ts904.415.2914310∗Ts967.916.8514811∗Ts104518.2115212∗Ts110720.8115815∗Ts122530.65214Figure 9
Time histories of control forces and torques expressed in the body-fixed frame with the prediction of 15 sampling instances.Figure 10
Time histories of control forces and torques expressed in the body-fixed frame with the prediction of 15 sampling instances.Figure 11
3-D overall trajectory of close proximity maneuver via MPC with prediction horizon with the prediction of 15 sampling instances.Figure 12
Time histories for the translational and rotational velocities in the body fixed frame with the prediction of 15 sampling instances.


## 6. 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.



---

*Source: 1305095-2021-11-12.xml*