Matlab系统辨识仿真CH2,CH5,CH7程序注释与剖析

1、第二部分程序注释及运行结果读者须知：为了便于读者理解，现将光盘上第一部分可直接在MATLAB6.I下运行的MATLAB程序的编号和书本上的内容对应如下，每个程序题目括号内的file.m是对应书本上的内容在光盘上第一部分的程序编号。第二章的随机序列产生程序例2.1用乘同余法产生随机数(见光盘FLch2sjxleg1.m)编程如下：A=6;N=100;x0=1;M=255;%初始化；fork=1:Nx2=A*x0;x1=mod(x2,M);%乘向余法递推100次；%x2和x0分力1J表示xi和xi-1；%将*2存储器的数除以M,取余数放x1(xi)中；v1=x1/256;%将*1存储器的数除以25

2、6得到小于1的随机数放v1中；v(:,k)=v1;%将v1中的数(i)存放在矩阵存储器V的第k列中，v(:,k)%表上行不变、列随递推循环次数变化；x0=x1;v0=v1;%xi-1=xi；endv2=v%递推100次结束；%该语句末无；，实现矩阵存储器v中随机数放在v2中，且可直接显示在MATLAB的window中；%grapher%以下是绘图程序；k=1:k1;plot(k,v,k,v,r);xlabel(k),ylabel(v);tktle(0-1)均匀分布的随机序列) 程序运行结果如图2.5所示。 3 B D7 .8 5 ,4 .3 CLN I 0010 2D 30 力口 50I507

3、80901DJk23均勺货布的随机产列图2.5采用MATLAB产生的(0,1)均匀分布的随机序列图 产生的（0-1）均匀分布的随机序列在程序运行结束后，产生的（0,1）均匀分布的随机序列，直接从MATLAB的window界面中copy出来如下（v2中每行存6个随机数）：v2=0.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.5508

4、0.31640.90230.43360.60940.66800.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.9023

5、0.43360.60940.66800.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.0820第二章的白噪声产生程序例2.2用乘同余法产生(见光盘FLch2bzsheg2.m)编程如下：A=6;x0=1;M=255;f=2;N=100;%初始化；x0=1;M=255;fork=1:Nx2=A*x0;x1=mod(x2,M);%乘向余法递推100次；%分别用x2和x0表小xi+1和xi-1;%取*2存储器的数除以M的

6、余数放x1(xi)中；v1=x1/256;%将*1存储器中的数除以256得到小于1的随机数放v1中;v(:,k)=(v1-0.5 )*f;x0=x1;v0=v1; endv2=v%将丫1中的数(i)减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环次数变化；%xi-1=xi；%递推100次结束；%该语句后无；，实现矩阵存储器v中随机数放在v2中,且可直接显示在MATLAB的window中；k1=k;%grapher%以下是绘图程序；k=1:k1;plot(k,v,k,v,r);xlabel(k),ylabel(v);tktle(-1,+1)均匀

7、分布的白噪声)程序运行结果如图2.6所示。图2.6采用MATLAB产生的（-1,+1）均匀分布的白噪声序列产生的（-1,1）均匀分布的白噪声序列在程序运行结束后，产生的（-1,1）均匀分布的白噪声序列，直接从MATLAB的window界面中copy出来如下（v2中每行存6个随机数）：v2=-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.453

8、1-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.921

9、90.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359*另外，书中图2.3白噪声的产生如下:显然，只要在例2.2程序的初始化部分中给N=300,f=6,运行程序就可以得到如图2.3所示的(-3,3)的白噪声过程.编程如下：A=6;x0=1;M=255;f=6;N

10、=300；%初始化；x0=1;M=255;fork=1:Nx2=A*x0;x1=mod(x2,M);%乘向余法递推100次；%分别用x2和x0表示xi+1和xi-1；%取*2存储器的数除以M的余数放x1(xi)中；v1=x1/256;%将*1存储器中的数除以256得到小于1的随机数放v1中；v(:,k)=(v1-0.5)*f;%将团中的数(i)减去0.5再乘以存储器f中的系数，存放在矩阵彳储器V的第k列中，v(:,k)表7K行小艾、列随递推循环次数变化；x0=x1;v0=v1;%xi-1=xi；endv2=v%递推100次结束；%该语句后无；，实现矩阵存储器v中随机数放在v2中，且可直接显示在

