Time-varying modal parameters identification of a spacecraft with rotating flexible appendage by recursive algorithm

Time-varying modal parameters identification of a spacecraft with rotating flexible appendage by recursive algorithm

Author’s Accepted Manuscript Time-varying modal parameters identification of a spacecraft with rotating flexible appendage by recursive algorithm Zhiy...

2MB Sizes 1 Downloads 51 Views

Author’s Accepted Manuscript Time-varying modal parameters identification of a spacecraft with rotating flexible appendage by recursive algorithm Zhiyu Ni, Ruinan Mu, Guangbin Xun, Zhigang Wu www.elsevier.com

PII: DOI: Reference:

S0094-5765(15)00363-X http://dx.doi.org/10.1016/j.actaastro.2015.10.001 AA5564

To appear in: Acta Astronautica Received date: 13 May 2015 Revised date: 17 August 2015 Accepted date: 1 October 2015 Cite this article as: Zhiyu Ni, Ruinan Mu, Guangbin Xun and Zhigang Wu, Timevarying modal parameters identification of a spacecraft with rotating flexible appendage by recursive algorithm, Acta Astronautica, http://dx.doi.org/10.1016/j.actaastro.2015.10.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Time-varying modal parameters identification of a spacecraft with rotating flexible appendage by recursive algorithm Zhiyu Ni a, *, Ruinan Mu a, c, Guangbin Xun b, Zhigang Wuc a

State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, China b

Department of Engineering Mechanics, Dalian University of Technology, China

c

School of Aeronautics and Astronautics, Dalian University of Technology, China

Corresponding author: Tel.: +86 0411 84706772; Fax: +86 0411 84706772. E-mail address: [email protected] (Zhiyu Ni)

Abstract: The rotation of spacecraft flexible appendage may cause changes in modal parameters. For this time-varying system, the computation cost of the frequently-used singular value decomposition (SVD) identification method is high. Some control problems, such as the self-adaptive control, need the latest modal parameters to update the controller parameters in time. In this paper, the projection approximation subspace tracking (PAST) recursive algorithm is applied as an alternative method to identify the time-varying modal parameters. This method avoids the SVD by signal subspace projection and improves the computational efficiency. To verify the ability of this recursive algorithm in spacecraft modal parameters identification, a spacecraft model with rapid rotational appendage, Soil Moisture Active/Passive (SMAP) satellite, is established, and the time-varying modal parameters of the satellite are identified recursively by designing the input and output signals. The results illustrate that this recursive algorithm can obtain the modal parameters in the high signal noise ratio (SNR) and it has better computational efficiency than the SVD method. Moreover, to improve the identification precision of this recursive algorithm in the low SNR, the wavelet de-noising technology is used to decrease the effect of noises. Keywords: Modal Parameter Identification, Time-varying System, Rotating Flexible Appendage, Recursive Algorithm, Spacecraft

1 Introduction For the spacecraft identification problem, most structures are still mainly processed as the time-invariant system, such as the in-orbit identification experiments of the Hubble space telescope [1] and the Engineering Test Satellite VI (ETS-VI) satellite [2]. If the structure is slow time-varying, caused by rotation of the solar panels in the identification experiment, the solar panels can be temporarily locked at certain angles in the identification process due to the relatively slow speed of solar panel

rotation [2-3]. Therefore, the time-varying system can be degraded into a time-invariant system. However, from the perspective of actual operations, we hope that in-orbit identification can be implemented without affecting the normal work conditions of the spacecraft, such as by using the input and output (I-O) signals (control torque signal/satellite attitude signal/appendage vibration signal, etc.) during spacecraft attitude adjustment to perform the corresponding in-orbit identification experiment. In this case, the identification of the time-varying parameters may be considered. In addition, because the rotation appendage carried by the spacecraft becomes increasingly large, the huge appendage may lead to the structural parameters varying with time [4]; for example, the solar panel rotation of the Engineering Test Satellite VIII (ETS-VIII) can cause a maximum 25% change in the structure parameter [5]. This type of effect is sometimes impossible to ignore. In this paper, a spacecraft with large rotational appendage, the Soil Moisture Active/Passive (SMAP) satellite, is selected as our research model. The SMAP mission is an observation mission proposed by the U.S. Jet Propulsion Laboratory to survey global soil moisture and freeze-thaw states. The launch of this satellite can help to improve the accuracy of weather forecasting, flood warning, and drought monitoring [6-8]. The Delta II rocket sent the satellite into a near-circular orbit at perigee 660 km and apogee 685 km on January 31, 2015, U.S. Eastern Standard Time. To maximize the scope of earth observation during the satellite’s in-orbit operation and obtain high-resolution global soil moisture maps, the on-board large deployable antenna can rotate at a speed of 4 s/r after deployment. Because the satellite structural configuration is changing during in-orbit operation, which has led to changes in the corresponding system modal parameters, the satellite can be studied as a time-varying system. The identification methods of modal parameters for some spacecrafts have been performed in other studies [1,9-10]. In these methods, the eigensystem realization algorithm (ERA) method has been successfully applied to the parameter identification experiments of time-invariant systems many times [11-13], such as the identification of the Galileo spacecraft and the Hubble space telescope by Juang [14] and the ETS-VI/VIII satellite in-orbit parameter identification [3-5,15]. The ERA method is a typical system realization method. By constructing a Hankel matrix corresponding to the I-O data and performing singular value decomposition (SVD) on the Hankel matrix at each time step, the system observability matrix is established to determine the state-space model parameters {A, B, C} , and the corresponding modal parameters are obtained from the state-space model [13]. If the system is time-varying, then the repeated-experiments methods can be used, such as the pseudo modal method

