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

系统辨识大作业最小二乘法及其相关估值方法应用学号：2012302259姓名：王家琦基于最小二乘法的多种系统辨识方法研究最小二乘法的引出在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。设单输入-单输出线性定长系统的差分方程为x (5.1.1)式中:uk为随机干扰；xk为理论上的输出值。xk只有通过观测才能得到，在观测过程中往往附加有随机干扰。xy (5.1.2)式中:nk为随机干扰。由式(5.1.2)x (5.1.3)将式(5.1.3)带入式(5.1.1)得y (5.1.4)我们可能不知道nk的统计特性，在这种情况下，往往把nk看做均值为0 设ξ( (5.1.5)则式(5.1.4)可写成y (5.1.6) 在观测uk时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定ξ(k)不仅包含了xk的测量误差，而且包含了uk的测量误差和系统内部噪声。假定ξ(k) 现分别测出n+N个随机输入值y1，yyy⋮y上述N个方程可写成向量-矩阵形式设P (5.3.5)于是θ (5.3.6) 如果再获得1组新的观测值un+N+1和y(n+N+1)，则又增加1y (5.3.7)式中yψ将式(5.3.1)和式(5.3.7)合并，并写成分块矩阵形式，可得Y (5.3.8)根据上式可得到新的参数估值θ (5.3.9)式中P根据矩阵求逆引理可以求得递推最小二乘法辨识公式θ (5.3.19)K (5.3.20)P (5.3.21)由于进行递推计算需要给出θN和PN初值P0和θ0，通过计算证明，可以取初值：θ0=0，P0辅助变量法辅助变量法是一种可克服最小二乘有偏估计的一种方法，对于原辨识方程y=Φθ+ξ (5.4.1)当ξ(k)是不相关随机序列时，最小二乘法可以得到参数向量θ的一致无偏估计。但是，在实际应用中ξ(k)往往是相关随机序列。假定存在着一个(2n+1)×N的矩阵Zlim (5.4.2)式中Q是非奇异的。用ZT乘以式(5.4.1)Z (5.4.3)由上式得θ= (5.4.4)如果取θ (5.4.5)作为θ估值，则称θIV为辅助变量估值，矩阵Z成为辅助变量矩阵，Z中的元素称常用的辅助变量法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。广义最小二乘法广义最小二乘法是能克服最小二乘法有偏估计的另一种方法，这种方法计算比较复杂但效果比较好。 下面直接介绍广义最小二乘法的计算步骤：(1)应用得到的输入和输出数据u(k)和yk(k=1a求出θ的最小二乘估计θ(2)计算残差ee (3)用残差e1k代替ξ (4)计算y1(k)yu (5)应用得到的y1(k)和a用最小二乘法重新估计θ，得到θ的第2次估值θ(2)。然后按步骤(2)计算残差e2(k)，按步骤(3)重新估计f，得到估值f2。再按照步骤(4)计算y2(k)和u2k，按照步骤(5)求一种交替的广义最小二乘法求解技术(夏式法)这种方法是夏天长提出来的，又称夏式法。以上讨论过的广义最小二乘法的特点在于系统的输入和输出信号反复过滤。一下介绍的夏式法是一种交替的广义最小二乘法求解技术，它不需要数据反复过滤，因而计算效率较高。这种方法可消去最小二乘估计中的偏差，而且由这种方法导出的计算方法也比较简单。基于以上的几种方法，有y=Φθ+ξ (5.7.1)因而有y= (5.7.2)应用最小二乘法可得到参数估值θ (5.7.3)可以推出θ (5.7.11)上式中的第1项是θ最小二乘估计θLS，第2项是偏差项θB，所以必须准确计算 为了准确计算θB专题解答设但输入-单输出系统的差分方程为yξ取真实值θTku(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，用最小二乘法估计参数θT用递推最小二乘法估计θT用辅助变量法估计参数θT设ξk+2用夏式法估计参数θT详细分析和比较所获得的参数辨识结果，并说明上述参数便是方法的优缺点。根据题目要求的解法，利用Matlab编程实现系统辨识的估值利用最小二乘法估计的结果如下：最小二乘法方差aabb0.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部分程序运行结果递推最小二乘法方差aabb0.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部分程序运行结果：辅助变量法方差aabb0.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部分程序运行结果：广义最小二乘法方差aabb0.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部分程序运行结果：夏式法方差aabb0.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源程序：%最小二乘估计%clearu=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];n=normrnd(0,sqrt(0.1),1,31);z=zeros(1,30);fork=3:31z(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);endh0=[-z(2)-z(1)u(2)u(1)]';HLT=[h0,zeros(4,28)];fork=3:30h1=[-z(k)-z(k-1)u(k)u(k-1)]';HLT(:,k-1)=h1;endHL=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)];%求出FAIc1=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.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];z(2)=0;z(1)=0;n=normrnd(0,sqrt(0.1),1,31);fork=3:31z(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);endc0=[0.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%取相对误差E=0.000000005c=[c0,zeros(4,30)];%被辨识参数矩阵的初始值及大小e=zeros(4,30);%相对误差的初始值及大小fork=3:30;%开始求Kh1=[-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;%求被辨识参数ce1=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;%给下次用ife2<=Ebreak;%如果参数收敛情况满足要求，终止计算endend%display('方差为0.0001递推最小二乘法辨识后的结果是：');a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%display('a1,a2,b1,b2经过递推最小二乘法辨识的结果是：');fori=3:31;if(c(1,i)==0)q1=c(1,i-1);break;endendfori=3:31;if(c(2,i)==0)q2=c(2,i-1);break;endendfori=3:31;if(c(3,i)==0)q3=c(3,i-1);break;endendfori=3:31;if(c(4,i)==0)q4=c(4,i-1);break;endenda1=q1;a2=q2;b1=q3;b2=q4;%clear%辅助变量递推最小二乘法估计na=2;nb=2;siitt=[1.6420.7150.390.35]';siit0=0.001*eye(na+nb,1);p=10^6*eye(na+nb);siit(:,1)=siit0;y(2)=0;y(1)=0;x(1)=0;x(2)=0;j=0;u=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];n=normrnd(0,sqrt(0.01),1,31);fork=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.6420.7150.390.35]')^2;e=max(abs((siit(:,k-1)-siit0)./siit0));ess(:,k-2)=siit(:,k-1)-siitt;siit0=siit(:,k-1);enda1=siit0(1);a2=siit0(2);b1=siit0(3);b2=siit0(4);clear%广义最小二乘估计clear;nn=normrnd(0,sqrt(0.5),1,31)';uk=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9580.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8901.1441.177-0.390];yk(1)=0;yk(2)=0;fori=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;fori=1:29;A(i,:)=[-yk(i+1)-yk(i)uk(i+1)uk(i)];endsiit=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);fori=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);endfori=1:29;fai(i,:)=[-e(i+1)-e(i)];endf=inv(fai'*fai)*fai'*e(3:31)';fori=3:31;yk(i)=yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);endyk(2)=yk(2)+f(1)*yk(1);fori=3:30;uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);forj=1:30fori=1:29;A(i,:)=[-yk(i+1)-yk(i)uk(i+1)uk(i)];endsiit=inv(A'*A)*A'*yk(3:31)';e(1)=yk(1);e(2)=yk(2)+siit(1)*(yk(1))-siit(3)*uk(1);fori=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);endfori=1:29;fai(i,:)=[-e(i+1)-e(i)];endf=inv(fai'*fai)*fai'*e(3:31)';k1(j)=f(1);k2(j)=f(2);fori=3:31;yk(i)=yk(i)+f(1)*(yk(i-1))+f(2)*(yk(i-2));endyk(2)=yk(2)+f(1)*yk(1);fori=3:30uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);endsiit';%用夏氏偏差修正法估计参数clear;u=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];n=normrnd(0,sqrt(0.5),1,31);z=zeros(1,30);fork=3:31z(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);endh0=[-z(2)-z(1)u(2)u(1)]';HLT=[h0,zeros(4,28)];fork=3:30h1=[-z(k