基于樣條插值的不完備信息系統參數估計

李 珍,魏利勝?,程運昌(.安徽工程大學電氣工程學院,安徽蕪湖,4000; .國網霍邱縣供電有限責任公司,安徽六安,37400)

基于樣條插值的不完備信息系統參數估計

李 珍1,魏利勝1?,程運昌2
(1.安徽工程大學電氣工程學院,安徽蕪湖,241000; 2.國網霍邱縣供電有限責任公司,安徽六安,237400)

研究一種不完備信息系統參數辨識插補法.首先將每個輸入信號加上時間戳,即將輸入信號進行時間順序排序,使數據在特定的時刻缺失;并采用樣條插補法,自適應補充缺失數據,即當缺失一個輸出數據時,應用線性樣條插值法填補缺失數據,當連續缺失兩個或兩個以上的數據時,應用三次樣條插值法進行補充;在此基礎上應用粒子群算法對插補后的不完備信息系統進行參數精確估計;最后,通過MATLAB仿真實驗驗證了該方法在不同缺失率下,其參數估計的精確性與有效性.

不完備信息;樣條插值;粒子群算法;參數估計

系統辨識已廣泛應用于各行各業[1].系統的參數估計是根據已知的輸入信號與輸出數據來估計此系統中控制器的參數值.在此過程中,往往會因為儀器或人為因素導致某些數據的丟失.如果直接使用存在丟失的數據來估計參數,結果顯然會存在偏差[2].因此,必須采取一種方法如插補法將缺失的數據補充完整,再進行參數的估計[3].目前,由于數據缺失使得很多系統研究都不能更順利地進行,因此插補法已成為廣大研究者們關注的焦點.加權估計法、似然插補法以及多重插補法等被廣泛運用[4];Cheti Nicoletti[5]對不同插補方法進行了比較;Hai Zhong[6]對多重插補法又做了新的探討;林東方[7]對期望最大化算法進行了探討;R Pintelon[8]等將缺失的數據當作未知的參數;而譚彥華[9]將樣條插值法應用到模糊系統中.

對于插補法的研究已取得了一定的成績,但是,當缺失多個數據時,回歸插補法效果不好;缺失單個數據時,多重插補法的效果不好;缺失率大的時候,似然插補法效果不好.加權插補法需要假設同一調整組內的單元具有相等的回答概率[2].而當缺失單一的數據時,采取線性的樣條插補法.當缺失兩個或兩個以上的數據時,就應用三次樣條插值法.樣條插值可以在缺失率相對較大時仍具有良好的插補效果,且不需要添加辨識的參數.而粒子群算法的程序較簡單,需要修改的參數很少.因此,采用樣條插值法與粒子群算法的結合來實現對不完備信息系統參數的較精確估計.

1 樣條插值法基本原理

根據完備信息系統參數估計的原理框圖來畫出不完備信息系統參數估計的原理框圖如圖1所示.在已知系統輸入信號u(t)與噪聲信號v(t)的前提下,應用樣條插補法將存在數據缺失的輸出補充完整,得到輸出值z(t).將插補后的系統輸出值z(t)與選擇的數學模型輸出值^Z(t)之間的誤差e(t)返回到粒子群算法,粒子群算法辨識出一組參數,將該組參數帶入數學模型,重新產生一個數學模型輸出值^Z(t).重復取兩者之間的誤差,直到達到事先規定的迭代次數或者誤差e(t)在事先規定的精度范圍內時,停止運算,此時對應的參數組就是該算法辨識的參數大小.

1.1 樣條插值法基本原理

根據已知的數據對,尋求一組較好的擬合多項式[10].運用該多項式來擬合每對相鄰數據點之間的函數曲線,然后求出缺失值[3].將區間[a,b]劃分成n等份,即a＝t1＜…＜tn＝b,稱為Ω.函數z(t)需要滿足:

線性樣條插值法:

當z(t)在每對[ti,ti＋1](i＝1,2,…,n－1)上的多項式次數不超過1次,并且滿足1個以上的子區間函數最高次數為1次,則稱z(t)是關于Ω的線性樣條函數.線性樣條函數的表達式如下:

