# Variable Mass Control and Parameter Identification of Spacecraft Orbit Refueling Process

Authors: Xu Han; Zhi Li; Hui Li; Longfei Huang; Yujia Pang
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

A maintainable refueling vehicle is the future development direction of the space system. In the process of fuel filling, capture docking, and configuration transformation, the most significant factor that affects the attitude control of the system is the continuous or sudden change of the angular torque of the system. In this paper, we study the control system of the variable mass body in an on-orbit service in the process of configuration transformation. The large errors of torque of inertia may make the narrow sense TEA (torque equilibrium attitude) of the system deviate greatly from the earth-oriented attitude. In order to avoid the error caused by the partial linearization of the system model, the feedback linearization method of the nonlinear system is used to design the controller to realize the tracking of the narrow sense TEA in the process of configuration transformation. Different from the traditional attitude control method, in order to avoid the high cost of the control system caused by the change of system mass characteristics and the change of system angular torque, CMG (control moment gyroscope) angular torque is introduced into the controller. We design a joint controller of attitude control and angular torque management, which can effectively stabilize the system and reduce the angular torque saturation of the attitude control system during the on-orbit service.


---

## Body

## 1. Introduction

Low cost, high reliability, fast response, and maintainability are the development direction of the future space system. The function of on-orbit refueling may become the design requirement of future spacecraft. The on-orbit refueling mission is a part of the space service support system, which takes the on-orbit spacecraft with insufficient fuel, exhausted fuel, or propulsion system failure as the application object. By means of cabin addition, fuel filling, or module replacement, the function of the spacecraft propulsion system can be supplemented or restored, and the mission capability of spacecraft on-orbit can be improved or extended. The on-orbit refueling mission consists of two parts: the replenishment vehicle and the replenished vehicle.The on-orbit refueling process usually includes (1) the capture and docking process of the replenishment vehicle and the replenished vehicle to realize the combination connection of the two vehicles and (2) the fuel filling and transfer process after the fuel tank connection. In the process of fuel filling, capture docking, and configuration change, the biggest influence on the attitude control of the system is the constant change or mutation of the angular torque of the system, especially when CMG is used as the actuator. The CMG needs to be unloaded by air jet. Because CMG absorbs the momentum of the system change quickly to saturation. And the control structure coupling is easily caused by the unloading of the jet during this period, which makes the control condition worse.The variable mass control problems of spacecraft during the refueling process can be divided into three categories: the variable mass control problems of the spacecraft docking process, the variable mass control problems of the refueling process, and the variable mass control problems of the spacecraft separation process. The change of system mass in the process of fuel transfer can be regarded as continuous, while the system mass changes suddenly in the process of capture docking and separation. Taking the system configuration change whose mass change rate is between the above two as an example, the control of the variable mass body is studied.Wu et al. [1] used the hybrid method to study the attitude stability of small satellites, in order to solve the problem of excessive output of the controller. Chu et al. [2] developed an approximate dynamic model with uncertain parameters considering the uncertainty of model parameters and proposed a robust adaptive control strategy to compensate for or reject these uncertainties, respectively. Huang and Uang [3] aimed at the problems of various disturbances and parameter changes encountered in the process of the spacecraft space mission, the sliding film control was applied to PID control, and the two were combined to complete the spacecraft attitude robust control. Qin [4] proposed a controller design method for singularly perturbed systems. Zhou and Zeng [5] proposed a new nonlinear robust H ∞ control method for spacecraft attitude maneuver problems with external disturbances and perturbed disturbances. Yuan et al. [6] proposed a decoupling control algorithm based on a robust adaptive method to solve the problem of spacecraft attitude control with disturbance and torque of inertia uncertainty. Tong and Li [7] proposed a static output feedback controller using multiobjective synthesis technology and studied the robust stability and disturbance suppression of spacecraft with parameter variation and control input saturation constraints. However, this method does not consider the uncertainties of the attachment frequency and rigid-flexible coupling matrix. Yang et al. [8] studied the orbit robust control of low earth orbit spacecraft under the condition of parameter variation during orbit transfer during rendezvous. Liang et al. [9] planned the angular trajectory of the satellite’s attitude maneuver around the Euler axis and designed a variable structure control law based on the error quaternion and error angular velocity between the actual and planned positions of the satellite. Based on the time-varying autoregressive sliding average model and the good local function fitting ability of the wavelet basis function, Lei et al. [10] used the Mexican cap wavelet function as the spatial base of the time-varying coefficients of the TARMA model and constructed a time-varying autoregressive sliding average model of the functional series based on the wavelet function. And decoupling estimation of time-varying coefficients is achieved.In this paper, from the angle of system angular torque, the angular torque of the system is introduced into the control system, and the angular torque and attitude of the system are jointly controlled to achieve a good balance between the attitude and angular torque of the system and ensure that the system is stable in a certain equilibrium attitude, and the angular torque of the system will not be saturated due to the change of the mass characteristics of the system. The joint control of attitude and angular torque is referred to as attitude control/angular torque management (ACMM).


## 2. Problem Formulation

### 2.1. Simplified System Model

The replenishment vehicle operates the replenished vehicle with a manipulator. Taking the configuration transformation process of the resupplied aircraft from “I” configuration to “L” configuration as an example, the model is established.Assuming that the solar panels are locked and the slow variables are ignored, the system dynamics model of the configuration transformation process is as follows:(1)Istω˙s+I˙stωs+RbIω˙bI+ωs×Istωs+ωs×RbIωbI+∑i=1nF¨siηai=Ts.(2)η¨ai+2ζaiΩaiη˙ai+Ωai2ηai+FtiTX¨+FsiTω˙s+FaiTω˙ai+FIiTω˙bI=0.Among them, formula (1) is the attitude dynamics equation of the central rigid body, and formula (2) is the vibration equation of the solar panel. Ist is the expression of the instantaneous torque of inertia of the combination in the body coordinate system (Fs). The specific expression is as follows:(3)Ist=Ibb+IIr−mrbIs×2,where Ibb is the torque of inertia of the replenishment vehicle relative to its body coordinate system. IIr is the representation of the torque of inertia of the replenished aircraft in its reference coordinate system. Let the torque of inertia of the replenished vehicle can be expressed as III in the body coordinate system, and the transformation matrix between the body coordinate system and the reference coordinate system is LIrIb, and then,(4)IIr=LIrIbIIILIrIbT.ωs is the component representation of the angular rate of rotation of the system relative to the inertial frame in the body assembly coordinate system. rbIb is the component representation of the vector from the centroid of the replenishment vehicle to the centroid of the replenished vehicle in the body coordinate system of the assembly. m is the reduced mass of the two-body system. If the mass of the replenishment vehicle and the replenished vehicle is mb and mI, respectively, the reduced mass can be expressed as follows:(5)m=mbmImb+mI.ωbI is the rotational angular velocity of the replenished vehicle relative to the replenished vehicle, that is, the relative attitude angular velocity, which is described in the body coordinate system of the assembly. RbI is the motion coupling coefficient of the replenishment vehicle and the replenished vehicle, which is expressed in the body assembly coordinate system. FIi is the coupling coefficient between the solar panel and the rotational motion of the replenished vehicle, which is expressed in the system of the replenished vehicle.FaiFsi are the coupling coefficient matrix of solar panel rotation and satellite rotation. η represents the displacement caused by the force. Fti is the principal vector array of system forces.When there is no confusion, the subscripts is omitted to identify the system parameter, and the coordinate system of the variable is clearly marked. Let ωAB be the angular velocity vector of a coordinate system relative to the B coordinate system. When B coordinate system is the inertial system, B is omitted. For example, XAB is the component representation of XA in the B coordinate system. The torque of inertia Is, the angular velocity of inertia ωs, and control torque Τs in the centroid coordinate system of the system are abbreviated as Ib, ωb, and Τb, respectively.In the process of configuration transformation, the dynamic equation and environmental torque model of CMG are consistent with that of long-term normal on-orbit flight. The final system model for controller design is as follows [11]:Attitude dynamics:(6)Ibtω˙b=Tcb+Tgvrob+Tdlb+Tjb+Tgb+Tdb,where ωb is the angular velocity of inertia; Ib is the torque of inertia; Tcb is the output torque of CMG; Tgvrob is the gyro coupling torque; Tdlb is the perturbation torque of the supplied vehicle; Tgb is the gravity gradient torque; Tdb is the atmospheric disturbance torque.LetFoo0x0y0z0 be the orbital coordinate system, and Fbobxbybzb be the body coordinate system. The attitude kinematics model can be expressed as follows:(7)ωxbωybωzb=φ˙cθ−ψ˙cφsθθ˙+ψ˙sφφ˙sθ+ψ˙cφcθ−ω0sψcθ+sφsθcψω0cψcφ−ω0sφcθcψ−sψsθ,where ψ, φ, θ are the yaw, roll, and pitch angles of the aircraft, respectively. “s” and “c” represent “sin” and “cos” functions, respectively. ωxb, ωyb, and ωzb represent the triaxial component of the absolute angular velocity ωb of the combination in the coordinate system Fb, and ω0 is the orbital angular velocity. It is assumed that ω0 is constant when the assembly runs in a circular orbit.The CMG kinetic model is as follows:(8)−h˙cb−ωb×hcb=Tcb,where hcb is the sum of the absolute angular torque of each CMG to its own center of mass, which is described in the body coordinate system of the combination; Tcb is the attitude control torque of CMG.When the attitude of the system satisfies the following equation, all kinds of torques can be balanced.(9)−ωb∗×Ibtωb∗+Tdlb+Tgb+Tdb=0.If the controller is designed to make the attitude of the combination track the TEA under this condition, the configuration transformation operation can be completed without unloading the angular torque of CMG. However, due to the complexity of configuration transformation, TEA is difficult to solve. And because of the rapid change of TEA in this process, large space structures are generally not suitable for rapid maneuver. So, the attitude control strategy of tracking dynamic TEA in the configuration transformation process is not feasible.The time of configuration transformation is shorter than that of the normal flight in orbit, and the angular torque accumulation caused by aerodynamic torque is much smaller than that caused by gravitational gradient torque and perturbation torque. Although the perturbation torque is one of the main reasons for CMG angular torque accumulation, the perturbation angular torque in different stages of the transfer process can cancel each other to a great extent after the path and velocity of the manipulator are reasonably planned. Therefore, in the preliminary design of the controller, the influence of aerodynamic torque and perturbation torque can be ignored temporarily. In this case, only the gravitational gradient torque and the orbital gyroscopic torque are considered, and the attitude satisfying the following form during configuration transformation is defined as “narrow sense TEA.”(10)−ωb∗×Ibtωb∗+Tgb=0.In the process of configuration transformation, the purpose of ACMM controller design is to make the system track the narrow sense TEA. The physical meaning of TEA in the orbit system is obvious when only considering the gravitational gradient torque and the orbital gyroscopic torque. Therefore, the controller design of this part is carried out in the orbit system, which is not only convenient for the derivation of feedback linearization control law but also can further study the characteristics of TEA in the orbit system [12].


