*Palle Andersen Rune
Brincker*

*Structural Vibration Solutions A/S *

*NOVI** Research
Park*

*Niels Jernes Vej 10, Sohngaardsholmsvej
57,*

*DK-9220 Aalborg East. **DK-9000
Aalborg*

*Denmark**Denmark*

In the traditional
input-output modal analysis the estimation of modal parameters have been
performed using a somewhat deterministic mathematical framework. One of the
major hurdles for people of this traditional modal community to overcome, when
turning to output-only modal analysis, is the switch of the mathematically
framework. In output-only modal analysis the mathematically framework involves
the use of statistics and introduction of concepts such as optimal prediction,
linear system theory and stochastic processes.

The two general assumptions
made in output-only modal analysis is that the underlying physical system
behaves linearly and time-invariant. The linearity imply that if an input with
a certain amplitude generates an output with a certain amplitude, then an input
with twice the amplitude will generate an output with twice the amplitude as
well. The time-invariance implies that the underlying physical system does not
change in time. One of the typical parametric model structures to use in
output-only modal analysis of linear and time-invariant physical systems is the
stochastic state space system.

_{} (1)

The first part of this model
structure is called the state equation and models the dynamic behavior of the
physical system. The second equation is called the observation or output
equation, since this equation controls which part of the dynamic system that
can be observed in the output of the model. In this model of the physical
system, the measured system response *y** _{t}*
is generated by two stochastic processes

The philosophy is that the
dynamics of the physical system is modeled by the *n*´*n* state matrix ** A**. Given an

The state space model (1) is
only applicable for linear systems that do not have time-varying changes of its
characteristics. However, this is not the only restriction. The only way to
obtain an optimal estimate of a state space model on the basis of measured
system response, is to require that the system response is a realization of a
Gaussian distributed stochastic process that has zero mean.

In other words, in the
applied stochastic framework the system response is modelled by a stochastic
process *y** _{t}* defined
as

_{}
(2)

and the principal assumption
is that the measured system response is a realization of this process. It is
seen that this process is completely described by its covariance function_{}**. **This means that if we can estimate a state space
model having the correct covariance function this model will completely
describe the statistically properties of the system response. An estimated
model fulfilling this is called covariance equivalent. The estimator that can
produce such model is called an optimal estimator.

Since the system response of
the linear state space model is a Gaussian stochastic process it implies that _{}**, **_{}** and **_{}**all are Gaussian stochastic processes as well. Since
the input processes **_{}** **and _{}**are unknown we make the simplest possible assumption
about their statistical properties, which is to assume that they are two
correlated zero-mean Gaussian white noise processes, defined by their
covariance matrices as**

_{}
(3)

The Gaussian stochastic
process describing the state_{}**is also zero-mean and completely described by its
covariance function**

_{}
(4)

Using (1) to (4) the
following relations can be established

_{} (5)

The matrix ** G**
is the covariance between system response

_{}
(6)

There are two additional
system matrices turns out to play an important role

_{}
(7)

These are the extended
observability matrix_{}**and the reversed extended stochastic controllability
matrix **_{}**.**

One of the most important parts of all estimation is the ability to predict the measurements optimally. In output only modal analysis this means to be able to predict the measured system response optimally. An optimal predictor is defined as a predictor that results in a minimum error between the predicted and measured system response. If the system response can be predicted optimally it implies that a model can be estimated in an optimal sense.

Recall
that the state vector ** _{}** completely describes the system dynamics at
time

_{}
(8)

In the
Gaussian case the optimal predictor of ** _{} **is then given by the conditional mean value

_{}
(9)

So, the
optimal predictor of ** _{} **is defined as the mean value of

_{}
(10)

This
error is the part of ** _{} **that cannot be predicted by

In order
to predict the system response a similar conditional mean can be formulated for_{}

_{} (11)

The last
part of this equation is obtained by inserting (1) and assuming that** _{} **and

The two predictors (9) and (11) are related through the so-called Kalman filter for linear and time-invariant systems, see e.g. Goodwin et al. [6]

_{}
(12)

