最小二乘法+相关分析法

关于最小二乘法+相关分析法第一页，共二十六页，编辑于2023年，星期一最小二乘法1.问题描述考虑如下线性系统：

其中，u(k)为系统激励信号，y(k)为系统输出，e(k)为模型噪声。第二页，共二十六页，编辑于2023年，星期一系统模型第三页，共二十六页，编辑于2023年，星期一

其中G(z-1)是系统函数模型，N(z-1)为有色噪声系统模型，e(k)为白噪声v(k)经过系统函数为N(z-1)的系统后的输出。通常第四页，共二十六页，编辑于2023年，星期一则系统可表示为

设样本和参数集为h(k)为可观测的量，差分方程可写为最小二乘形式第五页，共二十六页，编辑于2023年，星期一

如何系统噪声e(k)存在的情况下从该方程中正确的解出，即是系统辨识的任务。为了求出，我们面临三大问题：一是输入信号的选择，二是判决准则的选取，三是辨识算法的选择。选择输入

为了准确辨识系统参数，我们对输入信号有两大要求，一是信号要能持续的激励系统所有状态，二是信号频带能覆盖系统的频带宽度。除此之外还要求信号有可重复性，不能是不可重复的随机噪声，因此我们通常选择M序列或逆M序列作为输入。第六页，共二十六页，编辑于2023年，星期一准则函数

本次探讨最小二乘类辨识方法，在此选取准则函数使准则函数的估计值记做，称作参数的最小二乘估计值。令k=1,2,3,……L，可构成线性方程组：第七页，共二十六页，编辑于2023年，星期一式中：准则函数相应变为极小化,求得参数的估计值，将使模型更好的预报系统的输出。第八页，共二十六页，编辑于2023年，星期一最小二乘法设使得，则有展开上式，并根据以下两个向量微分公式：得正则方程：第九页，共二十六页，编辑于2023年，星期一当为正则阵时，有且有所以唯一使得，这种通过极小化式计算的方法称作最小二乘法。

第十页，共二十六页，编辑于2023年，星期一递推最小二乘法第十一页，共二十六页，编辑于2023年，星期一第十二页，共二十六页，编辑于2023年，星期一递推最小二乘法（RLS）步骤如下第十三页，共二十六页，编辑于2023年，星期一举例设某物理量Y与X满足关系式Y=aX2+bX+c，实验获得一批数据如下表，试辨识模型参数a,b和c。X 1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10Y 9.6 4.1 1.3 0.4 0.05 0.1 0.7 1.8 3.8 9.0分析：Y=aX2+bX+c+=[X(k)*X(k)X(k)1]'*[abc]+=+第十四页，共二十六页，编辑于2023年，星期一程序X=[1.01 2.033.024.0156.027.038.049.0310];Y=[9.64.11.30.40.050.1 0.7 1.8 3.8 9.0];%实验输入数据、实验输出数据symsabc%定义待辨识参数theta=[a;b;c];%theta包含待辨识参数a，b，ctheta1=zeros(3,1);%对象参数初始化P=10^6*eye(3);%构造初始P阵fork=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);1]K=P*phi/(1+phi'*P*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');%标示相应曲线对应的参数axis([110-1020]);%设置坐标轴范围第十五页，共二十六页，编辑于2023年，星期一实验结果theta1=0.4575-5.073413.3711第十六页，共二十六页，编辑于2023年，星期一相关分析法原理：一个单入单出线性定常系统的动态特性可用它的脉冲响应函数g(σ)来描述。第十七页，共二十六页，编辑于2023年，星期一令输入为白噪声，这时的自相关函数为：则根据维纳—霍夫积分方程可得：只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数。第十八页，共二十六页，编辑于2023年，星期一实例分析函数模型

其中K=120,T1=8.2Sec,T2=6.1Sec

输入信号采用6阶M序列，周期为63bit.幅度为1，采样周期为1s,利用相关分析法计算出系统的脉冲响应值,并与系统的脉冲响应理论值.实验步骤1：得到仿真的输入输出数据M序列2：利用相关分析法，求得系统的脉冲响应函数。3：利用相关分析法求得的系统脉冲响应函数与实际脉冲响应函数。4：对系统的输出加入噪声，画图比较利用相关分析法求得的系统脉冲响应函数与实际脉冲响应函数。第十九页，共二十六页，编辑于2023年，星期一程序%创建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;forn=1:Nptemp=xor(M(6),M(5));if(temp==0)M_XuLie(n)=a;elseM_XuLie(n)=-a;endM(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序列完毕第二十页，共二十六页，编辑于2023年，星期一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)=0fork=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);plot(u,'r');holdon;plot(x,'k');plot(y,'b');legend('u(k)','x(k)','y(k)');第二十一页，共二十六页，编辑于2023年，星期一%第二步，将白噪声添加入输出信号%产生白噪声信号vfangcha=0.5;%随意指定的方差v=fangcha*randn(1,63*4);%信号叠加,输出实际信号z(k)z=y+v;%figure(2);%plot(v,'m');%holdon;%plot(z,'y');figure(2);%打印无白噪声污染信号plot(y,'b');holdon;%打印白噪声信号plot(v,'m');%打印白噪声污染后的信号plot(z,'g');legend('y(k)','v(k)','z(k)');%计算Rmz(k)fork=1:NpRmz(k)=0;%初始化为0fori=(Np+1):((r+1)*Np)Rmz(k)=Rmz(k)+u(i-k)*z(i);endRmz(k)=Rmz(k)/(r*Np); end第二十二页，共二十六页，编辑于2023年，星期一%计算cc=-Rmz(Np-1);%计算脉冲响应估计值g1g1=Np*(Rmz+c)/((Np+1)*a^2*delta_T);%计算理论脉冲g0fork=1:Npg0(k)=K/(T1-T2)*(exp(-k*delta_T/T1)-exp(-k*delta_T/T2));end%计算脉冲响应估计误差delta_gdelta_g=sqrt(sum((