基于改進粒子群算法的分數(shù)階系統(tǒng)辨識方法

王 強，趙志誠，桑 博(太原科技大學電子信息工程學院，太原 030024)

分數(shù)階微積分與整數(shù)階微積分誕生于同一時期，距今已有三百多年的歷史。隨著計算機技術(shù)的不斷發(fā)展，分數(shù)階微積分理論應用到了許多研究領(lǐng)域，如：電化學[1]，熱力學[2]，生物學[3]，控制理論[4]等。研究表明，用分數(shù)階微積分對實際動態(tài)對象和過程建模的好處很明顯，分數(shù)階模型具有更低的階數(shù)，更少的參數(shù)，更高的建模精度等優(yōu)點[5]，而且傳統(tǒng)的整數(shù)階系統(tǒng)模型是無法代替分數(shù)階系統(tǒng)模型的[6]。目前對于分數(shù)階系統(tǒng)辨識的研究尚處于初始階段，沒有統(tǒng)一的辨識方法。由于分數(shù)階微積分的物理意義尚不能確定，目前幾乎所有的分數(shù)階系統(tǒng)辨識方法均不考慮模型結(jié)構(gòu)參數(shù)的辨識，均是通過先驗知識或者整數(shù)階模型假設(shè)分數(shù)階模型結(jié)構(gòu)參數(shù)。文獻[7]通過極大似然算法得到等效準則函數(shù)，將參數(shù)辨識問題轉(zhuǎn)化為參數(shù)優(yōu)化問題，模型誤差主要集中在低頻段；文獻[8]提出了基于主元分析的分數(shù)階系統(tǒng)子空間辨識方法，文獻[9]基于分數(shù)階Riemann-Liouville定義把調(diào)制函數(shù)方法擴展到分數(shù)階系統(tǒng)辨識，計算量較大。文獻[10]對含有一個或者兩個偽極點的分數(shù)階模型結(jié)構(gòu)提出一種辨識方法，把一個無限維問題轉(zhuǎn)化到有限維范圍，但不能向更多偽極點分數(shù)階模型推廣。

本文提出一種基于改進的粒子群算法實現(xiàn)分數(shù)階系統(tǒng)辨識的方法，可實現(xiàn)模型參數(shù)和階次同時辨識。辨識精度高，而且易于編程實現(xiàn)，在有輸出隨機干擾的情況下也能辨識到精確的結(jié)果。

1 分數(shù)階數(shù)學描述

分數(shù)階微積分的發(fā)展經(jīng)歷了一個漫長的過程。它可以看作是整數(shù)階微積分的推廣。

1.1 分數(shù)階微積分定義

在分數(shù)階微積分理論發(fā)展過程中，出現(xiàn)了多種分數(shù)階微積分的定義。常用的有Cauchy定義，Riemann-Liouville定義，Grumwald-Letnikov定義。

Cauchy分數(shù)階積分定義：

(1)

式中，q為分數(shù)階積分階次，C為f(t)單值與解析開區(qū)域的光滑曲線，Γ為伽馬函數(shù)。

Riemann-Liouville分數(shù)階積分定義為：

(2)

(3)

(4)

式中，h為步長，[X]表示對X的取整。當q>0，該式表示分數(shù)階微分；當q<0，該式表示分數(shù)階積分。

1.2 分數(shù)階系統(tǒng)的數(shù)學模型

分數(shù)階系統(tǒng)就是利用分數(shù)階微分方程來描述其動態(tài)特性的系統(tǒng)。線性分數(shù)階SISO系統(tǒng)的微分方程一般表達式：

anDαny(t)+an-1Dαn-1y(t)+…+a0y(t)=

bmDβmu(t)+bm-1Dβm-1u(t)+…+b0u(t)

(5)

式中，ai、bi為實數(shù)，微分階次αi、βi可以為小數(shù)或者分數(shù)，u(t)為系統(tǒng)輸入信號，y(t)為系統(tǒng)輸出信號。在零初始條件下，其傳遞函數(shù)表達式為：

(6)

目前，辨識一般的分數(shù)階系統(tǒng)是很困難的，因為系統(tǒng)的階數(shù)很難確定。文獻[11]提出連續(xù)階次分布的概念，把分數(shù)階系統(tǒng)一般表達式轉(zhuǎn)化成公因子階次表達式：

anDnαy(t)+an-1D(n-1)αy(t)+…+a0y(t)=

bmDmαu(t)+bm-1D(m-1)αu(t)+…+b0u(t)

(7)

式中，α為公因子階次。

2 權(quán)值自適應粒子群算法

2.1 基本粒子群算法

粒子群優(yōu)化算法采用“群體”和“進化”的概念，依據(jù)個體的目標函數(shù)值大小來確定最優(yōu)解。

