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

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

2、s a very wide 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 leastsquares method is applied to illustrate the importance of system identification.keyword

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

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

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

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

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

8、據(jù)原型是一 類十分重要的問題，最常用的逼近原則是讓實測數(shù)據(jù)和估計數(shù)據(jù)之間的距離平方 和最小，這即是最小二乘法。最小二乘法是一種經(jīng)典的數(shù)據(jù)處理方法。在系統(tǒng)辨 識領域中，最小二乘法是一種得到廣泛應用的估計方法 ，可用于動態(tài)系統(tǒng)，靜 態(tài)系統(tǒng)，線性系統(tǒng)，非線性系統(tǒng)?？捎糜陔x線估計，也可用于在線估計。這種 辨識方法主要用于在線辨識。在隨機的環(huán)境下，利用最小二乘法時，并不要求觀 測數(shù)據(jù)提供其概率統(tǒng)計方面的信息，而其估計結果，卻有相當好的統(tǒng)計特性。matlab是一套高性能數(shù)字計算和可視化軟件，它集成概念設計，算法開發(fā)，建模仿真,實時實現(xiàn)于一體，構成了一個使用方便、界面友好的用戶環(huán)境，其強大的擴展功能為各領域

9、的應用提供了基礎。對于比較復雜白生產(chǎn)過程，由于過程的輸入輸出信號一般總是可以測量 的，而且過程的動態(tài)特性必然表現(xiàn)在這些輸入輸出數(shù)據(jù)中，那么就可以利用輸入輸出數(shù)據(jù)所提供的信息來建立過程的數(shù)學模型。這種建模方法就稱為系統(tǒng)辨識。把辨識建模稱作“黑箱建?！?。2.1.2 最小二乘法系統(tǒng)辨識結構：本文把待辨識的過程看作“黑箱”。只考慮過程的輸入輸出特性，而不強調(diào) 過程的內(nèi)部機理。圖中，輸入u(k)和輸出z(k)是可以觀測的；g是系統(tǒng)模型，用來描述系統(tǒng)的輸 入輸出特性；n是噪聲模型，v(k)是白噪聲，e(k)是有色噪聲，根據(jù)表示定理： 可以表示為e(k) =n v(k)1、1g ( z 1) b(z j n

10、 ( z 1) d(z ja(z 1)c(z 1)a(z 1) 1 4z1 a2z2 lanaz nab(z1) dz1 b2z2 l 吐112c(z ) 1 cz c?zcnazd(z 1) d逐1 d2z 2dnbzn2.1.3準則函數(shù)設一個隨機序列z(k),k (1,2, ，l)的均值是參數(shù) 的線性函數(shù): ez(k) ht(k),其中h(k)是可測的數(shù)據(jù)向量，那么利用隨機序列的一個實現(xiàn)， 使準則函數(shù)：l2j( )z(k) ht(k)k 1(式 2-2)達到極小的參數(shù)估計值稱作的最小二乘估計。最小二乘格式：z(k) ht(k) e(k),為模型參數(shù)向量，ek為零均值隨機2.2廣義最小二乘法

11、2.2.1 廣義最小二乘數(shù)學模型a(z1)z(k)11b(z1)u(k)-v(k)c(z1)式中,u(k)和z(k)表示系統(tǒng)的輸入輸出；v(k)是均值為零的不相關的隨機 序列；且a(z 1)1a1z 1a2z2a%z nab(z 1)z1 b2z 2bnbz nbc(z 1)1gz 1c2z2cncz nc2.2.2 廣義最小二乘遞推算法如下?(k)?(k 1) kf(k)zf(k) hf(k)?(k 1)1kf(k) pf(k 1)hf(k)hf(k)pf(k 1)hf(k) 1 1pf (k)i kf(k)hf(k)pf(k 1)?e(k)?e(k 1) ke(k)e?(k) he(k)

12、?e(k 1)1ke(k)pe(k 1)he (k) he(k)pe(k 1)he(k) 1 1pe(k) ike(k)he(k)pe(k 1)式中hf (k)zf(k 1), ,zf(kna),uf (k1),uf (knb)he(k) e?(k 1), , e?(k nc)e?(k) z(k) h (k) ?(k)2.2.3 廣義最小二乘遞推算法的計算步驟：?(0)(充分小的實向量)1. 給定初始條件2 利用式zf (k)uf(k)p(0) a2i (a為充分大的數(shù))?e(0) 0pe(0) ic(z 1)z(k)i，計算 zf(k)和 uf(k);c(z1)z(k) 0 a ,a,b,

13、,b3利用式1 na 1nb，構造 hf(k) ；hf (k) zf(k 1), , zf (k na),uf(k 1), ,uf (k nb)?(k)?(k1) kf(k)zf(k)hf(k)?(k1)4利用式 kf(k) pf(k 1)hf(k)hf(k)pf(k 1)h,(k) 11 遞推計算 &k)；pf (k)ikf(k)hf(k)pf(k1)5利用 e(k) z(k) h(k)4k)和h(k) z(k 1), , z(k na),u(k 1), ,u(k nb) 計算 e?(k) ；6根據(jù) he(k) e?(k 1), e?(k nc) 來構造 he(k) ；?(k)7利用ke(k

14、)pe(k)?(k 1) ke(k)砥k) he(k)(k 1)_1pe(k 1)he(k) he(k)pe(k 1)h,(k) 1ike(k)he(k)pe(k 1)返回第2步進行迭代計算，直至獲得滿意的辨識結果3工程實例3.1 典型系統(tǒng)建模以某微循環(huán)流體系統(tǒng)模型的參數(shù)辨識為例.我們已經(jīng)得到該系統(tǒng)模型的差分 方程形式，取特定點的壓力波作為模型的輸入，以另一點的壓力波作為模型的輸 出.由于我們采集的數(shù)據(jù)是實時的，因此采用在線辨識方法.由于建立的微循環(huán)流 體系統(tǒng)模型是一個單輸入、單輸出的模型，為使參數(shù)估計的結果很好地跟蹤參數(shù) 真值的變化，我們采用漸消記憶的最小二乘法對系統(tǒng)模型參數(shù)進行辨識，即強調(diào)

