最小二乘法+相关分析法

1、关于最小二乘法+相关分析法现在学习的是第一页，共26页最小二乘法1.问题描述考虑如下线性系统： 其中，u(k)为系统激励信号，y(k)为系统输出，e(k)为模型噪声。 1111abnanbz kaz ka z k nbu kb u k ne k 现在学习的是第二页，共26页系统模型现在学习的是第三页，共26页 其中G(z-1)是系统函数模型，N(z-1)为有色噪声系统模型，e(k)为白噪声v(k)经过系统函数为N(z-1)的系统后的输出。通常111111, B zD zG zN zA zC z11212112121aabbnnnnA za za za zB zb zb zb z 1121211

2、2121ccddnnnnC zc zc zc zD zd zd zb z 现在学习的是第四页，共26页则系统可表示为 设样本和参数集为 1111BzDzzkukvkAzCz121 2( ) -( -1 ) , - ( -2), . - ( - ), ( -1 ), ( -2), ., ( - ) , ,., , , ,., TTnnhkzkzkzk n ukukuk na aa bbbh(k)为可观测的量， 差分方程可写为最小二乘形式 ( )Tz khke k现在学习的是第五页，共26页 如何系统噪声e(k)存在的情况下从该方程中正确的解出，即是系统辨识的任务。 为了求出，我们面临三大问题：一

3、是输入信号的选择，二是判决准则的选取，三是辨识算法的选择 。选择输入 为了准确辨识系统参数，我们对输入信号有两大要求，一是信号要能持续的激励系统所有状态，二是信号频带能覆盖系统的频带宽度。除此之外还要求信号有可重复性，不能是不可重复的随机噪声，因此我们通常选择M序列或逆M序列作为输入。 现在学习的是第六页，共26页准则函数 本次探讨最小二乘类辨识方法，在此选取准则函数 使准则函数的估计值记做，称作参数的最小二乘估计值。 令k=1,2,3,L，可构成线性方程组： 2211TkkJe kz khk TLLLzkHkek现在学习的是第七页，共26页式中：准则函数相应变为 极小化,求得参数的估计值，将

4、使模型更好的预报系统的输出。 1122, 0101121211LLababLabzezezez Le LzznuunzznuunHz Lz Lnu Lu Ln 2211LLTTLLLLkkJe kz khkzHzH现在学习的是第八页，共26页最小二乘法 设 使 得，则有展开上式，并根据以下两个向量微分公式：得正则方程： LS minJ 0LSTLLLLJzHzH2 TTTTa xaxx Axx AAx为对称阵TTLLLSLLH HH z现在学习的是第九页，共26页当 为正则阵时，有且有 所以唯一 使得 ，这种通过极小化式计算的方法称作最小二乘法。 TLLH H1TTLSLLLLH HH z 2

5、220LSTLLJH HLS minJ现在学习的是第十页，共26页递推最小二乘法现在学习的是第十一页，共26页现在学习的是第十二页，共26页递推最小二乘法（RLS）步骤如下 现在学习的是第十三页，共26页举例 设某物理量Y与X 满足关系式Y=aX2+bX+c，实验获得一批数据如下表，试辨识模型参数a,b和c。 X1.012.033.024.0156.027.038.049.0310 Y9.64.11.30.40.050.10.71.83.89.0 分析：Y=aX2+bX+c+ =X(k)*X(k) X(k) 1*a b c+ = +现在学习的是第十四页，共26页程序X=1.01 2.03 3.

6、02 4.01 5 6.02 7.03 8.04 9.03 10;Y=9.6 4.1 1.3 0.4 0.05 0.10.71.83.89.0;%实验输入数据、实验输出数据syms a b c % 定义待辨识参数theta=a;b;c; %theta包含待辨识参数a，b，ctheta1=zeros(3,1); %对象参数初始化P=106*eye(3); %构造初始P阵for k=1:10 %仿真步长范围1到10phi=X(k)*X(k);X(k);1;%y=aX*X+bX+c=phi*theta%theta=a；b；c；phi=X(k)*X(k);X(k);1K=P*phi/(1+phi*P*

7、phi); %递推最小二乘法K阵的递推公式theta=theta1+K*(Y(k)-phi*theta1); %theta的递推公式P=(eye(3)-K*phi)*P; %递推最小二乘法P阵的递推公式theta1=theta; %theta的最终估计向量theta2(:,k)=theta; %theta估计向量矩阵化，目的是为了%下面的plot仿真图像输出endtheta1 %输出参数估计值plot(1:10,theta2) %输出参数逐步递推估计的轨迹图像xlabel(k); %设置横坐标为步长kylabel(参数估计a，b，c); %纵坐标为估计参数a，b，clegend(a,b,c);

