递推最小二乘辨识.docx

实验4递推最小二乘法的实现

实验报告

哈尔滨工业大学

航天学院控制科学与工程系

专业：自动化

班级：1040101

姓名：

日期：2013年10月23日
1．实验题目：递推最小二乘法的实现
2．实验目的：
熟悉并掌握递推最小二乘法的算法原理。
3．实验主要原理
给定系统

（1）
其中，为待辨识的未知参数，是不相关随机序列。为系统的输出，为系统的输入。分别测出个输出、输入值，则可写出个方程，具体写成矩阵形式，有
（2）
设
，
则式（2）可写为
（3）
式中：y为N维输出向量；为N维噪声向量；为维参数向量；为测量矩阵。为了尽量减小噪声对估值的影响，应取，即方程数目大于未知数数目。
的最小二乘估计为
（4）
为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。设已获得的观测数据长度为N，将式（3）中的、和分别用来代替，即
（5）
用表示的最小二乘估计，则
（6）
令，则
（7）
如果再获得一组新的观测值和，则又增加一个方程
（8）
式中

将式（5）和式（8）合并，并写成分块矩阵形式，可得
（9）
于是，类似地可得到新的参数估值

（10）
式中

（11）
应用矩阵求逆引理，从求得与的递推关系式出发，经过一系列的推导，最终可求得递推最小二乘法辨识公式：
（12）
（13）
（14）
为了进行递推计算，需要给出和的初值和。
推荐取值方法为：假定，c是充分大的常数，I为单位矩阵，则经过若干次递推之后能得到较好的参数估计。

4．实验对象或参数
给定系统
（15）
即。假设实际系统的参数为，，，，，但是不已知，即不可测。取的零均值白噪声。输入信号取为
（16）
要求编制MATLAB程序，运用递推最小二乘法对这一系统的参数进行在线辨识，并将辨识结果与实际参数进行对比。

5．程序框图

6．程序代码
functionshiyan4
Y(1)=0;Y(2)=0;
n=2;theta=[2;1.3;0.4;0.88;2.2];
c=1000;
U(1)=1.5*sin(0.2*1);U(2)=1.5*sin(0.2*2);U(3)=1.5*sin(0.2*3);
Pn0=c*c*eye(5);
Tn0=zeros(5,1);
num=30;
data=zeros(num,5);
%生成白噪声
M=2147483647;
a=65539;
b=10000;
X=[];
R=[];
X(1)=35;
nn=12*(num+5);
forj=1:1:nn
x0=(a*X(j)+b);
X(j+1)=mod(x0,M);
R(j)=X(j+1)/M;
end
v=[];
forjj=0:1:num+4
jj1=jj*12+1;jj2=jj*12+11;x2=R(jj1);
forii=jj1:1:jj2
x2=x2+R(ii+1);
end
v(jj+1)=x2-6;
V(jj+1)=v(jj+1)/30;
end
%递推函数
while(num~=0)
HL(1)=-Y(n);HL(2)=-Y(n-1);
HL(3)=U(n+1);HL(4)=U(n);HL(5)=U(n-1);
Y(n+1)=HL*theta+V(n+1);
K=Pn0*HL'*(inv(1+HL*Pn0*HL'));
Pn=Pn0-K*HL*Pn0;
Tn=Tn0+K*(Y(n+1)-HL*Tn0);
data(n-1,:)=Tn;
Pn0=Pn;
Tn0=Tn;
n=n+1;
KK=[Tn(2)-16Tn(1)-4];
U(n+1)=KK*[Y(n);Y(n-1)]+[Tn(3)Tn(4)Tn(5)]*[1.5*sin(0.2*(n+1));1.5*sin(0.2*n);1.5*sin(0.2*(n-1))];
num=num-1;
end
data
%画图
a1=data(:,1);
a2=data(:,2);
b0=data(:,3);
b1=data(:,4);
b2=data(:,5);
k=0:1:length(a1)-1;
subplot(2,1,1);
stairs(k,a1,'r-');
xlabel('timek');
ylabel('a1');
subplot(2,1,2);
stairs(k,a2,'r-');
xlabel('timek');
ylabel('a2');
subplot(3,1,1);
stairs(k,b0,'r-');
xlabel('timek');
ylabel('b0');
subplot(3,1,2);
stairs(k,b1,'r-');
xlabel('timek');
ylabel('b1');
subplot(3,1,3);
stairs(