應(yīng)用最小二乘一次完成法和遞推最小二乘法算法的系統(tǒng)辨識(shí)

1、目錄1 引言 21.1 概述 21.2 辨識(shí)的基本步驟 22 系統(tǒng)辨識(shí)輸入信號(hào)的產(chǎn)生方法和理論依據(jù) 32.1 白噪聲序列 32.1.1 白噪聲序列的產(chǎn)生方法 32.2 M 序列的產(chǎn)生 42.2. 1 偽隨機(jī)噪聲 42.2.2 M 序列的產(chǎn)生方法 43 應(yīng)用經(jīng)典辨識(shí)方法的辨識(shí)方案。 63.1 經(jīng)典辨識(shí)方法概述 63.2 經(jīng)典辨識(shí)方法的實(shí)現(xiàn) 64 最小二乘法的理論基礎(chǔ) 74.1 最小二乘法 74.1.1 最小二乘法估計(jì)中的輸入信號(hào) 94.1.2 最小二乘估計(jì)的概率性質(zhì) 94.2 遞推最小二乘法 105 兩種算法的實(shí)現(xiàn)方案 115.1 最小二乘法一次完成算法實(shí)現(xiàn) 115.1.1 最小二乘一次完成算法

2、程序框圖 115.1.2 一次完成法程序 115.1.3 一次完成算法程序運(yùn)行結(jié)果 115.1.4 辨識(shí)數(shù)據(jù)比較 125.1.5 程序運(yùn)行曲線 125.2 遞推最小二乘法的實(shí)現(xiàn) 125.2.1 遞推算法實(shí)現(xiàn)步驟 125.2.2 程序編制思路： 135.2.3 遞推最小二乘法程序框圖 145.2.4 程序運(yùn)行曲線 155.2.5 測(cè)試結(jié)果 165.2.6 遞推數(shù)據(jù)表 166 結(jié)論 167 參考文獻(xiàn) 178 附錄 176應(yīng)用最小二乘一次完成法和遞推最小二乘法算法的系統(tǒng)辨識(shí)摘要 ：本題針對(duì)一個(gè)單輸入單輸出系統(tǒng)的便是問題，辨識(shí)的輸入信號(hào)采用的是偽隨機(jī)二位式序列（ M 序列），系統(tǒng)噪聲為獨(dú)立 同分布高斯

3、隨機(jī)向量序列 （白噪聲），辨識(shí)的算法是遞推最小二乘法和廣義最小二乘法， 本文簡(jiǎn)單描述應(yīng)用經(jīng)典辨識(shí)方法的辨識(shí) 方案，詳細(xì)描述了輸入信號(hào)、噪聲的產(chǎn)生方法及 matlab 程序，闡述了用兩種不同算法的辨識(shí)原理并對(duì)它們的推導(dǎo)過程及辨識(shí)程 序編制思路做了詳細(xì)的描述。最后結(jié)合真值與估計(jì)值對(duì)不同辨識(shí)算法的優(yōu)劣進(jìn)行了比較。關(guān)鍵詞 ：系統(tǒng)辨識(shí) M 序列 最小二乘法1 引言1.1 概述系統(tǒng)辨識(shí)是現(xiàn)代控制理論中的一個(gè)分支， 它是根據(jù)系統(tǒng)的輸入輸出時(shí)間函數(shù)來確定描述系統(tǒng)行為的數(shù)學(xué)模型。 通過辨識(shí)建立數(shù)學(xué)模型的目的是估計(jì)表征系統(tǒng)行為的重要參數(shù)，建立一個(gè)能模仿真實(shí)系統(tǒng)行為的模型，用當(dāng) 前可測(cè)量的系統(tǒng)的輸入和輸出預(yù)測(cè)系統(tǒng)輸

4、出的未來演變，以及設(shè)計(jì)控制器。對(duì)系統(tǒng)進(jìn)行分析的主要問題是根據(jù)輸入時(shí)間函數(shù)和系統(tǒng)的特性來確定輸出信號(hào)。對(duì)系統(tǒng)進(jìn)行控制的主要問題 是根據(jù)系統(tǒng)的特性設(shè)計(jì)控制輸入，使輸出滿足預(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è)方面并不是截然分開的，而是

5、可以交織在一起進(jìn)行 的。1.2 辨識(shí)的基本步驟 先驗(yàn)知識(shí)和建模目的的依據(jù)。先驗(yàn)知識(shí)指關(guān)于系統(tǒng)運(yùn)動(dòng)規(guī)律、數(shù)據(jù)以及其他方面的已有知識(shí)。這些知識(shí)對(duì) 選擇模型結(jié)構(gòu)、設(shè)計(jì)實(shí)驗(yàn)和決定辨識(shí)方法等都有重要作用。用于不同目的的模型可能會(huì)有很大差別。 實(shí)驗(yàn)設(shè)計(jì)。辨識(shí)是從實(shí)驗(yàn)數(shù)據(jù)中提取有關(guān)系統(tǒng)信息的過程，設(shè)計(jì)實(shí)驗(yàn)的目標(biāo)之一是要使所得到的數(shù)據(jù)能包 含系統(tǒng)更多的信息。主要包括輸入信號(hào)設(shè)計(jì)，采樣區(qū)間設(shè)計(jì)，預(yù)采樣濾波器設(shè)計(jì)等。 結(jié)構(gòu)辨識(shí)。即選擇模型類中的數(shù)學(xué)模型M的具體表達(dá)形式。除線性系統(tǒng)的結(jié)構(gòu)可通過輸入輸出數(shù)據(jù)進(jìn)行辨識(shí)外 ，一般的模型結(jié)構(gòu)主要通過先驗(yàn)知識(shí)獲得。 參數(shù)估計(jì)。知道模型的結(jié)構(gòu)后，用輸入輸出數(shù)據(jù)確定模型中的未知參

