最小二乘法+相关分析法

1、关于最小二乘法+相关分析法第一张，PPT共二十六页，创作于2022年6月最小二乘法1.问题描述考虑如下线性系统： 其中，u(k)为系统激励信号，y(k)为系统输出，e(k)为模型噪声。第二张，PPT共二十六页，创作于2022年6月系统模型第三张，PPT共二十六页，创作于2022年6月 其中G(z-1)是系统函数模型，N(z-1)为有色噪声系统模型，e(k)为白噪声v(k)经过系统函数为N(z-1)的系统后的输出。通常第四张，PPT共二十六页，创作于2022年6月则系统可表示为 设样本和参数集为 h(k)为可观测的量， 差分方程可写为最小二乘形式第五张，PPT共二十六页，创作于2022年6月 如

2、何系统噪声e(k)存在的情况下从该方程中正确的解出，即是系统辨识的任务。 为了求出，我们面临三大问题：一是输入信号的选择，二是判决准则的选取，三是辨识算法的选择 。选择输入 为了准确辨识系统参数，我们对输入信号有两大要求，一是信号要能持续的激励系统所有状态，二是信号频带能覆盖系统的频带宽度。除此之外还要求信号有可重复性，不能是不可重复的随机噪声，因此我们通常选择M序列或逆M序列作为输入。 第六张，PPT共二十六页，创作于2022年6月准则函数 本次探讨最小二乘类辨识方法，在此选取准则函数 使准则函数的估计值记做，称作参数的最小二乘估计值。 令k=1,2,3,L，可构成线性方程组： 第七张，PP

3、T共二十六页，创作于2022年6月式中：准则函数相应变为 极小化,求得参数的估计值，将使模型更好的预报系统的输出。第八张，PPT共二十六页，创作于2022年6月最小二乘法 设 使 得，则有展开上式，并根据以下两个向量微分公式：得正则方程： 第九张，PPT共二十六页，创作于2022年6月当 为正则阵时，有且有 所以唯一 使得 ，这种通过极小化式计算的方法称作最小二乘法。 第十张，PPT共二十六页，创作于2022年6月递推最小二乘法第十一张，PPT共二十六页，创作于2022年6月第十二张，PPT共二十六页，创作于2022年6月递推最小二乘法（RLS）步骤如下 第十三张，PPT共二十六页，创作于20

4、22年6月举例设某物理量Y与X 满足关系式Y=aX2+bX+c，实验获得一批数据如下表，试辨识模型参数a,b和c。X1.012.033.024.0156.027.038.049.0310Y9.64.11.30.40.050.10.71.83.89.0分析：Y=aX2+bX+c+ =X(k)*X(k) X(k) 1*a b c+ = +第十四张，PPT共二十六页，创作于2022年6月程序X=1.012.03 3.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

5、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*phi); %递推最小二乘法K阵的递推公式theta=theta1+K*(Y(k)-phi*theta1); %theta的递推公式P=(eye(3)-K*phi)*P; %递推最小二乘

6、法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); %标示相应曲线对应的参数axis(1 10 -10 20); %设置坐标轴范围第十五张，PPT共二十六页，创作于2022年6月实验结果theta1 =0.4575 -5.0734 13

7、.3711第十六张，PPT共二十六页，创作于2022年6月相关分析法原理：一个单入单出线性定常系统的动态特性可用它的脉冲响应函数g()来描述。第十七张，PPT共二十六页，创作于2022年6月令输入为白噪声，这时的自相关函数为：则根据维纳霍夫积分方程可得： 只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数。第十八张，PPT共二十六页，创作于2022年6月实例分析函数模型 其中K=120,T1=8.2Sec,T2=6.1Sec 输入信号采用6阶M序列，周期为63bit. 幅度为1，采样周期为1s, 利用相关分析法计算出系统的脉冲响应值, 并与系统的脉冲响应理论值. 实验

8、步骤 1 ：得到仿真的输入输出数据M序列 2：利用相关分析法，求得系统的脉冲响应函数。 3：利用相关分析法求得的系统脉冲响应函数与实际脉冲响应函数。 4：对系统的输出加入噪声，画图比较利用相关分析法求得的系统脉冲响应函数与实际脉冲响应函数。第十九张，PPT共二十六页，创作于2022年6月程序%创建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 =

9、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;end%stairs(M_XuLie)%生成M序列完毕第二十张，PPT共二十六页，创作于2022年6月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);%初

10、始化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/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);pl

11、ot(u,r);hold on;plot(x,k);plot(y,b);legend(u(k),x(k),y(k);第二十一张，PPT共二十六页，创作于2022年6月%第二步，将白噪声添加入输出信号%产生白噪声信号vfangcha = 0.5;%随意指定的方差v = fangcha * 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(

12、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第二十二张，PPT共二十六页，创作于2022年6月%计算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%计算脉冲响应估计误差delta_gdelta_g=sqrt(sum(g0-g1).2)/sum(g0.2);fi