系统辨识

系 统 辨 识 作 业系统辨识作业：•已知某系统为单输入/单输出系统，其测量噪声为有色噪声，分布未知。

现给出一个实验样本(如下表所示)，求该系统模型。

说明：可采用GLS ，ELS ，IV 等，要定阶，要比较仅用RLS 的计算结果 一、问题分析在估计模型参数时需要已知模型的阶数，但是由于本系统模型阶数也是未知的，所以本系统需要先由输入/输出数据通过辩识得出系统的阶数。

然后根据辨识的系统阶数再分析求解系统模型。

二、模型阶数的辨识按照品质指标“残差平方总和”定阶，如高阶系统模型相应的系数为零，则可退化成相应的低阶系统即低阶模型可视为高阶模型的特例。

理论上高阶模型的精度不低于低阶模型，但是考虑到计算机的舍入误差的影响，过高的阶数亦能引起模型精度的下降。

一般说低阶模型描述粗糙，高阶模型精度高，但是代价亦大。

根据逼近的观点，定阶往往是考虑多种因素的折衷。

定阶一般是按照假设——检验的步骤进行的，检验过程中往往带有主观成分。

一般说来低阶模型描述粗糙，高阶模型精度高。

残差平方总和J(n)是模型阶数的函数在不同的模型阶数的假设下，参数估计得到的J(n)值亦不同。

定阶的最简单办法是直接用J(n)。

设模型阶数的“真值”为n 0 ,当n < n 0 时随着n 的增加，J(n)值将明显的下降；而当n ≥ n 0 时随着n 的增加，J(n)值变化将不显著。

因此，由J(n)曲线随着n 的增加最后一次陡峭下降的n 值定做n 的估计值。

用数理统计的检验方法，判断n 的增加使得J(n)值改善是否明显。

