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

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

2、5.1.4辨識數(shù)據(jù)比較125.1.5程序運行曲線125.2遞推最小二乘法的實現(xiàn)125.2.1遞推算法實現(xiàn)步驟125.2.2程序編制思路：135.2.3遞推最小二乘法程序框圖145.2.4程序運行曲線155.2.5測試結(jié)果165.2.6遞推數(shù)據(jù)表166結(jié)論167參考文獻178附錄17應(yīng)用最小二乘一次完成法和遞推最小二乘法算法的系統(tǒng)辨識摘要：本題針對一個單輸入單輸出系統(tǒng)的便是問題，辨識的輸入信號采用的是偽隨機二位式序列（M序列），系統(tǒng)噪聲為獨立同分布高斯隨機向量序列（白噪聲），辨識的算法是遞推最小二乘法和廣義最小二乘法，本文簡單描述應(yīng)用經(jīng)典辨識方法的辨識方案，詳細描述了輸入信號、噪聲的產(chǎn)生方法及m

3、atlab程序，闡述了用兩種不同算法的辨識原理并對它們的推導(dǎo)過程及辨識程序編制思路做了詳細的描述。最后結(jié)合真值與估計值對不同辨識算法的優(yōu)劣進行了比較。關(guān)鍵詞：系統(tǒng)辨識 M序列 最小二乘法 1引言1.1概述系統(tǒng)辨識是現(xiàn)代控制理論中的一個分支，它是根據(jù)系統(tǒng)的輸入輸出時間函數(shù)來確定描述系統(tǒng)行為的數(shù)學(xué)模型。通過辨識建立數(shù)學(xué)模型的目的是估計表征系統(tǒng)行為的重要參數(shù)，建立一個能模仿真實系統(tǒng)行為的模型，用當(dāng)前可測量的系統(tǒng)的輸入和輸出預(yù)測系統(tǒng)輸出的未來演變，以及設(shè)計控制器。對系統(tǒng)進行分析的主要問題是根據(jù)輸入時間函數(shù)和系統(tǒng)的特性來確定輸出信號。對系統(tǒng)進行控制的主要問題是根據(jù)系統(tǒng)的特性設(shè)計控制輸入，使輸出滿足預(yù)先規(guī)

4、定的要求。而系統(tǒng)辨識所研究的問題恰好是這些問題的逆問題。通常，預(yù)先給定一個模型類M（即給定一類已知結(jié)構(gòu)的模型），一類輸入信號u和等價準則JL(y，yM)(一般情況下，J是誤差函數(shù)，是過程輸出y和模型輸出yM的一個泛函)；然后選擇使誤差函數(shù)J達到最小的模型，作為辨識所要求的結(jié)果。系統(tǒng)辨識包括兩個方面：結(jié)構(gòu)辨識和參數(shù)估計。在實際的辨識過程中，隨著使用的方法不同，結(jié)構(gòu)辨識和參數(shù)估計這兩個方面并不是截然分開的，而是可以交織在一起進行的。1.2辨識的基本步驟先驗知識和建模目的的依據(jù)。先驗知識指關(guān)于系統(tǒng)運動規(guī)律、數(shù)據(jù)以及其他方面的已有知識。這些知識對選擇模型結(jié)構(gòu)、設(shè)計實驗和決定辨識方法等都有重要作用。用于

5、不同目的的模型可能會有很大差別。實驗設(shè)計。辨識是從實驗數(shù)據(jù)中提取有關(guān)系統(tǒng)信息的過程，設(shè)計實驗的目標之一是要使所得到的數(shù)據(jù)能包含系統(tǒng)更多的信息。主要包括輸入信號設(shè)計，采樣區(qū)間設(shè)計，預(yù)采樣濾波器設(shè)計等。結(jié)構(gòu)辨識。即選擇模型類中的數(shù)學(xué)模型M的具體表達形式。除線性系統(tǒng)的結(jié)構(gòu)可通過輸入輸出數(shù)據(jù)進行辨識外  ，一般的模型結(jié)構(gòu)主要通過先驗知識獲得。參數(shù)估計。知道模型的結(jié)構(gòu)后，用輸入輸出數(shù)據(jù)確定模型中的未知參數(shù)。實際測量都是有誤差的，所以參數(shù)估計以統(tǒng)計方法為主。模型適用性檢驗。造成模型不適用主要有三方面原因：模型結(jié)構(gòu)選擇不當(dāng)；實驗數(shù)據(jù)誤差過大或數(shù)據(jù)代表性太差；辨識算法存在問題。檢驗方法主要

