系統(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)用，并通過(guò)實(shí)例分析說(shuō)明了最小二乘法應(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 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)識(shí)自然、利用自然的能力越來(lái)越強(qiáng)，對(duì)于未 知對(duì)象的探索也越來(lái)越深入.我們所研究的對(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ī) 律均一無(wú)所知的話，則把這樣的研究對(duì)象稱之為

4、“黑箱” 研究灰箱和黑箱時(shí)，將 研究的對(duì)象看作是一個(gè)系統(tǒng)，通過(guò)建立該系統(tǒng)的模型，對(duì)模型參數(shù)進(jìn)行辨識(shí)來(lái)確 定該系統(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í)的、廣泛的意義1. 1系統(tǒng)辨識(shí)簡(jiǎn)介系統(tǒng)辨識(shí)是根據(jù)系統(tǒng)的輸入輸出時(shí)間函數(shù)來(lái)確定描述系統(tǒng)行為的數(shù)學(xué)模型。 現(xiàn)代控制理論中的一個(gè)分支。通過(guò)辨識(shí)建立數(shù)學(xué)模型的目的是估計(jì)表征系統(tǒng)行為 的重要參數(shù)，建立一個(gè)能模仿真實(shí)系統(tǒng)行為的模型，用當(dāng)前可測(cè)量的系統(tǒng)的輸入 和輸出預(yù)測(cè)系統(tǒng)輸出的未來(lái)演變，以及設(shè)計(jì)控制器。對(duì)系統(tǒng)進(jìn)行分析的主要問(wèn)題 是根據(jù)輸入時(shí)間函數(shù)和系統(tǒng)的特性來(lái)

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

6、的目的是至關(guān)重要的。 它對(duì)模型類（模型結(jié)構(gòu)）、輸入信號(hào)和等價(jià)準(zhǔn)則的選擇都有很大的影響。通過(guò)辨 識(shí)建立數(shù)學(xué)模型通常有四個(gè)目的。 估計(jì)具有特定物理意義的參數(shù)有些表征系統(tǒng)行為的重要參數(shù)是難以直接測(cè)量的，例如在生理、生態(tài)、環(huán)境、經(jīng)濟(jì)等系統(tǒng)中就常有這種情況。這就需要 通過(guò)能觀測(cè)到的輸入輸出數(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ù)測(cè) 這是辨識(shí)的一個(gè)重要應(yīng)用方面，其目的是用迄今為止系統(tǒng)的可測(cè) 量的輸入和輸出去預(yù)測(cè)系統(tǒng)輸出的未來(lái)的演變。例如最常

7、見(jiàn)的氣象預(yù)報(bào)，洪水預(yù) 報(bào)，其他如太陽(yáng)黑子預(yù)報(bào)，市場(chǎng)價(jià)格的預(yù)測(cè)，河流污染物含量的預(yù)測(cè)等。預(yù)測(cè)模型辨識(shí)的等價(jià)準(zhǔn)則主要是使預(yù)測(cè)誤差平方和最小。只要預(yù)測(cè)誤差小就是好的預(yù)測(cè) 模型，對(duì)模型的結(jié)構(gòu)及參數(shù)則很少再有其他要求。這時(shí)辨識(shí)的準(zhǔn)則和模型應(yīng)用的 目的是一致的，因此可以得到較好的預(yù)測(cè)模型。 控制 為了設(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í)最小二乘方法簡(jiǎn)介對(duì)工程實(shí)踐中測(cè)得的數(shù)據(jù)進(jìn)行理論分析，用恰當(dāng)?shù)暮瘮?shù)去模擬數(shù)據(jù)原型是一 類十分重要的問(wèn)題，最常用的逼近原則是讓實(shí)測(cè)數(shù)

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

