系統(tǒng)辨識(shí)—最小二乘法

1、最小二乘法參數(shù)辨識(shí) 201403027 摘要:系統(tǒng)辨識(shí)在工程中的應(yīng)用非常廣泛,系統(tǒng)辨識(shí)的方法有很多種,最小二乘法是一種應(yīng)用極其廣泛的系統(tǒng)辨識(shí)方法.闡述了動(dòng)態(tài)系統(tǒng)模型的建立及其最小二乘法在系統(tǒng)辨識(shí)中的應(yīng)用,并通過實(shí)例分析說明了最小二乘法應(yīng)用于系統(tǒng)辨識(shí)中的重要意義.關(guān)鍵詞:最小二乘法;系統(tǒng)辨識(shí);動(dòng)態(tài)系統(tǒng)Abstract: System identification in engineering is widely used, system identification methods there are many ways， least squares method is a very wide

2、range of application of system identification method and the least squares method elaborated establish a dynamic system models in System Identification applications and examples analyzed by the least squares method is applied to illustrate the importance of system identification.Keywords: Least Squa

3、res; system identification; dynamic system引言隨著科學(xué)技術(shù)的不斷發(fā)展,人們認(rèn)識(shí)自然、利用自然的能力越來越強(qiáng),對(duì)于未知對(duì)象的探索也越來越深入.我們所研究的對(duì)象,可以依據(jù)對(duì)其了解的程度分為三種類型:白箱、灰箱和黑箱.如果我們對(duì)于研究對(duì)象的內(nèi)部結(jié)構(gòu)、內(nèi)部機(jī)制了解很深入的話,這樣的研究對(duì)象通常稱之為“白箱”;而有的研究對(duì)象,我們對(duì)于其內(nèi)部結(jié)構(gòu)、機(jī)制只了解一部分,對(duì)于其內(nèi)部運(yùn)行規(guī)律并不十分清楚,這樣的研究對(duì)象通常稱之為“灰箱”;如果我們對(duì)于研究對(duì)象的內(nèi)部結(jié)構(gòu)、內(nèi)部機(jī)制及運(yùn)行規(guī)律均一無所知的話,則把這樣的研究對(duì)象稱之為“黑箱”.研究灰箱和黑箱時(shí),將研究的對(duì)象看作是

4、一個(gè)系統(tǒng),通過建立該系統(tǒng)的模型,對(duì)模型參數(shù)進(jìn)行辨識(shí)來確定該系統(tǒng)的運(yùn)行規(guī)律.對(duì)于動(dòng)態(tài)系統(tǒng)辨識(shí)的方法有很多,但其中應(yīng)用最廣泛,辨識(shí)效果良好的就是最小二乘辨識(shí)方法,研究最小二乘法在系統(tǒng)辨識(shí)中的應(yīng)用具有現(xiàn)實(shí)的、廣泛的意義.11 系統(tǒng)辨識(shí)簡介系統(tǒng)辨識(shí)是根據(jù)系統(tǒng)的輸入輸出時(shí)間函數(shù)來確定描述系統(tǒng)行為的數(shù)學(xué)模型。現(xiàn)代控制理論中的一個(gè)分支。通過辨識(shí)建立數(shù)學(xué)模型的目的是估計(jì)表征系統(tǒng)行為的重要參數(shù)，建立一個(gè)能模仿真實(shí)系統(tǒng)行為的模型，用當(dāng)前可測量的系統(tǒng)的輸入和輸出預(yù)測系統(tǒng)輸出的未來演變，以及設(shè)計(jì)控制器。對(duì)系統(tǒng)進(jìn)行分析的主要問題是根據(jù)輸入時(shí)間函數(shù)和系統(tǒng)的特性來確定輸出信號(hào)。對(duì)系統(tǒng)進(jìn)行控制的主要問題是根據(jù)系統(tǒng)的特性設(shè)計(jì)控

