%系統辨識實驗--批量廣義最小二乘法
%實驗數據預處理
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（3）;
plot（Data（:1）‘--blue‘）;
hold?on;
plot（Data（:2）‘--red‘）;
hold?on;
title（‘批量增廣最小二乘法‘）;

%批量最小二乘法
D=3;????%純滯后
N=M-20;??%一次批處理采用的數據量
%第一次假設系統結構
%y（k）+a（1）y（k-1）+a（2）y（k-2）=b（1）x（k-1-D）+b（2）x（k-2-D）+w（k）
y_data=Data（:1）;
x_data=Data（:2）;
fai=zeros（N4）;
for?i=3:N+3
????y_f（i-2:）=y_data（i+D）;
????fai（i-2:）=[-y_data（i-1+D）-y_data（i-2+D）x_data（i-1）x_data（i-2）];
end
theta1=inv（fai‘*fai）*fai‘*y_f;


while（1）
????%利用前面得到的theta估計e（k）
????e_data=zeros（1M）;
????for?i=1:D+1
????????y_e（i）=y_data（i）;
????end;
????for?i=D+2:M
????????y_e（i）=-theta1（1）*y_e（i-1）-theta1（2）*y_e（i-2）+theta1（3）*x_data（i-D）+theta1（4）*x_data（i-1-D）;
????????e_data（i）=y_data（i）-y_e（i）;
????end
????%估計c（1）c（2）值
????%假設e的模型：e（k）=c（1）*e（k-1）+c（2）*e（k-2）+w（k）
????fai=zeros（N2）;
????for?i=3:N+3
????????e_f（i-2:）=e_data（i）;
????????fai（i-2:）=[-e_data（i-1）-e_data（i-2）];
????end
????theta_c=inv（fai‘*fai）*fai‘*e_f;
????disp（theta_c）;
????%用c（1）c（2）對xy進行濾波得到xbyb
????for?i=1:2
????????xb_data（i）=x_data（i）;
????????yb_data（i）=y_data（i）;
????end
????for?i=3:M
????????xb_data（i）=theta_c（1）*x_data（i-1）+theta_c（2）*x_data（i-2）;
????????yb_data（i）=theta_c（1）*y_data（i-1）+theta_c（2）*y_data（i-2）;
????end

????%用xbyb進行最小二乘估計
????fai=zeros（N4）;
????y_f=zeros（N1）;
????for?i=3:N+3
????????y_f（i-2:）=yb_data（i+D）;
????????fai（i-2:）=[-yb_data（i-1+D）-yb_data（i-2+D）xb_data（i-1）xb_data（i-2）];
????end
????theta=inv（fai‘*fai）*fai‘*y_f;
????
????theta_diff=（theta-theta1）./theta1;
????for?i=1:2
????????if?theta_diff（i）<0
????????????theta_diff（i）=-theta_diff（i）;
????????end
????end
????
???if?（max（abs（theta_diff））<0.0001）
????????break;
????????else
????????theta1=theta;
????end?????
end

disp（theta）;
%增廣最小二乘法與實際數據的比較
for?i=1:D+2
????y_e（i）=y_data（i）;
end;
%?ne=0.0025*randn（M1）;
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）+theta_c（1）*ne（i-11）+theta_c（2）*ne（i-21）;
????????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（3）;
hold?on;
plot（y_e‘--g‘）;
figure（1）;
plot（y_data-y_e‘‘g‘）;
title（‘誤差曲線‘）;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????5015??2009-01-09?21:06??RLS2.M

?????文件???????6178??2009-01-10?00:06??RLS_JXJY.m

?????文件???????2797??2009-01-10?14:41??GLSE.m

?????文件???????3839??2009-01-06?16:43??LS.m

?????文件???????1121??2009-04-23?17:01??nerual_xtbs.m

?????文件???????5778??2009-01-10?11:20??RELS.m

?????文件???????1774??2009-01-10?12:41??RELS_JXJY.m

-----------?---------??----------?-----??----

????????????????26502????????????????????7