The
matrix ** _{}** is called the non-steady state Kalman gain and

_{}
(13)

The last
of these equations is called the Ricatti equation. The Kalman filter predicts
the state ** _{}**on the basis of the previous
predicted state

_{}
(14)

Given
that the initial state prediction is** _{}**and the initial state prediction
covariance matrix

At start up the Kalman
filter (12) will experience a transient phase where the prediction of the state
will be non-steady. However, if the state matrix ** A** is stable the filter
will enter a steady state as time approach infinity. When this steady state is
reached the covariance matrix of the predicted state vector

_{} (15)

The steady state Kalman gain
is now calculated from

_{}
(16)

The last equation is now
called an algebraic Ricatti equation. Assuming all matrices but ** P**
is known this equation can be solved using eigenvalue decomposition, see Aoki
[2] and Overschee et al [1].

If the last equation in (15)
is rearranged the following state space system is obtained

_{}
(17)

This system is called the
innovation state space system. The major difference between this system and (1)
is that the state vector has been substituted with its prediction, and that the
two input processes of (1) have been converted into one input process – the
innovations. This state space system is widely used as model structure in
output only modal analysis, see e.g. Ljung [3] and Söderström et al. [4].

The Kalman filter defined in
the last section turns out to be the key element in the group of estimation
techniques known as the stochastic subspace techniques.

From (17) it is seen that if
sufficiently many states of (1), let’s say *j*
states, can be predicted, i.e. _{}** **and _{}**, then the **** A** and

_{}
(18)

This is a valid approach
since the innovations are assumed to be Gaussian white noise. Since ** A**
and

So the fundamental problem
to solve in the stochastic subspace identification technique is to extract the
predicted states from the measured data. To show how this is performed,
consider the state space system in (1) and take the conditional expectation on
both sides of both equations to yield

_{} (19)

Now assume that a recursion
is started at time step *q*. Inserting
the first equation in (19) recursively into itself *i* times and finally inserting the result the last of the equations
in (19) leads to the following formulation

_{}
(20)

This equation shows the
relation between the initial predicted state _{}**and the prediction of the free (noise free) response
of the system **_{}**. By stacking ***i*
equations on top of each other the following set of equations are obtained

_{} (21)

By introducing the vector _{}**as the left-hand side and noticing that the first
part of the right-hand side is equal to the extended observability matrix **_{}**we actually obtain the following expression for the
predicted states**

_{}
(22)

The matrix_{}**is actually the pseudo-inverse of **_{}**. This equation shows that if we can estimate **_{}**and **_{}**for several values of ***q*, we can in fact estimate the predicted states for several values
of *q* as well.

In this section we will focus on the estimation of the predicted free
response _{}**. We will
estimate a set of vectors **_{}**and gather
them column by column in a matrix **** O**.

In order to predict the system response a conditional mean similar to (11) can be formulated.

_{}
(23)

This conditional mean is the
prediction of the future system response _{}** **given
the past system response from time *t =
i+q-*1 down to *t = q*. This
conditional expectation is only an approximation of (11) since the conditioning
vector stops a time *t = q* and not *t =* 0. The approximation is only good if
*i* is sufficiently high. For zero-mean
Gaussian stochastic processes this conditional mean can be calculated by, see
e.g. Melsa et al. [5].

_{}
(24)

Since the error _{}**is zero-mean and uncorrelated and is independent of
the conditional mean **_{}**and the conditioning vector **_{}**the conditional mean (24) is also called the
orthogonal projection of the vector **_{}**onto the vector **_{}**.**

In order to estimate all
elements of ** _{}**we need to extend (24) to allow estimation of

_{}
(25)

In the last equation a new *i**p*´*ip*
matrix _{}** is
introduced for simplicity. This matrix is defined as **

_{}
(26)

Incidentally, the matrix _{}** is also
equal to **

_{}
(27)

As seen in (18) we need a
bank of predicted state estimates _{}**for ***q = i*
to *q = i+j-*1 for a sufficiently large
value of *j*. To
estimate these state in one operation based on the approach in (25) we need to
define the following two matrices _{}**and **_{}**as**

