A new method to identify the mass parameters of spacecraft

Abstract A two-step identification method is proposed, both the moment of inertia and the mass properties are identified. A new index parameter which is different from the commonly used condition number is first defined for designing the optimal excitation. A method is introduced based on the least squares algorithm. Detailed mathematical equations and numerical analysis are proposed. Simulation results show that the proposed method can identify the mass parameters of spacecraft.


Introduction
On orbit identification of mass properties is a classical research for spacecraft attitude control. Based on the conservation of angular momentum, Kun et al. (2017) presents a method to design and calculate the optimal excitation. A simple direct-search method is applied to calculate the optimal excitation. But the algorithm can only identify the moment of inertia. Huang et al. (2016) derive a dynamics model. Then, a full adaptive controller is introduced including a control law, a dynamic adaption law, and a kinematic adaption law. Zhang et al. (2015) presents a new method that the aircraft are modeled as rigid bodies, there are three steps to solve the problem. Xie et al. (2016) proposed a method which uses Lagrange method to establish dynamic equation. The PSO(particle swarm optimization) algorithm, Observer-Kalman filter identification and ERA(eigensystem realization algorithm) are combined to identify the mass property.
In 1987, Bergmann et al. (1987) first used second order filter and Kalman filter to identify the mass property of spacecraft. Bergmann and Dzielski (1990) also found that using only torque-producing actuator of the flight, such as the momentum wheel, the mass property of spacecraft can also be identified. Richfield et al. (1988) found symbolic manipulation is too complex for online use, they proposed a suboptimal direct strategy which can reduce the time of convergence. Yang et al. (2014) found a linear regression model through attitude dynamics equations, the traditionally used condition number was defined as the critical index. Tanygin and Williams (1997) uses coasting maneuvers as the excitation in mass parameters estimation process. This method improves the effectiveness of estimation. Palimaka and Burlton (1992) used a weighted-least-squares method and on-orbit gyro data to identify mass parameters.
On the basis of angular-momentum conservation, Norman et al. (2011) used the batch-least-squares-method and the Extended-Kalman-Filter to estimate the mass property in spite of sensor noise. Ma et al. (2008) presents a twostep method, 1st step is to identify the mass properties, 2nd step is to identify the inertia matrix. But the method is robotics-based, it requires measuring velocities which is a little difficult than forces and accelerations. Sekhavat et al. (2009) pointed out that the optimal excitation trajectory is found by optimizing the traditionally used condition number. An unscented dual state-parameter Kalman-Filter is adopted to estimate the mass property. Lam et al. (1991) uses extended Kalman filter and neuro-filter to do the simulation work. He found that the neuro-filter is more robust when the data obtained was not well. Different from the traditionally used condition number, a new index parameter which is different from the commonly used condition number is first defined for designing the optimal excitation.
In this paper we propose another method, both the inertia matrix and the position of the center of mass are iden-tified. The on-orbit mass parameter identification method is introduced based on the least squares algorithm. A detailed mathematical equation and numerical analysis are introduced. The remainder of this paper is organized as follows. Chapter 2 introduces the identification method of mass property of aircraft. Chapter 3 shows the adopted optimal excitation. The identification method under special circumstances is presented in Chapter 4. Numerical simulations that demonstrate performances of the proposed method are presented in Chapter 5. This paper is brought to a close with some concluding remarks in Chapter 6.

The Classical Least Squares Algorithm(LS)
The system measurement equation at time t k is Where X is the n×1 parameter vector to be estimated, φ k is the kth time 3×n measurement matrix, y k is the kth time output 3×1 vector n k is disturbance 3×1 vector. We can get φ k and y k from the telemetry data of sensors on orbit. When we have k sets of telemetry data, the system measurement equation is When the disturbance N k = 0, the unique solution is Usually the disturbance N k ≠ 0, so we get In order to minimize the influence of disturbance, we have to choose 3k>n.Ŷ WhereŶ k is the optimal estimation of Y k ,Xis the optimal estimation of X. Define the error vector e k Using the classical Least Squares Algorithm to find the minimum value of J.
Let ∂J ∂X = 0, we can get From the telemetry data, we know that Φ k is full rank ma- Equation (9) is the classical Least Squares Algorithm. But this algorithm is suitable for off-line estimation, we have to calculate pseudo-inverse matrix which is very large.

