Building an EDM Process Model by an Instrumental Variable Approach

Building an EDM Process Model by an Instrumental Variable Approach

Available online at www.sciencedirect.com Procedia CIRP 6 (2013) 456 – 462 The Seventeenth CIRP Conference on Electro Physical and Chemical Machinin...

498KB Sizes 1 Downloads 148 Views

Available online at www.sciencedirect.com

Procedia CIRP 6 (2013) 456 – 462

The Seventeenth CIRP Conference on Electro Physical and Chemical Machining (ISEM)

Building an EDM process model by an instrumental variable approach Ming Zhou*, Zhigang Chen, Jufen Niu Beijing Engineering Research Center of monitering for construction safety,school of Machtronic and Automobile Engineering, Beijing University of Civil Engineering and Architecture, Beijing, China, 100044 * Corresponding author. Tel.: 86-010-68322103 fax: 86-010-68322514.E-mail address: [email protected].

Abstract An electrical discharge machining (EDM) process model is critical in control system for stabilizing machining process and improving productivity. A two order AR model was firstly proposed for the process, but the deduced formulation of the model parameter estimates is proved biased. An instrumental variable approach is then proposed to correct the biased estimates. This approach was completed by two interactive Kalman filters. One filter provides instrumental variables, namely gap state estimates, through the latest estimated parameters; while the other uses the gap state estimates to recursively estimate a two order autoregressive (AR) model parameters. The two Kalman filters by way of interactively supporting each other forms a new feasible way of online modeling an EDM process. © 2013 The Authors. Published by Elsevier B.V. Selection and/or peer-review under responsibility of Professor Bert Lauwers

© 2013 The Authors. Published by Elsevier B.V. Selection and/or peer-review under responsibility of Professor Bert Lauwers Keywords: Electrical discharge machining; Modeling; Instrumental Variable; Kalman filter

1. Introduction Electrical discharge machining (EDM) method brings two advantages over traditional machining methods; one is its capability of machining conductive materials with no or little consideration of material hardness, the other is its untouched machining which generates little machining forces to the machined surfaces in a process. These two advantages extend its applications broadly in manufacturing world, especially in punches, dies, and micro machining, etc. However, EDM productivity is quite low because even a state-of-art EDM tool has only open-loop or half closed-loop control systems. The main problem to form a complete closed-loop control system is the difficulty of modelling the process precisely, because the metal removal process of EDM involves non-stationary, nonlinear and stochastic characteristics, which outdated classical linear control strategies in EDM control system. Adaptive control strategies from modern control theories seemed finding their applications in the process, since an adaptive controller can modify its behaviours in response to a rapidly varied

2212-8271© 2013 The Authors. Published by Elsevier B.V. Selection and/or peer-review under responsibility of Professor Bert Lauwers doi:10.1016/j.procir.2013.03.087

machining state in the process [1]. However, the success of a model-predictive-control performance largely depends on an accuracy of a process model [2]. Therefore, it is prerequisite to build a process model precise enough to reflect the on-going machining states in an EDM process. Initially, a feasible process model might be developed from a comprehensive study of its dynamical properties. It had been proved that the EDM process is a deterministic nonlinear process of chaotic behaviour [3]. To further understand the complexity of an EDM process, a chaotic analysis technique, called the method of correlation dimension, was employed to determine the fractal dimensions of the process, a quantification of its complexity. The computation turns out that the complexity of the process is slightly less than 2 correlation dimensions. This result makes clear that an EDM process is of a low dimensional structure [3]. In other words, it is enough to use a two-order system to approximate the process that has also been accidentally evidenced by a formerly built process model [4]. To further understand the process behavior, this paper builds a process model by an instrumental variable

457

Ming Zhou et al. / Procedia CIRP 6 (2013) 456 – 462

approach from another perspective, completed by two Kalman filters. One filter serves the functionality of estimating gap states from gap state identifications with the latest estimated process parameters; the other serves to recursively estimate the process (an AR model) parameters with the gap state estimates and then provides one-step gap state predictions.

machining state is identified to deteriorate, it is better to alter the current machining situation immediately than gradual recovery from the deteriorated states; otherwise, the workpiece surface may be damaged by arcing. In this sense, to achieve such an effect by control action, it is preferable to let m 1 in B(q) polynomial.

2. Formulating the process model

s(k ) -an estimate of a gap state s(k ) at time k

