Initialization of output-error identification methods – Comparison between ARX and RPM models

Proceedings of the 15th IFAC Symposium on System Identification Saint-Malo, France, July 6-8, 2009

Initialization of output-error identification methods – Comparison between ARX and RPM models Elie Tohme ∗,∗∗ , R´ egis Ouvrard ∗ , Thierry Poinot ∗ , Jean-Claude Trigeassou ∗ , Antoine Abche ∗∗ LAII, University of Poitiers, 40 avenue du Recteur Pineau 86022 Poitiers Cedex, France, (Tel: 33- 5 49 45 35 12; e-mail: [email protected]). ∗∗ Faculty of Engineering, University of Balamand, P.O. Box 100 Tripoli, Lebanon, (e-mail: [email protected])

∗

Abstract: The main disadvantage of the Output-Error (OE) identification methods is that they may converge to a secondary optimum. A good initialization converges to the global optimum. The ARX model is often selected as initialization step to OE algorithms. However, the ARX model may be too biased and may not lead to a good initialization. This paper presents an approach based on the Reinitialized Partial Moment (RPM) to obtain a good initialization for OE methods. The RPM model presents an implicit filter that replaces the necessary explicit filter required by the ARX model. The result are encouraging and they have shown that the RPM based approach has to lead to a better initialization than the conventional techniques. Keywords: ARX model, Global convergence, Initialization of output-error methods, Reinitialized partial moments, RPM model 1. INTRODUCTION

systems. They have pointed out a problem of bad initialization, i.e. a convergence to a secondary optimum.

An Output-Error (OE) method 1 can be summarized as follows: (1) an initial model is defined from the available a priori knowledge of the true system, (2) the (initial) model is simulated using the same excitation as that of the true system, (3) a quadratic criterion defined in terms of the difference between the system output and the model output is formed, (4) The parameters of the model are estimated by minimizing the criterion and a jump to the step (2) is performed until convergence is achieved.

Ljung [2003] has proposed to reduce the bias of the ARX model via a low-pass filter. The Reinitialized Partial Moment (RPM) model, introduced by Trigeassou [1987], presents the advantage of an implicit embedded filter. Hence, the bias of the parameters estimated by this EE approach is reduced, and the RPM model yields a good initialization of OE methods. The main interest of the RPM model is the adaptation to the nature of the noise.

The OE approach provides, under certain hypotheses, an asymptotically unbiased parameter estimation (Landau [1979]). On the other hand, the main disadvantage is its possible convergence to a secondary optimum. Subsequently, this would lead to inaccurate results. In this work, the initialization of OE identification methods in Discrete-Time (DT) domain is studied. The Equation-Error (EE) methods are commonly used for this initialization, e.g. the parameters of an ARX model are estimated to be applied as initial values of the OE algorithms. Rao and Garnier [2002, 2004] and Ljung [2003] have shown that the ARX model is very biased for certain 1 In this contribution, the terms of OE methods and Equation-Error (EE) methods denote the approaches introduced by Landau [1976, 1979], Ljung and S¨ oderstr¨ om [1983]. The EE methods are used to estimate parameters of linear-in-parameters model, and are based on the least squares. The OE methods can be applied to nonlinearin-parameters model, and consist in an iterative optimization procedure.

978-3-902661-47-0/09/$20.00 © 2009 IFAC

302

This paper is organized as follows. The system, the models and the approaches are presented in Section 2. Section 3 presents the simulation results. The disadvantage of the ARX model and the characteristics of the RPM model are highlighted. The bias of both models is studied in Sections 4 and 5. The conclusion is given in Section 6. 2. SYSTEM AND MODELS Consider a stable DT transfer-function system defined by y(k) = G(q, θ0 )u(k) + H(q, θ0 )e(k)

(1)

where {u(k), y(k)} are the DT input/output signals, {e(k)} is a white noise, G(q, θ0 ) and H(q, θ0 ) are DT transfer functions, and θ0 is the true parameter vector. In this work, four different models are studied: • the OE model G(q, θOE ) =

B(q, θOE ) , H(q, θOE ) = 1 F (q, θOE )

(2)

10.3182/20090706-3-FR-2004.0197

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 • the ARX model G(q, θARX ) =

B(q, θARX ) A(q, θARX )

, H(q, θARX ) =

1 (3) A(q, θARX )

• the ARX model with filter L(q) B(q, θARXf ilt ) A(q, θARXf ilt ) 1 ARXf ilt H(q, θ )= L(q)A(q, θARXf ilt )

