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

实验4递推昀小二乘法的实现实验报告专业：U自动化班级：姓名：日期：年月日11．实验题目：U递推昀小二乘法的实现2．实验目的熟悉并掌握递推昀小二乘法的算法原理。3．实验主要原理给定系统12()(1)(2)()nykaykaykaykn01()(1)()()nbukbukbuknk（1）其中a1，a2…an，b0，b1，b2…bn为待辨识的未知参数，()k是不相关随机序列。y为系统的输出，u为系统的输入。分别测出n+N个输出、n+N输入值y(1)，y(2)…y(n+N)，u(1)，u(2)…u(n+N)，则可写出N个方程，具体写成矩阵形式，有10(1)()(1)(1)(1)(1)(2)(1)(2)(2)(2)(2)()(1)()()()()nnaynynyununaynynyununbynNynNyNunNuNnNb（2）设10(1)(1)(2)(2),,()()nnaynnaynnybynNnNb，()(1)(1)(1)(1)(2)(2)(2)(1)()()()ynyunuynyunuynNyNunNuN则式（2）可写为y（3）2式中：y为N维输出向量；为N维噪声向量；为2n+1维参数向量；为N×(2n+1)测量矩阵。为了尽量减小噪声对估值的影响，应取N＞2n+1，即方程数目大于未知数数目。的昀小二乘估计为1()TTy（4）为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。设已获得的观测数据长度为N，将式（3）中的y、和分别用,,NNNY来代替，即NNNY（5）用N表示的昀小二乘估计，则1TTNNNNNY（6）令1TNNNP，则TNNNNPY（7）如果再获得一组新的观测值(1)unN和(1)ynN，则又增加一个方程111TNNNy（8）式中11(1),(1)NNyynNnN1()(1)(1)(1)TNynNyNunNuN将式（5）和式（8）合并，并写成分块矩阵形式，可得T111NNNNNYy（9）于是，类似地可得到新的参数估值31TT1TTT1111NNNNNNNNNYyT1T11111NNNNNTNNNNNYPyPYy（10）式中1T1TT111T11NNNNNTNNNNP11T11NNNP（11）应用矩阵求逆引理，从求得1NP与NP的递推关系式出发，经过一系列的推导，昀终可求得递推昀小二乘法辨识公式：T1111NNNNNNKy（12）1T11111NNNNNNKPP（13）1TT111111NNNNNNNNNPPPPP（14）为了进行递推计算，需要给出NP和N的初值0P和0。推荐取值方法为：假定2000,PcI，c是充分大的常数，I为(21)(21)nn单位矩阵，则经过若干次递推之后能得到较好的参数估计。44．实验对象或参数给定系统12012()(1)(2)()(1)(2)()ykaykaykbukbukbukk（15）即2n。假设实际系统的参数为12a，21.3a，00.4b，10.88b，22.2b，但是不已知，即不可测。取()[0.1,0.1]k的零均值白噪声。输入信号取为()1.5sin0.2ukk（16）要求编制MATLAB程序，运用递推昀小二乘法对这一系统的参数进行在线辨识，并将辨识结果与实际参数进行对比。5．程序框图开始结束程序初始化与赋初值用户输入初值设定常数c设定初值θ0和P0用户输入所需递推次数rec_num生成零均值白噪声ξ(k)由系统模型计算获得yN+1满足递推次数？由递推昀小二乘辨识公式得到KN+1、PN+1、θN+1打印结果并绘制相关图形是否由先前输入输出生成ΨN+1将θN+1存入final_dataθN+1与实际值得相对误差存入final_percent根据状态反馈控制律计算增益系数矩阵K及控制U56．程序代码function[]=Lab_4()%--------------------------------------------------------实验题目及初始化定义disp('')disp('实验4递推昀小二乘法的实现')disp('')disp('**本实验对象系统的辨识模型为：')disp('y(k)=-a1*y(k-1)-a2*y(k-2)+b0*u(k)+b1*u(k-1)+b2*u(k-2)+ξ(k).')disp('假定实际参数：a1=2,a2=1.3,b0=0.4，b1=0.88,b2=2.2.')disp('')disp('**本实验系统输入信号取为u(k)=1.5*sin(0.2*k);输出受到零均值白噪声ξ(k)∈[-0.1,0.1]污染.')disp('')a0=65539;M=2147483647;x=123456;b=10000;Y(1)=0;Y(2)=0;U(1)=1.5*sin(0.2*1);U(2)=1.5*sin(0.2*2);U(3)=1.5*sin(0.2*3);V=[];k_num=2;T_NT=[2;1.3;0.4;0.88;2.2];c=input('**为了进行推算，初值设定时，假定P0=c*c*I，请输入充分大的常数c=');P_N0=c*c*eye(5);T_N0=zeros(5,1);rec_num=input('请输入所需递推步数rec_num=');disp('')rec_NUM=rec_num;final_data=zeros(rec_num,5);final_percent=zeros(rec_num,5);%--------------------------------------------------产生白噪声作为干扰信号ξ(k)fori=1:rec_num+2x=mod(a0*x+b,M);V(i)=0.2*(x/M-0.5);end%------------------------------------------------------递推昀小二乘法在线计算while(rec_num~=0)HL(1)=-Y(k_num);HL(2)=-Y(k_num-1);HL(3)=U(k_num+1);HL(4)=U(k_num);HL(5)=U(k_num-1);Y(k_num+1)=HL*[2;1.3;0.4;0.88;2.2]+V(k_num+1);K_NN=P_N0*HL'*(inv(1+HL*P_N0*HL'));P_NN=P_N0-K_NN*HL*P_N0;T_NN=T_N0+K_NN*(Y(k_num+1)-HL*T_N0);final_data(k_num-1,:)=T_NN;T_ND1=abs(T_NN-T_NT);T_ND2=[T_ND1(1,1)/T_NT(1,1)T_ND1(2,1)/T_NT(2,1)T_ND1(3,1)/T_NT(3,1)T_ND1(4,1)/T_NT(4,1)T_ND1(5,1)/T_NT(5,1)];final_percent(k_num-1,:)=T_ND2;P_N0=P_NN;T_N0=T_NN;k_num=k_num+1;K=[T_N0(2)-16T_N0(1)-4];U(k_num+1)=K*[Y(k_num);Y(k_num-1)]+[T_N0(3)T_N0(4)T_N0(5)]*[1.5*sin(0.2*(k_num+1));1.5*sin(0.2*k_num);1.5*sin(0.2*(k_num-1))]6rec_num=rec_num-1;end%-----------------------------------------------显示识别结果并与实际值进行比较disp('**运用递推昀小二乘法对这一系统参数进行辨识，计算得到:')disp('')disp('相对于理论值[a1a2b0b1b2]=[21.30.40.882.2]，各次递推结果对应如final_data矩阵所示：')final_datadisp('实际值相对于理论值的误差如final_percent矩阵所示:')final_percent%-------------------------------------------------------------------相关绘图k=0:1:rec_NUM-1;subplot(3,2,1)stairs(k,final_data(:,1)','k-');holdonline([0rec_NUM-1],[T_NT(1,1)T_NT(1,1)])axis([0rec_NUM-1-0.52.5])title('a_1收敛情况图形')subplot(3,2,3)stairs(k,final_data(:,2)','k-');holdonline([0rec_NUM-1],[T_NT(2,1)T_NT(2,1)])axis([0rec_NUM-1-0.21.5])title('a_2收敛情况图形')subplot(3,2,2)stairs(k,final_data(:,3)','k-');holdonline([0rec_NUM-1],[T_NT(3,1)T_NT(3,1)])axis([0rec_NUM-10.11.4])title('b_0收敛情况图形')subplot(3,2,4)stairs(k,final_data(:,4)','k-');holdonline([0rec_NUM-1],[T_NT(4,1)T_NT(4,1)])axis([0rec_NUM-10.61.3])title('b_1收敛情况图形')subplot(3,2,6)stairs(k,final_data(:,5)','k-');holdonline([0rec_NUM-1],[T_NT(5,1)T_NT(5,1)])axis([0rec_NUM-1-0.52.5])title('b_2收敛情况图形')end7．实验结果及分析本实验给定系统为12012()(1)(2)()(1)(2)()ykaykaykbukbukbukk要求运用递推昀小二乘法对所给系统参数进行在线识别，其主要目的是设计状态反馈控制7律，因此首先采用状态空间分析法对系统进行分析。闭环系统原理图如下图所示。根据系统识别模型，对于该二阶系统，不妨设状态量1()xyk，21(1)xxyk，所以该系统状态空间表达式为：12(1)0(2)1xykXAXBVxykYCX其中：2101Aaa，01B，10C，012(2)(1)()ukVbbbukuk