系统辨识最小二乘法大作业-(2)

系统辨识大作业最小二乘法及其相关估值方法应用学院：自动化学院学号：姓名：日期：基于最小二乘法的多种系统辨识方法研究一、实验原理1.最小二乘法在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。设单输入-单输出线性定长系统的差分方程为(5.1.1)式中:为随机干扰；为理论上的输出值。只有通过观测才能得到，在观测过程中往往附加有随机干扰。的观测值可表示为(5.1.2)式中:为随机干扰。由式(5.1.2)得(5.1.3)将式(5.1.3)带入式(5.1.1)得(5.1.4)我们可能不知道的统计特性，在这种情况下，往往把看做均值为0的白噪声。设(5.1.5)则式(5.1.4)可写成(5.1.6)在观测时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定不仅包含了的测量误差，而且包含了的测量误差和系统内部噪声。假定是不相关随机序列(实际上是相关随机序列)。现分别测出个随机输入值，则可写成个方程，即上述个方程可写成向量-矩阵形式(5.1.7)设则式(5.1.7)可写为(5.1.8)式中：为维输出向量；为维噪声向量；为维参数向量；为测量矩阵。因此式(5.1.8)是一个含有个未知参数，由个方程组成的联立方程组。如果，方程数少于未知数数目，则方程组的解是不定的，不能唯一地确定参数向量。如果，方程组正好与未知数数目相等，当噪声时，就能准确地解出(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)为最小来确定估值。求对的偏导数并令其等于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.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)作为估值，则称为辅助变量估值，矩阵成为辅助变量矩阵，中的元素称为辅助变量。常用的辅助变量法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。5.广义最小二乘法广义最小二乘法是能克服最小二乘法有偏估计的另一种方法，这种方法计算比较复杂但效果比较好。下面直接介绍广义最小二乘法的计算步骤：(1)应用得到的输入和输出数据和，按模型求出的最小二乘估计(2)计算残差(3)用残差代替，计算(4)计算和(5)应用得到的和按模型用最小二乘法重新估计，得到的第2次估值。然后按步骤(2)计算残差，按步骤(3)重新估计，得到估值。再按照步骤(4)计算和，按照步骤(5)求的第3次估值。重复上述循环，之道的估值收敛为止。二、实验要求设但输入-单输出系统的差分方程为取真实值，输入数据：X1=1;X2=1;X3=1;X4=1;%移位寄存器输入Xi初T态（1111），Yi为移位寄存器各级输出n=2000;%置M序列总长度fori=1:nY4=X4;Y3=X3;Y2=X2;Y1=X1;X4=Y3;X3=Y2;X2=Y1;X1=xor(Y3,Y4);%异或运算U(i)=Y4;ifY4==0U(i)=1;elseU(i)=Y4*(-1);endendu=U;用的真实值利用查分方程求出作为测量值，为均值为0，方差为0.1，0.5的不相关随机序列。(1)用最小二乘法估计参数。(2)用递推最小二乘法估计。(3)用辅助变量法估计参数。(4)设，用广义最小二乘法估计参数。(5)详细分析和比较所获得的参数辨识结果，并说明上述参数便是方法的优缺点。根据题目要求的解法，利用Matlab编程实现系统辨识的估值三、实验结果实验方法方差试验次数误差最小二乘0.510.17080.20700.39200.20891.814820.16800.19830.40080.21191.774830.18200.17240.38430.23471.820740.22300.19690.38570.22341.82390.110.82810.13730.39080.03351.047120.80050.11910.38660.02071.082430.79630.11470.38990.01791.089040.78280.10300.38870.01601.10650.0111.62170.70080.39010.34180.026121.62240.70070.38930.34240.025431.62140.70080.38940.34100.026741.62250.70070.38980.34260.02530.00111.64190.71490.39000.34990.000221.64170.71480.38990.34990.000331.64170.71490.38990.34980.000441.64180.71480.39000.35000.0002实验方法方差试验次数误差递推最小二乘0.510.20150.20590.39510.21111.866020.16340.18940.39970.24101.866030.22590.17470.38670.20271.866040.17970.19800.41510.21481.86600.110.82720.13640.39060.03001.686020.76530.08920.39400.01461.866030.82360.13110.39380.03401.866040.78290.10650.38920.01131.68000.0111.62230.70070.39120.34300.023021.62360.70140.39080.34390.024531.62050.69990.39060.34190.023441.62240.70060.39050.34290.02360.00111.64180.71480.39000.35000.00121.64190.71490.39000.35000.00231.64180.71490.39000.35000.002141.64170.71480.38990.34980.0031实验方法方差试验次数误差辅助变量0.510.72380.09100.37110.00151.865220.00150.18310.40550.14761.866031.20770.42350.41090.17931.866040.71850.11550.38400.03811.85930.111.63520.70620.39610.35571.866021.59110.68530.39550.32841.866031.61270.69590.39620.34111.866041.63180.70660.39050.34671.86600.0111.64210.71500.38940.35000.235421.64200.71480.39030.35050.261231.64010.71440.39000.34870.243141.64230.71490.38920.35000.23450.00111.64210.71500.39000.35010.000121.64200.71490.39000.35000.000231.64210.71500.39000.35000.000141.64210.71500.39000.35010.0001实验方法方差试验次数误差广义最小二乘0.510.90190.23380.39850.04480.975020.97130.27980.39530.07250.902330.89120.24370.40140.02290.781440.99830.30690.40350.06930.93410.111.61580.69940.39150.33910.032421.61400.69800.39310.34030.034331.61410.69770.39050.33830.034941.61940.70140.39270.34080.02810.0111.64370.71620.39010.35090.002321.64440.71660.38980.35080.003031.64390.71640.39000.35050.002441.64410.71620.39000.35100.00260.00111.64200.71500.39000.35000.000021.64190.71500.39000.34990.000131.64200.71500.39000.35000.000141.64190.71500.39000.35000.00011）输入信号图像：2）递推最小二乘图像3）辅助变量法最小二乘四、实验结论1)在噪声方差比较小的情况下，各种方法所获得的估值比较理想，但随着噪声方差的增大，估值的偏差随之增大。当方差为0.5时，估计值与真实值之间偏差很大。当噪声方差小于0.01以后辨识结果非常接近真实值，且各种辨识方法直接差别很小。2)横向比较看来广义最小二乘法能够更好地还原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值(最小二乘是有偏估计，广义最小二乘是无偏估计)。3）从表格的误差以来也可以明显的看出，广义最小二乘效果最好。4）虽然从获得的图形看出递推最小二乘在很少的数据量的时候就可以达到收敛，但是实际运行时数据量太少辨识结果往往非常差。本实验就生成2000个M序列来代替原始的数据输入。5）由于噪声序列的随机性，往往导致实验结果存在某一时刻结果非常好，对实验噪声随机序列进行固定可以获得很好的结果。6）总的来说随着输入数据的增加便是结果会很接近真实值，而且广义最小二乘辨识结果一直很好。五、程序代码%A=179;N=2000;%初始化;%x0=11;M=2^35;%fork=1:N%乘同余法递推2000次;%x2=A*x0;%x2和x0分别表示xi和xi-1;%x1=mod(x2,M);%将x2存储器的数除以M，取余数放x1中（xi）;%v1=x1/2^35;%将x1存储器的数除以M得到小于1的随机数放v1中;%v(:,k)=v1;%将v1中的数存放在矩阵存储器的第k列中;%x0=x1;%xi-1=xi;%v0=v1;%end%递推2000次结束;%u=v;%周期为15的M序列X1=1;X2=1;X3=1;X4=1;%移位寄存器输入Xi初T态（1111），Yi为移位寄存器各级输出n=2000;%置M序列总长度fori=1:nY4=X4;Y3=X3;Y2=X2;Y1=X1;X4