11、MATLAB的window中；k1=k;%grapherk=1:k1;plot(k,v,k,v,r);%以下是绘图程序；xlabel(k),ylabel(v);tktle(-1,+1)均匀分布的白噪声)程序运行结果如图2.3所示。图2.3白噪声过程第二章的M序列产生程序例2.3用移位寄存器产生M序列的MATLAB软件实现（见光盘FLch2bzsheg3.m）编程如下：X1=1;X2=0;X3=1;X4=0;%移位寄存器输入Xi初T态（0101）,Yi为移位寄存器各级输出m=60;%置M序列总长度fori=1:m%1#Y4=X4;Y3=X3;Y2=X2;Y1=X1;X4=Y3;X3=Y2;X2=

12、Y1;X1=xor(Y3,Y4);%异或运算ifY4=0U(i)=-1;elseU(i)=Y4;endendM=U%绘图i1=ik=1:1:i1;plot(k,U,k,U,rx)xlabel(k)ylabel（M序列）title（移位寄存器产生的M序列）程序运行结果如图2.8所示。移位寄存器产生的M序列图2.8软件实现的移位寄存器产生的M序列图四级移位寄存器产生的M序列M =Columns 1 through 10-11-11111-1-1-1Columns 11 through 201-1-111-11-111Columns 21 through 3011-1-1-11-1-11Column

13、s 31 through 40-11-111-1-1-1Columns 41 through 501-1-111-11-111Columns 51 through 6011-1-1-11-1-111i1=60第五章的递推的极大似然法辨识程序例5.2系统模型如图5,5所示。试用递推的极大似然法对系统辨识的参数集入 v(k)e(k)图5,5 例5,2系统模型u u(k)v(k)随机信号，输入信号为幅值为1的M序列或随机信号，要求画出程序流程图，打印出程序(程序中带有注释)和辨识中的参数、误差曲线。解：首先解释编程所用的部分字母：由于在MATLAB语言中无法用希腊字母描述、P(k) Pi,(i 0,

14、1) ; K(k) K ;1, 2, 3, 4级寄存器的输出；无法用上标及下标，故依)用o和。1表不；令产生M序列时，a(i),b(i),c(i),d(i)表示四级移位寄存器的第编程如下(光盘上该程序：FLch5RMLeg2.m,可在MATLAB6.I下直接运行)：编程如下：clear%清零a(1)=1;b(1)=0;c(1)=1;d(1)=0;u(1)=d(1);z(1)=0;z(2)=0;%初始化fo门=2:1200%产生m序列u(i)a(i)=xor(c(i-1),d(i-1);b(i)=a(i-1);c(i)=b(i-1);d(i)=c(i-1)；u(i)=d(i)；endu;%若取去

15、；可以在程序运行中观测到m序列v=randn(1200,1);%产生正态分布随机数V=0;%计算噪声方差fo门=1:1200V=V+v(i)*v(i);endV1=V/1200;fork=3:1200%根据v和u计算zz(k)=1.2*z(k-1)-0.6*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);endo1=0.001*ones(6,1);p0=eye(6,6);%赋初值zf(1)=0.1;zf(2)=0.1;vf(2)=0.1;vf(1)=0.1;uf(2)=0.1;uf(1)=0.1;%迭代计算参数值和误差值fork=3:1200h=

16、-z(k-1);-z(k-2);u(k-1);u(k-2);v(k-1);v(k-2);hf=h;K=p0*hf*inv(hf*p0*hf+1);p=eye(6,6)-K*hf*p0;v(k)=z(k)-h*o1;o=o1+K*v(k);p0=p;o1=o;a1(k)=o(1);a2(k)=o(2);b1(k)=o(3);b2(k)=o(4);d1(k)=o(5);d2(k)=o(6);e1(k)=abs(a1(k)+1.2);e2(k)=abs(a2(k)-0.6);e3(k)=abs(b1(k)-1.0);e4(k)=abs(b2(k)-0.5);e5(k)=abs(d1(k)+1.0);

