广义递推最小二乘辨识

系统辨识上机实验报告广义递推最小二乘辨识一、实验目的1通过实验掌握广义最小二乘辨识算法；2运用MATLAB编程，掌握算法实现方法。二、实验原理广义最小二乘法的根本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。如果滤波模型选择得适宜，对数据进行了较好的白色化处理，那么直接利用普通最小二乘法就能获得无偏一致估计。广义最小二乘法所用的滤波模型实际上就是一种动态模型，在整个迭代过程中不断靠偏差信息来调整这个滤波模型，使它逐渐逼近于一个较好的滤波模型，以便对数据进行较好的白色化处理，使模型参数估计称为无偏一致估计。理论上说，广义最小二乘法所用的动态模型经过几次迭代调整后，便可对数据进行较好的白化处理，但是，当过程的输出噪信比比拟大或模型参数比拟多时，这种数据白色化处理的可靠性就会下降。此时，准那么函数可能出现多个局部收敛点，因而辨识结果可能使准那么函数收敛于局部极小点上而不是全局极小点上。这样，最终的辨识结果往往也会是有偏的。其收敛速度比拟慢，需要经过屡次迭代计算，才能得到较准确的参数估计值。一般情况下，经过屡次迭代后，估计值便会收敛到稳态值。但在某些情况下〔如噪声比拟低时〕存在局部极小值，估计值不一定收敛到准那么函数的全局极小值上。为了防止参数估计值收敛到局部极小值，最好选定初值接近最优解，一般可以用最小二乘法的批处理估计值作为初值。如果系统是时变的，或为了克服数据饱和现象，可以在两次RLS算法中分别引进遗忘因子。三、实验内容<1>数据获取 实验数据按照表10-1，为二阶线性离散系统的输入输出数据<2>数据处理 为了提高辨识精度，实验者必须对原始数据进行剔除坏数据、零均值化、工频滤波等处理。实验进行了白化滤波处理。<3>辨识算法 利用处理过的数据〔取适当的数据长度〕，选择某种辨识方法〔如RLS递推最小二乘法、RELS、RIV或RML等参数估计算法及F-检验或AIC定阶法〕，估计出模型参数和阶次，同时分析辨识结果。本实验采用广义递推最小二乘法进行系统辨识。三、广义递推最小二乘法〔RLS〕原理广义最小二乘法是用迭代的松弛算法对最小二乘估计的一种改良，它的根本思想是引入一个白化滤波器，把相关噪声转换为白噪声，基于对观测数据先进行一次滤波处理然后利用普通最小二乘法对滤波后的数据进行辨识。广义最小二乘法的计算步骤如下：1给定初始条件：包括给定的输入输出数据或者产生的数据序列，初始状态矩阵P0，被辨识参数的初始值〔取一个充分小的实向量〕，滤波器参数与矩阵初值。2利用式计算滤波后的输入输出序列。3对于二阶离散系统，利用式构造。4利用三个式子递推计算辨识矩阵5利用式计算，并根据构造。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(1);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.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(1);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];%用最小二乘递推算法辨识参数：a1,a2,b1,b2c0=[0.0010.0010.0010.001]';%被辨识参数的初始值直接给取，取一个充分小的实向量ce0=[0.001;0.001];%滤波器C(q-1)的参数初值p0=10^6*eye(4,4);%初始状态P0也采用直接取方式，直接给出初始状态P0，取一个充分大的实数单位矩阵pe0=10^6*eye(2,2);%计算滤波器的P矩阵初值yf=zeros(1,200);yf(1)=0;yf(2)=0;uf=zeros(1,200);uf(1)=0;uf(2)=0;e=zeros(1,200);%残差序列e(1)=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)+ce0(1)*y(k-1)+ce0(2)*y(k-2);%滤波后的输入序列uf(k)=u(k)+ce0(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=h1'*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%求出K的值d1=yf(k)-h1'*c0;c1=c0+k1*d1;%求被辨识参数ce1=c1-c0;%求参数当前值与上一次的值的差值e2=e1./c0;%求参数的相对变化e(:,k)=e2;%把当前相对变化的列向量参加误差矩阵的最后一列c0=c1;%新获得的参数作为下一次递推的旧参数c(:,k)=c1;%把当前所辨识参数的c1列向量参加辨识参数矩阵的最后一列p1=p0-k1*k1'*[h1'*p0*h1+1];%求p(k)值p0=p1;%把当前的p(k)值给下次用%%%递推算法的实质就是c和P在原来的根底上计算，因此要为下时刻的这两个值做准备e(k)=y(k)-[-y(k-1),-y(k-2),u(k-1),u(k-2)]*c1;%计算残差e(k),注意这里的h的组成he1=[-e(k-1);-e(k-2)];xe=1+he1'*pe0*he1;xe1=inv(xe);ke1=pe0*he1*xe1;%求当前的K矩阵de1=e(k)-he1'*ce0;ce1=ce0+ke1*de1;ce(:,k)=ce1;%保存滤波器参数ee1=ce1-ce0;ee2=ee1./ce0;%模型参数的相对误差Ee(:,k)=ee1;ce0=ce1;%为下一时刻做准备pe1=pe0-ke1*ke1'*[he1'*pe0*he1+1];%为下一时刻赋值pe0=pe1;%为下一时刻P矩阵赋值n=n+1;%完成一次递推，统计值加1end%大循环结束c%显示被辨识参数e%显示辨识结果的收敛情况n%显示收敛条件满足时算法最终递推次数a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);%别离参数A1=a1(200)A2=a2(200)B1=b1(200)B2=b2(200)figure(2);i=1:200;plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':')%直角坐标下的线性刻度曲线title('图2被辨识参数')figure(3);i=1:200;plot(i,e