三次樣條插值法:

當z(t)在每對區間[ti,ti＋1](i＝1,2,…,n－1)上都是次數不超過3次的多項式,并且最少在1個子區間上函數最高次數為3次,就稱z(t)是關于Ω劃分的三次樣條函數.令經過已知n個樣本數據點(ti, Yi)(i＝1,2,…,n),t1＜t2＜…＜tn的三次樣條函數為z(t)[9].令:其中,未知系數ai,bi,ci,di需要滿足條件:z(ti)＝Yi,z(ti－0)＝z(ti＋0)和z′(ti－0)＝z′(ti＋0),z″(ti－0)＝z″(ti＋0).此處添加的是第2類邊界條件即z″(a)＝Y″(a),z″(b)＝Y″(b).又令z″(ti)＝Ni,z″(ti＋1)＝Ni＋1,令hi＝ti＋1－ti[8].將z″(t)進行兩次積分得到如下的表達式:

根據z(ti)＝Yi,z(ti＋1)＝Yi＋1,可得:

保持z(t)在節點處的連續性,即z′(ti－0)＝z′(ti＋0)和z′i(ti－0)＝z′i＋1(ti＋0).令再把第一邊界條件z′(t1)＝Y′1,z′(tn)＝Y′n帶入z′(t)得出下列的方程組:

1.2 粒子群優化算法基本原理

設實函數F是定義在歐式空間ED的某一區域S上的.空間維數為D,粒子個數為L,那么第i個粒子在該區域S中的位置、速度以及該時刻的適應度分別表示為:X 0i＝(X 0i1,X 0i2,…,X 0iD)T∈S,Vi＝(Vi1,Vi2,…,ViD)和Fitnessi＝F(X 0i),其中i＝1,2,…,L.第i個粒子尋求到曾經的最佳位置和對應的最佳適應度分別是:Y 0i和p(i)[1].第i個粒子搜索到的全部粒子中的最優位置和對應的最佳適應度分別表示為:pg和fitness(D,pg).引入慣性權重w的粒子群算法的步驟如文獻[1]中所示.

其中d＝1,2,…,D,每個粒子的位置與速度按下式變化如下:

收斂性證明:

其中w∈(0,1)的隨機數,同時c1,c2∈(0,2)的隨機數.Vi、X 0i是相互獨立的,因此可以對其一維進行分析.由于每一個粒子的自身歷史最優位置Y 0與全部粒子中的最優位置pg是保持不變的[11],因此可以將Y 0、pg、w、c1及c2當作常數.化簡式(8)、式(9)得:

又

將式(10)、式(11)帶入式(12)得出:

運用特征方程法來求解式(13)c.＝c1＋c2,此一元二次方程的求解有3種情況:當Δ＝0時,λ＝λ1＝λ2＝此時X 0(i)＝(A0＋A1?i)?λi,其中A0,A1為待定的系數;當Δ＞0時,λ1,2＝此時X 0(i)＝A0＋A1?λi1＋A2?λi2;其中待定系數為A0,A1,A2;若Δ＜0時,那么X 0(i)＝A0＋A1?λi1＋A2?λi2,其中待定系數為A0,A1,A2.

2 基于樣條插值的不完備信息的參數辨識

對于不完備信息的一類單輸入單輸出線性系統[13]基本結構圖如圖2所示.由圖2可知,系統的輸入信號、噪聲信號分別為u1(k)、v1(k).v1(k)是均值為零,方差為σ2的高斯白噪聲.系統的傳遞函數用G(z－1)來表示,輸出用Z(k)來表示.存在數據缺失的輸出是Z0(k).u1(k)與Z0(k)是可以觀測的.z(k)是采用插補法插補后的系統輸出.一般情況下有下列表達式:

當滿足如式(17)關系時,該系統是滑動模型.

已知參數的階次na,nb和nd,令:

化簡式(17),得到該系統的最小二乘格式:

運用粒子群算法將辨識的問題轉化成優化問題來思考.由式(14)可表示出系統的偏差準則函數:

粒子群算法是從進化的方面來考慮問題的,Pbest(θ)是該算法的適應度函數.由于式(20)是在線的優化問題,因此運用PSO(Particle Swarm Optimization,粒子群優化)算法探求滿足適應度函數值最優(即值最小)的參數,這些參數就是需要估計的參數大小.

3 仿真實驗

仿真實驗采取MATLAB軟件來完成.首先在完整數據的條件下,采用最小二乘遺傳遞推算法、PSO算法產生輸入信號、噪聲信號、輸出信號、初始隨機位置矩陣和對應的速度矩陣以及一些輔助信號,在不同缺失率下產生不同的輸出,結合保存下來的數據,分別采用樣條插值粒子群算法和無插值最小二乘遺傳遞推算法來實現該類系統的參數辨識[13].

對于該系統的樣條插補粒子群優化算法和無插補最小二乘遺傳遞推算法的仿真結果如式(21)所示:

式中,v(k)是零均值的白噪聲,輸入u(k)采用6階M序列,幅度為1,數據個數L＝100,D＝5,加權陣取單位陣.給定初始條件:beita＝0.999,miu＝beita^2,P 0＝10^6?eye(D),theta 0＝10^(－3)?ones(D,1), P(1:D,1:D,1)＝P 0,theta(1:D,1)＝theta 0,A(k)＝1.學習的因子:c1＝1.496 2,c2＝1.496 2,其慣性權重w＝0.729 8,最大循環次數Max DT＝100,待估計參數個數為D.滿足w＜1,c＞0與2?w－c＋2＞0,所以算法收斂.最后比較兩種插補方法插補后系統的參數估計、性能函數.以下給出仿真的所有數據與圖形,為樣條插值法在系統辨識領域的應用打下基礎.

模型產生的隨機輸入信號、噪聲信號如圖3所示,圖3中噪聲信號的不同線條個數代表了不同的維數,噪聲信號是一個L×L維的信號.

不同缺失率下不同插補辨識法的參數估計值如表1所示.由表1可知,在相同的數據缺失率下,采取樣條插值粒子群算法使不完備系統的參數估計較準確、收斂速度較快、穩定性也較強.5%缺失率時,參數及性能函數的變化情況如圖4所示.10%缺失率時,參數及性能函數的變化情況如圖5所示.20%缺失率時,參數及性能函數的變化情況如圖6所示.30%缺失率時,參數及性能函數的變化情況如圖7所示.在圖4～圖7中,a1表示參數a1的真實值;a11表示參數a1在各缺失率下采用樣條插值粒子群算法進行辨識的估計值;a12表示參數a1在各缺失率下采用文獻[14]方法即無插值最小二乘遺傳遞推算法進行辨識的估計值;Pbest表示該系統在各缺失率下采用樣條插值粒子群算法進行辨識的性能函數;J表示該系統在各缺失率下采用無插值最小二乘遺傳遞推算法進行辨識的性能函數.圖4中其他參數的表示同a1的各種表示意義是一樣的.

在相同缺失率下,采用樣條插值粒子群算法估計的參數值更接近真實值,且其收斂的速度也更快.由表1可知,在缺失率為5%時,對于參數a1,采用樣條插值粒子群算法估計值a11為－1.499 3,而采用無插值最小二乘遺傳遞推算法的估計值a12為－1.107 1;對于參數a2,采用樣條插值粒子群算法估計值a21為0.698 4,而采用無插值最小二乘遺傳遞推算法的估計值a22為0.366 1;對于參數b1,采用樣條插值粒子群算法估計值b11為0.982 4,而采用無插值最小二乘遺傳遞推算法的估計值b12為1.106 0;對于參數b2,采用樣條插值粒子群算法估計值b21為0.495 8,而采用無插值最小二乘遺傳遞推算法的估計值b22為0.397 7;對于參數d1,采用樣條插值粒子群算法估計值d11為0.932 7,而采用無插值最小二乘遺傳遞推算法的估計值d12為0.004 5.從數字可以直觀看出,樣條插值粒子群算法的估計值更準確.表格中的其他缺失率下的數據同樣可以得出此結論.

表1 不同缺失率下不同插補辨識法的參數估計值(迭代100次)

