系统辨识最小二乘法大作业.docx

1、系统辨识最小二乘法大作业系统辨识大作业最小二乘法及其相关估值方法应用学院：自动化学院专业：信息工程学号：2021302171姓名：马志强日期：基于最小二乘法的多种系统辨识方法研究1.最小二乘法的引出在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。设单输入-单输出线性定长系统的差分方程为 (5.1.1)式中:为随机干扰；为理论上的输出值。只有通过观测才能得到，在观测过程中往往附加有随机干扰。的观测值可表示为 (5.1.2)式中:为随机干扰。由式(5.1.2)得 (5.1.3)将式(5.1.3)带入式(5.1.1)得 (5.1.4)我们可能不知道的统计特性，在这种情况下，往往把看做均值为0的

2、白噪声。 设 (5.1.5)那么式(5.1.4)可写成 (5.1.6) 在观测时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定不仅包含了的测量误差，而且包含了的测量误差和系统内部噪声。假定是不相关随机序列(实际上是相关随机序列)。 现分别测出个随机输入值，那么可写成个方程，即上述个方程可写成向量-矩阵形式 (5.1.7)设那么式(5.1.7)可写为 (5.1.8)式中：为维输出向量；为维噪声向量；为维参数向量；为测量矩阵。因此式(5.1.8)是一个含有个未知参数，由个方程组成的联立方程组。如果，方程数少于未知数数目，那么方程组的解是不定的，不能唯一地确定参数向量。如果，方程组

3、正好与未知数数目相等，当噪声时，就能准确地解出 (5.1.9)如果噪声，那么 (5.1.10)从上式可以看出噪声对参数估计是有影响的，为了尽量较小噪声对估值的影响。在给定输出向量和测量矩阵的条件下求系统参数的估值，这就是系统辨识问题。可用最小二乘法来求的估值，以下讨论最小二乘法估计。2.最小二乘法估计算法设表示的最优估值，表示的最优估值，那么有 (5.1.11)写出式(5.1.11)的某一行，那么有 (5.1.12)设表示与之差，即- (5.1.13)式中成为残差。把分别代入式(5.1.13)可得残差。设那么有 (5.1.14)最小二乘估计要求残差的平方和为最小，即按照指数函数 (5.1.15

4、)为最小来确定估值。求对的偏导数并令其等于0可得 (5.1.16) (5.1.17)由式(5.1.17)可得的最小二乘估计 (5.1.18)3.递推最小二乘法 为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。 设已获得的观测数据长度为，将式(5.1.8)中的和分别用来代替，即 (5.3.1)用的最小二乘估计，那么 (5.3.2)设 (5.3.5)于是 (5.3.6) 如果再获得1组新的观测值和，那么又增加1个方程 (5.3.7)式中将式(5.3.1)和式(5.3.7)合并，并写成分块矩阵形式，可得 (5.3.8)根据上式可得到新的参数估值 (5.3.9)式中根据矩阵求逆引理可

5、以求得递推最小二乘法辨识公式 (5.3.19) (5.3.20) (5.3.21)由于进行递推计算需要给出和初值和，通过计算证明，可以取初值：，c是充分大的常数，为单位矩阵，那么经过假设干次递推之后能够得到较好的参数估计。3.辅助变量法辅助变量法是一种可克服最小二乘有偏估计的一种方法，对于原辨识方程 (5.4.1)当是不相关随机序列时，最小二乘法可以得到参数向量的一致无偏估计。但是，在实际应用中往往是相关随机序列。假定存在着一个的矩阵满足约束条件 (5.4.2)式中是非奇异的。用乘以式(5.4.1)等号两边得 (5.4.3)由上式得 (5.4.4)如果取 (5.4.5)作为估值，那么称为辅助变

