# A Bias-Corrected Method for Fractional Linear Parameter Varying Systems

Authors: Zaineb Yakoub; Omar Naifar; Messaoud Amairi; Manel Chetoui; Mohamed Aoun; Abdellatif Ben Makhlouf
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This paper proposes an identification algorithm for the fractional Linear Parameter Varying (LPV) system considering noisy scheduling and output measurements. A bias correction technique is provided in order to compensate for the bias caused by the least squares algorithm. This approach was created to estimate either coefficients or fractional-order differentiation, and it has been proven to produce unbiased and reliable results. The suggested method’s performance is assessed by the identification of two fractional models and was compared with Nelder–Mead Simplex method.


---

## Body

## 1. Introduction

Fractional-order systems may be utilized to more effectively represent the dynamics of many real-world systems due to the memory characteristic and nonlocality of the fractional-order calculus [1–4]. As a result, the study of fractional-order systems has become a hotspot for research (see [5–7] and the references therein), and it has evolved in a variety of fields, including system identification [2, 4, 8, 9], stability [10, 11], robust stability [12, 13], control [14–16], and diagnosis [17–19]. For more than two decades, fractional-order system identification has been a major research issue. As a result, various approaches in the literature have been created [20]. This paragraph provides a quick overview of fractional system identification methods. Le Lay [21] developed the first study that dealt with fractional systems identification in 1998. The identification of open-loop systems with fractional models is established in this article. Battagalia et al. [6] established a fractional model which produced the transient thermal behavior of a system in 2000. Nearby this year, authors in [22] have discussed new ways to model a system with the fractional model. In this work, the estimation of the fractional dynamics of a lead-acid battery has been investigated. Subsequently, in 2001, the use of state variable filters (svf) was first extended by Cois et al. to fractional systems, with an instrumental variable (IV) procedure [23]. Later, in [24], the fractional orthogonal bases Laguerre functions have been synthesized. In 2006, the identification of lead-acid battery state of charge estimation was developed by Sabatier et al. in [1]. In 2008, Malti et al. have developed the fractional optimal IV method for model identification [25]. Then, in 2011, fractional thermal systems identification was illustrated by Gabano et al. [26]. In [27], the genetic algorithm is used for the continuous-time identification. More significantly, in 2015, and prompted by the reality that many technologies can only run in fractional closed-loop configurations due to stability and sensitivity issues, many systems can only work in fractional closed-loop configurations. The authors in [28] have proposed a bias compensated method for system identification of the fractional closed-loop system. In 2018, the nonlinear system identification of fractional Wiener models has been established by Sersour et al. [29]. Very recently, in 2020-2021, the fractional order system identification has attracted more and more attention: In [30] the authors have proposed an output error-based method to identify multi input single output (MISO) systems using fractional models. In paper [31], the authors consider the problem of identifying linear time invariant (LTI) fractional systems from noisy input and output measurements using the bias correction scheme. Moreover, in [32, 33], fractional derivatives based on the epidemic system have also been used to deal with some epidemic behaviors like the CORONA virus. In conclusion, the expansion of fractional calculus applications and other mathematical phenomena more closely reflects a real system than traditional integer models.According to the abovementioned studies, despite the integer LPV modeling field has emerged rapidly over the recent years and it was used in many control applications like aircrafts [34–36], only few literatures have been presented on fractional-order LPV systems. Therefore, the fractional LPV system identification was developed firstly in [37]. The authors developed the simplified refined instrumental variable continuous-time (SRIVC) method to identify fractional LPV systems with output error (OE) technique and a nonlinear programming. More recently, the fractional closed-loop LPV system identification has been investigated in [38].According to the abovementioned studies, the main contribution of this paper is to establish a parameter identification method for continuous-time fractional-order LPV systems under the Grünwald–Letnikov definition. This method is based on the bias compensated least squares algorithm known by its simplicity of concept and its power in the fields of identification and control theory. Firstly, the parameters vector estimation is obtained by the fractional-order least squares algorithm combined with the state variable filter approach. Then, to achieve the estimation consistency, the bias produced by the least squares is eliminated. Moreover, this algorithm combines with a nonlinear algorithm optimization to identify both fractional-order differentiations and linear coefficients in this work. In addition, this contribution seems interesting to deal with the more realistic case where the scheduling variables are affected by a measurement noise.The paper has the following structure: the fractional calculus and systems description is given in section2. In section 3, the identification problem is formulated. Section 4 establishes the creation of the fractional-order bias removed least squares technique for estimating both fractional-order and LPV coefficients. Section 5 contains two numerical examples as well as comparison simulations. In section 6, the closing comments are presented.


## 2. Preliminaries

In this section, we introduce the main notations, definitions, and theorems that will be used in the rest of this paper.


### 2.1. Notation

This section presents the used notations in this paper. Letℝ denote the set of all real numbers. ℝ+∗ be the set of all strictly positive real numbers. ℝn is defined as the set of real vectors of dimensions n. ℝn×m is the set of n×m matrices. The inverse and the transpose of the matrix X are defined as follows X−1 and XT, respectively. x2=xTx designates the squares of the 2-norm of x. D=p=d/dt is the differential operator. The Laplace operator is symbolized by s. The sampled period is denoted by T, and Nd defines the number of data.


### 2.2. Fractional-Order Differentiation

Unlike integer calculus, there are several definitions of fractional derivative. The known three kinds of definitions for fractional calculus are Riemann–Liouville (R-L), Grünwald–Letnikov (G-L), and Caputo definitions [39]. In this paper, the Grünwald–Letnikov definition is used.Definition 1 (see [40]).
The Grünwald–Letnikov fractional derivative approximation of the orderα of function χt is described as(1)χGLD0,tαt≃1Tα∑k=0t/T−1kαkχt−kT,where αk is the generalized form of the Newton’s binomial to non-integer orders.


### 2.3. Fractional LPV System Description

Consider the single-input-single-output (SISO) fractional-order LPV differential equation using the definition of Grünwald–Letnikov defined as follows:(2)∑n=0NanρtDαnyt=∑m=0MbmρtDβmut,where the input and the output signals are presented by ut and yt, respectively. The scheduling measurement is given by ρt. The linear coefficients an,bm∈ℝ0<n≤N;0≤m≤M and the non-integer orders are αn,βm∈ℝ+ and they verify the following equation:(3)α1<α2<⋅⋅⋅<αN;β0<β1<⋅⋅⋅<βM.The commensurate order differential equation is given by(4)∑i=0naaiρtDiνyt=∑j=0nbbjρtDjνut,with ν designates the commensurate order and where na=αn/υ and nb=βm/υ are integers and verify na>nb.Assuming null initial conditions, the SISO fractional commensurate order LPV system representation in the Laplace domain is described by(5)Gs=Bs,ρtAs,ρt=∑j=0nbbjρtsjν∑i=0naaiρtsiν,where ρt is the scheduling variable.


## 2.1. Notation

This section presents the used notations in this paper. Letℝ denote the set of all real numbers. ℝ+∗ be the set of all strictly positive real numbers. ℝn is defined as the set of real vectors of dimensions n. ℝn×m is the set of n×m matrices. The inverse and the transpose of the matrix X are defined as follows X−1 and XT, respectively. x2=xTx designates the squares of the 2-norm of x. D=p=d/dt is the differential operator. The Laplace operator is symbolized by s. The sampled period is denoted by T, and Nd defines the number of data.


## 2.2. Fractional-Order Differentiation

Unlike integer calculus, there are several definitions of fractional derivative. The known three kinds of definitions for fractional calculus are Riemann–Liouville (R-L), Grünwald–Letnikov (G-L), and Caputo definitions [39]. In this paper, the Grünwald–Letnikov definition is used.Definition 1 (see [40]).
The Grünwald–Letnikov fractional derivative approximation of the orderα of function χt is described as(1)χGLD0,tαt≃1Tα∑k=0t/T−1kαkχt−kT,where αk is the generalized form of the Newton’s binomial to non-integer orders.


## 2.3. Fractional LPV System Description