proposed by Liu [16] and the time-varying ERA (TV-ERA) method proposed by Majji [17]. These methods can be considered the improvement of time invariant identification methods, and the multiple sets of I-O data are applied in these repeated-experiments methods. Because the in-orbit spacecraft often has periodic characteristics, the repeated-experiments methods are reasonable for identifying the time-varying spacecraft system. However, because these methods require the SVD, the computational cost of the identification procedure is still very high. Some control problems, such as the self-adaptive control, often need to obtain the latest dynamical parameters for updating the controller parameters in real time. In addition, the identified system parameters can help to track and monitor the on-orbit working condition of spacecraft. Consequently, a faster identification method is required. To avoid the SVD and implement the on-line identification, a recursive subspace method called projection approximation subspace tracking (PAST) was proposed by Yang [18]. The PAST algorithm is based on the theory of signal subspace projection. By tracking the signal subspace matrix from the I-O data, this method avoids the SVD computation and thus improves the computational efficiency. The PAST method has already been applied to the modal parameter identification experiments of some time invariant/variant mechanical systems, but few studies have investigated in spacecraft parameter identification. The primary purpose of this paper is to study the recursive identification problem of the time-varying modal parameters of the SMAP satellite caused by the rotation of large flexible antenna. With the control torque during satellite attitude motion as the input signal, the attitude angle/angular velocity and appendages vibration signals as the output signal, the PAST algorithm is used to provide an alternative identification approach for the commonly used pseudo modal method to identify the time-varying modal parameters recursively. Because the SVD has been avoided and the required data decrease in the recursive identification procedure, the computational speed has been improved. In the numerical simulation, the time-varying frequencies are identified by designing the I-O signals. A comparison of the frequency values under different signal noise ratios (SNRs) proves that the recursive algorithm can identify the time-varying modal parameters of the SMAP satellite effectively under the higher SNR. Besides, when the SNR is lower, the wavelet de-noising technology is employed to process the output signal with noise to improve the identification precision. Finally, the numerical results also illustrate that the recursive algorithm has faster computation speeds than the frequently-used SVD methods for time-varying spacecraft system.

The contents of this paper are organized as follows: Section 2 reviews the spacecraft modelling principle for the rigid-flexible coupling structure, and the state-space equation of the SMAP satellite is provided. Section 3 briefly introduces the basic procedures of the PAST recursive algorithm, and the satellite time-varying modal parameters are identified using PAST method. By designing the I-O signals, the results of the numerical simulation verify the effectiveness of the recursive algorithm in Section 4. Finally, some conclusions are presented in Section 5.

2 Establishment of the Satellite State-space Model This section describes the establishment of the SMAP satellite model. First, the satellite structure is simplified, and the coordinate system is defined. Next, the problem of finite element modal analysis on appendages (such as the antenna and solar panels) is considered. Then, the system’s rigid-flexible coupling dynamical equation is established [19]. Finally, the system state-space equation is developed by considering the selection of the satellite I-O signals. 2.1 Description and Simplification of the Satellite Structure The SMAP satellite has a deployable antenna with a diameter of 6 m. The antenna consists of the front/rear net, tension tie, rim trusses, and metallic fittings. The antenna is connected to the body of the satellite by two fixed support rods. To increase the ground scanning surface area during in-orbit operations, the structure formed by the antenna and the two connection rods can rotate quickly at a speed of 4 s/r. The definition for the satellite coordinate system is shown in Figure 1. The Earth’s centre of mass is selected as the origin O of the inertial coordinate system O-XYZ. The X-axis is on the equatorial plane and points to the spring equinox. The Z-axis is perpendicular to the equatorial plane and coincides with the Earth’s rotational axis. The Y-axis is defined according to the right-hand rule. The entire satellite’s centre of mass is selected as the origin ob of the satellite coordinate system ob-xbybzb. The zb-axis points in the opposite direction along the centre of the Earth. The xb-axis points in the in-orbit flight direction of the satellite. The yb-axis is similarly defined according to the right-hand rule. The rigid body hinge attachment point P is selected as the origin oa of the appendage coordinate system oa-xayaza. The za-axis is along the direction of the connection rod. The structure formed by the antenna and the two connection rods rotate around the zb-axis during operation (the antenna rotational angle is defined as  ); therefore, the xa-axis and ya-axis in the appendage coordinate system will change as the

