www.pudn.com > guangyizuixiaofangcha.rar > guangyizuixiaofangcha.m

%y(t)=0.8y(t-1)+0.5u(t-2)+kesai(t) 
% J=E[(y(t+2)-yr(t))^2+0.1u(t)^2] 
%yr(t)为周期为200，幅度为20的方波，总步数t= 800 
clc 
clear 
Qe=0.81; 
randn('seed',1); 
M=800; 
kesai=sqrt(Qe)*randn(1,M+10); 
yr(1)=20;%*****参考输入******* 
% y(1)=0;y(2)=0;;%**********初值选取******** 
% u(1)=0;u(2)=0; 
y(1)=0.01;y(2)=0.01;u(1)=0.1;u(2)=0.1; 
 
sita(:,2)=[0.1 0.1]';%sita=[a1 b0]' 
sita(:,1)=[0.1 0.1]';%sita=[a1 b0]' 
P(:,5:6)=100000*eye(2); 
fai(:,1)=[0 0]';  %fai'=[-y(t-1) u(t-2)] 
fai(:,2)=[0 0]'; 
for t=1:M+10 
    if mod(t,100)==0 
        yr(t+1)=-yr(t); 
%         yr(t+2)=-yr(t); 
    else 
        yr(t+1)=yr(t); 
%         yr(t+2)=yr(t); 
    end 
end 
% Austrom 方法 
for t=3:M 
    y(t)=0.8*y(t-1)+0.5*u(t-2)+kesai(t); 
    P(:,2*t+1:2*t+2)=P(:,2*(t-1)+1:2*(t-1)+2)-P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2)*... 
                  (P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2))'/(1+fai(:,t-2)'*P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2)); 
    sita(:,t)=sita(:,t-1)+(P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2)*(y(t)-fai(:,t-2)'... 
                  *sita(:,t-1))/(1+fai(:,t-2)'*P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2))); 
%   f1'=sita(1,t); 
%   g0=sita(1,t)^2; 
% u(t)=(-sita(2,t)*sita(1,t)*u(t-1)+yr(t)-sita(1,t)*sita(1,t)*y(t))/(sita(2,t)+0.1/sita(2,t)); 
     u(t)=1/(0.1/sita(2,t)+sita(2,t))*(-sita(1,t)*sita(2,t)*u(t-1)+yr(t)-sita(1,t)^2*y(t)); 
     %fai(:,t-2)=[-y(t-1) u(t-2)]'; 
    fai(:,t-1)=[y(t) u(t-1)]';  
end 
% for i=3:M 
%     fai(:,i)=[y(i-1) u(i-2)]'; 
%     y(i)=0.8*y(i-1)+0.5*u(i-2)+e(i); 
%     sita(:,i)=sita(:,i-1)+(P*fai(:,i)*(y(i)-fai(:,i)'*sita(:,i-1)))/(1+fai(:,i)'*P*fai(:,i)); 
%     P=P-(P*fai(:,i)*fai(:,i)'*P)/(1+fai(:,i)'*P*fai(:,i)); 
%     f1(i)=sita(1,i); 
%     g1(i)=sita(1,i)*f1(i); 
%     u(i)=(-sita(2,i)*f1(i)*u(i-1)+yr(i)-g1(i)*y(i))/(sita(2,i)+0.1/sita(2,i)); 
% end 
% Bok-Jenkins 方法 
% for t=3:M 
%     y(t)=0.8*y(t-1)+0.5*u(t-2)+kesai(t); 
%     P(:,2*t+1:2*t+2)=P(:,2*(t-1)+1:2*(t-1)+2)-P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2)*... 
%                   (P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2))'/(1+fai(:,t-2)'*P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2)); 
%     sita(:,t)=sita(:,t-1)+(P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2)*(y(t)-fai(:,t-2)'... 
%                   *sita(:,t-1))/(1+fai(:,t-2)'*P(:,2*(t-1)+1:2*(t-1)+2)*fai(:,t-2))); 
% %   f1'=sita(1,t); 
% %   g0=sita(1,t)^2; 
%     u(t)=1/(0.1/sita(2,t)+0.5)*(yr(t)-sita(1,t)*sita(1,t)*y(t)-sita(1,t)*sita(2,t)*u(t-1)); 
%     %fai(:,t-2)=[-y(t-1) u(t-2)]'; 
%     fai(:,t-1)=[y(t) u(t-1)]';  
% end 
for t=1:M 
e(t)=y(t)-yr(t); 
end 
a=yr(800)-y(800) 
b=yr(730)-y(730) 
sita(:,800) 
for t=1:M 
    m(t)=0.8; 
    n(t)=0.5; 
    l1(t)=0; 
    l2(t)=0; 
    l3(t)=20; 
    l4(t)=-20; 
end 
figure 
t=1:200; 
subplot(2,2,1); 
plot(t,m(t),t,sita(1,t),'r'); 
subplot(2,2,2); 
plot(t,n(t),t,sita(2,t),'r'); 
figure 
t=1:800; 
subplot(2,2,1); 
plot(t,l1(t),t,u(t),'r'); 
subplot(2,2,2); 
plot(t,l2(t),t,y(t),'r',t,l3(t),t,l4(t)); 
figure 
t=1:800; 
subplot(2,2,1); 
plot(t,e(t),'r');