最小二乘法+相關(guān)分析法

關(guān)于最小二乘法+相關(guān)分析法第1頁，課件共26頁，創(chuàng)作于2023年2月最小二乘法1.問題描述考慮如下線性系統(tǒng)：

其中，u(k)為系統(tǒng)激勵信號，y(k)為系統(tǒng)輸出，e(k)為模型噪聲。第2頁，課件共26頁，創(chuàng)作于2023年2月系統(tǒng)模型第3頁，課件共26頁，創(chuàng)作于2023年2月

其中G(z-1)是系統(tǒng)函數(shù)模型，N(z-1)為有色噪聲系統(tǒng)模型，e(k)為白噪聲v(k)經(jīng)過系統(tǒng)函數(shù)為N(z-1)的系統(tǒng)后的輸出。通常第4頁，課件共26頁，創(chuàng)作于2023年2月則系統(tǒng)可表示為

設(shè)樣本和參數(shù)集為h(k)為可觀測的量，差分方程可寫為最小二乘形式第5頁，課件共26頁，創(chuàng)作于2023年2月

如何系統(tǒng)噪聲e(k)存在的情況下從該方程中正確的解出，即是系統(tǒng)辨識的任務(wù)。為了求出，我們面臨三大問題：一是輸入信號的選擇，二是判決準則的選取，三是辨識算法的選擇。選擇輸入

為了準確辨識系統(tǒng)參數(shù)，我們對輸入信號有兩大要求，一是信號要能持續(xù)的激勵系統(tǒng)所有狀態(tài)，二是信號頻帶能覆蓋系統(tǒng)的頻帶寬度。除此之外還要求信號有可重復(fù)性，不能是不可重復(fù)的隨機噪聲，因此我們通常選擇M序列或逆M序列作為輸入。第6頁，課件共26頁，創(chuàng)作于2023年2月準則函數(shù)

本次探討最小二乘類辨識方法，在此選取準則函數(shù)使準則函數(shù)的估計值記做，稱作參數(shù)的最小二乘估計值。令k=1,2,3,……L，可構(gòu)成線性方程組：第7頁，課件共26頁，創(chuàng)作于2023年2月式中：準則函數(shù)相應(yīng)變?yōu)闃O小化,求得參數(shù)的估計值，將使模型更好的預(yù)報系統(tǒng)的輸出。第8頁，課件共26頁，創(chuàng)作于2023年2月最小二乘法設(shè)使得，則有展開上式，并根據(jù)以下兩個向量微分公式：得正則方程：第9頁，課件共26頁，創(chuàng)作于2023年2月當為正則陣時，有且有所以唯一使得，這種通過極小化式計算的方法稱作最小二乘法。

第10頁，課件共26頁，創(chuàng)作于2023年2月遞推最小二乘法第11頁，課件共26頁，創(chuàng)作于2023年2月第12頁，課件共26頁，創(chuàng)作于2023年2月遞推最小二乘法（RLS）步驟如下第13頁，課件共26頁，創(chuàng)作于2023年2月舉例設(shè)某物理量Y與X滿足關(guān)系式Y(jié)=aX2+bX+c，實驗獲得一批數(shù)據(jù)如下表，試辨識模型參數(shù)a,b和c。X 1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10Y 9.6 4.1 1.3 0.4 0.05 0.1 0.7 1.8 3.8 9.0分析：Y=aX2+bX+c+=[X(k)*X(k)X(k)1]'*[abc]+=+第14頁，課件共26頁，創(chuàng)作于2023年2月程序X=[1.01 2.033.024.0156.027.038.049.0310];Y=[9.64.11.30.40.050.1 0.7 1.8 3.8 9.0];%實驗輸入數(shù)據(jù)、實驗輸出數(shù)據(jù)symsabc%定義待辨識參數(shù)theta=[a;b;c];%theta包含待辨識參數(shù)a，b，ctheta1=zeros(3,1);%對象參數(shù)初始化P=10^6*eye(3);%構(gòu)造初始P陣fork=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);1]K=P*phi/(1+phi'*P*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%輸出參數(shù)估計值plot([1:10],theta2)%輸出參數(shù)逐步遞推估計的軌跡圖像xlabel('k');%設(shè)置橫坐標為步長kylabel('參數(shù)估計a，b，c');%縱坐標為估計參數(shù)a，b，clegend('a','b','c');%標示相應(yīng)曲線對應(yīng)的參數(shù)axis([110-1020]);%設(shè)置坐標軸范圍第15頁，課件共26頁，創(chuàng)作于2023年2月實驗結(jié)果theta1=0.4575-5.073413.3711第16頁，課件共26頁，創(chuàng)作于2023年2月相關(guān)分析法原理：一個單入單出線性定常系統(tǒng)的動態(tài)特性可用它的脈沖響應(yīng)函數(shù)g(σ)來描述。第17頁，課件共26頁，創(chuàng)作于2023年2月令輸入為白噪聲，這時的自相關(guān)函數(shù)為：則根據(jù)維納—霍夫積分方程可得：只要記錄x(t)、y(t)的值，并計算它們的互相關(guān)函數(shù)，即可求得脈沖響應(yīng)函數(shù)。第18頁，課件共26頁，創(chuàng)作于2023年2月實例分析函數(shù)模型