antenna rotates. To simplify the discussion, this paper does not consider the effects of orbital factors and gravity gradients on spacecraft structural identification.

Figure 1 Simplified SMAP satellite model and its coordinate definition

The following definition and simplification are used in modelling: (1) With the exception of the two solar panels that extend outward from the satellite body, the rest are directly simplified into one cuboid rigid body model (abbreviated as central rigid body b, see Figure 2). (2) To ensure the stability of the entire satellite’s attitude during the antenna rotation, because the centre of mass of the antenna (abbreviated as antenna a) is not on the zb-axis of the satellite coordinate system, the mass of the two connection rods (abbreviated as rods PQ and QR) between antenna a and central rigid body b are adjusted so that the centre of mass of the entire satellite remains on the zb-axis after the composite structure c is deployed (the composite structure c consisting of the antenna and the two connection rods). (3) Point P is the hinge attachment point between the composite structure c and the central satellite body b . The angle between the connection rod PQ and the zb-axis is defined as  . The smaller angle between connection rods QR and PQ is defined as  . (4) Because the mass of the two solar panels (abbreviated as panels l and r) on the satellite is minute in comparison with the mass of the entire satellite, these two panels are assumed to have no effect on the position of the entire satellite’s centre of mass. That is, the entire satellite’s centre of mass remains unchanged during the rotation of composite structure c .

Figure 2 Defined plane configuration of the SMAP satellite

2.2 Finite Element Analysis of the Satellite Appendages Modelling and modal analysis of appendages such as the antenna and the solar panels are required to solve the corresponding appendage parameters such as the modal vibration shape, frequency, and radius vector of the node. These parameters will be applied for the rigid-flexible coupling dynamical equation of the satellite. The established antenna model is shown in Figure 3.

Figure 3 Finite element model of the deployable antenna a

2.3 Rigid-flexible Coupling Dynamical Model of the Satellite Through the aforementioned model simplification, this satellite can be considered as a rigid-flexible coupling structure with a rotatable appendage, where the appendages of the satellite have the composite structure c and panels l and r. This paper does not consider the satellite’s translational motion but is only concerned with the rotational motion. The attitude angle is defined as ψ  [ψx , ψ y , ψz ] . It is supposed that the angular velocity ψ of the satellite is very small, namely, ψ  0 , and the structure c is rotating with a uniform rotational speed of  , that is, the angular

acceleration  is   0 . When a 3 1 dimensional control torque u(t ) is applied to the satellite body, the rigid-flexible coupling equation of the satellite can be written as:

J (t )ψ   Fj (t )η j  u(t )

(1)

j

FjT (t )ψ  ηj  2ζ j Ω j ηji  Ω2j ηj  0 ( j  c, l , r )

(2)

where η j is defined as the arbitrary first p orders modal coordinate of appendage j. ζ j and Ω j are the damping ratio and natural frequency of appendage j, respectively, by the finite element analysis in Section 2.2. Matrix F is the coupling coefficient for the appendage vibration and spacecraft rotation.

J is the rotation inertia matrix of the entire satellite. See the references [20-21] for more details on spacecraft modelling theory. Define vector ξ  ψ T

ηcT

ηlT

T

ηrT  ; then, Equations (1) and (2) can be written in the form of

the following dynamical equation:

M (t )ξ  Eξ  Kξ  Lu(t )

(3)

where

 J (t ) Fc (t ) Fl (t ) Fr (t )  033  F T (t )   I 2 c Ωc , E M (t )   cT  Fl (t )   I  T   I   Fr (t )  033    2 Ω c  , L   I 33  K   0  2   Ωl     Ωr2  

2 l Ωl

     2 r Ωr 

and I represents a unit matrix. Define the n1 dimensional state vector x (t )   ξ T

T

ξ T  , and

rewrite Equation (3) in the form of a state equation:

x(t )  A(t ) x(t )  B(t )u(t )

(4)

where A(t ) and B(t ) are the n  n and n  r time-varying system and input matrices, respectively.

0  A(t )   1   M (t ) K

  0  , B(t )   1    M (t ) E   M (t ) L I

1

T

The m1 dimensional output vector y (t ) is selected as y(t )  ψ T , sT , ψ T , sT  , where s is the vibration displacement of the appendage. Then, the measurement equation can be written as:

H y(t )  Cx (t )   

 ξ    H   ξ 

(5)

where C (k ) is the m  n output matrix, and matrix H is

I  H=   ( j  c, l , r ) Φ j  Φ j ( j  c, l , r ) represents the corresponding modal matrix of appendage j. In the numerical simulation, the modal matrix is obtained by finite element modelling. In the actual in-orbit identification, the system output signals are measured by using the gyros on the centre rigid body and the sensors on the appendages.

