广义递推最小二乘辨识

系统辨识作业,广义递推最小二乘辨识

系统辨识上机实验报告

广义递推最小二乘辨识

一、实验目的

1 通过实验掌握广义最小二乘辨识算法；

2 运用MATLAB编程，掌握算法实现方法。

二、实验原理

广义最小二乘法的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。如果滤波模型选择得合适，对数据进行了较好的白色化处理，那么直接利用普通最小二乘法就能获得无偏一致估计。广义最小二乘法所用的滤波模型实际上就是一种动态模型，在整个迭代过程中不断靠偏差信息来调整这个滤波模型，使它逐渐逼近于一个较好的滤波模型，以便对数据进行较好的白色化处理，使模型参数估计称为无偏一致估计。理论上说，广义最小二乘法所用的动态模型经过几次迭代调整后，便可对数据进行较好的白化处理，但是，当过程的输出噪信比比较大或模型参数比较多时，这种数据白色化处理的可靠性就会下降。此时，准则函数可能出现多个局部收敛点，因而辨识结果可能使准则函数收敛于局部极小点上而不是全局极小点上。这样，最终的辨识结果往往也会是有偏的。

其收敛速度比较慢，需要经过多次迭代计算，才能得到较准确的参数估计值。一般情况下，经过多次迭代后，估计值便会收敛到稳态值。但在某些情况下（如噪声比较低时）存在局部极小值，估计值不一定收敛到准则函数的全局极小值上。为了防止参数估计值收敛到局部极小值，最好选定初值接近最优解，一般可以用最小二乘法的批处理估计值作为初值。如果系统是时变的，或为了克服数据饱和现象，可以在两次RLS算法中分别引进遗忘因子。

三、实验内容

<1> 数据获取 实验数据按照表10-1，为二阶线性离散系统的输入输出数据

<2> 数据处理 为了提高辨识精度，实验者必须对原始数据进行剔除坏数据、零均值化、工频滤波等处理。实验进行了白化滤波处理。

系统辨识作业,广义递推最小二乘辨识

<3> 辨识算法 利用处理过的数据（取适当的数据长度），选择某种辨识方法（如RLS递推最小二乘法、RELS、RIV或RML等参数估计算法及F-检验或AIC定阶法），估计出模型参数和阶次，同时分析辨识结果。

本实验采用广义递推最小二乘法进行系统辨识。

三、广义递推最小二乘法（RLS）原理

广义最小二乘法是用迭代的松弛算法对最小二乘估计的一种改进，它的基本思想是引入一个白化滤波器，把相关噪声转换为白噪声，基于对观测数据先进行一次滤波处理然后利用普通最小二乘法对滤波后的数据进行辨识。广义最小二乘法的计算步骤如下：

1 给定初始条件：包括给定的输入输出数据或者产生的数据序列，初始状态矩阵P0，被辨识参数的初始值（取一个充分小的实向量），滤波器参数与矩阵初值。

12 利用式zf(k) C(z 1)z(k) uf(k) C( z )u(计算滤波后的输入输出序列。k)

3 对于二阶离散系统，利用式hf(k) [ zf(k 1), zf(k 2),uf(k 1),uf(k 2)]T构造hf(k)。

4 利用

(k) (k 1) Kf(k)[zf(k) hT(k) (k 1)]