_{}
(28)

_{}
(29)

The index *p* in (29) signifies that the matrix
contains system response of the past compared to the system response we are
predicting. Since we assume that the system response is stationary, i.e. that ** _{}**, equation
(25) can easily be extended using (28) and (29) to yield

_{} (30)

With this equation the first
of the two major tasks in the stochastic subspace identification technique has
been fulfilled.

If the extension in (30) is
carried on to (22) we obtain the following relation

_{}
(31)

The matrix_{}** is a bank of predicted states and is defined as**

_{}
(32)

As seen the matrix_{}**only depends on system response and system response covariance,
and can therefore be estimated directly from the measured system response. In
Overschee et al. [1] a method based on the QR decomposition is presented (For
more on the QR decomposition, see e.g. Golub et al. [7]). This method estimates
**_{}** **directly from the
measurements without explicit need of the covariance estimates. By using that
method the stochastic subspace identification techniques can surely be called
data driven identification techniques.

In order to estimate ** A**
and

The only input we have for
the estimation is still only matrix _{}**, i.e. only information related to the system
response. The underlying system that has generated the measured response is
unknown, which means that we do not know the state space dimension of
underlying system. What this means can be seen from equation (31) that defines
the matrix **_{}**as the product **_{}**. The outer dimension of **_{}** and therefore
also of **_{}**is ***ip*´*j*. However, the
question is what the inner dimension of this product is. The inner dimension is
exactly the state space dimension of the underlying system.

So to find _{}**the first task is to determine this dimension.
We determine this dimension from **_{}**by using the Singular Value Decomposition or SVD,
see e.g. Golub et al. However, before taking the SVD we pre- and postmultiply **_{}**with the before mentioned weight matrices **_{}**and **_{}**which are user-defined. Now taking the SVD of the
resulting product yields**

_{}
(33)

Assuming that _{}**has full rank and that the rank of **_{}**is equal to the rank of **_{}**, the dimension of the inner product **_{}**is equal to the number of non-zero singular values,
i.e. number of diagonal elements of **_{}**. From the
last two equations of (33) we see that**_{}**is given by**

_{}
(34)

The non-singular *n*´*n*
matrix ** T** represents an arbitrary similarity transform. This means that
we have determined the extended observability matrix except for an arbitrary
similarity transformation, which merely means that we have no control over the
exact inner rotation of the state space system.

As seen the state space
dimension is determined as the number of diagnonal elements of _{}**, and **_{}**is found on the basis the reduced subspace of **_{}**. For these reasons it is no wonder why the
estimation techniques are called subspace identification techniques.**

Independently of the choice
of weight
matrices _{}**and **_{}**the estimation of the system matrices can be done in
the general way presented in this section. This approach presented here is not
the only one, but in the current context properly the most obvious choice. In
Overschee et al. [1] two other approaches are also described. The estimation
can be divided into three parts. **

Assuming that *N* samples of measured system response
are available the user needs to specify how many block rows *i* the matrix _{}** should have.
As seen from (33) the maximum state space dimension depends on the number of
block rows and will be ***ip*, where *p* is the dimension of the measured
system response vector _{}**. **

It should be remembered that
the maximum state space dimension corresponds to the maximum number of
eigenvalues that can be identified. It should also be remembered that *i* is the prediction horizon and as such
depends on the correlation length of the lowest mode to be identified.

_{}** are the
estimated using (30). However, in order to estimate the matrix**_{}**we also need to estimate the matrix**_{}**since**

_{}
(35)

This can be proven by proper
substitutions in the above equations, see also Overschee et al. [1]. _{}**is estimated
by deleting the first ***p* rows of _{}**.**

Now we have all the
information available that is needed in order to estimate a realization of the
innovation state space system defined in (17). Estimate the predicted states_{}**and **_{}**using (31) and (35), and set up the following matrix
of measured system response**

_{} (36)

Solve the least squares
problem

_{}
(37)

where _{}**is the pseudo inverse of **_{}**. The steady state Kalman gain **** K** is estimated by the
following relations. First estimate the reversed extended stochastic
controllability matrix