6、數(shù)。實(shí)際測(cè)量都是有誤差的，所以參 數(shù)估計(jì)以統(tǒng)計(jì)方法為主。 模型適用性檢驗(yàn)。造成模型不適用主要有三方面原因：模型結(jié)構(gòu)選擇不當(dāng)；實(shí)驗(yàn)數(shù)據(jù)誤差過大或數(shù)據(jù)代表性太差；辨識(shí)算法存在問題。檢驗(yàn)方法主要有利用先驗(yàn)知識(shí)檢驗(yàn)和利用數(shù)據(jù)檢驗(yàn)兩類。11.3設(shè)待辨識(shí)系統(tǒng)如圖 1 所示。二階系統(tǒng)為:參數(shù)真值為:1設(shè) f (k)1待辨識(shí)系統(tǒng)圖圖a(z 1) b(z1) f(z1)1 21aza?z1 2b1zb2z1 2f1zf2zai1.5; a20.7;b|1；b20.50; fo1；f11；f20.2z1)1, (k)(k),(k)為白色噪聲，簡(jiǎn)單描述應(yīng)用經(jīng)典辨識(shí)方法的辨識(shí)方案。12 (k)為有色噪聲，f(z )f

7、1z 1f2z11.00; f20.2 ,(k)為獨(dú)立同分布的高斯序列,0.4,2系統(tǒng)辨識(shí)輸入信號(hào)的產(chǎn)生方法和理論依據(jù)2.1白噪聲序列如果隨機(jī)序列 w t 均值為0,并且兩兩不相的關(guān)的，對(duì)應(yīng)自相關(guān)函數(shù)為Rw 12 l,l 0, 1, 2,ggg式中l(wèi)0； 00則稱這種隨機(jī)序列為白噪聲序列 2.1.1白噪聲序列的產(chǎn)生方法F面主要介紹(0, 1)均勻分布和正態(tài)分布隨機(jī)數(shù)的產(chǎn)生方法在計(jì)算機(jī)上產(chǎn)生(0, 1)均勻分布隨機(jī)數(shù)的方法很多，其中最簡(jiǎn)單、最方便的是數(shù)學(xué)方法。產(chǎn)生偽隨機(jī) 數(shù)的數(shù)學(xué)方法很多，其中最常用的是乘同余法和混合同余法。(1)乘同余法這種方法先用遞推同余式產(chǎn)生正整數(shù)序列Xi=Axi-1(mo

8、dM),i=1,2,3式中：M為2的方幕，k為大于2的整數(shù)；A = 3(mod8)或A = 5(mod8),且A不能太小；初值 x0取正奇數(shù)，例如取 x0=1.再令 i 十 1,2,.則是偽隨機(jī)序列，循環(huán)周期可達(dá)i(2)混合同余法混合同余法產(chǎn)生偽隨機(jī)數(shù)的遞推同余式為X Axi 1 c(mod M )式中 M 2k,K為大于2的整數(shù)；A三1(mod4),即A 2n 1，其中n為滿足關(guān)系式2W nW 34的整數(shù)。初Xk值X0為非負(fù)整數(shù)。令 i 丄，則 i是循環(huán)周期為2的偽隨機(jī)數(shù)序列M2.2 M序列的產(chǎn)生在進(jìn)行系統(tǒng)辨識(shí)時(shí)，選用白噪聲作為辨識(shí)輸入信號(hào)可以保證獲得較好的便是效果，但工程上難以實(shí)現(xiàn)。M序列

9、是一種很好的辨識(shí)輸入信號(hào)，它具有近似白燥聲的性質(zhì)，不僅可以保證有較好的辨識(shí)效果，而且工程上 易于實(shí)現(xiàn)。M序列是偽隨機(jī)二位式序列的一種形式。在介紹 M序列之前，先介紹一下偽隨機(jī)噪聲的概念。2.2. 1 偽隨機(jī)噪聲對(duì)白噪聲的一個(gè)樣本函數(shù)w(t)截取0,T時(shí)間內(nèi)一段，對(duì)其它時(shí)間段T,2T,2T,3T,，以 周期T延拖下去，這樣獲得的函數(shù) w(t)是周期T的函數(shù)，在0,T時(shí)間內(nèi)是白噪聲，在此時(shí)間之外是重復(fù)的白噪聲，它的自相關(guān) 函數(shù) RwE w t w t 的周期也是T。由于在0,T時(shí)間內(nèi)自相關(guān)函數(shù)R,就是白噪聲的自相關(guān)函數(shù)，它具有周期性，稱為w(t)為偽隨機(jī)噪聲。2.2.2 M序列的產(chǎn)生方法M序列是一