### 2.2. Mechanical Model

#### 2.2.1. Attitude Dynamics Equation

The derivative of coordinate transformation matrix from body coordinate system to orbit system is as follows:(11)L˙ob=ΩbooLob,where Ωboo=ωbo−η×.In the process of configuration transformation, the torque of inertia of the system changes constantly, but without confusion, the time mark is omitted. The relationship between the torque of inertia of the system in the body coordinate system and the orbit system is as follows:(12)Io=LobIbLbo,(13)I˙bo=ΩbooIbo−IboΩboo.The component of the absolute angular velocity of the combination in the orbit system is expressed as follows:(14)ωo=Lobωb.By deriving the time from equation (14) and using equation (10), the system dynamic equation described in the orbit system is obtained.(15)ω˙o=Io−1ωo×Ioωo+3ω02R^o×IoR^o+Tco+Tdo.


#### 2.2.2. Attitude Kinematics Equation

Note thatΘob=φobθobψobT is the Euler angle from the body coordinate system to the orbit system. In the case of no confusion, omit the lower corner mark indicating the direction of rotation. Using 3-1-2 transformation order, the kinematics equation of the system in the orbit system is as follows:(16)Θ˙ob=Rωoωo+η,where Rωo=−1/cosφcosφcosθ0cosφsinθsinφsinθcosφ−sinφsinθ−sinθ0cosθ.


#### 2.2.3. CMG Kinetic Equation

The dynamic equation of CMG angular torque in the orbit system is as follows:(17)h˙co=−ωo×hco−Tco.


#### 2.2.4. ACMM System Design Model

Equations (15)–(17) are written as state equations, and the system model for controller design is obtained as(18)Θ˙obω˙oh˙co=Rωoωo+ηΙo−1ωo×Ιoωo+3ω02R^o×ΙoR^o−ωo×hco+03Ιo−1−E3Tco+03Ιo−103Tdo.


## 2.1. Simplified System Model

The replenishment vehicle operates the replenished vehicle with a manipulator. Taking the configuration transformation process of the resupplied aircraft from “I” configuration to “L” configuration as an example, the model is established.Assuming that the solar panels are locked and the slow variables are ignored, the system dynamics model of the configuration transformation process is as follows:(1)Istω˙s+I˙stωs+RbIω˙bI+ωs×Istωs+ωs×RbIωbI+∑i=1nF¨siηai=Ts.(2)η¨ai+2ζaiΩaiη˙ai+Ωai2ηai+FtiTX¨+FsiTω˙s+FaiTω˙ai+FIiTω˙bI=0.Among them, formula (1) is the attitude dynamics equation of the central rigid body, and formula (2) is the vibration equation of the solar panel. Ist is the expression of the instantaneous torque of inertia of the combination in the body coordinate system (Fs). The specific expression is as follows:(3)Ist=Ibb+IIr−mrbIs×2,where Ibb is the torque of inertia of the replenishment vehicle relative to its body coordinate system. IIr is the representation of the torque of inertia of the replenished aircraft in its reference coordinate system. Let the torque of inertia of the replenished vehicle can be expressed as III in the body coordinate system, and the transformation matrix between the body coordinate system and the reference coordinate system is LIrIb, and then,(4)IIr=LIrIbIIILIrIbT.ωs is the component representation of the angular rate of rotation of the system relative to the inertial frame in the body assembly coordinate system. rbIb is the component representation of the vector from the centroid of the replenishment vehicle to the centroid of the replenished vehicle in the body coordinate system of the assembly. m is the reduced mass of the two-body system. If the mass of the replenishment vehicle and the replenished vehicle is mb and mI, respectively, the reduced mass can be expressed as follows:(5)m=mbmImb+mI.ωbI is the rotational angular velocity of the replenished vehicle relative to the replenished vehicle, that is, the relative attitude angular velocity, which is described in the body coordinate system of the assembly. RbI is the motion coupling coefficient of the replenishment vehicle and the replenished vehicle, which is expressed in the body assembly coordinate system. FIi is the coupling coefficient between the solar panel and the rotational motion of the replenished vehicle, which is expressed in the system of the replenished vehicle.FaiFsi are the coupling coefficient matrix of solar panel rotation and satellite rotation. η represents the displacement caused by the force. Fti is the principal vector array of system forces.When there is no confusion, the subscripts is omitted to identify the system parameter, and the coordinate system of the variable is clearly marked. Let ωAB be the angular velocity vector of a coordinate system relative to the B coordinate system. When B coordinate system is the inertial system, B is omitted. For example, XAB is the component representation of XA in the B coordinate system. The torque of inertia Is, the angular velocity of inertia ωs, and control torque Τs in the centroid coordinate system of the system are abbreviated as Ib, ωb, and Τb, respectively.In the process of configuration transformation, the dynamic equation and environmental torque model of CMG are consistent with that of long-term normal on-orbit flight. The final system model for controller design is as follows [11]:Attitude dynamics:(6)Ibtω˙b=Tcb+Tgvrob+Tdlb+Tjb+Tgb+Tdb,where ωb is the angular velocity of inertia; Ib is the torque of inertia; Tcb is the output torque of CMG; Tgvrob is the gyro coupling torque; Tdlb is the perturbation torque of the supplied vehicle; Tgb is the gravity gradient torque; Tdb is the atmospheric disturbance torque.LetFoo0x0y0z0 be the orbital coordinate system, and Fbobxbybzb be the body coordinate system. The attitude kinematics model can be expressed as follows:(7)ωxbωybωzb=φ˙cθ−ψ˙cφsθθ˙+ψ˙sφφ˙sθ+ψ˙cφcθ−ω0sψcθ+sφsθcψω0cψcφ−ω0sφcθcψ−sψsθ,where ψ, φ, θ are the yaw, roll, and pitch angles of the aircraft, respectively. “s” and “c” represent “sin” and “cos” functions, respectively. ωxb, ωyb, and ωzb represent the triaxial component of the absolute angular velocity ωb of the combination in the coordinate system Fb, and ω0 is the orbital angular velocity. It is assumed that ω0 is constant when the assembly runs in a circular orbit.The CMG kinetic model is as follows:(8)−h˙cb−ωb×hcb=Tcb,where hcb is the sum of the absolute angular torque of each CMG to its own center of mass, which is described in the body coordinate system of the combination; Tcb is the attitude control torque of CMG.When the attitude of the system satisfies the following equation, all kinds of torques can be balanced.(9)−ωb∗×Ibtωb∗+Tdlb+Tgb+Tdb=0.If the controller is designed to make the attitude of the combination track the TEA under this condition, the configuration transformation operation can be completed without unloading the angular torque of CMG. However, due to the complexity of configuration transformation, TEA is difficult to solve. And because of the rapid change of TEA in this process, large space structures are generally not suitable for rapid maneuver. So, the attitude control strategy of tracking dynamic TEA in the configuration transformation process is not feasible.The time of configuration transformation is shorter than that of the normal flight in orbit, and the angular torque accumulation caused by aerodynamic torque is much smaller than that caused by gravitational gradient torque and perturbation torque. Although the perturbation torque is one of the main reasons for CMG angular torque accumulation, the perturbation angular torque in different stages of the transfer process can cancel each other to a great extent after the path and velocity of the manipulator are reasonably planned. Therefore, in the preliminary design of the controller, the influence of aerodynamic torque and perturbation torque can be ignored temporarily. In this case, only the gravitational gradient torque and the orbital gyroscopic torque are considered, and the attitude satisfying the following form during configuration transformation is defined as “narrow sense TEA.”(10)−ωb∗×Ibtωb∗+Tgb=0.In the process of configuration transformation, the purpose of ACMM controller design is to make the system track the narrow sense TEA. The physical meaning of TEA in the orbit system is obvious when only considering the gravitational gradient torque and the orbital gyroscopic torque. Therefore, the controller design of this part is carried out in the orbit system, which is not only convenient for the derivation of feedback linearization control law but also can further study the characteristics of TEA in the orbit system [12].


