A bias correction method for fractional closed-loop system identification

Journal of Process Control 33 (2015) 25–36

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

A bias correction method for fractional closed-loop system identification Z. Yakoub, M. Chetoui, M. Amairi ∗ , M. Aoun University of Gabes, National Engineering School of Gabes, Tunisia Research Unit Modeling, Analysis and Control of Systems (MACS), 06/UR/11-12, Tunisia

a r t i c l e

i n f o

Article history: Received 27 November 2014 Received in revised form 16 April 2015 Accepted 30 May 2015 Keywords: Closed-loop system Commensurate-order Fractional differentiation Least squares State variable filter System identification

a b s t r a c t In this paper, the fractional closed-loop system identification using the indirect approach is presented. A bias correction method is developed to deal with the bias problem in the continuous-time fractional closed-loop system identification. This method is based on the least squares estimator combined with the state variable filter approach. The basic idea is to eliminate the estimation bias by adding a correction term in the least squares estimates. The proposed algorithm is extended, using a nonlinear optimization algorithm, to estimate both coefficients and commensurate-order of the process. Numerical example shows the performances of the fractional order bias eliminated least squares method via Monte Carlo simulations. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Recently, several researches show that the fractional calculus provides an excellent tool to describe the behavior of many complex physical systems due essentially to their long memory characteristic. This induces the use of fractional models in many applications involving various theoretical fields such as control [1,2], diagnosis [3] and system identification [4,5]. Concerning the time-domain fractional system identification, it has received great interest in the late nineties. In Cois et al. work [6] the fractional differentiation orders are assumed to be known a priori and only the coefficients are estimated by minimizing the quadratic criterion based on the equation-error. A synthesis of fractional Laguerre basis for system identification is developed in [7]. An optimal instrumental variable method for continuoustime fractional system identification has been proposed in [8]. In Chetoui et al. work [4] a new methods based on higherorder statistics are illustrated for the estimation of both fractional orders and coefficients of continuous-time errors-in-variables fractional models. In the bounded error context, Amairi et al. works extended several set-membership methods to the fractional system

∗ Corresponding author. Tel.: +216 23268539. E-mail addresses: [email protected] (Z. Yakoub), [email protected] (M. Chetoui), [email protected] (M. Amairi), [email protected] (M. Aoun). http://dx.doi.org/10.1016/j.jprocont.2015.05.005 0959-1524/© 2015 Elsevier Ltd. All rights reserved.

identification such as the outer bounding ellipsoid (OBE) [9] and the outer bounding parallelotope (OBP) [10]. The methods mentioned above are developed in open-loop conditions. However, for many industrial production process, the experimental data can only be obtained in closed-loop conditions for several reasons like stability, security, safety and performance constraints. In the literature, three approaches are proposed for the closed-loop system identification: the direct approach, the indirect approach and the joint input/output approach [11]. As for rational systems, (see [12,11,13] and the references therein for more details) the fractional closed-loop system identification has also attracted an attention recently. In [14], an unstable fractional first-order system with input time-delay has been identified using a graphical method based on the step response of the fractional closed-loop system. A more general indirect approach has been proposed in [15] where the least squares algorithm combined with the state variable filter is used. Simulation results have shown the presence of a bias on the estimates when the system is contaminated by a high level additive white noise. To solve this problem, the bias eliminated least squares (bels) method proposed in [16] is extended to the fractional case in [17]. This method is called the fractional order bias eliminated least squares (fobels). The fobels method is composed by three major steps. In the first step, the least squares method is used to estimate the closed-loop parameters. In the second step, the bias introduced by this method is estimated by an appropriate algorithm. The third step consists in computing the process parameters using the bias correction principle. The algorithm proposed in [17], besides its efficiency, it

26

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

requires a restriction on the controller order. To remove this restriction, like as the rational case (see [18] and [19] for more details), we propose in this paper an extended version of the fobels method called the prefileterd fobels (pfobels) which deals with an arbitrary order controller. The idea is to use a stable prefilter with an appropriate order connected to the external excitation of the closed-loop system. This paper is organized as follows: a mathematical background of fractional systems is presented in Section 2. The problem statement is illustrated in Section 3. In Section 4, the extension of the bias eliminated least squares method to the fractional closed-loop systems is developed. In Section 5, the algorithm of the prefiltered fobels method is described. Both coefficients and fractional orders are estimated in the Section 6. In Section 7, a numerical example shows the performances of the developed method. Finally, Section 8 concludes this paper. 2. Mathematical background Several definitions of the fractional differentiation have been proposed in the literature [20–22]. In this paper only the GrünwaldLetnikov definition is used [23]. Definition 1. The −Grünwald fractional differentiation of a continuous-time function f(t) relaxed at t = 0 (i.e. f (t) = 0, ∀t ≤ 0) is defined by 1  (t)   (−1)k h K

D f

   k

k=0

f (t − kh) ,

∀t ∈

R∗+

(1)

d

R∗+ denotes the set of all strictly positive real numbers. D = dt is the timedomain  differential operator, h is the sampling period, t = Kh and



is the Newton’s binomial generalized to fractional

k orders such as



   k

=

if k = 0  ( − 1) ( − 2) . . . ( − k + 1) if k > 0 k!

1

(2)

Consider a SISO commensurate-order1 fractional system described by na 

ai Di y (t) =

i=0

nb 

bj Dj yc (t)

(3)

j=0

where (ai , bj ) ∈ R2 are the linear coefficients of the differential equation,  ∈ R∗+ is the commensurate-order and where yc (t) and y (t) are respectively the input and the output signals. Applying the Laplace transform to the fractional differential Eq. (3), under zero initial conditions, yields the fractional transfer function B (s) G (s) = = A (s)

nb

bs j=0 j

na

as i=0 i

j

i

(4)

Stability condition for the fractional systems has been established in [21]. Theorem 1. A commensurate-order system described by (4) is Bounded Input Bonded Output (BIBO) stable iff 0<<2

