1、關于最小二乘法+相關分析法現在學習的是第一頁，共26頁最小二乘法1.問題描述考慮如下線性系統： 其中，u(k)為系統激勵信號，y(k)為系統輸出，e(k)為模型噪聲。 1111abnanbz kaz ka z k nbu kb u k ne k 現在學習的是第二頁，共26頁系統模型現在學習的是第三頁，共26頁 其中G(z-1)是系統函數模型，N(z-1)為有色噪聲系統模型，e(k)為白噪聲v(k)經過系統函數為N(z-1)的系統后的輸出。通常111111, B zD zG zN zA zC z11212112121aabbnnnnA za za za zB zb zb zb z 1121211

2、2121ccddnnnnC zc zc zc zD zd zd zb z 現在學習的是第四頁，共26頁則系統可表示為 設樣本和參數集為 1111BzDzzkukvkAzCz121 2( ) -( -1 ) , - ( -2), . - ( - ), ( -1 ), ( -2), ., ( - ) , ,., , , ,., TTnnhkzkzkzk n ukukuk na aa bbbh(k)為可觀測的量， 差分方程可寫為最小二乘形式 ( )Tz khke k現在學習的是第五頁，共26頁 如何系統噪聲e(k)存在的情況下從該方程中正確的解出，即是系統辨識的任務。 為了求出，我們面臨三大問題：一

3、是輸入信號的選擇，二是判決準則的選取，三是辨識算法的選擇 。選擇輸入 為了準確辨識系統參數，我們對輸入信號有兩大要求，一是信號要能持續的激勵系統所有狀態，二是信號頻帶能覆蓋系統的頻帶寬度。除此之外還要求信號有可重復性，不能是不可重復的隨機噪聲，因此我們通常選擇M序列或逆M序列作為輸入。 現在學習的是第六頁，共26頁準則函數 本次探討最小二乘類辨識方法，在此選取準則函數 使準則函數的估計值記做，稱作參數的最小二乘估計值。 令k=1,2,3,L，可構成線性方程組： 2211TkkJe kz khk TLLLzkHkek現在學習的是第七頁，共26頁式中：準則函數相應變為 極小化,求得參數的估計值，將

4、使模型更好的預報系統的輸出。 1122, 0101121211LLababLabzezezez Le LzznuunzznuunHz Lz Lnu Lu Ln 2211LLTTLLLLkkJe kz khkzHzH現在學習的是第八頁，共26頁最小二乘法 設 使 得，則有展開上式，并根據以下兩個向量微分公式：得正則方程： LS minJ 0LSTLLLLJzHzH2 TTTTa xaxx Axx AAx為對稱陣TTLLLSLLH HH z現在學習的是第九頁，共26頁當 為正則陣時，有且有 所以唯一 使得 ，這種通過極小化式計算的方法稱作最小二乘法。 TLLH H1TTLSLLLLH HH z 2

5、220LSTLLJH HLS minJ現在學習的是第十頁，共26頁遞推最小二乘法現在學習的是第十一頁，共26頁現在學習的是第十二頁，共26頁遞推最小二乘法（RLS）步驟如下 現在學習的是第十三頁，共26頁舉例 設某物理量Y與X 滿足關系式Y=aX2+bX+c，實驗獲得一批數據如下表，試辨識模型參數a,b和c。 X1.012.033.024.0156.027.038.049.0310 Y9.64.11.30.40.050.10.71.83.89.0 分析：Y=aX2+bX+c+ =X(k)*X(k) X(k) 1*a b c+ = +現在學習的是第十四頁，共26頁程序X=1.01 2.03 3.

6、02 4.01 5 6.02 7.03 8.04 9.03 10;Y=9.6 4.1 1.3 0.4 0.05 0.10.71.83.89.0;%實驗輸入數據、實驗輸出數據syms a b c % 定義待辨識參數theta=a;b;c; %theta包含待辨識參數a，b，ctheta1=zeros(3,1); %對象參數初始化P=106*eye(3); %構造初始P陣for k=1:10 %仿真步長范圍1到10phi=X(k)*X(k);X(k);1;%y=aX*X+bX+c=phi*theta%theta=a；b；c；phi=X(k)*X(k);X(k);1K=P*phi/(1+phi*P*

7、phi); %遞推最小二乘法K陣的遞推公式theta=theta1+K*(Y(k)-phi*theta1); %theta的遞推公式P=(eye(3)-K*phi)*P; %遞推最小二乘法P陣的遞推公式theta1=theta; %theta的最終估計向量theta2(:,k)=theta; %theta估計向量矩陣化，目的是為了%下面的plot仿真圖像輸出endtheta1 %輸出參數估計值plot(1:10,theta2) %輸出參數逐步遞推估計的軌跡圖像xlabel(k); %設置橫坐標為步長kylabel(參數估計a，b，c); %縱坐標為估計參數a，b，clegend(a,b,c);