Kf(k) P 1)h(k) [T

f1hf(kfk(fP) k(f 1 ( )]1h) k

Pf(k) [I Kf(k)hT

f(k)]Pf(k 1)

三个式子递推计算辨识矩阵 (k)

5 利用式e(k) z(k) hT(k) (k)计算e(k)，并根据he(k) [ e(k 1), e(k 2)]T构造he(k)。

6 利用

e(k) e(k 1) Ke(k)[e(k) hT(k) e(k 1)]

Ke(k) P)he(k)[1 heT(k)P)he(k)] 1 e(k 1e(k 1

TP(k) [I K(k)heee(k)]Pe(k 1)

系统辨识作业,广义递推最小二乘辨识

三个式子递推计算 e(k)。

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.744 0.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.277 0.983, -0.288, 0.846, 1.325, 0.723, 0.713, 0.643 0.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.201 0.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)

grid on

title('图1 输入信号')

y=[1.381 3.794 2.481 -0.280 -2.742 -1.554 2.129 2.691 3.427 2.199 ...

1.679 -1.249 1.371 0.637 3.131 -0.819 0.235 1.262 2.849 3.374...

2.346 0.664 3.015 0.561 2.271 3.650 0.625 2.305 0.364 1.857...

-0.912 -1.547 1.940 0.262 -0.379 -0.176 3.720 0.058 -0.752 1.983...

-0.923 3.361 4.240 -0.074 -0.481 3.780 2.137 0.086 0.638 -0.971...

-0.929 0.679 -0.664 -0.433 1.570 -2.785 -1.153 0.819 3.484 4.091...

-2.375 -2.561 -2.778 2.911 1.362 0.735 3.118 3.770 2.381 -0.812...

-1.635 0.589 1.550 -3.410 -1.249 -3.692 -2.358 2.552 -0.228 0.554...

2.178 2.471 0.743 -0.004 2.504 -3.204 -1.800 -1.284 0.159 0.426...

0.059 0.395 2.371 -0.157 2.248 3.297 2.329 2.780 2.375 1.873...

2.411 3.928 2.846 -2.215 -1.104 3.460 2.883 0.245 -0.231 -2.963...

2.072 -0.845 -0.074 1.037 2.468 -3.679 2.149 -0.081 1.639 -1.291...

2.548 -1.681 2.307 2.227 -1.558 0.008 2.055 -1.102 -1.427 0.350...

系统辨识作业,广义递推最小二乘辨识

2.736 2.965 -2.346 -1.510 0.809 -0.592 2.706 -1.941 2.275 2.802...

-1.337 2.091 -2.585 0.013 1.217 0.691 -0.491 2.114 0.333 -0.482...

3.388 2.082 3.797 4.079 5.036 1.250 -1.019 -0.160 -3.201 1.161...

3.926 1.789 -2.703 2.125 5.054 4.678 5.236 -0.241 2.152 0.356 ...

-3.519 2.213 1.527 -1.206 2.151 0.264 1.595 2.864 -0.539 1.982...

-3.104 -0.264 2.433 0.009 -1.360 -0.521 3.319 1.445 3.105 3.783...

-1.973 -0.138 -0.452 -0.586 -4.045 -1.743 2.577 3.849 0.367 1.324];

<2>初始值设置，包括被辨识参数的初始值、误差序列以及滤波器参数初值；

<3>迭代循环，辨识参数更新，根据误差调整滤波器参数，迭代计算被辨识参数，直至参数符合条件或迭代次数到。

<4>计算数据与图形显示，包括辨识参数辨识过程以及误差的收敛情况。

五、运行结果显示

1 输入数据：

2 被辨识参数：

系统辨识作业,广义递推最小二乘辨识

辨识结果：

ans =

-0.4251

0.0116

1.8869

-0.5464

3 被辨识参数误差的收敛情况：

系统辨识作业,广义递推最小二乘辨识

六、实验源程序

%递推的广义最小二乘法进行参数估计

clear;

close all;

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.744 0.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.277 0.983, -0.288, 0.846, 1.325, 0.723, 0.713, 0.643 0.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.201 0.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)

grid on

title('图1 输入信号')