5、制輸入，使輸出滿足預(yù)先規(guī)定的要求。而系統(tǒng)辨識(shí)所研究的問題恰好是這些問題的逆問題。通常，預(yù)先給定一個(gè)模型類=M（即給定一類已知結(jié)構(gòu)的模型），一類輸入信號(hào)u和等價(jià)準(zhǔn)則J=L(y，yM)(一般情況下，J是誤差函數(shù)，是過程輸出y和模型輸出yM的一個(gè)泛函)；然后選擇使誤差函數(shù)J達(dá)到最小的模型，作為辨識(shí)所要求的結(jié)果。系統(tǒng)辨識(shí)包括兩個(gè)方面：結(jié)構(gòu)辨識(shí)和參數(shù)估計(jì)。在實(shí)際的辨識(shí)過程中，隨著使用的方法不同，結(jié)構(gòu)辨識(shí)和參數(shù)估計(jì)這兩個(gè)方面并不是截然分開的，而是可以交織在一起進(jìn)行的。1.2系統(tǒng)辨識(shí)的目的在提出和解決一個(gè)辨識(shí)問題時(shí)，明確最終使用模型的目的是至關(guān)重要的。它對(duì)模型類（模型結(jié)構(gòu)）、輸入信號(hào)和等價(jià)準(zhǔn)則的選擇都有很大

6、的影響。通過辨識(shí)建立數(shù)學(xué)模型通常有四個(gè)目的。 估計(jì)具有特定物理意義的參數(shù)有些表征系統(tǒng)行為的重要參數(shù)是難以直接測量的，例如在生理、生態(tài)、環(huán)境、經(jīng)濟(jì)等系統(tǒng)中就常有這種情況。這就需要通過能觀測到的輸入輸出數(shù)據(jù)，用辨識(shí)的方法去估計(jì)那些參數(shù)。 仿真仿真的核心是要建立一個(gè)能模仿真實(shí)系統(tǒng)行為的模型。用于系統(tǒng)分析的仿真模型要求能真實(shí)反映系統(tǒng)的特性。用于系統(tǒng)設(shè)計(jì)的仿真，則強(qiáng)調(diào)設(shè)計(jì)參數(shù)能正確地符合它本身的物理意義。 預(yù)測這是辨識(shí)的一個(gè)重要應(yīng)用方面，其目的是用迄今為止系統(tǒng)的可測量的輸入和輸出去預(yù)測系統(tǒng)輸出的未來的演變。例如最常見的氣象預(yù)報(bào)，洪水預(yù)報(bào)，其他如太陽黑子預(yù)報(bào)，市場價(jià)格的預(yù)測，河流污染物含量的預(yù)測等。預(yù)測模

7、型辨識(shí)的等價(jià)準(zhǔn)則主要是使預(yù)測誤差平方和最小。只要預(yù)測誤差小就是好的預(yù)測模型，對(duì)模型的結(jié)構(gòu)及參數(shù)則很少再有其他要求。這時(shí)辨識(shí)的準(zhǔn)則和模型應(yīng)用的目的是一致的，因此可以得到較好的預(yù)測模型。 控制為了設(shè)計(jì)控制系統(tǒng)就需要知道描述系統(tǒng)動(dòng)態(tài)特性的數(shù)學(xué)模型，建立這些模型的目的在于設(shè)計(jì)控制器。建立什么樣的模型合適，取決于設(shè)計(jì)的方法和準(zhǔn)備采用的控制策略。 2最小二乘方法 2.1.1 系統(tǒng)辨識(shí)最小二乘方法簡介對(duì)工程實(shí)踐中測得的數(shù)據(jù)進(jìn)行理論分析，用恰當(dāng)?shù)暮瘮?shù)去模擬數(shù)據(jù)原型是一類十分重要的問題，最常用的逼近原則是讓實(shí)測數(shù)據(jù)和估計(jì)數(shù)據(jù)之間的距離平方和最小，這即是最小二乘法。最小二乘法是一種經(jīng)典的數(shù)據(jù)處理方法。在系統(tǒng)辨識(shí)領(lǐng)