## 2.2. Mechanical Model

### 2.2.1. Attitude Dynamics Equation

The derivative of coordinate transformation matrix from body coordinate system to orbit system is as follows:(11)L˙ob=ΩbooLob,where Ωboo=ωbo−η×.In the process of configuration transformation, the torque of inertia of the system changes constantly, but without confusion, the time mark is omitted. The relationship between the torque of inertia of the system in the body coordinate system and the orbit system is as follows:(12)Io=LobIbLbo,(13)I˙bo=ΩbooIbo−IboΩboo.The component of the absolute angular velocity of the combination in the orbit system is expressed as follows:(14)ωo=Lobωb.By deriving the time from equation (14) and using equation (10), the system dynamic equation described in the orbit system is obtained.(15)ω˙o=Io−1ωo×Ioωo+3ω02R^o×IoR^o+Tco+Tdo.


### 2.2.2. Attitude Kinematics Equation

Note thatΘob=φobθobψobT is the Euler angle from the body coordinate system to the orbit system. In the case of no confusion, omit the lower corner mark indicating the direction of rotation. Using 3-1-2 transformation order, the kinematics equation of the system in the orbit system is as follows:(16)Θ˙ob=Rωoωo+η,where Rωo=−1/cosφcosφcosθ0cosφsinθsinφsinθcosφ−sinφsinθ−sinθ0cosθ.


### 2.2.3. CMG Kinetic Equation

The dynamic equation of CMG angular torque in the orbit system is as follows:(17)h˙co=−ωo×hco−Tco.


### 2.2.4. ACMM System Design Model

Equations (15)–(17) are written as state equations, and the system model for controller design is obtained as(18)Θ˙obω˙oh˙co=Rωoωo+ηΙo−1ωo×Ιoωo+3ω02R^o×ΙoR^o−ωo×hco+03Ιo−1−E3Tco+03Ιo−103Tdo.


## 2.2.1. Attitude Dynamics Equation

The derivative of coordinate transformation matrix from body coordinate system to orbit system is as follows:(11)L˙ob=ΩbooLob,where Ωboo=ωbo−η×.In the process of configuration transformation, the torque of inertia of the system changes constantly, but without confusion, the time mark is omitted. The relationship between the torque of inertia of the system in the body coordinate system and the orbit system is as follows:(12)Io=LobIbLbo,(13)I˙bo=ΩbooIbo−IboΩboo.The component of the absolute angular velocity of the combination in the orbit system is expressed as follows:(14)ωo=Lobωb.By deriving the time from equation (14) and using equation (10), the system dynamic equation described in the orbit system is obtained.(15)ω˙o=Io−1ωo×Ioωo+3ω02R^o×IoR^o+Tco+Tdo.


## 2.2.2. Attitude Kinematics Equation

Note thatΘob=φobθobψobT is the Euler angle from the body coordinate system to the orbit system. In the case of no confusion, omit the lower corner mark indicating the direction of rotation. Using 3-1-2 transformation order, the kinematics equation of the system in the orbit system is as follows:(16)Θ˙ob=Rωoωo+η,where Rωo=−1/cosφcosφcosθ0cosφsinθsinφsinθcosφ−sinφsinθ−sinθ0cosθ.


## 2.2.3. CMG Kinetic Equation

The dynamic equation of CMG angular torque in the orbit system is as follows:(17)h˙co=−ωo×hco−Tco.


## 2.2.4. ACMM System Design Model

Equations (15)–(17) are written as state equations, and the system model for controller design is obtained as(18)Θ˙obω˙oh˙co=Rωoωo+ηΙo−1ωo×Ιoωo+3ω02R^o×ΙoR^o−ωo×hco+03Ιo−1−E3Tco+03Ιo−103Tdo.


## 3. Control System Design based on Adaptive Feedback Linearization

In the process of configuration transformation, the large change of torque of inertia may make the narrow sense TEA of the system deviate greatly from the earth-oriented attitude. In order to avoid the error caused by the partial linearization of the system model, the feedback linearization method of the nonlinear system is used to design the controller to realize the tracking of the narrow sense TEA. In the process of configuration transformation, the large change of torque of inertia may make the narrow sense TEA of the system deviate greatly from the earth-oriented attitude. In order to avoid the error caused by the partial linearization of the system model, the feedback linearization method of the nonlinear system is used to design the controller to realize the tracking of the narrow sense TEA in the process of configuration transformation.The core idea of feedback linearization is to make the nonlinear system realize accurate linearization of state or input/output under certain conditions through appropriate nonlinear state transformation and feedback transformation, so as to transform the nonlinear system design problem into a linear system design problem. It is different from the traditional method of local linearization using Taylor expansion [13], which does not ignore any nonlinear term in the linearization process, so this method is not only accurate but also holistic; that is, linearization is applicable to the whole region of the transformation. The limitation of this method is that it needs accurate information of the system. When the parameters are uncertain, it is unable to carry out accurate feedback linearization.In order to compensate for the dependence of the feedback linearization method on system parameters, an adaptive ACMM controller with online parameter identification is designed. The controller consists of two parts: online parameter identification loop and feedback linearization control loop. The structure diagram is shown in Figure1 [14].Figure 1
Structure diagram of adaptive feedback linearization controller.


### 3.1. Feedback Linearization Analysis

Under the assumption of a small angle, the coordinate transformation matrix fromFb to Fo can be approximately expressed as(19)Lob=I−Θob×,where E3 is the 3×3 unit matrix and Θob× is the antisymmetric matrix of Θob, Θob×=0−ψθψ0−φ−θφ0.The results are as follows:(20)Io=E3−Θob×IbE3+Θob×≅Ib−Θob×Ib+IbΘob×,Io−1≈Ib−1−Θob×Ib−1+Ib−1Θob×.Taking the system state variable asx=ΘobωohcoT, the state equation is written as a general nonlinear system(21)x˙≜fx+Gxu,where x∈R9,u∈R3,(22)fx=Rωoωo+ηIo−1ωo×Ioωo+3ω02R^o×IoR^o−ωo×hcoT,(23)Gx=03Io−1−E3T,u=Tco.The whole state of the system can be measured, and whenIo≠0, fx, Gx is continuous and smooth.Substituting the approximate expression equation (15) of Io−1 into the expression (23) of Gx, denote Gx=g1xg2xg3x. Definition:(24)G0x=spang1xg2xg3x,G1x=spang1xg2xg3xadfg1xadfg2xadfg3x,…Gix=spanadfkgjx:0≤k≤i,1≤j≤3,i=1,…,8,where span represents vector expansion; that is, G0x represents subspace formed by a linear combination of g1xg2xg3x, adfkgx=f,adfk−1gx. Where [] is Lie brackets, f,gx=∂g/∂xfx−∂f/∂xgx, hc0o=0hcyo0, x0=0ηhc0oT can be used to test.When1≤i≤8, the distribution Gix is a constant dimension in the neighborhood of x0.When1≤i≤7, the distribution Gix is involutivity.The distributionGix has dimension 9.According to the exact feedback linearization theorem, the system can be linearized by exact feedback. Therefore, three output variables with total relative order 9 can be used to define the state transformation, and the system can be transformed into a standard canonical form. Since the total relative order of the system is equal to the state dimension of the system, all states can be observed through the input-output relationship.


### 3.2. State Transition