17、e6(k)=abs(d2(k)-0.2);zf(k)=z(k)-d1(k)*zf(k-1)-d2(k)*zf(k-2);uf(k)=u(k)-d1(k)*uf(k-1)-d2(k)*uf(k-2);vf(k)=v(k)-d1(k)*vf(k-1)-d2(k)*vf(k-2);hf=-zf(k-1);-zf(k-2);uf(k-1);uf(k-2);vf(k-1);vf(k-2);endo1%若取去；可以在程序运行中观测到参数V1%绘图subplot(4,1,1)k=1:1200;plot(Ka1,kRa2,bRb1,TRb2,m:Rd1,gRd2,k);xlabel(k)ylabel(para

18、meter)Iegend(a1=-1.2,a2=0.6,b1=1.0,b2=0.5,d1=-1.0,d2=0.2);%图标炷title(TheparameteridendificationoftheRML);endsubplot(4,1,2)k=1:1200;plot(Ke1,kRe2,bRe3,TRe4,m,k,e5,gRe6,k);xlabel(k)ylabel(error)%title(误差曲线)endsubplot(4,1,3)k=1:1200;plot(k,u);xlabel(k)ylabel(input)%title(系统输入信号)endsubplot(4,1,4)k=1:1200

19、;plot(k,v);xlabel(k)ylabel(randomnoise)%title(系统所加的随机噪声)end程序运行结果如图5.7所示reTe moraprnFeThe parameter idendification of the RML1.5k10.502004008001000120001600 k一-50200400600 k80010001200eeFon moan?图5.7RML辨识参数曲线第七章的用改进的神经网络MBPT法辨识例7.1对具有随机噪声的二阶系统的模型辨识(光盘上编号：FLch7NNeg1)对具有随机噪声的二阶系统的模型辨识，进行标幺化以后系统的参考模型差分

20、方程为y(k)a1y(k-1)a2y(k-2)bu(k-1)(k)(7.90)式中，a10.3366,a20.6634,b0.68,(k)为随机噪声。由于神经网络的输出最大为1,所以，被辨识的系统应先标幺化，这里标幺化系数为5。利用图7.5正向建模(并联辨识)结构，神经网络选用3-9-9-1型，即输入层i,隐层j包寸2级，输出层k的节点个数分别为3、9、9、1个；由于神经网络的最大输出为1,因此在辨识前应对原系统参考模型标么化处理，辨识结束后再乘以标么化系数才是被辨识系统的辨识结果。1)编程如下：%w10ij表示第一隐层权值w1ij(k2),w11ij表示w/k1);w120ij表示第二隐层权

21、值w2ij(k2),w121ij表示w21ij(k1);w20j表示输出层权值w3ij(k2),w21j表示亚别尔1)；q表示隐层阈值；p表示输出层阈值；置标幺化系数f1=5等w10ij=.01.01.02;.1.11.02;.010.1;.11.01.02;.1.1.02;.11.1.1;.1.1.1;0.1.1;.10.1;w11ij=.1.2.11;.02.13.04;.09.08.08;.09.1.06;.1.11.02;.060.1;.1.1.1;0.10;.1.1.1;w20j=.01;.02;.1;.2;.1;.1;.1;.1;.1;w21j=0;0.1;.1;.02;0;.1;

22、.1;.1;.1;q0j=.9.8.7.6.1.2.1.1.1;q120j=q0j;q11j=.5.2.3.4.1.2.1.1.1;q12j=q11j;w121ij=w20j*q0j;w120ij=w20j*q11j;f1=5;q2j=0;%thresholdvaluep0=.2;k1=1;p1=.3;w=0;xj=111;%inputserror=0.0001;a1=1111;n=1;e1=0;e0=0;e2=0;e3=0;e4=0;yo=0;ya=0;yb=0;y0=0;y1=0;y2=0;y3=0;u=0;u1=0;u2=0.68;u3=.780;u4=u3-u2;k1=1;kn=28;