10、種離散二位式：隨機(jī)序列，所謂二位式”是指每個(gè)隨機(jī)變量只有 2種狀態(tài)?？捎枚嗉?jí)線性反饋 移位寄存器產(chǎn)生 M序列。M序列是最長線性反饋移存器序列的簡(jiǎn)稱，是由帶線性反饋的移存器產(chǎn)生的周期最長的一種序列。具有較強(qiáng)的抗干擾能力和較低的截獲概率，而且長的M序列更容易在一定的強(qiáng)噪聲中被提取，這樣就能夠充分保證數(shù)據(jù)的正常通信。通常產(chǎn)生偽隨機(jī)序列的電路為反饋移存器.一般說來，由n級(jí)移位寄存器產(chǎn)生的周期為N=2?-1的M序列，在一個(gè)循環(huán)周期內(nèi)，“0”出現(xiàn)的次數(shù)為 N 1/2, “1”出現(xiàn)的次數(shù)為 N 1/2。現(xiàn)在我們引入 M序列的本原多 項(xiàng)式的概念。若一個(gè) n次多項(xiàng)式f (x)滿足以下條件(1) f (X)為既約

11、的。(2) f (x)可整除(xm 1), m 2n 1。(3) f (x)除不盡(xq 1), q m則f (x)為本原多項(xiàng)式。一個(gè)4級(jí)M序列可以通過線性反饋移位寄存器產(chǎn)生，如下圖所示：圖2周期為15的偽隨機(jī)序列產(chǎn)生器圖每級(jí)移位寄存器由雙穩(wěn)態(tài)觸發(fā)器和門電路組成，稱為1位，分別以0和1來表示2種狀態(tài)。當(dāng)移位脈沖到來時(shí)，每位的內(nèi)容移至下一位，最后1位移出的內(nèi)容即為輸出。為了保持連續(xù)工作，將最后2級(jí)寄存器的內(nèi)容經(jīng)過適當(dāng)?shù)倪壿嬤\(yùn)算后反饋到第1級(jí)寄存器作為輸入。當(dāng)一個(gè)移位脈沖到來后，第1級(jí)寄存器的內(nèi)容送到第2級(jí)，第2級(jí)寄存器的內(nèi)容送到第3級(jí)，第3級(jí)寄存器的內(nèi)容送到第 4級(jí)，而第3級(jí)和第4級(jí)寄存器的內(nèi)容

12、作模和2相加后再反饋到第 1級(jí)寄存器。產(chǎn)生偽隨機(jī)序列時(shí)要求寄存器的初始狀態(tài)不全為0,因?yàn)槿?初始狀態(tài)將導(dǎo)致各級(jí)寄存器輸出永遠(yuǎn)是0。如果寄存器的初始內(nèi)容都是1,第1個(gè)移位脈沖來到后，4級(jí)寄存器的內(nèi)容變以 0111， 個(gè)周期的變化規(guī)律為：1111 0111 0011 0001 1000 0100 0010 1001 1100 01101011 0101 1010 1101 1110 1111 一個(gè)周期結(jié)束后，產(chǎn)生的15種不同的狀態(tài)。計(jì)算機(jī)模擬產(chǎn)生 M序列非常方便，先定義輸出序列長度和一個(gè)數(shù)組，數(shù)組個(gè)數(shù)等于移位寄存器的個(gè)數(shù)，通過 使用異或指令，再利用for循環(huán)指令，即可完成任意長度和級(jí)數(shù)M序列的產(chǎn)生

13、。圖3M序列產(chǎn)生的流程圖圖4MATLAB 仿真M序列圖形如圖3是產(chǎn)生M級(jí)長度為L的M序列程序流程圖，圖 4所示是Matlab仿真4級(jí)移位寄存器，長度 15的M 序列輸出圖形。23應(yīng)用經(jīng)典辨識(shí)方法的辨識(shí)方案。3.1經(jīng)典辨識(shí)方法概述用正弦信號(hào)、單位階躍信號(hào)和M序列辨識(shí)系統(tǒng)的方法都稱為經(jīng)典辨識(shí)方法。在經(jīng)典辨識(shí)方法中，用得最多的是脈沖響應(yīng)。用 M序列作為輸入信號(hào)，再用相關(guān)法處理測(cè)試結(jié)果，可以很方便地得到系統(tǒng)的脈沖響應(yīng)。 用白噪聲作為輸入信號(hào)確定系統(tǒng)脈沖響應(yīng)的方法，稱為相關(guān)法。相關(guān)法，就是根據(jù)維納-霍夫積分方程，禾U 用輸入信號(hào)的自相關(guān)函數(shù)和輸入與輸出的互相關(guān)函數(shù)確定系統(tǒng)脈沖響應(yīng)的方法。采用周期性的偽隨