G(q, θARXf ilt ) =

(4)

• the RPM model with filter M (q) B(q, θRP M ) A(q, θRP M ) 1 RP M H(q, θ )= M (q)A(q, θRP M )

G(q, θRP M ) =

(5)

The polynomials A, B and F are defined by A(q) = 1 + a1 q −1 + . . . + ana q −na , B(q) = b0 + . . . + bnb q −nb (6) F (q) = 1 + f1 q −1 + . . . + fnf q −nf

The ARX model with or without the filter L is well known and is described in (Ljung [1999]). The RPM model, introduced by Trigeassou [1987], is defined by yˆ(k) =

nb X

u bm β m (k) +

na X

(7)

where u βm (k) = q −m M (q)u(k), γ y (k) = (1 − M (q))y(k), ˆ k−n a

αyn (k)

= −q

−n

M (q)y(k), M (q) =

(i + 1)(i + 2) . . . (i + na − mi =

a (na − 1)!An ˆ

k

The selection of kˆ is not more difficult than the selection of the cut-off frequency and the order of the recommended prefilter of an ARX model. For a smaller value than the optimal one, the RPM model is nearer of the ARX model. For a larger value, the noise is not filtered and an accumulation of the error causes a divergence of the variance of the estimator yˆ(K). In practice, the value of kˆ is obtained as follows: kˆ is increased and a standard test, such as quadratic criterion or autocorrelation of the ˆ residuals, is evaluated to find the best k. Both ARX and RPM models are linear in parameters. The least square method is suitable to estimate the parameters. On the other hand, the OE model is nonlinear in parameters and an optimization method must be performed to estimate the parameters. In this context, the least square method is known as an equation-error method, and the optimization method as an output-error method (Landau [1979], Ljung and S¨ oderstr¨om [1983]). 3. RESULTS

an αyn (k) + γ y (k)

n=1

m=0

• If the perturbation is a coloured noise, the reinitialization interval must be within the interval ]na , kˆwn [.

X

i=0 a 1)An ˆ k−i

mi q −i , a = , An i

(8) i! (i − na )!