8、域中 ,最小二乘法是一種得到廣泛應(yīng)用的估計(jì)方法 ,可用于動(dòng)態(tài)系統(tǒng) ,靜態(tài)系統(tǒng) , 線性系統(tǒng) ,非線性系統(tǒng)。可用于離線估計(jì)，也可用于在線估計(jì)。這種辨識(shí)方法主要用于在線辨識(shí)。在隨機(jī)的環(huán)境下，利用最小二乘法時(shí)，并不要求觀測數(shù)據(jù)提供其概率統(tǒng)計(jì)方面的信息，而其估計(jì)結(jié)果，卻有相當(dāng)好的統(tǒng)計(jì)特性。MATLAB是一套高性能數(shù)字計(jì)算和可視化軟件 ,它集成概念設(shè)計(jì) ,算法開發(fā) ,建模仿真 ,實(shí)時(shí)實(shí)現(xiàn)于一體 ,構(gòu)成了一個(gè)使用方便、界面友好的用戶環(huán)境 ,其強(qiáng)大的擴(kuò)展功能為各領(lǐng)域的應(yīng)用提供了基礎(chǔ)。對(duì)于比較復(fù)雜的生產(chǎn)過程 ,由于過程的輸入輸出信號(hào)一般總是可以測量的 ,而且過程的動(dòng)態(tài)特性必然表現(xiàn)在這些輸入輸出數(shù)據(jù)中 ,那么就

9、可以利用輸入輸出數(shù)據(jù)所提供的信息來建立過程的數(shù)學(xué)模型。這種建模方法就稱為系統(tǒng)辨識(shí)。把辨識(shí)建模稱作“黑箱建?！?。2.1.2 最小二乘法系統(tǒng)辨識(shí)結(jié)構(gòu)：本文把待辨識(shí)的過程看作“黑箱”。只考慮過程的輸入輸出特性，而不強(qiáng)調(diào)過程的內(nèi)部機(jī)理。+e(k) 圖2-1 SISO系統(tǒng)辨識(shí)“黑箱”結(jié)構(gòu)圖y(k)u(k)z(k)v(k)圖中，輸入u(k)和輸出z(k)是可以觀測的；G是系統(tǒng)模型，用來描述系統(tǒng)的輸入輸出特性；N是噪聲模型，v(k)是白噪聲，e(k)是有色噪聲，根據(jù)表示定理：可以表示為e(k) =N v(k) 2.1.3 準(zhǔn)則函數(shù) 設(shè)一個(gè)隨機(jī)序列的均值是參數(shù)的線性函數(shù)： ，其中是可測的數(shù)據(jù)向量，那么利用隨機(jī)

10、序列的一個(gè)實(shí)現(xiàn)，使準(zhǔn)則函數(shù)： (式2-2)達(dá)到極小的參數(shù)估計(jì)值稱作的最小二乘估計(jì)。最小二乘格式：，為模型參數(shù)向量，為零均值隨機(jī)2.2 廣義最小二乘法2.2.1 廣義最小二乘數(shù)學(xué)模型式中，u(k)和表示系統(tǒng)的輸入輸出；v(k)是均值為零的不相關(guān)的隨機(jī)序列；且2.2.2 廣義最小二乘遞推算法如下式中2.2.3 廣義最小二乘遞推算法的計(jì)算步驟：1.給定初始條件 2利用式，計(jì)算和；3利用式，構(gòu)造；4利用式遞推計(jì)算；5利用和計(jì)算；6根據(jù)來構(gòu)造；7利用返回第2步進(jìn)行迭代計(jì)算，直至獲得滿意的辨識(shí)結(jié)果。3 工程實(shí)例3.1典型系統(tǒng)建模以某微循環(huán)流體系統(tǒng)模型的參數(shù)辨識(shí)為例.我們已經(jīng)得到該系統(tǒng)模型的差分方程形式,取

