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.