For more details about the origin of the implicit FIR filter M (q), see Tohme [2008]. This filter gives to the RPM model certain properties. The feature of a finite impulse response allows the removal of the transient effect of an infinite impulse response filter. However, the main property of the RPM model is the adaptation to the nature of the noise. A theoretical bias analysis, presented in Sections 4 and 5, illustrates the properties. The RPM model requires the selection of the design paˆ i.e. the reinitialization interval. The parameter rameter k, ˆ k allows the adaptation to the nature of the noise : • If the perturbation is an output-error white noise (i.e. w(k) = H(q, θ0 )e(k) is a white noise), an optimal (in terms of variance) reinitialization interval exists, namely kˆwn , for which the variance of the error is minimal and the bias is highly reduced. The parameter kˆ ˆ s ] (ts , should be selected such that the interval [0, kt sampling time) is equivalent to the double of the main time constant for an aperiodic system or the double of the (zero to 90%) rising time for an oscillating system. If kˆ = kˆwn , the discrete RPM model is close to an OE model as it is discussed in Section 5. • If the perturbation is an equation-error white noise, i.e. e(k) is a white noise, the reinitialization interval must be equal to na . Consequently, the RPM model is equivalent to an ARX model and the estimation is unbiased.

303

3.1 The Rao-Garnier test system The DT RPM method is a solution to the bias problems of the DT models described by Ljung [2003]. The paper of Ljung is motivated by the contributions of Rao and Garnier [2002, 2004] on the direct and indirect identification of continuous-time models. The proposed example by Rao and Garnier highlights the drawbacks of discretetime approaches, particularly, the ARX, OE, PEM and N4SID algorithms. Ljung has shown that these problems are due to the bias introduced by the ARX model. The ARX model is incorporated to initialize the OE and PEM algorithms, and can be considered as the basis of the N4SID subspace method. The solution proposed by Ljung consists of prefiltering the data. This is a commonly used approach. The DT RPM model does not encounter the same problems as other DT approaches, and consequently, yields a good estimation of the true parameters. Then, the RPM model can be considered as a good step that initializes other DT algorithms such as the OE algorithm. In this context, the ARX model and the RPM model are compared in terms of the initialization of the OE methods. The continuous-time transfer function to be evaluated is given by H(s) =

K (−T s + 1) 2 2 + 2ξ s/ω (s2 /ωn,1 + 2ξ1 s/ωn,1 + 1)(s2 /ωn,2 2 n,2 + 1)

(9)

where K = 1, T = 4s, ωn,1 = 20 rad/s, ωn,2 = 2 rad/s, ξ1 = 0.1 and ξ2 = 0.25. In this study, we compare only the estimated DT model, i.e. there is no comparison with the continuous-time models obtained by the direct or indirect identification (Rao and Garnier [2002, 2004]). The simulation protocol is the same as that considered by Ljung [2003]. The system is simulated for 10s with a sampling period ts = 0.01s. The number of measurements

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 is N = 1000. The system is excited by a zero-order hold multi-sine input signal defined as

4

10

(10)

A DT zero-mean white noise having a SN R = 10dB is added to the measured output. Fifty Monte Carlo experiments are simulated. Figure 1 shows segments of the input/output signals.

2

Amplitude

u(t) = sin(t) + sin(1.9t) + sin(2.1t) + sin(18t) + sin(22t)

10

0

10

−2

10

−4

10

−1

10

0

10

1

10

2

10

3

10

200

Output Phase (degrees)

40 20 0

0

−200

−400

−20

−1

−40

10 1

1.2

1.4

1.6

1.8

2 Time (sec)

2.2

2.4

2.6

2.8

0

10

3

1

10 Frequency (rad/s)

2

10

3

10

Input

Fig. 2. Bode plots of the ARX models (thin solid lines) and the true system (thick dashed line)

4 2 0 −2

4

10 1

1.2

1.4

1.6

1.8

2 Time (sec)

2.2

2.4

2.6

2.8

3

2

Amplitude

−4

10

−4

10

1−

0.001z −4

− + + + 5.774z −2 − 3.814z −3 + 0.951z −4

3.911z −1

(11)

and the same model structure is considered. In the sequel, the ”OE model” will denote the model identified using the OE algorithm to estimate the corresponding parameters.

0

10

1

10

2

10

3

10

200 Phase (degrees)

H(z) =

0.003z −3

−1

10

The corresponding DT system, obtained using zero order hold, is defined as follows 0.003z −2

0

10

−2

Fig. 1. Segments of the input/output signals. The noisefree output is represented by a dashed line

−0.001z −1

10

0 −200 −400 −600 −800 −1 10

0

10

1

10 Frequency (rad/s)

2

10

3

10

Fig. 3. Bode plots of the OE models (thin solid lines) and the true system (thick dashed line)

3.2 Results: ARX model Figure 2 illustrates the Bode plots of the estimated ARX models. As already stated earlier, they are the results of fifty Monte Carlo simulations. The estimated parameters of each ARX model are selected to initialize the OE algorithm. The Bode plots of the corresponding OE models are presented in Figure 3. As noted in (Ljung [2003]), the estimated ARX models are very biased. Hence, the OE models, obtained with the values of the ARX models as initialization of the OE algorithm, do not lead to an accurate estimation of the system. The initialization is too far from the optimal parameters and a convergence to a secondary optimum is observed. The same results are obtained with the PEM and N4SID algorithms (Ljung [2003]). 3.3 Results: RPM model

ARX models. Tohme [2008] shows that the design parameter selection of the RPM model is not too sensitive, and similar results are obtained for 60 ≤ kˆ ≤ 80. Finally, the OE algorithm is initialized using the estimated values of the RPM models. The corresponding Bode plots are shown in Figure 5. It is evident that the bias has disappeared. Thus the RPM model gives better initialization values than the ARX model. The variance observed at frequencies above 30rad/s is due to insufficient excitation (actually there is no excitation above 22rad/s). It is important to note that the estimated ARX model with filter L(q) are quite similar to the RMP model. They are not showed here and can be seen in (Ljung [2003]). 4. BIAS ANALYSIS

The Bode plots of the estimated RPM models are presented in Figure 4. The reinitialization interval kˆ is selected to be 75. It is clear that the bias of the estimated RPM models is much less than the bias of the corresponding

304

Suppose that the true system can be represented by an ARX model defined by (3). Introduce the following 1 notation w(k) = A(q,θ e(k). 0)

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 as an oversampling. This may be an additional reason of the bad results given by the ARX model.

2

10

The selection of a higher sampling period ts = 0.03s respects more the general rule of choosing a sampling period, i.e. 10 samples per period. The simulations were repeated by selecting a sampling period ts = 0.03s and the same number of samples N = 1000. The results of the estimation are almost unchanged. Only the estimated OE models show some improvement. However, they are still largely biased. The effect of a high sampling is a well known problem (Wahlberg and Ljung [1986], Ljung [1999]). The recommended solution is also the prefiltering of the data.

Amplitude

0

10

−2

10

−4

10

−1

10

0

1

10

2

10

3

10

10

Phase (degrees)

0

−200

−400

−600 −1 10

0

1

10

10 Frequency (rad/s)

2

4.2 The RPM model

3

10

10

Consider the limiting RPM model defined by θRP M ∗ = Fig. 4. Bode plots of the RPM models for kˆ = 75 (thin solid lines) and the true system (thick dashed line) 2

10

Let y0 (k) be the undisturbed output signal and let  T u φ(k) = αy10 (k) · · · αyn0a (k)β1u (k) · · · βn (k) b T  w w

(13)

+ α1 (k) · · · αna (k)0 · · · 0 = φy0 ,u (k) + φw,0 (k)

0

Amplitude

lim θRP M

N →∞

10

−2

Assuming the same conditions as in Subsection 4.1, the bias of the RPM model is defined by the following theorem. Theorem 1. (See Tohme [2008]) The asymptotic bias of the limiting RPM model is given by

10

−4

10

−1

10

0

1

10

2

10

3

10

10

Phase (degrees)

0

∗ −1 ∗ RP M ∗ θBIAS = θ0 − θRP M = RRP fRP M M ∗ T ¯ RRP M = Eφ(k)φ (k) ∗ T w ¯ ¯ fRP M = Eφ(k)w(k) − Eφ(k) φw,0 (k)θ0 + γ (k)

−200

where

−400




−600 −800 −1 10

0

1

10

10 Frequency (rad/s)

2

3

10

10

4.1 The ARX model The bias of the ARX model has been studied by Ljung [2003]. The bias of the ARX limiting model is given by ARX θBIAS = θ0 − θARX

∗

2 2 = σw (R + σw D)−1

2 σw



a1 . . . ana 0 . . . 0

T

(12)

where is the noise variance, θ0 is the vector of the true parameters, R is the regression matrix without noise, ∗ θARX is the limiting model and D is a matrix [I 0; 0 0] . This bias is the cause of the bad results of the ARX model (Figure 2) and the OE model (Figure 3). The recommended solution by Ljung is the prefiltering of the data by A(q, θ0 ) (the denominator of the transfer function representing the true system). Thus, the white output error is transformed into a white equation error. Consequently, the ARX model becomes unbiased. The sampling period in the above simulations is small, ts = 0.01s. The period of the oscillating mode ωn,1 = 20rad/s is 0.31s, i.e. 31 times the sampling period. This can be seen

305

.

2 If w(k) is a white noise with a variance σw , the bias becomes
−1 RP M 2 θBIAS =

Fig. 5. Bode plots of the OE models inititialized by the estimated RPM model values (thin solid lines) and the true system (thick dashed line)

(14)

,

RRP M + σw DRP M

FRP M

(15)

¯ y ,u (k)φT RRP M = Eφ y0 ,u (k) 0

DRP M 

= ˆ k−n a

X

2 mi

  i=0  ˆ  k−n  Xa  mi−1 mi   i=1  .  . .   k−n ˆ  Xa  mi−na +1 mi   i=na −1  0   . . . 0

ˆ k−n a

ˆ k−n a

X

X

mi−1 mi . . .

X

mi−na +1 mi 0 . . . 0

i=na −1 ˆ k−n a

i=1 ˆ k−n a

X

2

mi

i=0

. . mi−na +2 mi .

i=na −2

..

.

. . .

. . .

ˆ k−n a

...

...

X

2

mi

...

...

ˆ k−n a

X

2 mi

    .  . .     .  . .    0    0 .  .

...

0

0 ...

...

. . . 0

. . . . 0 ... 0

FRP M =



0 ...

i=0

ˆ k−n a

ˆ k−n a

X

X



mi−1 mi  mi−na +1 mi + + . . . + ana  a1     i=1 i=na −1 i=0   .   . .   ˆ ˆ ˆ   k−na k−na 2 X X Xa σw  k−n  2  a1 mi−na mi  mi + mi−na +1 mi + . . . + ana    i=na −1  i=na i=0   0   .   . . 0



15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 Consider that the system (9) is simulated with a white noise input having a unit variance. Figure 6 shows the ARX T ARX 0.5 norm, defined by (θBIAS θBIAS ) , of the bias of the limiting ARX model defined by (12), as well as the bias of the limiting RPM model defined by (15) for different ˆ It illustrates the values of the reinitialization interval k. 2 . dependence of the norms on the variance σw

DT form as given by (5) (Tohme [2008]). Consequently, by using the same procedure outlined in (Wahlberg and Ljung [1986]), the frequency distribution of the bias of the RPM model can be determined using the following weighting function

The results show that the bias of the parameters of the RPM model is much less than the bias of the parameters of the ARX model. It is to be noted that for kˆ = na = 4, the norm of the RPM model is superimposed over the norm of the ARX model. By progressively increasing the value ˆ the bias decreases until the norm of the bias reaches of k, the smallest value at kˆ ≈ 75. In fact, the norms are almost superimposed for 60 ≤ kˆ ≤ 90.

The analysis can be taken further by simulating the system that is defined by Equation (9). The input signal is assumed to have a spectrum φu (ω) ≡ 1. Figure 7 shows QARX (ω, θ0 ), QOE (ω, θ0 ) and QRP M (ω, θ0 ) for various conditions. The plot of QARX (ω, θ0 ) corresponds to the ARX model without prefiltering, i.e. L(eiω ) ≡ 1. It shows that the high frequencies have a weight factor 1012 times larger than the low frequencies. A low-pass filter L(q) allows a better distribution of the frequency weighting. In the context of a white output error, the ideal filter is L(q) = 1/A(q, θ0 ), i.e. the ideal filter of the algorithm of Steiglitz and McBride [1965]. It is equivalent to an OE model and the weighting function is QOE (ω, θ0 ) ≡ 1. The plot of QRP M (ω, θ0 ) for kˆ = 75 shows that the RPM model is closer to the OE model, i.e. the filter M (q) is closer to the ideal filter of Steiglitz-McBride.

2

10

ARX norm(θ BI AI S ) 0

10

RP M norm(θ BI AI S )

−2

10

kˆ = 10



∗

2



∗

2

QRP M (ω, θRP M ) = M (eiω ) φu (ω) A(eiω , θRP M )

(17)

−4

10

kˆ = 150

4

10

kˆ = 40

−6

10

QRP M (ω, θ 0 )

kˆ = 75 2

kˆ = 150

10

QOE (ω, θ 0 ) ≡ 1

−8

10

0

10

−10 −10

−8

10

−6

10

−4

10

−2

10

10

10

0

10

kˆ = 40

Power

10

kˆ = 75

−2

2 σw

¯ ¯2 QARX (ω, θ 0 ) = ¯A0 (ejω )¯

−4

10

ARX RP M Fig. 6. The dependence of the norms of θBIAS and of θBIAS 2 on σ for different kˆ

−6

10

kˆ = 10

w

−8

10

5. FREQUENCY DISTRIBUTION OF THE BIAS −10

10

Wahlberg and Ljung [1986], Ljung [1999] have proposed another approach to analyze the bias and its frequency distribution. In the case of the ARX model (4), the frequency distribution of the bias is determined using the weighting function ∗



2



∗

2

QARX (ω, θARXf ilt ) = L(eiω ) φu (ω) A(eiω , θARXf ilt ) (16)

where φu (ω) is the spectrum of the input signal, and L(q) is the recommended prefilter to solve the problems of the white output error and the fast sampling. The above equation is obtained for the one-step-ahead prediction. Actually, as it is pointed out in (Ljung [1999]), the lim∗ iting vector of parameters θARXf ilt is a result of a compromise between the minimization of the difference between the model G(eiω , θARXf ilt ) and the true system G(eiω , θ0 ) weighted by QARX (ω, θ) on one hand and fitting H(eiω ,θ0 ) 2 L(eiω ) to represent the spectrum of the output error on the other hand. Thus, it is of great importance to study the weighting function QARX (ω, θ). Consider a system that is represented by the RPM model output (7). The RPM model can be rewritten in a standard

306

−2

10

−1

0

10

10

1

10

Frequency (rad/s)

ˆ QARX (ω, θ0 ) using an Fig. 7. QRP M (ω, θ0 ) for different k, ARX model without filtering and QOE (ω, θ0 ) Figures 8 and 9 show the dependence of the sensitivity of the ARX and RPM models on the choice of their respective design parameters, i.e. the bandwidth of the filter L(q) for the ARX model and the reinitialization interval kˆ for the RPM model. The plot of QARX (ω, θARXf ilt ) is obtained using a 4th order low-pass Butterworth filter L(q), with a cutoff frequency ωc . It can be noticed that the frequency distributions of both weighting functions are similar. The choice of the design parameter of the RPM model is equivalent to that of the ARX model with prefiltering. The selection of kˆ between [60, 90] and ωc between [8rad/s, 12rad/s] yield similar results. The best values are ωc = 10rad/s and kˆ = 75. The order of the Butterworth filter L(q) is another design parameter to be taken into account in the case of the ARX model. Actually, the best results are obtained if the filter is assumed to have the same order as the true system.

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009

3

3

10

10

2

2

10

ωc = 8rad/s

1

10

Order6, ωc = 24rad/s Order5, ωc = 18rad/s

1

10

ωc = 10rad/s Power

Power

10

ωc = 6rad/s

0

10

−1

0

10

−1

10

10

Order4, ωc = 10rad/s ωc = 12rad/s −2

−2

10

10

ωc = 14rad/s

−3

10

−3

−2

−1

10

0

10

10

10

1

10

3

10

2

10

kˆ = 60 kˆ = 65

1

Power

kˆ = 75 kˆ = 80 kˆ = 85 kˆ = 90

−1

−2

10

−3 −2

10

0

10

1

10

Fig. 10. The dependence of QARX (ω, θARXf ilt ) on the order of the Butterworth filter The design parameter kˆ can be selected such that the RPM model is more efficient in the removal of the noise. For a white noise output error, the RPM model is close to the OE model. For a white noise equation error, the RPM model is equivalent to the ARX model. REFERENCES

kˆ = 70 0

10

10

−1

10

Frequency (rad/s)

Fig. 8. QARX (ω, θARXf ilt ) for various cut-off frequencies of a low-pass Butterworth filter

10

−2

10

Frequency (rad/s)

10

Order3, ωc = 1rad/s

−1

0

10

10

1

10

Frequency (rad/s)

Fig. 9. QRP M (ω, θRP M ) for various design parameters 60 ≤ kˆ ≤ 90 Figure 10 shows the dependence of QARX (ω, θARXf ilt ) on the order of the Butterworth filter. Since the true system is a 4th order, the best result is obtained using a 4th order filter. For each plot, the best cutoff angular frequency ωc has been chosen. 6. CONCLUSION The initialization of OE algorithms is crucial because it could lead to a convergence to the global optimum or to a secondary optimum. The EE algorithms are often used as an initialization procedure to OE algorithms. In DT domain, the ARX model is often selected. However, it must be used with precaution because it may be too biased to yield a good initialization. Hence, OE algorithms could converge to a secondary optimum. The RPM model is an alternative approach to achieve good initialization. The implicit filter embedded in the RPM model is characterized by better properties than the necessary explicit filter required by the ARX model. The Monte Carlo simulations show that the RPM model yields a good initialization for the OE algorithm. Consequently, the latter will converge to the global optimum.

307

I.D. Landau. Unbiased recursive identification using model reference adaptive techniques. IEEE Trans. Automat. Contr., AC-21(2):194–202, April 1976. I.D. Landau. Adaptive control : the model reference approach. Marcel Dekker, Control and systems theory, Vol. 8, 1979. L. Ljung. Initialisation aspects for subspace and outputerror identification methods. In European Control Conference, Cambridge, U.K., 2003. L. Ljung. System identification - Theory for the User. Prentice-Hall, Upper Saddle River, N.J., 2nd edition, 1999. L. Ljung and T. S¨ oderstr¨om. Theory and practice of recursive identification. MIT Press, London, England, 1983. G.P. Rao and H. Garnier. Numerical illustrations of the relevance of direct continuous-time model identification. In In Proceedings of 15th IFAC world congress, Barcelona, Spain, 2002. G.P. Rao and H. Garnier. Identification of continuoustime systems: direct or indirect? Systems Science, 30-3: 25–50, 2004. K. Steiglitz and L. McBride. A technique for the identification of linear systems. IEEE Transactions on Automatic Control, 10-4:461–464, 1965. E. Tohme. Initialization of Output Error Identification Algorithms. PhD thesis, Universit´e de Poitiers, France, http://laii.univ-poitiers.fr/THESES/elie tohme.pdf , 2008. J.C. Trigeassou. Contribution a ` l’extension de la m´ethode des moments en automatique. Application a ` l’identification des syst`emes lin´eaires. Th`ese d’etat, Universit´e de Poitiers, France, 1987. B. Wahlberg and L. Ljung. Design variables for bias distribution in transfer function estimation. IEEE Transactions on Automatic Control, AC-31:134–144, 1986.