系统辨识大作业报告

1、系 统 辨 识 大 作 业 报 告 班 级 13202-4 学 号 姓 名 亓子龙 报告日期 2013.12.03 一、辨识方法1.最小二乘法及其递推方法(1)最小二乘法：构造阵，利用公式计算；(2)递推最小二乘法：取前20个数据，利用基本最小二乘法给出和的初值和，然后利用公式：迭代计算。2.辅助变量法及其递推方法(1) 辅助变量法：首先利用基本最小二乘估计作为计算的初值，利用计算结果构造阵，依据公式迭代计算直至收敛。(2) 递推辅助变量法：前50个数据利用递推最小二乘估计和初值和，然后依据递推公式进行迭代计算：3.广义最小二乘法及其递推方法(1)广义最小二乘法：使用和按基本最小二乘求出估计值

2、，计算残差用残差代替并计算的估值，的计算公式为。利用进行数据滤波，再按最小二乘法重新估计，重复这些步骤直至收敛。(2)递推广义最小二乘法：前20个数据利用基本最小二乘求得递推的初值，然后按以下公式计算，根据所得结果对新的观测值进行数据滤波，然后进行的递推，递推公式如下：4.夏氏偏差修正法、夏氏改良法及递推夏氏法(1)夏氏偏差修正法：首先计算并作为的初值，然后计算残差并构造矩阵，同时计算矩阵和。其次由计算和，其中。利用公式获得的估计值，循环迭代直到基本保持不变。(2)夏氏改良法：基本同夏氏偏差修正法，只是的计算变为，减小了计算量。(3)夏氏递推算法：利用基本最小二乘获得递推初值和，递推公式如下：

3、5.增广矩阵法利用基本最小二乘构造阵，进而计算递推初值和。构造向量其中，然后根据递推方程递推计算：6.极大似然法 (1)选定初值，、由LS方法获得，可任意指定；(2)计算残差及指标；(3)计算梯度及海赛(Hassian)矩阵；(4)计算的新估值；(5)返回第(2)步直到收敛。7.辨识结果及噪声特性分析对于数据uy1.txt，辨识结果如表1：表1 uy1辨识结果辨识结果辨识方法基本最小二乘1.48550.78690.48370.1982递推最小二乘1.48550.78690.48370.1982辅助变量法1.49640.75880.48780.1910递推辅助变量法1.53530.83320.4

4、8730.2156广义最小二乘法1.50220.79980.48440.2034递推广义最小二乘法1.42570.26460.41290.1062夏氏修正法1.48860.79730.48220.2022夏氏改良法1.48860.79730.48220.2022夏氏递推法1.48250.79450.47970.2018增广矩阵法1.48360.78640.48280.1991-0.00060.0074极大似然法1.46300.78610.47930.19560.23820.3530对于数据uy2.txt，辨识结果如表2：表2 uy2辨识结果辨识结果辨识方法基本最小二乘1.11130.49630

5、.37910.1879递推最小二乘1.11130.49630.37910.1879辅助变量法1.38750.80360.36490.3289递推辅助变量法1.39650.70720.38020.2970广义最小二乘法1.23820.58840.37770.2368递推广义最小二乘法1.67000.90220.17310.0060夏氏修正法1.38340.71460.38840.2873夏氏改良法1.38340.71460.38840.2873夏氏递推法1.35740.69540.39980.2719增广矩阵法1.16740.55790.38620.21330.08340.0353极大似然法1.

6、81960.97110.38880.44781.14790.0574对于数据uy3.txt，辨识结果如表3：表3 uy3辨识结果辨识结果辨识方法基本最小二乘1.11580.48010.42540.1245递推最小二乘1.11580.48010.42540.1245辅助变量法4.34872.56280.53471.1997递推辅助变量法1.61720.86870.43250.3138广义最小二乘法1.27630.61170.42570.1898递推广义最小二乘法1.48920.61730.05710.5818夏氏修正法1.46280.77490.47110.2252夏氏改良法1.46280.77