(5)

1 All differentiation orders are exactly divisible by the same number, an integral number of times.

e(t)

yc (t) ys (t)

+

C(s)

-

+ +

G(s)

y0 (t)

+ +

y(t)

Fig. 1. Fractional closed-loop system.

and

∀sk ∈ C,

 

A sk = 0 such that

   arg s >   k

2

(6)

where sk is a pole of the commensurate transfer function. 3. Problem statement Consider the fractional closed-loop system depicted in Fig. 1. G(s) and C(s) describe respectively the process transfer function and the controller transfer function. The signals ys (t), yc (t) and y0 (t) are respectively the set-point input, the continuous-time noise-free input and the process output. The measurable output signal y(t) is eventually corrupted by an additive noise e(t) such as y (t) = y0 (t) + e (t)

(7)

The process transfer function is considered commensurate and given by B (s) G (s) = = A (s)

nb

bs j=0 j

na

as i=0 i

j

i

;

na ≥ nb

(8)

where (ai , bj ) ∈ R2 are the linear coefficients of the process transfer function and  ∈ R∗+ is the process commensurate-order. The controller transfer function is assumed to be known and described by a commensurate transfer function as C (s) =

Q (s) = P (s)

nq qr sr ; r=0 np l ps l=0 l

np ≥ nq

(9)

where (pl , qr ) ∈ R2 are the linear coefficients of the controller transfer function and  ∈ R∗+ is the controller commensurate-order. In this paper, the input signal yc (t) is considered perfectly known and the output signal y(t) is measured. The controller output (the control signal) is supposed unmeasurable. Thus, the indirect approach is required to identify the process with fractional model using a prior knowledge of the controller. Due to the feedback control, there exist a correlation between the unmeasurable output noise and the control signals [11]. In the case where the process and the controller are represented respectively by (8) and (9), the correlation is amplified due to the long memory aspect of the fractional differentiation. Recently, some works present contributions in the fractional closed-loop system identification context [15,17,24]. The simulation results presented in [15] have shown that for an important additive noise the estimated parameters are biased. To eliminate this bias, the optimal instrumental variable method combined with a nonlinear optimization algorithm is handled to identify both fractional transfer function coefficients and fractional orders [24]. An extension of the bias eliminated least squares (bels) method to fractional order case has been proposed in [17]. This method gives an unbiased estimation but, a restriction on the controller order is verified. The objective of this paper is to extend the bels method to fractional order case to identify the fractional closed-loop systems without noise modelling and without any restriction on the controller order.

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

4. Fractional system parameters estimation In the case of a zero set-point, the fractional closed-loop system presented in Fig. 1 can be described by y (t) = T (s) yc (t) + S (s) e (t)

(10)

27

where ˛ is the fractional order differentiation, the order Nf is an integer chosen such that Nf > ˛N  and wf denotes the cut-off frequency. Using the filtered input and output signals denoted, respectively, by ycf (t) and yf (t), equation (15) can be rewritten as 

where T (s) =

yf (t) +

G (s) 1 + C (s) G (s)

(11)

N 



An D

˛n

yf (t) =

n=1

M 



Bm Dˇm ycf (t) + f (t)

(18)

m=0

4.1. Ordinary least squares (ols) method

and S (s) =

1 1 + C (s) G (s)

are respectively the fractional closed-loop transfer function and the sensitivity function. In the rest of paper, the following assumptions are required to the the indirect closed-loop system identification [16,18,19]: A1–

The continuous-time noise-free input signal yc (t) is measurable, stationary and persistently exciting. The output noise e(t) and the input signal yc (t) are uncorrelated. The controller coefficients and the fractional commensurate-order  are known. The fractional closed-loop system is asymptotically stable. The polynomials A(s) and B(s) are coprimes, and so are P(s) and Q(s).

A2– A3– A4– A5–

Consider the closed-loop parameters vector  defined by

(12)

 = [A1 , . . ., AN  , B0 , . . ., BM  ]T

(19)

and the filtered regression vector ϕf described as



ϕfT (t) = −D

˛

1

yf (t) , . . ., −D

˛  N y f

(t) ,

D

ˇ

0

ycf (t) , . . ., D

ˇ  M yc f

(t) (20)

Eq. (18) can be reformulated as a linear regression form yf (t) = ϕfT (t)  + f (t)

(21)

The parameters vector  can be estimated by minimizing the ordinary least squares (ols) criterion defined as follows



JNs

ˆ ols 





1  2 ˆ ols f tk ,  Ns Ns −1

=



(22)

k=0

Substitution of expressions (8) and (9) into the equation (10) yields (A (s) P (s) + B (s) Q (s)) y (t) = B (s) P (s) yc (t) +  (t)

(13)

ˆ ols (Ns ) = Rˆ −1 Rˆ ϕ y  ϕf ϕf f f

where (t) is the equation-error defined by  (t) = A (s) P (s) e (t)

(14)

y(t) +



An D

˛n

y (t) =

n=1



Bm D

 ˇm

yc (t) +  (t)

where (An , Bm ) ∈ are the linear coefficients of the closed-loop transfer function and = (na + 1)(np + 1) − 1 = (nb + 1)(np + 1) − 1

A0 = 1,

 ˇ0 < ˇ1 < · · · < ˇM 



1+

s wf

Nf

1  ϕf (tk ) yf (tk ) Ns Ns −1

ˆ ols (Ns ) =  + Rˆ −1 Rˆ  ϕf ϕf ϕf f

⎡

−1 =  + Rˆ ϕ f ϕf

⎣

Rˆ yf f Rˆ yc

f

(25)

⎤ ⎦

(26)

f

ols

(16)

where Tˆ (s) denotes the estimate of the fractional closed-loop transfer function. A linear transformation can be applied to the input and the output signals to reduce the noise effect. It uses the extension of the Poisson’s filter to fractional case defined by [25] s˛



=  +  ˆ

The prediction of the noise-free output signal is given by

F˛ (s) =