Consider the single-input-single-output (SISO) fractional-order LPV differential equation using the definition of Grünwald–Letnikov defined as follows:(2)∑n=0NanρtDαnyt=∑m=0MbmρtDβmut,where the input and the output signals are presented by ut and yt, respectively. The scheduling measurement is given by ρt. The linear coefficients an,bm∈ℝ0<n≤N;0≤m≤M and the non-integer orders are αn,βm∈ℝ+ and they verify the following equation:(3)α1<α2<⋅⋅⋅<αN;β0<β1<⋅⋅⋅<βM.The commensurate order differential equation is given by(4)∑i=0naaiρtDiνyt=∑j=0nbbjρtDjνut,with ν designates the commensurate order and where na=αn/υ and nb=βm/υ are integers and verify na>nb.Assuming null initial conditions, the SISO fractional commensurate order LPV system representation in the Laplace domain is described by(5)Gs=Bs,ρtAs,ρt=∑j=0nbbjρtsjν∑i=0naaiρtsiν,where ρt is the scheduling variable.


## 3. Problem Formulation

### 3.1. Data Generating Fractional-Order LPV System

In this subsection, a brief description of Continuous-Time (CT) fractional order LPV systemG0 is presented. The generic form of the considered model is depicted in Figure 1 and governing mathematical relation is given by(6)G0:A0p,ρ0tχ0t=ℬ0p,ρ0tut,yt=χ0t+e0t.Figure 1
The fractional-order LPV model.Hence,p=d/dt denotes the differentiation operator. ut denotes the input signal of the plant. Let χ0t designates the noise-free output. yt is the available observation of the output signal that is affected by a noise sequence e0t which assumed to be a white noise.Remark 1.
The scheduling signal measurements are also contaminated by an additive measurement noisev0t(7)ρ0t=ρ0t+v0t,where ρ0t is the noise-free scheduling signal and v0t is an additive white noise with finite variance σv2 and independent of the output noise e0t.
The vector of scheduling parametersP0t∈ℝnρ+1 is expressed by(8)P0t=1ρ0tρ02t…ρ0nρtT.
A0p,ρ0t and ℬ0p,ρ0t are polynomials in p of degree na,nb expressed as follows:(9)A0p,ρ0t=1+∑i=1naai0ρ0tpiν;ℬ0p,ρ0t=∑j=0nbbj0ρ0tpjν,where ν is the commensurate order of the process.
The coefficientsai0 and bj0 are meromorphic and polynomial functions in ρ0t of maximum degree nρ and defined as follows:(10)ai0ρ0t=ai,00+∑l=1nρai,l0ρ0lt;bj0ρ0t=bj,00+∑l=1nρbj,l0ρ0lt,where ai0 and bj0 are the unknown parameters.


### 3.2. System Model: Fractional-Order LPV System

Consider the parameterized model presented by the following equation:(11)ℳ:yt=−∑i=1naaiρtDiνyt+∑j=0nbbjρtDjνut+εt.where εt is the residual term. The functions ai and bj which depend on the noisy observation ρt of the scheduling signal are polynomials parameterized as(12)aiρt=ai,0+∑l=1nρai,lρlt;bjρt=bj,0+∑l=1nρbj,lρlt.For convenience, the following notations are introduced in this subsection.The regression vectorφTt is expressed by(13)φTt=−Dνyt,…,−Dnaνyt,ut,…,Dnbνut.The findings of direct fractional derivative of noisy output are not precise. It is suggested that the state variable filter (svf) be used to offer fractional-order differentiation of both input and output signals to overcome this problem [41].The svf is presented by this equation as(14)Fαp=pαλFλF+pκ,where α∈ℝ+∗ is the order differentiation and λF is the cut-off frequency of the svf and κ∈ℕ selected so that κ>naν.We applyFαp to both sides of (11) and then(15)Ap,ρtyft=ℬp,ρtuft+εft.The filtered input/output signalsuft and yft are written as(16)Djνuft=Fjνput;Diνyft=Fiνpyt.Considereft represents the combined effects of the output observations noise et and the error introduced by the operation filtering. Hence, the error equation εft may be rewritten as(17)εft=Ap,ρteft.The filtered regression vector is expressed as follows:(18)φfTt=−Dνyft,…,−Dnaνyft,uft,…,Dnbνuft.The following linear regression form can be used to show the fractional system’s input and output behavior:(19)yft=ΦfTtθ+εft,where(20)Φft=φfTt⊗ρt,where ⊗ denotes the kronecker product.The parameters vectorθ is defined as(21)θ=a1T,…,anaT,b0T,…,bnbT∈ℝnρ+1na+nb+1,where(22)bj=bj,1bj,2…bj,nρT,ai=ai,1ai,2…ai,nρT.Certainly, the novelty made in this study is the use of the bias compensated method to identify the fractional LPV system fromNd data points of the input and noisy output signals.


## 3.1. Data Generating Fractional-Order LPV System

In this subsection, a brief description of Continuous-Time (CT) fractional order LPV systemG0 is presented. The generic form of the considered model is depicted in Figure 1 and governing mathematical relation is given by(6)G0:A0p,ρ0tχ0t=ℬ0p,ρ0tut,yt=χ0t+e0t.Figure 1
The fractional-order LPV model.Hence,p=d/dt denotes the differentiation operator. ut denotes the input signal of the plant. Let χ0t designates the noise-free output. yt is the available observation of the output signal that is affected by a noise sequence e0t which assumed to be a white noise.Remark 1.
The scheduling signal measurements are also contaminated by an additive measurement noisev0t(7)ρ0t=ρ0t+v0t,where ρ0t is the noise-free scheduling signal and v0t is an additive white noise with finite variance σv2 and independent of the output noise e0t.
The vector of scheduling parametersP0t∈ℝnρ+1 is expressed by(8)P0t=1ρ0tρ02t…ρ0nρtT.
A0p,ρ0t and ℬ0p,ρ0t are polynomials in p of degree na,nb expressed as follows:(9)A0p,ρ0t=1+∑i=1naai0ρ0tpiν;ℬ0p,ρ0t=∑j=0nbbj0ρ0tpjν,where ν is the commensurate order of the process.
The coefficientsai0 and bj0 are meromorphic and polynomial functions in ρ0t of maximum degree nρ and defined as follows:(10)ai0ρ0t=ai,00+∑l=1nρai,l0ρ0lt;bj0ρ0t=bj,00+∑l=1nρbj,l0ρ0lt,where ai0 and bj0 are the unknown parameters.


## 3.2. System Model: Fractional-Order LPV System

Consider the parameterized model presented by the following equation:(11)ℳ:yt=−∑i=1naaiρtDiνyt+∑j=0nbbjρtDjνut+εt.where εt is the residual term. The functions ai and bj which depend on the noisy observation ρt of the scheduling signal are polynomials parameterized as(12)aiρt=ai,0+∑l=1nρai,lρlt;bjρt=bj,0+∑l=1nρbj,lρlt.For convenience, the following notations are introduced in this subsection.The regression vectorφTt is expressed by(13)φTt=−Dνyt,…,−Dnaνyt,ut,…,Dnbνut.The findings of direct fractional derivative of noisy output are not precise. It is suggested that the state variable filter (svf) be used to offer fractional-order differentiation of both input and output signals to overcome this problem [41].The svf is presented by this equation as(14)Fαp=pαλFλF+pκ,where α∈ℝ+∗ is the order differentiation and λF is the cut-off frequency of the svf and κ∈ℕ selected so that κ>naν.We applyFαp to both sides of (11) and then(15)Ap,ρtyft=ℬp,ρtuft+εft.The filtered input/output signalsuft and yft are written as(16)Djνuft=Fjνput;Diνyft=Fiνpyt.Considereft represents the combined effects of the output observations noise et and the error introduced by the operation filtering. Hence, the error equation εft may be rewritten as(17)εft=Ap,ρteft.The filtered regression vector is expressed as follows:(18)φfTt=−Dνyft,…,−Dnaνyft,uft,…,Dnbνuft.The following linear regression form can be used to show the fractional system’s input and output behavior:(19)yft=ΦfTtθ+εft,where(20)Φft=φfTt⊗ρt,where ⊗ denotes the kronecker product.The parameters vectorθ is defined as(21)θ=a1T,…,anaT,b0T,…,bnbT∈ℝnρ+1na+nb+1,where(22)bj=bj,1bj,2…bj,nρT,ai=ai,1ai,2…ai,nρT.Certainly, the novelty made in this study is the use of the bias compensated method to identify the fractional LPV system fromNd data points of the input and noisy output signals.


## 4. Fractional LPV System Identification

The fractional parameters vector and fractional-order differentiation are both identified in this section. To begin, the best coefficient is calculated using the fractional-order bias compensated least squaresfobcls technique for initial values of the fractional-order differentiation. The fractional-order differentiation is then computed by combining the fobcls technique with a nonlinear algorithm.


