过程辨识仿真作业

1、例6 6仿真对象如图所示:Xv(k)l-zT+0.2z-21 - l.5z-i + 0.7z-2e(k) +/I u(k)l.Oz-1 + 05z-2y(k) + /1 - l.5zT + 0.7z-2z(k)图中，v(k)为服从N(OJ)分布的不相关随机噪声：控制X值，是数据的信噪比为23%o v(k) 经过噪声模型后，迭加在输出y(k)上；模型结构选用如下形式：z(k) + a1z(k -1)+ a2z(k -2) = b1n(k-1) + b2u(k -2) + e(k)输入信号u(k)采用4阶的M序列，特征多项式取F(s)二1sS4,幅度为1,循 环周期为31bit,数据长度L=480

2、,初始条件3(0) =0.001, P(0)=io6Io利用辅助 变量法辨识系统参数a、a2> “、b2o一、实验原理：设SISO过程采用如下数学模型&zT)z(k) = B(zT)ii(k) + e(k)描述，其中u(k)和z(k)表示过程的输入输出：e(k)是零均值的有色噪声；且f = l + a-1 +a2z_2 Hba/ z-lla1 B(z-1) = 1+tz-14-b:z-2 hnzf并且假定模型的阶次1、弭已知。由于e(k)是有色噪声，直接利用最小二乘法不能获得模型参数的无偏一致估计，这时可以利用辅助变量法，一起获得 无偏一致估计。1) 一次完成算法将模型化为最小二

3、乘格式z(k) = h。(k)e + e(k)或Zi = Hl +eL根据最小二乘知识，参数的最小估计值为A-11Tl抵=(H；Hl) HlZl = 0o+(HlHl) (Hlel)其中r II L.y7HLHl=7 Sh(k)hr(k).>E(h(k)l(k)LLk=l41 L1H£eL=- Xh(k)e(k) >Eh(k)e(k)LLk=l由Frechet定理知Ou 士 J = a+Eh(k)h7(k) 'Eh(k)e(k)如果e(k)是白噪声，则Eh(k)e(k)=0,那么如果e(k)不是白噪声，一般Eh(k)e(k)wO,为了使e(k)不是白噪声得情况下，

4、仍有北 "了仇，定义一个辅助变量矩阵 常、H；= h"(2)使之满足以下两个条件："H T H L L-s >Eh*(k)h(k) 是非奇异阵；I-/(1一H；%l>E(h*(k)e(k) = 0,即 h*(k)与 e(k)独立。其 L/中，K(k)称作辅助变量。如果所用辅助变量能满足上述两个条件，则有A* T1I = (HLrHL)HLZl = o(-HjJHl) (-HLreL)J/J/L-s > = &+Eh"(k)hT(k) Eh*(kXk)=仇 A其中。，称作辅助变量参数估计值。可见只要适当选择辅助变量，使之满A足以上

5、两个条件，参数就可以是无偏一致估计。2)辅助变量的选择选择辅助变量的基本原则就是上述两个条件必须满足。如果采用图a所示 的x(k)作为辅助变量，置图a辅助变量的选择h* (k) = -x(k -1),-x(k - nJ ,u(k -1),u(k -现)了当u(k)是持续激励信号时，必有Eh*(k)h工(k)是非奇异的，且因x(k) 只与u(k)有关，即h*(k)必与噪声无关；故有Eh*(k)e(k)二°。这样辅 助变量常有以下几种选择方法：1 .自适应滤波法把图a中辅助模型看成自适应滤波器时，辅助变量可按如下关系x(k)= i/(k)6(k)或C x(k) = h*r(k)(k)(*

6、)< O(k) = (l a)£(k-1) +闻k-d)a 取 0.0广0.1： d 取 010求得。2 .纯滞后法当图a中的辅助模型为纯滞后环节时，辅助变量取作x(k) = u(k-nb)其中，1%是多项式BQ-1)的阶次。此时，h'(k)可写成h*(k) = -u(k 取 - 1),-u(k- 1% 一%),u(k- l),-su(k-nb)c 当u(k)是持续激励信号，且与噪声无关，则辅助变量可满足上述两个条件。3 . Tally 原理如果噪声e(k)可看成下列模型的输出e(k) = D(z-)v(k)其中，v(k)是零均值的不相关噪声；且D(z-1 = 1 +

7、djZ-1 +d2z-2 + dlldZ-nd)则辅助变量取作x(k) = z(k-iid)即h* (k) = -z(k 4 1),z(k - 0 q ) ,u(k 1),u(k & ) 可满足上述两个条件，显然，只要输入信号u(k)是持续激励信号，条件即可 满足。另外，由于输入信号u(k)与噪声e(k)无关，故有Eu(k - i)e(k) = 0, i = 1,2,i1b且Ez(k-nd-i)e(k)= Ez(k-nd-i)e(k)D(z-1)v(k)=Ez(k- nd - i)v(k)+dEz(k-% -i)v(k- 1)+Ez(k- nd -i)v(k- nd) = 0其中,i