14、機(jī)信號(hào)作 為輸入信號(hào)可以使計(jì)算變得簡(jiǎn)單，所以用M序列辨識(shí)系統(tǒng)脈沖響應(yīng)是用得最廣泛的一種方法。3.2經(jīng)典辨識(shí)方法的實(shí)現(xiàn)由系統(tǒng)辨識(shí)圖可得如下關(guān)系統(tǒng)：y(k)= x(k)+ (k) g(k) (k)(3.2.1)其中：（k）為輸入白噪聲信號(hào)，x（k）為實(shí)際輸入信號(hào)，y（k）為實(shí)際輸出信號(hào)，g（ k ）為輸入信號(hào)為u（ k ）時(shí)的 脈沖響應(yīng)。根據(jù)維納-霍夫積分方程可得離散的維納霍夫積分方程：N 1k)(322)g( k)Rx(k 0令式中Ruxg(k)Rx(k)(323)再設(shè):RuxRux(0)Rux(1)Rux(2)Ru(0)Ru(1)Ru( 1)Ru(0)Ru( N 1)Ru( N 2)Rux(N

15、1)Ru(N 1)Ru(N 2)Ru(0)k則上式可寫成以下形式:再根據(jù)上式（3.2.3）可得:RuxRgg 1R 1Rux(3.2.4)由二電平M序列的自相關(guān)函數(shù)的計(jì)算公式可以求得112 R aN1N11N1NR 1111NN根據(jù)遞推公式求 m次觀測(cè)的脈沖響應(yīng)2 11N1 21(3.2.5)a N 11 12g(m)對(duì)系統(tǒng)進(jìn)行m次觀測(cè)，m由m次觀測(cè)得到的Rux()用Rux( ,m )來表示，則有如下遞推式:Rux(m,m)=m+ 1 k=oy(x) x(k- m)=Rux(m,m-1)+1m+ 1y(m) x(m- m) - Rux(m,m- 1)(326)綜合以上公式就可以得到以下遞推公式

16、：211x(m)1N121/ 、x(m- y(m)1)gm 1gm - gm-1 +.2m+ 1 a N+ 1 D112x(m-N+ 1)根據(jù)遞推式，可從 gm 1及新的觀測(cè)數(shù)據(jù)得到gm，隨著觀測(cè)數(shù)據(jù)的增加,(327)gm的精確度不斷的提高。34最小二乘法的理論基礎(chǔ)4.1最小二乘法設(shè)單輸入單輸出線性定長系統(tǒng)的差分方程表示為:y kak 1a2y k 2Lany k nb)u kb|Uk 1Lbnu knk其中nk n kk in(k)均值為0的白噪聲，現(xiàn)分另iJ測(cè)出n+N個(gè)輸出輸入值y(1),y(2), ,y( n+N),u(1),u(2 ,iu(n+N),則可寫出N個(gè)方程，寫成向量一矩陣形式

17、y n 1y nLy 1u n 1Lu1y n 2y n 1Ly 2u n 2Lu2MMMMy n Nyn N 1 Ly N u n NLuNa1Mn 1(4.1.1)ann 2boMMn Nbn7yn 1Mn1yyn 2ann2M，b°JMyn NMnNbny nLy1un 1Lu 1y n 1Ly2un 2Lu 2MMMyn N 1LyNun NLu N則式(4.1.1)可寫為y(4.1.2)式中：y為N維輸出向量；E為N為維噪聲向量；0為(2n+1)維參數(shù)向量；為N X(2n+1)測(cè)量矩陣。ai因此，式(4.1.1)是一個(gè)含有(2n+1)個(gè)未知參數(shù)，由 N個(gè)方程組成的聯(lián)立方程組

18、。在給定輸出向量y和測(cè)量矩陣的條件下求參數(shù)B的估計(jì)，這就是系統(tǒng)辨識(shí)問題。 設(shè) ？表示的估計(jì)值，？表示y的最優(yōu)估計(jì)，則有?式中：(4.1.3)? n? nMy nMb0Mbn設(shè) e(k)=y(k)- ?(k), e(k)稱為殘差，則有e=y- ?=y-最小二乘估計(jì)要求殘差的平方和最小，即按照指數(shù)函數(shù)Te eJ求j對(duì) ？的偏導(dǎo)數(shù)并令其等于 0可得：(4.1.4)由式(4.1.5)可得的最小二乘估計(jì):J為極小值的充分條件是:即矩陣T為正定矩陣，或者說是非奇異的。(4.1.5)(4.1.6)094.1.1最小二乘法估計(jì)中的輸入信號(hào)當(dāng)矩陣T的逆陣存在是，式(4.1.6)才有解。一般地，如果u(k)是隨機(jī)

19、序列或偽隨機(jī)二位式序列，則矩陣T是非奇異的，即( T)一1存在，式(4.1.6)有解。現(xiàn)在從T必須正定出發(fā)，討論對(duì)u(k)的要求。n N 1Tyyyuk nuyuu(4.1.7)當(dāng)N足夠大時(shí)有1 TRW.P.1yRuyRNRyuRu(4.1.8)如果矩陣T正定，則Ru是是對(duì)稱矩陣，并且各階主子式的行列式為正。當(dāng)N足夠大時(shí)，矩陣 Ru才是是對(duì)稱的。由此引出矩陣 T正定的必要條件是 u(k)為持續(xù)激勵(lì)信號(hào)。如果序列u(k)的n+1階方陣Ru是正定的，貝U稱 u(k)的n+1階持續(xù)激勵(lì)信號(hào)。下列隨機(jī)信號(hào)都能滿足 Ru正定的要求1. 有色隨機(jī)信號(hào)2. 偽隨機(jī)信號(hào)3. 白噪聲序列4.1.2最小二乘估計(jì)的

20、概率性質(zhì)最小二乘估計(jì)的概率性質(zhì)1) 無偏性由于輸出值y是隨機(jī)的，所以)是隨機(jī)的，但注意B不是隨機(jī)值。如果E ? E，則稱？是 無偏估計(jì)2) 致性估計(jì)誤差)的方差矩陣為2 彳1?1 TVarE(4.1.9)為有效值，參數(shù)(4.1.10)N N2limVar? lim R 10,w.p.1NN N上式表明當(dāng)nis時(shí)，)以概率1趨近于e o當(dāng)E (k)為不相關(guān)隨機(jī)序列時(shí)，最小二乘估計(jì)具有無偏性和一致性3) 有效性如果式(4.1.2)中的E的均值為0且服從正態(tài)分布的白噪聲向量，則最小二乘參數(shù)估計(jì)值 估計(jì)偏差的方差達(dá)到 Cramer-Rao不等式的下界，即Var? 2E t 1 M 1式中M為Fishe