y=[1.381 3.794 2.481 -0.280 -2.742 -1.554 2.129 2.691 3.427 2.199 ...

1.679 -1.249 1.371 0.637 3.131 -0.819 0.235 1.262 2.849 3.374...

2.346 0.664 3.015 0.561 2.271 3.650 0.625 2.305 0.364 1.857...

-0.912 -1.547 1.940 0.262 -0.379 -0.176 3.720 0.058 -0.752 1.983...

-0.923 3.361 4.240 -0.074 -0.481 3.780 2.137 0.086 0.638 -0.971...

-0.929 0.679 -0.664 -0.433 1.570 -2.785 -1.153 0.819 3.484 4.091...

-2.375 -2.561 -2.778 2.911 1.362 0.735 3.118 3.770 2.381 -0.812...

-1.635 0.589 1.550 -3.410 -1.249 -3.692 -2.358 2.552 -0.228 0.554...

2.178 2.471 0.743 -0.004 2.504 -3.204 -1.800 -1.284 0.159 0.426...

0.059 0.395 2.371 -0.157 2.248 3.297 2.329 2.780 2.375 1.873...

2.411 3.928 2.846 -2.215 -1.104 3.460 2.883 0.245 -0.231 -2.963...

2.072 -0.845 -0.074 1.037 2.468 -3.679 2.149 -0.081 1.639 -1.291...

2.548 -1.681 2.307 2.227 -1.558 0.008 2.055 -1.102 -1.427 0.350...

2.736 2.965 -2.346 -1.510 0.809 -0.592 2.706 -1.941 2.275 2.802...

-1.337 2.091 -2.585 0.013 1.217 0.691 -0.491 2.114 0.333 -0.482...

3.388 2.082 3.797 4.079 5.036 1.250 -1.019 -0.160 -3.201 1.161...

3.926 1.789 -2.703 2.125 5.054 4.678 5.236 -0.241 2.152 0.356 ...

-3.519 2.213 1.527 -1.206 2.151 0.264 1.595 2.864 -0.539 1.982...

-3.104 -0.264 2.433 0.009 -1.360 -0.521 3.319 1.445 3.105 3.783...

-1.973 -0.138 -0.452 -0.586 -4.045 -1.743 2.577 3.849 0.367 1.324];

%用最小二乘递推算法辨识参数：a1,a2,b1,b2

c0=[0.001 0.001 0.001 0.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组数据进行参数估计

for k=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 最小二乘法及其改进算法的性能和优缺点分析。

从最小二乘的定义中得知，最小二乘的思想就是寻找一个 的估计值 ，使

H 得各次测量的Zi(i 1, m)与由估计 确定的量测估计Zii之差的平方和最

小，由于此方法兼顾了所有方程的近似程度，使整体误差达到最小，因而对抑制误差是有利的。但是。最小二乘一般方法的估计精度不够高，这是由于对各个测量数据同等对待，而各次测量数据一般不会在相同的条件下获得，造成测量数据的置信度不变较大。因此，最小二乘法有很多改进算法，虽然没有一个是完美的，但是能够适应不同的情况、条件，对应选择不同的算法，其各自的性能及优缺点分析如下：

1) 基本的最小二乘法（LS，RLS）

对低噪声水平最有效，参数估计可很快收敛到真值，所需计算量相对较少。对于有色噪声的情况只能得到有偏估计。几乎不需要特殊的先验假设，RLS具有可靠的收敛性。

2）广义最小二乘法（GLS）

系统辨识作业,广义递推最小二乘辨识

一般可给出好的结果，但计算量较大，信噪比较小时收敛不到真值。可能出现局部最小值的情况。

3）多级最小二乘法（MSLS）

分三步获得系统模型参数和噪声模型参数的估计，计算量小，速度快，但精度难提高，工程应用不便。

4）增广最小二乘法（RELS）

算法简单，可同时得到参数和噪声模型的估计，工程应用效果很好。

5）加权最小二乘法

可对不同置信度的测量值采用加权的办法分别对待，置信度高的，权重取得大些；置信度低的，权重取的小些。但加权最小二乘法仅能用于事先能估计方程误差 对参数估计的影响。

6）递推最小二乘法

可以对计算量大、存储大的数据进行批处理，在线辨识，可以减少数据存储量，避免矩阵求逆，从而减少计算量。矩形窗（限定记忆）RLS方法需要保留一定的数据存储量，此存储量大小取决于矩形窗宽度，因而在应用范围上有一定程度的限制。