实验2---相关分析法辨识脉冲响应

实验2相关分析法辨识脉冲响应实验报告哈尔滨工业大学航天学院控制科学与工程系专业：自动化班级：1004102姓名：籍洋日期：2013年10月10日1．实验题目：相关分析法辨识系统脉冲响应2．实验目的通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。3．实验主要原理一个单入单出线性定常系统的动态特性可用它的脉冲响应函数g(σ)来描述。这样，只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数g(τ)。而在系统有正常输入的情形下，辨识脉冲响应的原理图如下图所示。0()()()ytgxtd则000()11lim()()(){lim()()}TTTTxtytxtdtgxtxtdtdTT上式两端同乘，进而取时间均值，有0()()()xyxRgRd则这就是著名的维纳霍夫积分方程。0()()(),()()()()()()()()xxxyxxyxtRkRkRgRdkgRgk如果输入是，这时的自相关函数为则根据维纳霍夫积分方程可得或者白噪声4．实验对象或参数下图为本实验的原理框图。系统的传递函数为)(sG,其中Sec26TSec,3812021..,TK；)()(kzku和分别为系统的输入和输出变量；)(kv为测量白噪声，服从正态分布，均值为零，方差为2v，记作),(~)(20vNkv；)(kg0为系统的脉冲响应理论值，)(ˆkg为系统脉冲响应估计值，)(~kg为系统脉冲响应估计误差。系统的输入采用M序列（采用实验1中的M序列即可），输出受到白噪声)(kv的污染。根据过程的输入和输出数据)(),(kzku，利用相关分析法计算出系统的脉冲响应值)(ˆkg，并与系统的脉冲响应理论值)(kg0比较，得到系统脉冲响应估计误差值)(~kg，当k时，应该有0)(~kg。相关分析法v(k)u(k)z(k)ckRtaNNkgMzPP)()()(ˆ2121210TtkTtkeeTTKkg//)()(ˆ)()(~kgkgkg0))(()(1121sTsTKsG1、模拟过程传递函数)(sG，获得过程的输入和输出数据)(),(kzku（采样时间取1秒）。(1)惯性环节其中，T为惯性环节的时间常数，K为惯性环节的静态放大倍数。若采样时间记作0T，则惯性环节的输出可写成：0011111000TkukuTeTTKkueTKkyekyTTTTTT)()())()()()()(///(2)传递函数)(sG仿真（串联）21211111TsTsTTKsG//)(令211TTKK＝，则)(sG的表达框图为：2、互相关函数的计算PPNrNiPMzizkiurNkR)()()()(111其中，r为周期数，1PNi表示计算互相关函数所用的数据是从第二个周期开始的，目的是等过程仿真数据进入平稳状态。（可分别令r=1、3，对比仿真结果）3、c的补偿补偿量c应取)(1PMzNR－，不能取)(PMzNR－。因为)(kRMz是周期函数，则有)()(0MzPMzRNR＝，故不能取)(PMzNR－。4、计算脉冲响应估计值111TsK/u(k)x(k)211Ts/y(k)TsK/1u(k)y(k)●脉冲响应估计值ckRtaNNkgMzPP)()()(ˆ21●脉冲响应估计误差PPNkNkgkgkgkg120120)()(ˆ)(＝5．程序框图TsK/1u(k)y(k)开始产生M序列和白噪声v求出系统在M序列作用下的输出y初始化参数T0,T1,T2,KZ=Y+V求Rmz求估计脉冲响应G求理论脉冲响应g0计算估计误差绘图结束6．程序代码function[sigma]=response(r)x=[0,1,0,1,1,0];%初始化Np=2^6-1;%M序列长度a=1;%振幅t=1;fori=1:Np*(r+1)y(i)=x(6);temp=xor(x(5),x(6));forj=5:-1:1x(j+1)=x(j);endx(1)=temp;endfori=1:Np*(r+1)if(y(i)==0)u(i)=a;elseu(i)=-a;endendK=120;T1=8.3;T2=6.2;T0=1;K1=K/T1/T2;x(1)=0;y(1)=0;fork=2:Np*(r+1)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)=exp(-T0/T2)*y(k-1)+T2*(1-exp(-T0/T2))*x(k-1)+T1*(T2*(exp(-T0/T2)-1)+T0)*(x(k)-x(k-1))/T0;%未经白噪声污染的输出endv=whitenoise(1,length(y));%产生白噪声z=y+v;%系统实际输出fork=1:Npsum=0;fori=Np+1:(r+1)*Npsum=u(i-k)*z(i)+sum;endRmz(k)=1/(r*Np)*sum;endc=-Rmz(Np-1);%补偿量c%计算脉冲响应估计值fork=1:NpG(k)=Np/((Np+1)*a^2*t)*(Rmz(k)+c);g0(k)=K/(T1-T2)*(exp(-k*t/T1)-exp(-k*t/T2));end%计算脉冲响应估计误差SUM1=0;SUM2=0;fork=1:Npe(k)=g0(k)-G(k);SUM1=e(k)^2+SUM1;SUM2=g0(k)^2+SUM2;endsigma=sqrt(SUM1/SUM2);step=0:Np-1;plot(step,[Rmz(63),Rmz(1:62)]);holdon;plot(step,[G(63),G(1:62)],'r');plot(step,[g0(63),g0(1:62)],'g');legend('互相关函数','脉冲响应估计值','脉冲响应理论值')end产生白噪声的函数：function[sig]=whitenoise(sigma,len)%白噪声产生函数，sigma为均方差，len为白噪声序列数据长度a=65539;M=2147483647;b=100;x(1)=12345;r(1)=x(1)/M;%第一部分为产生0-1的均匀分布随机数fori=1:12*lenx(i+1)=mod(a*x(i)+b,M);r(i+1)=x(i+1)/M;end%-----------------------------------------------------------------%n=12;fori=1:lensig(i)=0;forj=1:nsig(i)=sig(i)+r(n*(i-1)+j);%第二部分产生正态分布，方差为sigma的随机序列endendsig=(sig-12*0.5)*sigma;end7．实验结果及分析另r＝1,白噪声均方差sigma=0.5,运行命令：sigma=response(1)运行结果如下：得sigma=0.0416010203040506070-2-101234567互相关函数脉冲响应估计值脉冲响应理论值另r＝3,白噪声均方差sigma=0.5,运行命令：sigma=response(3,0.5)运行结果如下：得sigma=0.0373图像如下：010203040506070-2-101234567互相关函数脉冲响应估计值脉冲响应理论值比较发现r=1和r=3时产生的曲线基本相似,但是脉冲响应误差在r=3时更小。另r＝1,白噪声均方差sigma=1,运行命令：sigma=response(1,1)运行结果如下：得sigma=0.0594010203040506070-3-2-101234567互相关函数脉冲响应估计值脉冲响应理论值比较第一个与第三个结果,发现脉冲响应的估计误差是随着输入白噪声标准差的增大而增大的，白噪声标准差越小，对系统的输出干扰越小.8．结论在本次系统辨识的实验上机当中，在老师的指导之下，我利用相关分析法分析脉冲响应，得到r的值越大,得到的估计误差值越小;得到脉冲响应的估计误差是随着输入白噪声标准差的增大而增大的，带有白噪声污染的输出z，在白噪声标准差为0时与理想输出y是重合的，白噪声标准差越小，对系统的输出干扰越小.