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

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

2、ide range of applicati on of system ide ntificati on method and the least squares method elaborated establish a dyn amic system models in System Ide ntificati on applicati ons and examples an alyzed by the least squares method is applied to illustrate the importanee of system iden tificati on.Keywor

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

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

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

6、的目的是至關(guān)重要的。 它對模型類（模型結(jié)構(gòu)）、輸入信號和等價準(zhǔn)則的選擇都有很大的影響。通過辨 識建立數(shù)學(xué)模型通常有四個目的。 估計具有特定物理意義的參數(shù)有些表征系統(tǒng)行為的重要參數(shù)是難以直接測量的，例如在生理、生態(tài)、環(huán)境、經(jīng)濟(jì)等系統(tǒng)中就常有這種情況。這就需要 通過能觀測到的輸入輸出數(shù)據(jù)，用辨識的方法去估計那些參數(shù)。 仿真仿真的核心是要建立一個能模仿真實(shí)系統(tǒng)行為的模型。用于系統(tǒng) 分析的仿真模型要求能真實(shí)反映系統(tǒng)的特性。用于系統(tǒng)設(shè)計的仿真，則強(qiáng)調(diào)設(shè)計 參數(shù)能正確地符合它本身的物理意義。 預(yù)測 這是辨識的一個重要應(yīng)用方面，其目的是用迄今為止系統(tǒng)的可測 量的輸入和輸出去預(yù)測系統(tǒng)輸出的未來的演變。例如最常

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

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

9、信號一般總是可以測量 的，而且過程的動態(tài)特性必然表現(xiàn)在這些輸入輸出數(shù)據(jù)中，那么就可以利用輸入輸出數(shù)據(jù)所提供的信息來建立過程的數(shù)學(xué)模型。這種建模方法就稱為系統(tǒng)辨識。把辨識建模稱作“黑箱建模”。2.1.2 最小二乘法系統(tǒng)辨識結(jié)構(gòu)：本文把待辨識的過程看作“黑箱”。只考慮過程的輸入輸出特性，而不強(qiáng)調(diào) 過程的內(nèi)部機(jī)理。圖中，輸入u(k)和輸出z(k)是可以觀測的；G是系統(tǒng)模型，用來描述系統(tǒng)的輸 入輸出特性；N是噪聲模型，v(k)是白噪聲，e(k)是有色噪聲，根據(jù)表示定理： 可以表示為e(k) =N v(k)G(z)B ( z "I)A(z-)N ( z -1)D(z)C(z J)A(z'

10、;) = 1 a# azZ l| a“az 入B(z') = b才bzZ，III tnbnb；C(z-1) = 1 + ch1 + C2Z + cz-n D(z-1H dQ1 + dzz2 + dAbnb2.1.3準(zhǔn)則函數(shù)設(shè)一個隨機(jī)序列'z(k)，kE (1,2,丄)的均值是參數(shù)日的線性函數(shù)： Etz(k)'hT(k戶，其中h(k)是可測的數(shù)據(jù)向量，那么利用隨機(jī)序列的一個實(shí)現(xiàn)， 使準(zhǔn)則函數(shù)：L2J(R 二' z(k) -hT(kp 心(式2-2)達(dá)到極小的參數(shù)估計值二稱作二的最小二乘估計。最小二乘格式：z(k)二ht(k戶e(k)， 日為模型參數(shù)向量，ek)為零

11、均值隨機(jī)2.2廣義最小二乘法2.2.1廣義最小二乘數(shù)學(xué)模型A(z4)z(kB(zJ)u(kr-1v(k)C(z )式中，u(k)和z(k)表示系統(tǒng)的輸入輸出；v(k)是均值為零的不相關(guān)的隨機(jī) 序列；且'a(z)=1 +aZ +a2Z，* +anaz"aBCz'Hdz'+bzz'i+bnbzH'C(z') =1 +C1Z +C2Z " +Cnczf222廣義最小二乘遞推算法如下攻k)=畝 _1) + K f (k)z,k) hf(k)坯k 1)K ©) = Pf(k-1)hf(k)h?(k)Pf(k-1)hf(k)+1

12、卩 f (k)=l - K f(k)h«k)Pf(k-1):? (k) "(k -1) K e(k)g(k) - he(k)玖 k-1)K e(k) = Pe(k-1) h e(k) h e(k) P°(k - 1) h e(k)1丄Pe(k)珂I -K e(k)h e(k) Pe(k_1)式中hf(k)=乙(k-1),，-乙(k rOu(k 1),，q(k n)F he(k)=Gk-1),-§(k-nc)re(k)=z(k)- h Yk)k)2.2.3廣義最小二乘遞推算法的計算步驟：1.給定初始條件?(0)=；(充分小的實(shí)向量)Pf(0) = a2l

