Improve the identification convergence speed with decimation

In addition to the SSA-SSI-COV algorithm which pre-process the noisy data and allows to realize all identifiable modes, the convergence of time domain system identification algorithms can be easily improved by “decimation”, which is a two-step process: low-pass anti-aliasing filter and downsampling, thus, the choice of a proper sampling rate is an important issue in system identification [29,66]. The Nyquist frequency gives the lower bound for the sampling rate meanwhile the upper bound is determined by the numerical instability due to truncation and round-off errors in a digital computer. It is shown in [67] that the poles of the transfer function of a discrete-time system approach to one on the unit circle as the sampling interval becomes very small, which leads to numerical instability in the computation especially when the data is noise contaminated. In [68] a sampling rate between 10 times to 50 times the closed loop system bandwidth is recommended for the digital implementation of feedback systems.

The Canton Tower is a very flexible and long period structure with about 0.09 Hz

56

for its first mode, and there are 9 modes below 1 Hz. It is possible that a sampling rate of 50 Hz is too high for a fast convergence of these modes below 1 Hz through the use of stabilization diagram especially for the 1^st mode.

Three downsampling factors were chosen for comparison: 2, 5, and 10, which will reduce the original sampling rate (50 Hz) to 25 Hz, 10 Hz and 5 Hz respectively. To avoid aliasing, a low-pass filter must be used as an anti-aliasing filter to reduce the bandwidth of signal before it is downsampled, i.e., the signal are firstly low-passed by a Butterworth 10^th order filter with cut-off frequencies of 12.5 Hz, 5 Hz and 2.5 Hz corresponding to the respective Nyquist frequency. The same 400 seconds data length is downsampled for identification, thus, there are only 2000 points available for the case with sampling rate of 5 Hz, 4000 points for 10 Hz and 10000 points for 25 Hz. The system order is defined to be 60 when the data is downsampled to 5 Hz, and 90 when it is downsampled to 25 Hz and 10 Hz. The outcome stabilization diagram is shown in Figure 4-25 for the frequency ranging from 0 to 1 Hz. It is clear from a) that with the data downsampled to 25 Hz, it is still costly to reveal a stable diagram for the 1^st mode;

however, a stable diagram appears for all modes with few block rows when the data is downsampled to b) 10 Hz and c) 5 Hz, although for this latter not all modes were identified. Figure 4-26 shows the stability diagram for the frequency range between 1 to 6 Hz. Examples of complex modes shapes identified with different sampling rates are shown in Figure 4-27 for modes 1~10 (the wind modes are not shown), the same as those obtained in section 4.2 and 4.3 by SSI-COV and SSA-SSI-COV under a sampling rate of 50 Hz, these complex modes appear almost in a straight line, thus, the mode shape quality is not affected by a proper decimation.

Finally, the identified modal frequencies and damping ratios are summarized in

57

Table 4-2, one can note that the identified frequencies are the same as that obtained with the original sampling rate, only the identified damping ratio for 1^st mode (about 1%) is much lower than that identified through SSI-COV and SSA-SSI-COV using the original sampling rate of 50 Hz. From this analysis, it is concluded that a sampling rate of 10 Hz is suitable for identification and monitoring of Canton Tower which implies a significant reduction in computation effort. If only lower modes are interested, the sampling rate can be even reduced to 5 Hz and the convergence is even faster.

58

Chapter 5

Recursive Stochastic Subspace Identification algorithms

Different from the off-line analysis, the on-line system identification and damage detection, based on vibration data measured from the structural health monitoring system has received considerable attention recently. In this section, the Covariance-driven Recursive Stochastic Subspace Identification algorithm (RSSI-COV) is discussed, for later, to be used to estimate system modal parameters from the response measurements of time-varying systems. To consider the noise contaminated data, a recursive pre-processing technique called recursive singular spectrum analysis technique (rSSA) is introduced to enhance the accuracy and stability in the online tracking capability.

5.1 Recursive Covariance-driven Stochastic Subspace Identification algorithm (RSSI-COV)

For online application of SSI-COV, instead of arranging the block covariances in the so-called Toeplitz matrix as shown in (2.27), these must adopt the form of a Hankel Covariance matrix, which is the way it is outlined in NExT-ERA [10]. From the same stochastic properties shown in (2.25), the Hankel Covariance matrix has the following