8、 %標示相應曲線對應的參數axis(1 10 -10 20); %設置坐標軸范圍現在學習的是第十五頁，共26頁實驗結果 theta1 =0.4575 -5.0734 13.3711現在學習的是第十六頁，共26頁相關分析法 原理：一個單入單出線性定常系統的動態特性可用它的脈沖響應函數g()來描述。0 ( )( ) ()y tgx td則000()11lim( ) ()( )lim() () TTTTx ty t x tdtgx tx tdt dTT上式兩端同乘，進而取時間均值，有現在學習的是第十七頁，共26頁令輸入為白噪聲，這時的自相關函數為：則根據維納霍夫積分方程可得： 只要記錄x(t)、y(

9、t)的值，并計算它們的互相關函數，即可求得脈沖響應函數。dtttgRRxxxy0)()()()(kRxx)(1)(Rxykg現在學習的是第十八頁，共26頁實例分析函數模型 其中K=120,T1=8.2Sec,T2=6.1Sec 輸入信號采用6階M序列，周期為63bit. 幅度為1，采樣周期為1s, 利用相關分析法計算出系統的脈沖響應值, 并與系統的脈沖響應理論值. 實驗步驟 1 ：得到仿真的輸入輸出數據M序列 2：利用相關分析法，求得系統的脈沖響應函數。 3：利用相關分析法求得的系統脈沖響應函數與實際脈沖響應函數。 4：對系統的輸出加入噪聲，畫圖比較利用相關分析法求得的系統脈沖響應函數與實際脈

10、沖響應函數。現在學習的是第十九頁，共26頁程序%創建M序列Np=63;%循環周期delta_T = 1;%時鐘節拍a=1;%幅度%初始化M序列M(1)=1;M(2)=0;M(3)=0;M(4)=1;M(5)=1;M(6)=0;M_XuLie(Np) = 0;for n = 1 : Np temp = xor(M(6), M(5); if(temp = 0) M_XuLie(n) = a; else M_XuLie(n) = -a; end M(6) = M(5); M(5) = M(4); M(4) = M(3); M(3) = M(2); M(2) = M(1); M(1) = temp;e

11、nd%stairs(M_XuLie)%生成M序列完畢現在學習的是第二十頁，共26頁r = 3; % 周期數u=repmat(M_XuLie,1,r+1);%將M序列賦給輸入，作為輸入信號%第一步，從u(k)得到x(k),y(k)K = 120;T0 = 1;% 采樣時間T1 = 8.1;T2 = 6.2;K1=K/(T1*T2);%初始化X(k),Y(k)為0K2=1x(63)=0;y(63)=0for k = 2 : 63*4 %取得x(k)序列 x(k)=exp(-T0/T1)*x(k-1)+T1*K1*(1-exp(-T0/T1)*u(k-1)+T1*K1. *(T1*(exp(-T0/

12、T1)-1)+T0)*(u(k)-u(k-1)/T0 %取得y(k)序列 y(k)=exp(-T0/T2)*y(k-1)+T2*K2*(1-exp(-T0/T2)*x(k-1)+T2*K2. *(T2*(exp(-T0/T2)-1)+T0)*(x(k)-x(k-1)/T0end%獲取沒有白噪聲時候輸出完畢%作圖figure(1);plot(u,r);hold on;plot(x,k);plot(y,b);legend(u(k),x(k),y(k);現在學習的是第二十一頁，共26頁%第二步，將白噪聲添加入輸出信號%產生白噪聲信號vfangcha = 0.5;%隨意指定的方差v = fangcha

13、 * randn(1,63*4);%信號疊加,輸出實際信號z(k)z = y + v;%figure(2);%plot(v,m);%hold on;%plot(z,y);figure(2);%打印無白噪聲污染信號plot(y,b);hold on;%打印白噪聲信號plot(v,m);%打印白噪聲污染后的信號plot(z,g);legend(y(k),v(k),z(k);%計算Rmz(k)for k = 1 : Np Rmz(k)=0;%初始化為0 for i = (Np + 1) : (r+1)*Np) Rmz(k)=Rmz(k) + u(i-k)*z(i); end Rmz(k)=Rmz(k)/(r*Np);end現在學習的是第二十二頁，共26頁%計算cc=-Rmz(Np - 1);%計算脈沖響應估計值g1g1=Np*(Rmz+c)/(Np+1)*a2*delta_T);%計算理論脈沖g0for k = 1: Np g0(k)=K/(T1-T2)*(exp(-k*delta_T/T1)-exp(-k*delta_T/T2);end%計算脈沖響