9、信號(hào)一般總是可以測(cè)量 的，而且過(guò)程的動(dòng)態(tài)特性必然表現(xiàn)在這些輸入輸出數(shù)據(jù)中，那么就可以利用輸入輸出數(shù)據(jù)所提供的信息來(lái)建立過(guò)程的數(shù)學(xué)模型。這種建模方法就稱為系統(tǒng)辨識(shí)。把辨識(shí)建模稱作“黑箱建?！薄?.1.2 最小二乘法系統(tǒng)辨識(shí)結(jié)構(gòu)：本文把待辨識(shí)的過(guò)程看作“黑箱”。只考慮過(guò)程的輸入輸出特性，而不強(qiáng)調(diào) 過(guò)程的內(nèi)部機(jī)理。圖中，輸入u(k)和輸出z(k)是可以觀測(cè)的；G是系統(tǒng)模型，用來(lái)描述系統(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è)一個(gè)隨機(jī)序列'z(k)，kE (1,2,丄)的均值是參數(shù)日的線性函數(shù)： Etz(k)'hT(k戶，其中h(k)是可測(cè)的數(shù)據(jù)向量，那么利用隨機(jī)序列的一個(gè)實(shí)現(xiàn)， 使準(zhǔn)則函數(shù)：L2J(R 二' z(k) -hT(kp 心(式2-2)達(dá)到極小的參數(shù)估計(jì)值二稱作二的最小二乘估計(jì)。最小二乘格式：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廣義最小二乘遞推算法的計(jì)算步驟：1.給定初始條件?(0)=；(充分小的實(shí)向量)Pf(0) = a2l

13、(a為充分大的數(shù))?e(0)= 0fe(O)=l2利用式乙(kC(zA)z(k) (k)二 C(z 丄)z(k)計(jì)算 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戸遞推計(jì)算?(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)計(jì)算鍬);6根據(jù) he(k) =e(k 1),，e(k - nc)廣來(lái)構(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)行迭代計(jì)算，直至獲得滿意的辨識(shí)結(jié)果3工程實(shí)例3.1典型系統(tǒng)建模以某微循環(huán)流體系統(tǒng)模型的參數(shù)辨識(shí)為例我們已經(jīng)得到該系統(tǒng)模型的差分 方程形式，取特定點(diǎn)的壓力波作為模型的輸入，以另一點(diǎn)的壓

15、力波作為模型的輸 出由于我們采集的數(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)的差分方程為：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語(yǔ)言辨識(shí) 系統(tǒng)中的未知參數(shù)3、a2、“、b2。要求：用參數(shù)的真值利用差分方程求出z(k) 作為測(cè)量值，：(k

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

17、相同的步長(zhǎng)條件下，精度大體相同。算法優(yōu)化方案：（1）使用M序列（具有近似白噪聲的性質(zhì)）為輸入信號(hào)；（2）增加數(shù)據(jù)長(zhǎng)度去L；（3）減小噪聲信號(hào)v（k）的方差。4.1廣義最小二乘遞推算法的的 MATLA仿真結(jié)果及分析（1）、輸入選用題目給出的30個(gè)隨機(jī)數(shù)，即數(shù)據(jù)長(zhǎng)度去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.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ī)序列，辨識(shí)步長(zhǎng)均為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.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ī)序列得到的辨識(shí)效果更 好。噪聲均值相等時(shí)，方差越大，辨識(shí)效果越差，反之，方差越小辨識(shí)效果越 好?？梢酝ㄟ^(guò)增加步長(zhǎng)的方法提高辨識(shí)精度。下面給出以M序列作為輸入，噪聲均值為零，方差為 0.01的隨機(jī)序列，數(shù) 據(jù)長(zhǎng)度

20、取L=30,得到的變化曲線圖：下面給出以M序列作為輸入，噪聲均值為零，方差為 0.01的不相關(guān)隨機(jī)序 列，數(shù)據(jù)長(zhǎng)度去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)辨識(shí)理論及應(yīng)用，國(guó)防工業(yè)出版社，2006年2 方崇智，蕭德云，過(guò)程辨識(shí)，清華大學(xué)出版社，2002年3 賈秋玲，袁冬莉，欒云鳳，基于 MATLAB7.x/Simulink/Stateflow 系統(tǒng)仿真、分析及設(shè)計(jì)，西北工業(yè)大學(xué)

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

22、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序列的值為"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=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的前兩個(gè)初始值為零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); %輸出采樣信號(hào)(測(cè)量值)end%RGL曠義最小二乘辨識(shí)c0=0.0001 0.00

24、01 0.0001 0.0001' %直接給出被辨識(shí)參數(shù)的初始值,即一個(gè)充分小的實(shí)向量pf0=10A6*eye(4,4); %直接給出初始狀態(tài)P0,即一個(gè)充分大的實(shí)數(shù)單位矩陣ce0=0.001 0.001'pe0=eye(2,2);c=c0,zeros(4,299); %被辨識(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);

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

26、數(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'*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í)參數(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); %第二個(gè)圖形畫出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 的各次辨識(shí)結(jié)果title(' 參數(shù)變化曲線')%圖形標(biāo)題figure(3); %第三個(gè)圖形i=1:300; % 橫坐標(biāo)從1到300 plot(i,ea1,'r',i,ea2,'k:',i,eb1,'b',i,eb2,'k:')%畫出 al，a2,辨識(shí)結(jié)果的收斂情況ti

29、tle('誤差曲線')%圖形標(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(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+v(k)，試?yán)眠f推最小二乘法辨識(shí)參數(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