8、= 1,2,i1b那么Eh*(k)e(k)= 0则条件也满足。3)递推算法将(*)式一次完成算法写成Lt L北=(H f H 1y H 丁Zl = (Z h'h"(i) (E h*(i)z(i)i=li=l定义%T(k)= Sx(k)hYi)K(k) = P(k)h* (k)则和推导最小二乘递推算法一样，辅助变量法的递推算式(简称RIV)可以 写成6(k) = 0(k-l)+k(kX z(k)-hr(k)0(k -1) k(k)=P(k-l)h* (k)hr (k)P(k-l)h* (k-l)+l-1P(k)=I-K(k)hr (k)P(k-l)式中，h*(k)为辅助向量，如

9、果辅助变量x(k)选取(*吟式，则需要用最A小二乘法先递推若干步，以获取初步参数估计值。(k),作为辅助变量法的 递推初始状态。二、实验步骤1)输入信号U(k)的产生输入信号u(k)采用4阶M序列，特征多项式取F(s) = lss,幅 度为L循环周期为Np=31bit。2)噪声e(k)的产生如上图所示，e(k)由v(k)经过噪声模型后产生。要求v(k)为服从 N(M)分布的不相关随机噪声；控制人值，是数据的信噪比为23%。3)辅助变量的选择使用Tally原理给出辅助变量，辅助变量取做处)=2比一%),则 辅助变量可写成：h，(k) = -z(k - 为 1),_z(k ) ,u(k 1),u(

10、k 其中，11尸2； %=2；即h* (k) = -z(k - 3)-z(k - 4)9u(k -l)9u(k - 2)f4)利用辅助变量法进行参数估计初始化系统辨识初值：根据题目要求：3(0)=0.001, P(O)=1O6I(2)由于辅助最小：乘算法对初始值P(0)敏感，因此采用最小二乘算法(LS)进行辨识算法的起步,即用最小二乘辨识方法辨识前100步,然后用辅助最小二乘算法辨识。否则可能辨识没有可靠的收敛性；(3)计算残差£ (k)的统计性质;(5)计算阶跃响应。三、实验结果根据实验要求产生的4阶M序列如下图所示:图c幅度为1的4阶M序列图d RIV参数估计值的变化过程Is的辨

11、识结果_ h.1 1n s4 Ujr|1 J nj 一一 vx *I.一 jr<- DI-a2 b2n .v*->*Tu-0.5 -1 -1.5-V at50100 150 200 250 300 350 400 450 5001.5k图e Is参数估计值的变化过程图fRIV阶跃响应比较U!实验结果分析表1 RIV算法的辨识结果（噪信比23%）参数3132bib2静态增益噪声均值噪声方差真值-1.50.71.00.57.50.00.23估计值-1.51430.71000.90640.49757.53322.24E-20.2738表2阶跃响应比较k模型阶跃响应过程阶跃响应k模型阶跃响

12、应过程阶跃响应000166.40465.810310.05530.0553176.41815.877120.05430.0543186.54135.912931.29500.7868196.63126.089243.16852.6600206.54946.138155.44664.6122216.58665.951267.43506.7288227.39756.065178.94298.1230237.71347.262589.29008.9962247.72207.167499.46498.4569257.44956.9576108.90308.4759267.38406.5403118.5

13、5857.5047276.88606.6336127.774173810286.60505.9283137.05276.4403296.75915.8554146.56795.9047306.26736.2861156.42645.6819表3输出残差序列£&）的统计性质均值£（k）=229E-2P（0）1.0p（l）0 2621P（6）-0.1803P（n）0.3579p（16）03454P（2）04282P（7）-0.6902P（12）0.3255P（"）-0.4567P（3）-0.5079P（8）-0.7099p（13）0.7054p（18）-0.4

14、462P（4）-0.3423P（9）-0.9590p（14）1.2514p（19）-0.3574P（5）-0 1340PQO）-0.3243pQ5）0PQO）-0.0390由于本题的输入是有色噪声，因此用最小二乘辨识方法一般得不到无偏一致 估计。这里采用辅助变量法，递推计算480步后的辨识结果如表1所示。参数估 计的变化过程如图d所示，模型与过程的阶跃响应比较吻合，如图f与表2所示, 辨识结果是令人满意的。为进一步确认辨识结果，需要对所获得的模型进行检验。 计算输出残差£ （k）的均值和自相关系数，结果如表3所示。表3表明，由于 数值比较大，可以确信输出残差七（k）序列是有色噪声，而

