guangyizuixiaofangcha.m

这是黑龙江大学自动控制专业杨易鹏同学用matlab编写的广义最小二乘算法!

%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');