(完整word版)2003版系统辨识最小二乘法大作业

系统辨识最小二乘法作者：王景西北工业大学系统辩识大作业题目：最小二乘法系统辨识系统辨识最小二乘法作者：王景一、问题重述：用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数离散化有z^4-3.935z^3+5.806z^2-3.807z+0.9362----------------------------------------------=z^4-3.808z^3+5.434z^2-3.445z+0.8187噪声的成形滤波器离散化有4.004e-010z^3+4.232e-009z^2+4.066e-009z+3.551e-010-----------------------------------------------------------------------------=z^4-3.808z^3+5.434z^2-3.445z+0.8187采样时间0.01s要求：1.用Matlab写出程序代码；2.画出实际模型和辨识得到模型的误差曲线；3.画出递推算法迭代时各辨识参数的变化曲线；最小二乘法:在系统辨识领域中,最小二乘法是一种得到广泛应用的估计方法,可用于动态,静态,线性,非线性系统。在使用最小二乘法进行参数估计时,为了实现实时控制,必须优化成参数递推算法,即最小二乘递推算法。这种辨识方法主要用于在线辨识。MATLAB是一套高性能数字计算和可视化软件,它集成概念设计,算法开发,建模仿真,实时实现于一体,构成了一个使用方便、界面友好的用户环境,其强大的扩展功能为各领域的应用提供了基础。对4324326.51411.5320120232320YssssGUssss432120120232320ENWssss系统辨识最小二乘法作者：王景于一个简单的系统,可以通过分析其过程的运动规律,应用一些已知的定理和原理,建立数学模型,即所谓的“白箱建模”。但对于比较复杂的生产过程,该建模方法有很大的局限性。由于过程的输入输出信号一般总是可以测量的,而且过程的动态特性必然表现在这些输入输出数据中,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。这种建模方法就称为系统辨识。把辨识建模称作“黑箱建模”。系统辨识又分为参数辨识和阶次辨识,在本文中只讨论参数辨识问题最小二乘递推算法所用的模型：Z(k)=B()u(k)+v(k)最小二乘递推算法为:)(kv是服从N)1,0(分布的不相关随机噪声。)(1zG)()(11zAzB，)(1zN)()(11zCzD，考虑如图1示仿真对象,系统的差分方程为z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+)(kv(3.69)仿真对象选择如下的模型结构)()4(5)3()2(3)1()()4(4)3(3)2()1()(42121kvkubkubkubkubkubkzakzakzakzakz(3.68)++e(k)图1最小二乘递推算法辨识实例结构图y(k)u(k)z(k)v(k))(1zN)(1zG系统辨识最小二乘法作者：王景其中，)(kv是服从正态分布的白噪声N)1,0(。输入信号采用4位移位寄存器产生的M序列，幅度为1。按式(3.69)构造h(k)；加权阵取单位阵IΛL；利用式(3.61)计算K(k)、)(ˆkθ和P(k)，计算各次参数辨识的相对误差，精度满足要求式(3.67)后停机。最小二乘递推算法辨识的Malab6.0程序流程图：系统辨识最小二乘法作者：王景最小二乘递推算法辨识程序clear%清理工作间变量L=35;%M序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值fori=1:L;%开始循环，长度为Lx1=xor(y3,y4);%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”x2=y1;%第二个移位寄存器的输入是第一个移位寄存器的输出x3=y2;%第三个移位寄存器的输入是第二个移位寄存器的输出x4=y3;%第四个移位寄存器的输入是第三个移位寄存器的输出y(i)=y4;%取出第四个移位寄存器的幅值为0和1的输出信号，即M序列ify(i)0.5,u(i)=-1;%如果M序列的值为1,辨识的输入信号取“-1”elseu(i)=1;%如果M序列的值为0,辨识的输入信号取“1”end%小循环结束y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号ufigure(1);%第一个图形stem(u),gridon%显示出输入信号径线图并给图形加上网格z(2)=0;z(1)=0;z(3)=0;z(4)=0;%设z的前四个初始值为零fork=5:25;%循环变量从5到25z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4);%输出采样信号end%RLS递推最小二乘辨识c0=[0.0010.0010.0010.0010.0010.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(9,9);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%取相对误差E=0.000000005c=[c0,zeros(9,24)];%被辨识参数矩阵的初始值及大小e=zeros(9,25);%相对误差的初始值及大小fork=5:25;%开始求Kh1=[-z(k-1),-z(k-2),-z(k-3),-z(k-4),u(k),u(k-1),u(k-2),u(k-3),u(k-4)]';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)的值p0=p1;%给下次用ife2=Ebreak;%如果参数收敛情况满足要求，终止计算end%小循环结束end%大循环结束c,e%显示被辨识参数及其误差(收敛)情况%分离参数a1=c(1,:);a2=c(2,:);a3=c(3,:);a4=c(4,:);b1=c(5,:);b2=c(6,:);b3=c(7,:);b4=c(8,:);b5=c(9,:);ea1=e(1,:);ea2=e(2,:);ea3=e(3,:);ea4=e(4,:);eb1=e(5,:);eb2=e(6,:);eb3=e(7,:);eb4=e(8,:);eb5=e(9,:);figure(2);%第二个图形i=1:25;%横坐标从1到25plot(i,a1,'r',i,a2,':',i,a3,'r',i,a4,'k',i,b1,'y',i,b2,':',i,b3,'m',i,b4,':',i,b5,'g')%画出a1，a2，a3,a4,b1，b2,b3,b4,b5的各次辨识结果title('参数变化曲线')%图形标题figure(3);%第三个图形i=1:25;%横坐标从1到25plot(i,ea1,'r',i,ea2,'k:',i,ea3,'y',i,ea4,'m',i,eb1,'g',i,eb2,'c:',i,eb3,'r',i,eb4,':',i,eb5,'g')%画出a1，a2，a3,a4,b1，b2,b3,b4,b5的各次辨识结果的收敛情况title('误差曲线')%图形标题系统辨识最小二乘法作者：王景参数变化图相对误差图仿真结果表明，大约递推到第十五步时，参数辨识的结果基本达到稳定状态，即a1=-3.808,a2=5.434,a3=-3.445,a4=-0.8187,b1=1,b2=-3.935b3=5.806b4=-3.807b5=0.9362。此时，参数的相对变化量E0.000000005。从整个辨识过程来看，精度的要求直接影响辨识的速度。虽然最终的精度可以达到很小，但开始阶段的相对误差会很大，从相对误差图形中可看出，参数的最大相对误差会达到三位数系统辨识最小二乘法作者：王景增广最小二乘法算法所用的模型：Z(k)=B()u(k)+D()e(k)增广最小二乘法算法为：模型结构选用如下形式：)4)-v(k*4d3)-v(k*3d2)-v(k*2d1)-v(k*1d)4(5)3()2(3)1()()4(4)3(3)2()1()(42121kubkubkubkubkubkzakzakzakzakz系统辨识最小二乘法作者：王景增广最小二乘算法流程图如下页图所示系统辨识最小二乘法作者：王景增广最小二乘辨识程序clearL=55;%4位移位寄存器产生的M序列的周期y1=1;y2=1;y3=1;y4=0;%4个移位寄存器的输出初始值fori=1:L;x1=xor(y3,y4);%第一个移位寄存器的输入信号x2=y1;%第二个移位寄存器的输入信号x3=y2;%第三个移位寄存器的输入信号x4=y3;%第四个移位寄存器的输入信号y(i)=y4;%第四个移位寄存器的输出信号，M序列,幅值0和1，ify(i)0.5,u(i)=-1;%M序列的值为1时,辨识的输入信号取“-1”elseu(i)=1;%M序列的值为0时,辨识的输入信号取“1”endy1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号准备endfigure(1);%画第一个图形subplot(2,1,1);%画第一个图形的第一个子图stem(u),gridon%画出M序列输入信号v=randn(1,60);%产生一组60个正态分布的随机噪声subplot(2,1,2);%画第一个图形的第二个子图plot(v),gridon;%画出随机噪声信号u,v%显示输入信号和噪声信号z=zeros(13,55);%输出采样矩阵z(2)=0;z(1)=0;z(3)=0;z(4)=0;%输出采样的初值c0=[0.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(13,13);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.00000000005;%相对误差E=0.000000005c=[c0,zeros(13,54)];%被辨识参数矩阵的初始值及大小e=zeros(13,55);%相对误差的初始值及大小fork=5:55;%开始求Kz(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+4.004e-010*v(k-1)+4.232e-009*v(k-2)+4.066e-009*v(k-3)+3.551e-010*v(k-4);%系统在M序列输入下的输出采样信号h1=[-z(k-1),-z(k-2),-z(k-3),-z(k-4),u(k),u(k-1),u(k-2),u(k-3),u(k-4),v(k-1),v(k-