1、一、 问题描述考虑仿真对象： ()0.9(1).5(2)0.(3)0.7(1).5(2)(zkzkzkukukee4e,vN式中，u(k)和 z(k)是输入输出数据，v(k)是零均值、方差为 1 的不相关的随机噪声；u(k)采用与 e(k)不相关的随机序列。1. 设计实验，产生输入输出数据；2. 使用基本最小二乘法估计参数；3. 考虑其他适用于有色噪声的辨识方法估计参数；4. 模型验证。二、 问题分析对于单输入单输出系统（Single Input Single Output, SISO），如图 1所示，将待辨识的系统看作“灰箱”，它只考虑系统的输入输出特性，而不强调系统的内部机理。图 1 中，
2、输入 u(k)和输出 z(k)是可以测量的， 是1()Gz系统模型，用来描述系统的输入输出特性，y(k)是系统的实际输出。 是N噪声模型，v(k)是均值为零的不相关随机噪声，e(k)是有色噪声。图 1 SISO 系统的“灰箱”结构对于 SISO 随机系统，被辨识模型 为：()Gz12() nbbzyzuaa 其相应的差分方程为 11()()()nni iykykbuk若考虑被辨识系统或观测信息中含有噪声，被辨识模型可改写为 11()()()(nni izkaykbukv式中，z(k)为系统输出量的第 k 次观测值；y(k)为系统输出量的第 k 次真值，y(k-1)为系统输出量的第 k-1 次真
3、值，以此类推；u(k)为系统的第 k 个输入值，u(k-1)为系统的第 k-1 个输入值；v(k)为均值为 0 的不相关随机噪声。1. 数据生成本部分需要生成系统的输入输出数据以及噪声数据。1）白噪声的生成辨识数据通常包含有噪声，如果该噪声相关性较弱或者强度很小，可近似看作白噪声。本次实验问题中 ，即服从标准正态分布，可以()0,1)vtN:将噪声看作为服从正态分布的白噪声过程，在 Matlab 中可以由 randn 函数生成。2）输入数据的生成伪随机二进制序列（Pseudo Random Binary Sequence, PRBS）是广泛应用的一种伪随机序列，所谓“二进制”是指序列中每个随机
4、变量只有“0”或“1”两种逻辑状态。伪随机二进制序列可由多级线性反馈移位寄存器组成的随机信号发生器产生，其中具有最长循环周期的线性移位寄存器序列是伪随机二进制序列最常见的一种形式，简称 M 序列（Maximal Length Sequence）。M 序列由于具有近似白噪声的性质，而且工程上易于实现，能够保证较好的系统辨识精度，是普遍采用的一种辨识用输入信号。图 2 线性反馈移位寄存器产生伪随机二进制序列结构图以一个 4 级线性反馈移位寄存器产生伪随机二进制序列为例，如图 2所示。假设 4 个移位寄存器 输出的初态非全零，移位寄存器的工0123,a作原理是：一个移位脉冲来到后，每级移位寄存器的输
5、出移到下一级移位寄存器作为输入，最末一级移位寄存器的输出即为伪随机二进制序列。3）输出数据的生成根据给定的 SISO 系统，可以求出 z(k)的表达式：()0.9(1)0.5(2)0.(3)0.7(1).152v142zkzkzkzkukuee其理想系数值为.12312120.9,.,0.,.7,.5,.0,41aabc可以根据生成的白噪声序列和输入序列，以及必要的 0 初始值，带入表达式即可得到采样输出数据。2. 差分模型阶检验在实际场景中，辨识模型的阶数和纯时延往往是未知的，在很多情况下仅仅依靠猜测。在模型的阶数和纯时延不确定时，设系统模型为 11()()()(nni jytaytbyti
6、t其中 n 为模型的阶数， tCzet模型的阶估计可以采用多种方法，本实验采用比较简单易行的损失函数检验法。定义预报误差（噪声方差的估计值）的平方和为损失函数，即 1()NiJt当 n 从小增大时， 应随之减小，当 n 增大到某一值时， 应近似白噪声NJ NJ过程。采用以下的检验原则：在 n-1 这一点， 最后一次出现陡峭的下降，此NJ后就近似地保持不变或只有微小的下降，则取 。3. 参数辨识模型在系统辨识和参数估计领域中，最广泛的估计方法时最小二乘法和极大似然估计法。最小二乘法作为一种最基本的估计方法应用极为广泛，其他的大多数估计算法都与最小二乘法有关。它既可用于动态系统，也可用于静态系统；
7、可用于线性系统，也可用于非线性系统；可用于离线估计，也可用于在线估计。在随机环境下利用最小二乘法时，无须知道观测数据的概率统计信息，而这种方法获得的结果，却有相当好的统计性质。最小二乘参数估计方法来源于数理统计的回归分析，它能提供一个在最小方差意义上与实验数据最好的你和的模型。该估计在一定条件下有最佳的统计特性，即它们是一致的，无偏的和有效的。最小二乘法时一个经典的方法，概念简明，适应范围广，在一些情况下，可得到与极大似然法一样好的统计效果，它能很方便地与其他辨识算法建立关系。（1）递推最小二乘算法当获得一批数据后，可一次求得相应的参数估计值，这样处理问题的方法称为一次完成算法或批处理算法。它
8、在理论研究方面有许多方便之处，但当矩阵的维数增加时，矩阵求逆运算的计算量会急剧增加，将给计算机的计算速度和存储量带来负担，而且不适合在线辨识，无法跟踪参数随时间变化的情况。为了减少计算量，减少数据在计算机中所占的存储量，也为了实时地辨识出动态系统的特性，在用最小二乘法进行参数估计时，把它转化成参数递推的估计。参数递推估计是指被辨识的系统，每取得一次新的测量数据后，就在前一次估计结果的基础上，利用新引入的测量数据对前一次估计的结果进行修正，从而递推地得出新的参数估计值。这样，随着新测量数据的引入，一次接一次地进行参数估计，直到估计值达到满意的精确程度为止。最小二乘递推算法的基本思想可以概括为：当
9、前的估计值 =上次估计值 +修正项()k(1)k即新的估计值 是在旧的估计值 的基础上，利用新的观测数据对旧的()k1估计值进行修正而成的。考虑如下模型： ，其中 ， 分别是系统11()()()AzkBzukv()ukz的输入和输出； 是均值为零，方差为一的不相关白噪声。且满足：v112().nmzazazBbb令 1212()(),.(),.().Tnmhkzzkukab则使用加权最小二乘参数估计递推算法（Recursive Weighted Least Squares, RWLS）有：1 ()1)()()1()()()()1,.,.TT TkKkzhkKPhpIPhkzzknukm 图 3
10、 最小二乘递推辨识参数估计过程中信息的变换可以看出，取 的时候，加权最小二乘估计就退化成了最小二乘参数()1k估计的递推算法（Recursive Least Squares, RLS）。加权参数 可以在1范围内选择，如果选 ，所有的采样数据都是等同加权的，如果(0,1 1()k,则表示对新近获得的数据给予充分大的加权因子，而削弱历史观测()k数据的影响。实际计算时，需要首先确定初始参数 和 。(0)P一般说来选取 2(0)PI为 充 分 大 实 数为 充 分 小 的 向 量对于这样的系统，使用最小二乘法参数估计的递推算法进行辨识可以得到无偏估计，但是如果噪声模型必须用 表示时，此时就无法得到无
11、偏估1()Czvk计了，因为该方法没有把噪声模型考虑进去。使用递推最小二乘法进行系统辨识的过程如图 4 所示。（2）增广最小二乘法当噪声均值为 0 时，最小二乘参数估计算法为无偏估计；当噪声的均值不为 0 时，最小二乘参数估计算法为有偏估计。为了解决最小二乘参数估计的有偏性，将噪声模型的辨识同时考虑进去，因此称为增广最小二乘法。该算法可以看成是对一般最小二乘参数辨识算法的简单推广或扩充，因此又称为扩充最小二乘算法。考虑如下模型： ，其中 ， 分别111()()()AzkBzukCzvk()ukz是系统的输入和输出； 是均值为零，方差为一的不相关白噪声。且满足：v11212()nmvzazazb
12、bCcc在模型阶次 已经确定的情况下，令,nmv121212()(),.(),.(),(,.(). .TTnmvhkzzkukkvabc由于 是白噪声，所以利用增广最小二乘法这一形式上的变换，即可获()v得参数 的无偏估计。不过由于数据向量 中包含着不可测的噪声量()hk，因此要用相应的估计值代替。当 ；当(1),2,.()vkvk 0,()kv，即可得递推公式如下：0, 1Tzh1()()()()()(1),.,1,.,(),.()T TkKkvKPhPhIzhkzknukmvvk 一般说来选取初值 (0)PI增广最小二乘递推算法扩充了最小二乘法的参数向量 和数据向量 的()hk维数，把噪声
13、模型的辨识同时考虑进去。最小二乘法的许多结论对它都是适用的，但最小二乘法只能获得模型的参数估计，对于有色噪声，也就是噪声模型必须用 表示时，只能采用增广最小二乘递推算法，方可获得无偏估计。1()Czvk这是最小二乘参数估计的递推算法不可代替的。增广最小二乘算法的流程如图 5 所示。图 4 递推最小二乘参数辨识流程图 图 5 增广最小二乘算法流程图三、Matlab 仿真及运行结果在本次实验中， 取值为 1，则系统辨识对象为:()0.9().5(2)0.(3)0.7(1).5(2)v4zkzzkzkukkee1. 差分模型的阶检验在实际应用场景中，模型的阶数和纯时延往往是未知的，在此种情况下，一般
14、采用如下方法来确定模型的阶数：固定一个阶数 n，当 n 从小增大时，定义损失函数 应随之减小，当 n 增大到某一值时， 应近似白噪声过程。在NJ NJn-1 这一点， 最后一次出现陡峭的下降，此后就近似地保持不变或只有微小的下降，则取 。使用 Matlab 进行模型阶数确定代码如下：n%噪声方差的估计值和损失函数法定阶%clearclose all;L=100;y1=1;y2=1;y3=1;y4=0;for i=1:L+5;x1=xor(y3,y4); %第一个移位寄存器的输入信号x2=y1; %第二个移位寄存器的输入信号x3=y2; %第三个移位寄存器的输入信号x4=y3; %第四个移位寄存
15、器的输入信号y(i)=y4; %第四个移位寄存器的输出，即 M 序列，幅值为“0”或“1”if y(i)0.5;u(i)=-1; %M 序列的值为 1 时，辨识的输入取-1else u(i)=1; %M 序列为 0 时，输入取 1endy1=x1;y2=x2;y3=x3;y4=x4;%移位寄存器的输出序列endv=randn(L+5,1); %- ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)-%ksi(2) = 0;ksi(1) = 0; %设 ksi 的前两个初始值为 0lamda = 1;for k=3:L+5;ksi(k) = lamda*v(k) -
16、 ksi(k-1) - 0.41*ksi(k-2);endz(3)=0;z(2)=0;z(1)=0; %取 z 的前三个初始值为零for k=4:L+5; z(k) = -0.9*z(k-1) - 0.15*z(k-2) - 0.02*z(k-3) + 0.7*u(k-1) - 1.5*u(k-2) + ksi(k);%理想辨识输出采样信号end% 模型阶次 n=1 for i=1:L H1(i,1)=z(i); H1(i,2)=u(i); end estimate1=inv(H1*H1)*H1*z(2:L+1); D1=(z(2:L+1)-H1*estimate1)*(z(2:L+1)-H1
17、*estimate1)/L; %噪声方差的估计值AIC1=L*log(D1)+4*1; % 模型阶次 n=2 for i=1:L H2(i,1)=z(i+1); H2(i,2)=z(i); H2(i,3)=u(i+1); H2(i,4)=u(i); end estimate2=inv(H2*H2)*H2*z(3:L+2); D2=(z(3:L+2)-H2*estimate2)*(z(3:L+2)-H2*estimate2)/L; %噪声方差的估计值AIC2=L*log(D2)+4*2; % 模型阶次 n=3 for i=1:L H3(i,1)=z(i+2); H3(i,2)=z(i+1); H
18、3(i,3)=z(i); H3(i,4)=u(i+2); H3(i,5)=u(i+1); H3(i,6)=u(i); end estimate3=inv(H3*H3)*H3*z(4:L+3); D3=(z(4:L+3)-H3*estimate3)*(z(4:L+3)-H3*estimate3)/L; %噪声方差的估计值AIC3=L*log(D3)+4*3; % 模型阶次 n=4 for i=1:L H4(i,1)=z(i+3); H4(i,2)=z(i+2); H4(i,3)=z(i+1); H4(i,4)=z(i); H4(i,5)=u(i+3); H4(i,6)=u(i+2); H4(i,
19、7)=u(i+1); H4(i,8)=u(i); end estimate4=inv(H4*H4)*H4*z(5:L+4); D4=(z(5:L+4)-H4*estimate4)*(z(5:L+4)-H4*estimate4)/L; %噪声方差的估计值AIC4=L*log(D4)+4*4; % 模型阶次 n=5 for i=1:L H5(i,1)=z(i+4); H5(i,2)=z(i+3); H5(i,3)=z(i+2); H5(i,4)=z(i+1);H5(i,5)=z(i);H5(i,6)=u(i+4); H5(i,7)=u(i+3); H5(i,8)=u(i+2); H5(i,9)=u
20、(i+1);H5(i,10)=u(i);end estimate5=inv(H5*H5)*H5*z(6:L+5); D5=(z(6:L+5)-H5*estimate5)*(z(6:L+5)-H5*estimate5)/L; %噪声方差的估计值AIC5=L*log(D5)+4*5; plot(1:5,AIC1 AIC2 AIC3 AIC4 AIC5) grid ontitle(损失函数法判断模型阶次) xlabel(阶次) ylabel(损失函数) 程序运行结果如图 6 所示：图 6 模型阶次的判定由图 6 可知，当系统阶次从 1 到 5 变化时，在 n=3 及以后，损失函数只有微小的下降，在
21、n=4 以后又出现上升的趋势，所以取 n=3 作为模型的阶次。2. 参数辨识（1）递推最小二乘法根据前文所述递推最小二乘法算法流程，编写递推最小二乘法程序如下：% 递推最小二乘辨识clcclear all%采用 M 序列方法产生输入数据L=90;%90 个数据namenda=1;y1=1;y2=1;y3=1;y4=0;u_f=2;for i=1:Lx1=xor(y3,y4); x2=y1; x3=y2; x4=y3; y(i)=y4;if y(i)0.5u(i)=-u_f;elseu(i)=u_f;endy1=x1;y2=x2;y3=x3;y4=x4;endv=1*randn(L,1);eb(
22、2)=0;eb(1)=0;for k=3:Leb(k)=-eb(k-1)-0.41*eb(k-2)+namenda*v(k);endz(3)=0;z(2)=0;z(1)=0;for k=4:L;z(k)=-0.9*z(k-1)-0.15*z(k-2)-0.02*z(k-3)+0.7*u(k-1)-1.5*u(k-2)+eb(k);endfigure(1)i=1:30;subplot(2,1,1)stem(u(i);title(输入 M 序列)grid onsubplot(2,1,2)stem(z(i);title(输出序列)grid on;% c0=0.001 0.001 0.001 0.00
23、1 0.001;p0=103*eye(5,5);E= 0.0000005;c=c0,zeros(5,L-1);e=zeros(5,L);%for k=4:Lh1=-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2);k1=p0*h1*inv(h1*p0*h1+1);c1=c0+k1*(z(k)-h1*c0);e1=c1-c0;e2=e1./c0;e(:,k)=e2;c0=c1;c(:,k)=c1;p1(:,:,k)=p0-k1*k1*h1*p0*h1+1;p0=p1(:,:,k);if abs(e2)0.5;u(i)=-1; %M 序列的值为 1 时，辨识的输入为-1el
24、seu(i)=1; %M 序列的值为 0 时，辨识的输入为 1endy1=x1;y2=x2;y3=x3;y4=x4; %移位寄存器的输出序列,为下一次的输入信号做准备endfigure(1);subplot(3,1,1); stem(u),grid on %画出 M 序列输入信号径线图并给图形加上网格title(输入信号（M 序列）)%-随机噪声的产生-%v=randn(1,L); %产生一组 L=60 个正态分布的噪声subplot(3,1,2); plot(v),grid on %画出随机噪声信号title(随机噪声信号 e(k)%-实际噪声的产生-%- ksi(k)=lamda*v(k)
25、-ksi(k-1)-0.41*ksi(k-2)-%ksi(2) = 0;ksi(1) = 0; %设 ksi 的前两个初始值为 0lamda = 1;for k=3:L;ksi(k) = lamda*v(k) - ksi(k-1) - 0.41*ksi(k-2);endsubplot(3,1,3)plot(ksi), grid ontitle(实际噪声信号xi(k)%-辨识过程-%z(3)=0;z(2)=0;z(1)=0;zm(3)=0;zm(2)=0;zm(1)=0; %输出采样，不考虑噪声时系统输出，不考虑噪声时模型输出，模型输出的初值c0=0.001 0.001 0.001 0.001
26、0.001 0.001 0.001 0.001; %直接给出辨识被辨识参数的初始值，即一个充分小的实向量p0=106*eye(8,8); %直接给出初始向量 P0，即一个充分大的实数单位矩阵c=c0,zeros(8,L-1); %被辨识参数参数矩阵的初始值及大小e=zeros(8,L); %相对误差的初始值及大小for k=4:L; %开始求 Kz(k) = -0.9*z(k-1) - 0.15*z(k-2) - 0.02*z(k-3) + 0.7*u(k-1) - 1.5*u(k-2) + ksi(k);%系统在 M 序列下的输出采样信号h1=-z(k-1),-z(k-2),-z(k-3),
27、u(k-1),u(k-2),v(k),ksi(k-1),ksi(k-2); %为求 K（k）做准备x=h1*p0*h1+1;x1=inv(x);k1=p0*h1*x1; %Kd1=z(k)-h1*c0;c1=c0+k1*d1; %辨识参数 ce1=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; %给下次的计算使用end %整个循环结束c,e, %显示被辨识参数及参数收敛情况% 分离变量
28、a1=c(1,:);a2=c(2,:);a3=c(3,:);b1=c(4,:);b2=c(5,:); %分离出 a1,a2,a3,b1,b2d1=c(6,:);d2=c(7,:);d3=c(8,:); %分离出的 d1,d2,d3% 分离变量的收敛情况ea1=e(1,:);ea2=e(2,:);ea3=e(3,:);eb1=e(4,:);eb2=e(5,:); %分离出 a1,a2,b1,b2 的收敛情况ed1=e(6,:);ed2=e(7,:);ed3=e(8,:); %分离出 d1,d2,d3 的收敛情况figure(2); i=1:L;plot(i,a1,r,i,a2,r:,i,a3,r
29、+,i,b1,b,i,b2,b:,i,d1,g,i,d2,g:,i,d3,g+)grid onlegend(a1,a2,a3,b1,b2,c1,c2,c3)title(辨识曲线)figure(3); i=1:L;plot(i,ea1,r,i,ea2,r:,i,ea3,r+,i,eb1,b,i,eb2,b:,i,ed1,g,i,ed2,g:,i,ed3,g+) %画出各个参数的收敛情况grid onlegend(ea1,ea2,ea3,eb1,eb2,ec1,ec2,ec3)title(辨识误差曲线)figure(4);i=1:L;plot(i,z(i),g), grid on; title(辨
30、识系统的输出响应) %画出被辨识系统的输出响应程序运行结果如图 10 至图 13 所示。图 10 输入数据及噪声图 11 辨识曲线图 12 辨识误差曲线图 13 辨识的输出响应使用增广最小二乘算法进行系统辨识的结果如表 2 所示。表 2 增广最小二乘算法的辨识结果参数真值 0.9 0.15 0.02 0.7 -1.5 1.0 -1.0 -0.4估计值 0.9 0.15 0.02 0.7 -1.5 1.0 -1.0 -0.4从上表中可以看出，对于含有有色噪声的系统，采用增广最小二乘法进行系统辨识可以取得无偏估计。四、模型与误差分析对系统进行辨识后，需要对辨识所得的系统进行进一步的模型验证，以观察
31、所得的模型是否是对当前观测数据的最佳选择。对模型进行验证，可以采用事先生成的输入数据，分别输入到已知模型和经过系统辨识所得的模型中，检验所得模型是否存在问题。如果检验失败，可能存在的问题是：辨识所用的一组数据包含的信息部足或所选的模型类不合适。对本实验中的已知对象及经过系统辨识所得的模型，选取一组输入信号为： 1,1,1,u比较原参考对象的输出以及辨识模型的输出。具体代码如下：%模型检验%clear allclose allu=1,-1,1,-1,1,1,-1,1,1,-1,1,1,-1,-1,-1,1;%系统输入信号L=20z=zeros(1,L);x=zeros(1,L);y=zeros(
32、1,L);v=randn(1,L);%-实际噪声的产生-%- ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)-%ksi(2) = 0;ksi(1) = 0; %设 ksi 的前两个初始值为 0lamda = 1;for k=3:L;ksi(k) = lamda*v(k) - ksi(k-1) - 0.41*ksi(k-2);endfor k=4:16 z(k)=-0.9*z(k-1)-0.15*z(k-2)-0.02*z(k-3)+0.7*u(k-1)-1.5*u(k-2)+ksi(k) %真实模型的输出y(k)=-1.2317*y(k-1)-0.4735*y
33、(k-2)-0.0600*y(k-3)+0.8534*u(k-1)-1.2373*u(k-2)%递推最小二乘辨识模型输出x(k)=-0.9*x(k-1)-0.15*x(k-2)-0.02*x(k-3)+0.7*u(k-1)-1.5*u(k-2)+ksi(k) %增广最小二乘辨识模型的输出endfigure(1);subplot(3,1,1) stem(u) ,grid ontitle(系统输入信号)subplot(3,1,2)i=1:1:L; plot(i,z,r,i,y,b),grid onlegend(真实模型输出,递推最小二乘辨识模型)title(真实模型输出和递推最小二乘辨识模型输出值
34、)subplot(3,1,3)i=1:1:L; plot(i,z,r,i,x,g),grid onlegend(真实模型输出,增广最小二乘辨识模型)title(真实模型输出和增广最小二乘辨识模型输出值)figure(2);subplot(2,1,1)i=1:1:L; plot(i,z-y), grid on;title(真实模型与递推最小二乘辨识模型的差值)subplot(2,1,2)i=1:1:Lplot(i,z-x), grid on;title(真实模型与增广最小二乘辨识模型的差值)程序运行结果如图 14 和图 15 所示。图 14 不同辨识模型输出和真实系统输出对比图 15 真实模型与
35、辨识模型输出信号的差值由上图并结合表 1 和表 2，分析可得表 1 递推最小二乘算法的辨识结果参数真值 0.9 0.15 0.02 0.7 -1.5估计值 0.9367 0.1591 -0.0138 0.6275 -1.3558表 2 增广最小二乘算法的辨识结果参数真值 0.9 0.15 0.02 0.7 -1.5 1.0 -1.0 -0.4估计值 0.9 0.15 0.02 0.7 -1.5 1.0 -1.0 -0.4辨识给定的含有色噪声的系统时，最小二乘递推辨识算法的参数辨识结果是有偏差的，无法达到稳定值,但是这个误差是在-0.5,0.5范围内，还是可以接受的。而与最小二乘递推辨识算法相比，增广最小二乘递推辨识算法考虑了噪声模型，并且在递推到第九步时辨识结果就达到了稳定值，具有辨识速度快、辨识结果精确的特点。因此，对于含有有色噪声的系统，使用增广二乘辨识算法能够取得较好的效果。五、实验总结通过本次仿真实验，使我加深了对系统辨识算法的认识和了解，对于最小二乘参数辨识这种系统辨识领域中最基本最常用的方法有了较好的掌握，并通过 MATLAB 仿真熟悉了系统辨识的一般内容和步骤，对系统辨识的基本过程有了一个较为清晰的认识。最后，王立琦老师这一学期给我们所带的系统辨识这门课程使我受益良多，对我的科研工作有很大的启发借鉴意义。在此十分感谢王老师的指导和帮助。
