使用最小二乘法法进行系统辨识的两种方法

1、递推最小二乘法辨识与仿真现在有如下的辨识仿真对象：w(k)e(k)+y(k)u(k)图1递推最小二乘算法辨识G(z-1)N(z-)1)图中，v(k)是服从N(0,1)分布的不相关随机噪声。且G(z1)=B(7-1)，N(z1)=D(z-1)A(z1)C(z1)A(z1)=11.5a1z1+0.7z2=C(z1)B(z1)=1.0z1+0.5z2D(z1)=1选择上图所示的辨识模型。仿真对象选择如下的模型结构：y(k)+ay(k1)+ay(k2)=bu(k1)+bu(k2)+w(k)1212可得系统模型为：y(k)1.5y(k1)+0.7y(k2)=u(k1)+0.5u(k2)+w(k)递推最小

2、二乘法的推导公式如下：6(k)=6(k1)+K(k)z(k)he(k)6(k1)P(k)=P-1(k1)+h(k)he(k)-1=IK(k)ht(k)P(k1)K(k)=P(k1)h(k)he(k)P(k1)h(k)+1-1相关程序如下：%exp053%递推最小二乘法程序%clear%清理工作间变量L=55;%M序列的周期yl=l;y2=l;y3=l;y4=0;%四个移位寄存器的输出初始值fori=l:L;%开始循环，长度为Lxl=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或”x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出x3=y2;%第三

3、个移位积存器的输入是第2个移位积存器的输出x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出y(i)=y4;%取出第四个移位积存器幅值为0和1的输出信号，ify(i)0.5,u(i)=-0.03;%如果M序列的值为1时,辨识的输入信号取“-0.03”elseu(i)=0.03;%当M序列的值为0时,辨识的输入信号取“0.03”end%小循环结束yl=xl;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号uw=normrnd(0,sqrt(0.1),1,55);%加入白噪声figure(1);%第1个图形，伪随机序列stem(u),grid

4、on%以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;%取z的前两个初始值为零fork=3:55;%循环变量从3到55z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+w(k);%给出理想的辨识输出采样信号endc0=0.0010.0010.0010.001,0.001;%直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(5,5);%直接给出初始状态P0,即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=c0,zeros(5,54);%被辨识参数矩阵的初始值及大小e=ze

5、ros(5,55);%相对误差的初始值及大小fork=3:55;%开始求Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2),w(k);x=h1*p0*h1+1;x1=inv(x);%开始求K(k)k1=p0*h1*x1;%求出K的值d1=z(k)-h1*c0;c1=c0+k1*d1;%求被辨识参数ce1=c1-c0;%求参数当前值与上一次的值的差值e2=e1./c0;%求参数的相对变化e(:,k)=e2;%把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;%新获得的参数作为下一次递推的旧参数c(:,k)=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列p1=p0-k1

6、*k1*h1*p0*h1+1;%求出p(k)的值p0=p1;%给下次用ife2=Ebreak;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束c;%显示被辨识参数e;%显示辨识结果的收敛情况%分离参数a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);d1=c(5,:);ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);figure(2);%第2个图形i=l:55;%横坐标从1到55plot(i,a1,r,i,a2,:,i,b1,g,i,b2,:,i,b1,k)%画出a1，a2，b1，b2的各次辨识结果titl

7、e(系统辨识结果)图形标题A(z-1)=11.5az-1+0.7z-2=C(z-1)1B(z-1)=1.0z-1+0.5z-2D(z-1)=1一z-1+0.2z-2那么系统模型机构为：y(k)一1.5y(k一1)+0.7y(k一2)=u(k一1)+0.5u(k一2)+w(k)一w(k一1)+0.2w(k一2)相关程序如下：%exp054.m%递推最小二乘法编程%clear%清理工作间变量L=55;%M序列的周期yl=l;y2=l;y3=l;y4=0;%四个移位寄存器的输出初始值fori=1:L;%开始循环，长度为Lxl=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的

8、输出的“或”x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出y(i)=y4;%取出第四个移位积存器幅值为0和T的输出信号，ify(i)0.5,u(i)=-0.03;%如果M序列的值为1时,辨识的输入信号取“-0.03”elseu(i)=0.03;%当M序列的值为0时,辨识的输入信号取“0.03”end%小循环结束yl=xl;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号uw=normrnd(0,sqrt(0.1),1

9、,55);figure(1);%第1个图形，伪随机序列stem(u),gridon%以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;%取z的前两个初始值为零fork=3:55;%循环变量从3到55z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+w(k)-w(k-1)+0.2*w(k-2);%给出理想的辨识输出采样信号endc0=0.0010.0010.0010.001,0.001;%直接给出被辨识参数的初始值，即一个充分小的实向量p0=106*eye(5,5);%直接给出初始状态P0,即一个充分大的实数单位矩阵E=0.0000000

10、05;%相对误差E=0.000000005c=c0,zeros(5,54);%被辨识参数矩阵的初始值及大小e=zeros(5,55);%相对误差的初始值及大小fork=3:55;%开始求Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2),w(k),w(k-1),w(k-2);x=h1*p0*h1+1;x1=inv(x);%开始求K(k)k1=p0*h1*x1;%求出K的值d1=z(k)-h1*c0;c1=c0+k1*d1;%求被辨识参数ce1=c1-c0;%求参数当前值与上一次的值的差值e2=e1./c0;%求参数的相对变化e(:,k)=e2;%把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;%新获得的参数作为下一次递推的旧参数c(:,k)=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1*h1*p0*h1+1;%求出p(k)的值pO=pl;%给下次用fe2=Ebreak;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束c;%显示被辨识参数e;%显示辨识结果的收敛情况%分离参数a1=c(1