The state transformation is defined by the Lie derivative of the output variable. For the ACMM system, in order to avoid angular torque accumulation effectively under the premise of attitude stability, CMG angular torque and attitude stability information should be taken as output variables at the same time.(25)Hbo=hco+Ioωo.The second-order Lie derivative, the second-order Lie derivative, and the third-order Lie derivative, respectively, for the triaxial components of equationHbo are obtained as(26)Lf2Hb1ox=−ω02Hb1o+3ω02Iyzo,(27)Lf2Hb2ox=3ω02Iyzo,(28)Lf2Hb3ox=−ω03Hb1o−3ω02Iyzo.It can be seen from equation (12) that I˙o contains angular velocity information, and the control torque information will appear after further derivation. Therefore, according to equations (26)–(28), the third derivative of Hb1o, Hb2o and the fourth derivative of Hb3o will appear control torque information, which can be used as output variables.In order to meet the relative order requirements of the system, the magnitude after derivation is unified, the accuracy of numerical calculation is improved, and the output variables are selected as follows:(29)y1=Hb3ox,y2=ω0Hb2ox,y3=3ω03Ixyox.The nonlinear transformation of output variable definition is as follows:(30)z=Φx,where zT=z1Tz2Tz3T, xT=ΘobTωoThcoTT.The specific form of nonlinear transformation is as follows:(31)z11=Hb3ox,z12=LfHb3ox=−ω02Hb1o,z13=Lf2Hb3ox=−3ω03Iyzo−ω02Hb3o,z14=Lf3Hb3ox=−3ω03Iyo−Izoωx2+Ixyoωy2+ω0−Ixzoωz2+ω03Hb1o,(32)z21=ω0Hb2ox,z22=Lfω0Hb2ox=−3ω02Ixzo,z23=Lf2ω0Hb2ox=3ω03−Ixyoωx2+Izo−Ixoωy2+ω0+Iyzoωzo,(33)z31=3ω03Ixyox,z32=Lf3ω03IIxyox=3ω03−Ixzoωxo+Iyzoωyo+ω0−Ixo−Iyoωzo.The new equation of state is as follows:(34)z˙11=z12,z˙12=z13,z˙13=z14,z˙21=z22,z˙22=z23,z˙31=z32,(35)z˙14z˙23z˙32=Lf4y1Lf3y2Lf2y3+Lg1Lf3y1Lg2Lf3y1Lg3Lf3y1Lg1Lf2y2Lg2Lf2y2Lg3Lf2y2Lg1Lfy3Lg2Lfy3Lg3Lfy3u1ou2ou3o.Letz∗=z14z23z32T, (35) be written as follows:(36)z˙∗≜f∗x+Exuo.The definitions off∗x and decoupling matrix Ex are obvious.(37)Ex=−3ω03S1oIo−1,f∗x=c+3ω02S1oωo×η+S2oωo−η−ωo−η×S1oωo−η+S1oIo−1−ωo×Ioωo+3ω02R^o×IoR^o,where S1o=Iyo−IzoIxyo−Ixzo−IxyoIzo−IxoIyzoIxzo−IyzoIxo−Iyo， S2o=Iyo−Izo000Izo−Ixo000Ixo−Iyo, c=−ω02z100.


### 3.3. Input Transformation

WhenEx is reversible, the following input transformation is adopted for the system (25):(38)uo=E−1xv1−Lf4Hb3oxv2−Lf3ω0Hb2oxv2−Lf2ω03Hb2ox.System (25) will be equivalent to(39)z˙∗=v,where v=v1v2v3T.The system (39) has a linear input-output relationship. When its output is expected to track the desired trajectory zd without error, the corresponding linear control law can be designed.(40)v=Kzd−z+z˙d∗,where K is the control gain matrix in the form of (40) and z˙d∗ is the expected trajectory change rate.(41)K=k11k12k13k14000000000k21k22k23000000000k31k32.The characteristic equation of an equivalent linear system is as follows:(42)λ=s4+k14s3+k13s2+k12s+k11s3+k23s2+k22s+k21s2+k32s+k31.According to equations (36) and (41), the nonlinear control law for the linearized system (39) is obtained.(43)uo=E−1x−f∗x+v,(44)uo=Ioωo×η−ωo×Ioωo+3ω02R^o×IoRo−IoS1o−1ωo−ηS1oωo−η−S2oωo−η−13ω03IoS1o−1z˙d′−ε+Kzd−z,(45)uo=Lbouo.Equations (44) and (45) are the nonlinear control laws of the ACMM system in the track system. In this control law, it is assumed that the inertia characteristics of the system can be fully identified, and only the desired trajectory zd and the feedback gain matrix K need to be designed.


### 3.4. System Expected Trajectory

The stable working state of the ACMM system is a narrow sense TEA when only considering the gravitational gradient torque. At this point,Ixyo=Ixzo=Iyzo=0,Ixo=Ixp,Iyo=Iyo,Izo=Izp, and hcxo=0,hcyo=0. The total angular torque of the system at TEA is(46)Hbo=0hco0+Iip000Ijp000Ikp0−ω00=0hco−Ijpω00,where i≠j,j≠k,i≠k,Iip,Ijp,Ikp are the value of inertia along the three principal axes of inertia.After the total angular torque of the system at TEA is obtained, the desired trajectoryzd of the system can be determined. Although TEA does not require hcyo and can take any value within the capacity range of CMG, in order to simplify the design of the controller, the instruction value of hcyo is zero. The narrow sense of TEA requires the inertial principal axis to point along the orbit coordinate axis, but in the case of minimum attitude maneuver, the inertial principal axis in the y direction is generally pointed to the normal direction of the orbit; that is, Hbo=0−Iypω00.From the above analysis, we can get the expected trajectoryzd as follows:(47)zd=0000−Iypω00000T.


### 3.5. Equivalent Linear System Design

The closed-loop characteristic of the equivalent linear system is designed by the pole assignment method. In order to make the system have a certain stability margin and convergence speed, the closed-loop poles are placed in the sector area with an angle of±45° between the left side of the S plane S=−0.5ω0 and the real axis. Select 9 closed-loop poles as −0.5ω0, −0.707±0.1jω0, −1±0.5jω0, −1.414±0.5jω0, −2±jω0, and the specific distribution is shown in Figure 2.Figure 2
Pole assignment diagram of the equivalent linear system.For the fourth-order SISO system corresponding toz1, the closed-loop poles are assigned at −0.707±0.1jω0 and −1.414±0.5jω0. The corresponding fourth-order characteristic equation is as follows:(48)λ1s=s4+4.24ω0s3+6.75ω02s2+1.77ω03s+2.75ω04.Therefore,k11=2.75ω04,k12=1.77ω03,k13=6.75ω02,k14=4.24ω0.For the third-order single input single output system ofz2 choosing −0.5ω0 and −2±jω0 as its closed-loop poles, the corresponding feedback gain can be obtained in the same way k21=2.5ω03,k22=7.0ω02,k23=4.5ω0.For the second-order single input single output system ofz3, the closed-loop pole is −1±0.5jω0, and the feedback gain is calculated as k31=1.25ω02,k32=2.0ω0.


## 3.1. Feedback Linearization Analysis

Under the assumption of a small angle, the coordinate transformation matrix fromFb to Fo can be approximately expressed as(19)Lob=I−Θob×,where E3 is the 3×3 unit matrix and Θob× is the antisymmetric matrix of Θob, Θob×=0−ψθψ0−φ−θφ0.The results are as follows:(20)Io=E3−Θob×IbE3+Θob×≅Ib−Θob×Ib+IbΘob×,Io−1≈Ib−1−Θob×Ib−1+Ib−1Θob×.Taking the system state variable asx=ΘobωohcoT, the state equation is written as a general nonlinear system(21)x˙≜fx+Gxu,where x∈R9,u∈R3,(22)fx=Rωoωo+ηIo−1ωo×Ioωo+3ω02R^o×IoR^o−ωo×hcoT,(23)Gx=03Io−1−E3T,u=Tco.The whole state of the system can be measured, and whenIo≠0, fx, Gx is continuous and smooth.Substituting the approximate expression equation (15) of Io−1 into the expression (23) of Gx, denote Gx=g1xg2xg3x. Definition:(24)G0x=spang1xg2xg3x,G1x=spang1xg2xg3xadfg1xadfg2xadfg3x,…Gix=spanadfkgjx:0≤k≤i,1≤j≤3,i=1,…,8,where span represents vector expansion; that is, G0x represents subspace formed by a linear combination of g1xg2xg3x, adfkgx=f,adfk−1gx. Where [] is Lie brackets, f,gx=∂g/∂xfx−∂f/∂xgx, hc0o=0hcyo0, x0=0ηhc0oT can be used to test.When1≤i≤8, the distribution Gix is a constant dimension in the neighborhood of x0.When1≤i≤7, the distribution Gix is involutivity.The distributionGix has dimension 9.According to the exact feedback linearization theorem, the system can be linearized by exact feedback. Therefore, three output variables with total relative order 9 can be used to define the state transformation, and the system can be transformed into a standard canonical form. Since the total relative order of the system is equal to the state dimension of the system, all states can be observed through the input-output relationship.


## 3.2. State Transition