21、r矩陣，且114)漸近正態(tài)性In p y/M EIn p y /如果式(4.1.2)中的E是均值為0且服從正態(tài)分布的白噪聲向量，則最小二乘參數(shù)估計(jì)值2E T 14.2遞推最小二乘法為了實(shí)現(xiàn)實(shí)時(shí)控制，必須采用遞推算法，這種辨識(shí)方法主要用于在線辨識(shí)設(shè)已獲得的觀測(cè)數(shù)據(jù)長度為N，則Yn(4.1.11) 服從正態(tài)分布，(4.1.2)(4.2.1)估計(jì)方差矩陣為式中?NVarPnPnNyn如果再獲得1組新的觀測(cè)值，則又增加 1個(gè)方程yN得新的參數(shù)估計(jì)值PN 1式中應(yīng)用矩陣求逆引理，可得PN 1和PN的遞推關(guān)系式2Pn|Yn1yN矩陣求逆引理設(shè)A為nx n矩陣，B和C為nx m矩陣,則有矩陣恒等式bct得到

22、遞推關(guān)系式PnPn由于N PNN 1是標(biāo)量，因而上式可以寫成PN 1PnPn最后，得最小二乘法辨識(shí)公式?NPnPN 1PnPn1并且A + BCTI+CTA-1B1Bcta1B1cta 1KnPnyNN PNN PNn Pn1 Rl1 Fn1 N 1PNN 1 PN都是非奇異矩陣，(4.2.2)(4.2.3)1113(424)35有2種給出初值的辦法(1) 設(shè)NO ( N0>n)為N的初始值，可算出PN0(2) 假定? O,P0 C2|,c是充分大的常數(shù),PN 0noYN 0(425)I為(2n+1) x (2n+1)單位矩陣，則經(jīng)過若干次遞推之后能得到較好的參數(shù)估計(jì)。 5兩種算法的實(shí)現(xiàn)

23、方案5.1最小二乘法一次完成算法實(shí)現(xiàn)如果把式(4.1.6)中的'取好足夠的輸入一輸出數(shù)據(jù)以后一次計(jì)算出來，那么這種算法就式最小二乘法的一次完成法。5.1.1最小二乘一次完成算法程序框圖圖3.2最小二乘一次完成算法程序框圖5.1.2 一次完成法程序具體程序參見附錄45.1.3 一次完成算法程序運(yùn)行結(jié)果ans =-1.50000.70001.00000.5000a1 =-1.50000.7000bl =1.00005.1.4辨識(shí)數(shù)據(jù)比較?at?t?1.50000.70001.0000.5000參數(shù)真值1.5000.7001.000.50誤差00005.1.5程序運(yùn)行曲線o-1o o oo