23、e3=.055;z1=0;z12=0;q123j=0;t2j=0;o12j=0;r=0;r1=0;s=0.1;d2j=0;%+%calculatingoutputofthehiddenlayerv1=randn(40,1);form=1:40s1=0.1*v1(m)yn=.3366*y2+.6634*u1+s*s1;y1=y2;y2=yn;yp=yn;u0=u1;u1=u2;ya(m)=yn;fork=1:100%calculatingoutputofthehiddenlayer(1)fori=1:9x1=w11ij(i,1)*xj(:,1)+w11ij(i,2)*xj(:,2)+w11ij(

24、i,3)*xj(:,3);x=x1+q11j(:,i);o=1/1+exp(-x);o11j(i)=o；end%calculatingoutputofthehiddenlayer(2)fori=1:9forj=1:9z1=z1+w121ij(i,j)*o11j(:,j);endz=z1+q12j(:,i);o=1/1+exp(-x);o12j(i)=o;end%calculatingoutputoftheoutputlayerfori=1:9yb=yb+w21j(i,:)*o12j(:,i);endyi=yb+p1;y=1/1+exp(-yi);%calculatingerrorvaluebe

25、tweenaimandpracticevaluee0=e1;e1=e2;e2=(yp-y2/2;e(k)=e2;xj1=e2;xj2=e1;xj3=e0;xj=xj1xj2xj3;%revisingrightvalue(1)fori=1:9d1=o11j(:,i)*1-o11j(:,i)*d2j*w21j(i,:);%计算第1隐层误差反传信号do=o11j(:,i)*d1;qw=q11j(:,i)-q0j(:,i);q2j=q11j(:,i)+.8*do+.4*qw;q3j(:,i)=q2j;forj=1:3dw=w11ij(i,j)-w10ij(i,j);w12ij=w11ij(i,j)+.

26、8*do*xj(j)+.6*dw;w13ij(i,j)=w12ij;endendw10ij=w11ij;w11ij=w13ij;q0j=q11j;q11j=q3j;%revisingrightvalue(2)fori=1:9d1=o12j(:,i)*1-o12j(:,i)*d2j*w21j(i,:);%计算第2隐层误差反传信号do=o12j(:,i)*d1;qw=q12j(:,i)-q120j(:,i);t2j=q12j(:,i)+.8*do+.4*qw;q123j(:,i)=t2j;forj=1:9dw=w121ij(i,j)-w120ij(i,j);w122ij=w121ij(i,j)+.

27、8*do*o11j(j)+.6*dw;w123ij(i,j)=w122ij;endendw120ij=w121ij;w121ij=w123ij;q120j=q12j;q12j=q123j;%revisingrightvalue(3)ifm4,r=0.2;门=0.0001;elser=0.14;r1=0.005;endfori=1:9d2j=y*(1-y)*(yp-y);%计算输出误差反传信号dw=w21j(i,:)-w20j(i,:);w22j=w21j(i,:)+r*d2j*o12j(i)+.4*dw+n*e2;w23j(i,:尸w22j;endw20j=w21j;w21j=w23j;ph=

28、p1-p0;p2=p1+.96*(yp-y)+.58*ph+r1*e2;p0=p1;p1=p2;u=y;%k=k+1;ife2=.005break;elseendendya(m尸yp*f1;e3(m)=e2;ym(m尸y*f1;v(m)=s1;%m=m+1m6=mendw11ij=w13ijw121ij=w123ijw21j=w23jm1=m;%graphersubplot(3,1,1)m=1:m6;plot(m,ya,m,ym,rx),title(IdentifiedmodelbyMBPalgorithm),xlabel(k),ylabel(yaandym)legend(yaissystem

29、,ymisidentifiedmodel);%图标炷endsubplot(3,1,2)m=1:m6;plot(m,e3),xlabel(k),ylabel(error)endsubplot(3,1,3)m=1:m6;plot(m,v),xlabel(k),ylabel(randomnoise)end2）辨识结果IdentifiedmodelbyMBPalgorithmya is systemym is identified modelaa y图7.12MBP算法对具有随机噪声的二阶系统辨识结果3)各层权值wllij=w21j=0.23520.48550.2456-0.0589-0.10090.

30、15850.06900.10110.20930.19890.04940.04780.06070.23640.1214-0.15220.10050.12610.0211-0.1189-0.0157-0.15110.09940.04780.09930.09890.09940.0478-0.00070.0989-0.15060.04780.09930.24890.09940.0478w121ij=0.00230.00230.00230.01650.01820.01430.01040.00230.00330.02480.02980.02130.0129-0.0026-0.0008-0.0026-0.0

31、026-0.00270.14660.16720.12690.08660.00690.01680.00690.00700.00690.30430.34340.26380.18420.02380.04390.02380.02370.02380.15350.17280.13310.09340.01310.02310.01310.01300.01310.14660.16720.12690.08660.00690.01680.00690.00700.00690.14660.16720.12690.08660.00690.01680.00690.00700.00690.14660.16720.12690.

32、08660.00690.01680.00690.00700.00690.14660.16720.12690.08660.00690.01680.00690.00700.0069第七章的多维非线性辨识的MATLAB程序例7.2用神经网络学习二维非线性函数ypcos(2k1m2kzn)sin(2kzn)(光盘上编号：FLch7NNeg2)%初始化:w10ij表示Wij(k-2);w11ij表示wij(k-1);w20j表示Wjk(k-2);w21j表示wjk(k-1);q和p为阈值w10ij=.01.01.02;.01.01.02;.010.01;.001.001.002;.0010.002;.0

33、011.001.001;w11ij=-.1-.02.11;-.21.10-.19;-.14.15-.16;.14-.13.17;-.13.12.21;-.16-.23.13;w20j=.01;.02;.1;.2;.1;.1;w21j=.18;.9;.9;.7;.8;.9;q0j=.5.8.4.6.1.2;q1j=-.1.02.12.14-.02.02;q2j=0;%thresholdvaluep0=.2;k1=1;p1=.1;w23j=0;0;0;0;0;0;w22j=0;w=0;xj=0.50.30.2;%inputsya=000;yp=0;yy=0;m1=0;yam=0;yp1=0;qw=

34、0;yo=000;ya1=0;error=0.0001;n=1;q=0;e1=0;e0=0;e2=0;e3=0;e4=0;yo=0;ya=0;yb=0;y0=0;y1=0;y2=0;y3=0;u=0;u1=0;u2=0;k1=1;w0=1;kn=28;dj2=0.01;e3=.0055;a1=0.036;a2=0.036;a3=0.08;w0=1;dj2=0.01;%+初始化结束form=1:36forn=1:36q=pi*0.05*m;p=pi*0.02*n;yn=cos(2*q+2*p)*sin(2*p);yzs=0.1*yn;%xj=0.4+yzs0.30.2;%协调器使网络的输入随样本

35、动态变化yp=yn;fork=1:6%calculatingoutputofthehiddenlayerfori=1:6x1=w11ij(i,1)*xj(:,1)+w11ij(i,2)*xj(:,2)+w11ij(i,3)*xj(:,3);x=x1+q1j(:,i);o=tanh(x);o1j(i)=o;end%calculatingoutputoftheoutputlayerfori=1:3yb=yb+w21j(i,:)*o1j(:,i);endyi=yb+p1;y=tanh(yi);%calculatingerrorvaluebetweenaimandpracticevaluee0=e1;

36、e1=e2;e2=(yp-y).A2/2;%revisingrightvaluefori=1:6d1=1-o1j(:,i)*dj2*w23j(i,:);%计算隐层误差反传信号do=o1j(:,i)*d1;q3j(:,i尸q1j(i);forj=1:3dw=w11ij(i,j)-w10ij(i,j);ife20.004,a1=0;a2=0;a3=0;else,a1=0.022;a2=0.03;a3=0.003;end;w12ij=w11ij(i,j)+a1*do*xj(j)+a2*dw+a3*w;w13ij(i,j)=w12ij;endendw10ij=w11ij;w11ij=w13ij;q0j

37、=q1j;q1j=q3j;w=yp-y;w0=w;w=0.36*w0+(yp-y);%计算积累误差wife20.004,w=0.78*w;endfori=1:6d2j=y*(1-y)*(yp-y);%计算输出层误差反传信号dw=w21j(i,:)-w20j(i,:);w22j=w21j(i,:)+0.132*d2j*o1j(i)+0.26*dw;+0.026*w;w23j(i,:)=w22j;endw20j=w21j;w21j=w23j;ph=p1-p0;p2=p1+.3*(yp-y)+.102*ph;p0=p1;p1=p2;ife2=0.005break;elseend%判断误差end%对样

38、本的一次采样值训练k次结束ypp(m,n)=yn;%存储二维非线性样本yom(m,n尸y;%存储神经网络辨识结果e3(m,n)=e2;%存储训练误差end%变量n训练结束m2=mend%变量m训练结束w11ij=w13ij%句末无“；”，可直接在MATLAB主界面下观察隐层权值w21j=w23j%可直接在MATLAB主界面下观察输出层权值%graphersubplot(2,2,1);%plot3(m,n,ypp);%也可用此语句绘制样本的三维图mesh(ypp)%绘制样本的三维图xlabel(m),ylabel(n),zlabel(ypp),%三维坐标title(Identifiedmodel

39、byinp.algorithm),%图题subplot(2,2,2);mesh(yom)xlabel(m),ylabel(n),zlabel(yom);subplot(2,1,2);mesh(e3)xlabel(m),ylabel(n),zlabel(e3);3)辨识结果Identified model by inp. algorithm3 enm图7.13对二维非线性函数ypcos(2k1m2k2n)sin(2k2n)的辨识结果3)结论在程序设计和调试过程中，注意以下技巧：只要在程序中增加类似m或n的循环语句，便可用神经网对三维或三维以上的非线形函数训练，也可以用在复杂系统二个以上独立变量、

40、任意维非线性动态参数的辨识。程序中除了for,end,else和mesh语句外，一般应在语句末加分号。程序在m2=m,w11ij=w13ij和w21j=w23j不要加分号，目的有二个，一是用m2=m动态的显示程序执行的状态及程序每次循环所花的时间；二是在程序运行结束显示隐层和输出层的权值：w11ij=w21j=-0.00910.07580.20830.8293-0.12060.1949-0.09241.5274-0.01590.2657-0.04850.31840.2347-0.03190.2698-0.9167-0.01420.23080.31820.6995-0.0755-0.13810.

41、22571.5786调整修正因子a1a223和6b2b3,直到方差最小，收敛时间最短为止。一旦修正因子、循环次数m和n确定，调整第三层内循环k的次数，可以明显地改变训练精度和程序运行速度，k=20,最大方差e3max=0.01,程序运行时间ts=14s;k=6,e3max=0.032,ts=12s;k=3,e3max=0.08,ts=10s;k=2,e3=0.22,ts=9s.程序中，语句：ife20.05,a1=0;a2=0;a3=0;else,a1=0.02;a2=0.05;a3=0.005;end;协调器根据训练误差动态变化学习因子，如果在这条语句前加“”(即不要这条语句)，训练误差将参

42、加一倍多。这种改进BP网络比一般的BP网络在m,n,k相同的情况下平均收敛速度约提高两秒，且训练精度明显提高。第七章的模糊神经网络解耦MATLAB程序例7.3用隶属函数型神经网与模糊控制融合的解耦程序(FLch7FNNeg3.m)具有耦合的两个相邻子系统的差分方程为y(k)ay(k1)bhu(k1)bu*2)(k)yk)(7.122)y2(k)a21y(k1)b21u(k1)b22u(k2)2(k)yk)式中，i(k)随机噪声；yj(k)为两个相邻子系统之间的耦合；采用隶属函数型神经网与模糊控制融合的解耦方法(DMFFCNN方法实现解耦控制。1)开发的程序p1=0.5;q1=0.5;p2=0.