11、特定點(diǎn)的壓力波作為模型的輸入,以另一點(diǎn)的壓力波作為模型的輸出.由于我們采集的數(shù)據(jù)是實(shí)時(shí)的,因此采用在線辨識(shí)方法.由于建立的微循環(huán)流體系統(tǒng)模型是一個(gè)單輸入、單輸出的模型,為使參數(shù)估計(jì)的結(jié)果很好地跟蹤參數(shù)真值的變化,我們采用漸消記憶的最小二乘法對(duì)系統(tǒng)模型參數(shù)進(jìn)行辨識(shí),即強(qiáng)調(diào)新數(shù)據(jù)的作用,貶低老數(shù)據(jù)的作用.抽象出的SISO系統(tǒng)的差分方程為： 式（1-1）參數(shù)取真值為：，利用MATLAB 的M語言辨識(shí)系統(tǒng)中的未知參數(shù)、。要求：用參數(shù)的真值利用差分方程求出作為測量值，是均值為0，方差為0.1、0.5和0.01的不相關(guān)隨機(jī)序列。使用最小二乘算法辨識(shí)。3.2 廣義最小二乘遞推算法的MATLAB仿真（程序源代

12、碼見附錄）考慮仿真對(duì)象z(k)= -1.376z(k-1)-0.483z(k-2)+0.57u(k-1)+0.42u(k-2)+v(k)式中，v(k) 是均值為0，方差為0.01、0.1和0.5的不相關(guān)隨機(jī)序列。輸入信號(hào)采用4階M序列，幅度為1。選擇如下形式的辨識(shí)模型+y(k)u(k)e(k) z(k)+v(k)圖3.1 廣義最小二乘法辨識(shí)實(shí)例結(jié)構(gòu)圖其中取c1=0,c2=0.4結(jié)果分析及算法優(yōu)化由于辨識(shí)算法中輸入或噪聲信號(hào)為不相關(guān)隨機(jī)序列，所以每次辨識(shí)結(jié)果都不完全相同。但是，在相同輸入、相同的噪聲、相同的步長條件下，精度大體相同。算法優(yōu)化方案：（1）使用M序列（具有近似白噪聲的性質(zhì)）為輸入信號(hào)

13、；（2）增加數(shù)據(jù)長度去L；（3）減小噪聲信號(hào)v(k)的方差。4.1 廣義最小二乘遞推算法的的MATLAB仿真結(jié)果及分析(1)、輸入選用題目給出的30個(gè)隨機(jī)數(shù)，即數(shù)據(jù)長度去L=30，噪聲選用均值零，方差分別為0.5、0.1和0.01的隨機(jī)序列，辨識(shí)結(jié)果如表表4-1。表中給出了三種情況下辨識(shí)參數(shù)結(jié)果即表中的估計(jì)值，估計(jì)值與真值的相對(duì)誤差表4-1真值噪聲方差為0.5噪聲方差為0.1噪聲方差為0.01估計(jì)值相對(duì)誤差估計(jì)值相對(duì)誤差估計(jì)值相對(duì)誤差a11.3761.5534-2.01051.47300.12311.37230.0063a20.4830.6884-2.0054 0.58920.1061 0.4

14、8150.0231b10.570.7916-2.02160.67900.1821 0.57990.0510b20.420.6213-2.0013 0.53610.1161 0.42180.0234（2）、輸入均采用M序列，噪聲選擇均值為零，方差為0.5、0.1和0.01的隨機(jī)序列，辨識(shí)步長均為300步，辨識(shí)結(jié)果如表4-2。表中給出了三種情況下辨識(shí)參數(shù)結(jié)果即表中的估計(jì)值，估計(jì)值與真值的相對(duì)誤差. 表4-2真值噪聲方差為0.5噪聲方差為0.1噪聲方差為0.01估計(jì)值相對(duì)誤差估計(jì)值相對(duì)誤差估計(jì)值相對(duì)誤差a11.3761.3660-0.02801.37500.00791.37290.0005a20.4

15、830.4749-0.0701 0.48590.0292 0.4710-0.0068b10.570.6613-0.12490.58270.0582 0.57200.0051b20.420.4112-0.0823 0.43440.0983 0.4183-0.0049（3）數(shù)據(jù)結(jié)果分析：輸入采用M序列比采用隨機(jī)序列得到的辨識(shí)效果更好。噪聲均值相等時(shí)，方差越大，辨識(shí)效果越差，反之，方差越小辨識(shí)效果越好。可以通過增加步長的方法提高辨識(shí)精度。 下面給出以M序列作為輸入，噪聲均值為零，方差為0.01的隨機(jī)序列，數(shù)據(jù)長度取L=30，得到的變化曲線圖： 下面給出以M序列作為輸入，噪聲均值為零，方差為0.01的