7、490.47110.2252夏氏递推法1.30880.62960.46420.1864增广矩阵法1.25640.61390.44810.16990.19980.1012极大似然法1.65590.97700.41750.38791.06620.4991三组数据输入值都一样，输出值不同。对上面表1-表3的辨识结果进行分析，可以看出，uy1的噪声近似为白噪声，因为其他辨识方法的辨识结果和最小二乘法的结果很接近。而uy2和uy3的噪声是有色噪声，其他辨识方法的结果和最小二乘法的结果出入较大。二、辨识方法比较分析1.分析比较各种方法估计的精度1） 最小二乘法是成批处理观测数据，即离线辨识，在输入为白噪声

8、的情况下其辨识精度是最高的，但是在有色噪声的情况下偏差比较大，且估计的均方差随噪声均方差的增大而增大。2） 递推最小二乘法是在线辨识。理论上讲，其辨识精度应等于采用离线辨识的最小二乘法，但是由于在递推被辨识参数的初值中，P的取值相当大时，递推最小二乘法的结果很接近于最小二乘法的结果。其基本思想可以概括成：新的估计值=旧的估计值+修整项。它是一种实时控制算法。3） 辅助变量法在输入为有色噪声的情况下能克服上述两种方法的有偏估计缺点。但本题前一部分的输入是白噪声，因而辅助变量法的辨识精度差于上述两种方法。在输入为有色噪声的情况下，辨识结果较好。在计算时需构造辅助变量矩阵。4） 广义最小二乘法：如果

9、输入是有色噪声，则广义最小二乘法能克服估计的有偏性。但是广义最小二乘法是一种迭代方法，且收敛速度比较慢。在系统的噪声较大时，最小二乘法的指标J可能是多峰的，因此该算法未必收敛于真实参数，本题中，广义最小二乘法在白噪声情况下，其辨识的效果略差于最小二乘法和递推最小二乘法，尔在有色噪声情况下，其辨识的效果好于那两种方法。实际估计效果较好，工程上广泛应用。缺点是计算量大。5） 夏氏法也是为了克服输入为有色噪声而导致辨识的有偏性而提出的，改善广义最小二乘法的计算速度，提高计算效率。它还可以应用到多输入多输出系统。对于本题，它的优缺点基本同广义最小二乘法，其辨识结果同广义最小二乘法差不多。6） 增广矩阵

10、法：它是无偏估计，收敛性好，系统参数与噪声参数同时辨识。由于考虑了系统噪声的影响，采用递推方法，其精度和收敛性比递推最小二乘法有提高。7） 极大似然法以观测值的出现概率最大作为标准，为离线辨识方法，在输入为白噪声的情况下，结果的精度与最小二乘法不相上下。在输入为有色噪声的情况下，精度比最小二乘法能好。2.分析其计算量1） 最小二乘法：算法最简单，但是当输入输出数据量较大时，会出现维数较高的矩阵求逆和相乘。 所以数据量越大计算量将明显增加。2） 递推最小二乘法：递推最小二乘法在每一次迭代过程中都不需要矩阵求逆，只需要做矩阵或向量乘法，故每次迭代过程中，递推最小二乘法的计算量很小。3） 辅助变量法

11、：每一次迭代都需要求辅助模型的输出变量，并由此构造求辅助变量矩阵Z，做一次(2n+1) *(2n+1)的矩阵求逆，故计算量是最小二乘法的计算量整数倍。4） 广义最小二乘法：每一次迭代都需要利用上一次估计出的被辨识参数来求出残差，利用残差来构造(2n+1)*(2n+1)的残差阵 ，然后求得新的输入和输出序列，再一次利用最小二乘法重新估计被辨识参数 ，故每次迭代都需要利用两次最小二乘法，另外还要作数据滤波，故计算量远远大过辅助变量法。5） 夏氏法：每一次迭代过程中求残差矩阵，不需要数据的反复滤波，故计算量远小于广义最小二乘法，但明显大于辅助变量法的计算量。6） 增广矩阵法：增广矩阵法的估计参数中增