### 4.1. Fractional-Order Least Squaresfols Estimates

Firstly, the fractional commensurate orderν is meant to be known a priori, and the goal of this work is to identify only the coefficients.In this section, in order to identify the continuous-time fractional order LPV systems, the fractional-order least squaresfols method is detailed.Letθ^folsNd, the estimation of θ, be found by minimizing the following criterion:(23)CNdθ^folsNd=1Nd∑t=1Ndεf2t,where εft denotes the filtered equation error, which is provided by the equation as follows:(24)εft=yft−ΦfTtθ.As a result, the least squares estimatedθ is defined in the following way:(25)θ^folsNd=EΦfΦfT−1EΦfyf,where EΦfΦfT is the autocovariance matrix defined by(26)EΦfΦfT=∑t=0NdΦftΦfTt.AndEΦfyf is the cross-covariance vector given by(27)EΦfyf=∑t=0NdΦftyft.Proposition 1.
Under the assumption that the matrixEΦfΦfT is invertible, the difference between the fols estimates and the true parameters vector is given by(28)limNd⟶∞θ^folsNd−θ=Δ^θ,where Δ^θ is the bias yields by the fols method and given by the following expression:(29)Δ^θ=EΦfΦfT−1EΦfεf,=EΦfΦfT−1ΛEYfεf,where the matrix Λ=Ina×nρ+10nb+1×nρ∈ℝna+nb+1×nρ+1×na×nρ+1 and Yft=yft⊗ρt.Proof.
Substituting the filtered output signalyft with equations (18) in (24) produced the following equation:(30)θ^folsNd=∑t=0NdΦftΦfTt−1∑t=0NdΦftΦfTtθ+εft=θ+∑t=0NdΦftΦfTt−1∑t=0NdΦftεft=θ+EΦfΦfT−1EΦfεf.
Substituting equation (20) in the above equation yields(31)θ^folsNd=θ+EΦfΦfT−1Eφft⊗ρtεft,where(32)Eφft⊗ρtεft=E−Dνyft−Dνyftρt…−Dνyftρnρt⋮⋮⋮⋮−Dnaνyft−Dnaνyftρt…−Dnaνyftρnρtuftuftρt…uftρnρt⋮⋮⋮⋮DnbνuftDnbνuftρt…Dnbνuftρnρtεft.
Sinceut and vt are uncorrelated with the noise et, equation (32) can be written as follows:(33)Eφft⊗ρtεft=E−Dνyft⊗ρtεft⋮−Dnaνyft⊗ρtεft0⋮0.
Replacing equations (33) in (30), we get the following equation:(34)limNd⟶∞θ^folsNd=θ+EΦfΦfT−1ΛEyft⊗ρtεft=θ+Δ^θ.
Consequently, subtractingθ from limNd⟶∞θ^folsNd in equation (34) makes us obtain (29).
Equation (34) shows that the fols method fails to produce unbiased parameters and correct results because et and Φt are mutually correlated. As a result, the fractional-order bias compensated least squares fobcls approach is detailed in the next section to address the bias problem.


### 4.2. Fractional-Order Bias Compensated Least Squaresfobcls Algorithm

The bias compensated least squaresbcls algorithm was introduced firstly by Zheng [42]. In this study, the aim is to formulate an update version of the bcls algorithm for parameter estimation of a fractional LPV system from noisy scheduling and output observations.Thus, thefobcls estimator is expressed as(35)θ^fobclsNd=θ^folsNd−Δ^θ=θ^folsNd−EΦfΦfT−1ΛEYfεf.Thus, the assessment of the bias is important to the actual use of the unbiased identification method. Thefols approach, on the other hand, can give a single formula for computing the estimations. Therefore, further equations connecting the noise to other known variables are necessary. In order to achieve this aim, an enhanced system was used to raise the nominator parameters by na−nb dimensions which is necessary, where the introduced parameter bj=0j=nb+1,…,na+nb. Therefore, the augmented transfer function is given by(36)G¯p,ρt=B¯p,ρtAp,ρt=b0ρt+b1ρtpν+⋯+bna+nbρtpna+nbν1+a1ρtpν+⋯+anaρtpnaν,and the augmented numerator B¯p is written as(37)B¯p=b0+b1pν+⋯+bnbpnbν+⋯+bna+nbpna+nbν.θ¯ is presented by(38)θ¯=θb¯Tρt,where(39)b¯T=bnb+1T…bnb+naT=0…0∈ℝnanρ+1.Similarly to (18), the filtered regression form is described by(40)y¯ft=Φ¯fTtθ¯+εft,where(41)Φ¯ft=φ¯fTt⊗ρt,where φ¯fTt is expressed as follows:(42)φ¯fTt=−Dνyft…−Dnaνyftuft…Dna+nbνuft,=φfTtu¯fTt,where(43)u¯fTt=Dnb+1νuft…Dna+nbνuft.The estimatorθ¯^folsNd using the fols method is expressed as(44)θ¯^folsNd=EΦ¯fΦ¯fT−1EΦ¯ftyft=θ¯+EΦ¯fΦ¯fT−1Eφ¯fTt⊗ρtεft,where(45)EΦ¯fΦ¯fT=∑t=0NdΦ¯ftΦ¯fTt.Sincevt, et, and ut are mutually independent, we obtain(46)θ^folsNd=θ¯+EΦ¯fΦ¯fT−1E−Dνyft⊗ρtεft⋮−Dnaνyft⊗ρtεft0⋮0,=θ¯+EΦ¯fΦ¯fT−1Λ¯Eyft⊗ρtεft,=θ¯+EΦ¯fΦ¯fT−1Λ¯EYfεf.The arrayΛ¯ is given by the following equation:(47)Λ¯=Ina×nρ+10na+nb+1×nρ∈ℝ2na+nb+1×nρ+1×na×nρ+1,and the matrix EΦ¯fΦ¯fT may be described by(48)EΦ¯fΦ¯fT=EΦfΦfEΦfu¯fEΦfu¯fEu¯fu¯f.Consequently, its inverseEΦ¯fΦ¯fT−1 is given by the following expression:(49)EΦ¯fΦ¯fT−1=EΦfΦf−1+X−EΦfΦf−1EΦfu¯fD−1−D−1EΦfu¯fTEΦfΦf−1D−1.where X=EΦfΦfT−1EΦfu¯fD−1EΦfu¯fTEΦfΦfT−1 and D=Eu¯fu¯f−EΦfu¯fTEΦfΦfT−1EΦfu¯fCombining (43) with (48),(50)b¯^=−D−1EΦfu¯fTEΦfΦfT−1EΦfyf+D−1Eu¯fyf=−D−1EΦfu¯fTθ^folsNd−Eu¯fyf.On the other hand, it follows from (37), (38), (45), and (49) and since the augmented polynomial B¯p has known parameters vector b¯=0, it is evident that(51)b¯^=b¯−D−1EΦfu¯fTEΦfΦfT−1ΛEYfεf,=−D−1EΦfu¯fTEΦfΦfT−1ΛEYfεf.Thus, (50) and (52) give the following equality:(52)−D−1EΦfu¯fTEΦfΦfT−1ΛEYfεf=−D−1EΦfu¯fTθ^folsNd−Eu¯fyf.From (52), it is clear that(53)EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1EΦfu¯fTθ^folsNd−Eu¯fyf.Replacing (53) in (28) leads to(54)Δ^θ=EΦfΦfT−1ΛEΦfu¯fTEΦfΦfT−1Λ−1EΦfu¯fTθ^folsNd−Eu¯fyf.Finally, we obtain(55)θ^fobclsNd=θ^folsNd−Δ^θ,=θ^folsNd−EΦfΦfT−1ΛEΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTθ^folsNd−Eu¯fyf.To facilitate the use of this method, a detailed description is given in Algorithm 1.


#### 4.2.1. Convergence Analysis

