递推最小二乘法的应用

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 =，则计算可得

111111()/(112.36)/12.360.9191/34.12/12.36 2.7670

a T T T

b K T T =-=-=-===

用M 序列作为输入信号，得到输出信号，然后对参数进行辨识。程序如下：

%最小二乘的递推算法

%Z(k+1)=0.9191*Z(k)+2.7670*u(k)+v(k)

%======================================== clear clc

%==========400个产生 M 序列作为输入=============== x=[0 1 0 1 1 0 1 1 1]; %initial value n=403; %n 为脉冲数目 M=[]; %存放 M 序列 for i=1:n

temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j -1); end x(1)=temp; end

%===========产生均值为 0，方差为 1 的高斯白噪声========= v=randn(1,400);

%==============产生观测序列 z================= z=zeros(402,1); z(1)=-1; for i=2:401

z(i)=0.9191*z(i -1)+2.7670*M(i -1)+v(i -1); end %递推求解

P=10*eye(2); %协方差阵

Theta=zeros(2,401); %参数的估计值，存放中间过程估值

Theta(:,1)=[0.001;0.001];%参数估计值初始值a1=0.001,b1=0.001

u=0.99; %遗忘因子

for i=2:401

h=[-z(i-1);M(i-1)];

K=P*h/(h'*P*h+u);%增益矩阵

Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1));

P=(eye(2)-K*h')*P/u;

end

%==========================输出结果及作图============================= disp('参数a1 b1 的估计值:')

Theta(:,401)

a1=zeros(1,401);

a1(1,:)=-0.9191;

b1=zeros(1,401);

b1(1,:)=2.7670;

i=1:401;

figure(1)

plot(i,Theta(1,:),'m',i,a1,'g',i,Theta(2,:),'b',i,b1,'r');

e_a1=Theta(1,:)-a1;% a1的误差

e_b1=Theta(2,:)-b1;% b1的误差

figure(2)

plot(i,e_a1,'m',i,e_b1,'b');%误差曲线图

递推算法仿真曲线如图 1所示：

图 1递推算法仿真曲线

误差曲线如图 2所示：

图 2误差曲线

得到400个数据的a1 b1估计值为:11

ˆˆ0.9208, 2.6564a b =-=由估计值可得1ˆ12.63s T =，1

ˆ33.55K = 由递推算法仿真曲线图 1可知，参数估计收敛速度较快，占用内存较少。由误差曲线图 2可知，a1辨识效果较好，b1辨识误差较大。

1.2 一阶惯性迟延环节

炉内燃烧与传热过程可以简化为磨煤机动态和水冷壁动态两个部分，在此我

们将两者合并为一个一阶惯性迟延环节，研究表明此简化也能较好反映锅炉传热过程：

212()()1

s

K e G s B s T s τ-=⨯+

式中，B 为炉膛的燃料量，为已知量；D 为锅炉总有效吸热量；τ 为纯迟延时间；K2、T2为传递函数的系数。

针对纯迟延的参数应用改进的具有最小损失函数的递推最小二乘法辨识，具体辨识算法为：

12,

12()()()

()[(k 1),,(),(1),

,u(t d )][,,,,,,]

n m y k h k e k h k y y k n u t d m a a a b b b θθ=+=--------=

B 为已知量，令：

22()1

s

K e G s T s τ-=

+ 对τ、K2、T2进行辨识。假设采样时间1s T =，*d T τ=。将以上环节离散化，并写成差分方程的形式

11111111

()[(1)](1)(1)()/,/y k a y k b u k d v k a T T T b K T T =--+--+-=-=

相对于递推最小二乘法 h （k ）中多了一个时间常数 d ，其中，参数 h （k ）是各个参数的函数，输入输出观测向量是纯迟延时间d 的函数，残差 e （k ）取决于各参数的拟合误差，最小二乘估计是使目标函数J （k ，d ）=min ［J （k ，d ）］最小，由于 d 为离散值，当已知纯迟延时间在某一范围［dmin ，dmax ］内时，可以采用最小损失函数的方法将纯延迟时间 d 和其他参数向量一起辨识。具体实现方法可以分2 步来进行：

（1） 假设纯迟延时间已知，利用最小二乘法对其他参数进行估计

ˆˆ(,)min[(,)]J d J d θ

θ=； 例如假设d=2，则

11()[(1)](3)(1)y k a y k b u k v k =--+-+-