12、加了噪声的模型参数，故相应的矩阵及向量都比递推最小二乘法大，计算量也比递推最小二乘法大的多。7） 极大似然法: 极大似然法的计算量较大，但其参数估计量具有良好的渐进性质。3.分析噪声方差的影响由各种辨识方法可以看出，当噪声的方差增大时，估计的偏差明显增大，各项估计值得均方差也稍有增大。当噪声的方差为0，即系统输入输出数据不含噪声时，估计结果就是系统参数的真值。4.比较白噪声和有色噪声对辨识的影响噪声为白噪声时比噪声为有色噪声所引起的误差要小，且当噪声为白噪声时，最小二乘法估计具有无偏性和一致性。三、系统模型阶次的辨识1.用三种方法确定系统的阶次并辨识(1)按残差方差定阶：利用基本最小二乘法估计

13、各阶次时的， 按公式求得残差，然后计算。对某一系统，当时，随着的增大而减小。假定系统为阶，则在时出现最后一次陡峭的下降，从而判定系统阶次。(2)确定阶的Akaike信息准则：这个准则给出了一个定量的判断标准，它不要求建模人员主观的判断“陡峭的下降”。在一组可供选择的随机模型中，AIC最小的那个模型是一个可取的模型。白噪声情况下的AIC定阶公式为；有色噪声情况下的AIC定阶公式为。(3) 按残差白色定阶：计算残差的估值的自相关函数，检查其白色性，即可判断模型的阶次。残差的自相关函数为将规格化可得。（1）按残差方差和F检验法定阶（2）AIC定阶 （3）残差白色定阶 （4）辨识结果：残差方差F检验残

14、差白色数据uy1uy2uy3uy1uy2uy3uy1uy2uy3阶数222222222AIC白噪声有色噪声数据uy1uy2uy3uy1uy2uy3阶数2332322.噪声对定阶的影响由于噪声的存在，使得各种定阶方法以及参数辨识方法所得的结果与真值有所差异。当定阶模型中采用的是白噪声模型时，应对测试数据进行滤波处理，辨识结果的好坏很大程度上取决于有色噪声的白化程度，这点可从三组测试数据滤波前后的辨识效果中看出。3.三种方法的优劣及有效性残差方差定阶法：按估计误差方差最小来确定模型的阶。优点是算法较为简单；缺点是最后一次陡降依赖个人的主观判断；F检验和AIC定阶法：优点在于它是一个完全客观的模型。

15、残差白色定阶：这种方法通过检验残差的白色性来验证模型的估计是否合适。优点是概念清晰，算法简单；缺点是残差的白色性判别缺乏客观标准，存在人为因数。四、附录：Matlab程序1.最小二乘法%最小二乘法LSclear all;clc;format short;load uy1.txt;load uy2.txt;load uy3.txt;u=uy1(:,1),uy2(:,1),uy3(:,1);y=uy1(:,2),uy2(:,2),uy3(:,2);n=2;%阶次N=length(u)-n;%数据长度for i = 1:3 %三个系统一起进行辨识 Y = y(n+1:n+N, i); %输出向量 f

16、ai = -y(n:n+N-1, i), -y(1:N, i), u(n:n+N-1,i), u(1:N,i); %测量矩阵 sita = inv(fai * fai) * fai * Y; %参数向量 sitaLS( :, i) = sita;enddisp 最小二乘法三组数据辨识结果；disp(sitaLS);2.递推最小二乘法%递推最小二乘法RLSclear all%clcformat shortload uy1.txt;load uy2.txt;load uy3.txt;u=uy1( :,1),uy2( :,1),uy3( :,1);y=uy1( :,2),uy2( :,2),uy3(

17、 :,2);n = 2; %阶次 sita = zeros(4 , 1);sitaRLS = zeros(4, 3);for i = 1:3 p = (1e6)n * eye(n*n); for N = 0:length(u) - n - 1 psai = -y(n+N, i), -y(N+1, i), u(n+N,i), u(N+1,i); k = p*psai*inv(1 + psai*p*psai); p = p - p*psai*inv(1 + psai*p*psai)*psai*p; sita = sita + k*(y(n+N+1, i) - psai*sita); end sit

18、aRLS( :, i) = sita;enddisp 递推最小二乘法三组辨识结果数据;disp(sitaRLS);3.辅助变量法%辅助变量法IVLSclear all%clcformat shortload uy1.txt;load uy2.txt;load uy3.txt;u=uy1(:,1);y=uy2(:,2);y1=uy2(3:500,2);for i=3:500 Psi(i-2,:)=-y(i-1) -y(i-2) u(i-1) u(i-2);endsita=inv(Psi*Psi)*Psi*y1;Z=Psi;k=1;while k5 y_1=Z*sita; Z(1,:)=y(2)