Rˆ ϕf yf = E ϕf yf =

where E denotes the mathematical expectation. Replacing yf (t) by its expression (21) in (23) leads to

 ∈ R are For identifiability purpose, the orders ˛n ∈ R+ and ˇm + ordered such as

yˆ 0 (t) = Tˆ (s) yc (t)

(24)

k=0

k=0

˛0 = 0

˛1 < ˛2 < · · · < ˛N  ;

Ns −1

and Rˆ ϕf yf is the cross-covariance vector represented by:

R2

M

1  ϕf (tk ) ϕfT (tk ) Ns

(15)

m=0

N



Rˆ ϕf ϕf = E ϕf ϕfT =

M

(23)

where Rˆ ϕf ϕf is the auto-covariance matrix defined by



Eq. (13) can be rewritten as N

where the input and output signals are uniformly sampled at Ns instants {tk = kh}k=0 and Ns is the number of samples. ˆ ols is given by The optimal parameters vector 

(17)

where  ˆ

ols

is the bias of the ols estimator.

In the closed-loop, ϕf and  f are always correlated which justify the presence of the bias in the ols estimates ( = / 0) [15]. ˆ ols

4.2. Fractional order bias eliminated least squares method fobels In order to obtain an unbiased estimates, the bias eliminatedleast-squares (bels) method developed in [19] is extended to the fractional case and called the fractional order bias eliminated least squares (fobels) method.

28

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

Thus, the fobels estimator is given by

and

ˆ ols −  ˆ fobels =   ˆ 

ols

(27)

ˆ ols − Rˆ −1 Rˆ = ϕf ϕf ϕf f

Since ycf (t) and  f (t) are uncorrelated (see the assumption A2), the cross-covariance vector is given by:

 Rˆ ϕf f =

Rˆ yf f



0

= Rˆ yf f

(28)

where



=

IN 

0





∈ R(N +M +1)×N



(29)

and I is the identity matrix. The relation between the closed-loop and the process parameters vector is given by

 = M + 

(30)





q0

0

⎢. ⎢. ⎢. ⎢ ⎢ qnq ⎢ ⎢ ⎢ ⎢0 ⎢ ⎢. ⎢. ⎢. ⎢ ⎢ ⎢ ⎢ ⎢ Qc = ⎢ ⎢ ⎢ ⎢ .. ⎢. ⎢ ⎢ ⎢ .. ⎢. ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢. ⎣ ..

...

...

.. .

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ .. ⎥ ⎥ . ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 0 ⎥ 0(np −nq ) ⎥ ⎥ ⎥ ⎥ ⎥ q0 ⎥ ⎥ .. ⎥ ⎥ . ⎥ ⎥ qnq ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ .. ⎦ . .. .

0 0(np −nq ) ..

q0

.

.. .

..

.

qnq 0 ..

.

.. .

..

0

0

⎤

0

...

.

...

0



Pc1 ∈ RN ×na contains the last N rows and the na columns of Pc ,  Qc1 ∈ RN ×(nb +1) contains the last N rows and (nb + 1) columns of Qc  and Pc2 ∈ R(M +1)×(nb +1) contains the first (M + 1) rows and (nb + 1) columns of Pc . 

(31)

HT M = 0

W = I(N  +M  +1) − MM +

=

p0

...

pnp

0

...

...

0

T   ∈ R(N +M +1)

(32)

where by

M+

M=

Pc1

Qc1

0

Pc2

p0

⎢ ⎢ ⎢0 ⎢ ⎢ PcT = ⎢ .. ⎢. ⎢ ⎢ .. ⎣.

0

−1

MT

(38)

... ... ..

(33)

pnp 0

0 p0

.

 R(N +1)×(na +1)

... ...

... pnp

.

..

0

...

···

..

. 0

and Qc ∈

0

 R(N +1)×(nb +1)

⎤

⎥ ⎥ ⎥ ... ⎥ ⎥ .. ⎥ .. . . ⎥ ⎥ ⎥ 0 ⎦

.

p0

...

pnp



WM = I(N  +M  +1) M − M M T M

(39)

−1

MT M

(40)

=0 There exist a non-singular transformation matrix U such that the last columns of WU are linearly independent. So, the matrix H is defined as follows



H = WU



0

(41)

Inp (na +nb +2)



.. .

0

0 ..

rank(W ) ≥ (np (na + nb + 2)) and



Consider the matrices Pc ∈ defined respectively by

⎡

(37)

The rank of the matrix W verifies



M ∈ R(N +M +1)×(na +nb +1) is a known full column rank matrix described as follows



    R(N +M +1)×(N +M +1)

denotes the generalized inverse of the matrix M defined

M + = (M T M) 



(36)

Proof 1. Consider the square matrix W ∈ defined as follows

and



(35)

Proposition 1. There exist a matrix H ∈ R(N +M +1)×(np (na +nb +2)) with full column rank such that

where is the process parameters vector defined as follows

= a1 , . . ., ana , b0 , . . ., bnb

⎡

䊏 (34)

Proposition 2.



If np nb ≥ na − np the equation



−1 ˆ ols −  H T Rˆ ϕ Rˆ yf f = H T  f ϕf

has a unique solution.

(42)

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

Multiplying equation (30) by HT leads to

Proof 2. HT  = H

 T

M + 

Eq. (13) is rewritten as follows



(43)

Using Eq. (36), Eq. (43) is rewritten as H  = HT 

(44) P˜ (s) = P (s)

The ordinary least squares estimator is described by ˆ ols =  +   ˆ 

ols

=

(45)

Rˆ yf f

Multiplying Eq. (45) by HT and using Eq. (44) yields

H

T

−1 Rˆ ϕ f ϕf

(A (s) P (s) + B (s) Q (s)) y (t) = B (s) P˜ (s) y˜ c (t) +  (t)

Rˆ yf f = H

T



ˆ ols −  



Thus, if the matrix

D (s) N (s)

(52)