假設(shè)N維空間中，m個粒子組成一個群體，第i個粒子的位置記為Xi=(xi1，xi2，…，xiN)，其所經(jīng)歷的最好位置記為pbest=(pi1，pi2，…，piN)，運動速度記為Vi=(vi1，vi2，…，viN)，全局極值gbest=(gi1，gi2，…，giN)，則在迭代過程中，粒子按照式(8)、(9)更新自身速度與位置：

(8)

(9)

式中，i=1，2，3…，m；n=1，2，3…，N；c1、c2為學習因子，rand1、rand2是兩個相互獨立的0到1之間的隨機數(shù)。

式(8)表明粒子的速度更新由當前狀態(tài)、自身經(jīng)驗、和社會學習三部分組成，不僅平衡了局部和全局搜索能力，避免局部極小，同時使粒子間信息共享，最終到達最優(yōu)位置。將粒子的位置Xi帶入目標函數(shù)就得到相應的適應值fi=f(Xi).

2.2 權(quán)值自適應策略

基本粒子群算法存在容易陷入局部最優(yōu)的問題，而且ω的大小會影響到算法的全局和局部搜索能力，因此可以通過調(diào)整慣性權(quán)值ω來改進算法。

本文以所有粒子當前目標函數(shù)值的平均值作為粒子評價指標，把粒子群分為三個不同的子種群，每個子群根據(jù)粒子的目標函數(shù)值自適應的調(diào)整慣性權(quán)值ω，進而增強種群多樣性以提高算法的全局和局部搜索能力。

假設(shè)m個粒子的粒子群，當前所有粒子目標函數(shù)值的平均值記為favg：

(10)

式中，fi為粒子的當前適應值。

慣性權(quán)值ω值的自適應策略如下：

(11)

(12)

式中，ωmax為ω的最大值，iter為當前迭代次數(shù)，itermax為允許迭代最大次數(shù)。

3)當favg

(13)

式中，k1用來控制ω的上限，k1越大ω的上限越大，為使式(13)的ω值能夠大于1，本文取k1=1.5，則ω∈(0.5，1.1)；k2控制式(13)的調(diào)節(jié)能力，k2過大，ω會迅速越小，使得算法的全局搜索能力不足，k2過小，ω變化不明顯，容易導致算法陷入局部最優(yōu)，本文選取k2=0.9.

3 分數(shù)階系統(tǒng)辨識

文中將改進的粒子群算法應用到分數(shù)階系統(tǒng)辨識過程中，辨識步驟如下：

1)確定搜索向量，即待辨識的參數(shù)，包括模型參數(shù)和階次。

2)選擇目標函數(shù)。將單位階躍信號加入到實際系統(tǒng)中，得到輸出響應y(t)，將單位階躍信號加入到待辨識模型中，得到輸出響應為yM(t)，選取算法的適應度函數(shù)為，

3)算法初始化。初始化算法相關(guān)參數(shù)，產(chǎn)生隨機搜索向量。

4)確定截止條件。截止條件選擇最大迭代次或者目標函數(shù)值。

4 仿真實例

為檢驗所提分數(shù)階系統(tǒng)辨識方法的有效性，針對幾種不同模型做如下仿真，并考慮了無輸出干擾和有輸出干擾兩種情況。仿真取迭代次數(shù)100次為算法截止條件，即itermax=100.算法中參數(shù)設(shè)定如下：c1=c2=1.5，ωmin=0.4，ωmax=0.9.粒子初始化范圍在區(qū)間(0，3)內(nèi)。假設(shè)辨識的模型結(jié)構(gòu)為：

仿真過程中，將本文方法與文獻[12]的方法進行比較，為避免算法初始值對比較結(jié)果的影響，兩種辨識方法采用同一組初始值。

1)輸出響應無干擾，即在理想情況下通過上述方法對系統(tǒng)進行辨識。

利用本文方法得到的系統(tǒng)辨識結(jié)果為：目標函數(shù)值f=3.086 2e-28，a1=2.300 0，a0=1.800 0，β=0.800 0.模型的單位階躍響應如圖1，目標函數(shù)值的對數(shù)隨迭代次數(shù)變化曲線如圖2.由圖可見，辨識模型與實際系統(tǒng)的單位階躍響應曲線高精度擬合，算法迭代20次時，辨識的誤差值為0.001 6，文獻[12]辨識的誤差值為1.214 8，由此表明，本文方法僅需較少的迭代次數(shù)就可辨識到期望的系統(tǒng)模型。

圖1 模型Ⅰ及其辨識模型的階躍響應

圖2 模型Ⅰ目標函數(shù)值的對數(shù)變化曲線

該模型為成比例分數(shù)階系統(tǒng)模型，取α=2β，利用本文方法得到的系統(tǒng)辨識結(jié)果為：