6、有利用先驗知識檢驗和利用數(shù)據(jù)檢驗兩類。11.3設(shè)待辨識系統(tǒng)如圖1所示。 二階系統(tǒng)為：圖1待辨識系統(tǒng)圖參數(shù)真值為：1設(shè)，為白色噪聲，簡單描述應(yīng)用經(jīng)典辨識方法的辨識方案。2為有色噪聲，為獨立同分布的高斯序列，2系統(tǒng)辨識輸入信號的產(chǎn)生方法和理論依據(jù)2.1白噪聲序列如果隨機序列 均值為0，并且兩兩不相的關(guān)的，對應(yīng)自相關(guān)函數(shù)為式中則稱這種隨機序列為白噪聲序列2.1.1白噪聲序列的產(chǎn)生方法下面主要介紹（0，1）均勻分布和正態(tài)分布隨機數(shù)的產(chǎn)生方法 在計算機上產(chǎn)生（0，1）均勻分布隨機數(shù)的方法很多，其中最簡單、最方便的是數(shù)學(xué)方法。產(chǎn)生偽隨機數(shù)的數(shù)學(xué)方法很多，其中最常用的是乘同余法和混合同余法。（1）乘同余法這

7、種方法先用遞推同余式產(chǎn)生正整數(shù)序列Xi=Axi-1(modM),i=1,2,3式中：M為2的方冪，k為大于2的整數(shù)；A3(mod8)或A5(mod8),且A不能太??；初值x0取正奇數(shù)，例如取x0=1. 再令 則 是偽隨機序列，循環(huán)周期可達 。（2）混合同余法混合同余法產(chǎn)生偽隨機數(shù)的遞推同余式為式中 K為大于2的整數(shù)；A1(mod4),即，其中n為滿足關(guān)系式2n34的整數(shù)。初值x0為非負整數(shù)。令則 是循環(huán)周期為 的偽隨機數(shù)序列2.2 M序列的產(chǎn)生在進行系統(tǒng)辨識時，選用白噪聲作為辨識輸入信號可以保證獲得較好的便是效果，但工程上難以實現(xiàn)。M序列是一種很好的辨識輸入信號，它具有近似白燥聲的性質(zhì)，不僅可

8、以保證有較好的辨識效果，而且工程上易于實現(xiàn)。M序列是偽隨機二位式序列的一種形式。在介紹M序列之前，先介紹一下偽隨機噪聲的概念。2.2.1 偽隨機噪聲對白噪聲的一個樣本函數(shù)w(t)截取0,T時間內(nèi)一段，對其它時間段T,2T,2T,3T,以 周期 T延拖下去，這樣獲得的函數(shù)w(t)是周期T的函數(shù)，在0,T時間內(nèi)是白噪聲，在此時間之外是重復(fù)的白噪聲，它的自相關(guān)函數(shù) 的周期也是T。由于在0,T時間內(nèi)自相關(guān)函數(shù) 就是白噪聲的自相關(guān)函數(shù)，它具有周期性，稱為w(t)為偽隨機噪聲。2.2.2 M序列的產(chǎn)生方法M序列是一種離散二位式：隨機序列，所謂“二位式”是指每個隨機變量只有2種狀態(tài)?？捎枚嗉壘€性反饋移位寄存