13、(a為充分大的數(shù))?e(0)= 0fe(O)=l2利用式乙(kC(zA)z(k) (k)二 C(z 丄)z(k)計算 Zf(k)和 u«k)；3利用式滬q,，比力,bjhf(k)二-乙(k-1),廠z(k rQ,Uf(k1),，u© Ob)T構(gòu)造hf(k)；%k)二珂k -1) K f (k)zf(k) - hf(k)珂k - 1)4利用式 K ,k)二 Pf(k -1)hf(k)h；(k)Pf (k -1)h,k) 1戸遞推計算?(k)；P f(k)-I -K f(k)h ;(k) P f (k - 1)5 利用(k)二 z(k) - h (k)訛k)和h (k)十 z(

14、k-1), ，-z(k- ra),u(k-1), ,u(k-n)計算鍬);6根據(jù) he(k) =e(k 1),，e(k - nc)廣來構(gòu)造 he(k)；綬(k) = (?(k -1) + K e(k)&k) - hHk)8?(k - 1)7 利用 K e(k)二 Pe(k-1)h e(k) 'h e(k) Pg(1)h g(k)什Pe(k)=l - K e(k) hj(k) R(k-1)返回第2步進(jìn)行迭代計算，直至獲得滿意的辨識結(jié)果3工程實(shí)例3.1典型系統(tǒng)建模以某微循環(huán)流體系統(tǒng)模型的參數(shù)辨識為例我們已經(jīng)得到該系統(tǒng)模型的差分 方程形式，取特定點(diǎn)的壓力波作為模型的輸入，以另一點(diǎn)的壓

15、力波作為模型的輸 出由于我們采集的數(shù)據(jù)是實(shí)時的，因此采用在線辨識方法.由于建立的微循環(huán)流 體系統(tǒng)模型是一個單輸入、單輸出的模型，為使參數(shù)估計的結(jié)果很好地跟蹤參數(shù) 真值的變化，我們采用漸消記憶的最小二乘法對系統(tǒng)模型參數(shù)進(jìn)行辨識，即強(qiáng)調(diào)新數(shù)據(jù)的作用，貶低老數(shù)據(jù)的作用抽象出的SISO系統(tǒng)的差分方程為：z(k) Qz(k - 1) a2z(k - 2) =bi(k - 1) b2u(k - 2)(k) 式( 1-1)參數(shù)取真值為：J二376 0.483 0.57 0.42】，利用MATLAB的M語言辨識 系統(tǒng)中的未知參數(shù)3、a2、“、b2。要求：用參數(shù)的真值利用差分方程求出z(k) 作為測量值，：(k

16、)是均值為0，方差為0.1、0.5和0.01的不相關(guān)隨機(jī)序列。 使用最小二乘算法辨識。3.2廣義最小二乘遞推算法的 MATLA仿真(程序源代碼見附錄)考慮仿真對象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ī)序列。輸入信 號采用4階M序列，幅度為1。選擇如下形式的辨識模型圖3.1廣義最小二乘法辨識實(shí)例結(jié)構(gòu)圖其中取 c1=0,c2=0.4結(jié)果分析及算法優(yōu)化由于辨識算法中輸入或噪聲信號為不相關(guān)隨機(jī)序列，所以每次辨識結(jié)果都不 完全相同。但是，在相同輸入、相同的噪聲、

17、相同的步長條件下，精度大體相同。算法優(yōu)化方案：（1）使用M序列（具有近似白噪聲的性質(zhì)）為輸入信號；（2）增加數(shù)據(jù)長度去L；（3）減小噪聲信號v（k）的方差。4.1廣義最小二乘遞推算法的的 MATLA仿真結(jié)果及分析（1）、輸入選用題目給出的30個隨機(jī)數(shù)，即數(shù)據(jù)長度去L=30，噪聲選用均值 零，方差分別為0.5、0.1和0.01的隨機(jī)序列，辨識結(jié)果如表表 4-1。表中給出了三種情況下辨識參數(shù)結(jié)果即表中的估計值，估計值與真值的相對誤差表4-1真值噪聲方差為0.5噪聲方差為0.1噪聲方差為0.01估計值相對誤差估計值相對誤差估計值相對誤差a11.3761.5534-2.01051.47300.1231

