%渐消记忆最小二乘法
%实验数据预处理 xjirong
clear
filename='SI_Data.xls';
Data_org=xlsread(filename);
Max_y=max(Data_org(:,1));
Min_y=min(Data_org(:,1));
Max_x=max(Data_org(:,2));
Min_x=min(Data_org(:,2));
%对输入数据进行重新采样，减少点数
%每Q个点钟重新取出一个点来
Q=8;
M=floor(length(Data_org)/Q);
for i=1:M
Data_new(i,:)=Data_org(1+(i-1)*Q,:);
end
Data(:,1)=(Data_new(:,1)-Min_y)./(Max_y-Min_y);
Data(:,2)=(Data_new(:,2)-Min_x)./(Max_x-Min_x);
figure(1);
plot(Data(:,1),'--black');
hold on;
plot(Data(:,2),'--red');
hold on;
title('渐消记忆最小二乘法');
N=M; %递推出处理采用的数据量
y_data=Data(:,1);
x_data=Data(:,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
D=3; %纯滞后
delta=0.995; %遗忘因子
%假设系统结构
%y(k)=a(1)y(k-1)+a(2)y(k-2)+b(1)x(k-1-D)+b(2)x(k-2-D)+e(k)
theta_s=zeros(4,N);%[a1;a2;b1;b2]
P=100000*eye(4);
for k=D+2:N
X=[y_data(k-1);y_data(k-2);x_data(k-1-D);x_data(k-1-D)];
K=P*X/(delta+X'*P*X);
theta_s(:,k)=theta_s(:,k-1)+K*[y_data(k)-X'*theta_s(:,k-1)];
P=(P-K*X'*P)/delta;
end
%求最后的系数的平均值
sum=zeros(4,1);
for i=1:20
sum=sum+theta_s(:,N-i);
end
theta=sum/20;
% disp(theta);
%最小二乘法与实际数据的比较
for i=1:D+1
y_e(i)=y_data(i);
end;
y_e(D+2)=theta(1)*y_e(D+1)+theta(2)*y_e(D)+theta(3)*x_data(1);
for i=D+3:M
y_e(i)=theta(1)*y_e(i-1)+theta(2)*y_e(i-2)+theta(3)*x_data(i-1-D)+theta(4)*x_data(i-2-D);
end
%绘制预测曲线
figure(1);
hold on;
plot(y_e,':g');
disp('D=3,delta=0.995,J');
disp((y_data-y_e')'*(y_data-y_e'));
disp('二阶系数');
disp(theta);
%绘制误差曲线
figure(2);
hold on;
plot(y_data-y_e',':g');
%绘制系数收敛曲线
figure(3);
subplot(2,2,1);
plot(theta_s(1,:),'g');
title('a(1)系数收敛曲线');
hold on;
% figure(4);
subplot(2,2,2);
plot(theta_s(2,:),'g');
title('a(2)系数收敛曲线');
hold on;
% figure(5);
subplot(2,2,3);
plot(theta_s(3,:),'g');
title('b(1)系数收敛曲线');
hold on;
% figure(6);
subplot(2,2,4);
plot(theta_s(4,:),'g');
title('b(2)系数收敛曲线');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
D=6; %纯滞后
delta=0.995; %遗忘因子
theta_s=zeros(4,N);%[a1;a2;b1;b2]
P=100000*eye(4);
for k=D+2:N
X=[y_data(k-1);y_data(k-2);x_data(k-1-D);x_data(k-1-D)];
K=P*X/(delta+X'*P*X);
theta_s(:,k)=theta_s(:,k-1)+K*[y_data(k)-X'*theta_s(:,k-1)];
P=(P-K*X'*P)/delta;
end
%求最后的系数的平均值
sum=zeros(4,1);
for i=1:20
sum=sum+theta_s(:,N-i);
end
theta=sum/20;
% disp(theta);
%最小二乘法与实际数据的比较
for i=1:D+1
y_e(i)=y_data(i);
end;
y_e(D+2)=theta(1)*y_e(D+1)+theta(2)*y_e(D)+theta(3)*x_data(1);
for i=D+3:M
y_e(i)=theta(1)*y_e(i-1)+theta(2)*y_e(i-2)+theta(3)*x_data(i-1-D)+theta(4)*x_data(i-2-D);
end
%绘制预测曲线
figure(1);
hold on;
plot(y_e,':b');
disp('D=6,delta=0.995,J');
disp((y_data-y_e')'*(y_data-y_e'));
disp('二阶系数');
disp(theta);
%绘制误差曲线
figure(2);
hold on;
plot(y_data-y_e',':b');
%绘制系数收敛曲线
figure(3);
subplot(2,2,1);
plot(theta_s(1,:),'b');
% title('a(1)系数收敛曲线');
hold on;
% figure(4);
subplot(2,2,2);
plot(theta_s(2,:),'b');
% title('a(2)系数收敛曲线');
hold on;
% figure(5);
subplot(2,2,3);
plot(theta_s(3,:),'b');
% title('b(1)系数收敛曲线');
hold on;
% figure(6);
subplot(2,2,4);
plot(theta_s(4,:),'b');
% title('b(2)系数收敛曲线');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D=9; %纯滞后
delta=0.995; %遗忘因子
theta_s=zeros(4,N);%[a1;a2;b1;b2]
P=100000*eye(4);
for k=D+2:N
X=[y_data(k-1);y_data(k-2);x_data(k-1-D);x_data(k-1-D)];
K=P*X/(delta+X'*P*X);
theta_s(:,k)=theta_s(:,k-1)+K*[y_data(k)-X'*theta_s(:,k-1)];
P=(P-K*X'*P)/delta;
end
%求最后的系数的平均值
sum=zeros(4,1);
for i=1:20
sum=sum+theta_s(:,N-i);
end
theta=sum/20;
% disp(theta);
%最小二乘法与实际数据的比较
for i=1:D+1
y_e(i)=y_data(i);
end;
y_e(D+2)=theta(1)*y_e(D+1)+theta(2)*y_e(D)+theta(3)*x_data(1);
for i=D+3:M
y_e(i)=theta(1)*y_e(i-1)+theta(2)*y_e(i-2)+theta(3)*x_data(i-1-D)+theta(4)*x_data(i-2-D);
end
%绘制预测曲线
figure(1);
hold on;
plot(y_e,':r');
legend('实际输出','实际输入','D=3,delta=0.995','D=6,delta=0.995','D=9,delta=0.995',4);
disp('D=9,delta=0.994,J');
disp((y_data-y_e')'*(y_data-y_e'));
disp('二阶系数');
disp(theta);
%绘制误差曲线
figure(2);
hold on;
plot(y_data-y_e',':r');
legend('D=3,delta=0.995','D=6,delta=0.995','D=9,delta=0.995',1);
title('误差曲线');
%绘制系数收敛曲线
figure(3);
subplot(2,2,1);
plot(theta_s(1,:),'r');
legend('D=3,delta=0.995','D=6,delta=0.995','D=9,delta=0.995',4);
hold on;
% figure(4);
subplot(2,2,2);
plot(theta_s(2,:),'r');
hold on;
% figure(5);
subplot(2,2,3);
plot(theta_s(3,:),'r');
hold on;
% figure(6);
subplot(2,2,4);
plot(theta_s(4,:),'r');
hold on;
评论0