59 Hankel Covariance matrix can be constructed by arranging the output measurement data vectors as follows: matrix. One can find easily that the latter can be built by the following expression:

∑

obtained matrices A and C, _p^~can be set to 1 without further influence on the obtained models. Since the order of the Hankel Covariance matrix is “i” with data length N, then, the Hankel Covariance matrix (which is a square matrix since the same order is used for rows and columns) can be calculated as the summation over N-2i+1 rank-one matrices formed by y_k⁺y_k^−T , and later a new incoming data point N+1 will be converted into a

60

new rank-one matrix adding to the existing Hankel Covariance matrix.

There are three possibilities to formulate the adaptive Hankel matrix for RSSI-COV [17]:

1. Exponential forgetting: The old data is multiplied by a forgetting factor µ when a new data is added. incorporated (updating). Let the window length be L:

[ ] ^∑

3. Combined approach: It is the application of forgetting factor in updating as well as in downdating.

where the superscript N-2i+1 above the forgetting factor µ is because HN is formed by the summation over N-2i+1 rank-one matrices formed by y_k⁺y_k^−T, and since the Hankel Covariance matrix is updated first, the rank-one matrix “ N L i

(

N L i

)

^T

−− + +− +2 y 2

y ” to be downdated

should be multiplied by the corresponding forgetting factor to the power of N-2i+1.

The software library LAPACK is used in MATLAB to compute SVD of the Hankel Covariance matrix in an off-line manner, which uses the classical algorithms

61

like Householder reflections and QR algorithms [18], but it is very costly in the computation effort since it takes O(ab²) floating point operations to compute SVD[45], where a is the number of rows and b is number of columns of the matrix. It is not suitable for online applications. The need of a recursive fashion to update SVD was firstly found in the field of sensor array signal processing, in which the subspace estimation plays an important role and the online tracking of the direction of arrival (DOA) of the plane waves is the main issue. As a consequence, a new approach called Projection Approximation Subspace Tracking (PAST) was initially developed by Bin Yang [56], who takes advantage of a mathematical lemma to find the required column subspace as an unconstrained optimization problem. Later the algorithm is modified to its Extended Instrumental Variable version (EIV-PAST) by Gustafsson [20], which is a suitable algorithm for the structure of SSI-COV [17]. It is important to mention that PAST is not the unique algorithm to track the time varying subspace, it is classified in the category of the unconstrained quadratic optimization problem [34]. In the following, a brief description of PAST, its extension to EIV-PAST and implementation to recursive SSI-COV will be described.

5.1.1 Projection Approximation Subspace Tracking (PAST)

The PAST is originally a fast dominant-eigenvectors updating algorithm which is based on the following unconstrained cost function:

( )

^{W E}

( )

^{t WWH}

( )

^{t Tr}

{

^{Cz WHCzW WHCzWWHW}

}

V = z − z ² = −2 + (5.7)

where z

( )

t ∈ℜ^m×1 is a random vector, _⋅ denotes the Euclidean vector norm. E{·} and Tr{·} are the expectation and trace operator respectively. C_z is the signal covariance matrix defined as Cz ^{= E[z}( ) ( )^{t zT t ]}. The superscript ^H denotes Hermitian transpose, W

62

is a matrix with suitable dimensions. In our case, the covariance of vibration signals is real, the desired column subspace W is also real, therefore the Hermitian transpose can be treated as the usual matrix transpose. The cost function is known as Yang’s criterion, and there is a mathematical statement with respect to the matrix W:

Theorem 1[20]: The matrix W is a stationary point of V(W) if and only if T

U

W = , whereU ∈ ℜ ^m×r has orthonormal columns and contains any r distinct eigenvectors of Cz. All stationary points of V(W) are saddle points, except when U contains the r dominant eigenvectors, i.e., U = U′. In this case, V(W) attains the global minimum. Here,

T ∈ ℜ

^r×ris an arbitrary unitary matrix. Proof can be found in [56].

Avoiding the cumbersome mathematical derivations of the proof, more intuitively, the cost function is trying to minimize the error between the “transformed” or “filtered”