3 Recursive Identification of Satellite Time-varying Modal Parameters The PAST recursive algorithm is introduced in this section. Then, the identification of the time-varying modal parameters of the satellite is presented using the recursive method. 3.1 PAST Recursive Algorithm Consider adding a noise term to Equations (4) and (5). The state-space equation of the discrete time-varying system of the satellite can be described as:

x(k  1)  A(k ) x(k )  B(k )u(k )  wn (k )

(6)

y(k )  C (k ) x(k )  vn (k )

(7)

where wn (k ) and vn (k ) represent the process and measurement noise of the system, respectively. The I-O data sequences at time step k can be described as: T

u(k )  uT (k ) uT (k  1)

uT (k  M  1) 

y(k )   y T (k )

y T (k  M  1) 

y T (k  1)

T

(8) (9)

Equations (6) and (7) can be written as the following generalized equation:

y(k )  Γ (k ) x(k )  Δ(k )u(k )

(10)

where

C (k )   C (k  1) A(k ) Γ (k )     C (k  M  1) A(k  M  2) 0   C ( k  1) B(k ) Δ(k )     C (k  M  1) A(k  M  2)

     A(k )  0

B(k ) C (k  M  1) A(k  M  2)

B(k  1)

     0

where the matrix Γ (k ) is the observability matrix of the time-varying system. The procedures of the PAST algorithm are briefly introduced here. First, define a random data vector p and consider the following expectation function relationship:

J (W (k ))  E p(k )  Γ (k ) Γ H (k ) p(k )

2

(11)

where the superscript “ H ” represents the Hermitian transposition and the symbol “ E ” represents expectation. Suppose the matrix Γ (k ) is full rank, then the expectation for the covariance matrix of vector p(k ) can be written as Rpp (k )  E[ p(k ) pT (k )] . If the expectation function in Equation (11) is replaced by a scalar function, then Equation (11) can be written as: k

J (W (k ))   k i p(i )  Γ (k ) Γ H (k ) p(i )

2

(12)

i 1

where  is a forgetting factor. Then, Rpp (k )  E[ p(k ) pT (k )] can also be written as: k

Rpp    k i p(i) p H (i)

(13)

i 1

In Equation (12), the main goal of the PAST algorithm is to construct the projection relationship of data vector p(k ) onto the matrix Γ (k ) : q(k )  Γ H (k  1) p(k )

(14)

where vector q(k ) is the subspace projection of p(k ) . Substitute Equation (14) into Equation (12), and the modified value of Equation (12) is as follows: k

J ( Γ (k ))   k i p(i )  Γ (k )q(i )

2

(15)

i 1

To minimize the value of Equation (15), the following must be satisfied:

Γ (k )  Rpq (k ) Rqq1 (k )

(16)

Apply the matrix inversion lemma to Equation (16) and obtain the recursive form of the matrix

Γ (k ) , as follows: Γ (k )  Γ (k  1)  e(k ) g H (k )

(17)

where e(k )  p(k )  Γ (k  1)q(k ) , g(k )  Z (k  1)q(k )(  qH (k ) Z (k  1)q(k )) 1 , and Z (k ) is the recursive variable in the iterative process after the initial value is selected. Table 1 shows the procedure of PAST algorithm. See the references [18] by Yang for the specific derivation of the PAST algorithm.

Table 1 PAST algorithm Initialization:

 In  Γ (0)    , Z (0)  I n 0( mM n )n  For k=1, 2, 3, … :

Cost

q(k )  Γ H (k  1) p(k )

mMn

ν(k )  Z (k  1)q(k )

g(k )  ν(k )(  qH (k )v(k ))1

n2 2n

Z (k )  ( Z (k  1)  g(k )v H (k )) / 

2n2

e(k )  p(k )  Γ (k  1)q(k )

mMn

Γ (k )  Γ (k  1)  e(k ) g (k )

mMn

H

The computation complexity of PAST is 3mMn+O(n2) flops per iteration [18], whereas the complexities of the frequently-used identification method, which is based on SVD, is O(n3) at least. Therefore, the computational cost of PAST is lower than the SVD method. Particularly, when the order of spacecraft system is high, the advantage of computational efficiency for the recursive algorithm is obvious. In PAST algorithm, we need to apply the I-O data u(k ) and y(k ) to construct the data vector p(k ) . For this purpose, a data pre-processing way to construct and update the vector p(k ) is

introduced in Appendix A. 3.2 Identification of the Time-varying Modal Parameters After the matrix Γ (k ) is determined using the recursive algorithm, the discrete system matrix A(k ) can be obtained by the following equation:

A(k )  [ Γ1 (k  1)]† [ Γ2 (k )]

(18)

where Γ1 (k  1) and Γ 2 (k ) are the first m  (M  1) rows of Γ (k  1) and the last m  (M  1) rows of Γ (k ) , respectively. According to the pseudo modal analysis method by Liu [16], it is supposed that the eigenvalue of the time-varying system at each time instant is unchanged when the sampling time selected is sufficiently small. When the system matrix A(k ) is obtained from Equation (18), eigenvalue decomposition is performed to the identified matrix A(k ) at each time instant k , as follows: A(k )  Σ (k ) Λ(k ) Σ 1 (k )