The augmented fractional closed-loop system (51) can be identified using the identification data (˜yc (t), y(t)). So, Eq. (51) can be rewritten as a linear regression form ˜ + f (t) yf (t) = ϕ ˜ fT (t) 

(46)

−1 H T Rˆ ϕ f ϕf

(51)

where

T

−1  + Rˆ ϕ f ϕf

29

is a full column rank, the

estimate of the cross-covariance vector Rˆ yf f in the sense of least-

(53)

where ϕ ˜ f is the filtered regression vector defined as



ϕ ˜ fT (t) = −D

˛

1

yf (t) , . . ., −D

˛  N y f

(t) , D

ˇ

0



(54)

squares has a unique solution given by:



+

−1 Rˆ yf f = H T Rˆ ϕ f ϕf



ˆ ols −  HT 

(47)

Eq. (46) is determined by solving the above system of linear equations, in which there are (np (na + nb + 2)) equations (the −1 elements of the (np (na + nb + 2)) × N matrix H T Rˆ ϕ ) and N f ϕf unknowns (the elements of the Rˆ yf f ). Using definition of the matrix H and , the matrix

−1 H T Rˆ ϕ f ϕf

is a submatrix composed by the last N

columns of the non-singular np (na + nb + 2) rows and the first −1 . So, if np nb ≥ na − np (the number of rows is matrix (WU)T Rˆ ϕ f ϕf greater than the number of columns) this matrix is a full column T −1 T ˆ −1 rank. Using assumption A1, the matrix ((H T Rˆ ϕ ϕ ) (H Rϕ ϕ )) is f

f

f

f

and





M = np + nf + 1 (nb + 1) − 1

(55)

= (np + nf )(nb + 1) + nb

(56)

The augmented closed-loop parameters vector given by



˜ = A1 , . . ., AN  , B˜ 0 , . . ., B˜ M 

The process parameters are computed from the fobels estimates and the prior knowledge of the controller transfer function via the following expressions



ˆ fobels − 

ˆ fobels = M + 



= MT M

−1



ˆ fobels −  MT 

T

(57)

can be estimated by applying the ols estimator such as ˜ˆ ols (Ns ) =  ˜ +  ˜ˆ

ols

invertible. Under these conditions, an unique solution is obtained.  䊏

y˜ cf (t) , . . ., DˇM y˜ cf (t)

(58)

˜ + Rˆ −1 ˜ Rˆ y  = ˜ ϕ ˜ ϕ f f f

f

Since y˜ cf (t) and  f (t) are independents (see assumption A2), the cross-covariance vector is described by

 Rˆ ϕ˜ f f =



Rˆ yf f

˜ Rˆ y  = f f

0

(59)

(48) where



IN 

0

 

∈ R(N +M

+1)×N 

4.3. Prefiltered fractional order bias eliminated least squares method pfobels

˜ =

If np nb < na − np , the number of equations in (47) is not sufficient to uniquely determine Rˆ yf f (the number of equations is fewer

The relation between the augmented closed-loop and the process parameters vector is rewritten as follows

than the number of unknowns np (na + nb + 2) < N ). Thus, to solve this problem, an idea similar to Garnier’s idea proposed for rational systems in [19], is considered in this paper for fractional systems. The idea consists on using a stable prefilter of an order nf connected to the input signal and defined by



F  (s) =

N (s) = D (s)

nf

1 1+

s f

(49)

˜ = M

˜ + ˜ 

(61) (N  +M +1)×(na +nb +1)

˜ ∈R is a known full-column-rank where M matrix (rank(M) = na + nb + 1) defined by



˜ = M

Pc1

Qc1

0

P˜ c2



(62)

where

The filtered input signal is given by y˜ c (t) = F  (s) yc (t)

(60)

(50)

So, the objective is to obtain an augmented fractional transfer function and so a sufficient number of equations to solve the estimation problem.

⎡  ˜ =⎣

p0

...

pnp

0

...

...

0

⎤T ⎦ ∈ R(N  +M +1)

(63)

30

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

and P˜ c2 ∈ R(M

⎡p

0

⎢0 ⎢ ⎢ ⎢. ⎢ .. P˜ cT2 = ⎢ ⎢. ⎢. ⎢. ⎢ ⎣0

+1)×(n +1) b

is defined as follows

p˜ 1

...

p˜ np

0

...

...

0

p0

..

... p˜ 1 ..

.

...

0

...

0 ..

.

...

...

p˜ np

0

...

...

.

0 0

p0

p˜ 1

⎤

⎥ ⎥ ⎥ .. ⎥ .. . . ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ . . . p˜ np ⎦

• Step 1: design a stable prefilter F (s) as defined by Eq. (49) to obtain the prefiltered input signal y˜ c (t); ˜ using Eq. (66); • Step 2: compute the matrix H • Step 3: use y˜ c (t) and y(t) signals to estimate the augmented ˜ˆ ols and compute Rˆ ϕ˜ ϕ˜ ; closed-loop parameters vector  • Step 4: solve Eq. (68) to get Rˆ y  ; f f • Step 5: compute the asymptotic bias given by

ols

f

f

−1 ˜ ˆ Ryf f = Rˆ ϕ ˜ ˜ ϕ f

=

with = p0 +

pi , f

f

−1 ˆ = Rˆ ϕ Rϕ˜ f f ˜ ¯ ϕ

 ˜ˆ

(64)

p˜ i

f

f

Rˆ −1 ˜˜ ϕ ϕ f ˜f



+

˜ H T Rˆ −1 ˜ ϕ ˜ ϕ f

f



˜ˆ ols −  ˜T  H ˜

The augmented closed-loop unbiased parameters vector is given by

1 ≤ i ≤ np − 1

pnp

p˜ np =

˜ˆ fobels =  ˜ˆ ols −   ˜ˆ

f

ols

There exists (N + M + 1) × (N + M − na − nb ) full column-rank ˜ such that matrix H

• Step 6: compute the process parameters vector using Eq. (71).