signal ^z

( )

^t =^U′U′Tz

( )

^t and z(t), through a transformation matrix composed by the dominant eigenvectors of the signal covariance matrix Cz, i.e., U′U′^T . Usually the computed dominant eigenvectors are orthonormal, thus, U′U′^T is actually an orthogonal projection matrix as it is defined in linear algebra, and ^z

( )

^{t is the}

projection of the noisy signal into the signal subspace, thus, noise will be filtered out by this orthogonal projection. In fact, this “subspace filtering” concept was the original idea of Bin Yang to eliminate the measurement noise, which later will be implemented in recursive SSA.

Based on Theorem 1, instead of solving the Eigen-Decomposition problem through classical approaches, the unconstrained cost function only tries to update the dominant eigenvectors of the signal covariance matrix Cz. For time-varying systems, the dominant

63

where µdenotes the forgetting factor used in the summation to give different weights to the random vector z(k) in the summation, and the expectation in (5.7) is replaced by the summation. Also the definition of C_z has to be replaced by

( ) ∑ ( ) ( )

equation (5.7), after expand the cost function, it is a fourth order matrix equation to be solved. To adapt the solution to a Recursive Least Square (RLS) approach, an slow varying comparing to the sampling rate of data point.

With this assumption, since the dominant subspace W^H(k −¹) is already known

which is a typical optimization function in Least Square problems and can be minimized by:

64

where Czh is the covariance matrix formed by z(k) and h(k):

( )

^t

( ) ( )

^{k k} zh

( ) ( ) ( )

^{t t H t}

Thus, when there is a new incoming data at instant t, the matrix inversion lemma can be applied to (5.14) and the well-known RLS algorithm can be easily derived for updating W(t).

5.1.2 Instrumental Variable Projection Approximation Subspace Tracking (IV-PAST)

Since PAST by Bin Yang was formulated originally to treat antenna signals corrupted with additive noise through a subspace approach [56]. However, from the derivations and assumptions shown in section 2: in output-only SSI, the input source

uk

B is unknown, which together with system noise are assumed to be a stationary and spatially white noise, i.e., instead of a simple additive noise, it is rather the source of the system response. For this type of noise, it was proved in [43] that the normal least square formulation will lead to a biased solution and it is not appropriate to handle this type of problem; instead, an Instrumental Variable (IV) approach must be used. Since the instruments must be uncorrelated with noise, usually the same output measurement but with a time lag “i” can be chosen for this purpose. The modification of PAST to IV-PAST and later Extended IV-PAST was proposed by Gustafsson [20].

By introducing the instrument ^ξ

( )

^{t ∈ℜm×1}, the least square solution to the

65

objective function in (5.11) becomes the following equation:

( ) ( ) ( ) ( ) ( )

^]

( ) ( ) ( )

⁰

ξ weighted by the forgetting factor, can replace by the respective cross-covariance matrices Czξ

( )

^t and Chξ

( )

^t , which are defined and updated by:

The same matrix inversion lemma can be applied to (5.18) and the same RLS-like algorithm can be formulated for IV-PAST.

5.1.3 Extended Instrumental Variable Projection Approximation Subspace Tracking (EIV-PAST)

The inconvenience one may found in IV-PAST is that, to let the inverse of the cross-covariance matrix exist, the length of the instrument vector must have the same size as the measurement vector z(t). Another problem one may found in practice is that the cross-covariance matrix may be ill-conditioned for which is not amenable the inversion. The Extended Instrumental Variable was derived to solve this problem and make further stable the inversion process. Then, the cost function to be minimized can be replaced by its corresponding EIV formulation:

66

The least square solution of (5.19) is readily found to be:

( ) ^{t = ′} ( ) ^{t =} ( ) ( ) ( ) ( ) ^{t t} [

h

^t

^Th

^t ]

⁻¹

The difference between (5.15) and (5.20) is the effect of the extra-added Frobenius norm, which is able to fulfill the need of Recursive Least Square (RLS) via matrix inversion lemma, to treat the case of non-square matrices due to the use of instruments of different length, or, to obtain a better numerical stability in the recursion. A complete derivation and formulas for the EIV-RLS algorithm can be found in [43].

