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