9、器產(chǎn)生M序列。M序列是最長線性反饋移存器序列的簡稱，是由帶線性反饋的移存器產(chǎn)生的周期最長的一種序列。具有較強的抗干擾能力和較低的截獲概率，而且長的M序列更容易在一定的強噪聲中被提取，這樣就能夠充分保證數(shù)據(jù)的正常通信。通常產(chǎn)生偽隨機序列的電路為反饋移存器.一般說來, 由n級移位寄存器產(chǎn)生的周期為N=2-1的M序列，在一個循環(huán)周期內(nèi)，“0”出現(xiàn)的次數(shù)為,“1”出現(xiàn)的次數(shù)為?，F(xiàn)在我們引入M序列的本原多項式的概念。若一個n次多項式滿足以下條件（1）為既約的。（2）可整除,。（3）除不盡,則為本原多項式。 一個4級M序列可以通過線性反饋移位寄存器產(chǎn)生，如下圖所示：X1X2X3X4移位脈沖圖2 周期為15

10、的偽隨機序列產(chǎn)生器圖每級移位寄存器由雙穩(wěn)態(tài)觸發(fā)器和門電路組成，稱為1位，分別以0和1來表示2種狀態(tài)。當(dāng)移位脈沖到來時，每位的內(nèi)容移至下一位，最后1位移出的內(nèi)容即為輸出。為了保持連續(xù)工作，將最后2級寄存器的內(nèi)容經(jīng)過適當(dāng)?shù)倪壿嬤\算后反饋到第1級寄存器作為輸入。當(dāng)一個移位脈沖到來后，第1級寄存器的內(nèi)容送到第2級，第2級寄存器的內(nèi)容送到第3級，第3級寄存器的內(nèi)容送到第4級，而第3級和第4級寄存器的內(nèi)容作模和2相加后再反饋到第1級寄存器。產(chǎn)生偽隨機序列時要求寄存器的初始狀態(tài)不全為0，因為全0初始狀態(tài)將導(dǎo)致各級寄存器輸出永遠是0。如果寄存器的初始內(nèi)容都是1，第1個移位脈沖來到后，4級寄存器的內(nèi)容變以011

11、1，一個周期的變化規(guī)律為:1111 0111 0011 0001 1000 0100 0010 1001 1100 0110 1011 0101 1010 1101 1110 1111一個周期結(jié)束后，產(chǎn)生的15種不同的狀態(tài)。計算機模擬產(chǎn)生M序列非常方便，先定義輸出序列長度和一個數(shù)組，數(shù)組個數(shù)等于移位寄存器的個數(shù)，通過使用異或指令，再利用for循環(huán)指令，即可完成任意長度和級數(shù)M序列的產(chǎn)生。i<L停機輸入長度L,級數(shù)M置X初始狀態(tài)輸出，U(i)清零X(1)X(M)移位一次X(M) xor X(M-1)U(i)X(M) 圖3M序列產(chǎn)生的流程圖圖4MATLAB仿真M序列圖形如圖3是產(chǎn)生M級長度為

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

13、2經(jīng)典辨識方法的實現(xiàn)由系統(tǒng)辨識圖可得如下關(guān)系統(tǒng)： （3.2.1）其中：為輸入白噪聲信號，為實際輸入信號，為實際輸出信號，為輸入信號為時的脈沖響應(yīng)。根據(jù)維納-霍夫積分方程可得離散的維納霍夫積分方程： （3.2.2）令式中、則上式可寫成以下形式： （3.2.3）再設(shè)： 再根據(jù)上式（3.2.3）可得： （3.2.4）由二電平M序列的自相關(guān)函數(shù)的計算公式可以求得 (3.2.5)根據(jù)遞推公式求m次觀測的脈沖響應(yīng)對系統(tǒng)進行m次觀測，由m次觀測得到的用來表示，則有如下遞推式： （3.2.6）綜合以上公式就可以得到以下遞推公式： (3.2.7)根據(jù)遞推式，可從及新的觀測數(shù)據(jù)得到，隨著觀測數(shù)據(jù)的增加，的精確度不