f=7.684 6e-23，a2=2.300 0，a1=1.400 0，a0=1.000 0，β=1.100 0.如圖3所示，實際系統(tǒng)與辨識模型的單位階躍響應曲線幾乎完全擬合。同樣，由圖4可見，對于成比例分數(shù)階系統(tǒng)，本文方法優(yōu)于文獻[12]的方法。

該模型為不成比例分數(shù)階系統(tǒng)模型。利用本文方法得到的系統(tǒng)辨識結(jié)果為：目標函數(shù)值f=4.756 4e-20，a2=2.300 0，a1=1.400 0，a0=1.000，α=2.200 0，β=0.900 0，其階躍響應及目標函數(shù)值變化如圖5、6所示。由圖可見，實際模型與辨識模型的單位階躍響應曲線高精度擬合。對于不成比例分數(shù)階系統(tǒng)，本文方法優(yōu)于文獻[12]的方法，能夠得到很好的辨識效果。

圖3模型Ⅱ及其辨識模型的階躍響應

Fig.3StepresponseofmodelⅡanditsidentificationmodel

圖6 模型Ⅲ目標函數(shù)值的對數(shù)變化曲線

圖8 有干擾時模型Ⅲ目標函數(shù)值的對數(shù)變化曲線

2)輸出響應有干擾，在輸出響應y(t)上疊加-0.03～0.03的隨機干擾。下面僅對模型Ⅲ進行有隨機干擾的仿真分析。

利用本文方法得到的系統(tǒng)辨識結(jié)果為：目標函數(shù)值f=0.091 0，a2=2.287 0，a1=1.366 8，a0=0.999 3，α=2.187 6，β=0.891 5.如圖7，有隨機干擾的實際系統(tǒng)單位階躍響應與辨識模型的單位階躍響應曲線基本一致。同樣，由圖8表明，在輸出響應有干擾的情況下，本文方法明顯優(yōu)于文獻[12]方法，在迭代到10次之后，本文方法辨識的誤差值明顯小于文獻[12]方法辨識的誤差值。

5 結(jié)論

按照粒子目標函數(shù)值的不同，將粒子群分為三個子種群，引入慣性權(quán)值自適應律來改進基本粒子群算法，并利用改進的粒子群算法提出一種分數(shù)階系統(tǒng)辨識方法，實現(xiàn)對分數(shù)階系統(tǒng)的模型參數(shù)和階次參數(shù)同時辨識。方法選取實際系統(tǒng)與辨識系統(tǒng)的輸出誤差平方和為目標函數(shù)。仿真結(jié)果表明，本文方法可有效辨識成比例和不成比例分數(shù)階系統(tǒng)，算法辨識精度較高。

參考文獻：

[1] OLDHAM K B.Fractional differential equations in electro-chemistry[J].Advances in Engineering Software，2010，41(1)：9-12.

[2] GABANO J D，POINOT T.Fractional modeling and identification of thermal systems[J].Signal Processing，2011，91(3)：531-541.

[3] MAGIN R L.Fractional calculus models of complex dynamics in biological tissues[J].Computers & Mathematics with Applications，2010，59(5)：1586-1593.

[4] 桑海，趙志誠，孟香，等.交流伺服系統(tǒng)的一種PIλDμ控制器設(shè)計方法[J].太原科技大學學報，2013，34(6).

[5] 蘇密勇，譚永紅，王子民，等.同元次分數(shù)階模型的一種具有穩(wěn)定約束的頻域辨識算法[J].四川大學學報：工程科學版，2012，44(2)：130-134.

[6] ZHANG X，CHEN Y Q.Remarks on fractional order control systems[C]∥American Control Conference (ACC).Montreal，Canada，2012：5169-5173.

[7] LI W，PENG C，WANG Y.Identifying the fractional-order systems with frequency responses：A maximum likelihood algorithm[C]∥International Conference on Radar.Cheng Du，China，2011，2：1943-1948.

[8] ZENG L，CHENG P，YONG W.Subspace identification for commensurate fractional order systems using instrumental variables[C]∥Proceedings of the 30th Chinese Control Conference (CCC).Yantai，China，2011：1636-1640.

[9] LIU D Y，LALEG-KIRATI T M，GIBARU O，et al.Identification of fractional order systems using modulating functions method[C]∥American Control Conference (ACC).Washington DC，USA，2013：1679-1684.

[10] NAZARIAN P，HAERI M.Identification of fractional order model structures using conversion of infinite terms of a response to finite terms[C]∥7th International Conference on Electrical and Electronics Engineering (ELECO).Bursa，Turkey，2011:391-395.

[11] HARTLEY T T，LORENZO C F.Fractional-order system identification based on continuous order-distributions[J].Signal processing，2003，83(11)：2287-2300.

[12] MAITI D，CHAKRABORTY M，KONAR A.A novel approach for complete identification of dynamic fractional order systems using stochastic optimization algorithms and fractional calculus[C]∥International Conference on Electrical and Computer Engineering.Dhaka，Bengal，2008：867-872.