_{}
(38)

The covariance matrix ** G**
can then be extracted from the last

_{}
(39)

Estimate the Kalman gain ** K**
in (16) by solving the algebraic Ricatti equation in (16) first. Finally,
estimate the covariance matrix

As mentioned in section 3.2,
the only significant difference between the different stochastic subspace
algorithms is the choice of the weight matrices _{}**and **_{}**. In this chapter we will focus on three algorithms,
the Unweighted Principal Component algorithm, the Principal Component algorithm
and the Canonical Variate Analysis algorithm.**

The Unweighted Principal
Component algorithm is the most simple algorithm to incorporate into the
stochastic subspace frame work. As the name says it is an unweighted approach
which means that both weight matrices equals the unity matrix, see Overschee et
al. [1]

_{}
(40)

The reason is that this
algorithm determines the system order from the left singular vectors ** _{}**of the SVD of the following matrix

_{}
(41)

Since we have chosen the
weight to be unity the covariance of ** _{}**equals

_{}**
**(42)

This show that covariance
(42) is equal to (41) which means that (41) and (42) has the same left singular
vectors. From (34) we see that the extended observability matrix is determined
as

_{} (43)

This algorithm is also known
under the name N4SID.

The Principal Component
algorithm determines the system matrices from the singular values and the left
singular vectors of the matrix_{}**. This means that the singular values and left
singular vectors of **** _{}**must equal the singular values and left singular
vectors of

_{}
(44)

The covariance of ** _{}**now

_{}**
**(45)

This shows that the singular
values and the left singular vectors of (45) and_{}**are equal. From (34) we see that the extended
observability matrix is determined as**

_{}
(46)

Even though it appears that
(46) is equal to (43) they are not. We must remember that the SVD has been
taken on different matrices due to the different weights.

The Canonical Variate
Analysis algorithm, see Akaike 74,75, computes the principal angles and
directions between the row spaces of the matrix of past outputs_{}**and the matrix of future outputs **_{}**. The matrix of past outputs **_{}**is defined in (29), and the matrix **_{}**is defined in a similar manner**

_{}
(47)

The principal angles and
directions between the row spaces of ** _{}**and

_{}
(48)

Comparing with (30) we
obtain the same covariance matrix of (47) and ** _{}**, and therefore also the same principal angles and
directions between the row spaces of

_{}
(49)

The covariance of ** _{}**now

** _{}
**(50)

The covariance of (47) is
given by

_{} (51)

We see that the covariance
matrices of ** _{}**and (48) are equal. In the above it has been
used that

_{} (52)

Since ** _{} **is a
byproduct of the determination of

No matter what subspace
identification algorithm to use there are some common problems related to the
selection of the right state space dimension. The reason is that the singular
values in (33) never drop completely to 0. They merely drop significantly
compared to the largest singular values. So it is necessary provide some extra
tools that can aid in the selection process.

A widely used technique is
to make a frequency versus state space dimension plot and present the natural
frequency estimates of the poles of the estimated state space models.

1
Overschee, P. van & B. De Moor: *Subspace identification for linear systems –
Theory, Implementation, Applications.* Kluwer academic Publishers, ISBN
0-7923-9717-7, 1996.

2
Aoki, M.: *State
Space Modeling of Time Series.* Springer-Verlag, ISBN 0-387-52869-5, 1990.

3
Ljung, L.: *System
Identification – Theory for the user.* Prentice-Hall, ISBN 0-13-881640-9,
1987.

4
Söderström, T. & P. Stoica: *System Identification.* Prentice-Hall, ISBN 0-13-127606-9, 1989.

5
Melsa, J. L. & A.P. Sage: *An Introduction to Probability and Stochastic Processes.*
Prentice-Hall, ISBN 0-13-034850-3, 1973.

6
Goodwin, G.C. & K.S. Sin: *Adaptive Filtering, Prediction and Control.* Prentice-Hall, ISBN
0-13-004069-X, 1984.

7
Golub, G.H. & C.F. Van Loan: *Matrix Computations. *2^{nd} Ed.
The