where Λ(k )  diag(1 (k ),

,n (k )) is the diagonal eigenvalue matrix; i (k ) (i  1, 2,

(19)

, n) is the

time-varying conjugate complex eigenvalues; and Σ (k ) is the time-varying eigenvectors matrix. The ith pseudo eigenvalue is i (k )  exp( i (k )t  ji (k )t ) , j  1 , where  i (k ) and i (k ) are referred to as the ith system pseudo-damping ratio and pseudo-damped natural frequency, respectively, and t is the sampling time. By using the eigenvalue decomposition method to compute the pseudo natural frequency at each time instant k, the time-varying modal parameters of the satellite can be obtained. 3.3 Summary of the Identification Procedures The procedures of the PAST algorithm to identify the time-varying modal parameters are summarized below: Step 1 Using the I-O data sequences u(k ) and y(k ) in Equations (8) and (9), The data vector

p(k ) is determined recursively by the data pre-processing method described in Appendix A from Equations (A.5), (A.10), (A.12) and (A.13). Step 2 After the data vector p(k ) is determined, the PAST recursive algorithm in Section 3.1 is applied to identify the observability matrix Γ (k ) recursively by Equation (17). Step 3 From the expression form of the matrix Γ (k ) in Equation (10), the time-varying matrix

A(k ) is computed at each time instant by Equation (18). The time-varying parameters of the satellite can be determined by the eigenvalue decomposition of Equation (19).

4 Numerical Simulation A numerical simulation is applied in this section to identify the time-varying modal parameters of the SMAP satellite by the recursive algorithm. The numerical model of the satellite is established, and

the I-O signals are designed for identification. Then, the time-varying frequencies are obtained. 4.1 Satellite Model Parameters The size of the satellite model is illustrated in Figure 4. The values of some simulation parameters in Equations (1) and (2) are given in Appendix B. The rotational speed of composite structure c is

  90 deg/s. The first 4 orders and the first 8 orders of modal frequencies for the solar panels l / r and antenna a are respectively selected by finite element analysis. Therefore, the order n of state-space Equation (4) is n=(3+4×2+8)×2=38, the damping ratio of antenna a for each order is selected as 0.05, and the damping ratio of panels l and r for each order is selected as 0.1. The first 8 orders of the frequency values for antenna a are listed in Table 2. In the established satellite model, the 1st to 3rd orders of frequencies denote the rigid body motion of the satellite and are equal to zero, and the 4th-7th orders of the rigid-flexible coupling frequency values are shown in Figure 5. The standard zero-mean Gaussian random noise is the process noise and measurement noise, respectively, and the SNRs for process and measurement noise are selected as 50 dB and 40 dB. We only consider the open-loop system identification in this paper, and the situation of closed-loop identification will be discussed in a subsequent paper.

Figure 4 Size of the SMAP satellite model (unit: m)

Table 2 The first 8 orders of frequencies for antenna a, obtained by finite element analysis System orders

Frequency of the antenna (Hz)

1

2.4297

2

3.1047

3

8.3022

4

11.4596

5

14.2968

6

19.2184

7

22.9616

8

29.7757

Figure 5 The 4th to 7th original frequencies of the established satellite model

4.2 Design of the Input Signal and Collection of the Output Signal To simulate the loading and unloading process of the control torque signal produced by satellite reaction wheels, the input signal designed for the simulation is as shown in Figure 6.

Figure 6 Designed input torque signal of the satellite

In addition, the situation of the placement of the sensors on the appendages to collect the output signal is also considered. The out-of-plane vibration response signals of some nodes on each appendage are selected as the output signals (see Figure 7). The output matrix C of Equation (5) is constructed by extracting the corresponding elements from the modal matrix Φ j ( j  c, l , r ) . The initial condition

x (0) is zero. The response results of attitude angle ψ (t ) and attitude angular velocity ψ (t ) are as shown in Figure 8.

Figure 7 Sensor placements for each appendage

Figure 8 Open-loop responses of the attitude angles and attitude angular velocities of the SMAP satellite

4.3 Parameter Identification Results The following parameters are provided in the numerical simulation: forgetting factor  =0.98, sampling time t =0.01 s, and the Hankel parameter M=20 to ensure that the rank of the Hankel matrix is greater than the system order. The time-varying frequency of the satellite is identified using the proposed recursive method and pseudo modal method [16]. The results of the 4th to 7th orders of the rigid-flexible coupling frequency values are shown in Figure 9, and Table 3 shows the average relative errors of each frequency by these two methods. Figure 9 and Table 3 illustrate that these two methods can effectively identify the system frequencies. Although the relative error for recursive algorithm is slightly higher than the SVD method, the maximum error remains below 2%.