其中K=120,T1=8.2Sec,T2=6.1Sec

輸入信號采用6階M序列，周期為63bit.幅度為1，采樣周期為1s,利用相關(guān)分析法計算出系統(tǒng)的脈沖響應(yīng)值,并與系統(tǒng)的脈沖響應(yīng)理論值.實驗步驟1：得到仿真的輸入輸出數(shù)據(jù)M序列2：利用相關(guān)分析法，求得系統(tǒng)的脈沖響應(yīng)函數(shù)。3：利用相關(guān)分析法求得的系統(tǒng)脈沖響應(yīng)函數(shù)與實際脈沖響應(yīng)函數(shù)。4：對系統(tǒng)的輸出加入噪聲，畫圖比較利用相關(guān)分析法求得的系統(tǒng)脈沖響應(yīng)函數(shù)與實際脈沖響應(yīng)函數(shù)。第19頁，課件共26頁，創(chuàng)作于2023年2月程序%創(chuàng)建M序列Np=63;%循環(huán)周期delta_T=1;%時鐘節(jié)拍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;forn=1:Nptemp=xor(M(6),M(5));if(temp==0)M_XuLie(n)=a;elseM_XuLie(n)=-a;endM(6)=M(5);M(5)=M(4);M(4)=M(3);M(3)=M(2);M(2)=M(1);M(1)=temp;end%stairs(M_XuLie)%生成M序列完畢第20頁，課件共26頁，創(chuàng)作于2023年2月r=3;%周期數(shù)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)=0fork=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/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');holdon;plot(x,'k');plot(y,'b');legend('u(k)','x(k)','y(k)');第21頁，課件共26頁，創(chuàng)作于2023年2月%第二步，將白噪聲添加入輸出信號%產(chǎn)生白噪聲信號vfangcha=0.5;%隨意指定的方差v=fangcha*randn(1,63*4);%信號疊加,輸出實際信號z(k)z=y+v;%figure(2);%plot(v,'m');%holdon;%plot(z,'y');figure(2);%打印無白噪聲污染信號plot(y,'b');holdon;%打印白噪聲信號plot(v,'m');%打印白噪聲污染后的信號plot(z,'g');legend('y(k)','v(k)','z(k)');%計算Rmz(k)fork=1:NpRmz(k)=0;%初始化為0fori=(Np+1):((r+1)*Np)Rmz(k)=Rmz(k)+u(i-k)*z(i);endRmz(k)=Rmz(k)/(r*Np); end第22頁，課件共26頁，創(chuàng)作于2023年2月%計算cc=-Rmz(Np-1);%計算脈沖響應(yīng)估計值g1g1=Np*(Rmz+c)/((Np+1)*a^2*delta_T);%計算理論脈沖g0fork=1:Npg0(k)=K/(T1-T2)*(exp(-k*delta_T/T1)-exp(-k*delta_T/T2));end%計算脈沖響應(yīng)估計誤差delta_gdelta_g=sqrt