系统辨识大作业

系统辨识大作业——递推增广最小二乘法题目一：已知一系统为两输入单输出系统，观测数据受有色噪声污染，噪信比为N/S=0.1。系统经2000次采样，存放于文件T3T.TXT中。系统输入u1为7级M序列，u2为u1的63步移位序列。模型类可选为：A(q-1)y(k)=B1(q)u1(k)+B2(q)u2(k)+w(k)/C(q-1)。要求编制程序，辨识出该模型的结构及参数。（注：可以将模型变形，以适合算法）作业文档要求：描述问题；选择辨识方法并简单说明所选方法中的结构辨识原理和参数估计原理；程序流程图及程序清单；说明程序中用到的一些技术，如数据标准化、UD分解、稳定性判断等；结构搜索路线及各结构下的参数、残差；给出最终结果：A(q-1)=B1(q)=B2(q)=C(q-1)=并给出选择此最终结果的理由；用你的辨识结果来预报系统输出误差e(k)=y(k)-y(k)，并画出e(101)-e(400)的曲线图。1.问题描述已知一系统为两输入单输出系统，观测数据受有色噪声污染，噪信比为N/S=0.1。系统经2000次采样，存放于文件T3T.TXT中。系统输入u1为7级M序列，u2为u1的63步移位序列。模型类可选为：A(q-1)y(k)=B(q-1)u1(k)+C(q-1)u2(k)+D(q-1)w(k)把题目所选模型A(q-1)y(k)=B1(q)u1(k)+B2(q)u2(k)+w(k)/C(q-1)转化为问题描述的模型。2.辨识原理和参数估计原理2.1递推增广最小二乘法的辨识原理递推增广最小二乘法(简称RELS方法)是用于控制系统参数估计和结构辩识的一种方法,它基于岱(最小二乘)方法,使被控对象数学模型在误差信号平方和最小的意义上由实验数据拟合,同时把有色噪声看成是由白噪声合成的,从而解决了最小二乘法算法的有偏性和非一致性问题。其中nc=0；nd任意时是广义递推最小二乘法就是递推增广最小二乘法。2.3结构辨识在系统辨识中，系统的阶次是一个十分重要的参数，但在实际的情况中往往无法提前获知系统的阶数。这就需要在参数辨识之前先通过结构辨识，来确定阶数。一般来说经典的结构辨识方法主要有以下3种：F检验定阶法，AIC准则法，FPE准则法。（1）F检验定阶法假设系统的阶次为n（1,2,3...n），计算出模型输出和观测输出的残差平方和函数Jˆˆ()()TJzzzz(18)其中z为观测输出，ˆz为模型输出。一般来说J会随着模型阶次的增加而减少。实际进行系统辨识时，输入输出的信号采样的个数大大超过待辨识的参数个数，在这种情况下，随着模型阶次的增加J值先是显著下降，当模型的阶次大于真实系统的阶次时，J的显著下降现象就中止。利用这个原理可以判定模型应有的阶次。（2）AIC准则法赤池信息量准则，即Akaikeinformationcriterion、简称AIC，是衡量统计模型拟合优良性的一种标准，是由日本统计学家赤池弘次创立和发展的。赤池信息量准则建立在熵的概念基础上，可以权衡所估计模型的复杂度和此模型拟合数据的优良性。在一般的情况下，AIC可以表示为：22lnAICkL(19)其中：k是参数的数量，L是似然函数。假设条件是模型的误差服从独立正态分布。设n为观察数，RSS为剩余平方和，那么AIC变为：2ln(/)AICkRSSn增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个。赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。（3）FPE准则法对AR模型的一种定阶法.设tx是平稳AR(p)序列，12,,...Txxx是样本，ˆk是样本自协方差函数，ˆj(j=1，2，…，n)是拟合AR(n)参数的最小二乘估计.101ˆˆ()(1)(1)()niiinnFPEnTT（20）为最终预报误差.使FPE(n)达到最小的n称为AR(p)模型阶数p的估计.此方法是赤池弘次于1969年提出的所谓改进的残差方差图方法.最终预报误差准则名称的来源是对AR(p)序列观测值拟合k阶AR模型，使得一步预报均方误差达到最小，只有k取真阶p才能达到.所以定阶依据这种准则是合乎实际目的的。3程序流程图3.1结构辨识结果在参数辨识之前首先进行结构辨识，使用F检验法。取最后200个数据计算残差平方和函数值。具体数值如表1所示表1F检验法nJt（n-1，n）26.639830.8680645.005340.2363256.636550.23450.7292系统模型为111112()()()()()()()()AqzkBqukCqukDqvk。模型阶次每增加1，待辨识参数数目就增加4。查询F分布表可知(4,)2.3719F。观察表1中t的值，当阶次为4时，2.3719t，当阶次为5时，2.3719t。根据F检验法的思想，在模型阶次为5时，残差平方函数值显著下降的现象中止，所以实际系统的阶次为5。3.2程序流程图及程序清单模型完成定阶以后，就开始参数的辨识。由于系统的阶次为5，所以模型的递推表达式可以写成5554121110()()()()()iiiiiiiizkazkibukicukidvki（21）仿真程序的流程图如图1所示。图1递推输出误差辨识算法的程序流程图递推增广最小二乘法的辨识程序如下：%递推增广最小二乘法clearL=2000;n=5;size=5+4*(n-1);清零产生有色噪声信号画出有色噪声序列图加载T3t.txt文件并分别赋给z，u1，u2给辨识参数θ和P赋初值按照式(17)的第二式计算K(k)按照式(17)的第一式计算θ(k)按照式(17)的第三式计算P(k)计算被辨识参数的相对变化量参数是否满足要求？分离参数画出被辨识参数θ的各次递推估计式图形画出被辨识参数θ的相对误差的图形停机NYv=sqrt(0.1)*randn(2000,1);M=load('T3t.txt');z(:,1)=M(1:L,1);u1(:,1)=M(:,2);u2(:,1)=M(:,3);%递推增广最小二乘法c0=0.0001*ones(size,1);%直接给出被辨识参数的初始值,即一个充分小的实向量p0=1000000*eye(size,size);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.00005;%取相对误差Ec=[c0,zeros(size,L-1)];%被辨识参数矩阵的初始值及大小e=zeros(size,L);%相对误差的初始值及大小zm=zeros(L,1);ek=zeros(L,1);fork=6:L;%开始求Kh1=[-z(k-1,1),-z(k-2,1),-z(k-3,1),-z(k-4,1),-z(k-5,1),u1(k-1,1),u1(k-2,1),u1(k-3,1),u1(k-4,1),u1(k-5,1),u2(k-1,1),u2(k-2,1),u2(k-3,1),u2(k-4,1),u2(k-5,1),v(k),v(k-1),v(k-2),v(k-3),v(k-4),v(k-5)]';%为求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];%p1=p0-k1*h1'*p0;p0=p1;%给下次用end%循环结束fork=6:Lh2=[-zm(k-1,1),-zm(k-2,1),-zm(k-3,1),-zm(k-4,1),-zm(k-5,1),u1(k-1,1),u1(k-2,1),u1(k-3,1),u1(k-4,1),u1(k-5,1),u2(k-1,1),u2(k-2,1),u2(k-3,1),u2(k-4,1),u2(k-5,1),v(k),v(k-1),v(k-2),v(k-3),v(k-4),v(k-5)]';zm(k)=h2'*c(:,L);ek(k,1)=z(k,1)-zm(k);endJ=ek'*ek;J,c(:,L)%显示被辨识参数及参数收敛情况figure(1);k=1:L;plot(k,z,'r',k,zm,'b');xlabel('k');ylabel('幅值');title('观测输出与模型输出');figure(2);k=101:400;plot(k,ek(101:400,1),'b');xlabel('k');ylabel('幅值');title('残差');a1=c(1,:);a2=c(2,:);a3=c(3,:);a4=c(4,:);a5=c(5,:);figure(3);subplot(3,1,1);k=1:L;plot(k,a1(1,1:L),'b');xlabel('k');ylabel('幅值');title('a1');subplot(3,1,2);k=1:L;plot(k,a2(1,1:L),'b');xlabel('k');ylabel('幅值');title('a2');subplot(3,1,3);k=1:L;plot(k,a3(1,1:L),'b');xlabel('k');ylabel('幅值');title('a3');figure(4);subplot(2,1,1);k=1:L;plot(k,a4(1,1:L),'b');xlabel('k');ylabel('幅值');title('a4');subplot(2,1,2);k=1:L;plot(k,a5(1,1:L),'b');xlabel('k');ylabel('幅值');title('a5');b1=c(6,:);b2=c(7,:);b3=c(8,:);b4=c(9,:);b5=c(10,:);figure(5);subplot(3,1,1);k=1:L;plot(k,b1(1,1:L),'b');xlabel('k');ylabel('幅值');title('b1');subplot(3,1,2);k=1:L;plot(k,b2(1,1:L),'b');xlabel('k');ylabel('幅值');title('b2');subplot(3,1,3);k=1:L;plot(k,b3(1,1:L),'b');xlabel('k');ylabel('幅值');title('b3');figure(6);subplot(2,1,1);k=1:L;plot(k,b4(1,1:L),'b');xlabel('k');ylabel('幅值');title('b4');subplot(2,1,2);k=1:L;plot(k,b5(1,1:L),'b');xlabel('k');ylabel('幅值');title('b5');c1=c(11,:);c2=c(12,:);c3=c(13,:);c4=c(14,:);c5=c(15,:);figure(7);subplot(3,1,1);k=1:L;plot(k,c1(1,1:L),'b');xlabel('k');ylabel('幅值');title('c1');subplot(3,1,2);k=1:L;plot(k,c2(1,1:L),'b');xlabel('k');ylabel('幅值');title('c2');subplot(3,1,3);k=1:L;plot(k,c3(1,1:L),'b');xlabel('k');ylabel('幅值');title('c3');figure(8);subplot(2,1,1);k=1:L;plot(k,c4(1,1:L),'b');xlabel('k');ylabel('幅值');title(