讨论如下(1).当n=1时程序如下: clearu=zeros(100,1);%构造输入矩阵 z=zeros(100,1);%构造输出矩阵u=[-0.93249 0.34935 0.76165 -0.9964 -0.38894 -0.12288 0.021565 -0.49555 -0.61624 -1.912 0.22207 -0.31231 -0.17866 -1.8356 -0.26472 1.7642 -1.0418 1.1146 -2.0856 0.8152 1.5094 -0.5822 0.61097 0.35521 2.5907 1.5843 -0.9603 -0.27341 0.39947 0.17493 -1.7451 0.8112 1.2645 1.5682 0.63959 -0.47757 0.99697 0.058774 -0.16174 -1.2928 -0.04722 0.73182 -0.19644 0.091783 -1.1908 -0.90716 0.85388 0.33836 0.74074 0.54181 0.15676 -0.50569 -0.17521 1.3255 -2.488 0.50261 -1.1533 0.36407 0.65283 -0.05983∑=-=Nk T Kk y n J 12))(()(θϕ-1.1464 1.1406 1.3891 0.75736 -0.23474 0.1793 0.084435 0.35464 -0.21055 -1.0089 -2.5115 -1.3007 0.91626 1.9674 0.95495 0.499 0.36871 -0.12931 -1.572 -1.0032 1.0072 -0.20873 0.93346 -0.88953 -0.60081 -0.62593 1.4748 -1.2291 0.064272 -1.2139 0.15342 0.11332 0.23502 -1.0776 -0.32697 0.2859 0.67131 1.7297 0.70529 -1.2088]';%输入数据uz=[ 0.28412 -1.2053 -1.8813 -1.7496 -2.8686 -4.9332 -4.6691 -2.8167 -1.2759 -0.7668 -2.6371 -4.317 -6.3424 -7.5086 -7.8543 -7.7153 -3.1131 1.303 5.385 7.1704 6.7143 6.8642 6.8652 6.5386 4.9673 5.1895 5.4094 3.3101 -0.74817 -4.1961 -5.4271 -8.6193 -9.6902 -7.9423 -2.1012 4.2807 8.403 10.151 11.106 9.4447 6.4218 3.5482 2.3237 1.856 1.443 -0.24757 -4.6216 -6.7792 -4.6945 0.31967 4.1559 8.5023 9.8421 9.3436 9.3422 7.1959 3.2164 -1.2355 -4.1895 -4.4498 -4.2321 -4.5896 -4.4376 -1.9601 1.6745 4.572 6.75 8.2541 7.0644 3.8969 -0.72395 -7.4809 -13.613 -16.961 -16.283 -11.596 -4.368 3.5309 8.0703 7.7467 5.2576 4.7738 6.5092 6.221 2.3859 -4.1897 -9.2999 -9.2728 -9.1932 -10.059 -9.7339 -6.3149 -1.5968 3.5503 7.4666 9.883 11.627 11.234 10.671 10.235 ]';%输出数据zr=100;for p=1:(r-2) %利用循环生成观测矩阵h(p,:)=[-z(p+1) u(p+1)]; %endhl=h;for b=1:(r-2) %生成输出矩阵zl(b,:)=[z(b+2)];zl'endzl'%根据最小二乘法公式进行参数辩识c1=hl'*hl;c2=inv(c1);c3=hl'*zl;c=c2*c3;a1=c(1)a2=c(2)j=0;for k=4:100;hl=[-z(k-1);u(k-1)]';x=hl*c;y=z(k)-x;s=y*y;j=j+s;endj仿真结果如下a1 = -0.9280 a2 =1.0351 j = 755.1949(2)当 n=2时程序如下(输入输出数据同上,只给出不同于一阶系统的程序不同之处)其中U、Z分别是作业要求给出得的输入输出，数据输入同上。

r=100;%利用循环生成观测矩阵。

for p=1:(r-2)h(p,:)=[-z(p+1) -z(p) u(p+1) u(p)];endhl=h;%生成输出矩阵。

for b=1:(r-2)zl(b,:)=[z(b+2)];zl'endzl'%根据最小二乘法公式进行参数辩识c1=hl'*hl;c2=inv(c1);c3=hl'*zl;c=c2*c3;%Q ls%输出辩识参数a1=c(1)a2=c(2)b1=c(3)b2=c(4)j=0;%求J（n）for k=4:100; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';x=h1*c;y=z(k)-x;s=y*y；j=j+s;endj仿真结果如下a1 = -1.6585 b1 = 0.9656 a2 = 0.8266 b2 = 0.3972 j=98.9904(3)当n=3时程序如下(只给出与一阶系统的不同之处)r=100;for p=2:(r-2)h(p,:)=[-z(p+1) -z(p) -z(p-1) u(p+1) u(p) u(p-1)];Endhl=h;for b=2:(r-2)zl(b,:)=[z(b+2)];zl'endzl'c1=hl'*hl;c2=inv(c1);c3=hl'*zl;c=c2*c3;a1=c(1)a2=c(2)a3=c(3)b1=c(4)b2=c(5)b3=c(6)j=0;for k=4:100;hl=[-z(k-1);-z(k-2);-z(k-3);u(k-1);u(k-2);u(k-3)]';x=hl*c;y=z(k)-x;s=y*y;j=j+s;endj仿真结果如下a1 =-1.9030 a2 =1.2235 a3 = -0.1934 b1 = 0.9392 b2 = 0.1779 b3 = -0.3438 j = 87.6641数据分析如下利用LS法先对系统参数进行初步辩识并根据其确定残差平方总和J(n) 结果如下：阶数辩识n a i b i J (n) 1-0.9280 1.0351 755.1949 2 -1.6585 0.8266 0.96560.3972 98.99043 -1.9030 1.2235 -0.1934 0.93920.1779 -0.343887.6641根据公式，我们可以计算出相应的J （N ）变化率：T （1）=6.62897 T （2）=0.1292由此我可以看出，在由一阶到二阶的T （1）很大，说明J （N ）有一个很明显的变化，而T （2）很小，说明随着阶数的增加J （N ）的变化很小，据此我们有理由相信该系统应该是二阶的。

三、模型参数的辩识 1）、由上面的过程我们得到该系统是二阶的，又考虑到有色噪声的影响，所以采用递推增广最小二乘法（RELS ）进行辩识。

思路如下： ·考虑 CARMA 模型A (z -1) y (k) =B (z -1) u (k) +C (z -1) ε (k)其中：A (z -1) = 1 + a 1 z - 1 +…+ a n z -nB (z -1) = b 1 z - 1 +…+ b n z - nC (z -1) = 1 + c 1 z - 1 +…+ c r z - r还可表示为 ：以上两向量均由 2n 维扩展为 (2n+r) 维。

在式Kϕ中若ε ( k-1)、、、ε ( k-r) 是已知量，则可直接用 LS 法估计出 θ ，但是{ ε ( k) }是不可测的未知量。

一个简单可行的方法是用计算的ε ( k -1) 代替ε ( k -1) ,…. ε ( k -r) , { ε ( k) }由下式递推得出：ε (k) = y(k) - ϕk θ k -1 (计算残差) 其中：ϕk = [ -y(k-1),..,-y(k-n) , u(k-1),..,u(k-n) , ε (k-1),.., ε (k-r) ]T递推算式)1()1()()(++-=n J n J n J n T[][]1)2(111,...,,,...,,,...,)(),..,1(),(),...,1(),(),...,1(⨯+=--------=r n Tr n n TKc c b b a a r k k n k u k u n k y k y θεεϕ)()(k k y T kεθϕ+=θ N+1 = θ N + K N+1(y(N+1) – ϕN+1T θ N )ε (N+1) = y(N) - ϕN+1 θ NϕN = [-y(N-1),..,-y(N-n),u(N-1),,u(N-n),ε(N-1),..,ε (N-r) ]T 递推初值：θ 0= 0 ; P 0= 10 6 I; ε (0) = ε (-1) = … = ε (1-r) = 0用matlab 编程实现如下（仅列出程序中除输入输出部分输入输出同上一个程序） c0=[0.001 0.001 0.001 0.001 0.001 0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(6,6);%直接给出初始状态P0，即一个充分大的实数单位矩阵 E=5.0e-15;%取相对误差Ec=[c0,zeros(6,99)];%被辨识参数矩阵的初始值及大小 e=zeros(6,100);%相对误差的初始值及大小 v=zeros(1,100);for k=3:100; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2),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 e1=c1-c0;e2=e1./c0; %求参数的相对变化 e(:,k)=e2;c0=c1;%给下一次用v(k)=z(k)-h1'*c0;%预报噪声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, %显示被辨识参数及参数收敛情况a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);c1=c(5,:);c2=c(6,:); i=1:100;plot(i,a1,'r',i,a2,'r:',i,b1,'g',i,b2,'g:',i,c1,'k',i,c2,'k:') 仿真结果如下所示11111+++++=N N T N N N N P P Kϕϕϕ111111++++++-=N N TN NT N N N N NP P P P P ϕϕϕϕ参数收敛图形如下：图形分析:由图象可以看出在n 小于30时系统辨识参数误差比较大,但在n 大于20时随着n 的增大2）、当仅用RLS 做，程序实现如下：（输入输出部分同上）c0=[0.001 0.001 0.001 0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵 c=[c0,zeros(4,99)];%被辨识参数矩阵的初始值及大小 e=zeros(4,100);%相对误差的初始值及大小 for k=3:100; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(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; %求参数的相对变化 c0=c1;%给下一次用c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵二阶a 1= -1.5691b 1= 1.0045c 1= -0.4163 a 2= 0.7705 b 2= 0.5056c 2= 0.1563p1=p0-k1*k1'*[h1'*p0*h1+1];%find p(k) p0=p1;%给下次用 end%判断结束 end%循环结束c, e, %显示被辨识参数及参数收敛情况 a1=c(1,:) a2=c(2,:) b1=c(3,:) b2=c(4,:) i=1:100;plot(i,a1,'r',i,a2,'r:',i,b1,'g',i,b2,'g:') 仿真结果如下：参数收敛图形如下：三、结论：综上两种辨识做法，我们很容易得看出当n 小于10时，两种方法辨识的结果误差都很大，但两种方法对系统参数辨识的结果与理论值有着明显差别，当n 大于50时两种辨识方法所辨识的模型系统参数都以不同的收敛速度收敛于理论二阶a 1= -1.6585b 1= 0.9656 a 2= 0.8266b 2= 0.3971值，随着n的继续增大两种方法辨识结果的差别也越来越小，但n在20与50之间从图中可以看出递推增广最小二乘法的辨识结果明显好于标准递推最小二乘法。