最小二乘法系统辨识matlab 程序

用matlab实现的最小二乘法用于系统辨识，从文件中读取数据，然后辨识。里面包括各种基于最小二乘法的算法，普通最小二乘，广义最小二乘。。。等等

%系统辨识实验--批量广义最小二乘法
%实验数据预处理
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