˜ =H ˜ T ˜ T ˜ TM ˜ = 0 and H ˜ H

6. Fractional orders estimation

(65)

˜ is defined as follows The matrix H



˜ =W ˜U H



0

(66)

I(N+M −na −nb ) 

˜ ∈ R(N +M where the matrix W

+1)×(N  +M +1)

is defined by

˜ = I(N  +M +1) − M ˜M ˜+ W

˜ = 



+

f

f



˜ˆ ols −  ˜T  H ˜

(68)

The above expression represents a system of linear equations, in which there are N unknowns (elements of Rˆ yf f ) and ˜ So, in order to obtain an ˜ T Rˆ −1 ). (N + M + 1 − na − nb ) (rows of H ˜f ϕ ˜f ϕ

unique solution, it is required that (N + M + 1 − na − nb ) ≥ N . This condition is equivalent to: nf ≥

na − np nb + 1

(69)

na (nb +1)

− np , the expression (68) contains exactly N

If nf =

unknowns and

N

equations and the cross-covariance vector Rˆ yf f

is rewritten as follows: Rˆ yf f =



˜ T Rˆ −1 H ˜f ϕ ˜f ϕ

−1

˜

˜T H



−1

. . .,

AN  ,



B0 ,

. . .,

BM  ,

˛1 ,

. . .,



˜ˆ ols −   ˜



˜ˆ fobels −  ˜T  M ˜



˜ = A1 , 

...

5. Algorithm of the pfobels method The developed pfobels algorithm is composed by 6 steps:

(71)

. . .,

 ˇM



AN  ,

B0 ,

. . .,

BM ,





(73)

The main idea of the proposed algorithm consists in estimating the commensurate-order by minimizing the bias introduced by the ordinary least squares estimator combined with the state variable filter with respect to  [4,25]. This algorithm is called commensurate-order optimization combined to pfobels (copfobels). The quadratic criterion is defined as follows





=





1 2  ˜ˆ    2 ols

where the bias  ˜ˆ  ˜ˆ

ols

(74)

is given by:

−1 ˆ Rϕ˜ f f = Rˆ ϕ ˜ ˜ ϕ f

f

−1 = Rˆ ϕ Rˆ yf f ˜ ˜ ϕ f

ˇ0 ,

This vector is composed by 2(N + M + 1) elements and it can be reduced to 2(N + M + 1) elements if np nb ≥ na − np . Using the indirect approach, which supposes the knowledge of the controller transfer function, only the process commensurateorder  is estimated. In this case, the number of parameters can be reduced to N + M + 2. The closed-loop parameters vector is rewritten as follows

ols

(70)

˛N  ,

(72)

ols

˜ˆ fobels −  ˜+  ˜

ˆ fobels = M



A1 ,

J  ˜ˆ

Remark 1. If np nb ≥ na − np ; then, the prefilter is non langer needed (F (s) = 1 and y˜ c (t) = yc (t)) and the proposed algorithm is still applicable.

˜ TM ˜ = M



(67)

and U is a non-singular transformation matrix such that the last ˜ U are linearly independents. (N + M − na − nb ) columns of W ˜ T and using (65) leads to the estimate Multiplying Eq. (58) by H of the cross-covariance vector Rˆ yf f : ˜ ˜ T Rˆ −1 Rˆ yf f = H ˜ ϕ ˜ ϕ

In this section, the fractional differentiation orders and coefficients are supposed unknown and the closed-loop parameters vector is defined by

f



+

−1 −1 = Rˆ ϕ H T Rˆ ϕ ˜ ˜ ˜ ϕ ˜ ϕ f

f

f

f



˜ˆ ols −  HT  ˜

To solve this problem, the use of a nonlinear optimization technique is required, such as the Gauss–Newton algorithm. The co-pfobels algorithm is summarized in two steps: • Step 1: Initialize 0 (i=0)

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

31

Periodogram Power Spectral Density Estimate 0

Power/frequency (dB/Hz)

−10

−20

−30

−40

−50

−60

0

1

2

3

4

5 6 Frequency (Hz)

7

8

9

10

Fig. 2. The power spectral density of the external excitation.

– The differentiation commensurate-order is initialized in order to obtain an initial estimation of the fractional closed-loop 0 ˜ˆ , the bias introduced by the ols method parameters vector  ols

0ˆ

ols

˜  ols

• Step 2: Gauss–Newton optimization – Initialize  which is a positive constant. – Refine the commensurate-order by Gauss–Newton algorithm as follows



applying

the

T (s) =

∂J ∂T =  ∂ ∂

(76)

∂T ∂ ∂ ∂

(77)

– Estimate the fractional closed-loop parameters vector using the ols method. when – Calculate the bias introduced by this method i+1 ˆ ˜  ols

 = i+1 and compute the fractional closed-loop parameters i+1

˜ˆ vector  fobels . – Evaluate the criterion J(i+1 ). ˆ  2

˜  ols

and i = i + 1

 i+1 i    ˆ − ˜ˆ      – Repeat until  olsi+1 ols  <  or a maximum number of iter ˜ˆ  ations is reached.

C (s) = 9 1 +

1 s



(79)

(75) i

where the gradient ∂J and the approximated hessian H are ∂ respectively given by:

– =

(78)

where  = 1.2. The fractional closed-loop transfer function is defined as follows



H=

b0 a2 s2 + a1 s + 1

where the coefficients a2 , a1 and b0 are unknown. The commensurate-order is supposed known and equals to  = 0.25. The fractional controller transfer function is given by

)

 i+1 i −1 ∂J  = −  H ∂

G (s) =

0

˜ˆ and the closed-loop parameters vector  fobels .