6、量估值，矩阵成为辅助变量矩阵，中的元素称为辅助变量。常用的辅助变量法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。4.广义最小二乘法广义最小二乘法是能克服最小二乘法有偏估计的另一种方法，这种方法计算比拟复杂但效果比拟好。 下面直接介绍广义最小二乘法的计算步骤：(1)应用得到的输入和输出数据和，按模型求出的最小二乘估计(2)计算残差 (3)用残差代替，计算 (4)计算和 (5)应用得到的和按模型用最小二乘法重新估计，得到的第2次估值。然后按步骤(2)计算残差，按步骤(3)重新估计，得到估值。再按照步骤(4)计算和，按照步骤(5)求的第3次估值。重复上述循环，之道的估值收敛为止。5.一种交替的

7、广义最小二乘法求解技术(夏式法)这种方法是夏天长提出来的，又称夏式法。以上讨论过的广义最小二乘法的特点在于系统的输入和输出信号反复过滤。一下介绍的夏式法是一种交替的广义最小二乘法求解技术，它不需要数据反复过滤，因而计算效率较高。这种方法可消去最小二乘估计中的偏差，而且由这种方法导出的计算方法也比拟简单。基于以上的几种方法，有 (5.7.1)因而有 (5.7.2)应用最小二乘法可得到参数估值 (5.7.3)可以推出 (5.7.11)上式中的第1项是最小二乘估计，第2项是偏差项，所以必须准确计算。 为了准确计算，可采用迭代的方法。6.专题解答设但输入-单输出系统的差分方程为取真实值，输入数据如下所

8、示ku(k)ku(k)ku(k)111212122231323414245152561626717278182891929102030用的真实值利用查分方程求出作为测量值，为均值为0，方差为0.1，0.5的不相关随机序列。(1)用最小二乘法估计参数。(2)用递推最小二乘法估计。(3)用辅助变量法估计参数。(4)设，用广义最小二乘法估计参数。(5)用夏式法估计参数。(6)详细分析和比拟所获得的参数辨识结果，并说明上述参数便是方法的优缺点。根据题目要求的解法，利用Matlab编程实现系统辨识的估值利用最小二乘法估计的结果如下：最小二乘法方差1.6316局部程序运行结果递推最小二乘法方差局部程序运行

9、结果：辅助变量法方差局部程序运行结果：广义最小二乘法方差局部程序运行结果：夏式法方差局部程序运行结果：结论：通过编程计算，获得在噪声方差比拟小的情况下，各种方法所获得的估值比拟理想，但随着噪声方差的增大，估值的偏差随之增大，横向比拟看来夏式法与广义最小二乘法能够更好地复原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值。Matlab源程序：%最小二乘估计% clear u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.

10、086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177 -0.390; n=normrnd(0, sqrt(0.1), 1, 31); z=zeros(1,30);for k=3:31 z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);end h0=-z(2) -z(1) u(2) u(1); HLT=h0,zeros(4,28);for k=3:30 h1=-z(k) -

11、z(k-1) u(k) u(k-1); HLT(:,k-1)=h1;end HL=HLT; y=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16);z(17);z(18);z(19);z(20);z(21);z(22);z(23);z(24);z(25);z(26);z(27);z(28);z(29);z(30);z(31);%求出FAI c1=HL*HL; c2=inv(c1); c3=HL*y; c=c2*c3;%display(方差=0.1时，最小二乘法估计辨识参数如下：); a1=c(

12、1); a2=c(2); b1=c(3); b2=c(4); clear%递推最小二乘法估计 u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177 -0.390; z(2)=0; z(1)=0;n=normrnd(0, sqrt(0.1), 1, 31);for k=3:31 z

13、(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);end c0=0.001 0.001 0.001 0.001; %直接给出被辨识参数的初始值,即一个充分小的实向量 p0=106*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵 E=0.000000005; c=c0,zeros(4,30); %被辨识参数矩阵的初始值及大小 e=zeros(4,30); %相对误差的初始值及大小for k=3:30; %开始求K h1=-z(k-1),-z(k-2),

14、u(k-1),u(k-2); x=h1*p0*h1+1; x1=inv(x); %开始求K(k) 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; %把辨识参数c 列向量参加辨识参数矩阵的最后一列 p1=p0-k1*k1*h1*p0*h1+1; %求出 p(k)的值 p0=p1; %给下次用 if

