递推最小二乘估计及模型阶次辨识

实验二 递推最小二乘估计(RLS)

及模型阶次辨识(F-Test ）

1 实验方案设计

1.1 生成输入数据和噪声

用M 序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。 生成白噪声时，首先利用乘同余法生成U[0,1]均匀分布的随机数，再利用U[0,1]均匀分布的随机数生成标准正态分布的白噪声。 1.2 过程仿真

辨识模型的形式取)()()()()(11k e k u z B k z z A +=--，为方便起见，取n n n b a == 即

n

n n n z

b a b z b z B z a a a z a z A ------++++=++++=...1)(...1)(2

21

12211

用M 序列作为辨识的输入信号。 1.3 递推遗忘因子法

数据长度L 取534，初值????

??

?

????

???

?

??????==1000010000100001)0(001

.0)0(P θ 1.4 计算损失函数、噪声标准差

损失函数?????

?+---+-=μθ

μττ)()1()()]1(?)()([)1()(2k h k P k h k k h k z k J k J

噪声标准差θ

λ

dim )(?-=L L J

1.6 F-Test 定阶法计算模型阶次

统计量t

)22,2(~2

2

2)1()1()()1,(----++-=

+n L F n L n J n J n J n n t

其中，)(?J 为相应阶次下的损失函数值，L 为所用的数据长度，n 为模型

的估计阶次。若a t n n t >+)1,(，拒绝00:n n H >，若a t n n t <+)1,(，接受00:n n H >，其中αt 为风险水平α下的阀值。这时模型的阶次估计值可取1+n 。 1.6 计算噪信比和性能指标

噪信比2

2

y

e σση= 参数估计平方相对偏差i i i n

i i i θθθθθδ?~,~12

21-=???

? ??=∑= 参数估计平方根偏差i

i i n i i

n

i i

θθθθθδ?~,)()

~

(21

2212

2-==

∑∑== 2 编程说明

M 序列中，M 序列循环周期取15124=-=p N ，时钟节拍t ?=1Sec ，幅度1=a ，

特征多项式为1)(56⊕⊕=s s s F 。白噪声循环周期为32768215=。

)(s G 采样时间0T 设为1Sec ，5.0 ,1 ,7.0 ,5.12121===-=b b a a 。

3 源程序清单

3.1 正态分布白噪声生成函数

function v=noise(N)

%生成正态分布N(0,sigma)

%生成N个[0 1]均匀分布随机数

A=179; x0=11; M=2^15;

for k=1:N

x2=A*x0;

x1=mod(x2,M);

v1=x1/(M+1);

v(:,k)=v1;

x0=x1;

end

aipi=v;

sigma=1; %标准差

for k=1:length(aipi)

ksai=0;

for i=1:12

temp=mod(i+k,length(aipi))+1; ksai=ksai+aipi(temp);

end

v(k)=sigma*(ksai-6);

end

end

3.2 M序列生成函数

function [Np r M]=createM(n,a)

%生成长度为n的M序列,周期为Np，周期数为r

x=[1 1 1 1]; %初始化初态

for i=1:n

y=x;

x(2:4)=y(1:3);

x(1)=xor(y(1),y(4));

U(i)=2*y(4)-a;

end

M=U*a;

lenx=length(x);

Np=2^lenx-1;

r=n/Np;

end

3.3 加权最小二乘递推算法函数

function [Aes,Bes,Error]=RLS(na,nb,Z,U,f)

%Aes、Bes为参数估计值，na、nb为模型阶次，Z、U为输出输入数据，f为加权因子

N=na+nb;

n_max=length(Z);

X=0.001.*ones(N,1); %初始估计值

P=10^5.*eye(N); %初始P

e=0.0001; stop=1; %误差要求，循环停止信号

n=N;

Error=zeros(n_max,1);

while(stop==1&&n<=n_max)

H=[]; %新的数据向量

for i=1:na

H=[H;-Z(n-i)];

end

for j=1:nb

H=[H;U(n-j)];

end

K=P*H*inv(H'*P*H+f); %计算增益矩阵

X_past=X;

X=X+K*(Z(n)-H'*X); %计算新的估计值

P=P-K*K'*(H'*P*H+f); %计算下次递推用到的P

temp=abs((X-X_past)./X_past); %相对误差

stop=sum(temp)>=e; %判断精度

Error(n)=Z(n)-H'*X;

n=n+1;

end

Aes=X(1:na)';

Bes=X(na+1:N)';

3.4 主函数

clear %清理工作间变量

L=534; %M序列的周期,四级移位寄存器生成M序列,作为输入信号u(k)

ex=60; %在图像中展示的数据个数

a=1;

aa1=-1.5; aa2=0.7; bb1=1; bb2=0.5; %提前规定的a,b,c,d

[Np r u]=createM(L,a); %生成M序列

figure(1); %画第1个图形：u(k)

stem(u(1:ex)),grid; %以径的形式显示出部分输入信号并给图形加上网格

xlabel('k') %标注横轴变量

ylabel('输入信号') %标注纵轴变量

title(['四级移位寄存器生成M序列输入信号(前',int2str(ex),'位)']) %图形标题z(2)=0;z(1)=0; %取z的前两个初始值为零

y=z;

v=noise(L);% 生成白噪声

lamat=0.1;

for k=3:L; %循环变量从3到L

y(k)=-aa1*z(k-1)-aa2*z(k-2)+bb1*u(k-1)+bb2*u(k-2);

z(k)=y(k)+lamat*v(k); %给出辨识输出采样信号

end

ov=fangcha(v); %计算噪声方差

oy=fangcha(y); %计算信号方差

yita=sqrt(oy\ov); %计算噪信比

%用最小二乘递推算法辨识参数：a,b,c,d

e0=[0.001 0.001 0.001 0.001]'; %被辨识参数的初始值采用直接取方式，取一个充分小的实向量

p0=10^7*eye(4,4); %初始状态P0也采用直接取方式，取一个充分大的实数单位矩阵

E=1e-10; %相对误差E参考值取0.000000005

e=[e0,zeros(4,L-1)]; %被辨识参数矩阵的初始值及大小

eee=zeros(4,L); %相对误差的初始值及大小

n=0; %用于统计递推次数

for k=3:L; %开始递推运算

hk=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; %求h(k)

K=p0*hk*inv(hk'*p0*hk+1); %求K(k)

e1=e0+K*[z(k)-hk'*e0]; %求θ(k)

ee=(K*(z(k)-hk'*e0)); %求相对误差

eee(:,k)=ee; %把当前相对变化的列向量加入误差矩阵的最后一列

e0=e1; %新获得的参数作为下一次递推的旧参数

e(:,k)=e1; %把当前所辨识参数的c1列向量加入辨识参数矩阵的最后一列

pk=p0-K*K'*[hk'*p0*hk+1]; %求p(k)值

p0=pk; %把当前的p(k)值给下次用

n=n+1; %完成一次递推，统计值加1

if ee==E break; %若参数收敛满足要求，终止计算

end%小循环结束

end%大循环结束

a1=e(1,:); a2=e(2,:); b1=e(3,:); b2=e(4,:); ae1=eee(1,:); ae2=eee(2,:); be1=eee(3,:); be2=eee(4,:); %分离参数

figure(2); %画第2个图形

i=1:ex; %横坐标从1到L

plot(i,a1(i),'r',i,a2(i),'m',i,b1(i),'c',i,b2(i),'g') %画出a，b，c，d的各次辨识结果

grid on;

xlabel('k')

ylabel('辨识参数') %标注纵轴变量

title('最小二乘各次递推参数估计值') %图形标题

str1=' ,理论值为：';

legend(['a1',str1,num2str(aa1)],['a2',str1,num2str(aa2)],['b1',st r1,num2str(bb1)],['b2',str1,num2str(bb2)]);

figure(3); %画第3个图形

i=1:ex; %横坐标从1到L

plot(i,ae1(i),'r',i,ae2(i),'g',i,be1(i),'b',i,be2(i),'r:') %画出a,b,c,d的各次辨识结果的收敛情况

grid on;

xlabel('k') %标注横轴变量

ylabel('参数误差') %标注纵轴变量

legend('a1','a2','b1','b2');

title('参数的误差收敛情况') %图形标题

d=1;

Us=u(11:L); %舍弃前10个数据

Zs=v(11:L);

arfa=0.01; %置信度

if(arfa==0.05) %由自由度确定百分点

T_arfa=5.3;

elseif(arfa==0.05)

T_arfa=4.6;

elseif(arfa==0.025)

T_arfa=3.69;

elseif(arfa==0.05)

T_arfa=2.99;

elseif(arfa==0.1)

T_arfa=2.3;

else T_arfa=0.1;

end

norder_max=10; %最大阶数

J=zeros(1,norder_max); %各阶残差方差

norder=d; %起始阶数

stop=0; %终止信号

N=0;

PEs=zeros(norder_max,norder_max); %存储RLS估计结果

while(stop==0&&norder<=norder_max)

[Aes,Bes,Error]=RLS(norder,norder,Zs,Us,0.9); %RLS参数估计

theta=[Aes,Bes]';

thetalength=length(theta);

PEs(1:thetalength,norder)=theta;

J(norder)=Error'*Error; %残差方差

if norder>1

t=((J(norder-1)-J(norder))/J(norder))*((L-2*norder)/2);

if t<=T_arfa

N1=norder-1;

N2=N1-d;

end

end

norder=norder+1;

end

disp(['The estimated order by F_test is ',num2str(N1),',The η is ',num2str(yita)]);

4 曲线打印

图1 驱动序列：M序列

图2 递推函数估计值图像

图3 参数误差收敛图像

5 结果分析

根据输出The estimated order by F_test is 2,The ηis 0.20874 可以得出该模型估计阶数为2，噪信比为0.20874。

6 实验体会

通过本次试验，我不仅学到了噪声生成和系统辨识的方法，也复习了Matlab 编程的方法。在这次试验中，无论是专业课知识还是动手能力上都得到了更高层次提升。

递推最小二乘法推导(RLS)——全网最简单易懂的推导过程

递推最小二乘法推导（RLS）——全网最简单易懂的推导过程 作者：阿Q在江湖 先从一般最小二乘法开始说起 已知x和y的一系列数据，求解参数theta的估计。用矩阵的形式来表达更方便一些： 其中k代表有k组观测到的数据， 表示第i组数据的输入观测量，yi表示第i组数据的输出观测量。令： ，则最小二乘的解很简单， 等价于即参数解为：如果数据是在线的不断的过来，不停的采用最小二乘的解法来解是相当消耗资源与内存的，所

以要有一种递推的形式来保证对的在线更新。 进一步推导出递推最小二乘法（RLS） 我们的目的是从一般最小二乘法的解 推导出 的递推形式。一定要理解这里的下标k代表的意思，是说在有k组数据情况下的预测，所以k比k-1多了一组数据，所以可以用这多来的一组数据来对原本的估计进行修正，这是一个很直观的理解。下面是推导过程： 先看一般最小二乘法的解 下面分别对 和 这两部分进行推导变换，令

得到下面公式（1） 下面来变换得到公式（2） 下面再来，根据一般最小二乘法的解，我们知道下式成立，得到公式（3）（注：后续公式推导用到） 好了，有了上面最主要的三步推导，下面就简单了,将上面推导的结果依次代入公式即可：

至此，终于变成 的形式了。 通过以上推导，我们来总结一下上面RLS方程： 注：以上公式7中，左边其实是根据公式1，右边I为单位矩阵

公式（5）和（7）中，有些文献资料是用右边的方程描述，实际上是等效的，只需稍微变换即可。例如（5）式右边表达式是将公式（1）代入计算的。为简化描述，我们下面还是只讨论左边表达式为例。 上面第7个公式要计算矩阵的逆，求逆过程还是比较复杂，需要用矩阵引逆定理进一步简化。 矩阵引逆定理： 最终RLS的方程解为：

递推最小二乘法算法

题目： （递推最小二乘法） 考虑如下系统： )()4(5.0)3()2(7.0)1(5.1)(k k u k u k y k y k y ξ+-+-=-+-- 式中，)(k ξ为方差为0.1的白噪声。 取初值I P 610)0(=、00=∧ ）（θ。选择方差为1的白噪声作为输入信号)(k u ，采用PLS 法进行参数估计。 Matlab 代码如下： clear all close all L=400; %仿真长度 uk=zeros(4,1); %输入初值：uk(i)表示u(k-i) yk=zeros(2,1); %输出初值 u=randn(L,1); %输入采用白噪声序列 xi=sqrt(0.1)*randn(L,1); %方差为0.1的白噪声序列 theta=[-1.5;0.7;1.0;0.5]; %对象参数真值 thetae_1=zeros(4,1); %（）θ初值 P=10^6*eye(4); %题目要求的初值 for k=1:L phi=[-yk;uk(3:4)]; %400×4矩阵phi 第k 行对应的y(k-1),y(k-2),u(k-3), u(k-4) y(k)=phi'*theta+xi(k); %采集输出数据 %递推最小二乘法的递推公式 K=P*phi/(1+phi'*P*phi); thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1); P=(eye(4)-K*phi')*P; %更新数据 thetae_1=thetae(:,k); for i=4:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=2:-1:2 yk(i)=yk(i-1);

递推阻尼最小二乘法辨识算法公式的详细推导与说明

控制理论与控制工程 学位课程《系统辨识》考试报告 递推阻尼最小二乘法公式详细 推导 专业：控制理论与控制工程 班级：2011双控（研） 学生姓名：江南 学号：20110201016 任课教师：蔡启仲老师 2012年06月29 日

摘要 在参数辨识中，递推最小二乘法是用得最多的一种算法。但是，最小二乘法存在一些缺点，如随着协方差矩阵的减小，易产生参数爆发现象；参数向量和协方差矩阵的处置选择不当会使得辨识过程在参数收敛之前结束；在存在随机噪声的情况下，参数易产生漂移，出现不稳定等。为了防止参数爆发现象，Levenberg 提出在参数优化算法中增加一个阻尼项，以增加算法的稳定性。本文在一般的最小二乘法中增加了阻尼因子，构成了阻尼最小二乘法。又根据实时控制的要求，详细推到了递推阻尼最小二乘公式，实现在线辨识。 关键字：系统辨识，最小二乘法，递推算法 正文 1.题目的基本要求 已知单入单出系统的差分方程以及噪声，在应用最小二乘法进行辨识的时候，在性能指标中加入阻尼因子，详细推导阻尼最小二乘法的递推公式。 2.输入辨识信号和系统噪声的产生方法和理论依据 2.1系统辩识信号输入选择准则 （1）输入信号的功率或副度不宜过大，以免使系统工作在非线性区，但也不应过小，以致信噪比太小，直接影响辩识精度； （2）输入信号对系统的“净扰动”要小，即应使正负向扰动机会几乎均等； （3）工程上要便于实现，成本低。 2.2白噪声及其产生方法 （1） 白噪声过程 （2）白噪声是一种均值为0、谱密度为非0常数的平稳随机过程。 （3）白噪声过程定义：如果随机过程 () t ω的均值为0，自相关函数为 ()()2 R t t ωσδ= （2.2.1） 式中()t δ 为狄拉克（Dirac) 分布函数，即 (){ (),00,0 1t t t dt δδ∞ ∞=≠∞ ==? -且t （2.2.2） 则称该随机过程为白燥声过程。 2.3白噪声序列 (1) 定义 如果随机序列{() }w t 均值为0，并且是两两不相关的，对应的自相关函数为 ()2 ,0,1,2w l R l l σδ==±± 式中{1,0 0,0 l l l δ=≠=则称这种随机序列{()}w t 为白噪声序列。 2.4白噪声序列的产生方法 (1) （0，1）均匀分布随机数的产生 在计算机上产生（0，1）均匀分布随机数的方法很多，其中最简单、最方便的是数学方法。产生伪随机数的数学方法很多，其中最常用的是乘同余法和混合同余法。 ①乘同余法。

用matlab实现最小二乘递推算法辨识系统参数

自动化专业综合设计报告 设计题目：最小二乘递推算法辨识系统参数所在实验室：自动化系统仿真实验室 指导教师： 学生姓名 班级计082-2 班 学号 撰写时间：2012-3-1 成绩评定：

一．设计目的 1、学会用Matlab实现最小二乘法辨识系统参数。 2、进一步熟悉Matlab的界面及基本操作； 3、了解并掌握Matlab中一些函数的作用与使用； 二．设计要求 最小二乘递推算法辨识系统参数，利用matlab编程实现，设初始参数为零。z(k)-1.5*z(k-1)+0.7*z(k-2)=1*u(k-1)+0.5*u(k-2)+v(k); 选择如下形式的辨识模型： z(k)+a1*z(k-1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k); 三．实验程序 m= 3; N=100; uk=rand(1,N); for i=1:N uk(i)=uk(i)*(-1)^(i-1); end yk=zeros(1,N); for k=3:N yk(k)=1.5*yk(k-1)-0.7*yk(k-2)+uk(k-1)+0.5*uk(k-2); end %j=100;kn=0; %y=yk(m:j)'; %psi=[yk(m-1:j-1);yk(m-2:j-2);uk(m-1:j-1);uk(m-2:j-2)]'; %pn=inv(psi'*psi); %theta=(inv(psi'*psi)*psi'*y); theta=[0;0;0;0]; pn=10^6*eye(4); for t=3:N ps=([yk(t-1);yk(t-2);uk(t-1);uk(t-2)]); pn=pn-pn*ps*ps'*pn*(inv(1+ps'*pn*ps)); theta=theta+pn*ps*(yk(t)-ps'*theta); thet=theta'; a1=thet(1); a2=thet(2); b1=thet(3); b2=thet(4); a1t(t)=a1; a2t(t)=a2;b1t(t)=b1;b2t(t)=b2; end t=1:N; plot(t,a1t(t),t,a2t(t),t,b1t(t),t,b2t(t));

应用最小二乘一次完成法和递推最小二乘法算法的系统辨识讲解

1最小二乘法的理论基础 1.1最小二乘法 设单输入单输出线性定长系统的差分方程表示为： 其中δ(k)为服从N(0,1)的随机噪声,现分别测出n+N 个输出输入值y(1),y(2),…,y(n+N),u(1),u(2),…,u(n+N),则可写出N 个方程，写成向量－矩阵形式 (4.1.1) ()()()()()()()() 1201121n n y k a y k a y k a y k n b u k b u k b u k n k ξ=-------+ +-+ +-+()()()()()()101122,,n n a y n n y n a n y b y n N n N b ξξθξξ?? ??++????????????++? ???===??????????????++?????????? ???? ()()()()()()()()() () ()()()() ()( )()()10111121222112n n y n y n y u n u y n y n y u n u y n N y n N y N u n N u N a n a n b n N b ξξξ+--+???? ????+-+-+???? =?????????+-+--+???? ?? ???? ??+?? ??????+??+??????? ???+??????????

则式(1.1.1)可写为 （4.1.2） 式中：y 为N 维输出向量；ξ为N 为维噪声向量；θ为（2n+1）维参数向量；Φ为N ×（2n+1）测量矩阵。因此，式（4.1.1）是一个含有（2n+1）个未知参数，由N 个方程组成的联立方程组。 11y θφφξ--=- 在给定输出向量y 和测量矩阵Φ的条件下求参数θ的估计，这就是系统辨识问题。 设 表示 θ 的估计值，?表示y 的最优估计，则有 （4.1.3） 式中： ()()()10??1??2??,???n n a y n a y n y b y n N b θ???? +????????+????==????????+?????? ???? 设e(k)=y(k)- ?(k), e(k)称为残差，则有e=y- ?=y-Φθ 最小二乘估计要求残差的平方和最小，即按照指数函数 （4.1.4） 求Ｊ对 的偏导数并令其等于0可得： （4.1.5） 由式（4.1.5）可得的 θ 最小二乘估计： （4.1.6） J 为极小值的充分条件是： 即矩阵ΦT Φ为正定矩阵，或者说是非奇异的。 1.1.1最小二乘法估计中的输入信号 当矩阵ΦT Φ的逆阵存在是，式（1.1.6）才有解。一般地，如果u(k)是随机序列或伪随机二位式序列，则矩阵ΦT Φ是非奇异的，即（ΦT Φ）－1存在，式（1.1.6）有解。 现在从ΦT Φ必须正定出发，讨论对u(k)的要求。 y φθξ=+?θ??y θ=Φ()() ??T T J e e y y θ θ==-Φ-Φ?θ() ?20?T J y θ θ ?=-Φ-Φ=??T T y θ ΦΦ=Φ()1 ?T T y θ -=ΦΦΦ220?T J θ ?=ΦΦ>?1 n N yy yu T +-ΦΦ??

(完整word版)多种最小二乘算法分析+算法特点总结

第一部分：程序设计思路、辨识结果分析和算法特点总结 (2) 一：RLS遗忘因子法 (2) RLS遗忘因子法仿真思路和辨识结果 (2) 遗忘因子法的特点： (3) 二：RFF遗忘因子递推算法 (4) 仿真思路和辨识结果 (4) 遗忘因子递推算法的特点： (5) 三：RFM限定记忆法 (5) 仿真思路和辨识结果 (5) RFM限定记忆法的特点： (7) 四：RCLS偏差补偿最小二乘法 (7) 仿真思路和辨识结果 (7) RCLS偏差补偿最小二乘递推算法的特点： (9) 五：增广最小二乘法 (9) 仿真思路和辨识结果 (9) RELS增广最小二乘递推算法的特点： (11) 六：RGLS广义最小二乘法 (12) 仿真思路和辨识结果 (12) RGLS广义最小二乘法的特点： (14) 七：RIV辅助变量法 (14) 仿真思路和辨识结果 (14) RIV辅助变量法的特点： (16) 八：Cor-ls相关最小二乘法（二步法） (17) 仿真思路和辨识结果 (17) Cor-ls相关最小二乘法（二步法）特点： (18) 九：MLS多级最小二乘法 (19) 仿真思路和辨识结果 (19) MLS多级最小二乘法的特点： (22) 十：yule_walker辨识算法 (23) 仿真思路和辨识结果 (23) yule_walker辨识算法的特点： (24) 第二部分：matlab程序 (24) 一：RLS遗忘因子算法程序 (24) 二：RFF遗忘因子递推算法 (26) 三：RFM限定记忆法 (28) 四：RCLS偏差补偿最小二乘递推算法 (31) 五：RELS增广最小二乘的递推算法 (33) 六;RGLS 广义最小二乘的递推算法 (36) 七：Tally辅助变量最小二乘的递推算法 (39) 八：Cor-ls相关最小二乘法（二步法） (42) 九：MLS多级最小二乘法 (45) 十yule_walker辨识算法 (49)

递推最小二乘法

线性方程组的最优求解方法 一．递推最小二乘法 设线性方程组 b Ax = (1) 则有 k b k =x :A ),(, (n k Λ,2,1=) (2) 其中，[]kn k k a a a k ,,,:),(21Λ=A ，[]T n x x x ,,,21Λ=x 。 设 x :A ),()(k k f = (3) 下面采用基于递推最小二乘法(RLS)的神经网络算法来训练权值向量x ，以获得线性方程组（1）的解x 。由式（3）可知，若以)(k f 为神经网络输出，以k b 为神经网络训练样本，以x 为神经网络权值向量，[]kn k k a a a k ,,,:),(21Λ=A 为神经网络输入向量，则解线性方程组的神经网络模型如同1所示。 图1 神经网络模型 采用RLS 算法训练神经网络权值向量x ，其算法如下： （1）神经网络输出： x :A ),()(k k f = (4) （2）误差函数：

)()(k f b k e k -= (5) （3）性能指标： ∑==n k k e J 1 2)(21 (6) （4）使min =J 的权值向量x ，即为所求的神经网络权值向量x ，这是一个多变量线性优化问题，为此，由 0=??x J 可得最小二乘递推法（RLS ）： ]),([1k k k k k k b x :A Q x x -+=+ (7) ),(),(1),(:A P :A :A P Q k k k T k T k k += (8) k k k k P :A Q I P )],([1-=+ (9) ()n k ,,2,1Λ= 随机产生初始权值向量)1,(0 n rand =x ，设n n ?∈=R I P α0（α是足够大的正数（一般取10610~10=α），n n ?∈R I 是单位矩阵），通过对样本数据训练，即可获得神经网络权值 向量x ，此即为线性方程组（1）的解。 二．具有遗忘因子的递推最小二乘估计公式为： ]),([1k k k k k k b x :A Q x x -+=+ (10) ),(),(),(:A P :A :A P Q k k k T k T k k +=λ (11) k k k k P :A Q I P )],([11-= +λ (12) 式中，1:)],(:),([)(-=k A k A k T W P ，W 为加权对角阵：

递推最小二乘法的应用

1 递推最小二乘法在电厂模型辨识中的应用 电厂中大多数热工对象可以用一阶或二阶有迟延和非迟延的模型来表示，对这些模型中参数的辨识，递推最小二乘法是一种较好的方法。本文以火电厂部分典型一阶模型为例子，借助于某电厂现场数据，分别对以下几种环节进行辨识。 1.1 一阶惯性环节 火电厂中，来自锅炉的过热蒸汽，经高压调节汽门和导汽管道进入高压缸膨胀做功，高压缸的排汽回到锅炉再热器被重新加热，加热后的蒸汽经中压调节汽门进入中低压缸进一步膨胀做功，做功后的乏汽最终排入凝汽器变成凝结水，一般中压调节汽门的开度是高压调节汽门的3倍，即在机组负荷大于额定的30%或者滑压运行时，汽轮机的中压调门是完全开启的。因此，在简化模型中，汽机侧调速器一级压力与机组有功功率可以简化为一阶惯性环节如下： 1 11()1 K G s T s = + 将以上环节离散化，并写成差分方程的形式 11111111 ()[(1)](1)(1)()/,/y k a y k b u k v k a T T T b K T T =--+-+-=-= 其中 u 为调速器一级压力，y 为机组有功功率，()v k 为零均值方差为1的高斯白噪 声。 该论文依据递推最小二乘法原理，借助 MATLAB 工具编写程序，设定合适 的初始值和加权因子进行参数辨识，辨识结果为11?? 2.7547,0.9193a b ==-，由11??,a b 可得到134.12K =，112.36s T =进而得到系统的传递函数为： 134.12 ()12.361 G s s = + 下面运用递推最小二乘法对所得结果进行仿真：假设134.12K =，112.36s T =已知，采样时间为1s T =，则计算可得

几种最小二乘法递推算法的小结

一、 递推最小二乘法 递推最小二乘法的一般步骤： 1. 根据输入输出序列列出最小二乘法估计的观测矩阵?： ] )(u ... )1( )( ... )1([)(T b q n k k u n k y k y k ------=? 没有给出输出序列的还要先算出输出序列。 本例中， 2)]-u(k 1),-u(k 2),-1),-y(k -[-y(k )(T =k ?。 2. 给辨识参数θ和协方差阵P 赋初值。一般取0θ=0或者极小的数，取σσ,20I P =特别大，本例中取σ=100。 3. 按照下式计算增益矩阵G ： ) ()1()(1)()1()(k k P k k k P k G T ???-+-= 4. 按照下式计算要辨识的参数θ： )]1(?)()()[()1(?)(?--+-=k k k y k G k k T θ?θθ 5. 按照下式计算新的协方差阵P ： )1()()()1()(---=k P k k G k P k P T ? 6. 计算辨识参数的相对变化量，看是否满足停机准则。如满足，则不再递推；如不满足， 则从第三步开始进行下一次地推，直至满足要求为止。 停机准则：ε???<--) (?)1(?)(?max k k k i i i i 本例中由于递推次数只有三十次，故不需要停机准则。 7. 分离参数：将a 1….a na b 1….b nb 从辨识参数θ中分离出来。 8. 画出被辨识参数θ的各次递推估计值图形。 为了说明噪声对递推最小二乘法结果的影响，程序5-7-2在计算模拟观测值时不加噪 声， 辨识结果为a1 =，a2 = ，b1 = ，b2 =，与真实值a1 =， a2 = ， b1 = ，b2 =相差无几。

几种最小二乘法递推算法的小结

递推最小二乘法的一般步骤： 1. 根据输入输出序列列出最小二乘法估计的观测矩阵?： ] )(u ... )1( )( ... )1([)(T b q n k k u n k y k y k ------=? 没有给出输出序列的还要先算出输出序列。 本例中， 2)]-u(k 1),-u(k 2),-1),-y(k -[-y(k )(T =k ?。 2. 给辨识参数θ和协方差阵P 赋初值。一般取0θ=0或者极小的数，取σσ,20I P =特别 大，本例中取σ=100。 3. 按照下式计算增益矩阵G ： ) ()1()(1)()1()(k k P k k k P k G T ???-+-= 4. 按照下式计算要辨识的参数θ： )]1(?)()()[()1(?)(?--+-=k k k y k G k k T θ?θθ 5. 按照下式计算新的协方差阵P ： )1()()()1()(---=k P k k G k P k P T ? 6. 计算辨识参数的相对变化量，看是否满足停机准则。如满足，则不再递推；如不满足， 则从第三步开始进行下一次地推，直至满足要求为止。 停机准则：ε???<--) (?)1(?)(?max k k k i i i i 本例中由于递推次数只有三十次，故不需要停机准则。 7. 分离参数：将a 1….a na b 1….b nb 从辨识参数θ中分离出来。 8. 画出被辨识参数θ的各次递推估计值图形。 为了说明噪声对递推最小二乘法结果的影响，程序5-7-2在计算模拟观测值时不加噪声， 辨识结果为a1 =1.6417，a2 = 0.7148，b1 = 0.3900，b2 =0.3499，与真实值a1 =1.642， a2 = 0.715， b1 = 0.3900，b2 =0.35相差无几。 程序5-7-2-1在计算模拟观测值时加入了均值为0，方差为0.1的白噪声序列，由于噪声 的影响，此时的结果为变值，但变化范围较小，现任取一组结果作为辨识结果。辨识结果为a1 =1.5371， a2 = 0.6874， b1 = 0.3756，b2 =0.3378。 程序5-7-2-2在计算模拟观测值时加入了有色噪声，有色噪声为 E(k)+1.642E(k-1)+0.715E(k-2)，E(k)是均值为0，方差为0.1的白噪声序列，由于有色噪声的影响，此时的辨识结果变动范围远比白噪声时大，任取一组结果作为辨识结果。辨识结果为a1 =1.6676， a2 = 0.7479， b1 = 0.4254，b2 =0.3965。 可以看出，基本的最小二乘法不适用于有色噪声的场合。