Performance Assessment for a Class of Nonlinear Systems - Industrial

In this paper, a class of specified nonlinear systems are concerned. The nonlinear process model can be considered as a simplified Kolmogorov–Gabor ...
1 downloads 0 Views 2MB Size
ARTICLE pubs.acs.org/IECR

Performance Assessment for a Class of Nonlinear Systems Zhi Zhang, Li-Sheng Hu,* and Xue-Lian Zhang Department of Automation, Shanghai Jiao Tong University, and Key Laboratory of System Control and Information Processing, Ministry of Education of China, Shanghai, 200240, People’s Republic of China ABSTRACT: In this paper, a class of specified nonlinear systems are concerned. The nonlinear process model can be considered as a simplified KolmogorovGabor polynomial model, which is motivated by considering the actuator nonlinearity and the process output nonlinearity. In terms of this nonlinear system, the analytical performance benchmark and the associated benchmark controller are established by minimizing the variance of the measurement output. The existence conditions of feedback invariant are discussed. The calculation of performance index, based on correlation analysis, is proposed. To facilitate assessment, a filtering or whitening method is presented. The white noise sequence is estimated by dynamic principal component analysis (DPCA), based on the nonlinear autoregressive exogenous (NARX) model. The numerical examples are given to illustrate the effectiveness of the proposed assessment method.

1. INTRODUCTION In modern process industries, there are many control loops that are initially designed well to achieve the expected optimal performance. However, these designed controllers may gradually perform poorly, because of complex operating circumstances that are different from those envisaged at the design stage.1 Hence, it is necessary to offer an well-suited tool for control engineers to determine if specified performance targets and response characteristics are being met by the controlled process variables.2 Research on the control loop performance assessment has been increasing in the last two decades. By virtue of valuable contributions from Harris, time series analysis techniques could be used to find a suitable expression for the feedback controller invariant from routine operating data and the subsequent use of this as a benchmark to assess control loop performance.3 The associated minimum variance (MV) controller, as a benchmark controller, can also be obtained, which can be traced to the work of Åstr€om.4 The greatest advantage of this technique is that the unique performance benchmark can be analytically described by the parameters of the open-loop information, after model parameters are predetermined or estimated. A well-established performance benchmark can be used to test the deviation from the current control performance. Subsequently, Harris et al.,5 as well as Huang and Shah,6 extended the MV performance benchmark to the multivariable control loop performance assessment. Some excellent works on the control loop performance assessment can be found in the literature1,2,79 and references therein. However, it is well-known that most industrial processes are nonlinear in nature. The traditional methods based on linear systems may lead to biased results of performance assessment for nonlinear systems.10,11 The challenges of performance assessment for nonlinear systems are both in modeling and parameter estimation. These will bring difficulties for extending the traditional methods of linear systems to nonlinear ones. Recently, the pioneering contributions from Harris and Yu10 pointed out the direction of this issue. Harris and Yu10 gave many correct and interesting results. They proved that minimum variance feedback r 2011 American Chemical Society

invariant still exists for a class of nonlinear systems,12 which means that the minimum variance performance benchmark can be estimated from routine data. They also indicated that if the nonlinearities before the d-step (where d is the time delay in the control input channel) exist in the process model, the minimum variance feedback invariant does not exist. Based on these ideas, some excellent works of Yu et al.11,13,14 promoted the issue forward further when the control system is in the presence of nondifferentiable nonlinearities. For general nonlinear models, Yu et al.15 proposed a new performance index based on variance decomposition method. If the model information is available, Zhang et al.16 explored the feasibility of performance assessment for nonlinear control systems based on the T-S fuzzy model. Consider a class of specified nonlinear systems, this paper will explore establishing the analytical performance benchmark and the associated benchmark controller. The existence of feedback invariant for this nonlinear system is discussed. A whitening method using dynamic principal component analysis (DPCA), based on the nonlinear autoregressive exogenous (NARX) model, is also proposed. The major contribution of this paper is to provide the analytical expression of performance benchmark for a class of specified nonlinear systems and analyze the effects of nonlinearities before the d-step on the calculation of performance index in the higher-order statistics form. The remainder of this paper is organized as follows. Problem formulation on performance assessment for a class of nonlinear systems is presented in Section 2. Section 3 presents the performance benchmark and the associated benchmark controller for the nonlinear control system. The existence of the feedback invariant is discussed. In Section 4, the calculation of performance index based on correlation analysis is proposed. Section 5 presents a whitening method using DPCA based on NARX model for nonlinear measured data. In Section 6, the numerical Received: December 13, 2010 Accepted: August 6, 2011 Published: August 07, 2011 10557