16、不相關(guān)隨機(jī)序列，數(shù)據(jù)長度去L=300，得到的變化曲線圖：參考文獻(xiàn)1 李言俊，張科，系統(tǒng)辨識(shí)理論及應(yīng)用，國防工業(yè)出版社，2006年2方崇智，蕭德云，過程辨識(shí)，清華大學(xué)出版社，2002年3賈秋玲，袁冬莉，欒云鳳，基于MATLAB7.x/Simulink/Stateflow系統(tǒng)仿真、分析及設(shè)計(jì)，西北工業(yè)大學(xué)出版社，2006年4 侯媛彬, 汪梅, 王立琦,系統(tǒng)辨識(shí)及其MATLAB仿真,科學(xué)出版社,2004年附錄廣義最小二乘遞推算法的MATLAB仿真程序源代碼：clear %清理工作間變量L=300; % M序列的周期y1=1;y2=1;y3=1;y4=0; %四個(gè)移位寄存器的輸出初始值for i=1:

17、L;%開始循環(huán)，長度為L x1=xor(y3,y4); %第一個(gè)移位寄存器的輸入是第三個(gè)與第四個(gè)移位寄存器的輸出的“或” x2=y1; %第二個(gè)移位寄存器的輸入是第一個(gè)移位寄存器的輸出 x3=y2; %第三個(gè)移位寄存器的輸入是第二個(gè)移位寄存器的輸出 x4=y3; %第四個(gè)移位寄存器的輸入是第三個(gè)移位寄存器的輸出 y(i)=y4; %取出第四個(gè)移位寄存器的幅值為"0"和"1"的輸出信號(hào)，即M序列 if y(i)>0.5,u(i)=-1; %如果M序列的值為"1", 辨識(shí)的輸入信號(hào)取“-1” else u(i)=1; %如果M序列的

18、值為"0", 辨識(shí)的輸入信號(hào)取“1” end %小循環(huán)結(jié)束 y1=x1;y2=x2;y3=x3;y4=x4; %為下一次的輸入信號(hào)做準(zhǔn)備end %大循環(huán)結(jié)束，產(chǎn)生輸入信號(hào)ufigure(1); %第一個(gè)圖形stem(u),grid on %顯示出輸入信號(hào)M序列徑線圖并給圖形加上網(wǎng)格v=normrnd(0, sqrt(0.01), 1, 300);%均值為零的，方差為0.01或0.5或0.1不相關(guān)的隨機(jī)噪聲ze(2)=0;ze(1)=0;for k=3:301; ze(k)=0*ze(k-1)+0*ze(k-2)+v(k-1);%C(z1)=1,即取c1=0,c2=0endz

19、(2)=0;z(1)=0; %設(shè)z的前兩個(gè)初始值為零for k=3:301; %循環(huán)變量從3到301 z(k)=-1.376*z(k-1)-0.483*z(k-2)+57*u(k-1)+0.42*u(k-2)+ze(k-1); %輸出采樣信號(hào)（測量值） end%RGLS廣義最小二乘辨識(shí)c0=0.0001 0.0001 0.0001 0.0001' %直接給出被辨識(shí)參數(shù)的初始值,即一個(gè)充分小的實(shí)向量pf0=106*eye(4,4); %直接給出初始狀態(tài)P0，即一個(gè)充分大的實(shí)數(shù)單位矩陣ce0=0.001 0.001'pe0=eye(2,2);c=c0,zeros(4,299); %

20、被辨識(shí)參數(shù)矩陣的初始值及大小ce=ce0,zeros(2,299);e=zeros(4,300); %相對(duì)誤差的初始值及大小ee=zeros(2,300);s=0;%廣義最小二乘遞推算法的計(jì)算步驟 for k=3:300; zf(k)=z(k)+ce(1,k-2)*z(k-1)+ce(2,k-2)*z(k-2); uf(k)=u(k)+ce(1,k-2)*u(k-1)+ce(2,k-2)*u(k-2); hf1=-zf(k-1),-zf(k-2),uf(k-1),uf(k-2)' x=hf1'*pf0*hf1+1; x1=inv(x); %開始求K(k) k1=pf0*hf1*