19、y(1) u(2) u(1); Z(2,:)=-y_1(1) y(2) u(3) u(2); for i=5:500 Z(i-2,:)=-y_1(i-3) -y_1(i-4) u(i-1) u(i-2); end sitaIV=inv(Z*Psi)*Z*y1; if abs(sitaIV-sita)/sita)0.0001 break; end sita=sitaIV; k=k+1;enddisp(辅助变量法辨识结果：);disp(sita);4.递推辅助变量法%递推辅助变量法clear all%clcformat shortload uy1.txt;load uy2.txt;load uy3

20、.txt;n=2;L=length(uy1);N=L-n;U=uy1(:,1);Y=uy3(:,2);Y1=uy3(3:L,2);OL=-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2); S1=OL*OL;S2=inv(S1);S3=OL*Y1;sita=S2*S3;OF=OL;Fsita=sita,zeros(4,N-1);for i=1:1:100 YM=OF*sita; Z=OF; Z(2:N,1)=-YM(1:(N-1),1); Z(3:N,2)=-YM(1:(N-2),1); sita=inv(Z*OF)*Z*Y1; Fsita(:,i+1)=sita;en

21、ddisp(递推辅助变量法辨识结果：)disp(sita);5.广义最小二乘法%广义最小二乘法GYLSclear all;clc;format short;load uy1.txt;load uy2.txt;load uy3.txt;u = uy1( :, 1);y = uy1( :, 2), uy2( :, 2), uy3( :, 2);n = 2; %阶次%最小二乘法LSN = length(u) - n; %数据长度sitaLS = zeros(4, 3);for i = 1:3 %三个系统一起进行辨识 Y = y(n+1:n+N, i); %输出向量 fai = -y(n:n+N-1,

22、 i), -y(1:N, i), u(n:n+N-1), u(1:N); %测量矩阵 sita = inv(fai * fai) * fai * Y; %参数向量 sitaLS( :, i) = sita;end%广义最小二乘法GYLSsitaGYLS = zeros(4, 3);sita0 = zeros(4, 1); f = zeros(2, 1); Y = zeros(n+N, 1); U = zeros(n+N, 1);for i = 1:3 sita0 = sitaLS( :, i); fai = -y(n:n+N-1, i), -y(1:N, i), u(n:n+N-1), u(1

23、:N); while 1 e = zeros(N+n, 1); %残差 e(n+1:n+N) = y(n+1:n+N, i) - fai*sita0; omega = -e(n:n+N-1), -e(n-1:n+N-2); f = inv(omega * omega) * omega * e(n+1:n+N); Y(n+1:n+N) = y(n+1:n+N, i) + f(1)*y(n:N+1, i) + f(2)*y(n-1:N, i); U(n+1:n+N) = u(n+1:n+N) + f(1)*u(n:N+1) + f(2)*u(n-1:N); %使用最小二乘法再次估算sita fai