Recursive Least Squares Algorithm (RLS)
In Equation (35), the obtained moment of inertia is the best estimation using batch least squares algorithm. Batch least squares algorithm adopts the least squares algorithm to obtain the best estimation value once all the system input and output data are obtained. But this algorithm has the following disadvantages: The more data sampled, the higher the accuracy of the estimation. In order to obtain accurate identification results, the coefficient matrix is usually quite large, so the computational complexity is high, which cannot be used for online identification. The coefficient matrix and inverse matrix need to be recalculated every sampling interval. The algorithm can't be used if the coefficient matrix is not a full column rank matrix according to the matrix theory.
For these reason, the recursive least squares algorithm is introduced in 1950 by Placket. The recursive least squares algorithm is as follows: every time when it obtains new observation data, the algorithm uses the newly obtained observation data to correct the previous estimation result, then repeating estimation process recursively. The calculation is carried out one after another until a satisfactory accuracy. The algorithm renew the old parameter estimation result with the newly acquired data. The algorithm not only reduces the computational complexity and storage size, but also realizes online real-time identification.

The recursive least squares algorithm
Using the batch least square algorithm, the estimated value at time k isX Similarly, the estimation value at time k + 1 iŝ In order to get the recursive squares algorithm, we need to: 1. Recursive calculation of the inverse matrix Constructing information inverse matrix P k+1 (also known as covariance matrix) Then we get The matrix inversion formula is introduced Where A and C is the invertible matrix.
Recursive calculation of parameter estimatê 14) and (15) is the recursive least squares algorithm. Define K k as gain matrix.
Then (17) and (18) is another form of the recursive least squares algorithm. The selection of initial value of the recursive least squares algorithm is important. Traditionally, there are two ways to select the initial value: b) chooseX 0 as 0 vector; P 0 = αE, here the real number αis sufficiently large (10 5~1 0 6 ). b) The preliminary estimation ofX L and P L is based on initial sampling data. Then the value acquired is the initial value of the recursive least squares algorithm.

Weighted recursive-least-squares-algorithm
Suppose at time k, the weighting factor is λ k , the weighting matrix is Then the weighted recursive-least-squares-algorithm is

Forgetting Factor Recursive Least-Squares
Algorithm During the process of recursive calculation, when a new set of data is obtained, all the previous data are processed by a weighting factor ρ, 0 < ρ < 1.

Identification of mass property of aircraft
The dynamical equation of rigid body aircraft is The rotational dynamical equation of rigid body aircraft is Where Tc is control torque, T d is disturbance torque. If thruster is used as the actuator, then Where N is the number of thrusters, the value of λ j denoting the state of the thruster, λ j = 1 denotes the thruster is on, λ j = 0 denotes the thruster is off. Fc is the output force of thruster, D j is the direction of the output force of the jth thruster in spacecraft body coordinate system, r j is the position vector of the jth thruster in vehicle body coordinate system, rc is position vector of the center of mass of the spacecraft in vehicle body coordinate system. We need to identify the parameter M, rc, J.

Parameter identification of moment of inertia based on angular momentum conservation (using momentum wheels)
According to the angular momentum conservation Law, the dynamic equation of spacecraft is Where J is the moment of inertia matrix of the rigid spacecraft (without momentum wheels), ω is the angular velocity of spacecraft. Hw is the angular momentum of wheels, H 0 is the initial angular momentum of the combination of rigid spacecraft and momentum wheels.
∫︀ t f t0 T d is disturbance torque, t 0 is the initial time, t f is the finish time of the identification process. Let H d = ∫︀ t f t0 T d . Usually the disturbance torque is very small compared with the applied control torque of momentum wheel. In order to simplify the situation, let H d = 0.
Considering the special situation, when H 0 = 0, From equation (31), Then equation (31) is Equation (32) is a linear system of J v , at given time t i , we can get the angular velocity ω from sensor data, we can get Hw from the telemetry data of momentum wheels.
Suppose we get n sets of data, then If n≥ 2, Av is a full column rank matrix, The equation (34) is an over-determined linear system, we can used the least squares algorithm to get the value of J v , then we get the best estimate of J v .