15、新數(shù)據(jù)的作用，貶低老數(shù)據(jù)的作用.抽象出的siso系統(tǒng)的差分方程為：z(k) az(k 1) a?z(k 2) hu(k 1) b?u(k 2)(k) 式(1-1)參數(shù)取真值為：t 1.376 0.483 0.57 0.42 ,利用matla由勺m語言辨識 系統(tǒng)中的未知參數(shù)a1、a2、b1、b2。要求：用參數(shù)的真值利用差分方程求出z(k) 作為測量值，(k)是均值為0,方差為0.1、0.5和0.01的不相關隨機序列。 使用最小二乘算法辨識。3.2 廣義最小二乘遞推算法的 matla的真(程序源代碼見附錄)考慮仿真對象z(k)= -1.376z(k-1)-0.483z(k-2)+0.57u(k-1

16、)+0.42u(k-2)+v(k)式中，v(k)是均值為0,方差為0.01、0.1和0.5的不相關隨機序列。輸入信 號采用4階m序歹1，幅度為1。v(k)u(k)0.57z 1 0.42z 21 1.376 z 1 0.483z 2y(k)z(k)選擇如下形式的辨識模型圖3.1廣義最小二乘法辨識實例結構圖其中取 c1=0,c2=0.4結果分析及算法優(yōu)化由于辨識算法中輸入或噪聲信號為不相關隨機序列，所以每次辨識結果都不 完全相同。但是，在相同輸入、相同的噪聲、相同的步長條件下，精度大體相同。算法優(yōu)化方案：(1)使用m序列(具有近似白噪聲的性質(zhì))為輸入信號；(2)增加數(shù)據(jù)長度去l;(3)減小噪聲信

17、號v(k)的方差。4.1廣義最小二乘遞推算法的的 matla的真結果及分析(1)、輸入選用題目給出的30個隨機數(shù)，即數(shù)據(jù)長度去l=30,噪聲選用均值零,方差分別為0.5、0.1和0.01的隨機序列，辨識結果如表表 4-1。表中給出了三種情況下辨識參數(shù)結果即表中的估計值，估計值與真值的相對誤差表4-1真值噪聲方差為0.5噪聲方差為0.1噪聲方差為0.01傳計值相對誤差傳計值相對誤差情計值相對誤差a11.3761.5534-2.01051.47300.12311.37230.0063a20.4830.6884-2.00540.58920.10610.48150.0231b10.570.7916-2

18、.02160.67900.18210.57990.0510b20.420.6213-2.00130.53610.11610.42180.0234(2)、輸入均采用m序列，噪聲選擇均值為零，方差為0.5、0.1和0.01 的隨機序列，辨識步長均為300步，辨識結果如表4-2。表中給出了三種情況下辨識參數(shù)結果即表中的估計值，估計值與真值的相對誤差.表4-2真值噪聲方差為0.5噪聲方差為0.1噪聲方差為0.01傳計值相對誤差傳計值相對誤差情計值相對誤差a11.3761.3660-0.02801.37500.00791.37290.0005a20.4830.4749-0.07010.48590.029

19、20.4710-0.0068b10.570.6613-0.12490.58270.05820.57200.0051b20.420.4112-0.08230.43440.09830.4183-0.0049(3)數(shù)據(jù)結果分析：輸入采用 m序列比采用隨機序列得到的辨識效果更 好。噪聲均值相等時，方差越大，辨識效果越差，反之，方差越小辨識效果越 好?？梢酝ㄟ^增加步長的方法提高辨識精度。下面給出以m序列作為輸入，噪聲均值為零，方差為 0.01的隨機序列，數(shù) 據(jù)長度取l=30,得到的變化曲線圖：下面給出以m序列作為輸入，噪聲均值為零，方差為 0.01的不相關隨機序 列，數(shù)據(jù)長度去l=300,得到的變化曲線