The state transformation is defined by the Lie derivative of the output variable. For the ACMM system, in order to avoid angular torque accumulation effectively under the premise of attitude stability, CMG angular torque and attitude stability information should be taken as output variables at the same time.(25)Hbo=hco+Ioωo.The second-order Lie derivative, the second-order Lie derivative, and the third-order Lie derivative, respectively, for the triaxial components of equationHbo are obtained as(26)Lf2Hb1ox=−ω02Hb1o+3ω02Iyzo,(27)Lf2Hb2ox=3ω02Iyzo,(28)Lf2Hb3ox=−ω03Hb1o−3ω02Iyzo.It can be seen from equation (12) that I˙o contains angular velocity information, and the control torque information will appear after further derivation. Therefore, according to equations (26)–(28), the third derivative of Hb1o, Hb2o and the fourth derivative of Hb3o will appear control torque information, which can be used as output variables.In order to meet the relative order requirements of the system, the magnitude after derivation is unified, the accuracy of numerical calculation is improved, and the output variables are selected as follows:(29)y1=Hb3ox,y2=ω0Hb2ox,y3=3ω03Ixyox.The nonlinear transformation of output variable definition is as follows:(30)z=Φx,where zT=z1Tz2Tz3T, xT=ΘobTωoThcoTT.The specific form of nonlinear transformation is as follows:(31)z11=Hb3ox,z12=LfHb3ox=−ω02Hb1o,z13=Lf2Hb3ox=−3ω03Iyzo−ω02Hb3o,z14=Lf3Hb3ox=−3ω03Iyo−Izoωx2+Ixyoωy2+ω0−Ixzoωz2+ω03Hb1o,(32)z21=ω0Hb2ox,z22=Lfω0Hb2ox=−3ω02Ixzo,z23=Lf2ω0Hb2ox=3ω03−Ixyoωx2+Izo−Ixoωy2+ω0+Iyzoωzo,(33)z31=3ω03Ixyox,z32=Lf3ω03IIxyox=3ω03−Ixzoωxo+Iyzoωyo+ω0−Ixo−Iyoωzo.The new equation of state is as follows:(34)z˙11=z12,z˙12=z13,z˙13=z14,z˙21=z22,z˙22=z23,z˙31=z32,(35)z˙14z˙23z˙32=Lf4y1Lf3y2Lf2y3+Lg1Lf3y1Lg2Lf3y1Lg3Lf3y1Lg1Lf2y2Lg2Lf2y2Lg3Lf2y2Lg1Lfy3Lg2Lfy3Lg3Lfy3u1ou2ou3o.Letz∗=z14z23z32T, (35) be written as follows:(36)z˙∗≜f∗x+Exuo.The definitions off∗x and decoupling matrix Ex are obvious.(37)Ex=−3ω03S1oIo−1,f∗x=c+3ω02S1oωo×η+S2oωo−η−ωo−η×S1oωo−η+S1oIo−1−ωo×Ioωo+3ω02R^o×IoR^o,where S1o=Iyo−IzoIxyo−Ixzo−IxyoIzo−IxoIyzoIxzo−IyzoIxo−Iyo， S2o=Iyo−Izo000Izo−Ixo000Ixo−Iyo, c=−ω02z100.


## 3.3. Input Transformation

WhenEx is reversible, the following input transformation is adopted for the system (25):(38)uo=E−1xv1−Lf4Hb3oxv2−Lf3ω0Hb2oxv2−Lf2ω03Hb2ox.System (25) will be equivalent to(39)z˙∗=v,where v=v1v2v3T.The system (39) has a linear input-output relationship. When its output is expected to track the desired trajectory zd without error, the corresponding linear control law can be designed.(40)v=Kzd−z+z˙d∗,where K is the control gain matrix in the form of (40) and z˙d∗ is the expected trajectory change rate.(41)K=k11k12k13k14000000000k21k22k23000000000k31k32.The characteristic equation of an equivalent linear system is as follows:(42)λ=s4+k14s3+k13s2+k12s+k11s3+k23s2+k22s+k21s2+k32s+k31.According to equations (36) and (41), the nonlinear control law for the linearized system (39) is obtained.(43)uo=E−1x−f∗x+v,(44)uo=Ioωo×η−ωo×Ioωo+3ω02R^o×IoRo−IoS1o−1ωo−ηS1oωo−η−S2oωo−η−13ω03IoS1o−1z˙d′−ε+Kzd−z,(45)uo=Lbouo.Equations (44) and (45) are the nonlinear control laws of the ACMM system in the track system. In this control law, it is assumed that the inertia characteristics of the system can be fully identified, and only the desired trajectory zd and the feedback gain matrix K need to be designed.


## 3.4. System Expected Trajectory

The stable working state of the ACMM system is a narrow sense TEA when only considering the gravitational gradient torque. At this point,Ixyo=Ixzo=Iyzo=0,Ixo=Ixp,Iyo=Iyo,Izo=Izp, and hcxo=0,hcyo=0. The total angular torque of the system at TEA is(46)Hbo=0hco0+Iip000Ijp000Ikp0−ω00=0hco−Ijpω00,where i≠j,j≠k,i≠k,Iip,Ijp,Ikp are the value of inertia along the three principal axes of inertia.After the total angular torque of the system at TEA is obtained, the desired trajectoryzd of the system can be determined. Although TEA does not require hcyo and can take any value within the capacity range of CMG, in order to simplify the design of the controller, the instruction value of hcyo is zero. The narrow sense of TEA requires the inertial principal axis to point along the orbit coordinate axis, but in the case of minimum attitude maneuver, the inertial principal axis in the y direction is generally pointed to the normal direction of the orbit; that is, Hbo=0−Iypω00.From the above analysis, we can get the expected trajectoryzd as follows:(47)zd=0000−Iypω00000T.


## 3.5. Equivalent Linear System Design

The closed-loop characteristic of the equivalent linear system is designed by the pole assignment method. In order to make the system have a certain stability margin and convergence speed, the closed-loop poles are placed in the sector area with an angle of±45° between the left side of the S plane S=−0.5ω0 and the real axis. Select 9 closed-loop poles as −0.5ω0, −0.707±0.1jω0, −1±0.5jω0, −1.414±0.5jω0, −2±jω0, and the specific distribution is shown in Figure 2.Figure 2
Pole assignment diagram of the equivalent linear system.For the fourth-order SISO system corresponding toz1, the closed-loop poles are assigned at −0.707±0.1jω0 and −1.414±0.5jω0. The corresponding fourth-order characteristic equation is as follows:(48)λ1s=s4+4.24ω0s3+6.75ω02s2+1.77ω03s+2.75ω04.Therefore,k11=2.75ω04,k12=1.77ω03,k13=6.75ω02,k14=4.24ω0.For the third-order single input single output system ofz2 choosing −0.5ω0 and −2±jω0 as its closed-loop poles, the corresponding feedback gain can be obtained in the same way k21=2.5ω03,k22=7.0ω02,k23=4.5ω0.For the second-order single input single output system ofz3, the closed-loop pole is −1±0.5jω0, and the feedback gain is calculated as k31=1.25ω02,k32=2.0ω0.


## 4. Online Parameter Identification