24、 = -Y(n:n+N-1), -Y(1:N), U(n:n+N-1), U(1:N); sita = inv(fai * fai) * fai * Y(n+1:n+N); if (abs(sita0(1)- sita(1) 50) P2 = (1E6)n *eye(n*n); e(n+1:n+N) = y(n+1:n+N, i) - fai*sita; omega = -e(n:n+N-1), -e(n-1:n+N-2); f = inv(omega *omega) *omega * e(n+1:n+N); Y(n+1:n+N) = y(n+1:n+N,i) + f(1)*y(n:N+1,i

25、) + f(2)*y(n-1:N,i); U(n+1:n+N) = u(n+1:n+N) + f(1)*u(n:N+1) + f(2)*u(n-1:N); psi2 = -Y(n+N-1), -Y(N), U(n+N-1), U(N); K2 = P2*psi2*inv(1 + psi2*P2*psi2); P2 = P2 - P2*psi2*inv(1 + psi2*P2*psi2)*psi2*P2; sita = sita + K2*(Y(n+N) - psi2*sita); end end sitaRGYLS( :, i) = sita;enddisp 递推广义最小二乘法三组数据辨识结果

26、：;disp(sitaRGYLS);7.夏氏偏差修正法%夏式修正法XSCclear all;clc;format short;load uy1.txt;load uy2.txt;load uy3.txt;u=uy1(:,1);y=uy1(:,2),uy2(:,2),uy3(:,2);n=2;%阶次N = length(u) - n;sita = zeros(4, 1);sitaB0 = zeros(4, 1);sitaB = zeros(4, 1);sitaXSC = zeros(4, 3);omega = zeros(N, n);for i = 1:3 e = zeros(N+n, 1);

27、Y = y(n+1:n+N, i); fai = -y(n:n+N-1, i), -y(1:N, i), u(n:n+N-1), u(1:N); T = inv(fai*fai)*fai; M = eye(N) - fai*T; sitaLS( :, i) = T * Y; while (1) e(n+1:n+N) = Y - fai*sita; omega( :, 1) = -e(n:n+N-1); omega( :, 2) = -e(n-1:n+N-2); D = omega*M*omega; f = inv(D)*omega*M*Y; sitaB = T*omega*f; if (abs

28、(sitaB(1)- sitaB0(1)1e-6) break; end sitaB0 = sitaB; sita = sitaLS( :, i) - sitaB; end sitaXSC( :, i) = sita;enddisp 夏氏修正法三组数据辨识结果：;disp(sitaXSC);8.夏氏改良法%夏氏改良法XSIclear all;clc;format short;load uy1.txt;load uy2.txt;load uy3.txt;u=uy1(:,1);y=uy1(:,2),uy2(:,2),uy3(:,2);n=2;%阶次N=length(u)-n;%数据长度sita =

29、 zeros(4, 1);sita0 = zeros(4, 1);sitaXSI = zeros(4, 3);omega = zeros(N, n);e = zeros(N+n, 1);for i = 1:3 Y = y(n+1:n+N, i); fai = -y(n:n+N-1, i), -y(1:N, i), u(n:n+N-1), u(1:N); T = inv(fai*fai)*fai; sitaLS( :, i) = T * Y; while (1) e(n+1:n+N) = Y - fai*sita; E = e(n+1:n+N); omega( :, 1) = -e(n:n+N-

30、1); omega( :, 2) = -e(n-1:n+N-2); f = inv(omega*omega)*omega*E; sita = sitaLS( :, i) - T*omega*f; if (abs(sita(1)- sita0(1)1e-6) break; end sita0 = sita; end sitaXSI( :, i) = sita;enddisp 夏氏改良法三组辨识结果数据;disp(sitaXSI);9.递推夏氏法%递推夏氏算法RXSclear all;%clc;format short;load uy1.txt;load uy2.txt;load uy3.txt;

31、u=uy1(:,1);y=uy1(:,2),uy2(:,2),uy3(:,2);n=2;%阶次N=length(u)-n;%数据长度sitaRXS = zeros(4, 3);sita = zeros(4, 1);for i = 1:3 Y = y(n+1:n+N, i); %输出向量 fai = -y(n:n+N-1, i), -y(1:N, i), u(n:n+N-1), u(1:N); %测量矩阵 sita = inv(fai * fai) * fai * Y; %参数向量 sitaLS( :, i) = sita; P = (1e6)n * eye(3*n); beta = sitaL

32、S( :, i); 0; 0; for N = 1:length(u)-n e = zeros(n+N, 1); fai = -y(n:n+N-1, i), -y(1:N, i), u(n:n+N-1), u(1:N); e(n+1:n+N) = y(n+1:n+N, i) - fai*sita; psi = -y(n+N-1, i), -y(N, i), u(n+N-1), u(N), -e(n+N-1), -e(n+N-2); r = P*psi*inv(1+psi*P*psi); P = P - r*psi*P; beta = beta + r*(y(n+N, i) - psi*beta

33、); sita = beta(1:4); end sitaRXS( :, i) = sita;enddisp 夏氏递推法三组数据辨识结果：;disp(sitaRXS);10.增广矩阵法%增广矩阵法AMclear all;%clc;format short;load uy1.txt;load uy2.txt;load uy3.txt;u=uy1(:,1);y=uy1(:,2),uy2(:,2),uy3(:,2);n=2;%阶次N=length(u)-n;%数据长度sita = zeros(6 , 1);sitaAM = zeros(6, 3);e = zeros(N+n, 1);for i =

34、1:3 P = (1e6)n * eye(3*n); e(1) = y(1, i); e(2) = y(2, i); for N = 0:length(u) - n -1 psi = -y(n+N, i), -y(N+1, i), u(n+N), u(N+1), e(n+N), e(N+1); e(n+N+1) = y(n+N+1, i) - psi*sita; K = P*psi*inv(1 + psi*P*psi); P = P - K*psi*P; sita = sita + K*(y(n+N+1, i) - psi*sita); end sitaAM( :, i) = sita;end

35、disp 增广矩阵法三组数据辨识结果：;disp(sitaAM);11.极大似然法%极大似然辨识系统参数法MLclear all;clc;format short;load uy1.txt;load uy2.txt;load uy3.txt;u = uy1( :, 1);y = uy1( :, 2), uy2( :, 2), uy3( :, 2);n = 2; %阶次%最小二乘法LSN = length(u) - n; %数据长度sitaLS = zeros(4, 3);for i = 1:3 %三个系统一起进行辨识 Y = y(n+1:n+N, i); %输出向量 fai = -y(n:n+

36、N-1, i), -y(1:N, i), u(n:n+N-1), u(1:N); %测量矩阵 sita = inv(fai * fai) * fai * Y; %参数向量 sitaLS( :, i) = sita;endfor i = 1:3 sigma20 = 0.1; e = zeros(n+N, 1); dea1 = zeros(n+N, 1); dea2 = zeros(n+N, 1); deb1 = zeros(n+N, 1); deb2 = zeros(n+N, 1); dec1 = zeros(n+N, 1); dec2 = zeros(n+N, 1); sita = sitaL

37、S( :, i); 0.2; 0.2; fai = -y(n:n+N-1, i), -y(1:N, i), u(n:n+N-1), u(1:N); e(n+1:n+N) = y(n+1:n+N, i) - fai*sita(1:4); while (1) J = 0; dJ = 0; d2J = 0; sigma2 = 0; for k = n+1:n+N sigma2 = sigma2 + e(k)2/N; J = J + 0.5*e(k)2; dea1(k) = y(k-1, i) - sita(5)*dea1(k-1) - sita(6)*dea1(k-2); dea2(k) = y(k

38、-2, i) - sita(5)*dea2(k-1) - sita(6)*dea2(k-2); deb1(k) = -u(k-1) - sita(5)*deb1(k-1) - sita(6)*deb2(k-1); deb2(k) = -u(k-2) - sita(5)*deb2(k-1) - sita(6)*deb2(k-2); dec1(k) = -e(k-1) - sita(5)*dec1(k-1) - sita(6)*dec1(k-2); dec2(k) = -e(k-2) - sita(5)*dec2(k-1) - sita(6)*dec2(k-2); de = dea1(k), de

39、a2(k), deb1(k), deb2(k), dec1(k), dec2(k); dJ = dJ + e(k)*de; d2J = d2J + de*de; end sita = sita - inv(d2J)*dJ; if (sigma22 - sigma202)/sigma202)0.0001) break; end sigma20 = sigma2; end sitaML(:, i) = sita;enddisp 极大似然法三组数据辨识结果：;disp(sitaML);12.按残差方差定阶和残差白色定阶%残差方差定阶clear all;clc;format short;load uy

40、1.txt;load uy2.txt;load uy3.txt;u=uy1(:,1);y=uy3(:,2);Jie = 10;for n = 1:Jie N = length(u) - n; Y(1:N,1) = y(n+1:length(u); for i = 1:N for j = 1:n fai(i,j) = -y(n+i-j); fai(i,n+j) = u(n+i-j); end end sita = inv(fai*fai)*fai*Y; Err = Y - fai*sita; J(n) = Err*Err; R0=J(n)/N; R(n)=0; for k=1:1:N-n R(n)=R(n)+Err(k)*Err(k+n-1); end R(n)=R(n)/N; r(n)=R(n