Gap states are defined to be working indexes used to monitor EDM processes [4]. In this paper, we employed a way of sampling gap voltages and currents to identify a gap state with an A/D of 400k/s sampling rate, each 200k/s. It is well known that there are five sorts of discharging pulses in gap in machining, spark pulse, transient arc pulse, stable arc pulse, short pulse and open pulse. The time between two discharging pulses is called pulse-off. Each time as a gap voltage and a gap current are sampled, the state of this sampling from a discharging pulse or a pulse-off is immediately identified by pulse discrimination criterion [5]. Consequently, there are six sampling states, namely spark state, transient arc state, stable arc state, short state, open state, and pulse-off state. To obtain a gap state identification, calculate the ratio in Eq. (1) below.

y(k ) -a gap state identification at time k

y

(

spark

(1) where,

tran. arc

spark

,

) /(

tran.arc

spark

,

tran. arc

stab.arc

,

short

stab. arc

,

open

short

open

)

are respectively

accumulated numbers of spark state, transient arc state, stable arc state, short state, and open state in a time interval, y is a gap state identification. Gap state identification y is quantified machining state that can be used to monitor machining processes. Because an EDM process is of nonlinear characteristics and it has been proved that it is sufficient to use a two order system to approximate the process [3], we suggested a timely varied linear parameter model approximating the nonlinear behaviours of the process, namely tracking varied gaps states. Let a timely varied auto-regression with extra variable (ARX) model approximate an EDM process, defined by (2) A(q)s(k ) B(q)U (k ) w(k ) where A(q) 1 l a q i , i

i 1

B(q ) [

m j 1

b1 j q j ,

m j 1

b2 j q j ,

,

m

bnj q j ],

q

1

is a backward

j 1

shift operator, U T (k ) [u1 (k ), u2 (k ), , un (k )] is an input control variable vector containing power system and servo system variables u1 , u2 , , un , s(k ) is a gap

s(k ) -a gap state at time k

s (k ) -one-step gap state prediction at time 2 u

-the process noise variance

2 v

-the identification noise variance

2 1

,

2

-the process parameter noise variances

An EDM process signal in this paper, a series of gap state identifications, is not like the other signals obtained by manoeuvring the inputs and measuring the outputs. The state-of-art EDM tools are almost open-loop controlled. Open-loop controlled means that all EDM power parameters and servo system parameters once set are kept unchanged in machining and none of them are automatically adjusted to the changing gap states. From a mathematical view, this unchanged control variable which is in fact a constant, will not contribute to gap state variations. In such a case, we may neglect the influences of the control variables on the outputs of the model we will build; consequently, there will be no control variable in the model. Then a second order ARX model turns into an AR model, 2

s (k )

ai s(k i) w(k )

s(k 1) s(k 2)

i 1

s(k 1) s(k 2)

a1 a2

w(k ),

w(k )

(3)

T

where

a1 a2 is a vector containing two timely varied parameters a1 and a2 . Theoretically, a gap state identification y(k) may be assumed to be a corrupted gap state by an additive white noise n(k) with variance v2 , defined by y (k )

s(k ) n(k )

(4) where y(k ) , an gap state identification, is also a noisy observation of a gap state s(k ) .Substituting Eq. (3) in Eq. (4), y (k )

2

s ( k ) n( k )

ai s(k i) w(k ) n(k )

i 1

2

state at time k , w(k ) is assumed a zero-mean white noise disturbance with variance u2 . Because a two order system can approximate an EDM process, we took l 2 in A(q) polynomial. Practically, in a control scheme, if

2

k

ai y(k i)

i 1

y(k 1) with

2

ai n(k i) n(k ) w(k )

i 1

y (k

2)

(k )

(5)

458

Ming Zhou et al. / Procedia CIRP 6 (2013) 456 – 462

(k )

2

w(k ) n(k )

ai n(k i) .

enhancing a gap state identification to obtain a gap state estimate and the other providing estimated AR parameters.

(6)

i 1

Let Y p (k )

y (k )

T p 1) ,