5.1.4 Adaptation of EIV-PAST to RSSI-COV

Since covariance driven subspace can be considered as an SVD-enhanced Instrumental-Variable (IV) method [37], one may think that the IV-PAST algorithm is suitable to perform the SVD-updating task, or better said, the observability matrix O_i (column subspace of the Hankel covariance matrix) updating task, which is not true.

Let the random vector z(t) in IV-PAST formulation be replaced by the corresponding data vector in (5.2), which is y_k⁺∈ℜ^il×1 ; on the other side, the substitution of the instrument ^ξ

( )

^t is evidently y_k^−T_∈ℜ^1×il. IV Solution to the cost function will become:

67

window formulation (forgetting factor defined as the unit), and

∑ ( )

+

Similar to (5.18), the dominant eigenvectors can be found by

( )

^{t =U′IV}

( )

^{t =Hcovt}

( )

^{Htcov −1}

W (5.22)

However, either PAST or IV-PAST algorithm was derived to update dominant eigenvectors of the covariance matrix, the relationship between Eigen-Decomposition (ED) and SVD of the Hankel covariance matrix is preferred to be revised. The SVD of a Hankel Covariance matrix is defined as:

( )

_^

where U and V are orthonormal matrices, S is a diagonal matrix containing the singular values. But the column subspace U can be also obtained from the ED of the Hankel Covariance matrix multiplied by its transpose:

( ) ( )

¹

cov

covH = USV VS U =USS U = USS U−

H ^{T T T T T T T (5.24)}

From the relationships shown above, the desired observability matrix O_i is the same as the column subspace U1 extracted from Hankel Covariance matrix using SVD.

However, if Theorem 1 is reviewed, after solving (5.21) by least square, the obtained dominant eigenvector W

( )

^t =U′IV

( )

^t is the eigenvector of the Hankel Covariance

Symmetric matrix

68

matrix, it is NOT the desired column subspace. This latter must be computed via ED of the Hankel Covariance matrix multiplied by its transpose as shown in (5.24), i.e., the desired column subspace can be updated as the dominant eigenvectors of the

“Covariance of Hankel Covariance matrix”.

Fortunately there is EIV-PAST whose Frobenius norm is just satisfying this requirement. Again, substituting the random vector z(t) and the instrument ^ξ

( )

^{t by the}

corresponding data vector in (5.2), the objective function of EIV-PAST to be

where the moving window approach is adopted again. Similar to (5.20), the least square solution to (5.25) is the follows:

( ) ^{t = U}

¹

( ) ^{t =} ( ^H

^covt

^H

^{covt T}

)( ^H

^tcov

^H

^{covt T}

)

⁻¹

W

^(5.26)

Comparing (5.22) and (5.26), the first one is computing the dominant eigenvectors of the Hankel Covariance matrix, on the contrary, (5.26) is computing the dominant eigenvectors of H_t^covH_t^covTas that shown in (5.24), i.e., the desired columns subspace U1(t) of the Hankel Covariance matrix.

Hence, the so-called Extended Instrumental Variable Recursive Least Square (EIV-RLS) algorithm can be applied to solve the EIV-PAST problem, which fulfills the SVD-updating requirement of RSSI-COV to track the time-varying subspace U1(t). The explicit formulas to be implemented in RSSI-COV are shown below. Complete derivation of these formulas of EIV-RLS algorithm can be found in [43].

69 updated and downdated using the following algorithm:

Updating:

70

Several numerical studies of time varying systems subjected to white noise excitation were carried out to test the performance of RSSI-COV. Results show that for the SSI algorithms, forgetting factor leads to severe problems in the tracking of modal parameters. Since forgetting factor is introducing different weightings to the data meanwhile Hankel covariance matrix is being formed, the ability of this step to cancel out the random components no longer works. Moving window approach with forgetting factor set to one is adopted in this work since it better matches the assumptions of SSI-COV. The use of a moving window implies the same procedure shown above has to be done twice to complete the subspace updating for each new incoming data: after adding the new incoming data (updating), the oldest data has to be subtracted from the moving window (downdating). The same formulas shown above can be applied for downdating by setting forgetting factor µ equals to one, the data vectors y_t⁺₊₁ and _yt⁻₊₁^T should be replaced by the oldest data vector in the moving window, i.e., if the moving window length is L, these become: y_t⁺_−L+2i and _yt⁻_−L+2i^T. Moreover, several sign changes in the last four formulas have to be introduced for downdating:

Downdating:

71

( )

^{t+1 =}

[

^2Λ

( ) ( ) ( ) ( )

^{t+1 +ψt+1TPt+1 ψt+1}

]

^−1ψ

^{( ) ( )}

^t+1Pt+1*

K

µ

^{* (5.28e)}

( )

^{1 1}

( ) ( ) ( ) ( )

^1*

[

^{1 1 1 1}

] ^{( )}

¹

1t+ =U t+ − vt+ −U t+ ψt+ Kt+

U ^{* (5.28f)}

(

t +¹

)

= ¹2

[

P

(

t +¹

)

^* + P

(

t +¹

) (

^*ψ t +¹

) (

K t +¹

) ]

P µ ^(5.28g)

( ) (

t L

)

^T

t

t₊₁ = H^*₊₁ −h t +1 y⁻_{− +1}

H µ ^(5.28h)

(

^{t L}

)

^T

L t t

t

−− + +− +

+

+₁ = H^* ₁ − y ₁ y ₁

H µ ^(5.28i)

Since the updating task is done for the time-varying column subspace U1(t) from the formulas shown above, the system information can be then extracted from U1(t) as previously discussed in section 2.3, for each time instant.

5.2 Recursive Singular Spectrum Analysis (rSSA)

To be able to apply SSA in online filtering of vibration measurements, an on-line version of the algorithm that describe the current signal structure at each time instant is required. As mentioned before in section 2.7, the first step is to assemble the measurement data in a Hankel data matrix of N data points X(N) as shown in (2.47).

Although a forgetting factor 0<λ<1 can be applied to gives different weights to the data in terms of its age, through several simulation studies have been carried out, one conclude again that for subspace-based algorithms, the moving window approach should be adopted. The number of block rows is i’ and is kept constant meanwhile a new data point is added as a new column appended to the moving window Hankel matrix:

72 i.e., number of block rows of the Hankel data matrix; L’ is the length of moving window, K’=L’–i’+1 is the number of columns, with K’ > i’. For convenience, the subscript notation for the sliding window vector is different than that used in off-line SSA, which becomes now X_{N-L’+i’+j} , j =0,1,2, …, K’.

Similar to that shown in (5.24), by applying SVD to the Hankel data matrix

(

N+¹

)

X , the left singular vectors U can be computed via Eigen-Decomposition (ED) of the covariance of sliding window vectors:

T

where C_SSA(N+1) is the covariance matrix of the sliding window vector X_{N-L’+i’+j}.

To make rSSA algorithm possible, only the left singular vectors corresponding to the non-zero singular values of X(N+1) will be used and updated as on-line filters: since they correspond to the column subspace which span the range of X(N+1), i.e., any sliding window vector X_{N-+L’+i’+j} can be expressed as a linear combination of the computed column subspace from X(N+1).

Since the eigenvectors of the covariance matrix C_SSA(N+1) correspond to the desired column subspace, this is actually a typical rank-two modification of the symmetric eigen-problem, i.e., meanwhile a new data column is appended to the Hankel matrix (i.e., rank-one modification), an old data column is subtracted (i.e., rank-two

73

modification). The above-mentioned PAST algorithm is suitable to be implement to rSSA, because it is able to update in a recursive fashion the dominant eigenvectors of the signal covariance C_SSA(N+1), i.e., the subspace of X(N+1).

Adaptation of PAST to rSSA

From (5.8), the random vector z(k) can be substituted by the sliding window vector X_k:

( )

By introducing the same approximation:

( ) ( )

k

H k

k X

h′ = W −1 (5.32)

the original cost function is converted into a quadratic criterion:

( )

74

When the matrix inversion lemma is applied to (5.34), the well-known RLS algorithm can be easily derived, which is shown in the chart below.

1. From an initial SVD the recursive algorithm can be initialized with U′1

( )

^N , later

h can be updated using the following algorithm:

Updating:

在文檔中 協方差型隨機子空間識別法之應用 (頁 76-0)