It can be seen from the feedback control law (39) that the ACMM controller needs not only the angular rate and angular torque information of the combined body but also its torque of inertia information, which is provided by the online parameter identification unit, and its identification accuracy directly affects the control performance of the adaptive system. In this section, the control torque information and angular velocity information of the combination are used to identify the closed-loop control of the combination without additional excitation. The algorithm uses the least square method with a bounded gain forgetting factor, which can track the real parameters without continuous excitation.The standard form of the linear least squares [15] problem is Φx=y+ε or written as Φx≅y, where y is the measured vector, ε is the measured noise vector, x is the parameter to be identified, Φ is composed of known variables and parameters, and x^ is the solution of the least-squares algorithm, which minimizes the sum of squares of error Φx^-y.Firstly, the least square estimation model of the parameters to be identified is constructed according to the dynamic equation of the system, and the vector composed of six independent elements of the torque of inertia is taken as the unknown vectori=IxoIyoIzoIxyoIxzoIyzoT.(49)Ιoωo=D1ω˙oi,ωo×Ιoωo=D2ωoi,where(50)D1ω˙=ω˙xo000ω˙zoω˙yo0ω˙yo0ω˙zo0ω˙xo00ω˙zoω˙yoω˙xo0,D2ω˙=0−ωyoωzoωyoωzoωyo2−ωzo2ωxoωyo−ωxoωzoωxoωzo0−ωxoωzo−ωxoωyoωzo2−ωxo2ωyoωzo−ωxoωyoωxoωyo0ωxoωzo−ωyoωzoωxo2−ωyo2.The system dynamics equation in the orbit system can be reduced to the following equation:(51)D1ω˙oi+D2ωoi−3ω02D2R^i=Tco+Tdo.Since the angular acceleration information is not measurable and the angular velocity information is measurable, the left and right sides of equation (51) are integrated to obtain the following results:(52)D1ωo|t1t0+∫t0t1D2ωo−3ω02D2R^dti=∫t0t1Tco+Tdodt.If equation (52) is written in the form of standard linear least squares, then(53)Φ=D1ωo|t1t0+∫t0t1D2ωo−3ω02D2R^dt,y=∫t0t1Tco+Tdodt.In order to realize online real-time identification, the recursive form of least squares is used. Recursive least square estimation extracts the information of the estimated quantity from each measurement, which is used to modify the estimation obtained in the previous step [16]. The more times of measurement, the more times of correction, and the higher the accuracy of estimation.When the system is brought into TEA, the external force torques cancel each other, and the control torque required is small, which cannot meet the continuous excitation conditions required by the traditional least square method. Moreover, due to the “data saturation” phenomenon, when the observation data increase, the deviation between the estimated value obtained by the recursive least square parameter identification method and the real parameters will become larger and larger. Therefore, in parameter estimation, we should pay enough attention to the current data and gradually forget the old data which does not contain the current dynamic characteristics, so we use the recursive least squares estimator with bounded gain forgetting factor. The estimator still has bounded gain when the excitation is not continuous [17].SetX to be measured, Yj is the j measurement, and the measurement equation is(54)Yj=ΦjX+Vj,j=1,2,…,k,where Φj and Vj are the j measurement matrix and random measurement noise. The recursive least square estimation algorithm with the forgetting factor is as follows:(55)Kk=λE+Φk+1PkΦk+1T−1,Pk+1=1λPk−PkΦk+1TKkΦk+1Pk,X^k+1=X^k+Pk+1Φk+1TZk+1−Φk+1X^k.The choice of forgetting factor is of great significance to the stability of the system. When the continuous excitation is satisfied (such as in the dynamic process of tracking TEA), the forgetting factor of zero will lead to zero gain (i.e., it degenerates to the standard least squares method, resulting in the instability of tracking time-varying parameters) [18]. When the continuous excitation condition is not satisfied, the forgetting factor of the normal number will lead to a sharp increase in gain. The forgetting factor is chosen as the following bounded form:(56)λt=λ01−Pk0.λ0 and k0 is a normal number. The upper bounds of the maximum forgetting rate and the norm of the gain matrix are given, respectively; P is the gain matrix and represents the level of motivation. Formula (56) means when the norm of P is small (strong continuous incentive), the forgetting factor is λ0, forgetting is fast, and the system has a strong ability to track the changing parameters. When the norm of P increases, the forgetting speed decreases. And the forgetting speed is zero when the norm of P reaches a specified upper bound. In this case, λ0 is 0.95, and k0 is 109.


## 5. Controllability and Singularity

The nonlinear state transformation (34)–(36) are not global transformation. They are effective at the nonsingular points of the matrix in the new state equation. Nonsingular points are effective. In this case, the inverse matrix is(57)E−1t=−13ω03IoS1o−1.Formula (57) shows that the reversibility of E(x) depends on the reversibility of S1o. According to the definition of S1o, its reversibility is determined by the following formula:(58)Δ=detS1o=Ixo−IyoIyo−IzoIzo−Ixo+Ixo−IyoIxzo2+Izo−IxoIxzo2+Iyo−IzoIyzo2.Due to the existence of the torque of inertia of the inertial coordinate system of the system mass center and the track system, the following relations are(59)Io=LopIpLpo.Δ can be changed to(60)Δ=Ixp−IypIyp−IzpIzp−IxpΛφ∗,θ∗,Ψ∗.This is a function of the rotation Euler angle (narrow sense TEA) from the center of mass inertial principal axis coordinate system to the orbit system when the 3-1-2 rotation sequence is adopted, and the simplified expression ofΛφ∗,θ∗,Ψ∗ is as follows:(61)Λφ∗,θ∗,Ψ∗=cos2φ∗cos2θ∗cos2ψ∗+14sin2φ∗sin2θ∗sin2ψ∗1−3cos2φ∗.When the following two conditions are met,Δ≠0(62)Ixp≠Iyp≠Iyp≠Izp,Izp≠Ixp,Λ≠0.Among them, the first condition is a physical condition, which restricts the inertia characteristics of the system and can feedback linearization. The second condition limits the attitude range of the control law (43).


### 5.1. Torque of Inertia Constraint

Inertia constraintsIxp≠Iyp≠Iyp≠Izp,Izp≠Ixp are the inherent characteristic of the ACMM problem which only considers the gravity gradient torque, and even for the linear controller, the same constraint is needed. For example, assuming that the nonlinear dynamic equations of the system are linearized at zero attitude, the dynamic equations of pitch and roll/yaw in the ACMM system are decoupled. For pitch channel, take the state variable as x2=θ∗θ˙∗hcpT, the equation of state is(63)x˙2=A2x+b2u2p,where(64)A2=1Iyp0Iyp03ω02Izp−Ixp00000,b2=1Iyp0−1Iyp.The controllable matrix is:(65)S=1Iyp20−Iyp0−Iyp03ω02Izp−IxpIyp200.WhenIzp≠Ixp, rankS=2<3. The pitch axis is uncontrollable, and the roll/yaw axis has a similar conclusion.


### 5.2. Attitude Constraint

The nonlinear control law (39) requires Λ≠0. In order to avoid singularity, it is necessary to study the distribution of attitude angle when Λ=0. When using the 3-1-2 rotation sequence, according to the expression of Λ, and Ixp≠Iyp≠Iyp≠Izp,Izp≠Ixp, the singular surface near Θ∗=000 is obtained, as shown in Figure 3.Figure 3
Singular surface of the control law.The two surfaces in Figures4 and 5 will separate the TEA at θ∗=0,0,0 from other TEA. If the initial attitude of the system is located in the area surrounded by the surface in the figure, it cannot pass through the surface in the process of tracking TEA, so as to avoid the singularity of the nonlinear control law.Figure 4
Value of Λ whenφ∗=0°.Figure 5
Value of Λ whenφ∗=0° (close up view).Because Λ a function of three attitude angles, it is hard to observe Λ with the change of attitude angle. Considering that the change of yaw angle is more obvious in the process of configuration transformation, and the gravitational gradient torque is generally more obvious in the pitch axis, the roll angle is fixed here to investigate the influence of pitch angle and yaw angle on Λ.Figures6 and 7 show the value of Λ when φ∗=0, θ∗, ψ∗∈−180°,180°. As can be seen from the figure, φ∗=0, corresponding to the balanced attitude(66)Θ∗=mπ2nπ20,where the value of Λ is (−1)m + n, and it is local maximum 1 or local minimum −1. When φ∗=0, θ∗=±45°, ψ∗=±45°, Λ = 0, so if the initial time θ∗<45°, ψ∗<45°, φ = 0, the uniqueness TEA θ = (0, 0, 0) that can be achieved by a combination not passing through the singular surface of Λ = 0 in a graph. The value of Λ is 1.Figure 6
The value of Λ whenφ∗=−90°.Figure 7
The value of Λ whenφ∗=−45°.Figures6 to 13 show that when φ∗ changes from −90° to 90°, Λ changes with the θ∗, ψ∗. This set of curves shows that the value of Λ is symmetrically distributed with respect to θ∗=ψ∗ and θ∗=−ψ∗. When φ∗ increased from −90° to −45°, along θ∗=ψ∗, it is constant. But along the θ∗=−ψ∗, the amplitude of the fluctuation of Λ becomes smaller, as shown in the figure. When φ∗ increased to −45°, a single peak and trough appeared; when φ∗ increased from −45° to 0°, the amplitude of wave crest and wave trough becomes larger; when φ∗=0, the peak and trough values reach the maximum value of 1 and the minimum value of −1, respectively; when φ∗ increased from 0 to 45°, the amplitudes of wave crest and wave trough decrease, respectively; when φ∗ increased from 45° to 90°, Λ is symmetric about θ∗=ψ∗ and θ∗=−ψ∗ again. However, Λ rotates at 90 degrees. When θ∗=−ψ∗, Λ is a constant.Figure 8
The value of Λ whenφ∗=−60°.Figure 9
The value of Λ whenφ∗=−30°.Figure 10
Value of Λ whenφ∗=30°.Figure 11
Value of Λ whenφ∗=60°.Figure 12
Value of Λ whenφ∗=45°.Figure 13
Value of Λ whenφ∗=90°.It can be seen from the previous analysis that when only considering the gravitational gradient torque, the narrow sense TEA has the following form:(67)Θ∗=mπ2nπ2pπ2.Before applying the nonlinear control law, the distance between TEA and the singular point must be determined. The nonlinear control law is effective only when Λ at TEA is far from zero.Take the first and second partial derivatives of Λ:(68)∂Λ∂Θ∗=∂Λ∂φ∗∂Λ∂θ∗∂Λ∂ψ∗T,where(69)∂Λ∂φ∗=−2sin2φ∗cos2θ∗cos2ψ∗+14sin2θ∗sin2ψ∗cosφ∗1−3cos2φ∗+6sinφ∗sin2φ∗,∂Λ∂θ∗=−2cos2φ∗sin2θ∗cos2ψ∗+12cos2θ∗sin2ψ∗sinφ∗1−3cos2φ∗,∂Λ∂ψ∗=−2cos2φ∗cos2θ∗sin2ψ∗+12sin2θ∗cos2ψ∗sinφ∗1−3cos2φ∗,∂Λ∂Θ∗2=∂Λ∂φ∗2∂Λ∂φ∗∂θ∗∂Λ∂φ∗∂ψ∗∂Λ∂φ∗∂θ∗∂Λ∂θ∗2∂Λ∂θ∗∂ψ∗∂Λ∂φ∗∂ψ∗∂Λ∂θ∗∂ψ∗∂Λ∂ψ∗2T,where(70)∂Λ∂φ∗2=−4Λ+154sinφ∗sin2θ∗sin2ψ∗1+3cos2φ∗,∂Λ∂θ∗2=−4Λ,∂Λ∂ψ∗2=−4Λ,∂Λ∂φ∗∂θ∗=4sin2φ∗sin2θ∗cos2ψ∗+12cosφ∗cos2θ∗sin2ψ∗7−9cos2φ∗,∂Λ∂φ∗∂ψ∗=4sin2φ∗sin2θ∗cos2ψ∗+12cosφ∗cos2θ∗sin2ψ∗7−9cos2φ∗,∂Λ∂θ∗∂ψ∗=4cos2φ∗sin2θ∗sin2ψ∗+sinφ∗cos2θ∗cos2ψ∗1−3cos2φ∗.For TEA in the form of the formula (65), there are(71)∂Λ∂Θ∗|TEA=000T,(72)∂Λ∂Θ∗2=4−1n+p100010001mis even4−1n+p−1m000−1m−1m−1/20−1m−1/2−1m,mis odd.It can be seen that ifm is even, the eigenvalues of the second-order Hessian matrix are all 4(−1)n + p. When n + p is even, TEA is the minimum of Λ. When n + p is odd, TEA is the largest value of Λ. If m is odd, then the eigenvalue of the Hessian matrix is 8(−1)n + p, 4(−1)n + p, and 0. When n + p is an even number, θ∗=ψ∗±n−pπ/2φ∗<0 and θ∗=−ψ∗±n−pπ/2φ∗>0 (k is an integer) Λ are both minima −1; when n + p is an odd number, θ∗=ψ∗±n−pπ/2φ∗<0 and θ∗=−ψ∗±n−pπ/2φ∗>0 Λ are both maximum 1.It can be seen from the above analysis that a TEA is a maximum or a minimum. So, the nonlinear feedback control law given by formula (44) can avoid singular points in principle. But, when the initial attitude is not suitable, the system may still pass through the singular surface during the stabilization process. For example, if the initial angular rate of attitude is relatively large, the system is likely to reach the singular surface before stabilization. Therefore, whether the nonlinear control law can effectively drive the system to TEA is closely related to the initial attitude. If the reference input signal zd is selected reasonably, the singularity may be avoided in some programs, and the system will approach TEA along the ideal trajectory.