18、1.37230.0063a20.4830.6884-2.00540.58920.10610.48150.0231b10.570.7916-2.02160.67900.18210.57990.0510b20.420.6213-2.00130.53610.11610.42180.0234（2）、輸入均采用M序列，噪聲選擇均值為零，方差為0.5、0.1和0.01 的隨機(jī)序列，辨識步長均為300步，辨識結(jié)果如表4-2。表中給出了三種情況下辨識參數(shù)結(jié)果即表中的估計值，估計值與真值的相對誤差.表4-2真值噪聲方差為0.5噪聲方差為0.1噪聲方差為0.01估計值相對誤差估計值相對誤差估計值相對誤差a11.3

19、761.3660-0.02801.37500.00791.37290.0005a20.4830.4749-0.07010.48590.02920.4710-0.0068b10.570.6613-0.12490.58270.05820.57200.0051b20.420.4112-0.08230.43440.09830.4183-0.0049(3) 數(shù)據(jù)結(jié)果分析：輸入采用 M序列比采用隨機(jī)序列得到的辨識效果更 好。噪聲均值相等時，方差越大，辨識效果越差，反之，方差越小辨識效果越 好?？梢酝ㄟ^增加步長的方法提高辨識精度。下面給出以M序列作為輸入，噪聲均值為零，方差為 0.01的隨機(jī)序列，數(shù) 據(jù)長度

20、取L=30,得到的變化曲線圖：下面給出以M序列作為輸入，噪聲均值為零，方差為 0.01的不相關(guān)隨機(jī)序 列，數(shù)據(jù)長度去L=300,得到的變化曲線圖：File Edi t View Insert Toob Desktop Window HelpD Q S k貝Q的瓷盟h 螯皺變化曲線1. i戶L-L11 1 1 05Q100150200290300參考文獻(xiàn)1 李言俊，張科，系統(tǒng)辨識理論及應(yīng)用，國防工業(yè)出版社，2006年2 方崇智，蕭德云，過程辨識，清華大學(xué)出版社，2002年3 賈秋玲，袁冬莉，欒云鳳，基于 MATLAB7.x/Simulink/Stateflow 系統(tǒng)仿真、分析及設(shè)計，西北工業(yè)大學(xué)

21、出版社，2006年侯媛彬,汪梅,王立琦,系統(tǒng)辨識及其MATLAB真,科學(xué)出版社,2004年附錄廣義最小二乘遞推算法的MATLA仿真程序源代碼：clear %清理工作間變量L=300; % M序列的周期 y1=1;y2=1;y3=1;y4=0; %四個移位寄存器的輸出初始值for i=1:L;% 開始循環(huán)，長度為 Lx1=xor(y3,y4); %第一個移位寄存器的輸入是第三個與第四個移位寄存器的輸出的“或”x2=y1; %第二個移位寄存器的輸入是第一個移位寄存器的輸出x3=y2; %第三個移位寄存器的輸入是第二個移位寄存器的輸出x4=y3; %第四個移位寄存器的輸入是第三個移位寄存器的輸出y(

22、i)=y4; %取出第四個移位寄存器的幅值為"0"和"1"的輸出信號，即 M序列if y(i)>0.5,u(i)=-1; %如果M序列的值為"1",辨識的輸入信號取“ -1”else u(i)=1; % 如果M序列的值為"0",辨識的輸入信號取“ 1”end %小循環(huán)結(jié)束y1=x1;y2=x2;y3=x3;y4=x4; %為下一次的輸入信號做準(zhǔn)備end %大循環(huán)結(jié)束，產(chǎn)生輸入信號ufigure(1); % 第一個圖形stem(u),grid on %顯示出輸入信號 M序列徑線圖并給圖形加上網(wǎng)格v=normrn