We begin this part by demonstrating that the fobcls method is highly reliable.Theorem 1.
The proposedfobclsestimator for fractional order system identification gives consistent estimates. The optimal augmented estimator θfobclsNd verifies the following equation:(56)limNd⟶∞θ^fobclsNd=θ.Proof.
Firstly, we recall that the filtered output signal is expressed as(57)yft=ΦfTθ+εft.
Multiplying both side of equation (57) by u¯ft, we obtain(58)u¯ftyft=u¯ftΦfTθ
Because the noiseset and vt are mutually uncorrelated of each other and of ut, we get(59)Eu¯fyf=EΦfu¯ftθ.
So, using equation (53),(60)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTθ^folsNd−Eu¯fyf.
ReplacingEu¯fyf in the above equation by its expression given in equation (57), we obtain(61)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTθ^folsNd−EΦfu¯ftθ.
Thus, the result is shown as the following equation:(62)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦf−1Λ−1×EΦfu¯fTθ^folsNd−θ.
Now, replacing the termθ^folsNd−θ by the expression of Δ^θ defined in equation (54),(63)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTEΦfΦfT−1ΛEYfεf.
Consequently,(64)limNd⟶∞EYfεf=EYfεf.
Using equation (65) and equation (25), we find(65)limNd⟶∞θ^fobcls=θ^folsNd−EΦfΦfT−1ΛEYfεf=θ+EΦfΦf−1ΛEYfεf−EΦfΦfT−1ΛEYfεf.
Finally, lettingNd⟶∞, we obtain(66)limNd⟶∞θ^fobclsNd=θ.Algorithm 1: Fractional-order bias compensated least squares algorithm.
Step 1: Generate the data using the svfuft,yft|t=1,…,Nd.Step 2: Build the filtered regressor using equation (13)Step 3: Apply thefols method to obtain θ^fols using equation (43)Step 4: Compute the matrixΛ and the vector u¯ft and build the augmented regressor vector φf−Τt using equations (46), (42), and (41)Step 5: Compute the bias expression via equation (53)Step 6: Calculateθ^fobclsNd using equation (54)


### 4.3. Fractional-Order Estimation