14、斷的提高。134最小二乘法的理論基礎(chǔ)4.1最小二乘法設(shè)單輸入單輸出線性定長系統(tǒng)的差分方程表示為：其中 n(k)均值為0的白噪聲,現(xiàn)分別測出n+N個輸出輸入值y(1),y(2),y(n+N),u(1),u(2),u(n+N),則可寫出N個方程，寫成向量矩陣形式(4.1.1)則式(4.1.1)可寫為 （4.1.2）式中：y為N維輸出向量；為N為維噪聲向量；為（2n+1）維參數(shù)向量；為N×（2n+1）測量矩陣。因此，式（4.1.1）是一個含有（2n+1）個未知參數(shù)，由N個方程組成的聯(lián)立方程組。在給定輸出向量y和測量矩陣的條件下求參數(shù)的估計，這就是系統(tǒng)辨識問題。設(shè) 表示 的估計值，表示y的最

15、優(yōu)估計，則有 （4.1.3）式中：設(shè)e(k)=y(k)- (k), e(k)稱為殘差，則有e=y- =y- 最小二乘估計要求殘差的平方和最小，即按照指數(shù)函數(shù)（4.1.4）求對 的偏導(dǎo)數(shù)并令其等于0可得：（4.1.5）由式（4.1.5）可得的 最小二乘估計：（4.1.6）J為極小值的充分條件是：即矩陣T為正定矩陣，或者說是非奇異的。4.1.1最小二乘法估計中的輸入信號當(dāng)矩陣T的逆陣存在是，式（4.1.6）才有解。一般地，如果u(k)是隨機序列或偽隨機二位式序列，則矩陣T是非奇異的，即（T）1存在，式（4.1.6）有解?，F(xiàn)在從T必須正定出發(fā)，討論對u(k)的要求。（4.1.7）當(dāng)N足夠大時有（4.

16、1.8）如果矩陣T正定，則Ru是是對稱矩陣，并且各階主子式的行列式為正。當(dāng)N足夠大時，矩陣Ru才是是對稱的。由此引出矩陣T正定的必要條件是u(k)為持續(xù)激勵信號。如果序列u(k)的n+1階方陣Ru是正定的，則稱u(k)的n+1階持續(xù)激勵信號。下列隨機信號都能滿足Ru正定的要求1. 有色隨機信號2. 偽隨機信號3. 白噪聲序列4.1.2最小二乘估計的概率性質(zhì)最小二乘估計的概率性質(zhì)1） 無偏性由于輸出值y是隨機的，所以是隨機的，但注意不是隨機值。如果，則稱是無偏估計2）一致性估計誤差的方差矩陣為（4.1.9）上式表明當(dāng)N時，以概率1趨近于。當(dāng)(k)為不相關(guān)隨機序列時，最小二乘估計具有無偏性和一致性

17、3)有效性如果式（4.1.2）中的的均值為0且服從正態(tài)分布的白噪聲向量，則最小二乘參數(shù)估計值 為有效值，參數(shù)估計偏差的方差達到Cramer-Rao不等式的下界，即（4.1.10）式中M為Fisher矩陣，且（4.1.11）4）漸近正態(tài)性如果式（4.1.2）中的是均值為0且服從正態(tài)分布的白噪聲向量，則最小二乘參數(shù)估計值服從正態(tài)分布，即（4.1.2）4.2遞推最小二乘法為了實現(xiàn)實時控制，必須采用遞推算法，這種辨識方法主要用于在線辨識（4.2.1）設(shè)已獲得的觀測數(shù)據(jù)長度為N，則估計方差矩陣為式中于是如果再獲得1組新的觀測值，則又增加1個方程得新的參數(shù)估計值式中應(yīng)用矩陣求逆引理，可得和的遞推關(guān)系式矩陣

