广义递推最小二乘辨识系统辨识实验.docx

广义递推最小二乘辨识系统辨识实验.docx

《广义递推最小二乘辨识系统辨识实验.docx》由会员分享，可在线阅读，更多相关《广义递推最小二乘辨识系统辨识实验.docx（11页珍藏版）》请在冰点文库上搜索。

广义递推最小二乘辨识系统辨识实验

系统辨识上机实验报告

广义递推最小二乘辨识

一、实验目的

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列向量加入辨识参数矩阵的最后一列