广义递推最小二乘辨识.docx
《广义递推最小二乘辨识.docx》由会员分享，可在线阅读，更多相关《广义递推最小二乘辨识.docx（12页珍藏版）》请在冰豆网上搜索。
广义递推最小二乘辨识
系统辨识上机实验报告
广义递推最小二乘辨识
一、实验目的
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）
gridon
title（'图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.4251
0.0116
1.8869
-0.5464
3被辨识参数误差的收敛情况：
六、实验源程序
%递推的广义最小二乘法进行参数估计
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）
gridon
title（'图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,b2
c0=[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;%开始递推运算，求K
yf（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;%求被辨识参数c
e1=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;%完成一次递推，统计值加1
end%大循环结束
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,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:
'）%a1,a2,b1,b2的各次辨识结果的收敛情况
title（'图3被辨识参数误差的收敛情况'）
c（:
k）
七、辨识结果分析：
1通过辨识仿真表明：
广义最小二乘参数辨识精度高，尤其当噪声是有色噪声时。
最小二乘法是强烈有偏的，而广义却能得到参数的无偏估计。
同时看出Matlab具有强大的运算和分析功能，利用Matlab仿真进行系统辨识，可以大大提高辨识的速度和精度，并且辨识结果直观。
实验中所编写的基于Matlab语言编写的通用程序，对解决实际问题具有指导意义。
2最小二乘法及其改进算法的性能和优缺点分析。
从最小二乘的定义中得知，最小二乘的思想就是寻找一个
的估计值
，使得各次测量的
与由估计
确定的量测估计
之差的平方和最小，由于此方法兼顾了所有方程的近似程度，使整体误差达到最小，因而对抑制误差是有利的。
但是。
最小二乘一般方法的估计精度不够高，这是由于对各个测量数据同等对待，而各次测量数据一般不会在相同的条件下获得，造成测量数据的置信度不变较大。
因此，最小二乘法有很多改进算法，虽然没有一个是完美的，但是能够适应不同的情况、条件，对应选择不同的算法，其各自的性能及优缺点分析如下：
1）基本的最小二乘法（LS，RLS）
对低噪声水平最有效，参数估计可很快收敛到真值，所需计算量相对较少。
对于有色噪声的情况只能得到有偏估计。
几乎不需要特殊的先验假设，RLS具有可靠的收敛性。
2）广义最小二乘法（GLS）
一般可给出好的结果，但计算量较大，信噪比较小时收敛不到真值。
可能出现局部最小值的情况。
3）多级最小二乘法（MSLS）
分三步获得系统模型参数和噪声模型参数的估计，计算量小，速度快，但
升级会员 