y (k

(7)

where the subscript p denotes the number of vector rows [6]. Substituting the preceding Eq.(7) in Eq.(5) yields (8) y (k ) Y2T (k 1) (k ) . Extending the subscript p to N in Eq.(7) with consideration of Eq.(8), derive Y2T (k 1)

y (k )

(k )

YN (k ) y (k

Y2T (k

N 1) N

N)

(k

N 1)

(k ) .

(k )

(9)

Then, after simple manipulations, derive the least square (LS) estimate LS , LS

[

T N

(k )

N

(k )]

1

(10) with the estimate error [ NT (k ) LS

T N

(k )YN (k ) ,

N

(k )]

1

T N

(k ) (k ) .

3. Instrumental variable (IV) realization by two Kalman filters 3.1. Gap state estimates from gap state identifications To enhance a gap state identification y(k) to its gap state s(k), a Kalman filter algorithm is used to process the data to obtain the gap state estimate s(k ) , an enhancement of a gap state identification y(k). Because a gap state identification y(k) is a scalar, the state space representation of an EDM process is reduced to scalar expressions followed from Eqs. (3) and (4), s(k ) s (k ) K (k )( y (k ) s (k )) y (k ) s (k ) v(k ) (15) with i 2

s (k )

(11)

ai (k 1) s(k i)

i 1

From Eq.(5) and Eq.(6), it is found that (k ) is not a

(16)

white noise and y (k ) is correlated with

(k ) , so

where s(k ) is a gap state estimate at time k from a gap

(k ) (k ) will not tend to zero [6]. Consequently, the

state identification y(k), s (k ) is an one-step gap state

above estimate method in Eq.(10) will provide biased

prediction, v(k ) is an innovation sequence which is zero

parameter estimates.

mean white noise when the filter is optimal [6, 11],

T N

To obtain consistent estimates of AR parameters from noisy gap state identifications, the instrumental variable (IV) method has been used to overcome this biased

and K (k ) is called Kalman gain, given by K (k )

P (k )( P (k )

2

)

1

(17) where

2 v

is the variance of zero mean identification

estimate [6, 7]. The essential part of IV method is to

noise, and P (k ) is called the a priori error variance

replace

which is recursively estimated by

with a matrix

N (k )

with (k ) , that is, E[ TN (k ) (l )] 0 (12) Additionally,

T N

(k )

for

N

any

N (k )

k

uncorrelated and

l

.

(k ) is required to be non-singular.

The IV estimate of the AR parameters is hence given by , [Z NT (k ) N (k )] 1 Z NT (k )YN (k ) IV (13) and the estimate error is defined to be [Z NT (k ) N (k )] 1 Z NT (k ) (k ) . (14) IV IV From the above Eq.(12) and Eq.(14), it is deduced that IV estimate is unbiased and consistent. In this paper, the IV matrix is built by a filtered version of a gap state identification through a state estimator and a parameter estimator which allows the states and parameters to be estimated separately [8 10]. The two joint Kalman filters work in parallel, one

P (k )

a12 (k 1) P(k 1) a22 (k 1) P(k

2)

2 u

(18) with P( k )

(1 K (k )) P (k )

(19) where

2 u

is the variance of zero mean process noise,

and P(k ) is called the a posteriori error variance [6]. However, the AR parameters a1 (k ) , a2 (k ) and variance

2 u

,

2 v

are not known before hand and have to

be estimated or determined before running the Kalman filter for a gap state estimate. 3.2. Estimate of process parameters from a gap state estimate Substituting Eq.(16) in Eq.(15) leads to a formulation

459

Ming Zhou et al. / Procedia CIRP 6 (2013) 456 – 462

that a gap state estimate s(k ) can be expressed as a function of AR process model parameters, s (k )

2

ai (k 1) s (k i ) K (k )v(k )

i 1

[ s (k 1) s (k H (k ) (k 1)

where H (k )

2)]

a1 (k 1) a 2 (k 1)

K ( k )v ( k )

K (k )v(k ) ,

[s(k 1) s(k

(20) 2)] ,

(k 1) [a1 (k 1) a2 (k 1)]T . Because EDM process

is a non-stationary process, the process parameters could not be taken as invariants. However, the parameters are still assumed to be invariant in the time interval between two gap state identifications. Then, the process AR parameters are therefore defined by (21) (k ) (k 1) w (k ) where w (k ) is assumed zero-mean white noise with the 2 0 . covariance 1 Q

0

2 2

Then another Kalman filter can be used to estimate the process parameters. Followed from Eq.(20) and Eq.(21), the state space representation of parameters is (k ) (k 1) w (k ) . (22) y (k )

s (k )

H (k ) (k )

K ( k )v ( k )

The parameter measurement error covariance is 2 v

K 2 ( k )v 2 ( k )

K 2 (k )( y (k )

2

ai (k 1) s(k i)) 2

i 1

(23) Consequently, the estimate of parameters by Kalman filter algorithm is followed from preceding equations, (k ) (k 1) K (k )( y (k ) H (k ) (k 1)) y (k ) s(k ) H (k ) (k 1) v (k ) (24) where v (k ) is the innovation sequence, and K (k ) is a vector called Kalman filter gain, defined by the following equation [6, 11], K (k )

P (k ) H T (k )[ H (k ) P (k ) H T

2 v

]

1

(25) where P (k ) is called the a priori parameter error covariance matrix, which is recursively estimated by P (k )

P (k 1) Q

(26) where P (k ) is called the a posteriori parameter error covariance [6], given by P (k ) I 2 K (k ) H (k ) P (k ) (27) From the above Eq.s (21)-(27), the process parameters are estimated and fed to the Kalman filter to estimate a gap state for next cycle calculations.

Basically, the performance of a model is often evaluated by its adequacy in its applicable situations. The adequacy of an EDM model in this sense is its ability to track gap states with permitted errors. Thus, a criterion to testify the quality of an EDM model needs to be established. It should incorporate two kinds of errors: one is the permitted errors between gap state identifications and the corresponding gap states; the other is the permitted errors between one-step gap state predictions and the corresponding gap states. However, gap states are not known and only exist theoretically. A gap state estimate s(k ) is an enhancement from a gap state identification y(k) to the corresponding gap state. To be executable, we use gap state estimate s(k ) to substitute the gap state s(k). Then the model evaluation criterion is states: (1) The mean estimate errors between gap state estimates and the gap state identifications must be small within theoretical error bounds, calculated by [12] 2 ^ 1 N (28) J k1 y (k ) s(k ) , Nk1 . (29) Jk R12 N 2 ( 1 ) y (k ) N k 1

(2) The predictive errors between one-step gap state predictions and the gap state estimates must be small within theoretical error bounds, calculated by [12]

1 N

Jk2

2

N

s (k ) s (k ) ,

(30)

.

(31)

k 1

Jk

R22

N

( 1 ) s (k ) N k1

2

where s (k ) is an one-step gap state prediction. We employ a method, called grey correlation analysis, to obtain the optimized variances u2 , v2 , 21 , and

2 2

subject to the restrictions of the smallest of

estimate error and predictive error in a sense. Generally, strict restriction on prediction reduces prediction error while enlarging the estimate error and vice versa. To evaluate the model, we derived the optimized variances u2 =0.0325, v2 =0.30, 21 =0.01, and 22 =0.01 by grey correlation analysis with the restrictions of 70% weight on prediction and 30% on estimate. 4. Model validation Model validation was performed on a train of gap state identifications in Fig. 1. Fig. 1 shows the superimposed gap state estimates on the gap state identifications. Because the process is of non-stationary property and gap state in some parts varies abruptly, we specifically selected four representative parts in a short time interval (7.5 s) from the four different machining

Ming Zhou et al. / Procedia CIRP 6 (2013) 456 – 462

(c) Comparison of gap state estimates with gap state identificatons

1

Gap state identifications Gap state estimates

Gap state

0.8 0.6 0.4 0.2 0 11250

11251.5

11253

Time (Sec.)

11254.5

11256

11257.5

11256

11257.5

Errors between gap state estimates and gap state identifications

1

Theoretically positive bound

Errors

Error

0.5 0 -0.5 -1 11250

Theoretically negative bound 11251.5

11253

Time (Sec.)

11254.5

(d) Comparison of gap state estimates and gap state ientifications

1

Gap state identifications Gap state estimates

0.8

Gap state

situations to demonstrate the rational and stable estimates. It is easy to observe in Fig. 1(a) that gap state varies greatly. However, the developed model can still track the varied gap state identifications and the estimate errors fall in the theoretical bounds, indicating that the model works well in this phase. Fig. 1(b) and (c) shows parts of the estimates closely track the gap state identifications with very small errors in the restricted theoretical bounds. In deleterious machining process in Fig. 1(d), the gap state identifications appear erratic, but the estimates still track the identifications well with a few larger errors but still in the restricted theoretical bounds. These facts denote that the gap state estimates are stable and robust. Fig. 2 shows the one-step predictions based on the recovered or enhanced gap states, gap state estimates. In order to obtain clear impressions of precise one-step gap state predictions, we randomly select four parts from the

0.6 0.4 0.2 0 21000

21001.5

21003

Time (Sec.)

21004.5

21006

21007.5

21006

21007.5

Errors between gap state estimates and gap state identifications

1

Theoretically positive bound

Errors

0.5

Error

460

0 -0.5 Theoretically negative boound -1 21000

1.2

Gap state identifications Gap state estimates

1

a

b

c d

Gap State

0.8

0.6

21001.5

21003

Time (Sec.)

21004.5

Fig.1. Comparison of gap state estimates with gap state identifications in a non stationary process. (a) Comparison of gap state estimates with gap state identifications in phase A in 7.5 s. (b) Comparison of gap state estimates with gap state identifications in phase B in 7.5 s. (c) Comparison of gap state estimates with gap state identifications in phase C in 7.5 s. (d) Comparison of gap state estimates with gap state identifications in the deleterious machining process in 7.5 s.

0.4

0.2

0 0

3750

7500

11250 Time (Sec.)

15000

18750

22500

(a) Comparison of gap state esimates with gap state identifcations

1

Gap state identifications Gap state estimates

Gap state

0.8 0.6 0.4 0.2 0 1575

1576.5

1578

Time (Sec.)

1579.5

1581

1582.5

1581

1582.5

Errors between gap state estimates and gap state identifications

1

Theoretically positive bound

Errors

Error

0.5 0 -0.5 -1 1575

Theoretically negative bound 1576.5

1

1578 1579.5 Time (Sec.) (b) Comparison of gap state estimtes with gap state identifications

Gap state

0.6 0.4 0.2 5626.5

5628

Time (Sec.)

5629.5

5631

5632.5

5631

5632.5

Errors between gap state estimates and gap state identifications

1

Theoretically positive bound

Errors

Error

0.5 0 -0.5 -1 5625

5. Conclusion and discussion

Gap state identifications Gap state estimates

0.8

0 5625

plot, each in a short time interval (7.5 s) in Fig. 2(a) (d). From the observations of these figures, it is found that the predictions can track the trend of the gap state development; meanwhile, the accuracy of the predictions are acceptable because the prediction errors almost all fall in the restricted theoretical bounds except a tiny parts in the deleterious phase. Normally, it is recommended that a model be acceptable if the prediction errors in the theoretical bounds are greater than 68% of all the errors. From this criterion, the developed online model is rational, stable, and robust in operation.

Theoretically negative bound 5626.5

5628

Time (Sec.)

5629.5

This paper formulated a two order AR timely varied model for an EDM process, but whereafter proved that the parameter estimate formulation was biased. Then, an instrumental variable approach was proposed to correct the biased estimates of the model parameters. The obtained results are summarized as follows. (1) Two Kalman filters were formulated to implement the instrumental variable approach interactively, one for the estimates of gap states from gap state identifications, the other for the estimates of AR model parameters. (2) The process parameters in this paper are estimated based on the gap state estimates, not on the gap state identifications. Theoretically, this modelling approach is treated rationally and logically. (3) This paper provides not only the both comparisons of

461

Ming Zhou et al. / Procedia CIRP 6 (2013) 456 – 462

1.2

One-step gap state predictions Gap state estimates

1 C

Gap state

0.8

0.6 A 0.4

0.2

0 0

3750

7500

11250

Time (Sec.)

15000

18750

22500

Gap state estimates One-step predictions

Gap state

0.8 0.6 0.4 0.2 0 600

601.5

1

603 604.5 Time (Sec.) Errors between one-step predictions and gap state estimats Theoretically positive bound

Error

606

607.5

606

607.5

Errors

0.5 0 -0.5 -1 600

603

Time (Sec.)

604.5

(B) Comprisn of one-step pedictions with gap state estimtes

Gap state estimates One-step predictions

Gap state

0.8 0.6 0.4 0.2 4876.5

4878

Time (Sec.)

4879.5

1

Theoretically positive bound

Error

4881

4882.5

4881

4882.5

Errors between one-step predictions and gap state estimates

0.5

Errors

0 -0.5 Theoretically negative bound -1 4875

Gap state

0.4 0.2 0 11100

11101.5

11103

Time (Sec.)

11104.5

1

Theoretically positive bound 0.5

Error

11106

11107.5

11106

11107.5

Errors between one-step predictions and gap state estimates Errors

0 -0.5 Theoretically negative bound -1 11100

11101.5

11103

Time (Sec.)

11104.5

(D) Comprison of one-step pedictions wth gap state estimates

1

Gap state estimates One-step predictions

Gap state

0.8 0.6 0.4 0.2 0 20325

20326.5

20328

Time (Sec.)

20329.5

20331

20332.5

20331

20332.5

Errors between one-step preditions and gap state estimates

1

Theoretically positive bound

0.5

Errors

0 -0.5 Theoretically negative bound -1 20325

20326.5

20328

Time (Sec.)

20329.5

Acknowledgements The authors would like to thank the supports from Beijing Natural Science Foundation (No. 4122021), Doctoral Scientific Research Startup Foundation (No. 100901504), Beijing Education Committee Scientific Plan Foundation (No. 051101904) and partly from National Science Foundation (No.51004005).

Theoretically negative bound 601.5

1

0 4875

0.6

26250

(A) Comparison of one-step predictions with gp state estimates

1

Gap state estimates One-step predictions

Fig. 2. Comparison of one-step gap state predictions with gap state estimates in the EDM process. (a) Comparison of one-step gap state predictions with the gap state estimates in phase A in 7.5 s. (b) Comparison of one-step gap state predictions with the gap state estimates in phase B in 7.5 s. (c) Comparison of one-step gap state predictions with the gap state estimates in phase C in 7.5 s. (d) a Comparison of one-step gap state predictions with the gap state estimates in deleterious machining process in 7.5 s.

D B

(C) Comprison of one-step predictions with gap state estimtes

1 0.8

Error

gap state estimates with gap state identifications and one-step predictions with gap state estimates, but also the estimated errors and theoretically error bounds. That almost all the estimated and prediction errors of the models fall in the bounds shows the confidence of this modelling approach. (4) A validation test not only confirms the feasibility of the model but also reveals its robustness in any machining phase of an EDM process. The process model developed can not only be used in a conventional EDM process control, but also in a micro-EDM process control. In fact, the dynamical properties of conventional EDM process and micro-EDM process are almost the same. The only difference between them is that in micro-EDM process the discharging frequencies are higher than in conventional EDM, but with much smaller discharging energy. In other words, there is no problem in applying the developed model or the model in [4] in both situations.

4876.5

4878

Time (Sec.)

4879.5

References [1] Zhou M, Han F. Adaptive control for EDM process with a selftuning regulator. International Journal of Machine Tools & Manufacture 2009;49:462 9. [2] Wang YJ, Rawlings JB. A new robust model predictive control method I: theory and computation. Journal of Process Control 2004;14:231 47. [3] Zhou M, Han F, Tang B. The methodology of detecting dynamical properties of electrical discharge machining (EDM) process towards building a timely varied linear predictive model. In: 16th international symposium on electromachining (ISEM XVI). 2010. p. 73 8. [4] Zhou M, Han F, Soichiro I. A time-varied predictive model for EDM process. International Journal of Machine Tools and Manufacture 2008;48:1668 77. [5] Zhou M, Xianyi Meng et al. Building an EDM process model by an instrumental variable approach based on two interactive Kalman filters. Precision Engineering (2012), http://dx.doi.org/10.1016/j. precision engineering.2012.07.011 [6] Labarre D, Grivel E, Berthoumieu Y, Todini E, Najim M. Consistent estimate of autoregressive parameters from noisy observations based on two interacting Kalman filters. Signal Processing 2006;86:2863 76.

462

Ming Zhou et al. / Procedia CIRP 6 (2013) 456 – 462

[7] Wang KY, Polak E. Identification of linear discrete time systems using the instrumental variable method. IEEE Transactions on Automatic Control 1967;12(December (6)):707 18. [8] Nelson LW, Stear E. The simultaneous on-line estimate of parameters and states in linear systems. IEEE Transactions on Automatic Control 1976;21(February (2)):94 8. [9] Friedlander B, Stoica P, Soderström T. Instrumental variable methods for ARMA parameter estimate. In: IFAC, symposium on identification and system parameter estimation. 1985. p. 29 36. [10] Brown RG, Hwang PYC. Introduction to random signals and applied Kalman filtering. Third edition John Wiley & Sons, Inc.; 1997. [11] Anderson BDO, Moore JB. In: Kailath T, editor. Optimal filtering. Englewood Cliffs, NJ: Prentice-Hall; 1979. [12] Ljung L. System identification theory for the user. Second edition Prentice Hall PTR; 1999. p. 224 6.