Figure 9 Identification results of the 4th to 7th frequencies of the satellite calculated by the recursive method: (a) the 4th frequency; (b) the 5th frequency; (c) the 6th frequency; (d) the 7th frequency

Table 3 Error comparison of frequencies computed by the recursive and pseudo modal methods Frequency

Recursive algorithm

Pseudo modal method

4

1.5594%

1.0398%

5

1.0981%

1.0296%

6

1.4946%

1.4080%

7

1.0512%

1.0247%

Figure 10 shows the identification results of the 4th frequency using the recursive algorithm when the measurement noise is given as different SNRs. The error comparisons of the frequencies are presented in Table 3. As shown in Table 3, the results demonstrate that the computation accuracy is better when SNR≥30, and the relative error of the frequencies increases substantially when SNR=20 (or when the SNR is below 20). The simulation results prove that a lower SNR signal will dramatically affect the identified results because the SVD is avoided in the recursive algorithm.

Figure 10 Comparison of the 4th frequencies by recursive algorithm with different SNRs

Table 3 Error comparison of frequencies by recursive algorithm for different SNRs SNR (dB)

4

5

6

7

50

1.5594%

1.0981%

1.4946%

1.0512%

40

1.7109%

2.0256%

1.8578%

2.0246%

30

2.1987%

2.2484%

2.4784%

2.3543%

20

4.1449%

4.9258%

5.1278%

4.4857%

To decrease the effect of measurement noise, the MATLAB toolbox “wavelet” is applied to process the output signals y(t ) with a lower SNR before the identification. Figure 11 (a) and (c) show the original output signals ψ x and ψ x respectively for SNR=20. And the corresponding de-noising signals are shown in Figure 11 (b) and (d). The identified results after de-noising are shown in Figure 12 and Table 4. It can be seen that the identification precision in low SNR is improved obviously after the wavelet de-noising.

Figure 11 De-noising results of the attitude angles ψ x and attitude angular velocities ψ x using the wavelet toolbox (SNR=20)

Figure 12 Comparison of the 4th frequencies by recursive algorithm using the wavelet de-noising technology (SNR=20) Table 4 Error comparison of frequencies by recursive algorithm using the wavelet de-noising technology (SNR=20) Frequency

Before de-noising

After de-noising

4

4.1449%

2.1622%

5

4.9258%

2.6578%

6

5.1278%

3.0635%

7

4.4857%

2.5674%

Finally, the computational efficiencies of the two methods are compared. When the dimension of the Hankel matrices for the two methods is the same, 30 simulations are performed, and the average computation times are listed in Table 5. The simulation results also illustrate the conclusion of the

computation complexity in Section 3.1: because SVD is required in the pseudo modal method, the computation time is found to be much larger than that of the recursive algorithm. The results show that the recursive algorithm has a higher computational efficiency. Particularly, when the system order n is higher, the difference of computational time between these two methods is obvious.

Table 5 Comparison of the average computation time with M block rows of the Hankel matrix in MATLAB (unit: sec) Identification method

M=10

M=15

M=20

M=30

Recursive algorithm

17.43

28.56

43.64

68.89

Pseudo modal method

31.56

50.03

76.80

113.64

5 Conclusion In this paper, we attempt to provide an alternative recursive algorithm that differs from the common SVD methods [1, 13, 16-17] to identify the time-varying modal parameters of the spacecraft with rotational appendages. Because the recursive algorithm does not need to perform SVD at each time step and raises the computation speed, the latest identified results might be used for the on-line identification or the self-adaptive control system design in the future. The identification results of the numerical simulation demonstrate the validity of the recursive algorithm in identifying the time-varying modal parameters using the control torque signal and attitude/vibration signals. The relative error for recursive algorithm is only slightly higher than the SVD method. In addition, the simulation results also show that the recursive algorithm has a higher computational efficiency compared with the SVD-based methods, particularly when the system order is high. However, certain improvements are required for the recursive algorithm: if the SNR is too low in the identification (when SNR≤20, as shown in the example), then the results in Figure 10 and Table 3 illustrate that the identification precision becomes poor, and the method is easily disturbed by noise. For the recursive algorithm, because the SVD is avoided in determining the observability matrix, the noises have significant influence for the identified results. Thus, the signal de-noising technology is used in this paper to decrease the effect of noises as much as possible. Nevertheless, how to improve the ability of noise immunity for the recursive algorithm itself that needs to be studied further.

Appendix A: A data pre-processing way of I-O data for PAST algorithm For the establishment of the input vector p in PAST algorithm, there are several methods. Here we introduce a way to construct the vector p by the URV matrix decomposition. The detailed procedures of this method are described in reference [22]. First, the Hankel matrices corresponding to the I-O data at time step k can be written as:

u(2)  u(1)  u(2) u(3) U (k )      u( M ) u( M  1)

  y (1)    , Y (k )   y (2)     u(k  M  1)   y(M ) u(k )

y (2)

u(k  1)

y (3) y ( M  1)

y(k )  y (k  1)  (A.1)   y (k  M  1) 