15、 e2=E break; %如果参数收敛情况满足要求，终止计算 endend%display(方差为0.0001递推最小二乘法辨识后的结果是：); a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);%display(a1,a2,b1,b2经过递推最小二乘法辨识的结果是：);for i=3:31; if(c(1,i)=0) q1=c(1,i-1); break; endendfor i=3:31; if(c(2,i)=0) q2=c(2,i-1); break; endendfor i=3:31; if(c(3,i)=0) q3=c(3,i-1); break

16、; endendfor i=3:31; if(c(4,i)=0) q4=c(4,i-1); break; endend a1=q1; a2=q2; b1=q3 ; b2=q4;%clear%辅助变量递推最小二乘法估计 na=2; nb=2; siitt=1.642 0.715 0.39 0.35; siit0=0.001*eye(na+nb,1); p=106*eye(na+nb); siit(:,1)=siit0; y(2)=0;y(1)=0; x(1)=0;x(2)=0; j=0; u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573

17、0.626 0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177 -0.390; n=normrnd(0, sqrt(0.01), 1, 31);for k=3:31; h=-y(k-1),-y(k-2),u(k-1),u(k-2); y(k)=h*siitt+n(k)+1.642*n(k-1)+0.715*n(k-2); hx=-x(k-1),-x(k-2),u(k-1),u(k-2); kk=

18、p*hx/(h*p*hx+1); p=eye(na+nb)-kk*h*p; siit(:,k-1)=siit0+kk*y(k)-h*siit0; x(k)=hx*siit(:,k-1); j=j+(y(k)-h*1.642 0.715 0.39 0.35)2; e=max(abs(siit(:,k-1)-siit0)./siit0); ess(:,k-2)=siit(:,k-1)-siitt; siit0=siit(:,k-1); end a1=siit0(1); a2=siit0(2); b1=siit0(3); b2=siit0(4); clear%广义最小二乘估计 clear; nn =

19、 normrnd(0,sqrt(0.5),1,31); uk=1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 1.144 1.177 -0.390; yk(1)=0; yk(2)=0;for i=1:29; yk(i+2)=-1.642*yk(i+1)-0.715*yk(i)+0.39*uk(i+1)+

20、0.35*uk(i)+nn(i+2)+1.642*nn(i+1)+0.715*nn(i);end;for i=1:29; A(i,:)=-yk(i+1) -yk(i) uk(i+1) uk(i);end siit=inv(A*A)*A*(yk(3:31)+nn(2:30); e(1)=yk(1); e(2)=yk(2)+siit(1)*yk(1)-siit(3)*uk(1);for i=3:31; e(i)=yk(i)+siit(1)*yk(i-1)+siit(2)*yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2);endfor i=1:29; fai(i,:)=

21、-e(i+1) -e(i);end f=inv(fai*fai)*fai*e(3:31);for i=3:31; yk(i)=yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);end yk(2)=yk(2)+f(1)*yk(1);for i=3:30; uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);end uk(2)=uk(2)+f(1)*uk(1);for j=1:30for i=1:29; A(i,:)=-yk(i+1) -yk(i) uk(i+1) uk(i);end siit=inv(A*A)*A*yk(3:31); e(1)=yk(1);

22、e(2)=yk(2)+siit(1)*(yk(1)-siit(3)*uk(1);for i=3:31; e(i)=yk(i)+siit(1)*(yk(i-1)+siit(2)*(yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2);endfor i=1:29; fai(i,:)=-e(i+1) -e(i);end f=inv(fai*fai)*fai*e(3:31); k1(j)=f(1); k2(j)=f(2);for i=3:31; yk(i)=yk(i)+f(1)*(yk(i-1)+f(2)*(yk(i-2); endyk(2)=yk(2)+f(1)*yk(1);for i=3:30 uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);end uk(2)=uk(2)+f(1)*uk(1);end siit;%用夏氏偏差修正法估计参数 clear; u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 0.810