18、求逆引理 設(shè)A為n×n矩陣，B和C為n×m矩陣，并且A，ABCT和I+CTA-1B都是非奇異矩陣，則有矩陣恒等式（4.2.2）得到遞推關(guān)系式由于是標量，因而上式可以寫成（4.2.3）最后，得最小二乘法辨識公式（4.2.4）有2種給出初值的辦法（1）設(shè)N0（N0>n）為N的初始值，可算出（4.2.5）（2）假定C是充分大的常數(shù)，I為(2n+1) ×(2n+1)單位矩陣，則經(jīng)過若干次遞推之后能得到較好的參數(shù)估計。1245兩種算法的實現(xiàn)方案5.1最小二乘法一次完成算法實現(xiàn)如果把式（4.1.6）中的取好足夠的輸入輸出數(shù)據(jù)以后一次計算出來，那么這種算法就式最小二乘法的

19、一次完成法。5.1.1最小二乘一次完成算法程序框圖賦輸入信號初值u定義輸出觀測值的長度并計算系統(tǒng)的輸出值畫出輸入和輸出觀測值的圖形給樣本矩陣HL和zL賦值計算參數(shù)從中分離出并顯示出被辨識參數(shù)a1, a2, b1, b2停機圖3.2 最小二乘一次完成算法程序框圖5.1.2一次完成法程序具體程序參見附錄45.1.3一次完成算法程序運行結(jié)果ans = -1.5000 0.7000 1.0000 0.5000a1 = -1.5000a2 = 0.7000b1 =1.00005.1.4辨識數(shù)據(jù)比較1.50000.70001.0000.5000參數(shù)真值1.5000.7001.000.50誤差00005.1

20、.5程序運行曲線5.2遞推最小二乘法的實現(xiàn)5.2.1遞推算法實現(xiàn)步驟1）輸入M序列的產(chǎn)生程序，可以參見3.2當(dāng)中M序列產(chǎn)生方法（具體程序見附錄。）2）初始化。 一種簡單的方法是選擇和，其中C為盡量大的數(shù)，然后以它們?yōu)槠鹗贾颠M行計算(參照式4.2.3)?？梢宰C明，當(dāng)C時，按照遞推公式算得的將與上面用批處理算法求得的結(jié)果相同，當(dāng)N很大時，C的選擇便無關(guān)緊要了。3）遞推算法的停機標準：，是適當(dāng)小的數(shù)。4）最小二乘估計的遞推算法系統(tǒng)參數(shù)的步驟如下：（1）產(chǎn)生系統(tǒng)輸入信號M序列，輸入系統(tǒng)階次，計算輸入和輸出數(shù)據(jù)u（i），y（i），i=1,2,n;（2）置N=n，（I單位矩陣）；（3）形成行向量（4）計算

21、（5）輸入新測量數(shù)據(jù)y ( N+1 )和u（N+1）；（6）計算（7）計算（8）輸出和；（9）NN+1,轉(zhuǎn)（3）；5.2.2程序編制思路：（1）產(chǎn)生M序列，通過計算，依據(jù)算出輸出值，設(shè)定遞推初始值，采用4.2.4方法2設(shè)定辨識參數(shù)初值，為一個很小的量，P0為一個很大值的4×4矩陣，設(shè)定誤差量，作為停機標準。（2）依據(jù)設(shè)定，計算，可以算出為一標量，依據(jù)式4.2.5算出K1,K1為4×1矩陣。（3）依據(jù)式4.2.5計算出，和，加入估計參數(shù)矩陣；（4）計算誤差大小，求出相對變化率，加入誤差矩陣第一列。（5）再以和為已知參數(shù)，重復(fù)（2）（3）計算出參數(shù)，并加入估計參數(shù)矩陣。（6）如

