《系统辨识大作业汇总(共25页).doc》由会员分享，可在线阅读，更多相关《系统辨识大作业汇总(共25页).doc（25页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上一、 问题描述考虑仿真对象：式中，u(k)和z(k)是输入输出数据，v(k)是零均值、方差为1的不相关的随机噪声；u(k)采用与e(k)不相关的随机序列。1. 设计实验，产生输入输出数据；2. 使用基本最小二乘法估计参数；3. 考虑其他适用于有色噪声的辨识方法估计参数；4. 模型验证。二、 问题分析对于单输入单输出系统（Single Input Single Output, SISO），如图 1所示，将待辨识的系统看作“灰箱”，它只考虑系统的输入输出特性，而不强调系统的内部机理。图 1中，输入u(k)和输出z(k)是可以测量的，是系统模型，用来描述系统的输入输出特性，
2、y(k)是系统的实际输出。是噪声模型，v(k)是均值为零的不相关随机噪声，e(k)是有色噪声。图 1 SISO系统的“灰箱”结构对于SISO随机系统，被辨识模型为：其相应的差分方程为若考虑被辨识系统或观测信息中含有噪声，被辨识模型可改写为式中，z(k)为系统输出量的第k次观测值；y(k)为系统输出量的第k次真值，y(k-1)为系统输出量的第k-1次真值，以此类推；u(k)为系统的第k个输入值，u(k-1)为系统的第k-1个输入值；v(k)为均值为0的不相关随机噪声。1. 数据生成本部分需要生成系统的输入输出数据以及噪声数据。1） 白噪声的生成辨识数据通常包含有噪声，如果该噪声相关性较弱或者强度
3、很小，可近似看作白噪声。本次实验问题中，即服从标准正态分布，可以将噪声看作为服从正态分布的白噪声过程，在Matlab中可以由randn函数生成。2） 输入数据的生成伪随机二进制序列（Pseudo Random Binary Sequence, PRBS）是广泛应用的一种伪随机序列，所谓“二进制”是指序列中每个随机变量只有“0”或“1”两种逻辑状态。伪随机二进制序列可由多级线性反馈移位寄存器组成的随机信号发生器产生，其中具有最长循环周期的线性移位寄存器序列是伪随机二进制序列最常见的一种形式，简称M序列（Maximal Length Sequence）。M序列由于具有近似白噪声的性质，而且工程上易
4、于实现，能够保证较好的系统辨识精度，是普遍采用的一种辨识用输入信号。 图 2 线性反馈移位寄存器产生伪随机二进制序列结构图以一个4级线性反馈移位寄存器产生伪随机二进制序列为例，如图 2所示。假设4个移位寄存器输出的初态非全零，移位寄存器的工作原理是：一个移位脉冲来到后，每级移位寄存器的输出移到下一级移位寄存器作为输入，最末一级移位寄存器的输出即为伪随机二进制序列。3） 输出数据的生成根据给定的SISO系统，可以求出z(k)的表达式：其理想系数值为.可以根据生成的白噪声序列和输入序列，以及必要的0初始值，带入表达式即可得到采样输出数据。2. 差分模型阶检验在实际场景中，辨识模型的阶数和纯时延往往
5、是未知的，在很多情况下仅仅依靠猜测。在模型的阶数和纯时延不确定时，设系统模型为其中n为模型的阶数，模型的阶估计可以采用多种方法，本实验采用比较简单易行的损失函数检验法。定义预报误差（噪声方差的估计值）的平方和为损失函数，即当n从小增大时，应随之减小，当n增大到某一值时，应近似白噪声过程。采用以下的检验原则：在n-1这一点，最后一次出现陡峭的下降，此后就近似地保持不变或只有微小的下降，则取。3. 参数辨识模型在系统辨识和参数估计领域中，最广泛的估计方法时最小二乘法和极大似然估计法。最小二乘法作为一种最基本的估计方法应用极为广泛，其他的大多数估计算法都与最小二乘法有关。它既可用于动态系统，也可用于
6、静态系统；可用于线性系统，也可用于非线性系统；可用于离线估计，也可用于在线估计。在随机环境下利用最小二乘法时，无须知道观测数据的概率统计信息，而这种方法获得的结果，却有相当好的统计性质。最小二乘参数估计方法来源于数理统计的回归分析，它能提供一个在最小方差意义上与实验数据最好的你和的模型。该估计在一定条件下有最佳的统计特性，即它们是一致的，无偏的和有效的。最小二乘法时一个经典的方法，概念简明，适应范围广，在一些情况下，可得到与极大似然法一样好的统计效果，它能很方便地与其他辨识算法建立关系。（1）递推最小二乘算法当获得一批数据后，可一次求得相应的参数估计值，这样处理问题的方法称为一次完成算法或批处
7、理算法。它在理论研究方面有许多方便之处，但当矩阵的维数增加时，矩阵求逆运算的计算量会急剧增加，将给计算机的计算速度和存储量带来负担，而且不适合在线辨识，无法跟踪参数随时间变化的情况。为了减少计算量，减少数据在计算机中所占的存储量，也为了实时地辨识出动态系统的特性，在用最小二乘法进行参数估计时，把它转化成参数递推的估计。参数递推估计是指被辨识的系统，每取得一次新的测量数据后，就在前一次估计结果的基础上，利用新引入的测量数据对前一次估计的结果进行修正，从而递推地得出新的参数估计值。这样，随着新测量数据的引入，一次接一次地进行参数估计，直到估计值达到满意的精确程度为止。最小二乘递推算法的基本思想可以
8、概括为：当前的估计值=上次估计值+修正项即新的估计值是在旧的估计值的基础上，利用新的观测数据对旧的估计值进行修正而成的。考虑如下模型：，其中，分别是系统的输入和输出；是均值为零，方差为一的不相关白噪声。且满足：令则使用加权最小二乘参数估计递推算法（Recursive Weighted Least Squares, RWLS）有： 图 3 最小二乘递推辨识参数估计过程中信息的变换可以看出，取的时候，加权最小二乘估计就退化成了最小二乘参数估计的递推算法（Recursive Least Squares, RLS）。加权参数可以在范围内选择，如果选，所有的采样数据都是等同加权的，如果,则表示对新近获得
9、的数据给予充分大的加权因子，而削弱历史观测数据的影响。实际计算时，需要首先确定初始参数和。一般说来选取 对于这样的系统，使用最小二乘法参数估计的递推算法进行辨识可以得到无偏估计，但是如果噪声模型必须用表示时，此时就无法得到无偏估计了，因为该方法没有把噪声模型考虑进去。使用递推最小二乘法进行系统辨识的过程如图4所示。（2）增广最小二乘法当噪声均值为0时，最小二乘参数估计算法为无偏估计；当噪声的均值不为0时，最小二乘参数估计算法为有偏估计。为了解决最小二乘参数估计的有偏性，将噪声模型的辨识同时考虑进去，因此称为增广最小二乘法。该算法可以看成是对一般最小二乘参数辨识算法的简单推广或扩充，因此又称为扩
10、充最小二乘算法。考虑如下模型：，其中， 分别是系统的输入和输出；是均值为零，方差为一的不相关白噪声。且满足：在模型阶次已经确定的情况下，令由于是白噪声，所以利用增广最小二乘法这一形式上的变换，即可获得参数的无偏估计。不过由于数据向量中包含着不可测的噪声量，因此要用相应的估计值代替。当；当，即可得递推公式如下： 一般说来选取初值 增广最小二乘递推算法扩充了最小二乘法的参数向量和数据向量的维数，把噪声模型的辨识同时考虑进去。最小二乘法的许多结论对它都是适用的，但最小二乘法只能获得模型的参数估计，对于有色噪声，也就是噪声模型必须用表示时，只能采用增广最小二乘递推算法，方可获得无偏估计。这是最小二乘参
11、数估计的递推算法不可代替的。增广最小二乘算法的流程如图5所示。图4 递推最小二乘参数辨识流程图图5 增广最小二乘算法流程图三、Matlab仿真及运行结果在本次实验中，取值为1，则系统辨识对象为:1. 差分模型的阶检验在实际应用场景中，模型的阶数和纯时延往往是未知的，在此种情况下，一般采用如下方法来确定模型的阶数：固定一个阶数n，当n从小增大时，定义损失函数应随之减小，当n增大到某一值时，应近似白噪声过程。在n-1这一点，最后一次出现陡峭的下降，此后就近似地保持不变或只有微小的下降，则取。使用Matlab进行模型阶数确定代码如下：%噪声方差的估计值和损失函数法定阶%clearclose all;
12、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; %第四个移位寄存器的输入信号 y(i)=y4; %第四个移位寄存器的输出，即M序列，幅值为“0”或“1”if y(i)0.5;u(i)=-1; %M序列的值为1时，辨识的输入取-1else u(i)=1; %M序列为0时，输入取1 end y1=x1;y2=x2;y3=x3;y4=x4;%移位寄存器的输出序列endv=randn(L+5,1); %- ksi
13、(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) - 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=
14、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*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*estim
15、ate2)*(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); H3(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; % 模型阶
16、次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,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
17、,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(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 ont
18、itle(损失函数法判断模型阶次) xlabel(阶次) ylabel(损失函数) 程序运行结果如图6所示： 图6 模型阶次的判定由图6可知，当系统阶次从1到5变化时，在n=3及以后，损失函数只有微小的下降，在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:L x1=xor(y3,y4); x2=y1; x3=
19、y2; x4=y3; y(i)=y4; if y(i)0.5 u(i)=-u_f; else u(i)=u_f; end y1=x1;y2=x2;y3=x3;y4=x4;endv=1*randn(L,1);eb(2)=0;eb(1)=0;for k=3:L eb(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
20、(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.001 0.001;p0=103*eye(5,5);E= 0.;c=c0,zeros(5,L-1);e=zeros(5,L);%for k=4:L h1=-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;
21、 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时，辨识的输入为-1 else u(i)=1; %M序列的值为0时，辨识的输入为1 end y1=x1;y2=x2;y3=x3;y4=x4; %移位寄存器的输出序列,为下一次的输入信号做准备endfigure(1);subplot(3,1,1); stem(u),grid on %画出M序列输入信号径线图并给图形加上网格title(输入信号（M序列）) %-随机噪声的产生-%v=randn(1,L); %
22、产生一组L=60个正态分布的噪声subplot(3,1,2); plot(v),grid on %画出随机噪声信号title(随机噪声信号e(k) %-实际噪声的产生-%- 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);endsubplot(3,1,3)plot(ksi), grid ontitle(实际噪声信号xi(k) %-辨识过程-%z(3)=
23、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 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
24、*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),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; %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; %把递推出的辨识参数
25、c的列向量加入辨识参数矩阵 p1=p0-k1*k1*h1*p0*h1+1; % find p(k) p0=p1; %给下次的计算使用end %整个循环结束c,e, %显示被辨识参数及参数收敛情况% 分离变量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的收
26、敛情况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+,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 onleg
27、end(ea1,ea2,ea3,eb1,eb2,ec1,ec2,ec3)title(辨识误差曲线)figure(4);i=1:L;plot(i,z(i),g), grid on; title(辨识系统的输出响应) %画出被辨识系统的输出响应程序运行结果如图10至图13所示。 图10输入数据及噪声 图11 辨识曲线 图12辨识误差曲线图13辨识的输出响应使用增广最小二乘算法进行系统辨识的结果如表 2所示。表 2 增广最小二乘算法的辨识结果参数真值0.90.150.020.7-1.51.0-1.0-0.4估计值0.90.150.020.7-1.51.0-1.0-0.4从上表中可以看出，对于含有有色
28、噪声的系统，采用增广最小二乘法进行系统辨识可以取得无偏估计。四、模型与误差分析对系统进行辨识后，需要对辨识所得的系统进行进一步的模型验证，以观察所得的模型是否是对当前观测数据的最佳选择。对模型进行验证，可以采用事先生成的输入数据，分别输入到已知模型和经过系统辨识所得的模型中，检验所得模型是否存在问题。如果检验失败，可能存在的问题是：辨识所用的一组数据包含的信息部足或所选的模型类不合适。对本实验中的已知对象及经过系统辨识所得的模型，选取一组输入信号为：比较原参考对象的输出以及辨识模型的输出。具体代码如下： %模型检验%clear allclose all u=1,-1,1,-1,1,1,-1,1
29、,1,-1,1,1,-1,-1,-1,1;%系统输入信号L=20z=zeros(1,L);x=zeros(1,L);y=zeros(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);end for k=4:16 z(k)=-0.9*z(k-1)-0.15*z(k-2)-0.02*z(k-3)
30、+0.7*u(k-1)-1.5*u(k-2)+ksi(k) %真实模型的输出 y(k)=-1.2317*y(k-1)-0.4735*y(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,
31、z,r,i,y,b),grid onlegend(真实模型输出,递推最小二乘辨识模型)title(真实模型输出和递推最小二乘辨识模型输出值)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(真实模
32、型与增广最小二乘辨识模型的差值) 程序运行结果如图14和图15所示。 图14 不同辨识模型输出和真实系统输出对比 图15 真实模型与辨识模型输出信号的差值由上图并结合表1和表2，分析可得表 1 递推最小二乘算法的辨识结果参数真值0.90.150.020.7-1.5估计值0.93670.1591-0.01380.6275-1.3558表 2 增广最小二乘算法的辨识结果参数真值0.90.150.020.7-1.51.0-1.0-0.4估计值0.90.150.020.7-1.51.0-1.0-0.4辨识给定的含有色噪声的系统时，最小二乘递推辨识算法的参数辨识结果是有偏差的，无法达到稳定值,但是这个误
33、差是在-0.5,0.5范围内，还是可以接受的。而与最小二乘递推辨识算法相比，增广最小二乘递推辨识算法考虑了噪声模型，并且在递推到第九步时辨识结果就达到了稳定值，具有辨识速度快、辨识结果精确的特点。因此，对于含有有色噪声的系统，使用增广二乘辨识算法能够取得较好的效果。五、实验总结通过本次仿真实验，使我加深了对系统辨识算法的认识和了解，对于最小二乘参数辨识这种系统辨识领域中最基本最常用的方法有了较好的掌握，并通过MATLAB仿真熟悉了系统辨识的一般内容和步骤，对系统辨识的基本过程有了一个较为清晰的认识。最后，王立琦老师这一学期给我们所带的系统辨识这门课程使我受益良多，对我的科研工作有很大的启发借鉴意义。在此十分感谢王老师的指导和帮助。专心-专注-专业