## 5.1. Torque of Inertia Constraint

Inertia constraintsIxp≠Iyp≠Iyp≠Izp,Izp≠Ixp are the inherent characteristic of the ACMM problem which only considers the gravity gradient torque, and even for the linear controller, the same constraint is needed. For example, assuming that the nonlinear dynamic equations of the system are linearized at zero attitude, the dynamic equations of pitch and roll/yaw in the ACMM system are decoupled. For pitch channel, take the state variable as x2=θ∗θ˙∗hcpT, the equation of state is(63)x˙2=A2x+b2u2p,where(64)A2=1Iyp0Iyp03ω02Izp−Ixp00000,b2=1Iyp0−1Iyp.The controllable matrix is:(65)S=1Iyp20−Iyp0−Iyp03ω02Izp−IxpIyp200.WhenIzp≠Ixp, rankS=2<3. The pitch axis is uncontrollable, and the roll/yaw axis has a similar conclusion.


## 5.2. Attitude Constraint

The nonlinear control law (39) requires Λ≠0. In order to avoid singularity, it is necessary to study the distribution of attitude angle when Λ=0. When using the 3-1-2 rotation sequence, according to the expression of Λ, and Ixp≠Iyp≠Iyp≠Izp,Izp≠Ixp, the singular surface near Θ∗=000 is obtained, as shown in Figure 3.Figure 3
Singular surface of the control law.The two surfaces in Figures4 and 5 will separate the TEA at θ∗=0,0,0 from other TEA. If the initial attitude of the system is located in the area surrounded by the surface in the figure, it cannot pass through the surface in the process of tracking TEA, so as to avoid the singularity of the nonlinear control law.Figure 4
Value of Λ whenφ∗=0°.Figure 5
Value of Λ whenφ∗=0° (close up view).Because Λ a function of three attitude angles, it is hard to observe Λ with the change of attitude angle. Considering that the change of yaw angle is more obvious in the process of configuration transformation, and the gravitational gradient torque is generally more obvious in the pitch axis, the roll angle is fixed here to investigate the influence of pitch angle and yaw angle on Λ.Figures6 and 7 show the value of Λ when φ∗=0, θ∗, ψ∗∈−180°,180°. As can be seen from the figure, φ∗=0, corresponding to the balanced attitude(66)Θ∗=mπ2nπ20,where the value of Λ is (−1)m + n, and it is local maximum 1 or local minimum −1. When φ∗=0, θ∗=±45°, ψ∗=±45°, Λ = 0, so if the initial time θ∗<45°, ψ∗<45°, φ = 0, the uniqueness TEA θ = (0, 0, 0) that can be achieved by a combination not passing through the singular surface of Λ = 0 in a graph. The value of Λ is 1.Figure 6
The value of Λ whenφ∗=−90°.Figure 7
The value of Λ whenφ∗=−45°.Figures6 to 13 show that when φ∗ changes from −90° to 90°, Λ changes with the θ∗, ψ∗. This set of curves shows that the value of Λ is symmetrically distributed with respect to θ∗=ψ∗ and θ∗=−ψ∗. When φ∗ increased from −90° to −45°, along θ∗=ψ∗, it is constant. But along the θ∗=−ψ∗, the amplitude of the fluctuation of Λ becomes smaller, as shown in the figure. When φ∗ increased to −45°, a single peak and trough appeared; when φ∗ increased from −45° to 0°, the amplitude of wave crest and wave trough becomes larger; when φ∗=0, the peak and trough values reach the maximum value of 1 and the minimum value of −1, respectively; when φ∗ increased from 0 to 45°, the amplitudes of wave crest and wave trough decrease, respectively; when φ∗ increased from 45° to 90°, Λ is symmetric about θ∗=ψ∗ and θ∗=−ψ∗ again. However, Λ rotates at 90 degrees. When θ∗=−ψ∗, Λ is a constant.Figure 8
The value of Λ whenφ∗=−60°.Figure 9
The value of Λ whenφ∗=−30°.Figure 10
Value of Λ whenφ∗=30°.Figure 11
Value of Λ whenφ∗=60°.Figure 12
Value of Λ whenφ∗=45°.Figure 13
Value of Λ whenφ∗=90°.It can be seen from the previous analysis that when only considering the gravitational gradient torque, the narrow sense TEA has the following form:(67)Θ∗=mπ2nπ2pπ2.Before applying the nonlinear control law, the distance between TEA and the singular point must be determined. The nonlinear control law is effective only when Λ at TEA is far from zero.Take the first and second partial derivatives of Λ:(68)∂Λ∂Θ∗=∂Λ∂φ∗∂Λ∂θ∗∂Λ∂ψ∗T,where(69)∂Λ∂φ∗=−2sin2φ∗cos2θ∗cos2ψ∗+14sin2θ∗sin2ψ∗cosφ∗1−3cos2φ∗+6sinφ∗sin2φ∗,∂Λ∂θ∗=−2cos2φ∗sin2θ∗cos2ψ∗+12cos2θ∗sin2ψ∗sinφ∗1−3cos2φ∗,∂Λ∂ψ∗=−2cos2φ∗cos2θ∗sin2ψ∗+12sin2θ∗cos2ψ∗sinφ∗1−3cos2φ∗,∂Λ∂Θ∗2=∂Λ∂φ∗2∂Λ∂φ∗∂θ∗∂Λ∂φ∗∂ψ∗∂Λ∂φ∗∂θ∗∂Λ∂θ∗2∂Λ∂θ∗∂ψ∗∂Λ∂φ∗∂ψ∗∂Λ∂θ∗∂ψ∗∂Λ∂ψ∗2T,where(70)∂Λ∂φ∗2=−4Λ+154sinφ∗sin2θ∗sin2ψ∗1+3cos2φ∗,∂Λ∂θ∗2=−4Λ,∂Λ∂ψ∗2=−4Λ,∂Λ∂φ∗∂θ∗=4sin2φ∗sin2θ∗cos2ψ∗+12cosφ∗cos2θ∗sin2ψ∗7−9cos2φ∗,∂Λ∂φ∗∂ψ∗=4sin2φ∗sin2θ∗cos2ψ∗+12cosφ∗cos2θ∗sin2ψ∗7−9cos2φ∗,∂Λ∂θ∗∂ψ∗=4cos2φ∗sin2θ∗sin2ψ∗+sinφ∗cos2θ∗cos2ψ∗1−3cos2φ∗.For TEA in the form of the formula (65), there are(71)∂Λ∂Θ∗|TEA=000T,(72)∂Λ∂Θ∗2=4−1n+p100010001mis even4−1n+p−1m000−1m−1m−1/20−1m−1/2−1m,mis odd.It can be seen that ifm is even, the eigenvalues of the second-order Hessian matrix are all 4(−1)n + p. When n + p is even, TEA is the minimum of Λ. When n + p is odd, TEA is the largest value of Λ. If m is odd, then the eigenvalue of the Hessian matrix is 8(−1)n + p, 4(−1)n + p, and 0. When n + p is an even number, θ∗=ψ∗±n−pπ/2φ∗<0 and θ∗=−ψ∗±n−pπ/2φ∗>0 (k is an integer) Λ are both minima −1; when n + p is an odd number, θ∗=ψ∗±n−pπ/2φ∗<0 and θ∗=−ψ∗±n−pπ/2φ∗>0 Λ are both maximum 1.It can be seen from the above analysis that a TEA is a maximum or a minimum. So, the nonlinear feedback control law given by formula (44) can avoid singular points in principle. But, when the initial attitude is not suitable, the system may still pass through the singular surface during the stabilization process. For example, if the initial angular rate of attitude is relatively large, the system is likely to reach the singular surface before stabilization. Therefore, whether the nonlinear control law can effectively drive the system to TEA is closely related to the initial attitude. If the reference input signal zd is selected reasonably, the singularity may be avoided in some programs, and the system will approach TEA along the ideal trajectory.