23、d(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,即取 c仁0,c2=0endz(2)=0;z(1)=0; %設(shè)z的前兩個初始值為零for k=3:301; %循環(huán)變量從 3到301z(k)=-1.376*z(k-1)-0.483*z(k-2)+57*u(k-1)+0.42*u(k-2)+ze(k-1); %輸出采樣信號(測量值)end%RGL曠義最小二乘辨識c0=0.0001 0.00

24、01 0.0001 0.0001' %直接給出被辨識參數(shù)的初始值,即一個充分小的實(shí)向量pf0=10A6*eye(4,4); %直接給出初始狀態(tài)P0,即一個充分大的實(shí)數(shù)單位矩陣ce0=0.001 0.001'pe0=eye(2,2);c=c0,zeros(4,299); %被辨識參數(shù)矩陣的初始值及大小ce=ce0,zeros(2,299);e=zeros(4,300); %相對誤差的初始值及大小ee=zeros(2,300);s=0;%廣義最小二乘遞推算法的計算步驟for k=3:300;zf(k)=z(k)+ce(1,k-2)*z(k-1)+ce(2,k-2)*z(k-2);

25、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*x1;% 求出 K 的值d仁 zf(k)-hf1'*c0; c仁 c0+k1*d1; %求被辨識參數(shù) ce仁c1-c0; %求參數(shù)當(dāng)前值與上一次的值的差值e2=e1./c0; %求參數(shù)的相對變化e(:,k)=e2; %把當(dāng)前相對變化的列向量加入誤差矩陣的最后一列c0=c1; %新獲得的參數(shù)作為下一次遞推的舊參

26、數(shù)c(:,k)=c1; %把辨識參數(shù)c列向量加入辨識參數(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'*1.642 0.715 0.39 0.35')A2;%求準(zhǔn)則函數(shù)ee(k)=z(k)-h1'*c1;he1=-ee(k-1),-ee(k-2)'x=he1'*pe0*he1+1; x1=i nv(x);k仁 pe0*he1*x1;d仁 ee(k)-he1'*ce0;ce仁 ce0+

27、k1*d1;pe1=pe0-k1*he1'*pe0;ce0=ce1;ce(:,k)=ce1;pe0=pe1;end %大循環(huán)結(jié)束c%辨識參數(shù)變化矩陣%顯示被辨識參數(shù)及其誤差(收斂)情況%分離參數(shù)a1=c(1,1:300); a2=c(2,1:300);b 仁c(3,1:300);b2=c(4,1:300);c仁 ce(1,1:300);c2=ce(2,1:300);ea仁 e(1,1:300); ea2=e(2,1:300);eb 仁e(3,1:300); eb2=e(4,1:300); figure(2); %第二個圖形畫出al，b1，b2,的各次i=1:300; % 橫坐標(biāo)從1到3

28、00Plot(i,a1,'r',i,a2,'k',i,b1,'b',i,b2,'c',i,c1,'b',i,c2,'r') %a2,b1，b2,c1,c2 的各次辨識結(jié)果title(' 參數(shù)變化曲線')%圖形標(biāo)題figure(3); %第三個圖形i=1:300; % 橫坐標(biāo)從1到300 plot(i,ea1,'r',i,ea2,'k:',i,eb1,'b',i,eb2,'k:')%畫出 al，a2,辨識結(jié)果的收斂情況ti

29、tle('誤差曲線')%圖形標(biāo)題考慮仿真對象 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)。輸入信號采用4階M序列。選擇如下形式的辨識模型z(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+v(k)，試?yán)眠f推最小二乘法辨識參數(shù)al、a2、bl、b2。Np=15;r=4;X1=1;X2=1;X3=1;X4=1;men gth = r*Np;a=1;for i=1:1: m_len gthY4=X4;Y3=X3;Y2=X2;Y 仁X1;X4=Y3;X3=Y

30、2;X2=Y1;X仁 xor(Y3,Y4);if Y4=0M(i)=-a;elseM(i)=a;endendfigure;i=1:1:m_le ngth;plot(i,M);%白噪聲no ise = zeros(1,m_le ngth);for i=1:1: m_len gthtemp = no ise + 0.5*ra nds(1, m_len gth);no ise = temp;endno ise = no ise/12;%no ise = temp;% parameter of systemn=2;d=1;a 1= -1;a2=0.5;b1=1;b2=0.5;S_U0=0.2;S_Y0=0.2;% gen erate u,yuO=o nes(1,me ngth)*S_UO;U = M + u0 + noise;figure;i=1:1:m_le ngth;plot(i,U);%formulatio ny(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