where parameter M must satisfy the rank of the Hankel matrix to be greater than the system order n . Define the Hankel matrix at time step k  1 as:

U (k  1)  U (k ) u(k  1) , Y (k  1)  Y (k )

y(k  1)

(A.2)

The matrix U  (k ) that is normal to U (k ) is defined as: U  (k )  I  U T (k )[U (k )U T (k )]1U (k )

(A.3)

The recursive update of U (k ) may be expressed as

[U (k )U T (k )]1  [U ( k  1)U T ( k  1)  u( k ) uT ( k )]1  [U (k  1)U T (k  1)]1 

δ(k )uT (k )[U (k  1)U T (k  1)]1 1   (k )

(A.4)

where the parameters δ(k ) and  (k ) are defined as δ(k )  [U (k  1)U T (k  1)]1 u(k ) ,  (k )  uT (k )δ(k )

(A.5)

U T (k  1)  U T (k )[U (k )U T (k )]1U ( k )   T   u (k )  [(U † (k  1))T   (k ) T (k )U ( k  1) / (1  δ( k )) ( k ) / (1  δ( k ))]

(A.6)

leads to

where the symbol “ † ” represents the pseudo inverse. Therefore,

 I  U  (k  1) 0 U T (k )[U ( k )U T ( k )]U ( k )  I  U  ( k )    0 0  U T ( k ) δ( k ) δT ( k )U ( k  1) / (1   ( k )) U T ( k  1) δ( k ) / (1   ( k ))    δT (k )U ( k  1) / (1   ( k ))  ( k ) / (1   ( k ))   yields

(A.7)

U  (k  1) 0 U T (k  1)δ(k )  T 1 U  (k )        δ (k )U (k  1) 1 (1/ 1   ( k )) 0 0 1 1   (k )   

(A.8)

where the notation “   ” represents the double floor. Therefore,

Y ( k )U  ( k )  Y ( k  1) 

1 1   (k )

  U ( k  1) 0 y( k )     0 0  

 U T (k  1)δ(k )  T    δ (k )U (k  1) 1 (1/ 1   ( k ))  1   

(A.9)

 Y (k  1)U  (k  1) 0  z ( k )  wT ( k ) 1 where

z(k )  ([Y (k  1)U † (k  1)]u(k )  y(k )) / 1   (k )

(A.10)

wT (k )  δT (k )U T (k  1)  uT (k )[U (k  1)U T (k  1)]1U T (k  1) / 1   (k )

(A.11)

and the recursive forms for Y (k  1)U † (k  1) and U (k  1)U T (k  1) in Equations (A.10) and (A.11) are as follows:

[U (k )U T (k )]1  [U (k  1)U T (k  1)]1  δ(k )δT (k ) / 1   (k )

(A.12)

[Y (k )U † (k )]  [Y (k  1)U † (k  1)]  p(k )δT (k )

(A.13)

The vector z (k ) , which is derived from Equation (A.10), is what we need to the PAST algorithm as the vector p(k ) . When the data vector p(k ) is obtained, the matrix Γ (k ) can be derived by the recursive method in Section 3.1. See the reference [22] for a specific detailed description of the data pre-processing method.

Appendix B: Some selected simulation parameters for SMAP satellite The rotation inertia matrix J of the spacecraft is determined as follows:

J (t )  J b  mb rb rb  Ti JiTi T   mi rr (i  a, l , r, PQ, QR) i i

(B.1)

where m (  b, a, l , r, PQ, QR) is the mass (unit: kg); J  (  b, a, l , r, PQ, QR) denotes the rotation inertia (unit: kg m2); r (  b, a, l , r, PQ, QR) is the skew-symmetric matrix of the position vector from the centre of mass of the total system to the centre of mass of each component (unit: m) and

T (  a, l , r, PQ, QR) denotes the coordinate transformation matrix in the satellite coordinate system. mb =797.6, ma =64, ml = mr =31.2, mPQ =6, mQR =14

0 0 0 0   258.5553  14.1440    Jb   0 258.5553 0 3.7482 0   , Jl  J r   0  0 0 0 132.9333 0 10.4042   0  0 0  0 0   288 0 0.8650 10.5467      J a   0 288 0  , J PQ   0 0.04 0  , J QR   0 0.0933 0   0  0  0 0 576  0 0.8650 0 10.5467 0.3843 0 0.6343 1  0.6343 1   0  0  0      rb  -0.3843 0 0  , rl   0.6343 0 0.48 , rr   0.6343 0 0.48  0  1  1 0 0  0.48 0  0.48 0  4.6657 1.05 0.8835 0.4979   0  0    ra   4.6657 0 0  , rPQ   0.8835 0 0   1.05  0.4979 0 0  0 0  2.7301 1.4705  0 cos  rQR   2.7301 0 0  , Ta  TQR   sin   1.4705  0 0 0 

TPQ

cos    sin   0

 sin  cos  0

 sin  cos 0

0  1 0 0  0 sin  1  0  cos 

