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

系统辨识大作业最小二乘法及其相关估值方法应用学院：自动化学院专业：信息工程学号：2007302171姓名：马志强日期：2010.11.14基于最小二乘法的多种系统辨识方法研究1. 最小二乘法的引出在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。设单输入-单输出线性定长系统的差分方程为xk+a1xk-1+ank-n=b0uk+bnuk-n，k=1，2，3，(5.1.1)式中: uk为随机干扰；xk为理论上的输出值。xk只有通过观测才能得到，在观测过程中往往附加有随机干扰。xk的观测值yk可表示为yk=xk+nk(5.1.2)式中: nk为随机干扰。由式(5.1.2)得xk=yk-nk(5.1.3)将式(5.1.3)带入式(5.1.1)得yk+a1yk-1+anyk-n=b0uk+b1uk-1+bnuk-n+nk+i=1nai(k-i)(5.1.4)我们可能不知道nk的统计特性，在这种情况下，往往把nk看做均值为0的白噪声。设(k)=n(k)+i=1nai(k-i)(5.1.5)则式(5.1.4)可写成yk=-a1yk-1-a2yk-2-anyk-n+b0uk+b1uk-1+bnuk-n+(k)(5.1.6)在观测uk时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定(k)不仅包含了xk的测量误差，而且包含了uk的测量误差和系统内部噪声。假定(k)是不相关随机序列(实际上(k)是相关随机序列)。现分别测出n+N个随机输入值y1，y2，yn+N，u1，u2，u(n+N)，则可写成N个方程，即yn+1=-a1yn-a2yn-1-any1+b0un+1+b1un+bnu1+(n+1)yn+2=-a1yn+1-a2yn-any2+b0un+2+b1un+1+bnu2+(n+2)yn+N=-a1yn+N-1-a2yn+N-2-anyN+b0un+N+b1un+N-1+bnuN+(n+N)上述N个方程可写成向量-矩阵形式y(n+1)y(n+2)y(n+N)=-y(n)-y(n+1)-y(1)u(n+1)u(1)-y(2)u(n+2)u(2) -y(n+N-1)-y(N)u(n+N)u(N)a1anb0bn+(n+1)(n+2)(n+3)(5.1.7)设y=y(n+1)y(n+2)y(n+N)，=a1anb0bn，=(n+1)(n+2)(n+N)=-y(n)-y(n+1)-y(1)u(n+1)u(1)-y(2)u(n+2)u(2) -y(n+N-1)-y(N)u(n+N)u(N)则式(5.1.7)可写为y=+(5.1.8)式中：y为N维输出向量；为N维噪声向量；为(2n+1)维参数向量；为N(2n+1)测量矩阵。因此式(5.1.8)是一个含有(2n+1)个未知参数，由N个方程组成的联立方程组。如果N2n+1，方程数少于未知数数目，则方程组的解是不定的，不能唯一地确定参数向量。如果N=2n+1，方程组正好与未知数数目相等，当噪声=0时，就能准确地解出=-1y(5.1.9)如果噪声0，则=-1y-1(5.1.10)从上式可以看出噪声对参数估计是有影响的，为了尽量较小噪声对估值的影响。在给定输出向量y和测量矩阵的条件下求系统参数的估值，这就是系统辨识问题。可用最小二乘法来求的估值，以下讨论最小二乘法估计。2. 最小二乘法估计算法设表示 的最优估值，y表示 y的最优估值，则有y=(5.1.11)y=y(n+1)y(n+2)y(n+N)，=a1anb0bn写出式(5.1.11)的某一行，则有yk=-a1yk-1-a2yk-2-anyk-n+b0uk+b1uk-1+bnuk-n+k=-i=1naiy(k-i)+i=0nbiu(k-i)，k=n+1，n+2，n+N(5.1.12)设e(k)表示y(k)与y(k)之差，即ek= yk-yk=yk-yk=yki=1naiyk-i+i=0nbiuk-i=1+a1z-1+anz-nyk-b0+b1z-1+bnz-nuk=az-1yk-bz-1uk，k=n+1，n+2，n+N(5.1.13)式中az-1=1+a1z-1+anz-nbz-1=b0+b1z-1+bnz-nek成为残差。把k=n+1，n+2，n+N分别代入式(5.1.13)可得残差en+1，en+2，en+N。设e=en+1en+2en+NT则有ek= y-y=y-(5.1.14)最小二乘估计要求残差的平方和为最小，即按照指数函数J= eTe=(y-)T(y-)(5.1.15)为最小来确定估值。求J对的偏导数并令其等于0可得J=-2Ty-=0(5.1.16)T=Ty(5.1.17)由式(5.1.17)可得的最小二乘估计=(T)-1T y(5.1.18)3.递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。设已获得的观测数据长度为N，将式(5.1.8)中的y和分别用YN，N，N来代替，即YN=N+N(5.3.1)用N表示 的最小二乘估计，则N=(NTN)-1NT YN(5.3.2)设PN=（NTN）-1(5.3.5)于是N=PNNTYN(5.3.6)如果再获得1组新的观测值un+N+1和y(n+N+1)，则又增加1个方程yN+1=N+1T+N+1(5.3.7)式中yN+1=yn+N+1，N+1=(n+N+1)N+1T=-yn+N -yN+1 un+N+1 u(N+1)将式(5.3.1)和式(5.3.7)合并，并写成分块矩阵形式，可得YNyN+1=NN+1T +NN+1(5.3.8)根据上式可得到新的参数估值N+1=NN+1TTNN+1T-1NN+1TTYNyN+1=PN+1NN+1TTYNyN+1=PN+1(NTYN+N+1yN+1)(5.3.9)式中PN+1=NN+1TTNN+1T-1=（NTN+N+1N+1T）-1根据矩阵求逆引理可以求得递推最小二乘法辨识公式N+1=N+KN+1(yN+1-N+1TN)(5.3.19)KN+1=PN+1N+1(1+N+1TPNN+1)-1(5.3.20)PN+1=PN-PNN+1(1+N+1TPNN+1)-1N+1TPN(5.3.21)由于进行递推计算需要给出N和PN初值P0和0，通过计算证明，可以取初值：0=0，P0=c2I，c是充分大的常数，I为2n+1(2n+1)单位矩阵，则经过若干次递推之后能够得到较好的参数估计。3. 辅助变量法辅助变量法是一种可克服最小二乘有偏估计的一种方法，对于原辨识方程y=+(5.4.1)当(k)是不相关随机序列时，最小二乘法可以得到参数向量的一致无偏估计。但是，在实际应用中(k)往往是相关随机序列。假定存在着一个(2n+1)N的矩阵Z满足约束条件limN1NZT=0limN1NZT=Q(5.4.2)式中Q是非奇异的。用ZT乘以式(5.4.1)等号两边得ZTy=ZT+ZT (5.4.3)由上式得=(ZT)-1ZTy-(ZT)-1ZT(5.4.4)如果取IV=(ZT)-1ZTy(5.4.5)作为估值，则称IV为辅助变量估值，矩阵Z成为辅助变量矩阵，Z中的元素称为辅助变量。常用的辅助变量法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。4. 广义最小二乘法广义最小二乘法是能克服最小二乘法有偏估计的另一种方法，这种方法计算比较复杂但效果比较好。下面直接介绍广义最小二乘法的计算步骤：(1)应用得到的输入和输出数据u(k)和yk(k=1，2，3，n+N)，按模型az-1yk=bz-1uk+k求出 的最小二乘估计(1)=a1(1)an(1)b0(1)bn(1)(2)计算残差e1(k)e1k=a(1)z-1yk-b(1)z-1uk(3)用残差e1k代替 k，计算f(1)=(1)T(1)-1(1)Te1(4)计算y1(k)和u1(k)y1k=yk+f11yk-1+fm1y(k-m)u1k=uk+f11uk-1+fm1u(k-m)(5)应用得到的y1(k)和u1(k)按模型az-1y1k=bz-1u1k+(k)用最小二乘法重新估计，得到的第2次估值(2)。然后按步骤(2)计算残差e2(k)，按步骤(3)重新估计f，得到估值f 2。再按照步骤(4)计算y2(k)和u2k，按照步骤(5)求的第3次估值(3)。重复上述循环，之道的估值(i)收敛为止。5. 一种交替的广义最小二乘法求解技术(夏式法)这种方法是夏天长提出来的，又称夏式法。以上讨论过的广义最小二乘法的特点在于系统的输入和输出信号反复过滤。一下介绍的夏式法是一种交替的广义最小二乘法求解技术，它不需要数据反复过滤，因而计算效率较高。这种方法可消去最小二乘估计中的偏差，而且由这种方法导出的计算方法也比较简单。基于以上的几种方法，有y=+=f+(5.7.1)因而有y=+f+=f+(5.7.2)应用最小二乘法可得到参数估值f=TT-1TTy(5.7.3)可以推出=(T)-1T y-(T)-1f(5.7.11)上式中的第1项是最小二乘估计LS，第2项是偏差项B，所以必须准确计算B。为了准确计算B，可采用迭代的方法。6. 专题解答设但输入-单输出系统的差分方程为yk=-a1yk-1-a2yk-2+b1uk-1+b2uk-2+(k)k=k+a1k-1+a2k-2取真实值T=a1a2b1b2=1.6420.7150.390.35，输入数据如下所示ku(k)ku(k)ku(k)11.14711-0.958210.48520.201120.810221.6333-0.78713-0.044230.0434-1.159140.947241.3265-1.05215-1.474251.70660.86616-0.71926-0.34071.15217-0.086270.89081.57318-1.099281.14490.626191.450291.177100.433201.15130-0.390用的真实值利用查分方程求出yk作为测量值，k为均值为0，方差为0.1，0.5的不相关随机序列。(1) 用最小二乘法估计参数T=a1a2b1b2。(2) 用递推最小二乘法估计T=a1a2b1b2。(3) 用辅助变量法估计参数T=a1a2b1b2。(4) 设k+2k=k，用广义最小二乘法估计参数T=a1a2b1b2。(5) 用夏式法估计参数T=a1a2b1b2。(6) 详细分析和比较所获得的参数辨识结果，并说明上述参数便是方法的优缺点。根据题目要求的解法，利用Matlab编程实现系统辨识的估值利用最小二乘法估计的结果如下：最小二乘法方差a1a2b1b20.00011.62800.70280.39713.4491.63160.70590.39470.34941.63540.71200.39180.34631.63620.70820.39700.35271.63600.71650.39060.35341.62890.70460.39080.34370.0011.15430.67660.40640.33041.55770.63710.38680.32491.60500.68600.37370.32441.60600.68160.35830.31671.61950.70300.39070.33661.56700.65720.37520.31400.011.35380.50100.34860.17090.89560.16370.42370.06971.00080.20360.42360.12681.34030.47070.38260.26151.05740.22890.36820.13171.12310.29630.35920.15060.11.14240.27100.32840.22161.02550.17360.37660.18440.88960.11090.38930.08130.81820.11140.42980.09230.81000.01530.41220.13520.77150.13110.47140.06400.50.97510.10170.02710.07920.89380.07400.34840.16340.19270.01970.37620.15210.55060.03920.65100.05840.75600.04940.33720.14370.94590.13770.38150.1853部分程序运行结果递推最小二乘法方差a1a2b1b20.00011.67540.67870.52070.39241.30760.29000.00750.32841.51460.69631.14010.16391.57330.77820.41490.71101.16020.47530.67360.31591.20910.31920.52770.02750.0011.47670.40400.26790.55121.62590.75940.32530.37801.53930.47570.12680.43461.15480.17000.19260.82510.88580.07600.33850.04061.41290.31270.09920.83800.011.34850.34450.31940.37101.16390.32960.78130.21621.99461.23231.48520.03041.39240.35430.33190.45721.39820.36080.77730.31521.63460.72290.57800.39470.11.56240.71320.44220.41121.73350.71520.08440.63991.47630.53660.32550.31161.44770.34890.22180.22651.62160.70820.65950.42751.51050.40000.01130.22130.51.79270.94110.27300.34711.55560.88770.59720.12171.78681.25381.12480.21001.57330.74340.35890.13871.31930.60841.29710.30291.59590.53860.01410.6947部分程序运行结果：辅助变量法方差a1a2b1b20.00011.77990.85880.41470.40461.30760.29000.00750.32841.67350.75780.40600.34841.58120.65460.37710.36111.66570.74690.37720.35611.52810.65090.36450.33020.0011.62950.67750.40820.35011.64250.73050.39370.35431.55950.60520.35630.33621.41450.49250.40210.28731.63710.72700.36810.34181.35390.47330.39060.24890.011.34510.47000.38220.26801.36570.48930.43490.24961.37020.50090.43880.26111.18840.37070.34570.15211.36360.53300.41350.20641.31580.48360.45620.21690.11.55450.61670.41040.35121.59000.66480.40980.36641.66100.70290.40010.34951.51040.60140.39680.32561.56200.64960.39500.30971.44180.57060.41830.27830.51.49520.57040.37690.34831.55920.65410.43870.33301.36370.53020.38650.19161.55430.63240.32080.29291.53850.59940.36610.35761.35110.48390.34120.2208部分程序运行结果：广义最小二乘法方差a1a2b1b20.00011.64510.71820.38950.35031.64680.72050.38660.34921.63080.70380.38520.34571.64670.71880.38980.35441.64110.71560.38980.34771.64410.71850.38960.34880.0011.70270.76780.40590.37611.64770.71660.39730.36431.63910.71130.38360.35041.65860.71660.37910.36721.68310.74720.39940.37311.59640.68500.36550.32310.011.68260.74690.38590.36001.72450.76960.35650.36421.65770.71860.39100.38201.66560.72630.34250.34131.69030.73920.36330.38381.69420.74390.38340.38800.10.89960.15150.44350.06041.03800.30440.34100.04641.65710.65690.45890.42591.21800.45840.31520.10300.19570.09740.21640.10351.68590.72320.41440.44040.51.40530.67030.36720.13721.57830.61480.35120.41731.02330.06580.20350.40601.59730.81030.35440.21311.61910.67090.31110.48661.08040.22100.75000.5185部分程序运行结果：夏式法方差a1a2b1b20.00011.63670.70900.39100.35111.64510.71840.38760.34721.36180.21590.82050.58511.62910.70750.39110.34641.61730.69290.39550.34391.63070.70660.39120.34490.0011.53140.64070.40160.30921.56320.65750.42040.33191.61720.69460.39250.32581.56700.64890.40160.32331.53260.62430.39950.30701.62510.69460.36210.33600.010.88150.10890.37340.04401.37630.49720.41210.26541.50270.59520.38540.32771.34210.42680.39780.27901.24540.30240.36690.25251.33610.47410.49990.29000.10.66390.10740.68470.08271.14060.31540.42590.23670.17910.29160.94540.21261.20680.38270.74600.32200.64120.07040.58020.08771.03840.28780.74930.37930.50.93300.07520.97410.58820.62910.16860.46790.78690.59860.27920.26050.20080.44410.80330.74421.76670.62160.14400.58210.49060.65760.09991.00020.0263部分程序运行结果：结论：通过编程计算，获得在噪声方差比较小的情况下，各种方法所获得的估值比较理想，但随着噪声方差的增大，估值的偏差随之增大，横向比较看来夏式法与广义最小二乘法能够更好地还原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值。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.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) -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(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(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; %取相对误差E=0.000000005 c=c0,zeros(4,30); %被辨识参数矩阵的初始值及大小 e=zeros(4,30); %相对误差的初始值及大小for k=3:30; %开始求K h1=-z(k-1),-z(k-2),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 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; 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 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=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 = 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)+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,:)=-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); 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 -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