20、圖：j j參考文獻1 李言俊，張科，系統(tǒng)辨識理論及應用，國防工業(yè)出版社， 2006年2 方崇智，蕭德云，過程辨識，清華大學出版社， 2002年3 賈秋玲，袁冬莉，欒云鳳，基于 matlab7.x/simulink/stateflow 系統(tǒng)仿真、分析及設計，西北工業(yè)大學出版社， 2006年4侯媛彬，汪梅，王立琦，系統(tǒng)辨識及其matlab&真，科學出版社，2004年附錄廣義最小二乘遞推算法的matla防真程序源代碼：clear % 清理工作間變量l=300; % m 序列的周期y1=1;y2=1;y3=1;y4=0; % 四個移位寄存器的輸出初始值for i=1:l;% 開始循環(huán)，長度為 lx1=

21、xor(y3,y4); % 第一個移位寄存器的輸入是第三個與第四個移位寄存器的輸出 的“或”x2=y1; %第二個移位寄存器的輸入是第一個移位寄存器的輸出x3=y2; %第三個移位寄存器的輸入是第二個移位寄存器的輸出x4=y3; %第四個移位寄存器的輸入是第三個移位寄存器的輸出y(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)結束y1=x1;y2=x2;y3=x3;y4=x4; % 為下一次的

22、輸入信號做準備end %大循環(huán)結束，產(chǎn)生輸入信號ufigure(1); % 第一個圖形stem(u),grid on %顯示出輸入信號m序列徑線圖并給圖形加上網(wǎng)格v=normrnd(0, sqrt(0.01), 1, 300);% 均值為零的，方差為 0.01 或 0.5 或 0.1 不相 關的隨機噪聲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=0 end z(2)=0;z(1)=0; % 設 z 的前兩個初始值為零 for k=3:301; % 循環(huán)變量從3 到 301z(

23、k)=-1.376*z(k-1)-0.483*z(k-2)+57*u(k-1)+0.42*u(k-2)+ze(k-1); %輸出采樣信號(測量值)end%rglst義最小二乘辨識c0=0.0001 0.0001 0.0001 0.0001; % 直接給出被辨識參數(shù)的初始值, 即一個充分小的實向量pf0=10a6*eye(4,4); %直接給出初始狀態(tài) p0,即一個充分大的實數(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); % 相對誤差

24、的初始值及大小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);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 的值d1=zf(k)-hf1*c0; c1=c0+k1*d1; %求被辨識參數(shù)ce1=c1-c0; % 求參數(shù)當前值與

25、上一次的值的差值e2=e1./c0; %求參數(shù)的相對變化e(:,k)=e2; %把當前相對變化的列向量加入誤差矩陣的最后一列c0=c1; % 新獲得的參數(shù)作為下一次遞推的舊參數(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)人2;%求準則函數(shù)ee(k)=z(k)-h1*c1;he1=-ee(k-1),-ee(k-2);x=he1*pe0*he1+1

26、; 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)結束c%w識參數(shù)變化矩陣%顯示被辨識參數(shù)及其誤差( 收斂 ) 情況%分離參數(shù)a1=c(1,1:300); a2=c(2,1: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);f

27、igure(2); % 第二個圖形i=1:300; % 橫坐標從 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 的各次辨識結果title( 參數(shù)變化曲線) % 圖形標題figure(3); % 第三個圖形i=1:300; % 橫坐標從 1 到 300plot(i,ea1,r,i,ea2,k:,i,eb1,b,i,eb2,k:)%畫出a1， a2， b1， b2, 的各次辨識結果的收斂情況 title( 誤差曲線 ) % 圖形標題考慮仿真對象 z(k)+1.5z(k-1)+0.7z(

28、k-2)=u(k-1)+0.5u(k-2)+v(k),其中v(k)是服從正態(tài)分布的白噪聲n(0,1)。輸入信號采用4階m序歹u。選擇如下形式的辨識模型z(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+v(k),試利用遞推最小二乘法辨識參數(shù)al、a2、bl、b2。np=15;r=4;x1=1;x2=1;x3=1;x4=1;m length = r*np;a=1;for i=1:1:m_lengthy4=x4;y3=x3;y2=x2;y1=x1;x4=y3;x3=y2;x2=y1;x1=xor(y3,y4);if y4=0m(i)=-a;elsem(i)=a;end

29、endfigure;i=1:1:m_length;plot(i,m);%白噪聲noise = zeros(1,m length);for i=1:1:m lengthtemp = 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,y u0=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 y