Thefobcls technique can only be used in the scenario where the fractional orders are considered to be known a priori. The main contribution of this part is to design an algorithm extending the identification method studied previously to a more normal situation in which the fractional-order differentiation is supposed to be unknown and the coefficients are approximated. Fractional-order optimization-bias compensated least squares foo-bcls is the name of this algorithm. It is based on a combination of the fobcls coefficient estimate technique and a nonlinear differentiation order optimization algorithm. As an option to tackling the optimization problem, the Gauss–Newton method is offered.For that, consider the parameter vectorθ as follows:(67)θ=a1,…,ana,b0,…,bnb,ν.The challenge of parameter identification is framed as a functional minimization. As a result, the primary purpose of this method is to reduce the error described by equation (69) with regard to ν. As a result, the Gauss–Newton method is used to repeatedly update the commensurate order. As a result, the quadratic criteria are as follows:(68)Wθ^foo−bcls=12εt2,where the error εt is expressed by(69)εt=yt−y^t.At each the iterationj+1, this algorithm calculates the fractional differentiation order νj+1(70)νj+1=νj−γℋ−1∂W∂ν|νj.The gradient∂W/∂ν is described as follows:(71)∂W∂ν=∂εTt∂νεt.and the approximated hessian ℋ is given by(72)ℋ=∂εTt∂ν∂εt∂ν.The description of thefoo-bcls method is given by Algorithm 2.Data {uf (t), yf(t)t = 1, . . . , Nd}(1)
InitializationInitializeν0 .The initial estimation of theθˆ is obtained by Applying Algorithm 1.Compute W (θˆfoo–bcls) (0)Initializeλ ∈ R∗ .(2)
repeatApply the Gauss–Newton algorithm in order to update the commensurate order using equation (71)The parameters vector is estimated by applying Algorithm 1.Calculate the quadratic criterionW(θ^foo−bclsj+1 using equation (69)−λ=λ/2(3)
untilΔ^θi+1−Δ^θi/Δ^θi+1<σorj=MaxIter where MaxIter designates the maximum number of iteration(4)
returnθˆfoo−bcls.


## 4.1. Fractional-Order Least Squaresfols Estimates

Firstly, the fractional commensurate orderν is meant to be known a priori, and the goal of this work is to identify only the coefficients.In this section, in order to identify the continuous-time fractional order LPV systems, the fractional-order least squaresfols method is detailed.Letθ^folsNd, the estimation of θ, be found by minimizing the following criterion:(23)CNdθ^folsNd=1Nd∑t=1Ndεf2t,where εft denotes the filtered equation error, which is provided by the equation as follows:(24)εft=yft−ΦfTtθ.As a result, the least squares estimatedθ is defined in the following way:(25)θ^folsNd=EΦfΦfT−1EΦfyf,where EΦfΦfT is the autocovariance matrix defined by(26)EΦfΦfT=∑t=0NdΦftΦfTt.AndEΦfyf is the cross-covariance vector given by(27)EΦfyf=∑t=0NdΦftyft.Proposition 1.
Under the assumption that the matrixEΦfΦfT is invertible, the difference between the fols estimates and the true parameters vector is given by(28)limNd⟶∞θ^folsNd−θ=Δ^θ,where Δ^θ is the bias yields by the fols method and given by the following expression:(29)Δ^θ=EΦfΦfT−1EΦfεf,=EΦfΦfT−1ΛEYfεf,where the matrix Λ=Ina×nρ+10nb+1×nρ∈ℝna+nb+1×nρ+1×na×nρ+1 and Yft=yft⊗ρt.Proof.
Substituting the filtered output signalyft with equations (18) in (24) produced the following equation:(30)θ^folsNd=∑t=0NdΦftΦfTt−1∑t=0NdΦftΦfTtθ+εft=θ+∑t=0NdΦftΦfTt−1∑t=0NdΦftεft=θ+EΦfΦfT−1EΦfεf.
Substituting equation (20) in the above equation yields(31)θ^folsNd=θ+EΦfΦfT−1Eφft⊗ρtεft,where(32)Eφft⊗ρtεft=E−Dνyft−Dνyftρt…−Dνyftρnρt⋮⋮⋮⋮−Dnaνyft−Dnaνyftρt…−Dnaνyftρnρtuftuftρt…uftρnρt⋮⋮⋮⋮DnbνuftDnbνuftρt…Dnbνuftρnρtεft.
Sinceut and vt are uncorrelated with the noise et, equation (32) can be written as follows:(33)Eφft⊗ρtεft=E−Dνyft⊗ρtεft⋮−Dnaνyft⊗ρtεft0⋮0.
Replacing equations (33) in (30), we get the following equation:(34)limNd⟶∞θ^folsNd=θ+EΦfΦfT−1ΛEyft⊗ρtεft=θ+Δ^θ.
Consequently, subtractingθ from limNd⟶∞θ^folsNd in equation (34) makes us obtain (29).
Equation (34) shows that the fols method fails to produce unbiased parameters and correct results because et and Φt are mutually correlated. As a result, the fractional-order bias compensated least squares fobcls approach is detailed in the next section to address the bias problem.


## 4.2. Fractional-Order Bias Compensated Least Squaresfobcls Algorithm

The bias compensated least squaresbcls algorithm was introduced firstly by Zheng [42]. In this study, the aim is to formulate an update version of the bcls algorithm for parameter estimation of a fractional LPV system from noisy scheduling and output observations.Thus, thefobcls estimator is expressed as(35)θ^fobclsNd=θ^folsNd−Δ^θ=θ^folsNd−EΦfΦfT−1ΛEYfεf.Thus, the assessment of the bias is important to the actual use of the unbiased identification method. Thefols approach, on the other hand, can give a single formula for computing the estimations. Therefore, further equations connecting the noise to other known variables are necessary. In order to achieve this aim, an enhanced system was used to raise the nominator parameters by na−nb dimensions which is necessary, where the introduced parameter bj=0j=nb+1,…,na+nb. Therefore, the augmented transfer function is given by(36)G¯p,ρt=B¯p,ρtAp,ρt=b0ρt+b1ρtpν+⋯+bna+nbρtpna+nbν1+a1ρtpν+⋯+anaρtpnaν,and the augmented numerator B¯p is written as(37)B¯p=b0+b1pν+⋯+bnbpnbν+⋯+bna+nbpna+nbν.θ¯ is presented by(38)θ¯=θb¯Tρt,where(39)b¯T=bnb+1T…bnb+naT=0…0∈ℝnanρ+1.Similarly to (18), the filtered regression form is described by(40)y¯ft=Φ¯fTtθ¯+εft,where(41)Φ¯ft=φ¯fTt⊗ρt,where φ¯fTt is expressed as follows:(42)φ¯fTt=−Dνyft…−Dnaνyftuft…Dna+nbνuft,=φfTtu¯fTt,where(43)u¯fTt=Dnb+1νuft…Dna+nbνuft.The estimatorθ¯^folsNd using the fols method is expressed as(44)θ¯^folsNd=EΦ¯fΦ¯fT−1EΦ¯ftyft=θ¯+EΦ¯fΦ¯fT−1Eφ¯fTt⊗ρtεft,where(45)EΦ¯fΦ¯fT=∑t=0NdΦ¯ftΦ¯fTt.Sincevt, et, and ut are mutually independent, we obtain(46)θ^folsNd=θ¯+EΦ¯fΦ¯fT−1E−Dνyft⊗ρtεft⋮−Dnaνyft⊗ρtεft0⋮0,=θ¯+EΦ¯fΦ¯fT−1Λ¯Eyft⊗ρtεft,=θ¯+EΦ¯fΦ¯fT−1Λ¯EYfεf.The arrayΛ¯ is given by the following equation:(47)Λ¯=Ina×nρ+10na+nb+1×nρ∈ℝ2na+nb+1×nρ+1×na×nρ+1,and the matrix EΦ¯fΦ¯fT may be described by(48)EΦ¯fΦ¯fT=EΦfΦfEΦfu¯fEΦfu¯fEu¯fu¯f.Consequently, its inverseEΦ¯fΦ¯fT−1 is given by the following expression:(49)EΦ¯fΦ¯fT−1=EΦfΦf−1+X−EΦfΦf−1EΦfu¯fD−1−D−1EΦfu¯fTEΦfΦf−1D−1.where X=EΦfΦfT−1EΦfu¯fD−1EΦfu¯fTEΦfΦfT−1 and D=Eu¯fu¯f−EΦfu¯fTEΦfΦfT−1EΦfu¯fCombining (43) with (48),(50)b¯^=−D−1EΦfu¯fTEΦfΦfT−1EΦfyf+D−1Eu¯fyf=−D−1EΦfu¯fTθ^folsNd−Eu¯fyf.On the other hand, it follows from (37), (38), (45), and (49) and since the augmented polynomial B¯p has known parameters vector b¯=0, it is evident that(51)b¯^=b¯−D−1EΦfu¯fTEΦfΦfT−1ΛEYfεf,=−D−1EΦfu¯fTEΦfΦfT−1ΛEYfεf.Thus, (50) and (52) give the following equality:(52)−D−1EΦfu¯fTEΦfΦfT−1ΛEYfεf=−D−1EΦfu¯fTθ^folsNd−Eu¯fyf.From (52), it is clear that(53)EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1EΦfu¯fTθ^folsNd−Eu¯fyf.Replacing (53) in (28) leads to(54)Δ^θ=EΦfΦfT−1ΛEΦfu¯fTEΦfΦfT−1Λ−1EΦfu¯fTθ^folsNd−Eu¯fyf.Finally, we obtain(55)θ^fobclsNd=θ^folsNd−Δ^θ,=θ^folsNd−EΦfΦfT−1ΛEΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTθ^folsNd−Eu¯fyf.To facilitate the use of this method, a detailed description is given in Algorithm 1.


### 4.2.1. Convergence Analysis

We begin this part by demonstrating that the fobcls method is highly reliable.Theorem 1.
The proposedfobclsestimator for fractional order system identification gives consistent estimates. The optimal augmented estimator θfobclsNd verifies the following equation:(56)limNd⟶∞θ^fobclsNd=θ.Proof.
Firstly, we recall that the filtered output signal is expressed as(57)yft=ΦfTθ+εft.
Multiplying both side of equation (57) by u¯ft, we obtain(58)u¯ftyft=u¯ftΦfTθ
Because the noiseset and vt are mutually uncorrelated of each other and of ut, we get(59)Eu¯fyf=EΦfu¯ftθ.
So, using equation (53),(60)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTθ^folsNd−Eu¯fyf.
ReplacingEu¯fyf in the above equation by its expression given in equation (57), we obtain(61)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTθ^folsNd−EΦfu¯ftθ.
Thus, the result is shown as the following equation:(62)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦf−1Λ−1×EΦfu¯fTθ^folsNd−θ.
Now, replacing the termθ^folsNd−θ by the expression of Δ^θ defined in equation (54),(63)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTEΦfΦfT−1ΛEYfεf.
Consequently,(64)limNd⟶∞EYfεf=EYfεf.
Using equation (65) and equation (25), we find(65)limNd⟶∞θ^fobcls=θ^folsNd−EΦfΦfT−1ΛEYfεf=θ+EΦfΦf−1ΛEYfεf−EΦfΦfT−1ΛEYfεf.
Finally, lettingNd⟶∞, we obtain(66)limNd⟶∞θ^fobclsNd=θ.Algorithm 1: Fractional-order bias compensated least squares algorithm.
Step 1: Generate the data using the svfuft,yft|t=1,…,Nd.Step 2: Build the filtered regressor using equation (13)Step 3: Apply thefols method to obtain θ^fols using equation (43)Step 4: Compute the matrixΛ and the vector u¯ft and build the augmented regressor vector φf−Τt using equations (46), (42), and (41)Step 5: Compute the bias expression via equation (53)Step 6: Calculateθ^fobclsNd using equation (54)


## 4.2.1. Convergence Analysis

We begin this part by demonstrating that the fobcls method is highly reliable.Theorem 1.
The proposedfobclsestimator for fractional order system identification gives consistent estimates. The optimal augmented estimator θfobclsNd verifies the following equation:(56)limNd⟶∞θ^fobclsNd=θ.Proof.
Firstly, we recall that the filtered output signal is expressed as(57)yft=ΦfTθ+εft.
Multiplying both side of equation (57) by u¯ft, we obtain(58)u¯ftyft=u¯ftΦfTθ
Because the noiseset and vt are mutually uncorrelated of each other and of ut, we get(59)Eu¯fyf=EΦfu¯ftθ.
So, using equation (53),(60)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTθ^folsNd−Eu¯fyf.
ReplacingEu¯fyf in the above equation by its expression given in equation (57), we obtain(61)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTθ^folsNd−EΦfu¯ftθ.
Thus, the result is shown as the following equation:(62)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦf−1Λ−1×EΦfu¯fTθ^folsNd−θ.
Now, replacing the termθ^folsNd−θ by the expression of Δ^θ defined in equation (54),(63)limNd⟶∞EYfεf=EΦfu¯fTEΦfΦfT−1Λ−1×EΦfu¯fTEΦfΦfT−1ΛEYfεf.
Consequently,(64)limNd⟶∞EYfεf=EYfεf.
Using equation (65) and equation (25), we find(65)limNd⟶∞θ^fobcls=θ^folsNd−EΦfΦfT−1ΛEYfεf=θ+EΦfΦf−1ΛEYfεf−EΦfΦfT−1ΛEYfεf.
Finally, lettingNd⟶∞, we obtain(66)limNd⟶∞θ^fobclsNd=θ.Algorithm 1: Fractional-order bias compensated least squares algorithm.
Step 1: Generate the data using the svfuft,yft|t=1,…,Nd.Step 2: Build the filtered regressor using equation (13)Step 3: Apply thefols method to obtain θ^fols using equation (43)Step 4: Compute the matrixΛ and the vector u¯ft and build the augmented regressor vector φf−Τt using equations (46), (42), and (41)Step 5: Compute the bias expression via equation (53)Step 6: Calculateθ^fobclsNd using equation (54)


## 4.3. Fractional-Order Estimation

Thefobcls technique can only be used in the scenario where the fractional orders are considered to be known a priori. The main contribution of this part is to design an algorithm extending the identification method studied previously to a more normal situation in which the fractional-order differentiation is supposed to be unknown and the coefficients are approximated. Fractional-order optimization-bias compensated least squares foo-bcls is the name of this algorithm. It is based on a combination of the fobcls coefficient estimate technique and a nonlinear differentiation order optimization algorithm. As an option to tackling the optimization problem, the Gauss–Newton method is offered.For that, consider the parameter vectorθ as follows:(67)θ=a1,…,ana,b0,…,bnb,ν.The challenge of parameter identification is framed as a functional minimization. As a result, the primary purpose of this method is to reduce the error described by equation (69) with regard to ν. As a result, the Gauss–Newton method is used to repeatedly update the commensurate order. As a result, the quadratic criteria are as follows:(68)Wθ^foo−bcls=12εt2,where the error εt is expressed by(69)εt=yt−y^t.At each the iterationj+1, this algorithm calculates the fractional differentiation order νj+1(70)νj+1=νj−γℋ−1∂W∂ν|νj.The gradient∂W/∂ν is described as follows:(71)∂W∂ν=∂εTt∂νεt.and the approximated hessian ℋ is given by(72)ℋ=∂εTt∂ν∂εt∂ν.The description of thefoo-bcls method is given by Algorithm 2.Data {uf (t), yf(t)t = 1, . . . , Nd}(1)
InitializationInitializeν0 .The initial estimation of theθˆ is obtained by Applying Algorithm 1.Compute W (θˆfoo–bcls) (0)Initializeλ ∈ R∗ .(2)
repeatApply the Gauss–Newton algorithm in order to update the commensurate order using equation (71)The parameters vector is estimated by applying Algorithm 1.Calculate the quadratic criterionW(θ^foo−bclsj+1 using equation (69)−λ=λ/2(3)
untilΔ^θi+1−Δ^θi/Δ^θi+1<σorj=MaxIter where MaxIter designates the maximum number of iteration(4)
returnθˆfoo−bcls.


## 5. Numerical Example

The numerical examples provided in this section serve to verify the different aspects highlighted in the previous section of the paper.In order to give a quantitative evaluation of the proposed method in this paper, we adopted the performances indexes: the Best Fit Rate (BFR), which is defined as follows:(73)BFR=max1−∑t=1Ndyt−y^t2∑t=1Ndyt−y¯2,0,where y¯ is the mean of the output signal and y^ is the estimated model output and also the normalized relative quadratic error (NRQE) is expressed by(74)NRQE=1nmc∑k=1nmcθ^−θ2θ,where nmc designates the Monte Carlo realization number.The fractional-order optimization-least squaresfoo - ls and foo-bcls techniques are used to estimate fractional commensurate order and fractional LPV model coefficients. Moreover, these algorithms are compared with Nelder–Mead (NM) method which has been widely used in many scientific and engineering applications [43]. It is based on the Matlab fminsearch optimization function without constrained. In fact, the Matlab fminsearch is a nonlinear function for solving mathematical problems available on the matlab platform. It is based on the Simplex algorithm and used to optimize the solution by minimizing a given cost function [44]. It starts with the initial vector and tries to find the local minimum. The main reason for the popularity of this algorithm in practice is its simplicity.


### 5.1. Example 1

In this problem, the simulation studies are performed for the following system:(75)G1p,ρt=b0ρt1+a1ρtpν,where ν=0⋅400 is the dependence of the coefficients.The fractional LPV model input-output is given by(76)Ap,ρtyt=ℬp,ρtut+εt,where(77)Ap,ρt=1+a1ρtpα1,Bp,ρt=b0ρt,a1ρt=a1,0+a1,1ρt,b0ρt=b0,0+b0,1ρt,with(78)b0,0=1b0,1=1a1,0=1a1,1=1⋅5.Figure2 depicts the input signal ut as a pseudo random binary sequence (PRBS) with a uniform distribution between [−1 1]. With a sample time of T=0⋅05  seconds, the scheduling and output signals are monitored (see Figure 2). The sample count is Nd=6000. The output measurement is affected by a white noise et with varying signal-to-noise SNRy ratios, which is represented as(79)SNRy=10logvarytvaret,Figure 2
The fractional-order LPV system input and output measurement.Figure3 shows the observed scheduling signal defined by(80)ρt=0⋅5sin0⋅3π12t+vt,where vt is the white noise with a SNRρt of 20 [dB].Figure 3
The scheduling signal.


#### 5.1.1. The svf Parameters Selection

Our focus in this section is on the effects of the state variables parameters on thefoo-ls and foo-bcls algorithms. For the nmc=100 runs, a Monte Carlo simulation has been performed. In practice, finding the exact values of the svf parameters may be difficult, and iterative testing is typically required. NRQE will be used for this exam in this post. Accordingly, κ and λF can only be approximated based on the noise and order information.(1) The Choice ofλF. The choice of λF is too difficult so when it is too high, the noise is not filtered properly, and when it is too low, the system dynamics are filtered with the noise. A study on the effect of λF on the foo-ls and the foo-bcls methods is handled in this paragraph. Thus, the svf order κ=1 and λF is chosen as follows:(81)λF=0⋅01,0⋅5,1,10rad/seconds.The results are illustrated in Table1. The low NRQE values confirm that the foo-bcls performs good results when λF=0⋅5  rad/seconds.Table 1
Effect of the parameterλF on the foo-ls and foo-bcls methods (SNRρ=SNRy=20  [dB] and nmc=100  runs).
λFNRQEfoo−lsfoo-bcls0.010.5690.158 90.50.182 30.008 510.253 60.025 6100.985 00.269(2) Choice ofκ. At this level, λF=1  rad/seconds and the svf order κ is selected as(82)κ=1,2,3.Table2 shows that the simplicity of the implementation and the low value of NRQE is ensured by κ=1.Table 2
Effect of the parameterκ on the foo-ls and foo-bcls methods (SNRρ=SNRy=20  [dB] and nmc=100  runs).
κNRQEfoo−lsfoo−bcls10.182 30.008 520.152 80.004 530.137 60.00640.109 00.000 9


#### 5.1.2. Comparative Study

To show the relevance of the presented algorithm, a Monte Carlo simulation ofnmc=200 experiments is done. A Gaussian white noise is affected the output signal yt; the values SNRy1=12 [dB] and SNRy2=6 [dB] were considered.The fractional LPV model parameters vector is expressed as(83)θ=a1,0a1,1b0,0b0,1ν.λF=0⋅3 rad/seconds and κ=1.The results are recapitulated in Tables3 and 4 and illustrated graphically in Figure 4. It is evident from the values shown in Tables 3 and 4 that the foo-ls gives biased parameters estimates. On the other side, the foo-bcls method produces unbiased parameters estimates with a low NRQE, BFR ≃1, and low values of standard deviation. The observed good performance of this method was noticed even if in the case of a lower value of SNRy (SNRy=6[dB]).Table 3
Fractional system parameters estimation resulted byfoo-ls and the foo-bcls methods (SNRρ=20, SNRy=12, and 6 [dB] nmc=200  runs).
TrueSNR=12  dBSNR=6  dBfoo-lsfoo−bclsfoo-lsfoo−bclsa1,0=10⋅7740±0⋅01260⋅9989±0⋅01640⋅2764±0⋅01760⋅9960±0⋅0359a1,1=1⋅51⋅1256±0⋅00721⋅0007±0⋅00860⋅5842±0⋅01110⋅9980±0⋅0199b1,0=10⋅8706±0⋅00721⋅0007±0⋅00860⋅5842±0⋅01110⋅9980±0⋅0199b1,1=10⋅7683±0⋅01751⋅0001±0⋅00830⋅3315±0⋅03300⋅9979±0⋅0212ν=0⋅40⋅3068±0⋅01250⋅3999±0⋅00400⋅1513±0⋅02650⋅4053±0⋅0094Table 4
NRQE and BFR achieved byfoo-ls and the foo-bcls methods (SNRρ=20, SNRy=12 and 6 [dB], and nmc=200  runs).
MethodsSNR=12  dBSNR=6  dBNRQEBFRNRQEBFRfoo‐ls0.223 80.6280.658 00.473foo−bcls0.009 70.998 40.022 10.995 5Figure 4
The process estimates histograms obtained by thefoo-ls and the foo−bcls methods (SNRρ=20  [dB], SNRy=12and6  [dB], and nmc=200  runs). (a) foo-ls. (b) foo-bcls.
(a)(b)In the histograms below, we can clearly see that the obtained parameter estimates usingfoo-bcls are centered around the true ones, which confirms the consistency of this method.


### 5.2. Example 2

The model structure in this example is given by(84)G2p,ρt=b0ρt1+a1ρtpν+a2ρtp2ν.The fractional LPV model input-output is given as follows:(85)Ap,ρtyt=ℬp,ρtut+εt,where(86)Ap,ρt=1+a1ρtpν+a2ρtp2ν,ℬp,ρt=b0ρt,and (87)a1=a1,0+a1,1ρt,a2=a2,0+a2,1ρt,b0=b0,0+b0,1ρt,where ν=0⋅6 and the coefficients are described by the following equation:(88)a1ρt=0⋅5+0⋅25ρta2ρt=1+0⋅5ρtb0ρt=1+0⋅2ρt.In this example, the parameters vector of the fractional model is given by(89)θ=a1,0a1,1a2,0a2,1b0,0b0,1,ν.In the upcoming example, Figure5 shows that the input signal ut is, also, a PRBS distributed between [−1  1] and the output signal is observed with a sampling period T=0⋅1 seconds. The number of samples is Nd=4000. A Monte Carlo simulation of nmc=200 running with new noise realization is achieved at noise level SNRy=10  [dB].Figure 5
The input and the output signal.Figure6 presents the scheduling signal expressed as(90)ρt=1⋅10⋅5sin0⋅35π12t+vt,where vt is a white noise with a SNRρt of 20 [dB].Figure 6
The scheduling signal.We treated the same study developed in section (5.1.1), and the svf parameters are chosen according to the values of the NRQE, respectively, as follows:λF=0⋅3 rad/seconds and κ=2.Table5 and Table 6 summarize the results about means, standard deviation, BFR, and NRQE of the estimated LPV model parameters over 200 Monte Carlo simulations. Table 5 and Table 6 demonstrate that the foo-bcls method offers satisfying results in terms of unbiased parameters estimates with low value of the normalized relative quadratic error and standard deviation values.Table 5
Fractional system parameters estimation obtained by thefoo-bcls method and NM algorithm SNRρ=20,SNRy=10  [dB], and nmc=200  runs.
TrueSNR=10 dBfoo-bclsNM algorithma1,0  = 1.0001⋅0121±0⋅01570⋅8763±0⋅1124a1,1  = 0.5000⋅5057±0⋅02740⋅4590±0⋅0693a2,0  = 0.5000⋅5001±0⋅01230⋅6495±0⋅1173a2,1  = 0.2500⋅2514±0⋅03280⋅3224±0⋅0631b1,0  = 1.0001⋅0053±0⋅05050⋅9942±0⋅0125b1,1  = 0.2000⋅2023±0⋅11950⋅2091±0⋅0173ν  = 0.6000⋅5985±0⋅02270⋅5688±0⋅0286Table 6
NRQE and BFR obtained by thefoo-bcls method and NM algorithm (SNRρ=20, SNRy=10  [dB], and nmc=200 runs).
MethodsSNR=10 dBNRQEBFRfoo-bcls0.05860.9935NM algorithm0.07120.8986Figure7 illustrates the histograms of the estimates obtained by foo-bcls and NM algorithm. It checks that the estimates resulting using the foo-bcls method are close to the true parameters. On the other hand and despite its faster convergence, the NM algorithm gives biased estimates far for the true ones. Moreover, this algorithm requires an initial condition of a starting simplex, which can slow down the algorithm if it is badly chosen. In addition, the NM algorithm failed to converge in some of the Monte Carlo runs.Figure 7
The process estimates histograms of the estimates with thefoo-bcls method and NM algorithm (SNRρ=20  [dB], SNRy=10  [dB], and nmc=200  runs). (a) NM algorithm. (b) foo-bcls.
(a)(b)These results show that even with a large number of parameters, it is still possible to detect the coefficients as well as the fractional-order differentiation (case of the second example: seven parameters). This approach has the benefit of strong and faster convergence of the estimates to the true values when compared with other such methods.


## 5.1. Example 1

In this problem, the simulation studies are performed for the following system:(75)G1p,ρt=b0ρt1+a1ρtpν,where ν=0⋅400 is the dependence of the coefficients.The fractional LPV model input-output is given by(76)Ap,ρtyt=ℬp,ρtut+εt,where(77)Ap,ρt=1+a1ρtpα1,Bp,ρt=b0ρt,a1ρt=a1,0+a1,1ρt,b0ρt=b0,0+b0,1ρt,with(78)b0,0=1b0,1=1a1,0=1a1,1=1⋅5.Figure2 depicts the input signal ut as a pseudo random binary sequence (PRBS) with a uniform distribution between [−1 1]. With a sample time of T=0⋅05  seconds, the scheduling and output signals are monitored (see Figure 2). The sample count is Nd=6000. The output measurement is affected by a white noise et with varying signal-to-noise SNRy ratios, which is represented as(79)SNRy=10logvarytvaret,Figure 2
The fractional-order LPV system input and output measurement.Figure3 shows the observed scheduling signal defined by(80)ρt=0⋅5sin0⋅3π12t+vt,where vt is the white noise with a SNRρt of 20 [dB].Figure 3
The scheduling signal.


### 5.1.1. The svf Parameters Selection

Our focus in this section is on the effects of the state variables parameters on thefoo-ls and foo-bcls algorithms. For the nmc=100 runs, a Monte Carlo simulation has been performed. In practice, finding the exact values of the svf parameters may be difficult, and iterative testing is typically required. NRQE will be used for this exam in this post. Accordingly, κ and λF can only be approximated based on the noise and order information.(1) The Choice ofλF. The choice of λF is too difficult so when it is too high, the noise is not filtered properly, and when it is too low, the system dynamics are filtered with the noise. A study on the effect of λF on the foo-ls and the foo-bcls methods is handled in this paragraph. Thus, the svf order κ=1 and λF is chosen as follows:(81)λF=0⋅01,0⋅5,1,10rad/seconds.The results are illustrated in Table1. The low NRQE values confirm that the foo-bcls performs good results when λF=0⋅5  rad/seconds.Table 1
Effect of the parameterλF on the foo-ls and foo-bcls methods (SNRρ=SNRy=20  [dB] and nmc=100  runs).
λFNRQEfoo−lsfoo-bcls0.010.5690.158 90.50.182 30.008 510.253 60.025 6100.985 00.269(2) Choice ofκ. At this level, λF=1  rad/seconds and the svf order κ is selected as(82)κ=1,2,3.Table2 shows that the simplicity of the implementation and the low value of NRQE is ensured by κ=1.Table 2
Effect of the parameterκ on the foo-ls and foo-bcls methods (SNRρ=SNRy=20  [dB] and nmc=100  runs).
κNRQEfoo−lsfoo−bcls10.182 30.008 520.152 80.004 530.137 60.00640.109 00.000 9


### 5.1.2. Comparative Study

To show the relevance of the presented algorithm, a Monte Carlo simulation ofnmc=200 experiments is done. A Gaussian white noise is affected the output signal yt; the values SNRy1=12 [dB] and SNRy2=6 [dB] were considered.The fractional LPV model parameters vector is expressed as(83)θ=a1,0a1,1b0,0b0,1ν.λF=0⋅3 rad/seconds and κ=1.The results are recapitulated in Tables3 and 4 and illustrated graphically in Figure 4. It is evident from the values shown in Tables 3 and 4 that the foo-ls gives biased parameters estimates. On the other side, the foo-bcls method produces unbiased parameters estimates with a low NRQE, BFR ≃1, and low values of standard deviation. The observed good performance of this method was noticed even if in the case of a lower value of SNRy (SNRy=6[dB]).Table 3
Fractional system parameters estimation resulted byfoo-ls and the foo-bcls methods (SNRρ=20, SNRy=12, and 6 [dB] nmc=200  runs).
TrueSNR=12  dBSNR=6  dBfoo-lsfoo−bclsfoo-lsfoo−bclsa1,0=10⋅7740±0⋅01260⋅9989±0⋅01640⋅2764±0⋅01760⋅9960±0⋅0359a1,1=1⋅51⋅1256±0⋅00721⋅0007±0⋅00860⋅5842±0⋅01110⋅9980±0⋅0199b1,0=10⋅8706±0⋅00721⋅0007±0⋅00860⋅5842±0⋅01110⋅9980±0⋅0199b1,1=10⋅7683±0⋅01751⋅0001±0⋅00830⋅3315±0⋅03300⋅9979±0⋅0212ν=0⋅40⋅3068±0⋅01250⋅3999±0⋅00400⋅1513±0⋅02650⋅4053±0⋅0094Table 4
NRQE and BFR achieved byfoo-ls and the foo-bcls methods (SNRρ=20, SNRy=12 and 6 [dB], and nmc=200  runs).
MethodsSNR=12  dBSNR=6  dBNRQEBFRNRQEBFRfoo‐ls0.223 80.6280.658 00.473foo−bcls0.009 70.998 40.022 10.995 5Figure 4
The process estimates histograms obtained by thefoo-ls and the foo−bcls methods (SNRρ=20  [dB], SNRy=12and6  [dB], and nmc=200  runs). (a) foo-ls. (b) foo-bcls.
(a)(b)In the histograms below, we can clearly see that the obtained parameter estimates usingfoo-bcls are centered around the true ones, which confirms the consistency of this method.


## 5.1.1. The svf Parameters Selection

Our focus in this section is on the effects of the state variables parameters on thefoo-ls and foo-bcls algorithms. For the nmc=100 runs, a Monte Carlo simulation has been performed. In practice, finding the exact values of the svf parameters may be difficult, and iterative testing is typically required. NRQE will be used for this exam in this post. Accordingly, κ and λF can only be approximated based on the noise and order information.(1) The Choice ofλF. The choice of λF is too difficult so when it is too high, the noise is not filtered properly, and when it is too low, the system dynamics are filtered with the noise. A study on the effect of λF on the foo-ls and the foo-bcls methods is handled in this paragraph. Thus, the svf order κ=1 and λF is chosen as follows:(81)λF=0⋅01,0⋅5,1,10rad/seconds.The results are illustrated in Table1. The low NRQE values confirm that the foo-bcls performs good results when λF=0⋅5  rad/seconds.Table 1
Effect of the parameterλF on the foo-ls and foo-bcls methods (SNRρ=SNRy=20  [dB] and nmc=100  runs).
λFNRQEfoo−lsfoo-bcls0.010.5690.158 90.50.182 30.008 510.253 60.025 6100.985 00.269(2) Choice ofκ. At this level, λF=1  rad/seconds and the svf order κ is selected as(82)κ=1,2,3.Table2 shows that the simplicity of the implementation and the low value of NRQE is ensured by κ=1.Table 2
Effect of the parameterκ on the foo-ls and foo-bcls methods (SNRρ=SNRy=20  [dB] and nmc=100  runs).
κNRQEfoo−lsfoo−bcls10.182 30.008 520.152 80.004 530.137 60.00640.109 00.000 9


## 5.1.2. Comparative Study

To show the relevance of the presented algorithm, a Monte Carlo simulation ofnmc=200 experiments is done. A Gaussian white noise is affected the output signal yt; the values SNRy1=12 [dB] and SNRy2=6 [dB] were considered.The fractional LPV model parameters vector is expressed as(83)θ=a1,0a1,1b0,0b0,1ν.λF=0⋅3 rad/seconds and κ=1.The results are recapitulated in Tables3 and 4 and illustrated graphically in Figure 4. It is evident from the values shown in Tables 3 and 4 that the foo-ls gives biased parameters estimates. On the other side, the foo-bcls method produces unbiased parameters estimates with a low NRQE, BFR ≃1, and low values of standard deviation. The observed good performance of this method was noticed even if in the case of a lower value of SNRy (SNRy=6[dB]).Table 3
Fractional system parameters estimation resulted byfoo-ls and the foo-bcls methods (SNRρ=20, SNRy=12, and 6 [dB] nmc=200  runs).
TrueSNR=12  dBSNR=6  dBfoo-lsfoo−bclsfoo-lsfoo−bclsa1,0=10⋅7740±0⋅01260⋅9989±0⋅01640⋅2764±0⋅01760⋅9960±0⋅0359a1,1=1⋅51⋅1256±0⋅00721⋅0007±0⋅00860⋅5842±0⋅01110⋅9980±0⋅0199b1,0=10⋅8706±0⋅00721⋅0007±0⋅00860⋅5842±0⋅01110⋅9980±0⋅0199b1,1=10⋅7683±0⋅01751⋅0001±0⋅00830⋅3315±0⋅03300⋅9979±0⋅0212ν=0⋅40⋅3068±0⋅01250⋅3999±0⋅00400⋅1513±0⋅02650⋅4053±0⋅0094Table 4
NRQE and BFR achieved byfoo-ls and the foo-bcls methods (SNRρ=20, SNRy=12 and 6 [dB], and nmc=200  runs).
MethodsSNR=12  dBSNR=6  dBNRQEBFRNRQEBFRfoo‐ls0.223 80.6280.658 00.473foo−bcls0.009 70.998 40.022 10.995 5Figure 4
The process estimates histograms obtained by thefoo-ls and the foo−bcls methods (SNRρ=20  [dB], SNRy=12and6  [dB], and nmc=200  runs). (a) foo-ls. (b) foo-bcls.
(a)(b)In the histograms below, we can clearly see that the obtained parameter estimates usingfoo-bcls are centered around the true ones, which confirms the consistency of this method.


## 5.2. Example 2

The model structure in this example is given by(84)G2p,ρt=b0ρt1+a1ρtpν+a2ρtp2ν.The fractional LPV model input-output is given as follows:(85)Ap,ρtyt=ℬp,ρtut+εt,where(86)Ap,ρt=1+a1ρtpν+a2ρtp2ν,ℬp,ρt=b0ρt,and (87)a1=a1,0+a1,1ρt,a2=a2,0+a2,1ρt,b0=b0,0+b0,1ρt,where ν=0⋅6 and the coefficients are described by the following equation:(88)a1ρt=0⋅5+0⋅25ρta2ρt=1+0⋅5ρtb0ρt=1+0⋅2ρt.In this example, the parameters vector of the fractional model is given by(89)θ=a1,0a1,1a2,0a2,1b0,0b0,1,ν.In the upcoming example, Figure5 shows that the input signal ut is, also, a PRBS distributed between [−1  1] and the output signal is observed with a sampling period T=0⋅1 seconds. The number of samples is Nd=4000. A Monte Carlo simulation of nmc=200 running with new noise realization is achieved at noise level SNRy=10  [dB].Figure 5
The input and the output signal.Figure6 presents the scheduling signal expressed as(90)ρt=1⋅10⋅5sin0⋅35π12t+vt,where vt is a white noise with a SNRρt of 20 [dB].Figure 6
The scheduling signal.We treated the same study developed in section (5.1.1), and the svf parameters are chosen according to the values of the NRQE, respectively, as follows:λF=0⋅3 rad/seconds and κ=2.Table5 and Table 6 summarize the results about means, standard deviation, BFR, and NRQE of the estimated LPV model parameters over 200 Monte Carlo simulations. Table 5 and Table 6 demonstrate that the foo-bcls method offers satisfying results in terms of unbiased parameters estimates with low value of the normalized relative quadratic error and standard deviation values.Table 5
Fractional system parameters estimation obtained by thefoo-bcls method and NM algorithm SNRρ=20,SNRy=10  [dB], and nmc=200  runs.
TrueSNR=10 dBfoo-bclsNM algorithma1,0  = 1.0001⋅0121±0⋅01570⋅8763±0⋅1124a1,1  = 0.5000⋅5057±0⋅02740⋅4590±0⋅0693a2,0  = 0.5000⋅5001±0⋅01230⋅6495±0⋅1173a2,1  = 0.2500⋅2514±0⋅03280⋅3224±0⋅0631b1,0  = 1.0001⋅0053±0⋅05050⋅9942±0⋅0125b1,1  = 0.2000⋅2023±0⋅11950⋅2091±0⋅0173ν  = 0.6000⋅5985±0⋅02270⋅5688±0⋅0286Table 6
NRQE and BFR obtained by thefoo-bcls method and NM algorithm (SNRρ=20, SNRy=10  [dB], and nmc=200 runs).
MethodsSNR=10 dBNRQEBFRfoo-bcls0.05860.9935NM algorithm0.07120.8986Figure7 illustrates the histograms of the estimates obtained by foo-bcls and NM algorithm. It checks that the estimates resulting using the foo-bcls method are close to the true parameters. On the other hand and despite its faster convergence, the NM algorithm gives biased estimates far for the true ones. Moreover, this algorithm requires an initial condition of a starting simplex, which can slow down the algorithm if it is badly chosen. In addition, the NM algorithm failed to converge in some of the Monte Carlo runs.Figure 7
The process estimates histograms of the estimates with thefoo-bcls method and NM algorithm (SNRρ=20  [dB], SNRy=10  [dB], and nmc=200  runs). (a) NM algorithm. (b) foo-bcls.
(a)(b)These results show that even with a large number of parameters, it is still possible to detect the coefficients as well as the fractional-order differentiation (case of the second example: seven parameters). This approach has the benefit of strong and faster convergence of the estimates to the true values when compared with other such methods.


## 6. Conclusion

In this paper, a method on system identification using the least squares algorithm for fractional-order linear parameters varying systems is proposed. The estimates obtained by the least squares algorithm are biased. To compensate the measurement noise in the output signal, an effective method based on the well-known bias correction scheme is established. The calculation formula of the bias compensation is offered to obtain the satisfactory estimations of the model coefficients and fractional-order differentiation. The identification findings support the suggested method’s powerful capacity in managing with fractional-order LPV systems.



---

*Source: 7278157-2022-04-01.xml*