在相同缺失率下,由圖5可知,在缺失率為10%時,對于參數a1采用樣條插值粒子群算法估計值a11,在進化次數達到40次左右的時候就開始趨向收斂,且波動比較平緩,最終收斂于－1.473 4;而采用無插值最小二乘遺傳遞推算法的估計值a12,在進化次數達到60次左右的時候才開始趨向收斂,且波動不夠平緩,最終收斂于－0.807 1.在缺失率為10%時,對于參數a2采用樣條插值粒子群算法估計值a21,在進化次數達到30次左右的時候就開始趨向收斂,且波動比較平緩,最終收斂于0.665 5;而采用無插值最小二乘遺傳遞推算法的估計值a22,在進化次數達到40次左右的時候才開始趨向收斂,且波動不夠平緩,最終收斂于0.096 5.在缺失率為10%時,對于參數b1采用樣條插值粒子群算法估計值b11,在進化次數達到35次左右的時候就開始趨向收斂,且波動比較平緩,最終收斂于0.905 2;而采用無插值最小二乘遺傳遞推算法的估計值b12,在進化次數達到70次左右的時候才開始趨向收斂,且波動不夠平緩,最終收斂于0.987 6.在缺失率為10%時,對于參數b2采用樣條插值粒子群算法估計值b21,在進化次數達到35次左右的時候就開始趨向收斂,且波動比較平緩,最終收斂于0.479 2;而采用無插值最小二乘遺傳遞推算法的估計值b22在進化次數達到60次左右的時候才開始趨向收斂,且波動不夠平緩,最終收斂于0.713 4.在缺失率為10%時,對于參數d1采用樣條插值粒子群算法估計值d11,在進化次數達到35次左右的時候就開始趨向收斂,且波動比較平緩,最終收斂于0.876 6;而采用無插值最小二乘遺傳遞推算法的估計值d12,在進化次數達到60次左右的時候才開始趨向收斂,且波動不夠平緩,最終收斂于0.038 5.

從表1和參數仿真圖都可以看出,隨著缺失率的增大,辨識的效果越來越差.由表1可知,同樣采用樣條插值粒子群算法進行參數估計,對于參數a1,在缺失率為5%時,估計的值a11為－1.499 3;在缺失率為10%時,估計的值a11為－1.473 4;在缺失率為20%時,估計的值a11為－1.383 9;在缺失率為30%時,估計的值a11為－1.045 6.從數字的變化可以證明該結論.同樣的,由圖4b、圖5b、圖6b和圖7b可以看出,同樣采用樣條插值粒子群算法進行參數估計,對于參數b1,在缺失率為5%時,進化次數達到28次左右就開始趨向收斂,波動很平緩;在缺失率為10%時,進化次數達到35次左右開始趨向收斂,波動平緩;在缺失率為20%時,進化次數達到40次左右開始趨向收斂,波動相對較平緩;在缺失率為30%時,進化次數達到50次左右開始趨向收斂,波動不夠平緩.

在相同的缺失率下,采用樣條插值粒子群算法的系統性能函數收斂速度快得多.由圖4f、圖5f、圖6f和圖7f可知,在缺失率為5%時,采用樣條插值法粒子群算法的系統性能指標Pbest只需要20次左右就開始趨于收斂,性能函數的最大值只有110,收斂于6.420 2;采用無插值最小二乘遺傳遞推算法的系統性能指標J需要60次左右才開始趨于收斂,性能函數的最大值接近3 400,收斂于2.291 671 2e＋002.但是隨著缺失率的增大,不管是樣條插值粒子群算法還是無插補最小二乘遺傳遞推算法,插補后系統的性能指標的收斂速度都漸漸降低,性能指標的最大值都越來越大,同時最終收斂的值也隨之變大.如采用樣條插值粒子群算法的系統,在缺失率為10%時,系統的性能指標Pbest的最大值只有130,收斂于26.687 9;在缺失率為20%時,系統的性能指標J的最大值達到190,收斂于66.544 3.

4 結論

