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