– Compute J(0ˆ

The fractional process to be identified is described by

ols

7. Numerical example In this section, a numerical example is presented to illustrate the performances of the pfobels method when np nb < na − np .

B1 s 1 + A1 s

˛

1

ˇ

1

+ A2 s

˛

2

+ A3 s

(80)

˛

3

where ˛1 = 1.2; ˛2 = 1.45; ˛3 = 1.7 and ˇ1 = 1.2. Since na = 2 and nb = 0, the prefilter order equals to nf =



na nb +1

− np



= 1.

7.1. Influence of the prefilter The external excitation yc (t) is a Pseudo Random Binary Sequence (PRBS) where its power spectral density is presented in Fig. 2. The output signal is contaminated by an additive zero-mean white noise e(t) with a signal-to noise ratio (SNR) of 10dB. The input and the output signals are uniformly sampled with a sampling period h = 0.05sec. The number of samples is Ns = 6000. The state variable filter parameters are chosen as follows: Nf = 2 and ωf = 10 rad/sec. The prefilter cut-off frequency is arbitrarily chosen as



f =

5,

15,

25

rad/sec.

Thus, three prefilters are considered to illustrate the influence of the prefiltering on the quality of estimation: 

F1 (s) =

1 1+

s 5

(81)

32

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

Table 1 Parameters estimation of fractional closed-loop system obtained by ols and pfobels methods (SNR = 10 dB). ols f

f = 5 rad/sec

f = 15 rad/sec

f = 25 rad/sec

pfobels

True

˜ˆ 

(× 10

A1 = 1.222 A2 = −1.00 A3 = 0.889 B1 = 0.111 A1 = 1.222 A2 = −1.000 A3 = 0.889 B1 = 0.111 A1 = 1.222 A2 = −1.000 A3 = 0.889 B1 = 0.111

−0.1247 0.0574 −0.087 0.0007 −0.354 0.198 −0.0305 −0.0007 −0.3136 0.1683 −0.0246 −0.0005

0.18 0.10 0.01 0.008 0.15 0.08 0.02 0.001 0.12 0.07 0.02 0.0013

−2

)

NRQE

1.6408

1.827

1.7922

˜ˆ 

(× 10−2 )

1.2243 −1.0015 0.8902 0.1113 1.223 −1.0004 0.8895 0.1112 1.2215 −1.007 0.889 0.111

3.44 1.89 1.51 0.38 3.7 2.37 2.25 0.34 3.13 2.28 2.20 0.35

NRQE

0.0230

0.0243

0.0246

Table 2 Parameters estimation of fractional process obtained by ols and pfobels methods (SNR = 10 dB). ols f f = 5 rad/sec

f = 15 rad/sec

f = 25 rad/sec



F2 (s) =

True

 × 10

a1 = −4.5 a2 = 4 b0 = 0.5 a1 = −4.5 a2 = 4 b0 = 0.5 a1 = −4.5 a2 = 4 b0 = 0.5

0.2601 −0.0351 −0.1710 0.8944 −0.1384 −0.2848 0.7575 −0.1109 −0.2646

0.47 0.07 0.08 0.37 0.06 0.07 0.21 0.04 0.04

1 1+

pfobels

˜

ˆ

−2

(82)

s 15

NRQE 1.0338

1.136

1.1119

˜

ˆ

 × 10−2

−4.5069 4.0058 0.5011 −4.5096 4.0098 0.5014 −4.5030 4.0004 0.4997

8.52 6.75 1.72 10.18 9.83 1.53 10.26 9.89 1.56

NRQE 0.0233

0.0236

0.0237

7.2. Choice of the state variable filter parameters

(83)

In this section, the choice of the state variable filter parameters is studied. The same external excitation yc (t), where its spectral density is presented in Fig. 2, is applied to the fractional closed-loop system (80). The output y(t) is corrupted by an additive gaussian white noise e(t) (SNR=10 dB).

The performances of the proposed algorithms are assessed with the help of a Monte Carlo simulation of Nmc = 500 runs with white noise realizations (SNR=10 dB). For each Monte Carlo run, the fractional closed-loop system is estimated using respectively the ols and the pfobels algorithms. The mean of the closed-loop system parameters and the process parameters obtained by the ols and the pfobels, their standard deviation and the normalized relative quadratic error NRQE defined by

7.2.1. Choice of the state variable filter cut-off frequency In this section, the influence of the state variable filter cut-off frequency on the quality $ %of estimation is studied. The state variable filter order is Nf = ˛N  + 1 = 2. A Monte carlo simulation of 500runs, with SNR=10 dB, is carried out for different values of wf ranged between 0.01 and 30 as



F3 (s) =

1 1+

s 25

!  2 " ˆ  " Nmc  −  " 1  #   NRQE =  Nmc

(84)

k=1

are summarized respectively in Tables 1 and 2. The identification results show that the ols algorithm gives a biased estimates: all the estimated parameters are far from the true ones. In contrast, the pfobels estimator gives a consistent estimates. Indeed, estimated coefficients converge to the real ones whatever the prefilter cut-off frequency. So, the prefilter cut-off frequency can be chosen equals to 5 rad/sec. The histograms of the process estimates are presented in Fig. 3. It shows that all of the obtained parameters by the pfobels method are around the true values. This validates the consistency of the proposed method and its robustness against the choice of the prefilter cut-off frequency f .

ωf =

&

0.01,

0.1,

5, 10,

20,

30

'

rad/sec.

(85)

The identification results are presented in Table 3. They show that the NRQE obtained by the ols method remains high whatever the value of ωf . So, the ols method fails to give a good estimation of the parameters. However, Table 3 shows that the value of the NRQE remains small by applying the pfobels algorithm even if it becomes larger when wf moved away from the cut-off frequency of the fractional closed-loop system (ωf = 10 rad/sec).

Table 3 Variation of the NRQE according to the state variable filter cut-off frequency ωf (SNR = 10 dB). NRQE ωf (rad/sec)

0.01

0.1

5

10

20

30

ols pfobels

3.1274 1.0330

3.2105 0.0538

1.8252 0.0423

1.7813 0.0210

1.7075 0.0824

1.6785 0.1182

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

7.2.2. Choice of the state variable filter order The effect of the filter order Nf is investigated in this section. Nf is chosen to range between 1 and 4. The filter cut-off frequency is set to ωf = 10 rad/sec and the prefilter cut-off frequency is fixed to f = 5 rad/sec. A Monte carlo simulation of 500 runs with SNR=10 dB for different values of Nf is realized. Both the ols and the pfobels methods are applied. The obtained results are presented in Table 4. The simulation results show that the smallest value of the NRQE is obtained for Nf = 2. 7.2.3. Choice of data samples number In this section, the filter order and the cut-off frequency are fixed respectively to Nf = 2 and ωf = 10 rad/sec. The prefilter cut-off

(a)

33

Table 4 Variation of the NRQE according to the filter order Nf (SNR = 10 dB). NRQE Nf

1

2

3

4

ols pfobels

1.6878 0.0976

1.7814 0.0193

1.8780 0.0351

1.9634 0.0684

frequency is fixed to f = 5 rad/sec and the number of samples is chosen as Ns =

&

'

200, 1000, 3000, 6000

The identification results obtained for different values of Ns , by applying both ols and pfobels algorithms, are given in Table 5. A

(b)

(c) Fig. 3. The histograms of the process estimates with the pfobels method ( f = 5, 15, 25 rad/sec).

34

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

Table 5 Variation of the NRQE according to the data samples number Ns (SNR = 10 dB).

which validate the efficiency of the proposed algorithm. Nevertheless, for SNR=3 dB, the estimates are biased with a larger standard deviation. but this can be ameliorated with more data samples.

NRQE Ns

200

1000

3000

6000

ols pfobels

1.887 0.1465

1.785 0.1095

1.711 0.0625

1.641 0.0230

7.4. Fractional order estimation higher identification accuracy is obtained when Ns = 6000 which is expected because more information can be guaranteed with a large number of data samples. 7.3. The effect of noise level In this section, the influence of the SNR on the quality of estimation is studied. The signal excitation is the same used previously and the output signal is corrupted by an additive gaussian noise with SNR=10 dB and SNR=3 dB. The state variable filter cut-off frequency and the filter order are chosen respectively as: ωf = 10 rad/sec and Nf = 2. The prefilter cut-off frequency f = 5 rad/sec. The Monte Carlo simulation results are summarized in Tables 6 and 7. These tables show that the ols method has an important bias. and the pfobels estimator gives a consistent estimates. The histograms of the pfobels estimates are plotted in Fig. 4. They show that the pfobels estimator gives a satisfactory results

In this section, the process commensurate-order  is supposed unknown and estimated along with the coefficients (a1 , a2 and b0 ). The external excitation is the same signal used previously and the output signal y(t) is corrupted by an additive noise e(t) with SNR=5 dB. The state variable filter cut-off frequency is ωf = 10 rad/sec and the filter order is set to Nf = 2. The prefilter cut-off frequency is fixed to f = 5 rad/sec and the initial commensurateorder is fixed to 0 = 0.1. For each Monte Carlo run, coefficients and commensurate-order of the fractional closed-loop transfer function are estimated using the co-pfobels algorithm and the process parameters are determined using Eq. (71). The closed-loop and the process estimates are respectively presented in Tables 8 and 9. The histograms are plotted in Fig. 5. The identification results show that the estimated parameters obtained by the co-pfobels converge to the true values. Thus, in such a high-level noisy context, the co-pfobels algorithm is asymptotically unbiased.

Table 6 Parameters estimation of fractional closed-loop system obtained by ols and pfobels methods (SNR = 10 dB and SNR = 3 dB). ols SNR

10dB

3dB

pfobels

True

ˆ 

(× 10

A1 = 1.222 A2 = −1.00 A3 = 0.889 B1 = 0.111 A1 = 1.22 A2 = −1.00 A3 = 0.889 B1 = 0.111

−0.1247 0.0574 −0.087 0.0007 −0.2084 0.1137 −0.0184 −0.0012

0.18 0.10 0.01 0.008 0.22 0.08 0.02 0.001

(a)

−2

)

NRQE

1.6408

1.827

ˆ 

(× 10−2 )

1.2243 −1.0015 0.8902 0.1113 1.2267 −1.0044 0.8929 0.1116

3.44 1.89 1.51 0.38 6.38 5.40 5.01 0.712

(b) Fig. 4. The histograms of the process estimates obtained with the pfobels method.

NRQE

0.0230

0.0835

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

35

Table 7 Parameters estimation of fractional process obtained by the ols and pfobels methods (SNR = 10 dB and SNR = 3 dB). ols SNR SNR=10 dB

SNR=3 dB

pfobels

True

ˆ

(× 10−2 )

a1 = −4.5 a2 = 4 b0 = 0.5 a1 = −4.5 a2 = 4 b0 = 0.5

0.2601 −0.0351 −0.1710 0.5114 −0.0827 −0.2125

0.47 0.07 0.08 0.24 0.07 0.06

a1

NRQE 1.0338

1.145

a2 170

170

100

100

100

100

50

50

50

50

−3,5

0

3

4

−4.5069 4.0058 0.5011 −4.5196 4.018 0.5022

8.52 6.75 1.72 11.43 12.24 3.91

5

0 0.4

0.5

NRQE 0.0233

0.0840

nu

170

−4.5

(× 10−2 )

a3

170

0 −5.5

ˆ

0.6

0 0.24

0.25

0.26

Fig. 5. The histograms of the process estimates obtained with the co-pfobels method (SNR=5 dB).

Table 8 Parameters estimation of fractional closed-loop system by co-pfobels method (SNR = 5 dB). co-pfobels True

ˆ 

(× 10−2 )

A1 = 1.222 A2 = −1.00 A3 = 0.889 B1 = 0.111

1.224 −1.0209 0.8901 0.1114

6.28 7.12 6.61 0.70

 = 0.25

NRQE

0.0457

0.2491

Table 9 Parameters estimation of fractional process by the co-pfobels method (SNR = 5 dB).

method early developed for rational systems. In the case of a restriction in the controller order, a stable prefilter with an appropriate order is designed to process input data. When, The differentiation orders and the coefficients are considered unknown, a nonlinear optimization algorithm is combined to the fobels method in order to estimate both of them. The efficiency of the proposed estimator has been evaluated through a numerical example. The Monte Carlo simulations have shown that the estimator gives a good results and provides unbiased estimates. Future developments will investigate the relation between fobels and the instrumental variable estimator. Another interesting perspective is the extension of the proposed method for error in variables context.

co-pfobels True

ˆ

(× 10−2 )

NRQE

a1 = −4.5 a2 = 4 b0 = 0.5

−4.5630 4.0740 0.5011

11.18 10.02 0.314

0.0472

 = 0.25

0.2491

8. Conclusion In this paper, the fractional closed-loop system identification context is considered where an indirect approach is used to estimate system parameters using the input and the output signals. The developed fobels method is an extension of the well known bels

References [1] A. Ben Hmed, M. Amairi, M. Aoun, Stability and resonance conditions of the noncommensurate elementary fractional transfer functions of the second kind, Commun. Nonlinear Sci. Numer. Simul. 22 (1) (2015) 842–865. [2] B. Saidi, M. Amairi, S. Najar, M. Aoun, Bode shaping-based design methods of a fractional order PID controller for uncertain systems, Nonlinear Dyn. 80 (4) (2015) 1817–1838, http://dx.doi.org/10.1007/s11071-014-1698-1 [3] A. Aribi, C. Farges, M. Aoun, P. Melchior, S. Najar, M.N. Abdelkrim, Fault detection based on fractional order models: Application to diagnosis of thermal systems, Commun. Nonlinear Sci. Numer. Simul. 19 (10) (2014) 3679–3693. [4] M. Chetoui, M. Thomassin, R. Malti, M. Aoun, S. Najar, M.N. Abdelkrim, A. Oustaloup, New consistent methods for order and coefficient estimation of continuous-time errors-in-variables fractional models, Comput. Math. Appl. 66 (5) (2013) 860–872.

36

Z. Yakoub et al. / Journal of Process Control 33 (2015) 25–36

[5] M. Amairi, M. Aoun, S. Najar, M.N. Abdelkrim, Guaranteed frequency-domain identification of fractional order systems, Int. J. Model. Identif. Control 17 (1) (2012) 32–42. [6] O. Cois, A. Oustaloup, T. Poinot, J. Battaglia, Fractional state variable filter for system identification by fractional model, in: European Control Conference (ECC), 2001, pp. 2481–2486. [7] M. Aoun, R. Malti, F. Levron, A. Oustaloup, Synthesis of fractional laguerre basis for system approximation, Automatica 43 (9) (2007) 1640–1648. [8] R. Malti, S. Victor, A. Oustaloup, H. Garnier, et al., An optimal instrumental variable method for continuous-time fractional model identification, in: 17th IFAC World Congress, 2008, pp. 1–6. [9] M. Amairi, Recursive set-membership parameter estimation using fractional model, Circuits Sys. Signal Process. (CSSP) (2015) 1–28, http://dx.doi.org/10. 1007/s00034-015-0036-2, Springer. [10] M. Amairi, M. Aoun, S. Najar, M.N. Abdelkrim, Set membership parameter estimation of linear fractional systems using parallelotopes, in: 9th International Multi-Conference on Systems, Signals and Devices (SSD), IEEE, 2012, pp. 1–6. [11] U. Forssell, L. Ljung, Closed-loop identification revisited, Automatica 35 (7) (1999) 1215–1241. [12] I. Gustavsson, L. Ljung, T. Söderström, Identification of processes in closed loop identifiability and accuracy aspects, Automatica 13 (1) (1977) 59–75. [13] M. Gilson, P. Van den Hof, Instrumental variable methods for closed-loop system identification, Automatica 41 (2) (2005) 241–249. [14] M. Tavakoli-Kakhki, M.S. Tavazoei, Proportional stabilization and closed-loop identification of an unstable fractional order process, J. Process Control 24 (5) (2014) 542–549.

[15] Z. Yakoub, M. Chetoui, M. Amairi, M. Aoun, Indirect approach for closed-loop system identification with fractional models, in: International Conference on Fractional Differentiation and its Applications (ICFDA), IEEE, 2014, pp. 1–6. [16] W.X. Zheng, C.-B. Feng, A bias-correction method for indirect identification of closed-loop systems, Automatica 31 (7) (1995) 1019–1024. [17] Z. Yakoub, M. Amairi, M. Chetoui, M. Aoun, A bias-eliminated least squares method for continuous-time fractional closed-loop system identification, in: IEEE Conference on Control Applications (CCA), IEEE, 2014, pp. 128–133. [18] W.X. Zheng, Identification of closed-loop systems with low-order controllers, Automatica 32 (12) (1996) 1753–1757. [19] H. Garnier, M. Gilson, W. Zheng, A bias-eliminated least-squares method for continuous-time model identification of closed-loop systems, Int. J. Control 73 (1) (2000) 38–48. [20] A. Oustaloup, La dérivation non entière, Hermès (1995). [21] M. Aoun, Systemes linéaires non entiers et identification par bases orthogonales non entieres, Ph.D. thesis, Bordeaux 1, 2005. [22] S. Victor, R. Malti, Model order identification for fractional models, in: European Control Conference (ECC), IEEE, 2013, pp. 3470–3475. [23] K. Oldham, J. Spanier, The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order, 1974. [24] Z. Yakoub, M. Amairi, M. Chetoui, M. Aoun, On the closed-loop system identification with fractional models, Circuits Sys. Signal Process. (CSSP) (2015) 1–28, http://dx.doi.org/10.1007/s00034-015-0046-0, Springer. [25] S. Victor, R. Malti, H. Garnier, A. Oustaloup, et al., Parameter and differentiation order estimation in fractional models, Automatica 49 (4) (2013) 926–935.