通過仿真圖及仿真數據可以直觀地看出,對于不完備信息系統,在相同的數據缺失率下,樣條插值法可以實現較好的插補效果.在缺失率達到較高的數值時,樣條插值法依然能夠實現有價值的插補作用,從而使得該系統的參數得以較準確地估計.但是在缺失率增加到一定值時,由辨識參數的圖形和大小以及系統的性能函數可知,該插補法已無法較精確地模擬完整數據的系統.文章主要討論的是不同缺失率下樣條插值粒子群算法對系統參數的估計問題.樣條插值法是一種應用較為廣泛的插補法,但將其運用在需要辨識參數的系統中還較少.為了進一步探討樣條插值法的適用價值,可以在缺失率更高時采取其他的辨識方法,來實現更精確的參數估計.

[1] 劉黨輝,蔡遠文,蘇永芝,等.系統辨識方法及應用[M].北京:國防工業出版社,2010.

[2] 金勇進,邵軍.缺失數據的統計處理[M].北京:中國統計出版社,2009.

[3] 陳浩,華燈鑫,張毅坤,等.基于三次樣條函數的激光雷達數據可視化插值法[J].儀器儀表學報,2013,34(4):831-837.

[4] R J A Little,D B Rubin.Statistical analysis with missing data[M].U.S.A:Wiley&Sons,2002.

[5] U Nur,N T Longford,J E Cade,et al.The impact of handling missing data on alcohol consumption estimates in the UK wemon cohort study[J].Eur J Epidemiol,2009,24(1):589-595.

[6] H Zhong.The impact of missing data in the estimation of concentration index:a potential source of bias[J].Eur J Heslth Econ,2010,11(1):255-266.

[7] 林東方,宋迎春,金昊.不完全測量數據的EM處理算法[J].大地測量與地球動力學,2011,31(4):112-115.

[8] R Pintelon,J Schoukens.Identification of continuous-time systems with missing data[J].Instrumentation and Measurement,IEEE Transactions on,1999,48(3):736-740.

[9] 譚彥華,李洪興,馬秀娟,等.B樣條函數在模糊系統中的應用[J].控制理論與應用,2013,30(11):1 445-1 456.

[10]張曉丹,邵帥,劉欽圣.基于樣條函數的光滑支持向量機模型[J].北京科技大學學報,2012,34(6):718-725.

[11]高尚,湯可宗,蔣新姿.粒子群優化算法收斂性分析[J].科學技術與工程,2006,6(12):1 625-1 627.

[12]杜大軍,商立立,漆波,等.一種不完全信息下遞推辨識方法及收斂性分析[J].自動化學報,2015,41(8):1 502-1 515.

[13]侯媛彬,汪梅,王立琦.系統辨識及其MATLAB仿真[M].北京:科學出版社,2004.

[14]S Rhode,F Gauterin.Online estimation of vehicle driving resistance parameters with recursive least squares and recursive total least squares[J].Intelligent Vehicles Symposium(IV),2013 IEEE,2013,10(1109):269-276.

Parameter estimation for the incomplete information system based on spline interpolation

LI Zhen1,WEI Li-sheng1?,CHENG Yun-chang2
(1.College of Electrical Engineering,Anhui Polytechnic University,Wuhu 241000,China; 2.State Grid Huoqi County Electric Power Supply Company,Luan 237400,China)

In this issue,a recognition and interpolation method of parameter in an incomplete information system was discussed.Firstly,made the data lost on a certain time by sorting the input signals in the order of time after adding timestamp on each of them.The missing message was supplemented adaptively using spline interpolating,which means the data was supplemented adaptively using linear spline interpolating when one output signal was lost,using the cubic interpolating when two or more output signals were lost.Based on these,the parameter was precisely estimated by particle swarm optimization in the interpolation system.Finally,the precision and effectiveness of the parameter estimation in different missing rates was verified by the method of MATLAB simulation experiment.

incomplete information;spline interpolating;particle swarm optimization;parameter estimation

TP273

A

1672-2477(2015)05-0069-09

2015-09-24

國家自然科學基金資助項目(61203033、61172131、61271377)

李 珍(1989-),女,安徽樅陽人,碩士研究生.

魏利勝(1978-),男,安徽巢湖人,副教授,博士.