dx.doi.org/10.1021/ie200884s | Ind. Eng. Chem. Res. 2011, 50, 10557–10566

Industrial & Engineering Chemistry Research

ARTICLE

examples show the effectiveness of the proposed method. Finally, the conclusion remarks and the considerations of future works are given in Section 7.

2. PROBLEM FORMULATION Consider the following nonlinear system: ~ ðq1 Þut þ qd G ~ q ðq1 Þut 2 þ Gd ðq1 Þy t 2 y t ¼ qd G þ Ga ðq1 Þat

ð1Þ

where at ∈ R and ut ∈ R are the disturbance input and the control input of the system, respectively. yt ∈ Rp denotes the output of the system, and the quadratic terms ut2 and yt2 2 2 2 are defined as follows: ut2 := (ut,1 ,ut,2 ,...,ut,m )T, and yt2 := 2 2 2 2 T ~ (q1), ,yt,2 ,...,yt,p ) , respectively. Ga(q1), G(q1) := qdG (yt,1 1 d ~ 1 1 b ^ 1 Gq (q ) := q Gq (q ), and Gd (q ) := q Gd (q ) are the transfer functions for the associated input/output pairs. It should ~ q(q1), and G ^ d (q1) are the time~ (q1), G be noted here that G delay free parts of G(q1), Gq(q1), and Gd(q1), respectively. The disturbance input at is the white noise with zero mean, and covariance Σ2a. The time delay in the control input channel is represented by d (d g 1), and b denotes the time delay of the quadratic term yt2 (b g 1). It is necessary to know b and d. Because most industrial plants are characterized by nonlinear behaviors, the representations used to describe nonlinear processes are widely studied. One of the most popular classes of nonlinear model are block-oriented nonlinear models. Three of the more-common model structures are the Hammerstein model, the Wiener model, and the feedback block-oriented model.17 The Hammerstein model is a nonlinear model with a static nonlinear function followed by a linear dynamic block. The Wiener model is a nonlinear model with a linear dynamic block, followed by a static nonlinear function. The feedback blockoriented model consists of a static nonlinearity in the feedback path around a linear time-variant system. From the polynomial point of view, a nonlinear model without output feedback can be represented by the Volterra series model, and a nonlinear model with output feedback (e.g., NARMAX) can be represented by the KolmogorovGabor polynomial model.18 For example, a single-input single-output (SISO) nonlinear process can be described as a general KolmogorovGabor polynomial model with order k,18 m1

m2

in the inputs but nonlinear in the outputs. The Hammerstein model and Wiener model are described as mentioned in the preceding paragraph. It is interesting to note that the Hammerstein model, the Wiener model, and feedback block-oriented model have their corresponding NARMAX representations.17 In practical applications, the polynomial model is an useful modeling tool for the nonlinear systems, which is a good extension of the linear systems. The polynomial models have been receiving widespread attention. The process identifications of the distillation column and the heat-exchanger process using the Hammerstein model were studied in ref 20. The multivariable Hammerstein model and the multivariable Wiener model were used in ref 21 to describe the nonlinear behavior of binary system. Considering the phenomenon of output multiplicity which appears in some distillation processes, ref 22 used the feedback block-oriented model to represent the process, and that cannot be modeled by the Hammerstein or Wiener models. Ref 23 also reported that a form of the NARMAX model is used to model a real five-plate distillation column. The description of nonlinear process in this paper can be considered as a class of simplified KolmogorovGabor polynomial models. The stated model can be motivated by considering the actuator nonlinearity and the process output nonlinearity. Remark 1: In this paper, we only consider that the transfer function G(q1) from ut to yt has the decomposition, G(q1) = ~ (q1). Actually, it is a very simple case. For the more general qdG case (e.g., diagonal interactor and general interactor cases, following the procedure proposed in ref 6), the procedure developed in this paper is applicable to the case that the transfer function G(q1) has a more general time-delay property. The control performance assessment problem considered in this paper is to design a controller to minimize the variance Var{yt} of the measurement output yt. Furthermore, using the minimum variance Var{yt}|mv to assess control systems. In this paper, the controller to be designed is described as ut ¼  C1 ðq1 Þy t  C2 ðq1 Þy t 2  C3 ðq1 Þut 2 1

1

ð2Þ

1

where C1(q ), C2(q ), and C3(q ) are the controller transfer functions. The optimal problem then can be stated as Varfy t gjmv : ¼ min

C1 , C2 , C3

Varfy t g

subject to eq 1. To simplify the sequel development, q1 will be omitted, unless circumstances necessitate its presence.

where y is a mean value, and uti and ytj denote the input and output data at time instant t  i and t  j, respectively. According to the different motivations, the Kolmogorov Gabor polynomial model can be simplified as a parametric Volterra series model, a nonlinear differential equation model, a Hammerstein model, a Wiener model, and so on.19 The parametric Volterra series model realizes a linear feedback and models a nonlinearity only for the inputs. For example, PARX model in ref 10 can be considered as a class of parametric Volterra series models. The nonlinear differential equation model is the counterpart of the parametric Volterra series model, which is linear

3. BENCHMARKING TO NONLINEAR CONTROL SYSTEMS Before developing an algorithm for performance assessment of the control system, the performance benchmark and associated benchmark controller must be set up first. In this paper, we take minimum variance of the output yt as a benchmark for the performance assessment. In this section, a minimum variance controller will be designed and the associated minimum variance performance will be found for the nonlinear system that is described by eq 1. The closed-loop system with the controller described by eq 2 can be written as ~ ð  C1 y t  C2 y t 2  C3 ut 2 Þ þ Ga at y t ¼ qd G ~ q ut 2 þ G d y t 2 þ qd G 10558

dx.doi.org/10.1021/ie200884s |Ind. Eng. Chem. Res. 2011, 50, 10557–10566

Industrial & Engineering Chemistry Research

ARTICLE

Thus, eq 10 can be rewritten as

and

~ ut þ G ~ q ut 2 þ q d  b G ~ a at ^ d y t 2 þ G̅ a atþd þ G y tþd ¼ G

~ C1 Þ1 ½Ga at þ ðGd  qd G ~ C2 Þy t 2 y t ¼ ðI þ qd G

ð12Þ

~q  G ~ C3 Þut 2  þ qd ðG Decompose Ga and Gd as ~a Ga ¼ G̅ a þ qd G

ð3Þ

~d Gd ¼ G̅ d þ qd G

ð4Þ

We then consider the existence of the feedback invariant for two cases: ~ d, where Gd = 0 and Case 1: If b g d, we have Gd = Gd + qdG db ^ ~ Gd = q Gd. As proved in refs 10 and 11, the conditional expectation of yt+d can be obtained by ~ ut þ G ~ q ut 2 þ q d  b G ~ a at ^ d y 2t þ G ^y tþdjt ¼ Efy tþd jIt g ¼ G

~d = respectively. Note that if b g d, we have Gd = 0 and G ^ d. Furthermore, by applying the matrix inverse lemma,24 it qd-bG yields

ð13Þ where It is the information set.10 From eqs 1113, we have the prediction error

~ C1 ½ðG̅ a þ qd G ~ C 1 ÞG ~ a Þat y t ¼ ½I  qd ðI þ qd G ~q  G ~ C3 Þut 2  ~ d  qd G ~ C2 Þy t 2 þ ðG þ ðG̅ d þ qd G

y tþd  ^y tþdjt ¼ ð1 þ Ga, 1 q1 þ 3 3 3 þ Ga, d1 qd þ 1 Þatþd ¼ ̅Ga atþd

which directly leads to

Hence, under minimum variance control, the ^y t+d|t equals zero, the process output yt+d|mv will be dependent only on the recent d past disturbances:

y t ¼ G̅ a at þ G̅ d y t 2 ~a  G ~ C1 G̅ a Þat ~ C1 Þ1 ½ðG þ qd ðI þ qd G ~d  G ~ C2  G ~ C1 G̅ d Þy t 2 þ ðG ~q  G ~ C3 Þut 2  þ ðG

y tþd jmv ¼ ð1 þ Ga, 1 q1 þ 3 3 3 þ Ga, d1 qd þ 1 Þatþd ¼ G̅ a atþd

ð5Þ

If we want to get the feedback control invariant from eq 5, we may choose the controller transfer functions C1 as ~ 1 G ~ a G̅ 1 C1 ¼ G a

ð6Þ

C2 as ~ 1 G ~ 1 G ~d  G ~ a G̅ 1 C2 ¼ G ̅ d a G

In this case, eq 9 can be written as y t jmv ¼ G̅ a at

~ d. Since t + d  b > t, Case 2: If b < d, we have Gd = Gd + qdG according to Case 1, we have y tþdb ¼ ^y tþdbjt

ð7Þ

þ ð1 þ Ga, 1 q1 þ 3 3 3

C3 as ~ 1 G ~q C3 ¼ G

þ Ga, db1 qd þ b þ 1 Þatþdb

ð8Þ

ð9Þ

¼

y tþd

~ ut þ G ~ q ut 2 þ ̅Ga atþd þ G ~ a at G

h   i2 ^ d ^y tþdbjt þ 1 þ Ga, 1 q1 þ 3 3 3 þ Ga, db1 qd þ b þ 1 atþdb þG

However, the existence of the possibility that eq 9 can be called as the feedback invariant should be discussed. The discussions about the existence of the feedback invariant are stated as follows. Equation 1 can be rewritten as 2 2 ~ utd þ G ~ q utd ^ d y tb þ G þ Ga at yt ¼ G

2 ~ ut þ G ~ q ut 2 þ ̅Ga atþd þ G ~ a at þ G ^ d ð^y ¼G tþdbjt Þ

 i2 1 þ Ga, 1 q1 þ 3 3 3 þ Ga, db1 qd þ b þ 1 atþdb h   i ^ d ^y tþdbjt 1 þ Ga, 1 q1 þ 3 3 3 þ Ga, db1 qd þ b þ 1 atþdb þ 2G ^d þG

ð10Þ

We have the Diophantine identity as follows:

h

2 ~ q ut 2 þ ̅Ga atþd þ G ~ a at þ G ^ d ð^y ~ ut þ G ¼G tþdbjt Þ 2 ^ d atþdb þG þ

where Ga,0 (Ga,0 = 1) and Ga,i (for i = 1, ..., d  1) are Markov ~ a is the parameters of the disturbance transfer function, and G remaining rational and proper transfer function matrix. Let Dt = Gaat; we then have ~ a at Dtþd ¼ G̅ a atþd þ G

ð15Þ

and eq 12 can be rewritten as

respectively, we have y t jmv ¼ G̅ a at þ G̅ d y t 2

ð14Þ

þ2

10559

∑ i¼1

2 Ga,2 i atþdbi þ 2

db1



i¼1



j¼i þ 1



i¼1

Ga, i atþdb atþdbi



db2 db1

^ d ^y tþdbjt þ 2G

ð11Þ

db1

Ga, i Ga, j atþdbi atþdbj

atþdb þ

!

db1



i¼1

Ga, i atþdbi

ð16Þ

dx.doi.org/10.1021/ie200884s |Ind. Eng. Chem. Res. 2011, 50, 10557–10566

Industrial & Engineering Chemistry Research

ARTICLE

The conditional expectation of yt+d can be obtained by

Remark 2: To illustrate, consider a SISO system ! 1 2 at yt ¼ gut2 þ hyt1 þ 1  cq1

~ ut þ G ~ q ut 2 þ G ~ a at ^y tþdjt ¼ Efy tþd jIt g ¼ G db1



^ d ð^y tþdbjt Þ2 þ G ^ d Σa 2 ð1 þ þ G þ 2

db1



i¼1

Ga, i þ 2

i¼1

db2 db1





i¼1

j¼i þ 1

Ga,2 i

Ga, i Ga, j Þ ð17Þ

where at ∼ N(0,σa2), d = 2 and b = 1. We have ! 1 ytþ1 ¼ gut1 þ hyt 2 þ atþ1 1  cq1

where It is the information set.10 From eq 16 and eq 17 we have the prediction error 0 y tþd  ^y tþdjt ¼ G̅ a atþd þ þ

db1



i¼1

þ2

2 Ga,2 i atþdbi þ 2





j¼i þ 1

!



∑ i¼1



j¼i þ 1

Ga,2 i þ 2

db1



Ga, i

i¼1

!

ytþ2 ¼ gut þ hðgut1 þ hyt 2 þ cDt þ atþ1 Þ2 þ Dtþ2

¼ gut þ hðgut1 þ hyt 2 þ cDt þ atþ1 Þ2

Ga, i Ga, j

∑ i¼1

þ atþ2 þ catþ1 þ c2 Dt

!

db1

^ d ^y tþdbjt atþdb þ þ 2G

Ga, i atþdbi

ð18Þ

y tþd  ^y tþdjt ¼ G ̅ a atþd db1

∑ i¼1

þ2

^ytþ2jt ¼ Efytþ2 jIt g ¼ gut þ hðgut1 þ hyt 2 þ cDt Þ2



i¼1



j¼i þ 1

∑ i¼1

Ga, i atþdb atþdbi

!

¼ gut þ hðgut1 þ hyt 2 þ cDt Þ2 þ hσ a 2 þ c2 Dt

Ga, i Ga, j atþdbi atþdbj

db1



i¼1

db2 db1

∑ i¼1

þ hEfatþ1 2 jIt g þ 2hðgut1 þ hyt 2 þ cDt ÞEfatþ1 jIt g þ Efatþ2 jIt g þ cEfatþ1 jIt g þ c2 Dt

db1

db2 db1

^ d Σa 2 1 þ G þ2

The conditional expectation of yt+2 can be obtained by

2 ^ d @atþdb þ G

2 Ga,2 i atþdbi þ 2

¼ gut þ hðgut1 þ hyt 2 þ cDt Þ2 þ hatþ1 2 þ 2hðgut1 þ hyt 2 þ cDt Þatþ1 þ atþ2 þ catþ1 þ c2 Dt

Substituting eq 15 into eq 18 gives 0

þ

Ga,2 i þ 2

∑ Ga, iGa, j j¼i þ 1

!

db1



i¼1

Ga, i

  ^ d y tþdb  1 þ Ga, 1 q1 þ 2G

where It is the information set.10 Thus, under MV control, the ^yt+2|t equals zero. The MV controller can be obtained: ! ! c2 h h ut ¼  σ2 Dt  ðgut1 þ hy2t þ cDt Þ2  g g a g

þ 3 3 3 þ Ga, db1 qd þ b þ 1 Þatþdb ðatþdb  þ

db1



i¼1

Ga, i atþdbi

! 1 atþ2 1  cq1 ! #2 1 þ atþ1 1  cq1

Let Dt = [1/(1  cq1)]at; then,

db2 db1 i¼1

! 1 þ atþ2 1  cq1

Ga, i atþdb atþdbi

Ga, i Ga, j atþdbi atþdbj

db1

^ d Σa 2 1 þ G þ2



gut þ hytþ1 þ

¼ gut þ h gut1 þ hyt 2

db1 i¼1

¼

"

db2 db1 i¼1

ytþ2

2 ^ d @atþdb G

2

We have the prediction error ð19Þ

This equation indicates that, even under MV control, a feedback invariant does not exist, because of the presence of yt+db. Hence, if Case 1 (b g d) holds, the output yt and the quadratic outputs yt2 are dependent only on the disturbance before the d-step. Hence, the variance Var(yt|mv) of the output yt is invariant to the controller designing, which means that Var(yt) reaches its minimum. We then call the controller (described by eq 2) as a minimum variance controller for the nonlinear system (described by eq 1).

ytþ2  ^ytþ2jt ¼ hðatþ1 2  σa 2 Þ þ 2hðytþ1  atþ1 Þatþ1 þ atþ2 þ catþ1 This equation illustrates that, if b < d, even under MV control, a feedback invariant does not exist, because of the presence of yt+1. Remark 3: If b g d, eq 2 demonstrates the controller designing method for the nonlinear system. It may be a nonlinear controller, which can guarantee MV output of the nonlinear system. If b g d, the MV controller described by eq 2 is theoretical and may not be applicable in practice, because of its 10560

dx.doi.org/10.1021/ie200884s |Ind. Eng. Chem. Res. 2011, 50, 10557–10566

Industrial & Engineering Chemistry Research

ARTICLE

poor robustness and excessive control action. Nevertheless, as a benchmark, it does provide useful information. Remark 4: The process described by eq 1 is a direct extension of the linear systems. Without the nonlinearity in the process described by eq 1, the controller described by eq 2 will be reduced to the minimum variance controller specified by eq 6, which is also discussed in ref 6. With the nonlinearity in the process described by eq 1, the variation of the output yt is affected by the nonlinearity. If b < d, eq 9 shows that the nonlinearity will affect the minimum variance of the output. Unlike the case of linear systems, the yt|mv term in eq 9 is no longer Gaussian,25 and E{yt|mv} 6¼ 0, even if the disturbance input at is white noise. Furthermore, if b < d, the cross-correlation E{ytaTt } between the output and the disturbance input will be dependent not only on Ga, but also on Gd. All of these may indicate that, if b < d, the result of the performance assessment is false, since the feedback invariant does not exist. To facilitate the following development, Ga and Gd are rewritten as the Markov parameter representations

and we have d1



Efy t y Tt gjmv ¼ þ þ

d1

d1

∑ Gd, i j∑¼ 0 my y ði, jÞGTd, j i¼0 2 2

d1



i¼0

þ

Ga, j Σa 2 GTa, j

j¼0

d1



Gd, i

d1

j¼0

ð21Þ

my2 a ði, jÞGTa, j

d1

∑ Ga, i j∑¼ 0 may ði, jÞGTd, j i¼0 2

T where my2a(i,j) = may 2 (i,j) and my2y2(i,j) are higher-order moments and are defined as

2 aTt  j g my2 a ði, jÞ ¼ Efy ti

    T 2 may2 ði, jÞ ¼ E ati y tj

G̅ a ¼ Ga, 0 þ q1 Ga, 1 þ 3 3 3 þ qd þ 1 Ga, d1 and

and 1

G̅ d ¼ Gd, 0 þ q Gd, 1 þ 3 3 3 þ q

d þ 1



Gd, d1 m

respectively. Consider two cases as discussed above: ~ d = qdbG ~ d, where Gd = 0 and G ^ d. (1) If b g d, Gd = Gd + qdG Equation 9 can be rewritten as y t jmv ¼ G̅ a at ¼ Ga, 0 at þ q1 Ga, 1 at þ 3 3 3 þ qd þ 1 Ga, d1 at ¼

d1



j¼0

qj Ga, j at

y2 y2

ði, jÞ ¼ E

2 y ti



2 y tj

T 

Note that ( 3 )|mv corresponds to Case 2.

4. CALCULATION OF PERFORMANCE INDEX Consider the definition of the performance index for the control systems: 8 9 < n o n o1 = ð22Þ ηðdÞ ¼ trace E y t y Tt E y t y Tt : ; mv

and Efy t y Tt gjmv

¼

d1



j¼0

Ga, j Σa GTa, j

ð20Þ

2

Note that E{ 3 } denotes the expectation operator. As discussed in Case 1, the feedback invariant exists, and E{ytyTt }|mv is the minimum variance of the output with the minimum variance controller in the place. However, for real applications, the actual controller used may not be the minimum variance controller, which means the variance E{ytyTt } of the real process output will be different from E{ytyTt }|mv. Assessing this difference will provide a way to monitor the quality of the control systems. ~ d. Equation 9 can be (2) If b < d, we have Gd = Gd + qdG rewritten as y t jmv

¼

1

d þ 1

Ga, 0 at þ q Ga, 1 at þ 3 3 3 þ q

Clearly, when the feedback invariant exists and a process is under the minimum variance control, we have E{ytyTt }|mv = E{ytyTt } (that is, η(d) = 1; otherwise, E{ytyTt }|mv < E{ytyTt }, and 0 < η(d) < 1). Considering the cross-correlation between yt and at: (1) If b g d, the feedback invariant exists. We have yt|mv = Gaat, n o d1 T Ga, i Σa 2 i 0 and M estimate the matrix Gd,i, for i = 0, 1, ..., d  1, if E{ytaTti} 6¼ 0, using eq 27: Gd :¼ ½Gd, 0 , Gd, 1 ; :::; Gd, d1  ¼ Γd Md1 where Γi = [E{ytaTti},E{ytaTti1},...,E{ytaTtid+1}], for i = 0, 1, ..., ∞. Similarly, from eq 26, we have Ga :¼ ½Ga, 0 , Ga, 1 ; :::; Ga, d1  ¼ Γ0 Σ1  Gd M0 Σ1 ¼ Γ0 Σ1  Γd Md 1 M0 Σ1

where Fya(i) = E{ytyTt }1/2ΓiΣ1/2 and Fay(i) = Σ1/2ΓTi E{ytyTt }1/2 are the cross-correlation matrix between yt and at. Let Σy = diag{E{yyT}E{yyT}, ..., E{yyT}E{yyT}}. γya(i) = Σy1/2Mi Σ1/2 and γay(i) = Σ1/2MTi Σy1/2 are the cross-bicoherence matrices of ~ 0Σy1/2. the signals between yt and at. γyy(0) = Σy1/2M Note that γyy is not an exact auto-tricoherence matrix. The trace of every component on the primary diagonal is