过程辨识仿真作业.doc

例6.6仿真对象如图所示：图中，为服从分布的不相关随机噪声；控制值，是数据的信噪比为23%。经过噪声模型后，迭加在输出y(k)上； 模型结构选用如下形式：输入信号u(k)采用4阶的M序列，特征多项式取F(s)=1s，幅度为1, 循环周期为31bit，数据长度L=480, 初始条件=0.001, P(0)=I。利用辅助变量法辨识系统参数、。一、 实验原理：设SISO过程采用如下数学模型 描述，其中和表示过程的输入输出；是零均值的有色噪声；且 并且假定模型的阶次已知。由于是有色噪声，直接利用最小二乘法不能获得模型参数的无偏一致估计，这时可以利用辅助变量法，一起获得无偏一致估计。1) 一次完成算法 将模型化为最小二乘格式 或 根据最小二乘知识，参数的最小估计值为其中 由Frechet定理知如果是白噪声，则，那么如果不是白噪声，一般。为了使不是白噪声得情况下，仍有，定义一个辅助变量矩阵 = 使之满足以下两个条件：是非奇异阵； ，即与独立。其中，称作辅助变量。如果所用辅助变量能满足上述两个条件，则有 (*)其中称作辅助变量参数估计值。可见只要适当选择辅助变量，使之满足以上两个条件，参数就可以是无偏一致估计。2) 辅助变量的选择选择辅助变量的基本原则就是上述两个条件必须满足。如果采用图a所示的作为辅助变量，置图a 辅助变量的选择当是持续激励信号时，必有是非奇异的，且因只与有关，即必与噪声无关；故有。这样辅助变量常有以下几种选择方法：1 自适应滤波法 把图a中辅助模型看成自适应滤波器时，辅助变量可按如下关系或 (*) 取0.010.1；取010 求得。2 纯滞后法当图a中的辅助模型为纯滞后环节时，辅助变量取作其中，是多项式的阶次。此时，可写成当是持续激励信号，且与噪声无关，则辅助变量可满足上述两个条件。3 Tally原理如果噪声可看成下列模型的输出其中，是零均值的不相关噪声；且则辅助变量取作即可满足上述两个条件，显然，只要输入信号是持续激励信号，条件即可满足。另外，由于输入信号与噪声无关，故有且那么则条件也满足。3) 递推算法将(*)式一次完成算法写成定义则和推导最小二乘递推算法一样，辅助变量法的递推算式(简称RIV)可以写成式中，为辅助向量，如果辅助变量选取(*)式，则需要用最小二乘法先递推若干步，以获取初步参数估计值，作为辅助变量法的递推初始状态。二、 实验步骤1） 输入信号的产生输入信号采用4阶M序列，特征多项式取，幅度为1，循环周期为。2） 噪声的产生图b 仿真系统 如上图所示，由经过噪声模型后产生。要求为服从分布的不相关随机噪声；控制值，是数据的信噪比为23%。3） 辅助变量的选择使用Tally原理给出辅助变量，辅助变量取做，则辅助变量可写成：其中，=2；=2；即4） 利用辅助变量法进行参数估计(1)初始化系统辨识初值：根据题目要求：=0.001, P(0)=I(2)由于辅助最小二乘算法对初始值P(0)敏感，因此采用最小二乘算法（LS）进行辨识算法的起步，即用最小二乘辨识方法辨识前100步，然后用辅助最小二乘算法辨识。否则可能辨识没有可靠的收敛性；(3)计算残差(k)的统计性质；(5)计算阶跃响应。三、 实验结果 根据实验要求产生的4阶M序列如下图所示：图c 幅度为1的4阶M序列图d RIV参数估计值的变化过程图e ls参数估计值的变化过程图f RIV阶跃响应比较表1 RIV算法的辨识结果（噪信比=23%）参 数a1a2b1b2静态增益噪声均值噪声方差真 值-1.50.71.00.57.50.00.23估计值-1.51430.71000.90640.49757.53322.24E-20.2738表2 阶跃响应比较k模型阶跃响应过程阶跃响应k模型阶跃响应过程阶跃响应012345678910111213141500.0553 0.0543 1.2950 3.1685 5.4466 7.43508.9429 9.2900 9.4649 8.9030 8.5585 7.77417.0527 6.5679 6.4264 00.0553 0.0543 0.7868 2.6600 4.6122 6.72888.1230 8.9962 8.4569 8.4759 7.5047 7.38106.4403 5.9047 5.68191617181920212223242526272829306.4046 6.4181 6.54136.6312 6.5494 6.5866 7.3975 7.7134 7.72207.4495 7.3840 6.8860 6.6050 6.7591 6.26735.8103 5.8771 5.91296.0892 6.1381 5.9512 6.0651 7.2625 7.16746.9576 6.5403 6.6336 5.9283 5.8554 6.2861表3 输出残差序列(k)的统计性质均值(k)=-2.29E-2(0)1.0(1) 0.2621 (6)-0.1803(11)0.3579(16)0.3454(2)-0.4282 (7)-0.6902(12)0.3255(17)-0.4567(3)-0.5079 (8)-0.7099(13)0.7054(18)-0.4462(4)-0.3423 (9)-0.9590(14)1.2514(19)-0.3574(5)-0.1340 (10)-0.3243(15)0(20)-0.0390四、 实验结果分析由于本题的输入是有色噪声，因此用最小二乘辨识方法一般得不到无偏一致估计。这里采用辅助变量法，递推计算480步后的辨识结果如表1所示。参数估计的变化过程如图d所示，模型与过程的阶跃响应比较吻合，如图f与表2所示，辨识结果是令人满意的。为进一步确认辨识结果，需要对所获得的模型进行检验。计算输出残差(k)的均值和自相关系数，结果如表3所示。表3表明，由于数值比较大，可以确信输出残差(k)序列是有色噪声，而输入v(k)就是有色噪声，因此所获得的模型是可靠的。仿真实验表明，辅助变量法的计算量与最小二乘相当，但辨识结果却比最小二乘法好的多。因此辅助变量法是一种很有价值的辨识方法，尤其当噪声是有色的，而噪声的模型又不好确定时，辅助变量法就显示出了它的优势。但辅助变量法不能像增广最小二乘法或广义最小二乘法那样可以同时获得噪声模型的参数估计。附录 实验源程序n=4;%生成4阶M序列a=zeros(2n-1,n);a(1,n)=1;for i=2:1:2n-1a(i,1)=mod(a(i-1,n-1)+a(i-1,n),2);for p=2:1:na(i,p)=a(i-1,p-1);endendr=a(:,1);r1=r;r;r;r;r;r;r;r;r2=r1;r1;r1;r1;u=zeros(1,480);for i=1:480 u(i)=1*(1-2*r2(i);% 生成4阶M序列 endfigure(1)stairs(u);axis(1 480 -1.5 1.5)grid on;p=106*eye(4);theta=0.001 0.001 0.001 0.001;I=eye(4);a1=zeros(1,481);a2=a1;b1=a1;b2=a1;z=a1;y=a1;rou=zeros(1,20);%辅助变量%f=0 0 0 0;A=1 -1.5 0.7;B=0 1.0 0.5;C=1 -1 0.2;E=zeros(480,1);T=zeros(480,4);T(1,:)=theta;v=normrnd(0,1,480,1);%生成480均值，方差为1的正态分布噪声V=0.23*v;e=filter(C,A,V);%生成有色噪声z(1)=v(1);z(2)=v(2);%计算输出for k=3:1:480 z(k)=1.5*z(k-1)-0.7*z(k-2)+1.0*u(k-1)+0.5*u(k-2)+e(k);end%辅助变量法进行辨识for i=2:1:40 theta=theta+p*h*(h*p*h+1)(-1)*(z(i)-h*theta); p=(I-p*h*(h*p*h+1)(-1)*h)*p; h=-z(i) -z(i-1) u(i) u(i-1); T(i,:)=theta; E(i)=z(i)-h*T(i,:);endfor i=41:1:480 theta=theta+p*f*(h*p*f+1)(-1)*(z(i)-h*theta); p=(I-p*f*(h*p*f+1)(-1)*h)*p; h=-z(i) -z(i-1) u(i) u(i-1); f=-z(i-3) -z(i-4) u(i-1) u(i-2); %Tally原理 T(i,:)=theta; E(i)=z(i)-h*T(i,:);endk=1:1:480;figure(2);plot(k,T(:,1),b);axis(1 540 -1.6 1.6)hold on;grid on;k=1:1:480;plot(k,T(:,2),r);hold on;grid on;k=1:1:480;plot(k,T(:,3),g);hold on;grid on;k=1:1:480;plot(k,T(:,4),k);hold on;grid on;title(RIV);text(490,1.0,leftarrow b1);text(490,0.7,leftarrow a2);text(490,0.5,leftarrow b2);text(490,-1.5,leftarrow a1);xlabel(k);ylabel(theta);%输出残差的计算a=mean(E)for i=1:1:480 A0=E(i)*E(i);endR0=mean(A0);for i=1:1:479 A1=E(i)*E(i+1);endR(1)=mean(A1);for i=1:1:478 A2=E(i)*E(i+2);endR(2)=mean(A2);for i=1:1:477 A3=E(i)*E(i+3);endR(3)=mean(A3);for i=1:1:476 A4=E(i)*E(i+4);endR(4)=mean(A4);for i=1:1:475 A5=E(i)*E(i+5);endR(5)=mean(A5);for i=1:1:474 A6=E(i)*E(i+6);endR(6)=mean(A6);for i=1:1:473 A7=E(i)*E(i+7);endR(7)=mean(A7);for i=1:1:472 A8=E(i)*E(i+8);endR(8)=mean(A8);for i=1:1:471 A9=E(i)*E(i+9);endR(9)=mean(A9);for i=1:1:470 A10=E(i)*E(i+10);endR(10)=mean(A10);for i=1:1:469 A11=E(i)*E(i+11);endR(11)=mean(A11);for i=1:1:468 A12=E(i)*E(i+12);endR(12)=mean(A12);for i=1:1:467 A13=E(i)*E(i+13);endR(13)=mean(A13);for i=1:1:466 A14=E(i)*E(i+14);endR(14)=mean(A14);for i=1:1:465 A15=E(i)*E(i+15);endR(1)=mean(A1);for i=1:1:464 A16=E(i)*E(i+16);endR(16)=mean(A16);for i=1:1:463 A17=E(i)*E(i+17);endR(17)=mean(A17);for i=1:1:462 A18=E(i)*E(i+18);endR(18)=mean(A18);for i=1:1:461 A19=E(i)*E(i+19);endR(19)=mean(A19);for i=1:1:460 A20=E(i)*E(i+20);endR(20)=mean(A20);for i=1:1:20 rou(i)=R(i)/R0;endrou%阶跃响应Gain1=;Gain=;B=1.0 -1 0.2;A=1 -1.5 0.7;u=ones(1,30);v=normrnd(0,1,30,1); v=0.23*v;e=filter(A,B,v); junzhi=mean(e)%求噪声均值fangcha=std(e)%求噪声方差y=zeros(1,30);z=zeros(1,30);y(1)=v(1);z(1)=v(1);y(2)=v(2);z(2)=v(2);for i=3:1:30 y(i)=1.5*y(i-1)-0.7*y(i