实验二离散线性系统参数估计的递推最小二乘法 - Read

实验二离散线性系统参数估计的递推最小二乘法 - Read 实验二 离散线性系统参数估计的递推最小二乘法 实验目的 1) 掌握利用递推最小二乘法估计离散线性系统参数的基本 步骤 新产品开发流程的步骤课题研究的五个步骤成本核算步骤微型课题研究步骤数控铣床操作步骤 和使用要点。 2) 通过基本递推最小二乘(RLS)法，辅助变量(IV)法和增广最小二乘(ELS)法等三种常 用方法，了解最小二乘类型的递推估计算法的一般编程思路。 3) 熟悉估计算法中各个参数选择对辨识效果的影响。 实验原理 1) 辩识系统的基本情况和辩识工作的要求 待辩识系统的运行情况如图1所示，设其在静态工作点(U，Y)附近作局部线性化00 所得动态模型为 y+ay+…+ay=bu+…+bu，y= …= y=0(1) (k)1(k-1)n(k-n)1(k-1-d)n(k-n-d) (1-n)(0) 式中阶次n和时延d为结构性参数(一般需事先给出)，a，b为待估非结构性参数。 ii u (k) + u U0(k)Y=Y+y (k)0(k)待辩识系统 U 0 η (k) ζ (k) =U+u +ζ U’Y’=Y+y +η (k)0(k)(k) (k)0(k)(k) 图1 待辩识系统的基本运行情况 一般来说，采用伪随机信号M序列作为系统动态测试输入信号可以保证最小二乘法的存在并易于取得较好的参数估计结果。这是因为M序列的自相关函数很接近于周期性的函数，其功率谱密度中含有较宽范围的频率成分，族以激发起系统在各种频率下的动态行为。 本实验中取电平幅值V=1、周期为15的M序列信号(其波形见图2)作为动态测试信号u施加于对象输入端。本次实验的基本任务就是在含有噪声的输入、输出的实际观测值(k) U’=U+u +ζ和Y’=Y+y +η的基础上, 用最小二乘算法递推估计出系统参数(k)0(k)(k) (k)0(k)(k) { a, …,a ,b,…b}。 1n1n u (k) 01 2 3 4 56789101112131415 1 51 1 1 1 1 1 1 1 1 1 1 图2 M序列信号的波形 2) 基本关系式 取考虑噪声后的对象模型 ''''' y,,ay,？,ay，bu，？，bu，,(k)1(k,1)n(k,n)1(k,1,d)n(k,n,d)(k) , (2) ,,,，,k(k) 式中 '' (3) y,y，,,u,u，,(k)(k)(k)(k)(k)(k) ,'''',,[,y,？,,y|,u,？,,u] (4) ,,,,,,k(k1)(kn)(k1d)(knd),,,[,a,？,,a|b,？,b] (5) 1212 (6) ,,a,，？，a,,b,,？,b,(k)1(k,1)n(k,n)1(k,1,d)n(k,n,d) 3) 基本型最小二乘(LS)算法简介 (1) 批量算法(Batch algorithm): 取拟合残差 ˆˆˆˆˆe,,,,a，？，a,,b,,？,b, (7) (k)(k)1(k,1)n(k,n)1(k,1,d)n(k,n,d)的加权平方和作为参数估计的性能指标 NN 2',2ˆˆJ(,),we,w[y,,,] (8) ,,Nkkkkk()(),,kk11ˆw,0,k,1,2,？,N对参数估计量,作极小化可在条件下得到加权最小二乘估计 k ,,1,ˆ,,(,w,),wy (9) NWLSNNNNNN, 式中 ,',，，，，y11,,,,,',y22,,,,,,,y (10) ,NN,,,,？？,,,,,',y,,,,NN，，，， (11) W,diag[w,w,？,w]N12N (2) 递推算法(Recursive algorithm): 为了跟踪系统可能有的参数漂移性变化，按渐消记忆法取记忆因子为α=0.95~1。此时有 N,1 (12) W,diag[,,？,,,1]N且 ,,yW0，，，，，，NNn,Wy (13) ,,,,,N，1N，1N，1,,,,,,',,y01N，1N，1，，，，，，进一步可推得基本型最小二乘估计算法为 ',ˆˆˆ,,,,,P[y,,],，,N，1NN，1NN，1N，1N，1N ,,,[1,,P,],，N，1N，1NN，1, (14) ,,,P[P,P,,P],,,N，1NN，1NN，1N，1N, 初值设定为 2ˆ (15) ,,0,P,cI002n ''由于不是白噪声，且与和有关，以上基本型最小二乘估计结果是有偏的。,uy(k)(k)(k) '为了克服这种有偏性，在(即)情况下可采用辅助变量法(IV)，在较,,0u,u(k)(k)(k) 一般情况下可采用增广最小二乘法(ELS)。 4) 辅助变量(IV)估计算法简介 辅助变量(IV)估计批量算法为 ,,1,ˆ,,(ZW,)ZWy (16) NIVNNNNNN, 式中 ,,Zzzz,[,,？,]12NN,'',zyyuu,,,[,？,|,？,] (17) ,(1)()(1)()kvk,vk,nk,,dk,n,d ,yz,,()vkkk, ˆ采用,的延迟多步的低通滤波值，目的是使辅助变量与方,y,i,1,2,？,nk,IV(k)v(k,i)程误差尽量不相关。IV估计的递推算法在形式上与RLS算法相似，且因对初值设定敏,(k) 感，一般宜先以RLS算法启动，待约50步后在切换到IV法。辅助变量(IV)法估计的递推算 法框图见图4。 5) 增广最小二乘(ELS)估计算法简介 增广最小二乘(ELS)估计是根据考虑噪声后的对象模型的如下变形: ''''' y,,ay,？,ay，bu，？，bu,,,,,,(k)1(k1)n(kn)1(k1d)n(knd) ，,，,,，？，,,(k)1(k,1)n，d(k,n,d) , (18) ,,,，,k(k) 式中为白噪声， ,(k)'''', (19) ,,[,y,？,,y|u,？,u|,,？,,]k(k,1)(k,n)(k,1,d)(k,n,d)(k,1)(k,n,d) (20) ,,[,a,？,,a|b,？,b|,,？,,]1n1n1n，d相应的批量算法为 ,,1,ˆ,,(,W,),Wy (21) NELSNNNNNN, 式中 , (22) ,,[,,,,？,,]N12n ,但包含噪声，只能用递推方式估算如下: ,,i,0k(k,i) 'ˆ,,,,ˆ,,y,(k)(k)kk (23) ,'''',ˆ,,[,y,？,,y|u,？,u|,,？,,],k(k,1)(k,n)(k,1,d)(k,n,d)(k,1)(k,n,d), 增广最小二乘(ELS)估计的递推算法框图见图5。 实验步骤 1) 辩识对象运行情况的仿真 本实验取系统(1)的参数真值为。 n,2,d,1,a,,1,a,0.5,b,1,b,0.51212 u将M序列作为，即可按(1)式求得。再设工作点(U,Y)，并y,k,1,2,3,？00(k)(k) ’,,和，则对象的实际输入、输出U和Y可按(2)式求得。其中正态白生成随机噪声(k)(k)kk ,,噪声序列和根据中心极限定理采用12个在[-0.5,0,5]中均匀分布的随机数叠加的方kk 法生成，其方差可由噪声幅度因子调节。 在待辩识对象在大范围内是非线性的情况下，通过辩识得出的是对象在工作点附近的局部线性化动态模型，因此还应去处工作点Y，以得出输出端动态分量的实测值。但工作点0 无法单独测量，只能根据最新得到的推算: N,,N(,,3~4)个Y'p(k) N,1N,111ˆˆU,U'Y,Y' , (24) ,,0(k)(k,1)0(k)(k,1)NNi,0i,0 或递推计算如下: 1ˆˆ,U,U，U',U',,,0(k)0(k,1)(k)(k,N)N (25) ,1ˆˆ,,YYY'Y',，,,0(k)0(k,1)(k)(k,N)N, ''''式中的初值均可设为0。 U,YU,Ykk0()0()0(0)0(0) 于是，输出端动态分量的实测值可近似估算为 ''',u,U,U,(k)(k)0(k) (26) ,'''yYY,,,(k)(k)0(k), 2) 模型的阶和时延检查 为简单起见，可只用基本LS法来检验系统模型的阶次n和时延d。一个有效的作法是比较不同时延d和不同阶次n情况下得到的模型对于观测数据的拟和优良度。拟和优良度可用递腿参数估计已进入收敛的第N步所得残差平方和 2 N'2,ˆ (27) J,[y,,,](),1,kkk,，11kN 来衡量。为此，需作出不同时延d值下J随阶次n变化的折线如图6所示。一般来说，随阶次n增大而减小。能使J取最小值的那条折线所对应的d值，即可认为是合适的时延估计值ˆd。在此折线上，当n变得比真正的阶次n大时J几乎停止减小。因此，确定阶次的方法0 ˆ就是取最后一次明显地陡降后的n值作为合适的阶次估计值。 n 时延d也可通过在系统已稳定运行于稳态工作点()情况下，在输入端叠加阶跃或冲激信号，观测输出端开始出现响应的延迟时间来直接确定。 J ˆ d,d ˆ d,d ˆ d,d 01 2 3 4 5 6 n 图6 阶次和时延的检验 3)递推最小二乘算法，辅助变量(IV)估计算法，增广最小二乘(ELS)估算法流程图 基本递推最小二乘(RLS)算法的框图见图3。 初始化 ’’生成数据 u, Y, k=1-n-d,…,-1,0 (k-d) (k) 2ˆ设定初值为 ,,0,P,cI002n 1 ? k ,'' ,,[,y,？,,y|构造向量 ,,k(k1)(kn) '',u ,？,,u] ,,,,(k1d)(knd) ? k +1 k ''生成当前新数据 u,yk,dk()() ,,,Pkk,1k ,,,1,[,，,] kkk ' ˆˆˆ,,,,，,,[y,,,] kk,1kk(k)kk,1 tP,[P,,,,],kk,1kkk k-k>iT? 1p 打印结果 图3 基本递推最小二乘算法流程图 初始化 ’’生成数据 u, Y, k=1-n-d,…,-1,0 (k-d) (k) 2ˆ设定初值为(c充分大) ,,0,P,cI002n ,,0,y,0,k,00v(k) 1 ? k ˆ滤波 ,,(1,,),，,,,,,0.01~0.1,q,0~10kk,1k,1,q N k>50? Y '',''构造向量 ,u,？,,u],,[,y,？,,y|,,,,(k1d)(knd),,k(k1)(kn) '', ,u,？,,u]z,[,y,？,,y|,,,,(k1d)(knd)kv(k,1)v(k,n) ? k +1 k, ''u,y,y,z,生成当前新数据 k(,d)k()vk()kk 同RLS算 法中的Q ψ,P,P段程序 ,,,,kk1kkk1k 计算，zτγ,1[α，ψ]kkk, 'τˆˆˆθ,θ，γψ[y,θ],,kk1kk(k)kk1, tP,[P,γψψ]α,kk1kkk N k-k>iT? 1p 打印结果 图4 递推辅助变量(IV)估计算法流程图 初始化 ’’生成数据 u, Y, k=1-n-d,…,-1,0 (k-d) (k) 3 ˆ设定初值为(c充分大) ,,0,P,cI003n，d 1 ? k ,'' ,,[,y,？,,y|构造向量 ,,k(k1)(kn)'',u,？,,u|v,？,v] (k,1,d)(k,n,d)(k,1)(k,n,d) k? k +1 ''', ˆˆ生成当前新数据uyv,y,,, ,,(k,d)(k)(k)(k)kk,1 ,,,P kk,1k ,,1,[,，,]kkk ,ˆˆ ,,,，,,vkk,1kk(k) t P,[P,,,,],kk,1kkk k-k>iT? 1p 打印结果 图5 递推增广最小二乘(ELS)估计算法流程图 4)实验结果 递推最小二乘算法仿真图 递推辅助变量(IV)估计算法仿真图 递推增广最小二乘(ELS)估计算法仿真图 用基本递推最小二乘法(RLS)作递推估计，设定阶次及延时(n=2，d=1)，观测遗忘因子α、P(i,i)的初值C的不同取值对系统收敛速度和估计精度的影响。 1 α的选取对系统辩识的影响: 选取α=0.95，p(i,i)=10000时的仿真曲线 结论:α的增加有助于加快收敛速度 2 观察p(i,i)的改变对辨识结果的影响 当α=0.995，p(i,i)=8000时的曲线 结论:p(i,i)的增大有助于提高收敛速度。 (2) 在有噪声(aa=0.1，bb=0.4)，和正确设定阶次及延时(n=2，d=1)的情况下，观察当系统参数时变时，采用RLS法，比较不同的遗忘因子α对参数变化的跟踪效果的影响。 在aa=0.1，bb=0.4的情况下，α,0.95时的曲线 在α,0.98时的曲线 遗忘因子较大时对参数的跟踪效果较好，因为遗忘因子越大时各参数权值越相近，越接 近真实情况。 附录源程序 1. 递推最小二乘算法 clear all; close all; clc; % M array construction Np=15;r=4; X1=1;X2=1;X3=1;X4=1; m_length = r*Np; a=1; for i=1:1:m_length Y4=X4;Y3=X3;Y2=X2;Y1=X1; X4=Y3;X3=Y2;X2=Y1; X1=xor(Y3,Y4); if Y4==0 M(i)=-a; else M(i)=a; end end figure; i=1:1:m_length; plot(i,M); % 白噪声 noise = zeros(1,m_length); for i=1:1:m_length temp = noise + 0.5*rands(1,m_length); noise = temp; end noise = noise/12; %noise = temp; % parameter of system n=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5; S_U0=0.2;S_Y0=0.2; % generate u,y u0=ones(1,m_length)*S_U0; U = M + u0 + noise; figure; i=1:1:m_length; plot(i,U); %formulation y(1)=0;y(2)=0;y(3)=0; Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3); for k=4:m_length y(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2); Y(k)=S_Y0+y(k)+noise(k); end figure; i=1:1:m_length; plot(i,Y); %辨识 % premitive value c=0.5; P = (c^2)*eye(2*n); sita(:,3) = [0,0,0,0]'; alf = 0.995; % M %compute U0,Y0 sum_U = 0;sum_Y = 0; for k=1:1:m_length sum_U = sum_U + U(k); sum_Y = sum_Y + Y(k); end t_U0 = sum_U/m_length;t_Y0 = sum_Y/m_length; for k=1:1:m_length t_u(k) = U(k) - t_U0; t_y(k) = Y(k) - t_Y0; end figure; i=1:1:m_length; plot(i,t_u); figure; i=1:1:m_length; plot(i,t_y); for k=4:1:m_length %sita_f = sita; fai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d)]'; W = P*fai; dan = 1/(alf + fai'*W); sita(:,k) = sita(:,k-1) + dan*W*(t_y(k)-fai'*sita(:,k-1)); P = ( P - dan * W*W')/alf; end figure; i = 3:1:m_length; plot(i,sita(1,3:m_length),i,sita(2,3:m_length),i,sita(3,3:m_length),i,sita(4,3: m_length)); 2.递推辅助变量(IV)估计算法 clear all; close all; clc; % M array construction Np=15;r=7; X1=1;X2=1;X3=1;X4=1; m_length = r*Np; a=1; for i=1:1:m_length Y4=X4;Y3=X3;Y2=X2;Y1=X1; X4=Y3;X3=Y2;X2=Y1; X1=xor(Y3,Y4); if Y4==0 M(i)=-a; else M(i)=a; end end figure; i=1:1:m_length; plot(i,M); % 白噪声 noise = zeros(1,m_length); for i=1:1:m_length temp = noise + 0.5*rands(1,m_length); noise = temp; end noise = noise/12; %noise = temp; % parameter of system n=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5; S_U0=0.2;S_Y0=0.2; % generate u,y u0=ones(1,m_length)*S_U0; U = M + u0 + noise; figure; i=1:1:m_length; plot(i,U); %formulation y(1)=0;y(2)=0;y(3)=0; Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3); for k=4:m_length y(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2); Y(k)=S_Y0+y(k)+noise(k); end figure; i=1:1:m_length; plot(i,Y); %辨识 % premitive value c=2; P = (c^2)*eye(2*n); sita(:,3) = [0,0,0,0]'; alf = 0.95; % M %compute U0,Y0 sum_U = 0;sum_Y = 0; for k=1:1:m_length sum_U = sum_U + U(k); sum_Y = sum_Y + Y(k); end t_U0 = sum_U/m_length;t_Y0 = sum_Y/m_length; for k=1:1:m_length t_u(k) = U(k) - t_U0; t_y(k) = Y(k) - t_Y0; end figure; i=1:1:m_length; plot(i,t_u); figure; i=1:1:m_length; plot(i,t_y); % ro = 0.01; % 0.01 ~ 0.1 q = 1; % 0 ~ 10 l_sita = [0,0,0,0]'; yv(1)=0;yv(2)=0; for k=4:1:q+4 %sita_f = sita; l_sita = (1-ro)*l_sita + ro*sita(k-1); fai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d)]'; W = P*fai; dan = 1/(alf + fai'*W); sita(:,k) = sita(:,k-1) + dan*W*(t_y(k)-fai'*sita(:,k-1)); P = ( P - dan * W*W')/alf; end for k=q+5:1:m_length l_sita = (1-ro)*l_sita + ro*sita(:,k-1); fai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d)]'; if k<51 W = P*fai; dan = 1/(alf + fai'*W); sita(:,k) = sita(:,k-1) + dan*W*(t_y(k)-fai'*sita(:,k-1)); P = ( P - dan * W*W')/alf; else z = [-yv(k-48-1),-yv(k-48-2),t_u(k-1-d),t_u(k-2-d)]'; yv(k-48)=z'*l_sita; W = P*fai; Wk = P*z; dan = 1/(alf + fai'*W); dank = 1/(alf + fai'*P*z); sita(:,k) = sita(:,k-1) + dank * Wk * (t_y(k)-fai'*sita(:,k-1)); P = ( P - dank * Wk * Wk')/alf; end end figure; i = 3:1:m_length; plot(i,sita(1,3:m_length),i,sita(2,3:m_length),i,sita(3,3:m_length),i,sita(4,3: m_length)); 3.递推增广最小二乘(ELS)估计算法 clear all; close all; clc; % M array construction Np=15;r=4; X1=1;X2=1;X3=1;X4=1; m_length = r*Np; a=1; for i=1:1:m_length Y4=X4;Y3=X3;Y2=X2;Y1=X1; X4=Y3;X3=Y2;X2=Y1; X1=xor(Y3,Y4); if Y4==0 M(i)=-a; else M(i)=a; end end figure; i=1:1:m_length; plot(i,M); % 白噪声 noise = zeros(1,m_length); for i=1:1:m_length temp = noise + 0.5*rands(1,m_length); noise = temp; end noise = noise/12; %noise = temp; % parameter of system n=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5; S_U0=0.2;S_Y0=0.2; % generate u,y u0=ones(1,m_length)*S_U0; U = M + u0 + noise; figure; i=1:1:m_length; plot(i,U); %formulation y(1)=0;y(2)=0;y(3)=0; Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3); for k=4:m_length y(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2); Y(k)=S_Y0+y(k)+noise(k); end figure; i=1:1:m_length; plot(i,Y); %辨识 % premitive value c=2; P = (c^3)*eye(3*n+d); sita(:,3) = [0,0,0,0,0,0,0]'; alf = 0.995; % M %compute U0,Y0 sum_U = 0;sum_Y = 0; for k=1:1:m_length sum_U = sum_U + U(k); sum_Y = sum_Y + Y(k); end t_U0 = sum_U/m_length;t_Y0 = sum_Y/m_length; for k=1:1:m_length t_u(k) = U(k) - t_U0; t_y(k) = Y(k) - t_Y0; end figure; i=1:1:m_length; plot(i,t_u); figure; i=1:1:m_length; plot(i,t_y); v(1)=0;v(2)=0;v(3)=0; for k=4:1:m_length fai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d),v(k-1),v(k-2),v(k-3)]'; v(k) = t_y(k) - fai'*sita(:,k-1); W = P*fai; dan = 1/(alf + fai'*W); sita(:,k) = sita(:,k-1) + dan*W*v(k); P = ( P - dan * W*W')/alf; end figure; i = 3:1:m_length; plot(i,sita(1,3:m_length),'r',i,sita(2,3:m_length),'g',i,sita(3,3:m_length),'b' ,i,sita(4,3:m_length),'y',i,sita(5,3:m_length),'m',i,sita(6,3:m_length),'c',i,s ita(7,3:m_length),'k');