43、5;q2=0.5;p3=0.5;q3=0.5;p4=0.95;q4=0.5;p5=0.95;q5=0.5;p6=0.95;q6=0.5;p7=0.5;q7=0.5;p8=0.95;q8=0.5;p9=0.5;q9=0.5;p10=0.95;q10=0.5;p11=0.95;q11=0.5;p12=0.95;q12=0.5;p13=0.65;q13=0.05;%隶属函数型神经网的中心值、尺度因子和权向量初始化部分a10=-202;a11=-202;a20=-202;a21=-202;b10=1.51.51.5;b11=1.51.51.5;b20=1.51.51.5;b21=1.51.51.5;v

44、0=-1-0.5-0.5;-0.500.5;0.50.51;v1=-1-0.5-0.5;-0.500.5;0.50.51;%系统部分的初值yp0=0;yp1=0;ep0=0;ep1=0;up0=1.05;up1=1.39;x1=0.12;x2=0.24;y0=0;y1=0;u1=0;e0=0;e1=0;sp=10;k=1;se=0.02;sd=0.02;su=0.0522;%开始kp=0.284;ki=0.03;kd=0.626;%thesecondchannelreferencep12=0.5;q12=0.5;p22=0.5;q22=0.5;p32=0.5;q32=0.5;p42=0.95;