0  cos  sin  

0  1 0 0  1 0 0     0 0 sin(   ) cos(   )  , Tl  Tr  0 1 0 0 0 1  1  0  cos(   ) sin(   ) 

Acknowledgments The first author would like to thank the team of State Key Laboratory of Structural Analysis for Industrial Equipment for their helpful suggestions to improve the structure model. This work was supported by the National Natural Science Foundation of China (11372056).

References [ 1 ] J.N. Juang, M. Phan, L.G. Horta, et al., Identification of observer/Kalman filter Markov parameters-theory and experiments, J Guid. Control Dyn. 16(2) (1993) 320-329. [2] I. Yamaguchi, T. Kida, K. Komatsu, et al., ETS-VI on-orbit system identification experiments, JSME Int. J. 40(4) (1997) 623-629. [3] S. Adachi, I. Yamaguchi, T. Kida, et al., On-orbit system identification experiments on Engineering Test Satellite-VI, Control Eng. Pract. 7(7) (1999) 831-841. [4] A. Meguro, K. Shintate, M. Usui, et al., In-orbit deployment characteristics of large deployable

antenna reflector onboard Engineering Test Satellite VIII, Acta Astronaut. 65(9) (2009) 1306-1316. [5] T. Nagashio, T. Kida, T. Ohtani, et al., Design and implementation of robust symmetric attitude controller for ETS-VIII spacecraft, Control Eng. Pract. 18(12) (2010) 1440-1451. [6] D. Entekhabi, E. Njoku, P. O'Neill, et al., The soil moisture active/passive mission (SMAP), in: International Geoscience and Remote Sensing Symposium, 2008. [7] J.Y. Liu, Space-based large spinning sensor pointing and control design and its application to NASA’s SMAP spacecraft, in: AIAA Guidance, Navigation, and Control Conference, National Harbor, Maryland, 2014. [8] D. Entekhabi, E. Njoku, P. O'Neill, et al., The soil moisture active passive (SMAP) mission, Proc. of the IEEE. 98(5) (2010) 704-716. [9] I. Yamaguchi and T. Kida, System identification experiments of large space structures, Trans. of the Society of Instrument and Control Engineers. E-1(1) (2001) 312-321. [10] L. H. Geng, D. Y. Xiao, Q. Wang, et al., Attitude-control model identification of on-orbit satellites actuated by reaction wheels, Acta Astronaut. 66(5) (2010) 714-721. [11] I. Yamaguchi, T. Kida, T. Kasai, Experimental demonstration of LSS system identification by eigensystem realization algorithm, Proceedings of American Control Conference (Seattle, WA), American Automatic Control Council, Evanston, IL, 1995. [12] I. Yamaguchi, T. Kasai, H. Igawa, Multi-input multi-output system identification using impulse responses, J. Space Eng. 1(1) (2008) 1-11. [13] J.E. Cooper and J.R. Wright, Spacecraft in-orbit identification using eigensystem realization methods, J Guid. Control Dyn. 15(2) (1992) 352-359. [14] R.S. Pappa and J.N. Juang, Galileo spacecraft modal identification using an eigensystem realization algorithm, J. Astronaut. Sci. 33(1) (1985) 95-118. [15] T. Kasai, I. Yamaguchi, H. Igawa, et al., On-orbit system identification experiments of the Engineering Test Satellite-VIII, Trans. of Space Tech. Japan 7 (2009) 79-84. [16] K.F. Liu and L. Deng, Identification of pseudo-natural frequencies of an axially moving cantilever beam using a subspace-based algorithm, Mech. Syst. Signal Process. 20(1) (2006) 94-113. [17] M. Majji, J.N. Juang, J.L. Junkins, Time-varying eigensystem realization algorithm, J Guid. Control Dyn. 33(1) (2010) 13-28. [18] B. Yang, Projection approximation subspace tracking, IEEE Trans. Signal Process. 43(1) (1995)

95-107. [19] R.R. Craig, Coupling of substructures for dynamic analyses: an overview, Collection of the 41st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference and Exhibit, AIAA Paper, Atlanta, USA, April 2000. [20] P.C. Hughes, Modal identities for elastic bodies, with application to vehicle dynamics and control, J Appl. Mech. 47(1) (1980) 177-184. [21] H.B. Hablani, Constrained and unconstrained modes-Some modeling aspects of flexible spacecraft, J Guid. Control Dyn. 5(2) (1982) 164-173. [22] F. Tasker, A. Bosse and S. Fisher, Real-time modal parameter estimation using subspace methods: applications, Mech. Syst. Signal Process. 12(6) (1998) 797-808.

Highlights

The rotation of the large flexible appendage of spacecraft is considered. The time-varying modal parameters of this spacecraft structure are identified. Apply recursive algorithm without SVD to determine the time-varying modal parameters. The recursive algorithm has faster computational speed than common SVD-based methods. The wavelet de-noising is employed to improve identification precision for lower SNR.