24、o OJK( 岀輸 )岀輸散離5.2遞推最小二乘法的實(shí)現(xiàn)5.2.1遞推算法實(shí)現(xiàn)步驟1）輸入M序列的產(chǎn)生程序，可以參見3.2當(dāng)中M序列產(chǎn)生方法（具體程序見附錄。2）初始化。一種簡(jiǎn)單的方法是選擇2（參照式Cl，其中C為盡量大的數(shù)，然后以它們?yōu)槠鹗贾颠M(jìn)行計(jì)算4.2.3）?？梢宰C明，當(dāng)Cis時(shí)，按照遞推公式算得的將與上面用批處理算法求得的結(jié)果相同，當(dāng)N很大時(shí)，C的選擇便無關(guān)緊要了。3）遞推算法的停機(jī)標(biāo)準(zhǔn):max?i k ?! k 1? k&是適當(dāng)小的數(shù)。4）最小二乘估計(jì)的遞推算法系統(tǒng)參數(shù)的步驟如下:u (i), y (i), i=1,2,n;1)產(chǎn)生系統(tǒng)輸入信號(hào) M 序列，輸入系統(tǒng)階次，計(jì)算輸

25、入和輸出數(shù)據(jù)2)置 N=n,?N0,PN 1010II 單位矩陣)3)形成行向量1 y(n)L y(Nn)u(N) Lu(N1n4)計(jì)算 K N 1PNN 1 PN5)輸入新測(cè)量數(shù)據(jù)y ( N+1 )和u( N+1 )；6)計(jì)算?N?1NKN 1 y(N 1) N?N7)計(jì)算PN1IKN 1 N 1 PN ;8)輸出?N1和 PN9)N N+1,轉(zhuǎn)(3);5.2.2 程序編制思路 ：(1)產(chǎn)生M序列，通過計(jì)算，依據(jù)算出輸出值，設(shè)定遞推初始值，采用424方法2設(shè)定辨識(shí)參數(shù)初值，?為一個(gè)很小的量，P0為一個(gè)很大值的4X 4矩陣，設(shè)定誤差量，作為停機(jī)標(biāo)準(zhǔn)。( 2)依據(jù)設(shè)定,計(jì)算 1y(2)y(1)

26、u(1) u(2) T , T1 P01 1可以算出為一標(biāo)量,依據(jù)式425算出K1,K1為4X 1矩陣。( 3)依據(jù)式 4.2.5 計(jì)算出 ?1,和 P1 ,加入估計(jì)參數(shù)矩陣；( 4)計(jì)算誤差大小,求出相對(duì)變化率,加入誤差矩陣第一列。(5)再以?和R為已知參數(shù)，重復(fù)(2)-( 3)計(jì)算出參數(shù) ?，并加入估計(jì)參數(shù)矩陣。6) 如此循環(huán),計(jì)算出參數(shù)估計(jì)?n 。7) 判斷誤差變化率滿足要求否,如果滿足,則退出循壞,不再進(jìn)行計(jì)算。8) 分離估計(jì)參數(shù),繪出圖形。523遞推最小二乘法程序框圖如下所示:圖5最小二乘遞推算法辨識(shí)的Malab程序流程圖具體程序參見附錄524程序運(yùn)行曲線圖6M序列信號(hào)的波形Para

27、meter Ide ntificati on with Recursive Least Squares Method-5000-l'l11t111-II-rrrr1IrIdentification Precision3002001000-100-200-300-4001020304050607080901005.2.5測(cè)試結(jié)果見附錄65.2.6地退數(shù)據(jù)表遞推次數(shù)a1召1誤差a1 ?a2召2誤;差a2?b誤差b, bb2b2誤;差d b1-1.500.0011.4990.70.0010.6991.00.0010.9990.50.0010.49992-1.501.50.700.71.001

28、.00.500.53-1.50.0011.4990.70.0010.6991.00.750.250.50.750.254-1.5-1.500.70.0010.6991.00.750.250.50.750.255-1.5-1.500.70.701.00.750.250.50.750.256-1.5-1.500.70.701.01.000.50.507-1.51.500.70.701.01.000.50.50從以上數(shù)據(jù)可以看出，隨著遞推次數(shù)的增加，參數(shù)估計(jì)的誤差逐漸減少，最終趨于穩(wěn)定。下面是辨識(shí)誤差的分布趨勢(shì)。（誤差=新估計(jì)舊估計(jì)值古計(jì)值）56結(jié)論最小二乘法是應(yīng)用廣泛的辨識(shí)方法之一?？捎糜趧?dòng)態(tài)系統(tǒng)

29、，也可用于靜態(tài)系統(tǒng)；可用于線性系統(tǒng)，也可用 于非線性系統(tǒng)。通過本課程設(shè)計(jì)，了解和掌握了基于最小二乘法原理的一次性完成算法和遞推算法的實(shí)現(xiàn)方 法，完成了 Matlab 下的仿真計(jì)算，能夠精確的估計(jì)辨識(shí)參數(shù)。最小二乘一次性完成算法是離線算法，需要采 集大量數(shù)據(jù)，一次性完成計(jì)算，因此，數(shù)據(jù)計(jì)算量大，當(dāng)數(shù)據(jù)量很大時(shí)，數(shù)據(jù)輸入不方便，但在本課程設(shè)計(jì) 過程當(dāng)中，考慮到了此問題，運(yùn)用相應(yīng)的辦法，解決了矩陣輸入的問題。遞推算法適合于在線算法，利用原 有參數(shù)估計(jì)進(jìn)行下一步估計(jì)，可以做到運(yùn)算量小，實(shí)時(shí)進(jìn)行估計(jì)，根據(jù)仿真結(jié)果圖示，可以看出計(jì)算一定量 的數(shù)據(jù)后，參數(shù)估計(jì)趨于穩(wěn)定，只要滿足誤差要求，不必計(jì)算完所有數(shù)據(jù)。