22、此循環(huán)，計算出參數(shù)估計。（7）判斷誤差變化率滿足要求否，如果滿足，則退出循壞，不再進行計算。（8）分離估計參數(shù)，繪出圖形。5.2.3遞推最小二乘法程序框圖如下所示：Y工作間清零產(chǎn)生輸出采樣信號給被辨識參數(shù)和P賦初值按照式(4.2.4)的第三式計算P(k)計算被辨識參數(shù)的相對變化量參數(shù)收斂滿足要求？停機按照式(4.2.4)的第二式計算K(k)按照式(4.2.4)的第一式計算(k)第四個移位寄存器的輸出取反，并將幅值變?yōu)?.03得到辨識系統(tǒng)的輸入信號樣本值給M序列的長度L和移位寄存器的輸入賦初始值畫出被辨識參數(shù)的各次遞推估計值圖形分離參數(shù)畫出被辨識參數(shù)的相對誤差的圖形畫出辨識的輸入信號徑線圖形圖5

23、 最小二乘遞推算法辨識的Malab程序流程圖具體程序參見附錄5.2.4程序運行曲線圖6M序列信號的波形5.2.5測試結(jié)果見附錄65.2.6地退數(shù)據(jù)表遞推次數(shù)誤差誤差誤差誤差1-1.500.0011.4990.70.0010.6991.00.0010.9990.50.0010.49992-1.501.50.700.71.001.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.

24、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ù)估計的誤差逐漸減少，最終趨于穩(wěn)定。下面是辨識誤差的分布趨勢。（誤差=）56結(jié)論最小二乘法是應(yīng)用廣泛的辨識方法之一?？捎糜趧討B(tài)系統(tǒng)，也可用于靜態(tài)系統(tǒng)；可用于線性系統(tǒng)，也可用于非線性系統(tǒng)。通過本課程設(shè)計，了解和掌握了基于最小二乘法原理的一次性完成算法和遞推算法的實現(xiàn)方法，完成了Matlab下的仿真計算，能夠精確的估計辨識參數(shù)。最小二乘一次性完成算法是離線算法，需要采集大量數(shù)據(jù)，一次性完成計算，因此，數(shù)據(jù)計算量大，當(dāng)數(shù)據(jù)量很

25、大時，數(shù)據(jù)輸入不方便，但在本課程設(shè)計過程當(dāng)中，考慮到了此問題，運用相應(yīng)的辦法，解決了矩陣輸入的問題。遞推算法適合于在線算法，利用原有參數(shù)估計進行下一步估計，可以做到運算量小，實時進行估計，根據(jù)仿真結(jié)果圖示，可以看出計算一定量的數(shù)據(jù)后，參數(shù)估計趨于穩(wěn)定，只要滿足誤差要求，不必計算完所有數(shù)據(jù)。兩種算法不足之處，在考慮噪聲干擾時，看成了白噪聲，具有不相關(guān)性。如果噪聲引起的方程誤差是相關(guān)的，而不是白噪聲，可以通過下面兩種方法進行估計先估計未受噪音污染的系統(tǒng)輸出，然后再估計模型的參數(shù)，就是輔助變量法；使用一種濾波器，將相關(guān)的方程的誤差轉(zhuǎn)換為白噪聲再進行估計，就是增廣矩陣法。7參考文獻1李言俊 張科,系統(tǒng)

