辨識(shí)MATLAB編程：二階系統(tǒng)與漸消記憶實(shí)現(xiàn))
簡介這是一篇面向控制系統(tǒng)、信號(hào)處理與機(jī)器學(xué)習(xí)等領(lǐng)域?qū)W習(xí)者的最小二乘參數(shù)辨識(shí)程序說明文檔。資源共包含一個(gè)docx格式文檔壓縮包大小七十八KB內(nèi)容圍繞遞推最小二乘算法和漸消記憶最小二乘算法展開先構(gòu)造一個(gè)極點(diǎn)位于單位圓內(nèi)的穩(wěn)定系統(tǒng)生成五百點(diǎn)M序列作為輸入并疊加白噪聲擾動(dòng)隨后通過增益矩陣與估計(jì)方差矩陣的遞推更新實(shí)現(xiàn)系統(tǒng)參數(shù)的在線估計(jì)與辨識(shí)同時(shí)借助繪圖功能將參數(shù)估計(jì)過程可視化。文檔還針對漸消記憶法給出關(guān)鍵改進(jìn)通過引入遺忘因子削弱舊數(shù)據(jù)影響使估計(jì)結(jié)果更接近真實(shí)值其中包含完整的遞推公式、參數(shù)初值設(shè)定、循環(huán)迭代流程以及仿真圖形輸出便于讀者直接對照學(xué)習(xí)或改造成其他辨識(shí)對象。已有兩千一百五十八人學(xué)習(xí)下載適合需要掌握最小二乘類參數(shù)辨識(shí)原理、遞推算法實(shí)現(xiàn)及仿真驗(yàn)證方法的初學(xué)者和相關(guān)技術(shù)開發(fā)人員。1. 最小二乘參數(shù)辨識(shí)一套能直接改著用的 MATLAB 遞推程序做控制系統(tǒng)或信號(hào)處理的人大概率都遇到過這種事系統(tǒng)結(jié)構(gòu)大致清楚但里面的系數(shù)拿不準(zhǔn)。傳遞函數(shù)寫不出來狀態(tài)空間更是無從下手。這時(shí)候最直接的辦法就是給系統(tǒng)一個(gè)已知的激勵(lì)信號(hào)記錄輸出再用參數(shù)辨識(shí)把模型里的未知系數(shù)一點(diǎn)一點(diǎn)“逼”出來。這篇筆記要拆的是一份操作型很強(qiáng)的 MATLAB 程序——它用最小二乘遞推算法RLS辨識(shí)一個(gè)二階系統(tǒng)的四個(gè)參數(shù)然后給出漸消記憶最小二乘的改進(jìn)版本順帶把 M 序列激勵(lì)、白噪聲擾動(dòng)、增益矩陣更新這些細(xì)節(jié)全部攤開來講。適合剛接觸辨識(shí)理論但想盡快跑通一個(gè)完整例子的從業(yè)者也適合那些已經(jīng)在用工具箱、但想搞明白內(nèi)部迭代邏輯的人。這程序里辨識(shí)的對象是Z(k2)1.8Z(k1)-0.9Z(k)0.75u(k1)0.5u(k)v(k)。四個(gè)待辨識(shí)參數(shù)分別是 -1.8、0.9、0.75、0.5注意程序里 h 向量帶負(fù)號(hào)所以實(shí)際迭代目標(biāo)是 -1.8 和 0.9 這兩個(gè)帶符號(hào)的量。極點(diǎn)都在單位圓內(nèi)系統(tǒng)穩(wěn)定所以數(shù)據(jù)是收斂可用的。M 序列做輸入、白噪聲做擾動(dòng)這是辨識(shí)里最常見的測試組合。下面從算法原理到每一段代碼怎么改一步步來過。2. 遞推最小二乘的數(shù)學(xué)內(nèi)核為什么是 P、K、Theta 三個(gè)量在滾2.1 最小二乘的核心思想與遞推化的必要性最小二乘參數(shù)辨識(shí)的基本思路很簡單把系統(tǒng)寫成 y(k) h(k) * θ v(k) 的形式其中 h(k) 是回歸向量由過去的輸入輸出組成θ 是待辨識(shí)參數(shù)向量v(k) 是噪聲。我們的目標(biāo)是找一組 θ讓殘差平方和最小。批量最小二乘可以直接用矩陣求逆算出 θ(HH)^(-1)HZ但問題是如果數(shù)據(jù)是一批一批到的或者系統(tǒng)參數(shù)本身可能緩慢漂移每次都重新做一次矩陣求逆計(jì)算量會(huì)隨數(shù)據(jù)量增長而爆炸而且老數(shù)據(jù)會(huì)把新數(shù)據(jù)的信息“稀釋”掉。遞推最小二乘就是為了解決這個(gè)矛盾。它把批處理公式拆成一個(gè)迭代式子每來一個(gè)新數(shù)據(jù)就用增益矩陣 K(k) 對當(dāng)前的參數(shù)估計(jì) Theta(k-1) 做一個(gè)修正。修正量等于 K(k) 乘以預(yù)測誤差實(shí)際輸出減去當(dāng)前參數(shù)預(yù)測的輸出。這個(gè)過程中有三個(gè)量在同步滾參數(shù)估計(jì)值 Theta、誤差協(xié)方差矩陣 P、增益矩陣 K。P 矩陣本質(zhì)上記錄了參數(shù)估計(jì)的不確定度K 則是把 P 和當(dāng)前回歸向量 h 組合起來算出的修正步長。三者配合理論上可以從任意初值出發(fā)逐步逼近真實(shí)參數(shù)。這個(gè)程序用的系統(tǒng)方程是 Z(k2)1.8Z(k1)-0.9Z(k)0.75u(k1)0.5u(k)v(k)。把它改寫成最小二乘的標(biāo)準(zhǔn)形式回歸向量 h(k)[-z(k-1); -z(k-2); u(k-1); u(k-2)]程序中循環(huán)從 i3 開始所以 h 里取的是 z(i-1)、z(i-2)、M(i-1)、M(i-2)對應(yīng)前面的 z(k1)、z(k)、u(k1)、u(k)待辨識(shí)參數(shù) θ[a1; a2; b1; b2]其中 a1-1.8、a20.9、b10.75、b20.5。注意這里回歸向量里的 z 和 u 都帶了負(fù)號(hào)或正號(hào)符號(hào)問題在后面看代碼時(shí)要特別留神——很多人自己寫的時(shí)候就是在這里把符號(hào)搞反導(dǎo)致辨識(shí)出來的參數(shù)符號(hào)全錯(cuò)。2.2 RLS 迭代公式與程序中的變量一一對應(yīng)RLS 的標(biāo)準(zhǔn)迭代公式有三個(gè)增益矩陣K(k) P(k-1) * h(k) / (h(k) * P(k-1) * h(k) λ)。λ 是遺忘因子標(biāo)準(zhǔn) RLS 取 λ1漸消記憶法取 λ1如 0.998。協(xié)方差更新P(k) (I - K(k) * h(k)) * P(k-1) / λ。參數(shù)更新Theta(k) Theta(k-1) K(k) * (z(k) - h(k) * Theta(k-1))。程序里的P100*eye(4)對應(yīng) P 的初始值Theta(:,1)[1;1;1;1]對應(yīng)參數(shù)初值K[1;1;1;1]只是初始化一個(gè)列向量在循環(huán)里會(huì)被真正計(jì)算的值覆蓋。h[-z(i-1);-z(i-2);M(i-1);M(i-2)]就是當(dāng)前時(shí)刻的回歸向量。循環(huán)體里三行代碼和上面三個(gè)公式一一對應(yīng)。P 的初值取得比較大100 倍單位陣表示一開始對參數(shù)估計(jì)非常不確定這樣 K 的初始修正步長比較大參數(shù)能快速從初值 [1;1;1;1]這是一個(gè)故意給偏的初值收斂到真值附近。這里有個(gè)值得注意的細(xì)節(jié)程序里 hPh 1 中的 1 是噪聲方差假設(shè)。如果實(shí)際系統(tǒng)噪聲方差不是 1這個(gè) 1 應(yīng)該相應(yīng)調(diào)整。但在這個(gè)例子里vrandn(1,N) 產(chǎn)生的白噪聲方差恰好是 1所以不用改。很多人直接把工具箱里的 RLS 拿過來用根本沒意識(shí)到這個(gè)數(shù)字和噪聲方差掛鉤等換了系統(tǒng)換了噪聲辨識(shí)結(jié)果變差了卻找不到原因——這就是其中一條隱蔽的坑。3. 程序逐段拆解M 序列、白噪聲、觀測數(shù)據(jù)和遞推主循環(huán)3.1 M 序列產(chǎn)生移位寄存器與異或反饋的配合程序開頭這段是產(chǎn)生 M 序列偽隨機(jī)二進(jìn)制序列作為輸入激勵(lì)x[0 1 0 1 1 0 1 1 1]; % 寄存器初值9級(jí) n503; % 脈沖總數(shù) N500; % 實(shí)際使用的點(diǎn)數(shù) for i1:n tempxor(x(4),x(9)); % 第4級(jí)與第9級(jí)異或作為反饋 M(i)x(9); % 輸出取第9級(jí) for j9:-1:2 x(j)x(j-1); % 寄存器右移一位 end x(1)temp; % 反饋值寫入第1級(jí) end邏輯說明這是一個(gè) 9 級(jí)線性反饋移位寄存器LFSR反饋多項(xiàng)式取第 4 級(jí)和第 9 級(jí)的異或。每來一個(gè)時(shí)鐘脈沖寄存器內(nèi)容右移一位最高位丟棄最低位補(bǔ)上反饋值。輸出的 M 序列取的是第 9 級(jí)的值。這個(gè)序列是周期性的周期為 2^9 - 1 511。程序產(chǎn)生 503 個(gè)點(diǎn)取前 500 個(gè)使用所以實(shí)際用的是一整個(gè)周期內(nèi)的連續(xù)片段。參數(shù)說明寄存器初值 x 不能全零否則序列永遠(yuǎn)輸出 0這一點(diǎn)要特別注意。n 設(shè)為 503 而 N 設(shè)為 500多出來的 3 個(gè)點(diǎn)是留給后面生成觀測數(shù)據(jù)時(shí)用的——觀測數(shù)據(jù) z 需要從 i3 開始循環(huán)所以輸入序列必須多準(zhǔn)備幾個(gè)點(diǎn)才能讓下標(biāo)對齊。反饋抽頭選擇第 4 和第 9 級(jí)不是隨意的9 級(jí)寄存器用這對抽頭能保證序列周期達(dá)到最大長度 511。如果你要改寄存器級(jí)數(shù)抽頭位置也要跟著變否則周期會(huì)縮短激勵(lì)信號(hào)的頻譜覆蓋就不夠均勻。3.2 白噪聲與觀測數(shù)據(jù)生成驗(yàn)證模型時(shí)的關(guān)鍵路徑vrandn(1,N); % 方差為1的白噪聲 zzeros(N2,1); z(1)-1; z(2)0; for i3:N2 z(i)1.8*z(i-1)-0.9*z(i-2)0.75*M(i-1)0.5*M(i-2)v(i-2); end邏輯說明z 是系統(tǒng)的觀測輸出序列。z(1)-1、z(2)0 是兩個(gè)初始狀態(tài)不是為了湊數(shù)而是為了讓循環(huán)從 i3 開始時(shí)有前兩個(gè)輸出值可用。第 i 個(gè)觀測值由前兩個(gè)輸出、前兩個(gè)輸入和當(dāng)前噪聲共同決定。公式里的系數(shù) 1.8、-0.9、0.75、0.5 就是我們要辨識(shí)的真值。注意噪聲項(xiàng)是 v(i-2)和輸出 z(i) 對齊——因?yàn)榛貧w向量里的輸入輸出都是滯后項(xiàng)噪聲也要對應(yīng)滯后兩個(gè)時(shí)刻才合理。參數(shù)說明randn(1,N) 產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)均值 0 方差 1正好匹配 RLS 公式里分母上的“1”。如果你把噪聲方差調(diào)大比如改成 2*randn(1,N)那么 RLS 公式里的那個(gè) 1 也要相應(yīng)改成 4方差值否則增益矩陣的步長和噪聲水平不匹配參數(shù)估計(jì)的波動(dòng)會(huì)明顯變大。程序里 z 的長度是 N2比 M 序列多兩個(gè)點(diǎn)也是為了讓循環(huán)下標(biāo)從 3 開始后能恰好用到 M 序列的全部 500 個(gè)點(diǎn)。3.3 遞推主循環(huán)P、K、Theta 三個(gè)量的更新順序是成敗關(guān)鍵P100*eye(4); % 估計(jì)方差初值給大 Pstorezeros(4,N1); Pstore(:,1)[P(1,1),P(2,2),P(3,3),P(4,4)]; Thetazeros(4,N1); % 參數(shù)估計(jì)值存儲(chǔ) Theta(:,1)[1;1;1;1]; % 故意給偏的參數(shù)初值 K[1;1;1;1]; % 增益矩陣初值占位 for i3:N2 h[-z(i-1);-z(i-2);M(i-1);M(i-2)]; % 回歸向量 KP*h*inv(h*P*h1); % 增益矩陣 Theta(:,i-1)Theta(:,i-2)K*(z(i)-h*Theta(:,i-2)); % 參數(shù)更新 P(eye(4)-K*h)*P; % 協(xié)方差更新 Pstore(:,i-1)[P(1,1),P(2,2),P(3,3),P(4,4)]; % 記錄方差對角元 end邏輯說明循環(huán)里三步依次是算 K、更新 Theta、更新 P順序是固定的——必須先用舊的 P 算 K再用 K 修正 Theta最后才更新 P。換順序會(huì)導(dǎo)致 P 和 K 不匹配參數(shù)收斂軌跡會(huì)亂掉。存儲(chǔ)數(shù)組 Pstore 和 Theta 的列數(shù)都是 N1從 1 到 501但循環(huán)體里 Theta(:,i-1) 的下標(biāo)范圍是 2 到 501所以初值放在第 1 列最后一個(gè)估計(jì)值放在第 501 列。畫圖時(shí) i1:N1 正好覆蓋全部估計(jì)點(diǎn)。參數(shù)說明P 的初值取 100eye(4)表示初始協(xié)方差為 100即參數(shù)不確定度很高。這個(gè)值如果取太小比如 1參數(shù)收斂會(huì)明顯變慢取太大比如 10^6則前幾步的參數(shù)估計(jì)會(huì)劇烈震蕩。Theta 初值取 [1;1;1;1] 離真值 [-1.8;0.9;0.75;0.5] 有不小距離這樣能直觀看到收斂過程。回歸向量 h 里的符號(hào)很關(guān)鍵因?yàn)椴罘址匠淌?z(k2)1.8z(k1)-0.9*z(k)...移到等式左邊后z 項(xiàng)的系數(shù)分別變成 -1.8 和 0.9所以 h 里對應(yīng)位置是 -z(i-1) 和 -z(i-2)這樣辨識(shí)出來的 θ(1) 和 θ(2) 會(huì)收斂到 -1.8 和 0.9。4. 漸消記憶最小二乘遺忘因子 0.998 背后的原理和改法4.1 為什么標(biāo)準(zhǔn) RLS 會(huì)有穩(wěn)態(tài)偏差標(biāo)準(zhǔn) RLS 的推導(dǎo)里P 矩陣會(huì)隨著迭代次數(shù)增加而不斷縮小。這意味著后期增益矩陣 K 的模越來越小參數(shù)的修正步長也隨之變小。理論上如果系統(tǒng)是線性定常的、噪聲是零均值的參數(shù)會(huì)收斂到真值附近但實(shí)際仿真中你往往會(huì)發(fā)現(xiàn)估計(jì)值和真值之間始終存在一個(gè)小的偏差而且這個(gè)偏差不隨迭代次數(shù)增加而消失。原因在于RLS 的代價(jià)函數(shù)對所有歷史數(shù)據(jù)等權(quán)重看待越老的數(shù)據(jù)積累越多在總代價(jià)里占比越大新數(shù)據(jù)對參數(shù)修正的貢獻(xiàn)相對被稀釋。如果系統(tǒng)本身有輕微時(shí)變或者噪聲不是純白噪聲比如有色噪聲這個(gè)偏差會(huì)更明顯。把系統(tǒng)輸出和參數(shù)估計(jì)代入程序畫出的圖上你會(huì)發(fā)現(xiàn)四條參數(shù)曲線都收斂了但每條曲線和真值之間都有一條肉眼可見的“縫”。這個(gè)縫的大小取決于噪聲水平、P 初值和數(shù)據(jù)長度。P 初值取 100、數(shù)據(jù)取 500 點(diǎn)時(shí)偏差通常較小但不為零——這正是程序注釋里寫的“辨識(shí)結(jié)果會(huì)與其精確值始終存在著偏差”。4.2 改成漸消記憶法的關(guān)鍵兩個(gè)語句而已漸消記憶法的思路是給歷史數(shù)據(jù)加一個(gè)指數(shù)衰減的權(quán)重越老的數(shù)據(jù)權(quán)重越小這樣新數(shù)據(jù)能在代價(jià)函數(shù)中占據(jù)更大比例從而持續(xù)修正參數(shù)。實(shí)現(xiàn)上只需要改兩處% 原語句 KP*h*inv(h*P*h1); P(eye(4)-K*h)*P; % 改為漸消記憶法遺忘因子 λ0.998 KP*h*inv(h*P*h0.998); P(eye(4)-K*h)*P/0.998;邏輯說明分母上的 1 改成 0.998對應(yīng)遺忘因子 λ。P 更新時(shí)除以 0.998即乘以 1/0.998 ≈ 1.002效果是讓 P 矩陣不會(huì)像標(biāo)準(zhǔn) RLS 那樣快速收縮到接近零。這樣每步結(jié)束后 P 維持在一個(gè)“有彈性”的量級(jí)下一個(gè)數(shù)據(jù)點(diǎn)到來時(shí) K 仍然有足夠大的模參數(shù)修正步長不會(huì)衰減到微乎其微。綜合效果就是估計(jì)值能持續(xù)跟隨真值穩(wěn)態(tài)偏差明顯縮小。參數(shù)說明遺忘因子 λ 的取值范圍通常是 0.951.0。λ 越接近 1記憶越長參數(shù)越穩(wěn)定但跟蹤越慢λ 越小比如 0.95記憶越短跟蹤時(shí)變參數(shù)的能力越強(qiáng)但穩(wěn)態(tài)波動(dòng)也越大。0.998 是一個(gè)比較保守的值適合參數(shù)基本不變但希望消除穩(wěn)態(tài)偏差的場景。如果你的系統(tǒng)參數(shù)是緩變的可以考慮取 0.980.99如果參數(shù)變化較快可能需要 0.950.97但這時(shí)要考慮噪聲放大的副作用。注意在這里噪聲方差假設(shè)從 1 變成了 λ0.998所以分母的含義不只是噪聲方差而是噪聲方差與遺忘因子的組合——這個(gè)細(xì)節(jié)很多人會(huì)忽略。4.3 兩種算法的結(jié)果對比怎么看對比標(biāo)準(zhǔn) RLS 和漸消記憶法的參數(shù)估計(jì)曲線最直觀的差異在穩(wěn)態(tài)段。標(biāo)準(zhǔn) RLS 的曲線雖然收斂但每條曲線和真值之間始終有偏差漸消記憶法的曲線幾乎貼近真值線。P 矩陣的對角元變化也解釋了這一點(diǎn)標(biāo)準(zhǔn) RLS 的 P 對角元持續(xù)衰減到接近零漸消記憶法的 P 對角元穩(wěn)定在一個(gè)較小的非零值附近。這個(gè)差異在你需要高精度參數(shù)值時(shí)非常關(guān)鍵——比如你辨識(shí)出來的參數(shù)接下來要用于設(shè)計(jì)控制器偏差可能導(dǎo)致閉環(huán)極點(diǎn)偏出期望位置。5. 避坑指南從初值設(shè)定到噪聲匹配四個(gè)最容易翻車的地方5.1 現(xiàn)象P 初值給太小導(dǎo)致收斂極慢現(xiàn)象把 P100eye(4) 改成 Peye(4) 或 P0.1eye(4)參數(shù)估計(jì)曲線需要很長時(shí)間才能逼近真值500 個(gè)數(shù)據(jù)點(diǎn)可能不夠用。原因P 初值代表對參數(shù)初始不確定度的估計(jì)給小了等于默認(rèn)“初始參數(shù)已經(jīng)很準(zhǔn)”增益矩陣 K 的初始步長就小參數(shù)只能慢慢爬過去。解決保持相對較大的初值如 100 倍單位陣或者根據(jù)參數(shù)可能的變化范圍手動(dòng)估計(jì)一個(gè)合適的數(shù)量級(jí)。我一般的做法是先跑一次標(biāo)準(zhǔn) RLS 看 P 對角元收斂到多少再反推一個(gè)合理的初值一般取收斂值的 10100 倍。5.2 現(xiàn)象遺忘因子太小導(dǎo)致參數(shù)估計(jì)波動(dòng)大現(xiàn)象把遺忘因子從 0.998 改成 0.95參數(shù)曲線雖然在真值附近波動(dòng)但波動(dòng)的幅度明顯變大看起來“毛刺”很多。原因遺忘因子越小有效數(shù)據(jù)窗口越短等效于用更少的數(shù)據(jù)做估計(jì)方差自然增大。解決先試 0.990.998 這個(gè)區(qū)間確認(rèn)穩(wěn)態(tài)偏差確實(shí)被消除后再嘗試適當(dāng)減小。如果偏差還在先檢查是不是噪聲方差沒匹配而不是一味調(diào)小遺忘因子。我自己的習(xí)慣是只在對參數(shù)有明確時(shí)變需求時(shí)才把遺忘因子調(diào)到 0.98 以下否則盡量保持在 0.99 附近。5.3 現(xiàn)象回歸向量的符號(hào)寫反辨識(shí)出的參數(shù)符號(hào)全錯(cuò)現(xiàn)象把 h[-z(i-1);-z(i-2);M(i-1);M(i-2)] 寫成 h[z(i-1);z(i-2);M(i-1);M(i-2)]結(jié)果 θ(1) 收斂到 1.8 而不是 -1.8θ(2) 收斂到 -0.9 而不是 0.9。原因差分方程移項(xiàng)后系數(shù)符號(hào)變了但回歸向量沒有同步調(diào)整。解決先把差分方程寫成標(biāo)準(zhǔn)回歸形式 y(k)h(k)*θv(k)逐項(xiàng)確認(rèn)符號(hào)再寫代碼。這一步最基礎(chǔ)但最容易錯(cuò)而且錯(cuò)的隱蔽——最終參數(shù)的絕對值是對的正負(fù)號(hào)卻是反的。5.4 現(xiàn)象M 序列初值全零導(dǎo)致激勵(lì)信號(hào)為零辨識(shí)完全失敗現(xiàn)象寄存器初值改成 x[0 0 0 0 0 0 0 0 0]程序運(yùn)行后參數(shù)曲線完全不動(dòng)或者發(fā)散到異常值。原因全零初值的 LFSR 永遠(yuǎn)輸出 0系統(tǒng)完全沒有被激勵(lì)觀測數(shù)據(jù)里不含任何參數(shù)信息。解決確保寄存器初值至少有一個(gè) 1最省事的就是保持程序原來的初值 [0 1 0 1 1 0 1 1 1]。如果你自己重新生成 M 序列記得在代碼里加一個(gè)判斷如果 sum(x)0就報(bào)錯(cuò)提示重新初始化。5.5 現(xiàn)象噪聲方差與 RLS 公式中的常數(shù)不匹配估計(jì)值偏離真值現(xiàn)象把噪聲改成 v2*randn(1,N)但 RLS 公式里的常數(shù)保持 1最后參數(shù)估計(jì)偏離真值且曲線波動(dòng)變大。原因公式里的常數(shù)同時(shí)承擔(dān)了噪聲方差設(shè)定和遺忘因子兩個(gè)角色噪聲實(shí)際方差和假設(shè)值不一致時(shí)增益矩陣的步長不是最優(yōu)的。解決如果你改了噪聲方差記得同步修改inv(h*P*h常數(shù))里的常數(shù)改成噪聲方差的數(shù)值。這是最容易被忽略的一個(gè)參數(shù)因?yàn)樗卦诰仃嚽竽娴谋磉_(dá)式里不是單獨(dú)一句話。6. 把辨識(shí)結(jié)果落地殘差校驗(yàn) 利用帶遺忘因子的 RLS 做在線跟蹤當(dāng)你把漸消記憶法的程序跑通、參數(shù)曲線也收斂到真值附近之后還有一個(gè)值得做的驗(yàn)證步驟檢查模型輸出和實(shí)際觀測之間的殘差。做法很簡單把辨識(shí)出的參數(shù)代回系統(tǒng)方程用同一段輸入重新仿真一遍然后算殘差% 辨識(shí)完成后用估計(jì)參數(shù)計(jì)算模型輸出 z_est zeros(N2,1); z_est(1) -1; z_est(2) 0; for i3:N2 z_est(i) -Theta(1,end)*z_est(i-1) - Theta(2,end)*z_est(i-2) ... Theta(3,end)*M(i-1) Theta(4,end)*M(i-2); end residual z(3:N2) - z_est(3:N2); % 殘差序列 figure(2) subplot(2,1,1) plot(1:N, residual) title(模型殘差序列) subplot(2,1,2) histogram(residual, 30) title(殘差直方圖)邏輯說明殘差等于實(shí)際輸出減模型輸出。如果辨識(shí)效果好殘差應(yīng)該接近白噪聲——均值接近零、直方圖呈鐘形。如果殘差里還有明顯的周期性或趨勢說明模型結(jié)構(gòu)沒對比如階次不夠高或者輸入信號(hào)沒有充分激勵(lì)出系統(tǒng)的所有模態(tài)。參數(shù)說明注意這里模型輸出用 Theta(1,end)、Theta(2,end) 取的是最后一列的估計(jì)值。Theta(1,end) 是辨識(shí)出的 -a1因?yàn)榛貧w向量里 h(1)-z(i-1)所以計(jì)算模型輸出時(shí)要寫成-Theta(1,end)*z_est(i-1)。這個(gè)符號(hào)細(xì)節(jié)和前面回歸向量里的符號(hào)是呼應(yīng)的別在這里搞反。一個(gè)我自己常用的習(xí)慣是在線辨識(shí)場景里不只要畫參數(shù)曲線還要同時(shí)畫 P 矩陣的跡或?qū)窃鸵驗(yàn)樗苤庇^反映“當(dāng)前參數(shù)可信度”。如果跡在迭代后期還在持續(xù)下降說明系統(tǒng)還沒進(jìn)入穩(wěn)態(tài)如果跡已經(jīng)穩(wěn)住但有波動(dòng)說明遺忘因子在起作用。從那以后我每次跑辨識(shí)實(shí)驗(yàn)都會(huì)強(qiáng)制走一遍“觀察 P 跡 → 檢查殘差 → 確認(rèn)符號(hào)”這三步再把結(jié)果用于后續(xù)設(shè)計(jì)。這一步看著多花幾分鐘實(shí)際能擋掉大量后期的低級(jí)返工。希望這份程序和拆解思路能幫到你尤其是那些剛開始接觸遞推辨識(shí)、想盡快把理論落到代碼上的朋友——按上面的步驟改參數(shù)、看曲線、檢查殘差你很快就能摸清這套算法的手感。本文還有配套的精品資源點(diǎn)擊獲取 遞推最小二乘參數(shù)辨識(shí)MATLAB編程：二階系統(tǒng)與漸消記憶實(shí)現(xiàn)
簡介這是一篇面向控制系統(tǒng)、信號(hào)處理與機(jī)器學(xué)習(xí)等領(lǐng)域?qū)W習(xí)者的最小二乘參數(shù)辨識(shí)程序說明文檔。資源共包含一個(gè)docx格式文檔壓縮包大小七十八KB內(nèi)容圍繞遞推最小二乘算法和漸消記憶最小二乘算法展開先構(gòu)造一個(gè)極點(diǎn)位于單位圓內(nèi)的穩(wěn)定系統(tǒng)生成五百點(diǎn)M序列作為輸入并疊加白噪聲擾動(dòng)隨后通過增益矩陣與估計(jì)方差矩陣的遞推更新實(shí)現(xiàn)系統(tǒng)參數(shù)的在線估計(jì)與辨識(shí)同時(shí)借助繪圖功能將參數(shù)估計(jì)過程可視化。文檔還針對漸消記憶法給出關(guān)鍵改進(jìn)通過引入遺忘因子削弱舊數(shù)據(jù)影響使估計(jì)結(jié)果更接近真實(shí)值其中包含完整的遞推公式、參數(shù)初值設(shè)定、循環(huán)迭代流程以及仿真圖形輸出便于讀者直接對照學(xué)習(xí)或改造成其他辨識(shí)對象。已有兩千一百五十八人學(xué)習(xí)下載適合需要掌握最小二乘類參數(shù)辨識(shí)原理、遞推算法實(shí)現(xiàn)及仿真驗(yàn)證方法的初學(xué)者和相關(guān)技術(shù)開發(fā)人員。1. 最小二乘參數(shù)辨識(shí)一套能直接改著用的 MATLAB 遞推程序做控制系統(tǒng)或信號(hào)處理的人大概率都遇到過這種事系統(tǒng)結(jié)構(gòu)大致清楚但里面的系數(shù)拿不準(zhǔn)。傳遞函數(shù)寫不出來狀態(tài)空間更是無從下手。這時(shí)候最直接的辦法就是給系統(tǒng)一個(gè)已知的激勵(lì)信號(hào)記錄輸出再用參數(shù)辨識(shí)把模型里的未知系數(shù)一點(diǎn)一點(diǎn)“逼”出來。這篇筆記要拆的是一份操作型很強(qiáng)的 MATLAB 程序——它用最小二乘遞推算法RLS辨識(shí)一個(gè)二階系統(tǒng)的四個(gè)參數(shù)然后給出漸消記憶最小二乘的改進(jìn)版本順帶把 M 序列激勵(lì)、白噪聲擾動(dòng)、增益矩陣更新這些細(xì)節(jié)全部攤開來講。適合剛接觸辨識(shí)理論但想盡快跑通一個(gè)完整例子的從業(yè)者也適合那些已經(jīng)在用工具箱、但想搞明白內(nèi)部迭代邏輯的人。這程序里辨識(shí)的對象是Z(k2)1.8Z(k1)-0.9Z(k)0.75u(k1)0.5u(k)v(k)。四個(gè)待辨識(shí)參數(shù)分別是 -1.8、0.9、0.75、0.5注意程序里 h 向量帶負(fù)號(hào)所以實(shí)際迭代目標(biāo)是 -1.8 和 0.9 這兩個(gè)帶符號(hào)的量。極點(diǎn)都在單位圓內(nèi)系統(tǒng)穩(wěn)定所以數(shù)據(jù)是收斂可用的。M 序列做輸入、白噪聲做擾動(dòng)這是辨識(shí)里最常見的測試組合。下面從算法原理到每一段代碼怎么改一步步來過。2. 遞推最小二乘的數(shù)學(xué)內(nèi)核為什么是 P、K、Theta 三個(gè)量在滾2.1 最小二乘的核心思想與遞推化的必要性最小二乘參數(shù)辨識(shí)的基本思路很簡單把系統(tǒng)寫成 y(k) h(k) * θ v(k) 的形式其中 h(k) 是回歸向量由過去的輸入輸出組成θ 是待辨識(shí)參數(shù)向量v(k) 是噪聲。我們的目標(biāo)是找一組 θ讓殘差平方和最小。批量最小二乘可以直接用矩陣求逆算出 θ(HH)^(-1)HZ但問題是如果數(shù)據(jù)是一批一批到的或者系統(tǒng)參數(shù)本身可能緩慢漂移每次都重新做一次矩陣求逆計(jì)算量會(huì)隨數(shù)據(jù)量增長而爆炸而且老數(shù)據(jù)會(huì)把新數(shù)據(jù)的信息“稀釋”掉。遞推最小二乘就是為了解決這個(gè)矛盾。它把批處理公式拆成一個(gè)迭代式子每來一個(gè)新數(shù)據(jù)就用增益矩陣 K(k) 對當(dāng)前的參數(shù)估計(jì) Theta(k-1) 做一個(gè)修正。修正量等于 K(k) 乘以預(yù)測誤差實(shí)際輸出減去當(dāng)前參數(shù)預(yù)測的輸出。這個(gè)過程中有三個(gè)量在同步滾參數(shù)估計(jì)值 Theta、誤差協(xié)方差矩陣 P、增益矩陣 K。P 矩陣本質(zhì)上記錄了參數(shù)估計(jì)的不確定度K 則是把 P 和當(dāng)前回歸向量 h 組合起來算出的修正步長。三者配合理論上可以從任意初值出發(fā)逐步逼近真實(shí)參數(shù)。這個(gè)程序用的系統(tǒng)方程是 Z(k2)1.8Z(k1)-0.9Z(k)0.75u(k1)0.5u(k)v(k)。把它改寫成最小二乘的標(biāo)準(zhǔn)形式回歸向量 h(k)[-z(k-1); -z(k-2); u(k-1); u(k-2)]程序中循環(huán)從 i3 開始所以 h 里取的是 z(i-1)、z(i-2)、M(i-1)、M(i-2)對應(yīng)前面的 z(k1)、z(k)、u(k1)、u(k)待辨識(shí)參數(shù) θ[a1; a2; b1; b2]其中 a1-1.8、a20.9、b10.75、b20.5。注意這里回歸向量里的 z 和 u 都帶了負(fù)號(hào)或正號(hào)符號(hào)問題在后面看代碼時(shí)要特別留神——很多人自己寫的時(shí)候就是在這里把符號(hào)搞反導(dǎo)致辨識(shí)出來的參數(shù)符號(hào)全錯(cuò)。2.2 RLS 迭代公式與程序中的變量一一對應(yīng)RLS 的標(biāo)準(zhǔn)迭代公式有三個(gè)增益矩陣K(k) P(k-1) * h(k) / (h(k) * P(k-1) * h(k) λ)。λ 是遺忘因子標(biāo)準(zhǔn) RLS 取 λ1漸消記憶法取 λ1如 0.998。協(xié)方差更新P(k) (I - K(k) * h(k)) * P(k-1) / λ。參數(shù)更新Theta(k) Theta(k-1) K(k) * (z(k) - h(k) * Theta(k-1))。程序里的P100*eye(4)對應(yīng) P 的初始值Theta(:,1)[1;1;1;1]對應(yīng)參數(shù)初值K[1;1;1;1]只是初始化一個(gè)列向量在循環(huán)里會(huì)被真正計(jì)算的值覆蓋。h[-z(i-1);-z(i-2);M(i-1);M(i-2)]就是當(dāng)前時(shí)刻的回歸向量。循環(huán)體里三行代碼和上面三個(gè)公式一一對應(yīng)。P 的初值取得比較大100 倍單位陣表示一開始對參數(shù)估計(jì)非常不確定這樣 K 的初始修正步長比較大參數(shù)能快速從初值 [1;1;1;1]這是一個(gè)故意給偏的初值收斂到真值附近。這里有個(gè)值得注意的細(xì)節(jié)程序里 hPh 1 中的 1 是噪聲方差假設(shè)。如果實(shí)際系統(tǒng)噪聲方差不是 1這個(gè) 1 應(yīng)該相應(yīng)調(diào)整。但在這個(gè)例子里vrandn(1,N) 產(chǎn)生的白噪聲方差恰好是 1所以不用改。很多人直接把工具箱里的 RLS 拿過來用根本沒意識(shí)到這個(gè)數(shù)字和噪聲方差掛鉤等換了系統(tǒng)換了噪聲辨識(shí)結(jié)果變差了卻找不到原因——這就是其中一條隱蔽的坑。3. 程序逐段拆解M 序列、白噪聲、觀測數(shù)據(jù)和遞推主循環(huán)3.1 M 序列產(chǎn)生移位寄存器與異或反饋的配合程序開頭這段是產(chǎn)生 M 序列偽隨機(jī)二進(jìn)制序列作為輸入激勵(lì)x[0 1 0 1 1 0 1 1 1]; % 寄存器初值9級(jí) n503; % 脈沖總數(shù) N500; % 實(shí)際使用的點(diǎn)數(shù) for i1:n tempxor(x(4),x(9)); % 第4級(jí)與第9級(jí)異或作為反饋 M(i)x(9); % 輸出取第9級(jí) for j9:-1:2 x(j)x(j-1); % 寄存器右移一位 end x(1)temp; % 反饋值寫入第1級(jí) end邏輯說明這是一個(gè) 9 級(jí)線性反饋移位寄存器LFSR反饋多項(xiàng)式取第 4 級(jí)和第 9 級(jí)的異或。每來一個(gè)時(shí)鐘脈沖寄存器內(nèi)容右移一位最高位丟棄最低位補(bǔ)上反饋值。輸出的 M 序列取的是第 9 級(jí)的值。這個(gè)序列是周期性的周期為 2^9 - 1 511。程序產(chǎn)生 503 個(gè)點(diǎn)取前 500 個(gè)使用所以實(shí)際用的是一整個(gè)周期內(nèi)的連續(xù)片段。參數(shù)說明寄存器初值 x 不能全零否則序列永遠(yuǎn)輸出 0這一點(diǎn)要特別注意。n 設(shè)為 503 而 N 設(shè)為 500多出來的 3 個(gè)點(diǎn)是留給后面生成觀測數(shù)據(jù)時(shí)用的——觀測數(shù)據(jù) z 需要從 i3 開始循環(huán)所以輸入序列必須多準(zhǔn)備幾個(gè)點(diǎn)才能讓下標(biāo)對齊。反饋抽頭選擇第 4 和第 9 級(jí)不是隨意的9 級(jí)寄存器用這對抽頭能保證序列周期達(dá)到最大長度 511。如果你要改寄存器級(jí)數(shù)抽頭位置也要跟著變否則周期會(huì)縮短激勵(lì)信號(hào)的頻譜覆蓋就不夠均勻。3.2 白噪聲與觀測數(shù)據(jù)生成驗(yàn)證模型時(shí)的關(guān)鍵路徑vrandn(1,N); % 方差為1的白噪聲 zzeros(N2,1); z(1)-1; z(2)0; for i3:N2 z(i)1.8*z(i-1)-0.9*z(i-2)0.75*M(i-1)0.5*M(i-2)v(i-2); end邏輯說明z 是系統(tǒng)的觀測輸出序列。z(1)-1、z(2)0 是兩個(gè)初始狀態(tài)不是為了湊數(shù)而是為了讓循環(huán)從 i3 開始時(shí)有前兩個(gè)輸出值可用。第 i 個(gè)觀測值由前兩個(gè)輸出、前兩個(gè)輸入和當(dāng)前噪聲共同決定。公式里的系數(shù) 1.8、-0.9、0.75、0.5 就是我們要辨識(shí)的真值。注意噪聲項(xiàng)是 v(i-2)和輸出 z(i) 對齊——因?yàn)榛貧w向量里的輸入輸出都是滯后項(xiàng)噪聲也要對應(yīng)滯后兩個(gè)時(shí)刻才合理。參數(shù)說明randn(1,N) 產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)均值 0 方差 1正好匹配 RLS 公式里分母上的“1”。如果你把噪聲方差調(diào)大比如改成 2*randn(1,N)那么 RLS 公式里的那個(gè) 1 也要相應(yīng)改成 4方差值否則增益矩陣的步長和噪聲水平不匹配參數(shù)估計(jì)的波動(dòng)會(huì)明顯變大。程序里 z 的長度是 N2比 M 序列多兩個(gè)點(diǎn)也是為了讓循環(huán)下標(biāo)從 3 開始后能恰好用到 M 序列的全部 500 個(gè)點(diǎn)。3.3 遞推主循環(huán)P、K、Theta 三個(gè)量的更新順序是成敗關(guān)鍵P100*eye(4); % 估計(jì)方差初值給大 Pstorezeros(4,N1); Pstore(:,1)[P(1,1),P(2,2),P(3,3),P(4,4)]; Thetazeros(4,N1); % 參數(shù)估計(jì)值存儲(chǔ) Theta(:,1)[1;1;1;1]; % 故意給偏的參數(shù)初值 K[1;1;1;1]; % 增益矩陣初值占位 for i3:N2 h[-z(i-1);-z(i-2);M(i-1);M(i-2)]; % 回歸向量 KP*h*inv(h*P*h1); % 增益矩陣 Theta(:,i-1)Theta(:,i-2)K*(z(i)-h*Theta(:,i-2)); % 參數(shù)更新 P(eye(4)-K*h)*P; % 協(xié)方差更新 Pstore(:,i-1)[P(1,1),P(2,2),P(3,3),P(4,4)]; % 記錄方差對角元 end邏輯說明循環(huán)里三步依次是算 K、更新 Theta、更新 P順序是固定的——必須先用舊的 P 算 K再用 K 修正 Theta最后才更新 P。換順序會(huì)導(dǎo)致 P 和 K 不匹配參數(shù)收斂軌跡會(huì)亂掉。存儲(chǔ)數(shù)組 Pstore 和 Theta 的列數(shù)都是 N1從 1 到 501但循環(huán)體里 Theta(:,i-1) 的下標(biāo)范圍是 2 到 501所以初值放在第 1 列最后一個(gè)估計(jì)值放在第 501 列。畫圖時(shí) i1:N1 正好覆蓋全部估計(jì)點(diǎn)。參數(shù)說明P 的初值取 100eye(4)表示初始協(xié)方差為 100即參數(shù)不確定度很高。這個(gè)值如果取太小比如 1參數(shù)收斂會(huì)明顯變慢取太大比如 10^6則前幾步的參數(shù)估計(jì)會(huì)劇烈震蕩。Theta 初值取 [1;1;1;1] 離真值 [-1.8;0.9;0.75;0.5] 有不小距離這樣能直觀看到收斂過程。回歸向量 h 里的符號(hào)很關(guān)鍵因?yàn)椴罘址匠淌?z(k2)1.8z(k1)-0.9*z(k)...移到等式左邊后z 項(xiàng)的系數(shù)分別變成 -1.8 和 0.9所以 h 里對應(yīng)位置是 -z(i-1) 和 -z(i-2)這樣辨識(shí)出來的 θ(1) 和 θ(2) 會(huì)收斂到 -1.8 和 0.9。4. 漸消記憶最小二乘遺忘因子 0.998 背后的原理和改法4.1 為什么標(biāo)準(zhǔn) RLS 會(huì)有穩(wěn)態(tài)偏差標(biāo)準(zhǔn) RLS 的推導(dǎo)里P 矩陣會(huì)隨著迭代次數(shù)增加而不斷縮小。這意味著后期增益矩陣 K 的模越來越小參數(shù)的修正步長也隨之變小。理論上如果系統(tǒng)是線性定常的、噪聲是零均值的參數(shù)會(huì)收斂到真值附近但實(shí)際仿真中你往往會(huì)發(fā)現(xiàn)估計(jì)值和真值之間始終存在一個(gè)小的偏差而且這個(gè)偏差不隨迭代次數(shù)增加而消失。原因在于RLS 的代價(jià)函數(shù)對所有歷史數(shù)據(jù)等權(quán)重看待越老的數(shù)據(jù)積累越多在總代價(jià)里占比越大新數(shù)據(jù)對參數(shù)修正的貢獻(xiàn)相對被稀釋。如果系統(tǒng)本身有輕微時(shí)變或者噪聲不是純白噪聲比如有色噪聲這個(gè)偏差會(huì)更明顯。把系統(tǒng)輸出和參數(shù)估計(jì)代入程序畫出的圖上你會(huì)發(fā)現(xiàn)四條參數(shù)曲線都收斂了但每條曲線和真值之間都有一條肉眼可見的“縫”。這個(gè)縫的大小取決于噪聲水平、P 初值和數(shù)據(jù)長度。P 初值取 100、數(shù)據(jù)取 500 點(diǎn)時(shí)偏差通常較小但不為零——這正是程序注釋里寫的“辨識(shí)結(jié)果會(huì)與其精確值始終存在著偏差”。4.2 改成漸消記憶法的關(guān)鍵兩個(gè)語句而已漸消記憶法的思路是給歷史數(shù)據(jù)加一個(gè)指數(shù)衰減的權(quán)重越老的數(shù)據(jù)權(quán)重越小這樣新數(shù)據(jù)能在代價(jià)函數(shù)中占據(jù)更大比例從而持續(xù)修正參數(shù)。實(shí)現(xiàn)上只需要改兩處% 原語句 KP*h*inv(h*P*h1); P(eye(4)-K*h)*P; % 改為漸消記憶法遺忘因子 λ0.998 KP*h*inv(h*P*h0.998); P(eye(4)-K*h)*P/0.998;邏輯說明分母上的 1 改成 0.998對應(yīng)遺忘因子 λ。P 更新時(shí)除以 0.998即乘以 1/0.998 ≈ 1.002效果是讓 P 矩陣不會(huì)像標(biāo)準(zhǔn) RLS 那樣快速收縮到接近零。這樣每步結(jié)束后 P 維持在一個(gè)“有彈性”的量級(jí)下一個(gè)數(shù)據(jù)點(diǎn)到來時(shí) K 仍然有足夠大的模參數(shù)修正步長不會(huì)衰減到微乎其微。綜合效果就是估計(jì)值能持續(xù)跟隨真值穩(wěn)態(tài)偏差明顯縮小。參數(shù)說明遺忘因子 λ 的取值范圍通常是 0.951.0。λ 越接近 1記憶越長參數(shù)越穩(wěn)定但跟蹤越慢λ 越小比如 0.95記憶越短跟蹤時(shí)變參數(shù)的能力越強(qiáng)但穩(wěn)態(tài)波動(dòng)也越大。0.998 是一個(gè)比較保守的值適合參數(shù)基本不變但希望消除穩(wěn)態(tài)偏差的場景。如果你的系統(tǒng)參數(shù)是緩變的可以考慮取 0.980.99如果參數(shù)變化較快可能需要 0.950.97但這時(shí)要考慮噪聲放大的副作用。注意在這里噪聲方差假設(shè)從 1 變成了 λ0.998所以分母的含義不只是噪聲方差而是噪聲方差與遺忘因子的組合——這個(gè)細(xì)節(jié)很多人會(huì)忽略。4.3 兩種算法的結(jié)果對比怎么看對比標(biāo)準(zhǔn) RLS 和漸消記憶法的參數(shù)估計(jì)曲線最直觀的差異在穩(wěn)態(tài)段。標(biāo)準(zhǔn) RLS 的曲線雖然收斂但每條曲線和真值之間始終有偏差漸消記憶法的曲線幾乎貼近真值線。P 矩陣的對角元變化也解釋了這一點(diǎn)標(biāo)準(zhǔn) RLS 的 P 對角元持續(xù)衰減到接近零漸消記憶法的 P 對角元穩(wěn)定在一個(gè)較小的非零值附近。這個(gè)差異在你需要高精度參數(shù)值時(shí)非常關(guān)鍵——比如你辨識(shí)出來的參數(shù)接下來要用于設(shè)計(jì)控制器偏差可能導(dǎo)致閉環(huán)極點(diǎn)偏出期望位置。5. 避坑指南從初值設(shè)定到噪聲匹配四個(gè)最容易翻車的地方5.1 現(xiàn)象P 初值給太小導(dǎo)致收斂極慢現(xiàn)象把 P100eye(4) 改成 Peye(4) 或 P0.1eye(4)參數(shù)估計(jì)曲線需要很長時(shí)間才能逼近真值500 個(gè)數(shù)據(jù)點(diǎn)可能不夠用。原因P 初值代表對參數(shù)初始不確定度的估計(jì)給小了等于默認(rèn)“初始參數(shù)已經(jīng)很準(zhǔn)”增益矩陣 K 的初始步長就小參數(shù)只能慢慢爬過去。解決保持相對較大的初值如 100 倍單位陣或者根據(jù)參數(shù)可能的變化范圍手動(dòng)估計(jì)一個(gè)合適的數(shù)量級(jí)。我一般的做法是先跑一次標(biāo)準(zhǔn) RLS 看 P 對角元收斂到多少再反推一個(gè)合理的初值一般取收斂值的 10100 倍。5.2 現(xiàn)象遺忘因子太小導(dǎo)致參數(shù)估計(jì)波動(dòng)大現(xiàn)象把遺忘因子從 0.998 改成 0.95參數(shù)曲線雖然在真值附近波動(dòng)但波動(dòng)的幅度明顯變大看起來“毛刺”很多。原因遺忘因子越小有效數(shù)據(jù)窗口越短等效于用更少的數(shù)據(jù)做估計(jì)方差自然增大。解決先試 0.990.998 這個(gè)區(qū)間確認(rèn)穩(wěn)態(tài)偏差確實(shí)被消除后再嘗試適當(dāng)減小。如果偏差還在先檢查是不是噪聲方差沒匹配而不是一味調(diào)小遺忘因子。我自己的習(xí)慣是只在對參數(shù)有明確時(shí)變需求時(shí)才把遺忘因子調(diào)到 0.98 以下否則盡量保持在 0.99 附近。5.3 現(xiàn)象回歸向量的符號(hào)寫反辨識(shí)出的參數(shù)符號(hào)全錯(cuò)現(xiàn)象把 h[-z(i-1);-z(i-2);M(i-1);M(i-2)] 寫成 h[z(i-1);z(i-2);M(i-1);M(i-2)]結(jié)果 θ(1) 收斂到 1.8 而不是 -1.8θ(2) 收斂到 -0.9 而不是 0.9。原因差分方程移項(xiàng)后系數(shù)符號(hào)變了但回歸向量沒有同步調(diào)整。解決先把差分方程寫成標(biāo)準(zhǔn)回歸形式 y(k)h(k)*θv(k)逐項(xiàng)確認(rèn)符號(hào)再寫代碼。這一步最基礎(chǔ)但最容易錯(cuò)而且錯(cuò)的隱蔽——最終參數(shù)的絕對值是對的正負(fù)號(hào)卻是反的。5.4 現(xiàn)象M 序列初值全零導(dǎo)致激勵(lì)信號(hào)為零辨識(shí)完全失敗現(xiàn)象寄存器初值改成 x[0 0 0 0 0 0 0 0 0]程序運(yùn)行后參數(shù)曲線完全不動(dòng)或者發(fā)散到異常值。原因全零初值的 LFSR 永遠(yuǎn)輸出 0系統(tǒng)完全沒有被激勵(lì)觀測數(shù)據(jù)里不含任何參數(shù)信息。解決確保寄存器初值至少有一個(gè) 1最省事的就是保持程序原來的初值 [0 1 0 1 1 0 1 1 1]。如果你自己重新生成 M 序列記得在代碼里加一個(gè)判斷如果 sum(x)0就報(bào)錯(cuò)提示重新初始化。5.5 現(xiàn)象噪聲方差與 RLS 公式中的常數(shù)不匹配估計(jì)值偏離真值現(xiàn)象把噪聲改成 v2*randn(1,N)但 RLS 公式里的常數(shù)保持 1最后參數(shù)估計(jì)偏離真值且曲線波動(dòng)變大。原因公式里的常數(shù)同時(shí)承擔(dān)了噪聲方差設(shè)定和遺忘因子兩個(gè)角色噪聲實(shí)際方差和假設(shè)值不一致時(shí)增益矩陣的步長不是最優(yōu)的。解決如果你改了噪聲方差記得同步修改inv(h*P*h常數(shù))里的常數(shù)改成噪聲方差的數(shù)值。這是最容易被忽略的一個(gè)參數(shù)因?yàn)樗卦诰仃嚽竽娴谋磉_(dá)式里不是單獨(dú)一句話。6. 把辨識(shí)結(jié)果落地殘差校驗(yàn) 利用帶遺忘因子的 RLS 做在線跟蹤當(dāng)你把漸消記憶法的程序跑通、參數(shù)曲線也收斂到真值附近之后還有一個(gè)值得做的驗(yàn)證步驟檢查模型輸出和實(shí)際觀測之間的殘差。做法很簡單把辨識(shí)出的參數(shù)代回系統(tǒng)方程用同一段輸入重新仿真一遍然后算殘差% 辨識(shí)完成后用估計(jì)參數(shù)計(jì)算模型輸出 z_est zeros(N2,1); z_est(1) -1; z_est(2) 0; for i3:N2 z_est(i) -Theta(1,end)*z_est(i-1) - Theta(2,end)*z_est(i-2) ... Theta(3,end)*M(i-1) Theta(4,end)*M(i-2); end residual z(3:N2) - z_est(3:N2); % 殘差序列 figure(2) subplot(2,1,1) plot(1:N, residual) title(模型殘差序列) subplot(2,1,2) histogram(residual, 30) title(殘差直方圖)邏輯說明殘差等于實(shí)際輸出減模型輸出。如果辨識(shí)效果好殘差應(yīng)該接近白噪聲——均值接近零、直方圖呈鐘形。如果殘差里還有明顯的周期性或趨勢說明模型結(jié)構(gòu)沒對比如階次不夠高或者輸入信號(hào)沒有充分激勵(lì)出系統(tǒng)的所有模態(tài)。參數(shù)說明注意這里模型輸出用 Theta(1,end)、Theta(2,end) 取的是最后一列的估計(jì)值。Theta(1,end) 是辨識(shí)出的 -a1因?yàn)榛貧w向量里 h(1)-z(i-1)所以計(jì)算模型輸出時(shí)要寫成-Theta(1,end)*z_est(i-1)。這個(gè)符號(hào)細(xì)節(jié)和前面回歸向量里的符號(hào)是呼應(yīng)的別在這里搞反。一個(gè)我自己常用的習(xí)慣是在線辨識(shí)場景里不只要畫參數(shù)曲線還要同時(shí)畫 P 矩陣的跡或?qū)窃鸵驗(yàn)樗苤庇^反映“當(dāng)前參數(shù)可信度”。如果跡在迭代后期還在持續(xù)下降說明系統(tǒng)還沒進(jìn)入穩(wěn)態(tài)如果跡已經(jīng)穩(wěn)住但有波動(dòng)說明遺忘因子在起作用。從那以后我每次跑辨識(shí)實(shí)驗(yàn)都會(huì)強(qiáng)制走一遍“觀察 P 跡 → 檢查殘差 → 確認(rèn)符號(hào)”這三步再把結(jié)果用于后續(xù)設(shè)計(jì)。這一步看著多花幾分鐘實(shí)際能擋掉大量后期的低級(jí)返工。希望這份程序和拆解思路能幫到你尤其是那些剛開始接觸遞推辨識(shí)、想盡快把理論落到代碼上的朋友——按上面的步驟改參數(shù)、看曲線、檢查殘差你很快就能摸清這套算法的手感。本文還有配套的精品資源點(diǎn)擊獲取