- Email: [email protected]
Multivariable Impulse Response Estimation via Correlation Analysis and Its Application to Automated System Identification
input and the output variables is a classical method to determine whether two variables are related by a transfer function (Box and Jenkins, 1976) . To extend this method to multivariable systems, one must assume that input and measured disturbance sequences in each channel are statistically independent. In process industrial applications, this assumption is usually too restricting for a number of reasons. First, it is usually not realistic in practice to assume that input and measured disturbance sequences are independent. Secondly, this requirement may exclude orthogonal multifrequency signals such as the Schroeder-phased input (Rivera et al. , 1997) which are desirable in a process enviroment due to their "plant-friendly" character. In this paper, we propose a method based on the singular value decomposition of the information matrix from the correlation analysis to determine the potential existence of a transfer function relationship, in spite of correlation between measurements.
cussed in Section 4. Section 5 contains an example based on industrially-provided data. Section 6 provides a brief summary and some conclusions.
2. CORRELATION ANALYSIS OF MULTIVARlABLE SYSTEMS Consider a LTI, MIMO system with ni inputs, no outputs and nd measured disturbances: y(t) = p(q)u(t)
+ pd(q)d(t) + 1I(t)
(1)
where p and Pd are transfer function matrices of proper dimensions, 11 is a zero-mean unmeasured disturbance / noise vector. Sequences u(t), d(t) and 1I(t) are assumed to be quasi-stationary (Ljung, 1987). We shall allow u(t) and/ or d( t) in each channel to be correlated to each other. But the unmeasured disturbance 1I(t) is assumed to be independent to both d(t) and u(t) . Given measurements of sequences {u(t)} , {d(t)} and {y(t)} , the objective is to estimate P and Pd.
Clearly, many different identification methods can be applied to estimate system transfer functions. The quality of final estimates is determined by the choice of design variables involved in an identification process; these include the design of input perturbation, selection of model structures and selection of model parameter estimation algorithms etc. In the context of prediction error methods, the problem of how those design variables influence the model performance is well understood (Ljung, 1987), although in an asymptotic sense. But how to choose those design variables in practice is still an art and often demands very good a priori system information. Because reliable a priori system information is not always available and is supposed to be supplied by the identification process itself, the identification process is often an iterative procedure. In this context, totally automated, "single pass" identification may be an unrealistic expectation. But the incentive to reduce iterations, particularly as it relates to experimental testing, is clear. As such, nonparametric methods such as correlation analysis can serve a useful purpose in identification process monitoring, to assist in determining the "goodness" of a data set during identification testing, and to provide necessary system information to gUide the choice of design variables in the parametric estimation step.
In equation (1), entry (ij) of transfer function matrices p and Pd can be characterized by their finite impulse responses (FIR): k
k
. . = ~9(ih)q-1
Pt]1
~
1
j(i12 ) -I Pd i12 - ~ ~ 1 q -
,
1=0
(2)
1=0
i = 1, . .. , no; j1 = 1, . . . , n i ; j2 = 1, ... ,nd
The system (1) can then be represented as:
(3) j=11=0
In the correlation analysis framework, multiplying both sides of equation (3) by Uj(t - s) and dj(t - s) respectively, and taking the expectation operation gives n i
' UjYi (S) = L
k
iW L9i )'ujU,.,(s -l)+
w=11=0 nd
L
The paper is organized as follows . In Section 2 correlation analysis of multivariable systems and its properties are discussed. In section 3 the classical variable selection method and its limitations are explained. A transfer function detection/variable selection method based on the singular value decomposition of the information matrix is proposed. The recursive estimation method and the online experiment monitoring problem is dis-
k
(4)
Lf/iW ),u,d,.,(s -l)
w=11=0
j
= 1, . .. , ni ; i = 1, .. . , no ; S = 0, .. . , k n,
,djYi(S) = L
k
iW L9i ),dju",(S -l)+
w=11=0 nd k
L
Lf?w)'d,d,u(S -l)
w=11=0
1400
(5)
j = 1, ... ,nd; i = 1, ... , no; s = 0, .. . ,k
obtained for the estimator (7) ; and (b) the problem of selecting proper model structures encountered in parametric methods is avoided. The maximum lag k in estimator (7) is the only parameter which must be specified by users . Hence (7) is very appropriate as a preliminary step providing necessary a priori system knowledge to a parametric estimation method. For most parametric estimation methods, having an estimate of the system step response available a priori is most helpful, which can be obtained straightly by accumulation of proper parts of O. Furthermore, by rearranging 0, the Hankel matrix of the system can be obtained and from the Hankel matrix, a state-space model can be generated using a realization method. The procedure is one of subspace identification methods (Viberg, 1994).
where '"Y denotes the auto or cross correlation function between corresponding sequences. In matrix form , equations (4) and (5) are given by: '"Yi =
where
'"Yi ,
'"Yi
fo Oi,
(6)
i = 1, ... ,no
Oi are vectors representing =
['"YUIYi (0) ,··" '"YUIYi (k) , 1'u2Yi (0),···, '"Yu ni Yi (k),
1'd,Yi (0), . .. , '"YdndYi (k )]T,
Oi -- [90(il) , ... , 9k(il) , 90(i2) ,"', 9k(ini) ' j(il) f,(i2) . .. j(ind)]T f,o(il) ' • .. ' k '0 , , k ' and f
Other important features of the estimator (7) include:
is the information matrix
0
[1] The auto- and cross-correlation functions satisfy '"Yxx(k) = '"Yxx( -k) and 1'xy(k) = '"Yyx( -k). It follows that f ii is a Toeplitz matrix and f ij = f~i ' Hence f 0 and f are symmetric matrices and can be generated effectively. Also some efficient algorithms can be sought to solve the estimator (7). [2] A recursive algorithm can be easily built based on the estimator (6) . An attractive application of the recursive estimator is for on-line monitoring of the identification experiment execution. The problem will be further discussed in Section 4.
fo = [fij](ni+nd)x(ni+nd) with
'"Yxy(O) 1'xy( -1) f
ij
... '"Yxy( -k) ::: ~xy( -k + 1)
~XY(l) ~xY(O)
=
[
1
'"Yxy(k) 1'xy(k -1) .. . 1'xy(O)
and x
= {
y= {
if i :::; ni d i - ni otherwise
Ui
if j :::; ~i d j - ni otherwIse
Uj
It is noted that in practice, all the correlation functions in (7) must be estimated from an estimator. To use the conventional correlation function estimator (Bendat and Piersol, 1986), some standard data pretreatment procedures are necessary. For example, the measurement sequences {u(t)}, {d(t)} and {y(t)} should be differenced to remove possible trends (unstationarity) . Also to avoid ill-conditioning of the f matrix, all the measurement sequences should be properly scaled.
Equation (6) can be further put into a whole matrix form as: '"Y =
f 0,
(7)
with
TT
l' = ['"Yl '"Y2
T]T '"Yn o '
and
f
=
f01
1 is an identity matrix of proper dimension. Equation (7) is a [no(k+1)(ni+nd)] x [no(k+1)(ni+nd)] square linear equation set. It forms the correlation analysis estimator for a multivariable system.
3. APPLICATION TO INPUT jOUTPUT VARIABLE SELECTION In chemical engineering applications there are often many seemingly related measurements such as the composition, temperature, pressure and flowrate etc. But only some of them are really related by transfer functions and some are just measurement noise or process disturbance. To realize the identification process automation, it is necessary to establish some criterion to detect which measurements are really related by transfer functions. In other words, the task of transfer function detection is to determine whether the estimated parameter vector 0 do represent impulse responses among different in-
When all the input and measured disturbance sequences have persistent excitation of at least order k, the matrix f is guaranteed to be nonsingular and a unique solution to 0 exists. In the above problem setting, the estimate 0 is consistent. But due to the large amount of parameters involved with FIR models, it is usually associated with a large variance error, compared to the other parametric estimation methods. Two major advantages using the correlation analysis method are: (a) a simple closed-form solution can be
1401
put/measured disturbance and output channels or some of it just represent the nOise/disturbance. For a SISO system, the problem can be examined via the following null hypothesis test: the output sequence y is uncorrelated with the input sequence u . Under this hypothesis, we have (Bartlett, 1955)
var[-yuy(l)] ~ O"~O";(N _l)-I
lated to their variance through the following expressions respectively:
(8)
where the input sequence u should be prewhitened using a prewhitening filter and 0"2 denotes the variance. Strictly speaking, relation (8) only applies to the case where u is Gaussian distributed. But it is a common practice to use the relation for input sequences of other probabilic distributions. Because usually, N » l , (8) is sometimes approximated as:
Substituting (12) and (13) into (11), the following bounds on the standard deviation of 8 are obtained:
(9) where 'Yi is the ith element of the vector 'Y which can represents 'YUjYk or 'YdjYk for some j and k. Under the hypothesis that a output k is uncorrelated with a input or a measured disturbance j, varbd is approximated by:
Hence, with 99% confidence level, 'Yuy should stay in the range set by the following expression: -2.58
(N _ l) ~ y~
'Yuy () l
~ 2.58
(N _ l) y~
(10)
(15) when all the inputs and output i are scaled by their variance. To use relation (15) , all the measurement sequences should be prefiltered by a prewhitening filter of sufficiently high order. Relations (14) and (15) can be used as a criterion to detect the existence of transfer functions . If the estimated Bi from (6) goes beyond the bounds set by (14), existence of a transfer function between the input j and the output i is suggested.
If 'Yuy goes out of the range in (10), a transfer function between the output y and the input u is suggested. For a MIMO system, if sequences {u(t)}, {d(t)} and {v(t)} are uncorrelated to each other, the above method can be directly applied to each input or measured disturbance/output channel. But assuming that {u(t)} in each channel is statistically uncorrelated to each other will be too restricting in practice. For example, some processfriendly input sequences such as many multi-frequency signals will be excluded in consideration. Furthermore, it is usually not so realistic to assume that {d( t)} and { d( t )} / {u (t)} are uncorrelated to each other. When the input sequences and/or the measured disturbance sequences are correlated, using criterion (10) will lead to incorrect conclusions. For example, assume that Ui and Uj are correlated, and there exists a transfer function between Yk and Ui but not between Yk and Uj . Then amplitude of 'YUjYk is determined by 'YU j Ui and does not necessarily stay in a band set by (10) .
4. APPLICATION TO IDENTIFICATION EXPERIMENT MONITORlNG The correlation analysis estimator (6) forms a simple linear equation set. A recursive algorithm can be readily built. Denote
cp(t)=[uI(t-1) . .. uI(t-k)U2(t-1) . .. uni(t-k) dl (t-1) . .. dl(t-k) .. . dnd(t-k)]T
Another application of the correlation analysis estimator (6) we consider in this section is to establish a transfer function detection criterion. The criterion is a MIMO system extension to the singular value decomposition method of SISO systems discussed by Rivera (Rivera, 1992) . From equation (6) , we have Q:2(fo)1I8ill~ ~ lbi l1 2 ~ 0'2(fo)IIBill~
(6) can then be computed recursively as
i = 1, ... , no
L(t) _ P(t - l)cp(t) -1/a+cpT(t)P(t-l)cp(t)'
(11)
where O'(fo) and Q:(fo) are the maximum and the minimum singular values of ro · 118ill~ and Ibill~ can be re-
P(t)
= P(t _
) 17
1) _ P(t - l)cp(t) cpT(t)p(t - 1) (18)
l/a
1402
(
+ cpT(t)P(t - l)cp(t)
cussed in Section 4. Section 5 contains an example based on industrially-provided data. Section 6 provides a brief summary and some conclusions.
2. CORRELATION ANALYSIS OF MULTIVARlABLE SYSTEMS Consider a LTI, MIMO system with ni inputs, no outputs and nd measured disturbances: y(t) = p(q)u(t)
+ pd(q)d(t) + 1I(t)
(1)
where p and Pd are transfer function matrices of proper dimensions, 11 is a zero-mean unmeasured disturbance / noise vector. Sequences u(t), d(t) and 1I(t) are assumed to be quasi-stationary (Ljung, 1987). We shall allow u(t) and/ or d( t) in each channel to be correlated to each other. But the unmeasured disturbance 1I(t) is assumed to be independent to both d(t) and u(t) . Given measurements of sequences {u(t)} , {d(t)} and {y(t)} , the objective is to estimate P and Pd.
Clearly, many different identification methods can be applied to estimate system transfer functions. The quality of final estimates is determined by the choice of design variables involved in an identification process; these include the design of input perturbation, selection of model structures and selection of model parameter estimation algorithms etc. In the context of prediction error methods, the problem of how those design variables influence the model performance is well understood (Ljung, 1987), although in an asymptotic sense. But how to choose those design variables in practice is still an art and often demands very good a priori system information. Because reliable a priori system information is not always available and is supposed to be supplied by the identification process itself, the identification process is often an iterative procedure. In this context, totally automated, "single pass" identification may be an unrealistic expectation. But the incentive to reduce iterations, particularly as it relates to experimental testing, is clear. As such, nonparametric methods such as correlation analysis can serve a useful purpose in identification process monitoring, to assist in determining the "goodness" of a data set during identification testing, and to provide necessary system information to gUide the choice of design variables in the parametric estimation step.
In equation (1), entry (ij) of transfer function matrices p and Pd can be characterized by their finite impulse responses (FIR): k
k
. . = ~9(ih)q-1
Pt]1
~
1
j(i12 ) -I Pd i12 - ~ ~ 1 q -
,
1=0
(2)
1=0
i = 1, . .. , no; j1 = 1, . . . , n i ; j2 = 1, ... ,nd
The system (1) can then be represented as:
(3) j=11=0
In the correlation analysis framework, multiplying both sides of equation (3) by Uj(t - s) and dj(t - s) respectively, and taking the expectation operation gives n i
' UjYi (S) = L
k
iW L9i )'ujU,.,(s -l)+
w=11=0 nd
L
The paper is organized as follows . In Section 2 correlation analysis of multivariable systems and its properties are discussed. In section 3 the classical variable selection method and its limitations are explained. A transfer function detection/variable selection method based on the singular value decomposition of the information matrix is proposed. The recursive estimation method and the online experiment monitoring problem is dis-
k
(4)
Lf/iW ),u,d,.,(s -l)
w=11=0
j
= 1, . .. , ni ; i = 1, .. . , no ; S = 0, .. . , k n,
,djYi(S) = L
k
iW L9i ),dju",(S -l)+
w=11=0 nd k
L
Lf?w)'d,d,u(S -l)
w=11=0
1400
(5)
j = 1, ... ,nd; i = 1, ... , no; s = 0, .. . ,k
obtained for the estimator (7) ; and (b) the problem of selecting proper model structures encountered in parametric methods is avoided. The maximum lag k in estimator (7) is the only parameter which must be specified by users . Hence (7) is very appropriate as a preliminary step providing necessary a priori system knowledge to a parametric estimation method. For most parametric estimation methods, having an estimate of the system step response available a priori is most helpful, which can be obtained straightly by accumulation of proper parts of O. Furthermore, by rearranging 0, the Hankel matrix of the system can be obtained and from the Hankel matrix, a state-space model can be generated using a realization method. The procedure is one of subspace identification methods (Viberg, 1994).
where '"Y denotes the auto or cross correlation function between corresponding sequences. In matrix form , equations (4) and (5) are given by: '"Yi =
where
'"Yi ,
'"Yi
fo Oi,
(6)
i = 1, ... ,no
Oi are vectors representing =
['"YUIYi (0) ,··" '"YUIYi (k) , 1'u2Yi (0),···, '"Yu ni Yi (k),
1'd,Yi (0), . .. , '"YdndYi (k )]T,
Oi -- [90(il) , ... , 9k(il) , 90(i2) ,"', 9k(ini) ' j(il) f,(i2) . .. j(ind)]T f,o(il) ' • .. ' k '0 , , k ' and f
Other important features of the estimator (7) include:
is the information matrix
0
[1] The auto- and cross-correlation functions satisfy '"Yxx(k) = '"Yxx( -k) and 1'xy(k) = '"Yyx( -k). It follows that f ii is a Toeplitz matrix and f ij = f~i ' Hence f 0 and f are symmetric matrices and can be generated effectively. Also some efficient algorithms can be sought to solve the estimator (7). [2] A recursive algorithm can be easily built based on the estimator (6) . An attractive application of the recursive estimator is for on-line monitoring of the identification experiment execution. The problem will be further discussed in Section 4.
fo = [fij](ni+nd)x(ni+nd) with
'"Yxy(O) 1'xy( -1) f
ij
... '"Yxy( -k) ::: ~xy( -k + 1)
~XY(l) ~xY(O)
=
[
1
'"Yxy(k) 1'xy(k -1) .. . 1'xy(O)
and x
= {
y= {
if i :::; ni d i - ni otherwise
Ui
if j :::; ~i d j - ni otherwIse
Uj
It is noted that in practice, all the correlation functions in (7) must be estimated from an estimator. To use the conventional correlation function estimator (Bendat and Piersol, 1986), some standard data pretreatment procedures are necessary. For example, the measurement sequences {u(t)}, {d(t)} and {y(t)} should be differenced to remove possible trends (unstationarity) . Also to avoid ill-conditioning of the f matrix, all the measurement sequences should be properly scaled.
Equation (6) can be further put into a whole matrix form as: '"Y =
f 0,
(7)
with
TT
l' = ['"Yl '"Y2
T]T '"Yn o '
and
f
=
f01
1 is an identity matrix of proper dimension. Equation (7) is a [no(k+1)(ni+nd)] x [no(k+1)(ni+nd)] square linear equation set. It forms the correlation analysis estimator for a multivariable system.
3. APPLICATION TO INPUT jOUTPUT VARIABLE SELECTION In chemical engineering applications there are often many seemingly related measurements such as the composition, temperature, pressure and flowrate etc. But only some of them are really related by transfer functions and some are just measurement noise or process disturbance. To realize the identification process automation, it is necessary to establish some criterion to detect which measurements are really related by transfer functions. In other words, the task of transfer function detection is to determine whether the estimated parameter vector 0 do represent impulse responses among different in-
When all the input and measured disturbance sequences have persistent excitation of at least order k, the matrix f is guaranteed to be nonsingular and a unique solution to 0 exists. In the above problem setting, the estimate 0 is consistent. But due to the large amount of parameters involved with FIR models, it is usually associated with a large variance error, compared to the other parametric estimation methods. Two major advantages using the correlation analysis method are: (a) a simple closed-form solution can be
1401
put/measured disturbance and output channels or some of it just represent the nOise/disturbance. For a SISO system, the problem can be examined via the following null hypothesis test: the output sequence y is uncorrelated with the input sequence u . Under this hypothesis, we have (Bartlett, 1955)
var[-yuy(l)] ~ O"~O";(N _l)-I
lated to their variance through the following expressions respectively:
(8)
where the input sequence u should be prewhitened using a prewhitening filter and 0"2 denotes the variance. Strictly speaking, relation (8) only applies to the case where u is Gaussian distributed. But it is a common practice to use the relation for input sequences of other probabilic distributions. Because usually, N » l , (8) is sometimes approximated as:
Substituting (12) and (13) into (11), the following bounds on the standard deviation of 8 are obtained:
(9) where 'Yi is the ith element of the vector 'Y which can represents 'YUjYk or 'YdjYk for some j and k. Under the hypothesis that a output k is uncorrelated with a input or a measured disturbance j, varbd is approximated by:
Hence, with 99% confidence level, 'Yuy should stay in the range set by the following expression: -2.58
(N _ l) ~ y~
'Yuy () l
~ 2.58
(N _ l) y~
(10)
(15) when all the inputs and output i are scaled by their variance. To use relation (15) , all the measurement sequences should be prefiltered by a prewhitening filter of sufficiently high order. Relations (14) and (15) can be used as a criterion to detect the existence of transfer functions . If the estimated Bi from (6) goes beyond the bounds set by (14), existence of a transfer function between the input j and the output i is suggested.
If 'Yuy goes out of the range in (10), a transfer function between the output y and the input u is suggested. For a MIMO system, if sequences {u(t)}, {d(t)} and {v(t)} are uncorrelated to each other, the above method can be directly applied to each input or measured disturbance/output channel. But assuming that {u(t)} in each channel is statistically uncorrelated to each other will be too restricting in practice. For example, some processfriendly input sequences such as many multi-frequency signals will be excluded in consideration. Furthermore, it is usually not so realistic to assume that {d( t)} and { d( t )} / {u (t)} are uncorrelated to each other. When the input sequences and/or the measured disturbance sequences are correlated, using criterion (10) will lead to incorrect conclusions. For example, assume that Ui and Uj are correlated, and there exists a transfer function between Yk and Ui but not between Yk and Uj . Then amplitude of 'YUjYk is determined by 'YU j Ui and does not necessarily stay in a band set by (10) .
4. APPLICATION TO IDENTIFICATION EXPERIMENT MONITORlNG The correlation analysis estimator (6) forms a simple linear equation set. A recursive algorithm can be readily built. Denote
cp(t)=[uI(t-1) . .. uI(t-k)U2(t-1) . .. uni(t-k) dl (t-1) . .. dl(t-k) .. . dnd(t-k)]T
Another application of the correlation analysis estimator (6) we consider in this section is to establish a transfer function detection criterion. The criterion is a MIMO system extension to the singular value decomposition method of SISO systems discussed by Rivera (Rivera, 1992) . From equation (6) , we have Q:2(fo)1I8ill~ ~ lbi l1 2 ~ 0'2(fo)IIBill~
(6) can then be computed recursively as
i = 1, ... , no
L(t) _ P(t - l)cp(t) -1/a+cpT(t)P(t-l)cp(t)'
(11)
where O'(fo) and Q:(fo) are the maximum and the minimum singular values of ro · 118ill~ and Ibill~ can be re-
P(t)
= P(t _
) 17
1) _ P(t - l)cp(t) cpT(t)p(t - 1) (18)
l/a
1402
(
+ cpT(t)P(t - l)cp(t)