MATLAB是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。对于一个简单的系统 ,可以通过分析其过程的运动规律 ,应用一些已知的定理和原理,建立数学模型 ,即所谓的“白箱建模 ”。但对于比较复杂的生产过程 ,该建模方法有很大的局限性。由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。这种建模方法就称为系统辨识。把辨识建模称作“黑箱建模”。本章列举了在系统辨识中常用的基本算法的MATLAB程序。
8.1 白噪声的产生
如果在计算机上比较经济的产生统计上理想的各种不同分布的白噪声序列，则对系统辨识仿真研究提供了极大方便。为了简单起见，常把各种不同分布的白噪声序列称为随机数。从理论上讲，只要有了一种具有连续分布的随机数，就可以通过函数变换的方法产生其他任意分布的随机数。显然，在具有连续分布的随机数中，（0,1）均匀分布的随机数是最简单、最基本的一种，有了（0,1）均匀分布的随机数，便可以产生其他任意分布的随机数和白噪声。
本文首先产生（0,1）均匀分布的随机序列，然后对程序稍加改动，将产生（0,1）均匀分布的随机数统统减去0.5，相当于将原随机序列图的横坐标向上平移0.5，原随机序列变成了（-0.5,0.5）的白噪声，然后乘以存储器f的系数，这里取法f=2，则可得到产生（-1,1）均匀分布的白噪声。
8.1.1 乘同余法的介绍
用如下递推同余式产生正整数序列{}
（8.1）
式中，mod为取M的余数，如：M为2的方幂，即，k为k>2是整数；若k=3，则（mod8）或k=8，则（mod256），且A应取适中的值，如3 (8.2) 可以证明序列{}是伪随机数序列；同时还可以证明伪随机序列{}的循环周期达到最大值。 将式（8.1）和（8.2）合并为 （8.3） 式中，初值为。 8.1.2 举例 利用乘同余法，选R=2,A=6,k=8, ,递推100次，采用MATLAB的仿真语言编程，产生（-1,1）均匀分布的白噪声。 用乘同余法产生(见配套MATLAB程序 algorithms1.m) 1编程如下： A=6; x0=1; M=255; f=2; N=100； %初始化； x0=1; M=255; for k=1: N %乘同余法递推100次； x2=A*x0; %分别用x2和x0表示xi+1和xi-1； x1=mod (x2,M); %取x2存储器的数除以M的余数放x1（xi）中； v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中； v(:,k)=(v1-0.5 )*f; %将v1中的数（）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环次数变化； x0=x1; % xi-1= xi； v0=v1; end %递推100次结束； v2=v %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MATLAB的window中； k1=k; %grapher %以下是绘图程序； k=1:k1; plot(k,v,k,v,'r'); xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声') ② 程序运行结果如图8.1所示。 图8.1 采用MATLAB产生的(-1,+1)均匀分布的白噪声序列 ③ 产生的(-1，1)均匀分布的白噪声序列 在程序运行结束后，产生的(-1，1)均匀分布的白噪声序列，直接从MATLAB的window界面中copy出来如下（v2中每行存6个随机数）： v2 = -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 8.2 M序列的产生 在进行系统辨识时，选用白噪声作为辨识输入信号可以保证获得较好的辨识效果，但在工程上难以实现。M序列是一种很好的辨识输入信号，它具有近似白噪声的性质，不仅可以保证有较好的辨识效果，而且工程上又易于实现。 8.2.1 M序列的产生方式 M序列是一种离散二位式随机序列，所谓“二位式”是指每个随机变量只有两种状态。离散二位式随机序列是按照确定的方式产生的，实际上是一种确定序列。可用多级线性反馈移位寄存器产生M序列。每级移位寄存器由双稳态触发器和门电路组成，称为1位，分别以0和1来表示2中状态。当移位脉冲来到时，每位的内容（0或1）移到下一位，最后1位（即第n位）移出的内容即为输出。为了保持连续工作，将最后2级寄存器的内容经过适当的逻辑运算后反馈到第1级寄存器作为输入。 8.2.2 用四级移位寄存器产生M序列 用移位寄存器产生M序列的MATLAB软件实现(见配套MATLAB程序algorithms2.m) ① 编程如下： X1=1;X2=0;X3=1;X4=0; %移位寄存器输入Xi初T态（0101）， Yi为移位寄存器各级输出 m=60; %置M序列总长度 for i=1:m %1# Y4=X4; Y3=X3; Y2=X2; Y1=X1; X4=Y3; X3=Y2; X2=Y1; X1=xor(Y3,Y4); %异或运算 if Y4==0 U(i)=-1; else U(i)=Y4; end end M=U %绘图 i1=i k=1:1:i1; plot(k,U,k,U,'rx') xlabel('k') ylabel('M序列') title('移位寄存器产生的M序列') ② 程序运行结果如图8.2所示。 图8.2 软件实现的移位寄存器产生的M序列图 ③ 四级移位寄存器产生的M序列 M = Columns 1 through 10 -1 1 -1 1 1 1 1 -1 -1 -1 Columns 11 through 20 1 -1 -1 1 1 -1 1 -1 1 1 Columns 21 through 30 1 1 -1 -1 -1 1 -1 -1 1 1 Columns 31 through 40 -1 1 -1 1 1 1 1 -1 -1 -1 Columns 41 through 50 1 -1 -1 1 1 -1 1 -1 1 1 Columns 51 through 60 1 1 -1 -1 -1 1 -1 -1 1 1 i1 = 60 8.3 最小二乘一次完成算法的产生 考虑仿真对象 (8.4)其中，是服从正态分布的白噪声N(0,1)。输入信号采用4阶M序列，幅度为1。 选择如下形式的辨识模型 (8.5) 设输入信号的取值是从k =1到k =16的M序列，则待辨识参数为=。其中，被辨识参数、观测矩阵z L、H L的表达式为 , ， (8.6) 程序框图如图8.3所示。 Matlab仿真程序如下： %二阶系统的最小二乘一次完成算法辨识程序，（见配套MATLAB程序algorithms3.m) u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画输入信号u的径线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线 subplot(3,1,3) %画三行一列图形窗口中的第三个图形 stem(z),grid on %画出输出观测值z的径线图形，并显示坐标网格 u,z %显示输入信号和输出观测信号 %L=14 %数据长度 HL=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14)] %给样本矩阵HL赋值 ZL=[z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)] % 给样本矩阵z L赋值 %Calculating Parameters c1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c=c2*c3 %计算并显示 %Display Parameters a1=c(1), a2=c(2), b1=c(3),b2=c(4) %从中分离出并显示a1 、a2、 b1、 b2 %End 程序运行结果： >> u =[ -1，1，-1，1，1，1，1，-1，-1，-1，1，-1，-1，1，1] z =[ 0，0，0.5000，0.2500，0.5250，2.1125， 4.3012，6.4731，6.1988，3.2670，-0.9386， -3.1949，-4.6352，6.2165，-5.5800，-2.5185] HL = 0 0 1.0000 -1.0000 -0.5000 0 -1.0000 1.0000 -0.2500 -0.5000 1.0000 -1.0000 -0.5250 -0.2500 1.0000 1.0000 -2.1125 -0.5250 1.0000 1.0000 -4.3012 -2.1125 1.0000 1.0000 -6.4731 -4.3012 -1.0000 1.0000 -6.1988 -6.4731 -1.0000 -1.0000 -3.2670 -6.1988 -1.0000 -1.0000 0.9386 -3.2670 1.0000 -1.0000 3.1949 0.9386 -1.0000 1.0000 4.6352 3.1949 -1.0000 -1.0000 6.2165 4.6352 1.0000 -1.0000 5.5800 6.2165 1.0000 1.0000 ZL =[ 0.5000，0.2500，0.5250，2.1125，4.3012，6.4731，6.1988，3.2670，-0.9386，-3.1949， -4.6352，-6.2165，-5.5800，-2.5185]T c =[ -1.5000，0.7000，1.0000，0.5000]T a1 = -1.5000 a2 = 0.7000 b1 = 1.0000 b2 =0.5000 >> 从仿真结果表8.4可以看出，由于所用的输出观测值没有任何噪声成分，所以辨识结果也无任何误差。 8.4 最小二乘递推算法的产生 考虑图8.5所示的仿真对象， 图中，是服从N（0,1）分布的不相关随机噪声。且 ， ， 选择图8.5所示的辨识模型。仿真对象选择如下的模型结构 (8.7) 其中，是服从正态分布的白噪声N（0,1）。输入信号采用4位移位寄存器产生的M序列，幅度为0.03。按式 (8.8) 构造h (k)；加权阵取单位阵；计算K(k)、和P(k)，计算各次参数辨识的相对误差，精度满足要求后停机。 最小二乘递推算法辨识的Malab程序流程如图8.6所示。 图8.6 最小二乘递推算法辨识的Malab6.0程序流程图 下面给出具体程序。 %最小二乘递推算法辨识程序, （见配套MATLAB程序algorithms4.m) clear %清理工作间变量 L=15; % M序列的周期 y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值 for i=1:L;%开始循环，长度为L x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或” x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出 x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出 x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出 y(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列 if y(i)>0.5,u(i)=-0.03; %如果M序列的值为"1", 辨识的输入信号取“-0.03” else u(i)=0.03; %如果M序列的值为"0", 辨识的输入信号取“0.03” end %小循环结束 y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备 end %大循环结束，产生输入信号u figure(1); %第一个图形 stem(u),grid on %显示出输入信号径线图并给图形加上网格 z(2)=0;z(1)=0; %设z的前两个初始值为零 for k=3:15; %循环变量从3到15 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %输出采样信号 end %RLS递推最小二乘辨识 c0=[0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量 p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵 E=0.000000005; %取相对误差E=0.000000005 c=[c0,zeros(4,14)]; %被辨识参数矩阵的初始值及大小 e=zeros(4,15); %相对误差的初始值及大小 for k=3:15; %开始求K h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; x=h1'*p0*h1+1; x1=inv(x); %开始求K(k) k1=p0*h1*x1;%求出K的值 d1=z(k)-h1'*c0; c1=c0+k1*d1; %求被辨识参数c e1=c1-c0; %求参数当前值与上一次的值的差值 e2=e1./c0; %求参数的相对变化 e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数 c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1'*[h1'*p0*h1+1]; %求出 p(k)的值 p0=p1; %给下次用 if e2<=E break; %如果参数收敛情况满足要求，终止计算 end %小循环结束 end %大循环结束 c，e %显示被辨识参数及其误差(收敛)情况 %分离参数 a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); figure(2); %第二个图形 i=1:15; %横坐标从1到15 plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出a1，a2，b1，b2的各次辨识结果 title('Parameter Identification with Recursive Least Squares Method') %图形标题 figure(3); %第三个图形 i=1:15; %横坐标从1到15 plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况 title('Identification Precision') %图形标题 程序运行结果： c = 0.0010 0 0.0010 -0.4984 -1.2328 -1.4951 -1.4962 -1.4991 -1.4998 -1.4999 0.0010 0 0.0010 0.0010 -0.2350 0.6913 0.6941 0.6990 0.6998 0.6999 0.0010 0 0.2509 1.2497 1.0665 1.0017 1.0020 1.0002 0.9999 0.9998 0.0010 0 -0.24 0.7500 0.5668 0.5020 0.5016 0.5008 0.5002 0.5002 -1.5000 -1.5000 -1.5000 -1.4999 -1.4999 0.7000 0.7000 0.7000 0.7000 -0.7000 0.9999 0.9999 0.9999 0.9999 0.9999 0.5000 0.5000 0.5000 0.5000 0.5000 e = 0 0 0 -499.4200 1.4734 0.2128 0.0007 0.0020 0.0004 0.0000 0 0 0 0 -235.9916 -3.9416 0.0042 0.0070 0.0012 0.0001 0 0 249.8612 3.9816 -0.1466 -0.0607 0.0003 -0.0018 -0.0003 -0.0001 0 0 -249.8612 -4.0136 -0.2443 -0.1143 -0.0007 -0.0016 -0.0012 -0.0001 0.0001 0.0000 -0.0000 -0.0000 0.0000 0.0001 -0.0000 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.0000 -0.0000 -0.0004 0.0000 -0.0000 0.0000 -0.0000 图8.7 最小二乘递推算法的参数辨识数 图8.8 最小二乘递推算法的参数辨识结果 图8.9 最小二乘递推算法的参数辨识结果收敛情况 表8.2 最小二乘递推算法的辨识结果 参 数 a1 a2 b1 b2 真 值 -1.5 0.7 1.0 0.5 估计值 -1.4999 0.7 0.9999 0.5000 仿真结果表明，大约递推到第十步时，参数辨识的结果基本达到稳定状态，即a1=-1.4999, a2= 0.7000, b1=0.9999, b2=0.5000。此时，参数的相对变化量E 0.000000005。从整个辨识过程来看，精度的要求直接影响辨识的速度。虽然最终的精度可以达到很小，但开始阶段的相对误差会很大，从图8.7(3)图形中可看出，参数的最大相对误差会达到三为数 8.5 广义最小二乘算法的产生 广义最小二乘法所用的模型是： （8.9） （8.10） 其中我们选定模型参数：a1=1.5,a2=-0.7,b1=1,b2=0.5,c1=0,c2=0 广义最小二乘递推算法的计算步骤： 1. 给定初始条件 2. 利用式，计算和； 3. 利用式，构造； 4. 利用式递推计算； 5. 利用和，计算； 6. 根据来构造； 7. 利用 广义最小二乘法程序代码（见配套MATLAB程序algorithms5.m) clear %清理工作间变量 L=55; % M序列的周期 y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值 for i=1:L;%开始循环，长度为L x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或” x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出 x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出 x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出 y(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列 if y(i)>0.5,u(i)=-1; %如果M序列的值为"1", 辨识的输入信号取“-1” else u(i)=1; %如果M序列的值为"0", 辨识的输入信号取“1” end %小循环结束 y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备 end %大循环结束，产生输入信号u z(2)=0;z(1)=0; %设z的前四个初始值为零 for k=3:45; %循环变量从5到45 z(k)=1.5*z(k-1)-0.7*z(k-2)+1*u(k-1)+0.5*u(k-2); %输出采样信号 end %RGLS广义最小二乘辨识 c0=[0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量 pf0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵 ce0=[0.001 0.001]'; pe0=eye(2,2); c=[c0,zeros(4,44)]; %被辨识参数矩阵的初始值及大小 ce=[ce0,zeros(2,44)]; e=zeros(4,45); %相对误差的初始值及大小 ee=zeros(2,45); for k=3:45; %开始求K zf(k)=z(k)+0*z(k-1)+0*z(k-2); uf(k)=u(k)+0*u(k-1)+0*u(k-2); hf1=[-zf(k-1),-zf(k-2),uf(k-1),uf(k-2)]'; x=hf1'*pf0*hf1+1; x1=inv(x); %开始求K(k) k1=pf0*hf1*x1;%求出K的值 d1=zf(k)-hf1'*c0; c1=c0+k1*d1; %求被辨识参数c e1=c1-c0; %求参数当前值与上一次的值的差值 e2=e1./c0; %求参数的相对变化 e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数 c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列 pf1=pf0-k1*hf1'*pf0; %求出 p(k)的值 pf0=pf1; %给下次用 h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; ee(k)=z(k)-h1'*c1; he1=[-ee(k-1),-ee(k-2)]'; x=he1'*pe0*he1+1; x1=inv(x); k1=pe0*he1*x1; d1=ee(k)-he1'*ce0; ce1=ce0+k1*d1; pe1=pe0-k1*he1'*pe0; ce0=ce1; ce(:,k)=ce1; pe0=pe1; end %大循环结束 %显示被辨识参数及其误差(收敛)情况 %分离参数 a1=c(1,:); a2=c(2,:);b1=c(3,:);b2=c(4,:); c1=ce(1,:);c2=ce(2,:); ea1=e(1,:); ea2=e(2,:);eb1=e(3,:); eb2=e(4,:); figure(2); %第二个图形 i=1:45; %横坐标从1到25 plot(i,a1,'r',i,a2,':',i,b1,'b',i,b2,':',i,c1,'b',i,c2,':') %画出a1，a2,b1，b2,c1,c2的各次辨识结果 title('参数变化曲线') %图形标题 figure(3); %第三个图形 i=1:45; %横坐标从1到25 plot(i,ea1,'r',i,ea2,'k:',i,eb1,'g',i,eb2,'c:') %画出a1，a2，b1，b2,的各次辨识结果的收敛情况 title('误差曲线') %图形标题 图8.10 广义最小二乘法辨识参数变化曲线 图8.11 广义最小二乘法辨识误差变化曲线 结果分析：广义最小二乘法的基本思想是基于对数据先进行一次滤波预处理，然后利用普通的最小二乘法对滤波后的数据进行辨识。辨识结果a1=1.5016,a2=-0.7010,b1=1.0120,b2=0.4986,c1=0.0040,c2=0.0027与最小二乘法递推算法相比较，可以发现它们的结果基本上是一致的。这是因为本例的仿真对象相当于 =1。这种对象利用最小二乘法已可获得无偏一致估计，但广义最小二乘法同时又能给出噪声模型的参数估计。 8.6 增广最小二乘递推算法的产生程序 考虑图8.12所示的仿真对象，图中，是服从N（0,1）分布的不相关随机噪声。 且 ， ， 模型结构选用如下形式： 增广最小二乘递推辨识程序流程如图8.13所示。 图8.13 增广最小二乘递推算法辨识的Malab6.0程序流程图 Matlab程序如下。 %增广最小二乘递推辨识程序，（见配套MATLABalgorithms6.m) L=60; %4位移位寄存器产生的M序列的周期 y1=1;y2=1;y3=1;y4=0; %4个移位寄存器的输出初始值 for i=1:L; x1=xor(y3,y4); %第一个移位寄存器的输入信号 x2=y1; %第二个移位寄存器的输入信号 x3=y2; %第三个移位寄存器的输入信号 x4=y3; %第四个移位寄存器的输入信号 y(i)=y4; %第四个移位寄存器的输出信号，M序列, 幅值"0"和"1"， if y(i)>0.5,u(i)=-1; %M序列的值为"1"时,辨识的输入信号取“-1” else u(i)=1; %M序列的值为"0"时,辨识的输入信号取“1” end y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号准备 end figure(1); %画第一个图形 subplot(2,1,1); %画第一个图形的第一个子图 stem(u),grid on %画出M序列输入信号 v=randn(1,60); %产生一组60个正态分布的随机噪声 subplot(2,1,2); %画第一个图形的第二个子图 plot(v),grid on; %画出随机噪声信号 u ，v %显示输入信号和噪声信号 z=zeros(7,60);zs=zeros(7,60);zm=zeros(7,60);zmd=zeros(7,60); %输出采样矩阵、不考虑噪声时系统输出矩阵、不考虑噪声时模型输出矩阵、模型输出矩阵的大小 z(2)=0; z(1)=0; zs(2)=0; zs(1)=0; zm(2)=0; zm(1)=0; zmd(2)=0; zmd(1)=0; %输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出的初值 c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量 p0=10^6*eye(7,7); %直接给出初始状态P0，即一个充分大的实数单位矩阵 E=0.00000000005; %相对误差E=0.000000005 c=[c0,zeros(7,14)]; %被辨识参数矩阵的初始值及大小 e=zeros(7,15); %相对误差的初始值及大小 for k=3:60; %开始求K z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2); %系统在M序列输入下的输出采样信号 h1=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]'; %为求K(k)作准备 x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %K d1=z(k)-h1'*c0; c1=c0+k1*d1; %辨识参数c zs(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2); %系统在M序列的输入下的输出响应 zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[c1(1);c1(2);c1(3);c1(4)]; %模型在M序列的输入下的输出响应 zmd(k)=h1'*c1; %模型在M序列的输入下的输出响应 e1=c1-c0; e2=e1./c0; %求参数误差的相对变化 e(:,k)=e2; c0=c1; %给下一次用 c(:,k)=c1; %把递推出的辨识参数c 的列向量加入辨识参数矩阵 p1=p0-k1*k1'*[h1'*p0*h1+1]; %find p(k) p0=p1; %给下次用 if e2<=E break; %若收敛情况满足要求，终止计算 end %判断结束 end %循环结束 c, e, %显示被辨识参数及参数收敛情况 z, zs, zm %显示输出采样值、系统实际输出值、模型输出值 %分离变量 a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); %分离出a1、 a2、 b1、 b2 d1=c(5,:); d2=c(6,:); d3=c(7,:); %分离出d1、 d2、 d3 ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); %分离出a1、 a2、 b1、 b2的收敛情况 ed1=e(5,:); ed2=e(6,:); ed3=e(7,:); %分离出d1 、d2 、d3的收敛情况 figure(2); %画第二个图形 i=1:60;plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:',i,d3,'g+') %画出各个被辨识参数 title('Parameter Identification with Recursive Least Squares Method') %标题 figure(3); i=1:60; %画出第三个图形 plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:',i,ed2,'r+') %画出各个参数收敛情况 title('Identification Precision') %标题 figure(4); subplot(4,1,1); %画出第四个图形，第一个子图 i=1:60;plot(i,zs(i),'r') %画出被辨识系统在没有噪声情况下的实际输出响应 subplot(4,1,2); i=1:60;plot(i,z(i),'g') %第二个子图，画出被辨识系统的采样输出响应 subplot(4,1,3); i=1:60;plot(i,zm(i),'b') %第三个子图，画出模型含有噪声的输出响应 subplot(4,1,4); i=1:60;plot(i,zs(i),'b') %第四个子图，画出模型去除噪声后的输出响应 程序执行结果： >> u = [1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,1,-1,-1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1, -1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1] v = [1.1418,1.5519,1.3836,-0.7581,0.4427,0.9111,-1.0741,0.2018,0.7629,-1.2882,-0.9530,0.7782,-0.0063,0.5245, 1.33,0.4820,-0.7871,0.7520,-0.1669,-0.8162,2.0941,0.0802,-0.9373,0.6357,1.6820,0.5936,0.7902,0.1053, -0.1586,0.8709,-0.1948,0.0755,-0.5266,-0.6855,-0.2684,-1.1883,0.2486,-0.1025,-0.0410,-2.2476,-0.5108, 0.2492 0.3692,0.1792,-0.0373,-1.6033,0.3394,-0.1311,0.4852,0.5988,-0.0860,0.3253,-0.3351,-0.3224,-0.3824,-0.9534, 0.2336,1.2352,-0.5785,-0.5015] c = 0.0010 0 0.0010 -0.27 -0.9025 -0.9390 -0.8620 -1.1236 -1.5000 -1.5000 -1.5000 0.0010 0 0.0010 0.0010 -0.0734 -0.2661 -0.1199 0.3259 0.7000 0.7000 0.7000 0.0010 0 0.0592 0.4560 0.5300 0.5293 0.9458 0.5402 1.0000 1.0000 1.0000 0.0010 0 0.0572 0.8186 0.8500 0.8639 1.1157 1.2957 0.5000 0.5000 0.5000 0.0010 0 0.0796 0.7342 0.5417 0.5516 0.3847 0.2459 1.0000 1.0000 1.0000 0.0010 0 -0.04 -0.5981 -0.3420 -0.4456 -0.4044 -0.8443 -1.0000 -1.0000 -1.0000 0.0010 0 -0.0655 -0.7795 -0.8572 -0.7412 -0.4506 -0.1974 0.2000 0.2000 0.2000 e = 0 0 0 -279.9426 2.2355 0.0404 -0.0820 0.3035 0.3350 0.0000 0.0000 0 0 0 0 -74.37 2.6263 -0.5495 -3.71 1.1479 -0.0000 0.0000 0 0 58.2212 6.7007 0.1622 -0.0013 0.7868 -0.42 0.8512 0.0000 0.0000 0 0 -58.2212 -15.3052 0.0384 0.0163 0.2915 0.1613 -0.6141 -0.0000 -0.0000 0 0 -80.5566 -10.2283 -0.2621 0.0183 -0.3027 -0.3607 3.0663 0.0000 0.0000 0 0 -90.3556 5.6936 -0.4282 0.3029 -0.0925 1.0880 0.1844 -0.0000 0.0000 0 0 -66.4754 10.9059 0.0996 -0.1354 -0.3920 -0.5619 -2.0130 0.0000 0.0000 z = [0,0,-0.4400,-3.9913,-5.7014,-6.9415,-7.8178,-3.9097,1.4543,2.4075,3.5810,6.6598,6.0079,3.53,2.4376, -0.09,-2.3472,-2.3178,-4.4100,-6.9915,-6.0233,-5.8181,-3.6094,1.7476,5.5068,6.5756,8.0416,6.3933,2.3550, 0.6078,-2.3342,-2.9824,-3.9807,-5.5271,-6.6924,-8.7267,-6.5221,-2.5583,2.1343,2.3062,4.1939,6.4870, 6.3125,3.2878,0.8703,-3.0262,-2.7133,-3.2428,-3.7807,-4.8137,-6.6618,-5.5921,-2.9025,1.1385,3.1125,3.7363,6.0362,6.7499,2.6324,0.0478] zs = [0,0,-0.5000,-0.8401,4.17,4.2583,6.9212,8.3677,1.20,-5.4182,-2.0932,-2.1863,-7.9830,-5.8500,-0.5991, -1.6810,2.3509,2.9533,0.3337,3.4925,5.9002,4.09,6.0108,2.8415,-5.80,-6.5369,-4.5087,-7.9595,-5.4608, 1.4428,0.2369,4.4268,2.3396,2.3834,4.0042,4.6696,8.9054,5.1745,0.7719,-5.4923,-1.4653,-3.1765,-7.2947,-6.4279,-0.0129,0.4961,5.86,1.4516,1.49,1.9010,3.0741,7.1232,5.2248,1.9393,-4.2395,-3.3718,-1.9258, 6.93,-7.3995,1.2763 ] zm = [0,0,-0.4400,-3.9913,-5.7014,-6.9415,-7.8178,-3.9097,1.4543,2.4075,3.5810,6.6598,6.0079,3.53,2.4376, -0.09,-2.3472,-2.3178,-4.4100,-6.9915,-6.0233,-5.8181,-3.6094,1.7476,5.5068,6.5756,8.0416,6.3933,2.3550, 0.6078,-2.3342,-2.9824,-3.9807,-5.5271,-6.6924,-8.7267,-6.5221,-2.5583,2.1343,2.3062,4.1939,6.4870,6.3125, 3.2878,0.8703,-3.0262,-2.7133,-3.2428,-3.7807,-4.8137,-6.6618,-5.5921,-2.9025,1.1385,3.1125,3.7363, 6.0362,6.7499,2.6324,0.0478] 图8.14 辨识系统在没有噪声情况下的实际输出响应 图8.15被辨识系统的采样输出响应 图8.16 模型含有噪声的输出响应 图8.17 模型去除噪声后的输出响应 仿真结果表明，递推到第9步时，参数辨识的结果基本达到稳定状态，即a1= -1.5000, a2= 0.7000, b1=1.0000, b2=0.5000, d1=1.0000, d2=-1.0000, d3=0.2000。此时，辨识参数的相对变化量E 0.000000005。与最小二乘递推算法相比，增广最小二乘递推算法虽然考虑了噪声模型，同样具有速度快、辨识结果精确的特点。 8.7 递推极大似然算法的产生程序 递推的极大似然法系统模型如图8.18所示。 u 其中v(k)为随机信号，输入信号为幅值为的M序列或随机信号。 ① 编程如下(见配套MATLAB程序algorithms7): 编程如下： clear %清零 a(1)=1;b(1)=0;c(1)=1;d(1)=0;u(1)=d(1);z(1)=0;z(2)=0; %初始化 for i=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); end u; %若取去‘；’可以在程序运行中观测到m序列 v=randn(1200,1); %产生正态分布随机数 V=0; %计算噪声方差 for i=1:1200 V=V+v(i)*v(i); end V1=V/1200; for k=3:1200 %根据v和u计算z z(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); end o1=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; %迭代计算参数值和误差值 for k=3:1200 h=[-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); 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)]; end o1 %若取去‘；’可以在程序运行中观测到参数 V1 %绘图 subplot(4,1,1) k=1:1200; plot(k,a1,'k:',k,a2,'b',k,b1,'r',k,b2,'m:',k,d1,'g',k,d2,'k'); xlabel('k') ylabel('parameter') legend('a1=-1.2,','a2=0.6','b1=1.0','b2=0.5','d1=-1.0','d2=0.2'); %图标注 title('The parameter idendification of the RML'); end subplot(4,1,2) k=1:1200; plot(k,e1,'k',k,e2,'b',k,e3,'r',k,e4,'m',k,e5,'g',k,e6,'k'); xlabel('k') ylabel('error') %title('误差曲线') end subplot(4,1,3) k=1:1200; plot(k,u); xlabel('k') ylabel('input') %title('系统输入信号') end subplot(4,1,4) k=1:1200; plot(k,v); xlabel('k') ylabel('random noise') %title('系统所加的随机噪声') end ② 程序运行结果如图8.19 所示 图8.19 RML辨识参数曲线 8.8 辅助变量算法的产生 辨识模型为 （8.11） 其中噪声为有色噪声，其产生过程为 （8.12） v(k) 为零均值的不相关随机噪声。参照 Tally 法选取辅助变量，为误差传递函数的阶数，此处为 2。则有 。 辅助变量法的递推公式可写成 初始条件，。 经过编译计算，各个参数的估计值为 MATLAB程序如下： %Tally 辅助变量最小二乘的递推算法，（见配套MATLAB程序algorithms8） %Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+e(k)，e(k)为有色噪声 %e(k)=v(k)+0.5*v(k-1)+0.2*v(k-2),v(k)为零均值的不相关随机噪声 %======================================== clear clc %==========400 个产生 M 序列作为输入=============== x=[0 1 0 1 1 0 1 1 1]; %initial value n=403; %n 为脉冲数目 M=[]; %存放 M 序列 for i=1:n temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end %===========产生均值为 0，方差为 1 的高斯白噪声========= v=randn(1,400); e=[]; e(1)=0.3; e(2)=0.5; for i=3:400 e(i)=v(i)+0.5*v(i-1)+0.2*v(i-2); end %==============产生观测序列 z================= z=zeros(402,1); z(1)=-1; z(2)=0; for i=3:400 z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+e(i); end P=100*eye(4); %估计方差 Pstore=zeros(4,400); Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; Theta=zeros(4,400); %参数的估计值，存放中间过程估值 Theta(:,1)=[3;3;3;3]; Theta(:,2)=[3;3;3;3]; Theta(:,3)=[3;3;3;3]; Theta(:,4)=[3;3;3;3]; % K=zeros(4,400); %增益矩阵 K=[10;10;10;10]; for i=5:400 h=[-z(i-1);-z(i-2);M(i-1);M(i-2)]; hstar=[-z(i-2-1);-z(i-2-2);M(i-1);M(i-2)]; %辅助变量 K=P*hstar*inv(h'*P*hstar+1); Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1)); P=(eye(4)-K*h')*P; Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; end %==================结果输出及作图=================== disp('参数 a1 a2 b1 b2 的估计结果：') Theta(:,400) i=1:400; figure(1) plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title('待估参数过渡过程') figure(2) plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title('估计方差变化过程') 待估参数过渡过程 图8.20 辅助变量法参数过渡过程 估计方差变化过程 图8.21 辅助变量法方差变化过程下载本文