26、辨識理論及應(yīng)用,國防工業(yè)出版社,2006.7 .2劉宏才,系統(tǒng)辨識與參數(shù)估計,北京,冶金工業(yè)出版社,1999.3黃道平 ,MATLAB與控制系統(tǒng)的數(shù)字仿真及CAD化學(xué)工業(yè)出版社, 2004.104侯媛彬 汪梅 王立奇 系統(tǒng)辨識及其MATLAB仿真科學(xué)出版社 2004.25蔡季冰，系統(tǒng)辨識，北京，北京理工大學(xué)出版社，1991.8附錄附錄1：用乘同余法產(chǎn)生白噪聲 編程如下：A=6; x0=1; M=255; f=2; N=100； %初始化；x0=1; M=255;for k=1: N %乘同余法遞推100次； x2=A*x0; %分別用x2和x0表示xi+1和xi-1； x1=mod (x2,M

27、); %取x2存儲器的數(shù)除以M的余數(shù)放x1（xi）中； v1=x1/256; %將x1存儲器中的數(shù)除以256得到小于1的隨機數(shù)放v1中； v(:,k)=(v1-0.5 )*f; %將v1中的數(shù)（）減去0.5再乘以存儲器f中的系數(shù)，存放在矩陣存儲器v的第k列中，v(:,k)表示行不變、列隨遞推循環(huán)次數(shù)變化； x0=x1; % xi-1= xi； v0=v1;end %遞推100次結(jié)束；v2=v %該語句后無；，實現(xiàn)矩陣存儲器v中隨機數(shù)放在v2中，且可直接顯示在MATLAB的window中；k1=k;%grapher %以下是繪圖程序；k=1:k1;plot(k,v,k,v,r);xlabel(k

28、), ylabel(v);tktle( (-1,+1)均勻分布的白噪聲)1.2程序運行結(jié)果如圖6所示。圖6 采用MATLAB產(chǎn)生的(-1,+1)均勻分布的白噪聲序列附錄2 產(chǎn)生的(-1，1)均勻分布的白噪聲序列在程序運行結(jié)束后，產(chǎn)生的(-1，1)均勻分布的白噪聲序列，直接從MATLAB的window界面中copy出來如下（v2中每行存6個隨機數(shù)）：v2 = -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.92190.5703 0.4531 -0.2500 -0.4844 0.1016 -0.36720.8047 -0.1328 0.2188 0.3359 -0

29、.9531 -0.71880.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.48440.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2

30、500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.71880.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.13280.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.48440.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531

31、 -0.7188 0.6875 -0.8359附錄3 M序列產(chǎn)生程序function Out=Mgenerator(length,plength,a)%plength=4;length=15; %M序列周期=2plength-1 ,%置M序列總長度和級數(shù) ，a是幅度 for n=1:plength %移位寄存器輸入X(i)初T態(tài)（1111） X(n)=1;endfor i=1:length for j=plength:-1:2 X(j)=X(j-1); end X(1)=xor(X(plength-1),X(plength); %異或運算 if X(plength)=0 Out(i)=-1;

32、else Out(i)=X(plength); end end%save ('mdata','Out');%繪圖i1=ik=1:1:length;plot(k,Out,k,Out,'rx')xlabel('k')ylabel('M序列')title('移位寄存器產(chǎn)生的M序列')end附錄4程序運行結(jié)果如下圖所示圖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

33、;A=1;n=2; %output number; u=mseries(A,L);%u=Mserise(L,M,A); %系統(tǒng)辨識的輸入信號為一個周期的M序列,信號長度L；幅度A；z=zeros(1,16); %定義輸出觀測值的長度1*16dd=zeros(1,4);gg=zeros(14,1);for k=n+1:L+1 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想輸出值作為觀測值endsubplot(3,1,1) %畫三行一列圖形窗口中的第一個圖形stem(u) %畫出輸入信號u的經(jīng)線圖形ylabel('輸入u(k)'

34、)subplot(3,1,2) i=1:1:L+1; plot(i,z) ylabel('輸出z(k)')subplot(3,1,3) %畫三行一列圖形窗口中的第三個圖形stem(z),grid on %畫出輸出觀測值z的經(jīng)線圖形，并顯示坐標網(wǎng)格ylabel('離散輸出z(k)')%u,z; %顯示輸入信號和輸出觀測信號%L=L-1%數(shù)據(jù)長度for i=n+1:L+1 h=-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)

35、 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) -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) %給樣

