广义递推最小二乘辨识.docx

系统辨识上机实验报告广义递推最小二乘辨识一、实验目的1通过实验掌握广义最小二乘辨识算法；2运用MATLAB编程，掌握算法实现方法。二、实验原理广义最小二乘法的根本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。如果滤波模型选择得适宜，对数据进行了较好的白色化处理，那么直接利用普通最小二乘法就能获得无偏一致估计。广义最小二乘法所用的滤波模型实际上就是一种动态模型，在整个迭代过程中不断靠偏差信息来调整这个滤波模型，使它逐渐逼近于一个较好的滤波模型，以便对数据进行较好的白色化处理，使模型参数估计称为无偏一致估计。理论上说，广义最小二乘法所用的动态模型经过几次迭代调整后，便可对数据进行较好的白化处理，但是，当过程的输出噪信比比拟大或模型参数比拟多时,这种数据白色化处理的可靠性就会下降。此时，准那么函数可能出现多个局部收敛点，因而辨识结果可能使准那么函数收敛于局部极小点上而不是全局极小点上。这样，最终的辨识结果往往也会是有偏的。其收敛速度比拟慢，需要经过屡次迭代计算，才能得到较准确的参数估计值。一般情况下，经过屡次迭代后，估计值便会收敛到稳态值。但在某些情况下（如噪声比拟低时）存在局部极小值，估计值不一定收敛到准那么函数的全局极小值上。为了防止参数估计值收敛到局部极小值，最好选定初值接近最优解，一般可以用最小二乘法的批处理估计值作为初值。如果系统是时变的，或为了克服数据饱和现象，可以在两次RLS算法中分别引进遗忘因子。三、实验内容<1>数据获取实验数据按照表10-1,为二阶线性离散系统的输入输出数据<2>数据处理为了提高辨识精度，实验者必须对原始数据进行剔除坏数据、零均值化、工频滤波等处理。实验进行了白化滤波处理。<3>辨识算法利用处理过的数据（取适当的数据长度），选择某种辨识方法（如RLS递推最小二乘法、RELS>RlV或RML等参数估计算法及F-检验或AlC定阶法），估计出模型参数和阶次，同时分析辨识结果。本实验采用广义递推最小二乘法进行系统辨识。三、广义递推最小二乘法（RLS）原理广义最小二乘法是用迭代的松弛算法对最小二乘估计的一种改良，它的根本思想是引入一个白化滤波器，把相关噪声转换为白噪声，基于对观测数据先进行一次滤波处理然后利用普通最小二乘法对滤波后的数据进行辨识。广义最小二乘法的计算步骤如下：1给定初始条件：包括给定的输入输出数据或者产生的数据序列，初始状态矩阵PO,被辨识参数的初始值（取一个充分小的实向量），滤波器参数与矩阵初值。2利用式Z（Q=C（z->z（Q=C（z"）（Q计算滤波后的输入输出序列。3对于二阶离散系统，利用式号优）=-2/（4-1）,-2/（2-2）,勺伏一1）,勺（2-2）7'构造（女）。4利用三个式子递推计算辨识矩阵*）5利用式ee）=z（幻-力（幻仇幻计算贸&）,并根据小伏）=-«也-1）,-6（&-2）7构造4（幻。6利用三个式子递推计算。（幻。7返回第二步进行迭代计算，直至获得满意的辨识结果。四、实验步骤<1>输入输出数据：u=1.147,0.201,-0.787,-1.584-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,0.433,-1.177,-0.390.- 0.982,1.435,-0.119,-0.769,-0.899,0.882,-1.008,-0.844,0.628,-0.679.1.541,1.375,-0.984,-0.582,1.609,0.090,-0.813,-0.428,-0.848,-0.410.0.048,-1.099,-1.108,0.259,-1.627,-0.528,0.203,1.204,1.691,-1.235.- 1.228,-1.267,0.309,0.043,0.043,1.461,1.585,0.552,-0.601,-0.319.0.7440.829,-1.626,-0.127,-1.578,-0.822,1.469,-0.379,-0.212,0.178.0.493-0.056,-0.1294,1.228,-1.606,-0.382,-0.229,0.313,-0.161,-0.810.- 0.2770.983,-0.288,0.846,1.325,0.723,0.713,0.6430.463,0.786.1.161,0.850,-1.349,-0.596,1.512,0.795,-0.713,0.453,-1.604,0.889.-0.938,0.056,0.829,-0.981,-1.232,1.327,-0.681,0.114,-1.135,1.284.-1.2010.758,0.590,-1.007,0.390,0.836,-1.52,-1.053,-0.083,0.619.0.840-1.258,-0.354,0.629,-0.242,1.680,-1.236,-0.803,0.537,-1.100.1.417,-1.024,0.671,0.688,-0.123,-0.952,0.232,-0.793,-1.138,1.154.0.206,1.196,1.013,1.518,-0.553,-0.987,0.167,-1.445,0.630,1.255.0.311,-1.726,0.975,1.718,1.360,1.667,1.111,1.018,0.078,-1.665.-0.760,1.184,-0.614,0.994,-0.089,0,947,1.706,-0.395,1.222,-1.351.0.231,1,425,0.114,-0.689,-0.704,1.070,0.262,1.610,1.489,-1.602.0.020,-0.601,-0.020,-0.601,-0.235,1.245,1.226,-0.204,0.926,-1.297;figure(l);stem(u)gridontitle('图1输入信号')y=1.3813.7942.481-0.280-2.742-1.5542.1292.6913.4272.199.1.679-1.2491.3710.6373.131-0.8190.2351.2622.8493.374.2.3460.6643.0150.5612.2713.6500.6252.3050.3641.857.- 0.912-1.5471.9400.262-0.379-0.1763.7200.058-0.7521.983.- 0.9233.3614.240-0.074-0.4813.7802.1370.0860.638-0.971.- 0.9290.679-0.664-0.4331.570-2.785-1.1530.8193.4844.091.- 2.375-2.561-2.7782.9111.3620.7353.1183.7702.381-0.812.- 1.6350.5891.550-3.410-1.249-3.692-2.3582.552-0.2280.554.2.1782.4710.743-0.0042.504-3.204-1.800-1.2840.1590.426.0.0590.3952.371-0.1572.2483.2972.3292.7802.3751.873.2.4113.9282.846-2.215-1.1043.4602.8830.245-0.231-2.963.2.072-0.845-0.0741.0372.468-3.6792.149-0.0811.639-1.291.2.548-1.6812.3072.227-1.5580.0082.055-1.102-1.4270.350.2.7362.965-2.346-1.5100.809-0.5922.706-1.9412.2752.802.- 1.3372.091-2.5850.0131.2170.691-0.4912.1140.333-0.482.3.3882.0823.7974.0795.0361.250-1.019-0.160-3.2011.161.3.9261.789-2.7032.1255.0544.6785.236-0.2412.1520.356.- 3.5192.2131.527-1.2062.1510.2641.5952.864-0.5391.982.- 3.104-0.2642.4330.009-1.360-0.5213.3191.4453.1053.783.- 1.973-0.138-0.452-0.586-4.045-1.7432.5773.8490.3671.324;<2>初始值设置，包括被辨识参数的初始值、误差序列以及滤波器参数初值；<3>迭代循环，辨识参数更新，根据误差调整滤波器参数，迭代计算被辨识参数，直至参数符合条件或迭代次数到。<4>计算数据与图形显示，包括辨识参数辨识过程以及误差的收敛情况。五、运行结果显示1输入数据：2被辨识参数：辨识结果：ans=-0.42510.01161.8869-0.54643被辨识参数误差的收敛情况：六、实验源程序%递推的广义最小二乘法进行参数估计clear;closeall;display。广义递推最小二乘辨识);u=1.147,0.201,-0.787,-1.584-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,0.433,-1.177,-0.390.- 0.982,1.435,-0.119,-0.769,-0.899,0.882,1.008,-0.844,0.628,-0.679.1.541,1.375,-0.984,-0.582,1.609,0.090,-0.813,-0.428,-0.848,-0.410.0.048,-1.099,-1.108,0.259,1.627,-0.528,0.203,1.204,1.691,-1.235.- 1.228,-1.267,0.309,0.043,0.043,1.461,1.585,0.552,-0.601,-0.319.0.7440.829,-1.626,-0.127,-1.578,0.822,1.469,0.379,0.212,0.1780.493-0.056,-0.1294,1.228,-1.606,-0.382,-0.229,0.313,-0.161,-0.810.-0.2770.983,-0.288,0.846,1.325,0.723,0.713,0.6430.463,0.786.1.161,0.850,-1.349,-0.596,1.512,0.795,-0.713,0.453,-1.604,0.889.-0.938,0.056,0.829,-0.981,-1.232,1.327,-0.681,0.114,-1.135,1.284.- 1.2010.758,0.590,-1.007,0.390,0.836,-1.52,-1.053,-0.083,0.619.0.840-1.258,-0.354,0.629,-0.242,1.680,-1.236,-0.803,0.537,-1.100.1.417,-1.024,0.671,0.688,-0.1235-0.952,0.232,-0.793,-1.138,1.154.0.206,1.196,1,013,1.518,-0.553,-0.987,0.167,-1.445,0.630,1.255.0.311,-1.726,0.975,1.718,1.360,1.667,1.111,1.018,0.078,-1.665.- 0.760,1.184,-0.614,0.994,-0.089,0.947,1.706,-0.395,1.222,-1.351.0.231,1.425,0.114,-0.689,-0.704,1,070,0.262,1.610,1.489,-1.602.- .020,-0.601,-0.020,-0.601,-0.235,1.245,1.226,-0.204,0.926,-1.297;figure(1);stem(u)gridontitlef图1输入信号')y=1.3813.7942.481-0.280-2.742-1.5542.1292.6913.4272.199.1.679-1.2491.3710.6373.131-0.8190.2351.2622.8493.374.2.3460.6643.0150.5612.2713.6500.6252.3050.3641.857.- 0.912-1.5471.9400.262-0.379-0.1763.7200.058-0.7521.983.- 0.9233.3614.240-0.074-0.4813.7802.1370.0860.638-0.971.- 0.9290.679-0.664-0.4331.570-2.785-1.1530.8193.4844.091.- 2.375-2.561-2.7782.9111.3620.7353.1183.7702.381-0.812.- 1.6350.5891.550-3.410-1.249-3.692-2.3582.552-0.2280.554.2.1782.4710.743-0.0042.504-3.204-1.800-1.2840.1590.426.0.0590.3952.371-0.1572.2483.2972.3292.7802.3751.873.2.4113.9282.846-2.215-1.1043.4602.8830.245-0.231-2.963.2.072-0.845-0.0741.0372.468-3.6792.149-0.0811.639-1.291.2.548-1.6812.3072.227-1.5580.0082.055-1.102-1.4270.350.2.7362.965-2.346-1.5100.809-0.5922.706-1.9412.2752.802.- 1.3372.091-2.5850.0131.2170.691-0.4912.1140.333-0.4823.3882.0823.7974.0795.0361.250-1.019-0.160-3.2011.161.3.9261.789-2.7032.1255.0544.6785.236-0.2412.1520356.- 3.5192.2131.527-1.2062.1510.2641.5952.864-0.5391.982.- 3.104-0.2642.4330.009-1.360-0.5213.3191.4453.1053.783- 1.973-0.138-0.452-0.586-4.045-1.7432.5773.8490.3671.324;%用最小二乘递推算法辨识参数：al,a2,bl,b2c0=0.0010.0010.0010.001,;%被辨识参数的初始值直接给取，取一个充分小的实向量ce0=0.001;0.001;%滤波器C(q-1)的参数初值PO=IOA6*eye(4,4);%初始状态Po也采用直接取方式，直接给出初始状态P0,取个充分大的实数单位矩阵pe0=106*eye(2,2);%计算滤波器的P矩阵初值yf=zeros(1,200);yf(l)=0;yf(2)=0;uf=zeros(1,200);uf(l)=0;uf(2)=0;e=zeros(1,200);%残差序列e(l)=0;e(2)=0;c=c0,zeros(4,199);%被辨识参数矩阵的初始值及大小e=zeros(4,200);%相对误差-残差的初始值及大小ce=zeros(2,200);%滤波器参数每次辨识的结果n=0;%用于统计递推次数%对200组数据进行参数估计fork=3:200;%开始递推运算，求Kyf(k)=y(k)+ce(1)*y(k-1)+ce0(2)*y(k-2);%滤波后的输入序列uf(k)=u(k)+ce(1)*u(k-1)+ce0(2)*u(k-2);%滤波后的输出序列h1=-yf(k-1),-yf(k-2),uf(k-1),uf(k-2),;%求h(k)%滤波后的输出输出组成的辨识数据x=hl'*p*hl+l;xl=inv(x);kl=p*hl*xl%求出K的值dl=yf(k)-hl,*c;CI=CO+kl*dl;%求被辨识参数Cel=clcO;%求参数当前值与上一次的值的差值e2=el./c0;%求参数的相对变化e(:,k)=e2;%把当前相对变化的列向量参加误差矩阵的最后一列CO=C1;%新获得的参数作为下一次递推的旧参数c(:,k)=cl;%把当前所辨识参数的CI列向量参加辨识参数矩阵的最后一列p1=p-kl*k*h*pO*h1+1;%求p(k)值PO=P1;%把当前的p(k)值给下次用%递推算法的实质就是C和P在原来的根底上计算，因此要为下时刻的这两个值做准备e(k)=y(k)-y(k-l),-y(k-2),u(k-l),u(k-2)*cl;%计算残差e(k),注意这里的h的组成he1=-e(k-l);-e(k-2);xe=1+hel'*pe*he1;xel=inv(xe);kel=pe*hel*xel;%求当前的K矩阵del=e(k)-hel,*ce;cel=ce+kel*del;ce(:,k)=cel;%保存滤波器参数eel=cel-ce;ee2=eel.ce;%模型参数的相对误差Ee(:,k)=eel;ce=cel;%为下一时刻做准备pe1=pe-ke1*ke1,*he1,*pe*he1+1;%为下一时刻赋值pe=pel;%为下一时刻P矩阵赋值n=n+l;%完成一次递推，统计值加1end%大循环结束c%显示被辨识参数e%显示辨识结果的收敛情况n%显示收敛条件满足时算法最终递推次数al=c(l);a2=c(2);bl=c(3);b2=c(4,:);eal=e(l,:);ea2=e(2,:);ebl=e(3,:);eb2=e(4,:);%别离参数Al=al(200)A2=a2(200)Bl=bl(2OO)B2=b2(2OO)figure(2);i=1:200;plot(i,a1fi,a2,T,i,b1,g,i,b2,T)%直角坐标下的线性刻度曲线title。图2被辨识参数)figure(3);i=1:200;plot(i,eal,T,i,ea2,g,i,ebl,b,i,eb2,M)%al,a2,bl,b2的各次辨识结果的收敛情况title。图3被辨识参数误差的收敛情况)c(:,k)七、辨识结果分析：1通过辨识仿真说明：广义最小二乘参数辨识精度高，尤其当噪声是有色噪声时。最小二乘法是强烈有偏的，而广义却能得到参数的无偏估计。同时看出Matlab具有强大的运算和分析功能，利用MaUab仿真进行系统辨识，可以大大提高辨识的速度和精度，并且辨识结果直观。实验中所编写的基于Matlab语言编写的通用程序，对解决实际问题具有指导意义。2最小二乘法及其改良算法的性能和优缺点分析。从最小二乘的定义中得知，最小二乘的思想就是寻找一个。的估计值讥使得各次测量的Zja=I与由估计。确定的量测估计Z=Hj)之差的平方和最小，由于此方法兼顾了所有方程的近似程度，使整体误差到达最小，因而对抑制误差是有利的。但是。最小二乘一般方法的估计精度不够高，这是由于对各个测量数据同等对待，而各次测量数据一般不会在相同的条件下获得，造成测量数据的置信度不变较大。因此，最小二乘法有很多改良算法，虽然没有一个是完美的，但是能够适应不同的情况、条件，对应选择不同的算法，其各自的性能及优缺点分析如下：1）根本的最小二乘法（LS,RLSJ对低噪声水平最有效，参数估计可很快收敛到真值，所需计算量相对较少。对于有色噪声的情况只能得到有偏估计。几乎不需要特殊的先验假设，RLS具有可靠的收敛性。2）广义最小二乘法（GLS）一般可给出好的结果，但计算量较大，信噪比拟小时收敛不到真值。可能出现局部最小值的情况。3）多级最小二乘法（MSLS）分三步获得系统模型参数和噪声模型参数的估计，计算量小，速度快，但精度难提高,工程应用不便。4）增广最小二乘法（RELS）算法简单，可同时得到参数和噪声模型的估计，工程应用效果很好。5）加权最小二乘法可对不同置信度的测量值采用加权的方法分别对待，置信度高的，权重取得大些；置信度低的，权重取的小些。但加权最小二乘法仅能用于事先能估计方程误差对参数估计的影响。6）递推最小二乘法可以对计算量大、存储大的数据进行批处理，在线辨识，可以减少数据存储量，防止矩阵求逆，从而减少计算量。矩形窗（限定记忆）RLS方法需要保存一定的数据存储量,此存储量大小取决于矩形窗宽度，因而在应用范围上有一定程度的限制。