45、q42=0.5;p52=0.95;q52=0.5;p62=0.95;q62=0.5;p72=0.5;q72=0.5;p82=0.95;q82=0.5;p92=0.5;q92=0.5;p102=0.95;q102=0.5;p112=0.95;q112=0.5;p122=0.95;q122=0.5;p132=0.65;q132=0.05;%系统部分的初值yp02=0;yp12=0;ep02=0;ep12=0;up02=1.05;up12=1.39;x12=0.12;x22=0.24;y02=0;y12=0;u12=0;e02=0;e12=0;sp2=5;k=1;se2=0.02;sd2=0.02;

46、su2=0.0522;kp2=0.284;ki2=0.03;kd2=0.626;%+%子系统1(SUB1)的J循环开始forJ=1:50ep1=10-yp1;pid=kp*(ep1-ep0)+ki*ep1;up2=up1+pid;yp2=0.5*yp1+2.5*up2+2.5*up1;%yp2=0.8333*yp1-0.1667*yp0+1.6667*up2+1.6667*up1;yp(:,J尸yp2;up0=up1;up1=up2;ep0=ep1;yp0=ypi;ypi=yp2;endtime=1:1:50;n=0.28*rand(size(time);%产生的SUB1随机噪声n1=0.3*

47、rand(size(time);%产生的SUB2随机噪声while(k=50)%functiondealwithfor(T=1:20)X1=x1*se;X2=x2*sd;if(X1=2)X1=2;endif(X2=2)X2=2;end%thehiddeninglayeroutputofFNNfori=1:3forj=1:3A=(X1-a11(:,i)/b11(:,i).A2;B=(X2-a21(:,j)/b21(:,j)A2;h(i,j尸exp(-(A+B)/2);endend%theoutputsum=0;fori=1:3forj=1:3sum=sum+h(i,j)*v1(i,j);ende

48、ndot=sum;cu=su*ot;u2=u1+cu;disp(u2);if(u2=1)u2=1;end%+n1(:,k)+0.01*y12表示随机噪声和子系y2=0.5*y1+2.5*u2+2.5*u1+n1(:,k)+0.01*y12;统间的相互耦合disp(theoutputynumberisint2str(T);disp(y2);Y(：,k尸y2;E=0.5*(sp-y2).A2;e2=sp-y2;Es(:,k)=e2;x1=e2;x2=e2-e1;e0=e1;e1=e2;delot=(sp-y2)*2.5*su;fori=1:3forj=1:3dv=v1(i,j)-v0(i,j);v

49、2=v1(i,j)+p13*delot+q13*dv;v3(i,j尸v2;endend%refreshingthecenterandwidth%a11s1=0;forj=1:3s1=s1+v1(1,j)*h(1,j);endpa11=s1;a110=a10(:,1);a111=a11(:,1);dela11=pa11*(X1-a111)/b11(:,1).A2;da11=a111-a110;a112=a111+p1*dela11+q1*da11;a10(:,1)=a111;a11(:,1)=a112;%a12s2=0;forj=1:3s2=s2+v1(2,j)*h(2,j);endpa12=s

50、2;a120=a10(:,2);a121=a11(:,2);dela12=pa12*(X1-a121)/b11(:,2).A2;da12=a121-a120;a122=a121+p2*dela12+q2*da12;a10(:,2)=a121;a11(:,2)=a122;%a13s3=0;forj=1:3s3=s3+v1(3,j)*h(3,j);endpa13=s3;a130=a10(:,3);a131=a11(:,3);dela13=pa13*(X1-a131)/b11(:,32;da13=a131-a130;a132=a131+p3*dela13+q3*da13;a10(:,3)=a131;a11(:,3)=a132;%b11pb11=pa11;b110=b10(:,1);b111=b11(:,1);delb11=pb11*(X1-a111).A2/b111.A3;db11=b111-b110;b112=b111+p4*delb11+q4*db1