Identification of moment of inertia using momentum wheels
1. When T d = 0 Whereẇ is the derivative of w, Tc is the control torque.ω is the cross product matrix. Usually the disturb torque T d is relatively small compared with the applied control torque of momentum wheel, suppose T d = 0. At given time t i , we can get the angular velocity w from gyroscope data, we can getẇ from angular accelerometer data. We can get Tc from telemetry data of momentum wheel. If there is n sets of data, If Av is a full column rank matrix, The equation (39) is an over-determined linear system, we can used the least squares algorithm to get the value of J v , then we get the best estimation of J v .
2. When T d ≠ 0 During two successive sampling time t i and t i+1 , from equation (37), Usually the time between two successive sampling time t i and t i+1 is very short, suppose T di+1 − T di = 0, then If we get (n(n is even number) sets of data, then If Av is a full column rank matrix, The equation (46) is an over-determined linear system, we can used the least squares algorithm to get the value of J v , then we get the best estimate of J v .

Identification of moment of inertia using thrusters According to equation (28)
Where rc is position vector of the center of mass. rm is position vector of the accelerometer. f j is the vector of output force of thrusters.
(︀ r j − rc )︀ is the cross product matrix. J and rc are the parameters needed to be identified.
If we only want to identify the moment of inertia J, the thruster can be installed in symmetrical position. Then Subsequent deduction is similar to section 2.1.2.1.

Parameter identification of mass 2.2.1 Identification of mass and center of mass at the same time using thrusters
Where is the velocity data of the point where the accelerometer is stalled, υc is the velocity of the center of mass, rc is position vector of the center of mass. rm is position vector of accelerometer.ω is the cross product matrix. f j is the vector of output force of thrusters.
Take the derivative of equation (50), Whereυm is the output of the accelerometer, ωc is the output of gyroscope data, rm is position vector of accelerometer. f j is the force vector of thrusters, M is the mass of spacecraft. Define Then equation (51) is BMv =υm (52) When there is n sets of data, If Bv is a full column rank matrix, The equation (54) is an over-determined linear system, we can used the least squares algorithm to get the value ofM v then get the best

Combined estimation of moment of inertia, mass and center of mass using thrusters
Change the form of equation (48), Change the form of equation (51) [ Then we get Equation (57) can be used to estimate the moment of inertia, mass and center of mass at the same time.
3 Optimal excitation for identification of mass 3.1 Identification method using angular momentum conservation Kun et al. (2017) found that excitation design is one of the important aspect in the identification of inertia parameters and the form of excitation has a great influence on the identification result. A performance index which is similar to but not the condition number is first defined as the benchmark for designing the optimal excitation. But the article only can identify the moment of inertia, there is still mass and center of mass need to be identified. This paper proposes a new method based on the article above-mentioned.
Because of the conservation of angular momentum, At time t k , Where Suppose there is N sets of data, Then Ψ d,N is usually impossible to measured, equation (65) is an over-determined linear system, we can used the least squares algorithm to get the best estimate.
From equation (65) Ψ The identification error is In equation (67) , in order to get high accuracy estimation, on one hand, we can decrease the value of ⃦ ⃦ Ψ d,N ⃦ ⃦ , e.g. adopting the input shaper to eliminate the measurement noise, on the other hand, we can decrease the value of κ, e.g. choose the optimal excitation.
When considering the measurement noise of ω, equa- Where ∆Ω N denoting the measurement noise of ω.
The novel idea of this paper is that we found κ is different from the cond (︀ Ψ w,N )︀ . Then the identification error is In order to find κ, first we write down the It is a 6-order equation about λ, we can't get an analytical solution, we can only solve it numerically. When we get the λ. The optimal excitation is obtained by solving the following optimal control problem.
Subject toḢw = Tw The equation (81) can only ensure finding the optimal excitation at final time duration. In order to maintain κ at minimum value during the whole identification process, the optimal excitation is obtained by solving the following optimal control problem.
Subject toḢw = Tw Subject toḢw = Tw We find that in equation (81), equation (82) and equation (83), the optimal excitation is not the same.
Because the subsequent value of κ is not bigger than the previous value of κ. Then Then We also can use the dual merit-index to find the optimal excitation. Firstly, minimizing the scalar variable κ, secondly, minimizing the change of angular velocity during the identification process. Which means the initial angular velocity is equal to the final angular velocity. The optimal excitation is obtained by solving the following optimal control problem.
Subject toḢw = Tw If we choose an arbitrary time t f , in order to minimize κ (κ is getting smaller as the increase of sampling data), the difference between the initial angular velocity and the terminal angular velocity is small. The optimal excitation is obtained by solving the following optimal control problem.
Subject to f (Hw) = min Where κe and ωe denoting the expected condition number and expected initial angular velocity deviation.

Identification method using rotational dynamical equation
Above-mentioned identification algorithm uses the angular momentum conservation. There is another algorithm from the perspective of rotational dynamics.

Identification of mass and center of mass using momentum wheels
According to angular momentum conservation ]︁ T (100) Then When there are N sets of data Change form into N is measurement noise and is impossible to estimate, equation (105) is an over-determined linear system, we can used the least squares algorithm to get the best estimate.
The identification error is (107) is The process of identification of moment of inertia and identification of mass and center of mass are independent. We can identify them separately.
In this paper, a two-step identification method is proposed. Firstly, applying the optimal excitation to identify the moment of inertia. Secondly, based on the identified moment of inertia, estimating the mass and center of mass.
Another novel idea of this article is a new method to find the optimal excitation to identify the mass and center of mass. A new scalar variable κ is proposed as the benchmark index (different from the condition-number) to identify the optimal excitation.
When we get the moment of inertia, we use this equation Where If measurement noise is relatively small, the least squares algorithm is used to get the best estimate of from equation (110).
If measurement noise is relatively large, then we can used the least squares algorithm to get x from equation (110).
Where ∆A m,N is measurement noise of acceleration matrix ∆A m,N . Then the estimation error is , the optimal excitation is the trajectory to minimize κ. Because is dependent on angular velocity and angular acceleration. So dynamical equation have to be involved in the algorithm. Then If dynamical equation is involved in the process, it is difficult to solve the problem. To avoid that, another method is proposed.
Adopting the numerical solution method The Runge-Kutta method is: Where f is derivative function, h is integration step. If we adopt Euler method, then From equation (117) From equation (118) From equation (101) )︃ Where Ωc = J −1 Tc + J −1 Hw J −1 Hw When there are N sets of data, then we can used the least squares algorithm to get the best estimatê The identification error is , the optimal trajectory is the trajectory that minimize κ.
Where κ is related to the output angular momentum and output torque of momentum wheels. Let