## 6. Numerical Examples of Configuration Transformation Process

It is assumed that the mass characteristics of the two sections are shown in the table (Table1).Table 1
Configuration transformation cabin and configuration parameters.
Replenishment vehicleReplenished vehicleOne lineL configurationMass (kg)22000220004400044000Torque of inertia (105 kgm2)Ix1.301.302.6015.23Iy6.176.1739.5614.27Iz6.866.8640.9426.90Ixy0006.934Ixz0000Iyz0000Suppose that the configuration transformation process takes 5000 seconds, and in the first 1000 seconds, the replenished vehicle is pushed out 5 meters along the longitudinal axis of the replenishment vehicle with a mechanical arm at a constant speed; when 1000 seconds to 4000 seconds, rotate the replenished vehicle 90° around the negative direction of the yaw axis of the replenishment vehicle. Considering the transfer time and the speed limit of the end of the manipulator, adopt the rotation process of acceleration uniform deceleration. When 1000 seconds to 2000 seconds, rotate the replenished vehicle with constant acceleration, when 2000 seconds to 3000 seconds, rotate with uniform speed, and when 3000 seconds to 4000 seconds, rotate with uniform deceleration. The ideal situation is that when the replenished vehicle rotates 90° relative to the replenishment vehicle, the angular velocity of the relative replenishment vehicle can just be reduced to zero. From 4000 seconds to 5000 seconds, the replenished vehicle is pulled to the lateral interface of the node cabin at a constant speed along they-axis direction of the replenishment vehicle by the mechanical arm to complete the redocking. From the process of configuration transformation, the relative attitude of the replenishment vehicle and the replenished vehicle is mainly reflected in the yaw axis.In the simulation model, the replenishment vehicle is a rigid body dynamic model with flexible appendages. The solar panel of the replenished vehicle is retracted during the transfer process, regardless of the influence of its flexible mode, the flexibility of the manipulator, and its dynamic relationship. It is assumed that it can move according to the planned ideal trajectory, and the influence of atmospheric disturbance is considered during the transfer process. In the process of configuration transformation, the replenished vehicle to be transferred is in an uncontrolled state, and the whole system is controlled by a set of CMGs of pentagonal pyramid configuration. The nominal angular torque capacity of each CMG is 1000 nms, and the minimum angular torque on the angular torque envelope of each CMG is 4200 nms. It is still assumed that the combined body runs in a circular orbit 400 km away from the earth’s surface, and the orbital angular velocity is lowerω0 is a constant value of 0.0011 rad/s, and the sampling period is 200 ms. The initial attitude angle and angular velocity are assumed to be 1 θ = [0.0007 0.0008 0.001]°, θ = [0.001–0.001 0.001]°/s. The angular torque of the initial CMG is calculated from the initial frame angle of each CMG in the pentagonal pyramid configuration.Without considering the influence of the manipulator, the relative attitude of the replenishment vehicle and the replenished vehicle in the transfer process is shown in Figure14. If the replenishment vehicle is in the uncontrolled free-floating state during the configuration transformation, the attitude motion of the replenishment vehicle during the transfer process is shown in Figure 15. If the angular torque management of CMG is not carried out in the process of configuration transformation, and the replenishment vehicle maintains directional stability to the ground, the angular torque that CMG needs to absorb in the process of configuration transformation is shown in Figure 16. As can be seen from Figure 17, if the replenishment vehicle is not controlled during the transfer process, the attitude angle drift is too large. As shown in Figures 16 and 18, if the replenishment vehicle is strictly oriented to the ground during the transfer process, the disturbance angular torque absorbed by CMG will be large and reach saturation quickly, which will affect the accuracy of the transfer process.Figure 14
Relative attitude angles of replenishment vehicle and replenished vehicle.Figure 15
Angular torque of CMG under directional stability.Figure 16
Attitude angle of replenishment vehicle in free floating.Figure 17
Main inertia identification results.Figure 18
CMG singularity measure for earth-oriented stability.After torque management is adopted, the identification results of the parameter identification unit are shown in Figures17 and 19 when supplying the narrow TEA of the aircraft attitude tracking system during the transfer process. The closed-loop identification effect of the least square method with bounded gain forgetting factor is ideal, and it can track the changing torque of inertia of the system in the whole configuration transformation process.Figure 19
Identification result of inertia product.The angular torque and attitude angle stored in CMG are shown in Figures20–25. After the angular torque planning, the angular torque of CMG is strictly controlled within its capacity range during the whole transfer process (Figures 20 and 22), and the singular measure of CMG is far away from zero (Figure 23). The attitude maneuver of the replenishment vehicle is mainly reflected on the yaw axis (Figures 21 and 24), which can be verified by the relative motion of the replenishment vehicle and the replenished vehicle during the transfer process. In the controller design, the change history of nine states of the equivalent linear system is shown in Figure 25. Except the fifth state needs to track the time-varying reference input and the eighth rolling/pitching inertia product changes constantly during the transfer process, the other amplitudes fluctuate less, and the dynamic effect is ideal.Figure 20
CMG angular torque (narrow TEA tracking).Figure 21
Attitude angle of replenishment vehicle.Figure 22
Control torque (narrow TEA tracking).Figure 23
CMG singular measure (narrow TEA tracking).Figure 24
Attitude angular velocity of replenishment vehicle.Figure 25
Equivalent linear system state variables.


## 7. Conclusion

This paper studies the control problem of the variable mass body in the process of on-orbit service. Different from the traditional method, the feedback linearization method of the nonlinear system is used to design the controller. By realizing the tracking of narrow TEA in the process of configuration transformation, the error caused by the angular deviation between narrow TEA and earth orientation attitude caused by the large change of moment of inertia is avoided, and the ACMM controller is obtained. In order to improve the input information accuracy of the ACMM controller, the least square method with a bounded gain forgetting factor is introduced to track the real variable parameters. This method does not need continuous excitation. A joint attitude control/angular momentum management controller is designed to effectively stabilize the system and reduce the angular momentum saturation of the attitude control system during on-orbit service. In the numerical simulation of the configuration transformation process, we compared the changes in CMG angular momentum management in the configuration transformation process. If the supply aircraft is kept strictly oriented to the ground during the transfer process, the disturbance angular momentum absorbed by CMG will quickly reach saturation and need to be unloaded, which will affect the accuracy of the transfer process. If CMG management is carried out and the least square method with bounded gain forgetting factor is used for closed-loop parameter identification, CMG saturation is effectively controlled. This shows that the ACMM controller can track the changing moment of inertia of the system, which provides a solution to the problem of variable quality control in on-orbit service.



---

*Source: 9660218-2022-06-13.xml*