8、 %标示相应曲线对应的参数axis(1 10 -10 20); %设置坐标轴范围现在学习的是第十五页，共26页实验结果 theta1 =0.4575 -5.0734 13.3711现在学习的是第十六页，共26页相关分析法 原理：一个单入单出线性定常系统的动态特性可用它的脉冲响应函数g()来描述。0 ( )( ) ()y tgx td则000()11lim( ) ()( )lim() () TTTTx ty t x tdtgx tx tdt dTT上式两端同乘，进而取时间均值，有现在学习的是第十七页，共26页令输入为白噪声，这时的自相关函数为：则根据维纳霍夫积分方程可得： 只要记录x(t)、y(

9、t)的值，并计算它们的互相关函数，即可求得脉冲响应函数。dtttgRRxxxy0)()()()(kRxx)(1)(Rxykg现在学习的是第十八页，共26页实例分析函数模型 其中K=120,T1=8.2Sec,T2=6.1Sec 输入信号采用6阶M序列，周期为63bit. 幅度为1，采样周期为1s, 利用相关分析法计算出系统的脉冲响应值, 并与系统的脉冲响应理论值. 实验步骤 1 ：得到仿真的输入输出数据M序列 2：利用相关分析法，求得系统的脉冲响应函数。 3：利用相关分析法求得的系统脉冲响应函数与实际脉冲响应函数。 4：对系统的输出加入噪声，画图比较利用相关分析法求得的系统脉冲响应函数与实际脉

10、冲响应函数。现在学习的是第十九页，共26页程序%创建M序列Np=63;%循环周期delta_T = 1;%时钟节拍a=1;%幅度%初始化M序列M(1)=1;M(2)=0;M(3)=0;M(4)=1;M(5)=1;M(6)=0;M_XuLie(Np) = 0;for n = 1 : Np temp = xor(M(6), M(5); if(temp = 0) M_XuLie(n) = a; else M_XuLie(n) = -a; end M(6) = M(5); M(5) = M(4); M(4) = M(3); M(3) = M(2); M(2) = M(1); M(1) = temp;e

11、nd%stairs(M_XuLie)%生成M序列完毕现在学习的是第二十页，共26页r = 3; % 周期数u=repmat(M_XuLie,1,r+1);%将M序列赋给输入，作为输入信号%第一步，从u(k)得到x(k),y(k)K = 120;T0 = 1;% 采样时间T1 = 8.1;T2 = 6.2;K1=K/(T1*T2);%初始化X(k),Y(k)为0K2=1x(63)=0;y(63)=0for k = 2 : 63*4 %取得x(k)序列 x(k)=exp(-T0/T1)*x(k-1)+T1*K1*(1-exp(-T0/T1)*u(k-1)+T1*K1. *(T1*(exp(-T0/

12、T1)-1)+T0)*(u(k)-u(k-1)/T0 %取得y(k)序列 y(k)=exp(-T0/T2)*y(k-1)+T2*K2*(1-exp(-T0/T2)*x(k-1)+T2*K2. *(T2*(exp(-T0/T2)-1)+T0)*(x(k)-x(k-1)/T0end%获取没有白噪声时候输出完毕%作图figure(1);plot(u,r);hold on;plot(x,k);plot(y,b);legend(u(k),x(k),y(k);现在学习的是第二十一页，共26页%第二步，将白噪声添加入输出信号%产生白噪声信号vfangcha = 0.5;%随意指定的方差v = fangcha

13、 * randn(1,63*4);%信号叠加,输出实际信号z(k)z = y + v;%figure(2);%plot(v,m);%hold on;%plot(z,y);figure(2);%打印无白噪声污染信号plot(y,b);hold on;%打印白噪声信号plot(v,m);%打印白噪声污染后的信号plot(z,g);legend(y(k),v(k),z(k);%计算Rmz(k)for k = 1 : Np Rmz(k)=0;%初始化为0 for i = (Np + 1) : (r+1)*Np) Rmz(k)=Rmz(k) + u(i-k)*z(i); end Rmz(k)=Rmz(k)/(r*Np);end现在学习的是第二十二页，共26页%计算cc=-Rmz(Np - 1);%计算脉冲响应估计值g1g1=Np*(Rmz+c)/(Np+1)*a2*delta_T);%计算理论脉冲g0for k = 1: Np g0(k)=K/(T1-T2)*(exp(-k*delta_T/T1)-exp(-k*delta_T/T2);end%计算脉冲响