36、本矩陣HL賦值for j=n+1:L+1 zz=z(j); gg(j-2,:)=zz;end %c=inv(dd'*dd)*dd'*gg; %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)% 給樣本矩陣zL賦值%calculating parameters%計算參數(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

37、*c3 %計算并顯示 %DISPLAY PARAMETERSa1=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;%四個移位積存器的輸出初始值for i=1:L;%開始循環(huán)，長度為L x1=xor(y3,y4);%第一個移位積存器的輸入是第3個與第4個移位積存器的輸出的“或” x2=y1;%第二個移位積存器的輸入是第3個移位積存器的輸出 x3=y2;%第三個移位積存器的輸入是第2個移位積存器

38、的輸出 x4=y3;%第四個移位積存器的輸入是第3個移位積存器的輸出 y(i)=y4;%取出第四個移位積存器幅值為"0"和"1"的輸出信號， if y(i)>0.5,u(i)=-0.5;%如果M序列的值為"1"時,辨識的輸入信號取“-0.03” else u(i)=0.5;%當(dāng)M序列的值為"0"時,辨識的輸入信號取“0.03” end%小循環(huán)結(jié)束 y1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號做準備end%大循環(huán)結(jié)束，產(chǎn)生輸入信號u figure(1);%第1個圖形stem(u),gri

39、d on%以徑的形式顯示出輸入信號并給圖形加上網(wǎng)格z(2)=0;z(1)=0;%取z的前兩個初始值為零for k=3:L;%循環(huán)變量從3到25 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%給出理想的辨識輸出采樣信號 end%RLS遞推最小二乘辨識c0=0.001 0.001 0.001 0.001'%直接給出被辨識參數(shù)的初始值,即一個充分小的實向量p0=106*eye(4,4);%直接給出初始狀態(tài)P0，即一個充分大的實數(shù)單位矩陣E=0.000000005;%相對誤差E=0.000000005c=c0,zeros(4,99);%被辨識參數(shù)矩

40、陣的初始值及大小e=zeros(4,100);%相對誤差的初始值及大小for k=3:L; %開始求K h1=-z(k-1),-z(k-2),u(k-1),u(k-2)' x=h1'*p0*h1+1; x1=inv(x); %開始求K(k) k1=p0*h1*x1;%求出K的值 d1=z(k)-h1'*c0; c1=c0+k1*d1;%求被辨識參數(shù)c e1=c1-c0;%求參數(shù)當(dāng)前值與上一次的值的差值 e2=e1./c0;%求參數(shù)的相對變化 e(:,k)=e2; %把當(dāng)前相對變化的列向量加入誤差矩陣的最后一列 c0=c1;%新獲得的參數(shù)作為下一次遞推的舊參數(shù) c(:,k

41、)=c1;%把辨識參數(shù)c 列向量加入辨識參數(shù)矩陣的最后一列 p1=p0-k1*k1'*h1'*p0*h1+1;%求出 p(k)的值 p0=p1;%給下次用 if e2<=E break;%若參數(shù)收斂滿足要求，終止計算 end%小循環(huán)結(jié)束end%大循環(huán)結(jié)束c%顯示被辨識參數(shù)e%顯示辨識結(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個圖形i=1:L;%橫坐標從1到15plot(i,a1,&#

42、39;r',i,a2,':',i,b1,'g',i,b2,':') %畫出a1，a2，b1，b2的各次辨識結(jié)果title('Parameter Identification with Recursive Least Squares Method')%圖形標題figure(3); %第3個圖形i=1:L; %橫坐標從1到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %畫出a1，a2，b1，b2的各次辨識結(jié)果的收

43、斂情況title('Identification Precision') %圖形標題附錄6 遞推最小二乘法程序測試結(jié)果c = Columns 1 through 8 0.0010 0 0.0010 -0.4987 -1.2329 -1.5000 -1.5000 -1.5000 0.0010 0 0.0010 0.0010 -0.2350 0.7000 0.7000 0.7000 0.0010 0 0.2510 1.2503 1.0668 1.0000 1.0000 1.0000 0.0010 0 -0.2490 0.7503 0.5668 0.5000 0.5000 0.5000 Columns 9 through 16 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 0 0.7000 0.7000 0.7000 0.7000 0.7000 0.7000 0.7000 0 1.0000 1.0000 1.0000 1.0000 1.0000