21、x1;%求出K的值 d1=zf(k)-hf1'*c0; c1=c0+k1*d1; %求被辨識(shí)參數(shù)c e1=c1-c0; %求參數(shù)當(dāng)前值與上一次的值的差值 e2=e1./c0; %求參數(shù)的相對(duì)變化 e(:,k)=e2; %把當(dāng)前相對(duì)變化的列向量加入誤差矩陣的最后一列 c0=c1; %新獲得的參數(shù)作為下一次遞推的舊參數(shù) c(:,k)=c1; %把辨識(shí)參數(shù)c 列向量加入辨識(shí)參數(shù)矩陣的最后一列 pf1=pf0-k1*hf1'*pf0; %求出 p(k)的值 pf0=pf1; %給下次用 h1=-z(k-1),-z(k-2),u(k-1),u(k-2)' s=s+(z(k)-h1

22、'*1.642 0.715 0.39 0.35')2;%求準(zhǔn)則函數(shù) ee(k)=z(k)-h1'*c1; he1=-ee(k-1),-ee(k-2)' x=he1'*pe0*he1+1; x1=inv(x); k1=pe0*he1*x1; d1=ee(k)-he1'*ce0; ce1=ce0+k1*d1; pe1=pe0-k1*he1'*pe0; ce0=ce1; ce(:,k)=ce1; pe0=pe1; end %大循環(huán)結(jié)束c%辨識(shí)參數(shù)變化矩陣%顯示被辨識(shí)參數(shù)及其誤差(收斂)情況%分離參數(shù)a1=c(1,1:300); a2=c(2,1

23、:300);b1=c(3,1:300);b2=c(4,1:300); c1=ce(1,1:300);c2=ce(2,1:300); ea1=e(1,1:300); ea2=e(2,1:300);eb1=e(3,1:300); eb2=e(4,1:300);figure(2); %第二個(gè)圖形i=1:300; %橫坐標(biāo)從1到300plot(i,a1,'r',i,a2,'k',i,b1,'b',i,b2,'c',i,c1,'b',i,c2,'r') %畫出a1，a2,b1，b2,c1,c2的各次辨識(shí)結(jié)果t

24、itle('參數(shù)變化曲線') %圖形標(biāo)題figure(3); %第三個(gè)圖形i=1:300; %橫坐標(biāo)從1到300plot(i,ea1,'r',i,ea2,'k:',i,eb1,'b',i,eb2,'k:') %畫出a1，a2，b1，b2,的各次辨識(shí)結(jié)果的收斂情況title('誤差曲線') %圖形標(biāo)題考慮仿真對(duì)象z(k)+1.5z(k-1)+0.7z(k-2)=u(k-1)+0.5u(k-2)+v(k)，其中v(k)是服從正態(tài)分布的白噪聲N(0,1)。輸入信號(hào)采用4階M序列。選擇如下形式的辨識(shí)模型z(

25、k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+v(k)，試?yán)眠f推最小二乘法辨識(shí)參數(shù)a1、a2、b1、b2。Np=15;r=4;X1=1;X2=1;X3=1;X4=1;m_length = r*Np;a=1;for i=1:1:m_length Y4=X4;Y3=X3;Y2=X2;Y1=X1; X4=Y3;X3=Y2;X2=Y1; X1=xor(Y3,Y4); if Y4=0 M(i)=-a; else M(i)=a; endendfigure;i=1:1:m_length;plot(i,M);% 白噪聲noise = zeros(1,m_length);for

26、 i=1:1:m_length temp = noise + 0.5*rands(1,m_length); noise = temp;endnoise = noise/12;%noise = temp;% parameter of systemn=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5;S_U0=0.2;S_Y0=0.2;% generate u,yu0=ones(1,m_length)*S_U0;U = M + u0 + noise;figure;i=1:1:m_length;plot(i,U);%formulationy(1)=0;y(2)=0;y(3)=0;Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3);for k=4:m_length y(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2); Y(k)=S_Y0+y(k)+noise(k);endfigure;i=1:1:m_length;plot(i,Y);