30、兩種算法不足之處，在考慮噪聲 干擾時(shí)，看成了白噪聲，具有不相關(guān)性。如果噪聲引起的方程誤差是相關(guān)的，而不是白噪聲，可以通過下面 兩種方法進(jìn)行估計(jì)先估計(jì)未受噪音污染的系統(tǒng)輸出，然后再估計(jì)模型的參數(shù)，就是輔助變量法；使用一 種濾波器，將相關(guān)的方程的誤差轉(zhuǎn)換為白噪聲再進(jìn)行估計(jì)，就是增廣矩陣法。7 參考文獻(xiàn)1 李言俊 張科,系統(tǒng)辨識(shí)理論及應(yīng)用 ,國防工業(yè)出版社 ,2006.7 .2 劉宏才 ,系統(tǒng)辨識(shí)與參數(shù)估計(jì) ,北京,冶金工業(yè)出版社 ,1999.3 黃道平，MATLAB與控制系統(tǒng)的數(shù)字仿真及 CAD化學(xué)工業(yè)出版社,2004.104 侯媛彬汪梅王立奇系統(tǒng)辨識(shí)及其MATLAB仿真科學(xué)出版社 2004.25

31、 蔡季冰，系統(tǒng)辨識(shí)，北京，北京理工大學(xué)出版社，1991.8 附錄附錄 1：用乘同余法產(chǎn)生白噪聲 編程如下：A=6; x0=1; M=255; f=2;N=100 ；%初始化；x0=1; M=255;for k=1: N x2=A*x0;x1=mod （x2，M）; v1=x1/256;%乘同余法遞推 100次；%分別用x2和x0表示xi+i和xi-1 ；%取x2存儲(chǔ)器的數(shù)除以M的余數(shù)放x1(xi)中；%將x1存儲(chǔ)器中的數(shù)除以256得到小于1的隨機(jī)數(shù)放v1中；v（:，k）=（v1-0.5 ）*f;%將v1中的數(shù)(Q減去0.5再乘以存儲(chǔ)器f中的系數(shù)，存放在矩陣存儲(chǔ)器v的第k列中，v(:,k)表示行

32、不變、列隨遞推循環(huán)次數(shù)變化；x0=x1;v0=v1;% xi-1= xi；end v2=v%遞推100次結(jié)束；%該語句后無；'實(shí)現(xiàn)矩陣存儲(chǔ)器 v中隨機(jī)數(shù)放在v2中，且可直接顯示在MATLAB 的 window 中；k1=k; %grapher k=1:k1;%以下是繪圖程序；plot(k，v，k，v， 'r ');xlabel(k')， ylabel( v');tktle( (-1，+1)均勻分布的白噪聲 ') 1.2 程序運(yùn)行結(jié)果如圖 6所示。圖6采用MATLAB產(chǎn)生的（-1,+1）均勻分布的白噪聲序列附錄2產(chǎn)生的（-1 , 1）均勻分布的白噪

33、聲序列在程序運(yùn)行結(jié)束后，產(chǎn)生的（-1, 1）均勻分布的白噪聲序列，直接從 MATLAB的window界面中copy出來如下（v2中每行存6個(gè)隨機(jī)數(shù)）：v2 =-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.718

34、80.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-

35、0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359附錄3 M序列產(chǎn)生程序function Out=Mge nerator(le ngth,ple ngth,a)%plength=4;length=15; %M序列周期=2yiength-1 ,% 置M序列總長度和級(jí)數(shù)，a是幅度for n=1:plength%移位寄存器輸入 X(i)初T態(tài)(1111 )X( n)=1;endfor i=1:

36、le ngthfor j=ple ngth:-1:2X(j)=X(j-1);endX(1)=xor(X(ple ngth-1),X(ple ngth);%異或運(yùn)算if X(ple ngth)=OOut(i)=-1;elseOut(i)=X(ple ngth);endend%save ('mdata','Out');%繪圖i1=ik=1:1:le ngth;plot(k,Out,k,Out,'rx')xlabel('k')ylabel('M 序列')title('移位寄存器產(chǎn)生的M序列')end附錄4

37、程序運(yùn)行結(jié)果如下圖所示移位寄存器產(chǎn)生的ri序列圖7移位寄存器產(chǎn)生的M序列附錄4最小二乘一次完成法%Oncefinish %model is y(k)-1.5y(k-1)+0.7(k-2)=u(k-1)+0.5(k-2)+v(k); L=50;A=1;n=2; %output number; u=mseries(A,L);%u=Mserise(L,M,A); %系統(tǒng)辨識(shí)的輸入信號(hào)為一個(gè)周期的M序列，信號(hào)長度L ;幅度A;z=zeros(1,16);%定義輸出觀測(cè)值的長度 1*16dd=zeros(1,4); gg=zeros(14,1);for k=n+1:L+1z(k)=1.5*z(k-1)-

38、0.7*z(k-2)+u(k-1)+0.5*u(k-2); % 用理想輸出值作為觀測(cè)值end subplot(3,1,1) %畫三行一列圖形窗口中的第一個(gè)圖形 stem(u) %畫出輸入信號(hào) u 的經(jīng)線圖形 ylabel(' 輸入 u(k)')subplot(3,1,2) i=1:1:L+1;plot(i,z)ylabel(' 輸出 z(k)')subplot(3,1,3)%畫三行一列圖形窗口中的第三個(gè)圖形stem(z),grid on%畫出輸出觀測(cè)值 z 的經(jīng)線圖形，并顯示坐標(biāo)網(wǎng)格ylabel(' 離散輸出 z(k)')%u,z;%顯示輸入信號(hào)

39、和輸出觀測(cè)信號(hào)%L=L-1% 數(shù)據(jù)長度for i=n+1:L+1h=-z(i-1) -z(i-2) u(i-1) u(i-2);dd(i-2,:)=h;%dd=(L-1)*2nend%HL=-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11

40、) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14) % 給樣本矩陣 HL 賦值 for j=n+1:L+1zz=z(j);gg(j-2,:)=zz;end%c=inv(dd'*dd)*dd'*gg;給樣本矩陣 zL 賦值%ZL=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)%ca

41、lculating parameters% 計(jì)算參數(shù)c1=dd'*dd; c2=inv(c1); c3=dd'*gg; c=c2*c3;c'%c1=dd'*dd; c2=inv(c1); c3=dd'*ZL; c=c2*c3 % 計(jì)算并顯示 %DISPLAY PARAMETERS a1=c(1), a2=c(2), b1=c(3), b2=c(4); % 從 中分離出并顯示 a1 、 a2、 b1、 b2 %End 附錄 5 遞推最小二乘算法程序 clear;clc;% 清理工作間變量L=100;% M 序列的周期 y1=1;y2=1;y3=1;y4=0

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

43、 -0.03else u(i)=0.5;% 當(dāng)M序列的值為"0"時(shí),辨識(shí)的輸入信號(hào)取“ 0.03 ” end% 小循環(huán)結(jié)束y1=x1;y2=x2;y3=x3;y4=x4;% 為下一次的輸入信號(hào)做準(zhǔn)備en d%大循環(huán)結(jié)束，產(chǎn)生輸入信號(hào) u figure(1);% 第 1 個(gè)圖形stem(u),grid on% 以徑的形式顯示出輸入信號(hào)并給圖形加上網(wǎng)格 z(2)=0;z(1)=0;% 取 z 的前兩個(gè)初始值為零 for k=3:L;% 循環(huán)變量從 3 到 25z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%給出理想的辨識(shí)輸出采樣信號(hào)e

44、nd %RLS遞推最小二乘辨識(shí)c0=0.001 0.001 0.001 0.001'% 直接給出被辨識(shí)參數(shù)的初始值 ,即一個(gè)充分小的實(shí)向量p0=10A6*eye(4,4);% 直接給出初始狀態(tài)P0,即一個(gè)充分大的實(shí)數(shù)單位矩陣E=0.000000005;% 相對(duì)誤差 E=0.000000005c=c0,zeros(4,99);% 被辨識(shí)參數(shù)矩陣的初始值及大小e=zeros(4,100);% 相對(duì)誤差的初始值及大小for k=3:L; % 開始求 Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2)' x=h1'*p0*h1+1; x1=inv(x); %開始

45、求 K(k)k仁pO*h1*x1;%求出K的值d1=z(k)-h1'*c0; c1=c0+k1*d1;%求被辨識(shí)參數(shù) ce1=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ù)矩陣的最后一列p1=p0-k1*k1'*h1'*p0*h1+1;% 求出 p(k) 的值p0=p1;% 給下次用if e2<=E break;% 若參數(shù)收斂滿足要求，終止計(jì)算end

46、% 小循環(huán)結(jié)束end%大循環(huán)結(jié)束c%顯示被辨識(shí)參數(shù)e%顯示辨識(shí)結(jié)果的收斂情況%分離參數(shù)a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); figure(2);% 第 2 個(gè)圖形i=1:L;% 橫坐標(biāo)從 1 到 15 plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %畫出 a1， a2， b1， b2 的各次辨識(shí)結(jié)果title('Parameter Identification with Recursive Least Squares Method')%圖形標(biāo)題figure(3); % 第 3 個(gè)圖形 i=1:L; % 橫坐標(biāo)從 1 到 15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %畫出 a1，a2，b1，b2 的各次辨識(shí)結(jié)果的收斂情況title('Identification Precision') %圖形標(biāo)題附錄 6 遞推最小二乘法程序測(cè)試結(jié)果c =Columns 1 t