系统辨识最小二乘参数估计matlab

1、精选优质文档-倾情为你奉上最小二乘参数估计摘要：最小二乘的一次性完成辨识算法（也称批处理算法），他的特点是直接利用已经获得的所有（一批）观测数据进行运算处理。这种算法在使用时，占用内存大，离线辨识，观测被辨识对象获得的新数据往往是逐次补充到观测数据集合中去的。在应用一次完成算法时，如果要求在每次新增观测数据后，接着就估计出系统模型的参数，则需要每次新增数据后要重新求解矩阵方程。最小二乘辩识方法在系统辩识领域中先应用上已相当普及，方法上相当完善，可以有效的用于系统的状态估计，参数估计以及自适应控制及其他方面。关键词：最小二乘（Least-squares），系统辨识（System Identifi

2、cation）目录： 1.目的1.1掌握系统辨识的理论、方法及应用1.2熟练Matlab下最小二乘法编程1.3掌握M序列产生方法2.设备 PC机1台（含Matlab软件） 3引言3.1 课题背景最小二乘理论是有高斯（K.F.Gauss）在1795年提出：“未知量的最大可能值是这样一个数值，它使各次实际观测值和计算值之间的差值的平方乘以度量其精度的数值以后的和最小。”这就是最小二乘法的最早思想。最小二乘辨识方法提供一个估算方法，使之能得到一个在最小方差意义上与实验数据最好拟合的数学模型。递推最小二乘法是在最小二乘法得到的观测数据的基础上，用新引入的数据对上一次估计的结果进行修正递推出下一个参数估

3、计值，直到估计值达到满意的精确度为止。4数学模型的结构辨识4.1根据汉格尔矩阵估计模型的阶次设一个可观可控的SISO过程的脉冲响应序列为个g(1),g(2),g(L),可以通过汉格尔（Hankel）矩阵的秩来确定系统的阶次。令Hankel阵为：，其中决定阵地维数，k可在1至间任意选择。则有。如果（过程的真实阶次），那么Hankel阵的秩等于。因此可以利用Hankel阵的奇异性来确定系统的阶次。4.2根据残差平方和估计模型的阶次SISO过程的差分方程模型的输出残差为，数据长度L，为阶时的数据矩阵，为阶时的参数的估计量，为模型阶次估计值，为真实阶次，则残差平方和函数：残差平方和有这样的性质：当L足

4、够大时，随着增加先是显著地下降，当>时，值显著下降的现象就终止。这就是损失函数法来定阶的原理。图1-1 损失函数法原理。 4.3根据AIC准则估计模型的阶次.图1-2 具体的定阶用法是：对不同阶次首先用极大似然法估计参数，然后计算似然函数值及值，找到使的作为。5 程序%待辨识系统 z(k)=0.1*z(k-1)-0.5*z(k-2)+0.6*z(k-3)+u(k-1)+0.5*u(k-2)-0.2*u(k-3)+v(k)/800%clcclear %清理工作间变量L=300; % M序列的周期x1=1;x2=1;x3=1;x4=0;x5=1;x6=0; %四个移位积存器的输出初始值for

5、 k=1:L; %开始循环，长度为L u(k)=xor(x3,x4); %第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或” x6=x5;x5=x4;x4=x3;x3=x2;x2=x1;x1=u(k);end %大循环结束，产生输入信号u plot(u) %绘图M序列v=randn(300,1); %随机误差干扰z=zeros(1,300);for k=4:300 z(k)=0.10*z(k-1)+0.55*z(k-2)+0.40*z(k-3)+0.70*u(k-1)+0.90*u(k-2)-3.50*u(k-3)+v(k)/400; %用理想输出值作为观测值endH=zeros(

6、300,6); %定义一个H“0”矩阵for i=4:300 H(i,:)=-z(i-1) -z(i-2) -z(i-3) u(i-1) u(i-2) u(i-3);%用循环产生H矩阵 z1(i,:)=z(i); %用循环产生z矩阵end%计算参数%c=inv(H'*H)*H'*z1%带入公式书上3.1.23a1=c(1),a2=c(2),a3=c(3),b1=c(4),b2=c(5),b3=c(6)%辨识出参数%系统阶次辨识AIC算法%bb=zeros(5,1);n=1; %假设为1阶for i=2:300 H1(i,:)=-z(i-1) u(i-1); zz1(i,:)=z

7、(i);endaa1=inv(H1'*H1)*H1'*zz1bb(1)=(zz1-H1*aa1)'*(zz1-H1*aa1)/L;AIC(1)=L*log(bb(1)+4*n;n=2; %假设为2阶for i=3:300 H2(i,:)=-z(i-1) -z(i-2) u(i-1) u(i-2); zz2(i,:)=z(i);endaa2=inv(H2'*H2)*H2'*zz2bb(2)=(zz2-H2*aa2)'*(zz2-H2*aa2)/L;AIC(2)=L*log(bb(2)+4*n;n=3; %假设为3阶for i=4:300 H3(i,

8、:)=-z(i-1) -z(i-2) -z(i-3) u(i-1) u(i-2) u(i-3); zz3(i,:)=z(i);endaa3=inv(H3'*H3)*H3'*zz3bb(3)=(zz3-H3*aa3)'*(zz3-H3*aa3)/L;AIC(3)=L*log(bb(3)+4*n;n=4; %假设为4阶for i=5:300 H4(i,:)=-z(i-1) -z(i-2) -z(i-3) -z(i-4) u(i-1) u(i-2) u(i-3) u(i-4); zz4(i,:)=z(i);endaa4=inv(H4'*H4)*H4'*zz4b

9、b(4)=(zz4-H4*aa4)'*(zz4-H4*aa4)/L;AIC(4)=L*log(bb(4)+4*n;n=5; %假设为5阶for i=6:300 H5(i,:)=-z(i-1) -z(i-2) -z(i-3) -z(i-4) -z(i-5) u(i-1) u(i-2) u(i-3) u(i-4) u(i-5); zz5(i,:)=z(i);endaa5=inv(H5'*H5)*H5'*zz5bb(5)=(zz5-H5*aa5)'*(zz5-H5*aa5)/L;AIC(5)=L*log(bb(5)+4*n;x=min(AIC)for i=1:5 if(AIC(i)=x) N=i %所辨识出的阶次N endendplot(1:5,AIC(1) AIC(2) AIC(3) AIC