Identification of mass and center of mass using thrusters
Identification error is , in order to get high accuracy estimation, on one hand, minimizing the value of ⃦ ⃦ Ψ d,N ⃦ ⃦ , e.g. adopt the input shaper to eliminate the measurement noise, on the other hand, minimizing the value of, e.g. choose the optimal excitation.
When considering the measurement noise of ω equation (66) Where ∆Ω N denoting the measurement noise of ω.

Identification of inertia parameters of rigid spacecraft
In the simulation, Kun et al. (2017) already used the optimal excitation to identify the moment of inertia of spacecraft. Based on the identified moment of inertia of spacecraft. An optimal excitation is generated to estimate the mass and center of mass. The values of κ are derived from two excitations (optimal excitation and sinusoidal excitation). When the value of κ at the last exciting duration is close to zero, it means that the optimal excitation has strong insensitivity to the measurement noise and disturbance torque. Figure 2 shows values of κ calculated from the optimal excitation and the sinusoidal excitation. Obviously, the κ  value of the optimal excitation is getting smaller than the κ value of the sinusoidal excitation. Table 1 shows the actual mass parameters and the identified mass parameter.
Simulation results show that the proposed method can be used to get high-accuracy identification results on orbit.

Identification of parameters of flexible spacecraft
The method mentioned above can also be applied to the identification of flexible body. As shown in Figure 2, OX b Y b Z b is the body coordinate system of the hub. L is the solar panels, EJ is the so called flexural rigidity, ρ is line density of the solar panels.   Table 2 show the physical parameter. Table 3 show the related parameters of the model.
J v is the true value of the moment of inertia,Ĵ v is the identified moment of inertia, defining FM:  Figure 3 shows values of κ of 4 excitations. The corresponding values of FM are displayed in Figure 4. We can get the optimal excitation design for identifying inertia parameters of spacecraft in Kun et al. (2017). The excitation of case 1 is the optimal excitation and it has the minimum value of κ. Also the identification error is the minimum. For detailed analysis, please reference Kun et al. (2017).   The relative errors are too much at time 42.5 s. But at time 100 s, the relative error has decreased sharply to meet the accuracy requirements of the identification.

Conclusion
In this paper, a two-step identification method is proposed, both the moment of inertia and the mass properties are identified. A new index parameter which is different from the commonly used condition number is first defined for designing the optimal excitation. A method is introduced based on the least squares algorithm. A detailed mathematical equations and numerical analysis are proposed. Simulation results show that the method can identify the mass parameters of spacecraft. There are 3 novelty idea: 1. In my previously published work (Kun et al. 2017), only the moment of inertia is identified. In this paper, both the moment of inertia and mass property are identified. 2. A new index parameter which is different from the commonly used condition number is first defined for designing the optimal excitation. 3. Not only rigid spacecraft, the identification method could be used to identify the parameters of flexible spacecraft.