15、输入v（k）就是有色 噪声，因此所获得的模型是可靠的。仿真实验表明，辅助变量法的计算量与最小 二乘相当，但辨识结果却比最小二乘法好的多.因此辅助变量法是一种很有价值 的辨识方法，尤其当噪声是有色的，而噪声的模型乂不好确定时，辅助变量法就 显示出了它的优势。但辅助变量法不能像增广最小二乘法或广义最小二乘法那样可以同时获得噪声模型的参数估计。附录实验源程序n=4;*牛成4阶M字歹lja=zeros (2人nT, n);a (1, n) =1;for i=2:l:2"n-la (i, 1) =mod (a (i-l, n-l) +a (i-l, n) , 2);for p=2:1:na (

16、i, p) =a (i-1, p-1)；endendr=a(: r 1);rl=r; r;r;r;r;r;r;r;r2=rl;rl;rl;rl;u=zeros(1,480);for i=l:480u (i) =1* (l-2*r2 (i) ;% 生成 4 阶M 序列endfigure(1)stairs(u);axis(1 480 -1.5 1.5)grid on;p=10人6*eye(4);theta=0.001 0.001 0. 001 0. 001 1 ;I=eye (4);al=zeros (lf 481);a2=al;bl=al;b2=al;z=al;y=al;rou=zeros (1

17、, 20);f=0 0 0 0 f ;A=L -1.5 0.7;B=0 1.0 0.S;C=l -1 0.2;E=zeros(480, 1);T=zeros(480, 4);T (1, :) =thetaf;v=normrnd (0, 1, 480, I) ; *，上成48。均值，方不为 1的R i v=0.Z3*v;e=filter (C, A, V) ; %上成有色噪声z(l)=v(l);z(2)=v (2);%计芽输出for k=3:1:480z (k)=1.5*z (k-l)-0.7*z(k-2)+l. 0*u (k-1) +0. 5*u (k-2) +e (k); end务辅助变量法

18、进行辨识 for i=2:1:40theta=theta+p*h* ( (hf *p*h+l) A (-1) ) (z (i)-hf *theta);p=(I-p*h* (hf *p*h+l) (-1) *hf) *p; h=-z(i) -z(i-l) u(i) u(i-l) 1 ; T (i, : ) =thetaf ;E(i)=z(i)-hf*T(ir : ) f ;endfor i=41: 1: 480theta=theta+p*f * ( (hf *p*f+l)人(-1) ) * (z (i) -hf *theta);p=(hf *p*f+l) /s(-l)*hf )*p;h=-z (

19、i) -z (i-1) u(i) u(i-l) 1 ;f=-z (i-3) -z(i-4) u(i-l) u(i-2) f ; *Tally原理 T (i, : ) =thetaf ;E (i)=z (i)-hf *T (i, : ) f ; endk=l:1: 480;figure (2);plot (k,T(:, 1), »bT);axis (1 540 -1. 6 1.6)hold on;grid on;k=l:1:480;plot(k,T(:,二)，)；hold on;grid on;k=l:1:480;plot(k,T (:, 3), 'g)；hold on;gri

20、d on;k=l:1:480;plot (k,T(:, 4),;hold on;grid on;title (f RIVpA±«E3*ia xu1);text (490/ 1. 0/ 1 Meftarrovz blf);text(490r 0.7r 1leftarrow aZf);text(490r 0.5r 1leftarrow b2 f);text(490,-1.5, fieftarrow alf);xlabel(fkf);ylabel('theta');8输出残差的计算a=mean(E)for i=l:1：480A0=E(i)*E(i);endR0=m

21、ean(AO);for i=l:1：479A1=E (i + 1);endR (1) =mean (Al);for i=l: 1： 478A2=E(i)*E(i+2);endR (2) =mean (A2);for i=l:1：477A3=E (i)*E(i+3);endR (3) =mean (A3);for i=l: 1： 476A4=E (i)*E (i+4);endR (4 ) =mean (A4);for i=l: 1: 475A5=E(i)*E(i + 5);endR =mean (A5);for i=l: 1： 474A6=E (i)*E (i + 6);endR (6) =me

22、an (A6);for i=l:1：473A7=E(i)*E(i+7);endR =mean (A7);for i=l: 1： 472A8=E(i)*E(i + 8);endR (8)=mean (A8);for i=l:1：471A9=E(i)*E(i + 9);R (9) =mean (A9);for i=l: 1: 470A10=E (i)*E(i + 10); endR (10) =mean (A10);for i=l:1:469A11=E (i)*E(i + ll); endR (11) =mean (All);for i=l:1:468A12=E(i)*E(i + 12);endR

23、 (12) =mean (A12);for i=l:1:467A13=E(i)*E(i + 13); endR (13) =mean (A13);for i = l: 1: 466A14=E(i)*E(i + 14);endR (14) =mean (A14);for i = l: 1: 465A15=E(i)*E(i + 15); endR (1) =mean (Al);for i=l:1:464A16=E(i)*E(i + 16);endR (16) =mean (A16);for i=l: 1: 463A17=E(i)*E(i+17);endR (17) =mean (A17);for i=l: 1:462A18=E (i)*E(i + 18);endR (18)