限定記憶最小二乘法

%限定記憶最小二乘的遞推算法%Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k)%========================================clearclc%==========400產(chǎn)生M序列作為輸入===============x=[010110111];%initialvaluen=403;%n為脈沖數(shù)目M=[];%存放M序列fori=1:ntemp=xor(x(4),x(9));M(i)=x(9);forj=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========產(chǎn)生均值為0，方差為1的高斯白噪聲=========v=randn(1,402);%==============產(chǎn)生觀測序列=================z=zeros(402,1);z(1)=-1;z(2)=0;fori=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i);end%遞推求解P_a=100*eye(4);%估計方差Theta_a=[3;3;3;3];Pstore=zeros(4,381);%存放方差中間過程Theta_Store=zeros(4,381);%參數(shù)的估計值，存放中間過程估值L=20;%記憶長度fori=3:L-1h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];K=P_a*h*inv(h'*P_a*h+1);Theta_a=Theta_a+K*(z(i)-h'*Theta_a);P_a=(eye(4)-K*h')*P_a;endfork=0:380hL=[-z(k+L-1);-z(k+L-2);M(k+L-1);M(k+L-2)];K_b=P_a*hL*inv(1+hL'*P_a*hL);Theta_b=Theta_a+K_b*(z(k+L)-hL'*Theta_a);P_b=(eye(4)-K_b*hL')*P_a;hk=[-z(k+L);-z(k+L-1);M(k+L);M(k+L-1);];K_a=P_b*hk*inv(1+hk'*P_b*hk);Theta_a=Theta_b-K_a*(z(k+L+1)-hk'*Theta_b);P_a=(eye(4)+K_a*hk')*P_b;Pstore(:,k+1)=[P_a(1,1),P_a(2,2),P_a(3,3),P_a(4,4)];Theta_Store(:,k+1)=Theta_a;end%========================輸出結(jié)果及作圖===========================disp('參數(shù)ala2blb2的估計值為：')Theta_ai=1:381;figure(1)plot(i,Theta_Store(1,:),i,Theta_Store(2,:),i,Theta_Store(3,:),i,Theta_Store(4,:))title('待估參數(shù)過渡過程')figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:))title('估計方差過渡過程')運行結(jié)果如下：參數(shù)a1a2b1b2的估計值為：Theta_a=-1.54260.72971.00690.3365彳MATLAB7.11.0(R2010b)FileEditDebugParallelDesktopWindowHelp:DtS|#■?町O|擲國薊I?I[P:\ProgramFiles\MATLAB\R2010b\binShortcutsHowtoAdd創(chuàng)What'sNewCurrentFolder ?x衛(wèi)；bin? ▼Ql￡|Name▼勺yiwangyinNi.m巧CurrentFolder ?x衛(wèi)；bin? ▼Ql￡|Name▼勺yiwangyinNi.m巧yasuochuli2.m巧yasuochulil.m勺yasuochuli.m*1xindianxinhaol.m*1xindianxinhao.m勺xiaozaochuli.mxiaoguoduibi.m巧xiandingjiyi.mmworker.bat巧tuxiang.m巧tubianxinhao.m勺tubiandian.m_|tubiandian.asv*1qiyidian.m_,qiyidian.asvDProductRootsHImw_mpiexec.bat一msvc_modules_installer.pm,,mexutils.pm_|mexsetup.pm鬲mexext.bat.,mex.pl1(強1mayDetailsWorkspace *口？'x肓i到甩■11國叨“”JHddJHJHddJHddJdJHdd■Value[0.0080;0.0(*[-0.0068;0.0[0.0035;-0.0三20<1x403doi—<4x4doubl<4x4doubl<4x381doi<4x381doi[-1.5426;0.7▼CommandHistory -? image(IJ:”…colormap(map):nbcol=size(map,1)'?…colorbar:；%—2011/6/3019:08—%；%—2011/6/3020:00—%i%—