CN105808929A - 一种高炉十字测温装置中心温度估计方法 - Google Patents

Abstract

一种高炉十字测温装置中心温度估计方法属于高炉冶炼自动化控制领域，特别涉及一种基于多输出支持向量回归机(M‑SVR)的高炉十字测温装置中心温度估计方法。本发明构造一个具有输出自反馈结构，并考虑上一时刻输入数据的十字测温中心温度估计模型，对十字测温中心五个测温点温度进行在线估计。保证在十字测温中心温度点无法正常测量时，操作人员能够根据模型的温度估计值估计炉内煤气流分布，及时调整布料矩阵，保证高炉稳定顺行。本发明包括以下步骤：(1)模型输入变量的选择；(2)建模策略的确定；(3)模型训练样本的确定。

Description

一种高炉十字测温装置中心温度估计方法

技术领域

本发明属于高炉冶炼自动化控制领域，特别涉及一种基于多输出支持向量回归机(M-SVR)的高炉十字测温装置中心温度估计方法。

背景技术

十字测温作为广泛使用的高炉炉顶温度监测装置，其温度分布是高炉稳定顺行的重要指标。准确地监测十字测温温度，估计炉内煤气流分布，及时调整布料制度，能够辅助高炉稳定顺行，达到高产高效的效果。然而，高炉中心带温度过高，导致十字测温装置传感器的寿命较短，特别是高炉中心的五个测温点(测温点5[T5]、测温点6[T6]、测温点16[T16]、测温点15[T15]、测温点17[T17])传感器极易损坏，并且十字测温传感器更换周期较长，因而无法对炉内温度进行实时测量，严重影响高炉操作人员对高炉内部煤气流分布的判断，导致无法及时的调整布料矩阵，降低了高炉的生产效率。为此，需要建立十字测温装置中心温度估计模型，实时估计十字测温装置中心温度，保证在十字测温中心温度点传感器有损坏时，操作人员可以根据温度的估计值来判断炉内煤气流分布，及时调整布料矩阵。

发明内容

为了解决十字测温装置中心温度点传感器极易损坏，更换周期长，经常不能连续测温的不足，本发明基于文献[M.Sanchez-Fernandez,M.de-Prado-Cumplido,J.Arenas-Garcia and F.Perez-Cruz.SVM multiregression for nonlinear channel estimationin multiple-input multiple-output systems[J].IEEE Transactions on SignalProcessing,2004:106-108.]中提出的多输出支持向量回归机(M-SVR)算法，构造一个具有输出自反馈结构，并考虑上一时刻输入数据的十字测温中心温度估计模型，对十字测温中心五个测温点温度(T5,T6,T16,T15,T17)进行在线估计。保证在十字测温中心温度点无法正常测量时，操作人员能够根据模型的温度估计值估计炉内煤气流分布，及时调整布料矩阵，保证高炉稳定顺行。

为实现上述目的，本发明采用如下技术方案，本发明包括以下步骤：

(1)模型输入变量的选择

需要估计的高炉十字测温中心温度点为T5,T6,T16,T15,T17，根据工艺机理和变量之间的相关性，共选取以下9个变量作为输入：十字测温点3、十字测温点4、十字测温点8、十字测温点10、十字测温点19、顶温东南、顶温西北、顶温西南、顶温西北；

(2)建模策略的确定

考虑到高炉炼铁过程的动态特性以及输入输出之间的时序和滞后关系，引入输出自反馈，并将上一时刻的输入也纳入模型的输入变量中，构成如下23输入5输出估计模型；

$\{{\tilde{y}}_1\left(t\right),\mathrm{...},{\tilde{y}}_5\left(t\right)\}=f\{x_1\left(t\right),\mathrm{...},x_9\left(t\right),x_1(t-T),\mathrm{...},x_9(t-T),{\tilde{y}}_1(t-T),\mathrm{...},{\tilde{y}}_5(t-T)\}---\left(1\right)$

(3)模型训练样本的确定

从厂方数据库分别读入正常生产情况下一段时间的输入和输出样本数据集，根据(2)建模策略的确定中所提建模策略，将上一时刻的输入输出采样数据也纳入模型的输入中，得训练模型的输入变量为x_i＝{x₁(t),...,x₉(t),x₁(t-T),...,x₉(t-T),y₁(t-T),...,y₅(t-T)},输出变量为y_i＝{y₁(t),...,y₅(t)}；

式中，y₁(t),...,y₅(t)为初始训练模型时中心五点温度(T5,T6,T15,T16,T17)当前时刻的实际输出；

首先采用噪声尖峰滤波算法用于剔除高炉生产的噪声尖峰跳变数据；然后采用移动平均滤波算法用于剔除尖峰跳变滤波后训练数据中较小的高频测量噪声波动干扰；

将输入输出数据归一化，缩至[-1,1]区间；

得初始训练模型所需的输入数据集为X＝{x_i|i＝1,2,...,L},输出数据集Y＝{y_i|i＝1,2,...,L}；

式中，x_i＝{x₁(t),...,x₉(t),x₁(t-T),...,x₉(t-T),y₁(t-T),...,y₅(t-T)}，y_i＝{y₁(t),...,y₅(t)}，L为训练样本组数，y₁(t),...,y₅(t)为初始训练模型时中心五点温度(T5,T6,T15,T16,T17)的实际输出；

(4)M-SVR估计模型的训练和使用

(A)开始：初始化所有变量；

(B)如果选择为十字测温温度估计模型训练，转至(C)，读取估计模型训练所需的样本数据集；如果选择为中心温度估计，转至(I)，读取已完成训练的十字测温中心温度估计模型；

(C)读取模型训练数据集：读入(3)模型训练样本的确定中处理完毕的样本输入数据集X和样本输出数据集Y；

(D)M-SVR模型相关待定参数确定：包括目标函数的惩罚因子C，允许偏差ε，核函数以及核函数的相关参数；根据本发明所提出的估计模型的实际情况，核函数选择Sigmoid函数如式(2)；

式中，d为核函数的伸缩变量；

(E)M-SVR模型的训练过程及模型预设参数的确定：

基于模型训练样本集以及(D)中预先设定的模型相关待定参数，惩罚因子C和允许偏差ε进行估计模型训练；本发明所依据的M-SVR算法通过最小化公式(3)中性能指标L_p(W,b)来进行估计模型的训练，求出与最优分类超平面对应的M-SVR估计模型参数向量 $W^{*}={\[w_i^{*}\]}_{i=1}^N,b^{*}={\[b_i^{*}\]}_{i=1}^N,$ 从而得到十字测温中心温度估计模型 ${\tilde{y}}_i=\mathrm{\Φ}\left(x_i\right)W^{*}+b^{*},$ 式中x_i为模型输入数据，为模型输出，N为模型输出向量个数；

$\begin{array}{c}{L}_p(W,b)=\frac{1}{2}\underset{j=1}{\overset{N}{\Σ}}\left|\right|w_i|{|}^2+C\left(\underset{i=1}{\overset{L}{\Σ}}L\right(u_i)+\frac{dL\left(u_i\right)}{du}{|}_{u_i}\frac{{u_i}^2-{\left(u_i\right)}^2}{2u_i})\\ =\frac{1}{2}\underset{j=1}{\overset{N}{\Σ}}\left|\right|w_i|{|}^2+\frac{1}{2}\underset{i=1}{\overset{L}{\Σ}}{a}_i{u_i}^2+CT\end{array}---\left(3\right)$

式中， $u_i=\left|\right|y_i-{\tilde{y}}_i\left|\right|$ 为估计误差， $a_i=\frac{C}{u_i}\frac{dL\left(u_i\right)}{du}{|}_{u_i}=\left\{\begin{array}{cc}0,& u_i<\mathrm{\&epsiv;}\\ \frac{2C(u_i-\mathrm{\&epsiv;})}{u_i},& u_i\&GreaterEqual;\mathrm{\&epsiv;}\end{array}\right.,$ C为目标函数惩罚因子， $L\left(u_i\right)=\left\{\begin{array}{cc}0,& u_i<\mathrm{\&epsiv;}\\ {(u_i-\mathrm{\&epsiv;})}^2,& u_i\&GreaterEqual;\mathrm{\&epsiv;}\end{array}\right.,T=\underset{i=1}{\overset{L}{\&Sigma;}}L\left(u_i\right)-\frac{dL\left(u_i\right)}{du}{|}_{u_i}\frac{u_i^2}{2u_i}$ 为不依赖 $W^{*}={\&lsqb;w_i^{*}\&rsqb;}_{i=1}^N,b^{*}={\&lsqb;b_i^{*}\&rsqb;}_{i=1}^N$ 的常量ε为允许偏差；

采用最小二乘迭代法求解 $W^{*}={\&lsqb;w_i^{*}\&rsqb;}_{i=1}^N,b^{*}={\&lsqb;b_i^{*}\&rsqb;}_{i=1}^N,$ 再依据模型 ${\tilde{y}}_i=\mathrm{\&Phi;}\left(x_i\right)W^{*}+b^{*}$ 得到估计值；

(F)建模效果评估：引入公式(4)相对误差δ(％)对估计模型的建模效果进行评估，如果满足实际工况标准条件，则结束本次M-SVR模型训练过程，转至(G)；如果不符合，则重新训练，转(D)；

$\mathrm{\&delta;}=\frac{1}{L}\underset{i=1}{\overset{L}{\&Sigma;}}\frac{|y_i-{\tilde{y}}_i|}{y_i}\&times;100\%---\left(4\right)$

式中，为模型的估计值，yi为模型的实际输出值，L为样本个数；

(G)保存M-SVR模型：估计模型训练结束，符合实际工况标准条件，得到的M-SVR模型保存，用于高炉十字测温中心温度的估计；

(H)读取M-SVR模型：调出训练完成的高炉十字测温中心温度估计模型 ${\tilde{y}}_i=\mathrm{\&Phi;}\left(x_i\right)W^{*}+b^{*};$

(I)读取估计模型输入的数据：

根据(2)建模策略中所述的输入变量，从厂方数据库中读取相应的输入数据；首次使用模型进行温度估计时，上一时刻温度的估计值使用上一时刻实际值代替；

(J)判断数据是否异常或者缺失；判断模型23个输入数据是否有噪声尖峰跳变以及是否有数据缺失情况；若有则转(K)进行数据处理，否则转(L)进行估计运算；

(K)数据处理：若有噪声尖峰跳变数据，则利用数据库的历史数据进行噪声尖峰滤波处理；若有数据缺失情况，则用上一时刻的相应变量数据进行替换；

(L)温度在线估计运算：将输入变量数据进行归一化处理，调用(G)中保存的M-SVR估计模型进行中心温度在线估计；

(M)十字测温中心温度估计结果显示：在估计系统人机界面上显示本次十字测温中心温度的在线估计结果；

公式中涉及的符号表示含义如下：

当前采样时刻测温点3温度x₁(t)(℃)；

当前采样时刻测温点4温度x₂(t)(℃)；

当前采样时刻测温点8温度x₃(t)(℃)；

当前采样时刻测温点10温度x₄(t)(℃)；

当前采样时刻测温点19温度x₅(t)(℃)；

当前采样时刻东南顶温温度x₆(t)(℃)；

当前采样时刻西北顶温温度x₇(t)(℃)；

当前采样时刻东北顶温温度x₈(t)(℃)；

当前采样时刻西南顶温温度x₉(t)(℃)；

上一采样时刻测温点3温度x₁(t-T)(℃)；

上一采样时刻测温点4温度x₂(t-T)(℃)；

上一采样时刻测温点8温度x₃(t-T)(℃)；

上一采样时刻测温点10温度x₄(t-T)(℃)；

上一采样时刻测温点19温度x₅(t-T)(℃)；

上一采样时刻东南顶温温度x₆(t-T)(℃)；

上一采样时刻西北顶温温度x₇(t-T)(℃)；

上一采样时刻东北顶温温度x₈(t-T)(℃)；

上一采样时刻西南顶温温度x₉(t-T)(℃)；

上一采样时刻测温点5温度y₁(t-T)(℃)；

上一采样时刻测温点6温度y₂(t-T)(℃)；

上一采样时刻测温点16温度y₃(t-T)(℃)；

上一采样时刻测温点15温度y₄(t-T)(℃)；

上一采样时刻测温点17温度_y5(t-T)(℃)；

上一采样时刻测温点5估计值

上一采样时刻测温点6估计值

上一采样时刻测温点16估计值

上一采样时刻测温点15估计值

上一采样时刻测温点17估计值

作为一种优选方案，本发明所述步骤(3)一段时间的输入和输出样本数据集是指：一小时的输入和输出样本数据集。

作为另一种优选方案，本发明所述步骤(3)剔除尖峰跳变滤波后训练数据中较小的高频测量噪声波动干扰是指：剔除尖峰跳变滤波后训练数据中温度值小于10℃的高频测量噪声波动干扰。

作为另一种优选方案，本发明所述式(2)中d＝0.5；式(3)中C＝0.6，ε＝0.1。

作为另一种优选方案，本发明所述采用最小二乘迭代法求解具体步骤如下：

Step1：初始化相关参数，设定k＝0，模型参数值W(k)＝0，b(k)＝0，计算相应的估计误差u_i和支持向量权值a_i；

Step2：根据上一步骤中的初始化条件，由公式(3)求出对应的估计模型参数极值点W_s,b_s，得到线性搜索的下降方向 $Dir\left(k\right)=\left[\begin{array}{c}{W}_s-W\left(k\right)\\ b_s-b\left(k\right)\end{array}\right];$

Step3：确定线性搜索步长λ＝1，并由上一步骤中得到的线性搜索的下降方向Dir(k)，根据 $\left[\begin{array}{c}W(k+1)\\ b(k+1)\end{array}\right]=\left[\begin{array}{c}W\left(k\right)\\ b\left(k\right)\end{array}\right]+\mathrm{\&lambda;}Dir\left(k\right)$ 求出下一步迭代的解W(k+1),b(k+1)；若得到性能指标L_p(W,b)_k+1≥L_p(W,b)_k，则减小搜索步长λ，重新进行 $\left[\begin{array}{c}W(k+1)\\ b(k+1)\end{array}\right]=\left[\begin{array}{c}W\left(k\right)\\ b\left(k\right)\end{array}\right]+\mathrm{\&lambda;}Dir\left(k\right),$ 直到性能指标减小为止，得到下一个迭代点，其中，L_p(W,b)_k+1,L_p(W,b)_k分别表示迭代k+1次和迭代k次的式(3)定义的性能指标取值；

Step4：由迭代解W(k+1),b(k+1)计算模型误差u_i和支持向量权值a_i，令k＝k+1，跳转至Step2，继续执行以上步骤，直到性能指标L_p(W,b)不再减小，跳出循环，此时得到最优的模型参数W^*＝W(k),b^*＝b(k)。

其次，本发明所述满足实际工况标准条件是指：δ≤5％。

另外，本发明还包括以下步骤：

(N)数据保存：将本次估计的相关输入输出数据保存到相应的历史数据库，供之后的估计系统评估、修正以及查询；

(O)在线估计结束？若需继续进行十字测温中心温度的在线估计，则转至(I)；否则转至(P)；

(P)结束。

本发明有益效果。

本发明采用高炉炼铁实际生产数据，并考虑输入输出变量之间的时滞关系，将上一时刻的输入输出数据也纳入模型的输入，应用智能建模算法M-SVR，建立高炉十字测温中心温度估计模型，实现了十字测温中心温度的在线估计。当十字测温装置有损坏时，操作人员可根据本发明估计的温度值判断炉内煤气流分布，及时调整布料矩阵，保证高炉稳定顺行。

附图说明

下面结合附图和具体实施方式对本发明做进一步说明。本发明保护范围不仅局限于以下内容的表述。

图1为高炉十字测温装置示意图。

图2为基于M-SVR的十字测温中心温度估计软件流程框图。

图3为基于M-SVR的十字测温中心温度估计效果图。

图1中：1高炉，2十字测温装置，3顶温测量装置，4数据采集装置，5运行测温软件的计算机系统。

图3所用标记符号如下：

测温点5——T5

测温点6——T6

测温点15——T15

测温点16——T16

测温点17——T17

具体实施方式

如图1所示，十字测温装置和顶温测量装置分别安装在炉喉和炉顶相应位置。数据采集器连接测温系统，并通过通信总线连接运行M-SVR估计软件的计算机系统。

本发明的实现方法包括，(1)模型输入变量的确定(2)建模策略的确定(3)模型训练样本的确定(4)M-SVR估计模型的训练和使用。

为方便描述，对本发明使用的符号和术语定义如下：

当前采样时刻测温点3温度x₁(t)(℃)；

当前采样时刻测温点4温度x₂(t)(℃)；

当前采样时刻测温点8温度x₃(t)(℃)；

当前采样时刻测温点10温度x₄(t)(℃)；

当前采样时刻测温点19温度x₅(t)(℃)；

当前采样时刻东南顶温温度x₆(t)(℃)；

当前采样时刻西北顶温温度x₇(t)(℃)；

当前采样时刻东北顶温温度x₈(t)(℃)；

当前采样时刻西南顶温温度x₉(t)(℃)；

上一采样时刻测温点3温度x₁(t-T)(℃)；

上一采样时刻测温点4温度x₂(t-T)(℃)；

上一采样时刻测温点8温度x₃(t-T)(℃)；

上一采样时刻测温点10温度x₄(t-T)(℃)；

上一采样时刻测温点19温度x₅(t-T)(℃)；

上一采样时刻东南顶温温度x₆(t-T)(℃)；

上一采样时刻西北顶温温度x₇(t-T)(℃)；

上一采样时刻东北顶温温度x₈(t-T)(℃)；

上一采样时刻西南顶温温度x₉(t-T)(℃)；

上一采样时刻测温点5温度y₁(t-T)(℃)；

上一采样时刻测温点6温度y₂(t-T)(℃)；

上一采样时刻测温点16温度y₃(t-T)(℃)；

上一采样时刻测温点15温度y₄(t-T)(℃)；

上一采样时刻测温点17温度y₅(t-T)(℃)；

上一采样时刻测温点5估计值

上一采样时刻测温点6估计值

上一采样时刻测温点16估计值

上一采样时刻测温点15估计值

上一采样时刻测温点17估计值

(1)模型输入变量的确定

如图1所示，需要估计的高炉十字测温中心温度点为十字测温点5、十字测温点6、十字测温点15、十字测温点16、十字测温点17(T5,T6,T15,T16,T17)。根据工艺机理，影响高炉十字测温中心点温度的因素众多，主要包括十字测温相邻温度点温度、顶温、雷达料线、顶压等30个变量，若将这些变量同时作为模型的输入，必然会增加模型的复杂度，降低模型的精度。因此本发明首先对影响十字测温中心温度的相关变量进行约减，选择出与中心五点温度相关性最高并能代表其他变量的作为模型的输入。

本发明选取高炉正常生产情况下5000组数据，包含十字测温中心五点和其余30个变量(采样频率均为10s)，利用多元统计分析对高炉采集到的实际数据进行因子分析和相关性分析。对五个输出变量(T5,T6,T15,T16,T17)进行因子分析得到一个可以代表其他五个输出的主因子，再将得到的主因子与其余30个变量进行相关性分析，选出估计模型的输入变量。考虑工程实际，选取相关性在0.5以上的作为输入变量，共选取以下9个变量作为输入：十字测温点3、十字测温点4、十字测温点8、十字测温点10、十字测温点19、顶温东南、顶温西北、顶温西南、顶温西北。

(2)建模策略的确定

考虑到高炉炼铁过程的动态特性以及输入输出之间的时序和滞后关系，引入输出自反馈，并将上一时刻的输入也纳入模型的输入变量中，构成如下23输入5输出估计模型。

$\{{\tilde{y}}_1\left(t\right),\mathrm{...},{\tilde{y}}_5\left(t\right)\}=f\{x_1\left(t\right),\mathrm{...},x_9\left(t\right),x_1(t-T),\mathrm{...},x_9(t-T),{\tilde{y}}_1(t-T),\mathrm{...},{\tilde{y}}_5(t-T)\}---\left(2\right)$

(3)模型训练样本的确定

根据以往研究，M-SVR在小样本情况下也具有较好的建模和泛化能力，为增加模型的运算速度，减少模型的复杂度，从历史数据中选择正常生产一小时的数据为样本数据进行模型的训练。实际情况中十字测温温度和顶温温度的采样频率均为10s，因此无需对样本进行采样频率的调整。

从厂方数据库分别读入正常生产情况下一小时的输入和输出样本数据集，根据(2)建模策略的确定中所提建模策略，将上一时刻的输入输出采样数据也纳入模型的输入中，得训练模型的输入变量为x_i＝{x₁(t),...,x₉(t),x₁(t-T),...,x₉(t-T),y₁(t-T),...,y₅(t-T)},输出变量为y_i＝{y₁(t),...,y₅(t)}；式中，y₁(t),...,y₅(t)为初始训练模型时中心五点温度(T5,T6,T15,T16,T17)当前时刻的实际输出。

由于测量数据可能存在随机误差，首先采用噪声尖峰滤波算法用于剔除高炉生产的噪声尖峰跳变数据；然后采用移动平均滤波算法用于剔除尖峰跳变滤波后训练数据中较小(温度值小于10(℃))的高频测量噪声波动干扰；

为消除数据量纲对模型的影响，将输入输出数据归一化，缩至[-1,1]区间。

至此，M-SVR模型的训练样本均已处理完毕得初始训练模型所需的输入数据集为X＝{x_i|i＝1,2,...,L},输出数据集Y＝{y_i|i＝1,2,...,L}；式中，x_i＝{x₁(t),...,x₉(t),x₁(t-T),...,x₉(t-T),y₁(t-T),...,y₅(t-T)}，y_i＝{y₁(t),...,y₅(t)}，L为训练样本组数，y₁(t),...,y₅(t)为初始训练模型时中心五点温度(T5,T6,T15,T16,T17)的实际输出。

(4)M-SVR估计模型的训练和使用

本发明的所提方法是以软件的方式实现的，其流程框图如图2所示，分为M-SVR估计模型训练过程和M-SVR估计模型测量过程，其详细步骤如下：

(A)开始：初始化所有变量；

(B)如果选择为十字测温温度估计模型训练，转至(C)，读取估计模型训练所需的样本数据集；如果选择为中心温度估计，转至(I)，读取已完成训练的十字测温中心温度估计模型；

(C)读取模型训练数据集：读入(3)模型训练样本的确定中处理完毕的样本输入数据集X和样本输出数据集Y；

(D)M-SVR模型相关待定参数确定：包括目标函数的惩罚因子C，允许偏差ε，核函数以及核函数的相关参数；根据本发明所提出的估计模型的实际情况，当设定目标函数的惩罚因子C＝0.6，允许偏差ε＝0.1，核函数选择Sigmoid函数(如式2)，核函数参数伸缩变量d＝0.5时，模型的估计精度最好；

式中，d为核函数的伸缩变量。

(E)M-SVR模型的训练过程及模型预设参数的确定：

基于模型训练样本集以及(D)中预先设定的模型相关待定参数，惩罚因子C和允许偏差ε进行估计模型训练；本发明所依据的M-SVR算法通过最小化公式(3)中性能指标L_p(W,b)来进行估计模型的训练，求出与最优分类超平面对应的M-SVR估计模型参数向量 $W^{*}={\&lsqb;w_i^{*}\&rsqb;}_{i=1}^N,b^{*}={\&lsqb;b_i^{*}\&rsqb;}_{i=1}^N,$ 从而得到十字测温中心温度估计模型 ${\tilde{y}}_i=\mathrm{\&Phi;}\left(x_i\right)W^{*}+b^{*},$ 式中x_i为模型输入数据，为模型输出，N为模型输出向量个数。

$\begin{array}{c}{L}_p(W,b)=\frac{1}{2}\underset{j=1}{\overset{N}{\&Sigma;}}\left|\right|w_i|{|}^2+C\left(\underset{i=1}{\overset{L}{\&Sigma;}}L\right(u_i)+\frac{dL\left(u_i\right)}{du}{|}_{u_i}\frac{{u_i}^2-{\left(u_i\right)}^2}{2u_i})\\ =\frac{1}{2}\underset{j=1}{\overset{N}{\&Sigma;}}\left|\right|w_i|{|}^2+\frac{1}{2}\underset{i=1}{\overset{L}{\&Sigma;}}{a}_i{u_i}^2+CT\end{array}---\left(3\right)$

式中， $u_i=\left|\right|y_i-{\tilde{y}}_i\left|\right|$ 为估计误差， $a_i=\frac{C}{u_i}\frac{dL\left(u_i\right)}{du}{|}_{u_i}=\left\{\begin{array}{cc}0,& u_i<\mathrm{\&epsiv;}\\ \frac{2C(u_i-\mathrm{\&epsiv;})}{u_i},& u_i\&GreaterEqual;\mathrm{\&epsiv;}\end{array}\right.,$ C为目标函数惩罚因子， $L\left(u_i\right)=\left\{\begin{array}{cc}0,& u_i<\mathrm{\&epsiv;}\\ {(u_i-\mathrm{\&epsiv;})}^2,& u_i\&GreaterEqual;\mathrm{\&epsiv;}\end{array}\right.,T=\underset{i=1}{\overset{L}{\&Sigma;}}L\left(u_i\right)-\frac{dL\left(u_i\right)}{du}{|}_{u_i}\frac{u_i^2}{2u_i}$ 为不依赖 $W^{*}={\&lsqb;w_i^{*}\&rsqb;}_{i=1}^N,b^{*}={\&lsqb;b_i^{*}\&rsqb;}_{i=1}^N$ 的常量ε为允许偏差。

采用最小二乘迭代法求解具体步骤如下：

Step1：初始化相关参数，设定k＝0，模型参数值W(k)＝0，b(k)＝0，计算相应的估计误差u_i和支持向量权值a_i；

Step2：根据上一步骤中的初始化条件，由公式(3)求出对应的估计模型参数极值点W_s,b_s，得到线性搜索的下降方向 $Dir\left(k\right)=\left[\begin{array}{c}{W}_s-W\left(k\right)\\ b_s-b\left(k\right)\end{array}\right];$

Step3：确定线性搜索步长λ＝1，并由上一步骤中得到的线性搜索的下降方向Dir(k)，根据 $\left[\begin{array}{c}W(k+1)\\ b(k+1)\end{array}\right]=\left[\begin{array}{c}W\left(k\right)\\ b\left(k\right)\end{array}\right]+\mathrm{\&lambda;}Dir\left(k\right)$ 求出下一步迭代的解W(k+1),b(k+1)。若得到性能指标Lp(W,b)_k+1≥L_p(W,b)_k，则减小搜索步长λ，重新进行 $\left[\begin{array}{c}W(k+1)\\ b(k+1)\end{array}\right]=\left[\begin{array}{c}W\left(k\right)\\ b\left(k\right)\end{array}\right]+\mathrm{\&lambda;}Dir\left(k\right),$ 直到性能指标减小为止，得到下一个迭代点，其中，L_p(W,b)_k+1,L_p(W,b)_k分别表示迭代k+1次和迭代k次的式(3)定义的性能指标取值；

Step4：由迭代解W(k+1),b(k+1)计算模型误差u_i和支持向量权值a_i，令k＝k+1，跳转至Step2，继续执行以上步骤，直到性能指标L_p(W,b)不再减小，跳出循环，此时得到最优的模型参数W^*＝W(k),b^*＝b(k)。

再依据模型 ${\tilde{y}}_i=\mathrm{\&Phi;}\left(x_i\right)W^{*}+b^{*}$ 得到估计值；

(F)建模效果评估：引入公式(4)相对误差δ(％)对估计模型的建模效果进行评估，如果满足实际工况标准条件δ≤5％，则结束本次M-SVR模型训练过程，转至(G)；如果不符合，则重新训练，转(D)；

$\mathrm{\&delta;}=\frac{1}{L}\underset{i=1}{\overset{L}{\&Sigma;}}\frac{|y_i-{\tilde{y}}_i|}{y_i}\&times;100\%---\left(4\right)$

式中，为模型的估计值，y_i为模型的实际输出值，L为样本个数。

(G)保存M-SVR模型：估计模型训练结束，符合实际工况标准条件，将得到的M-SVR模型保存，用于高炉十字测温中心温度的估计；

(H)读取M-SVR模型：调出训练好的高炉十字测温中心温度估计模型

(I)读取估计模型输入的数据：

根据(2)建模策略中所述的输入变量，从厂房数据库中读取相应的输入数据。首次使用模型进行估计时，上一时刻温度的估计值使用上一时刻实际值代替；

(J)判断数据是否异常或者缺失；判断模型23个输入数据是否有噪声尖峰跳变以及是否有数据缺失情况；若有则转(K)进行数据处理，否则转(L)进行估计运算；

(K)数据处理：若有噪声尖峰跳变数据，则利用数据库的历史数据进行噪声尖峰滤波处理；若有数据缺失情况，则用上一时刻的相应变量数据进行替换；

(L)温度在线估计运算：将输入变量数据进行归一化处理，调用(G)中保存的M-SVR估计模型进行中心温度在线估计；

(M)十字测温中心温度估计结果显示：在估计系统人机界面上显示本次十字测温中心温度的在线估计结果；

(N)数据保存：将本次估计的相关输入输出数据保存到相应的历史数据库，供之后的估计系统评估、修正以及查询。

(O)在线估计结束？若需继续进行十字测温中心温度的在线估计，则转至(I)；否则转至(P)；

(P)结束。

由于高炉炼铁系统具有显著的非线性动态特性以及易受原料的影响，为了确保估计模型的精度，当温度估计值相对误差大于5％时需要利用新的样本数据对模型进行重新训练。

本发明采用高级语言C#进行估计方法的软件实现。该软件界面实现了数据显示、查询、估计结果显示以及查询等功能，可以方便地让操作人员获得其所需要的信息。另外，安装该估计软件计算机上装有OPC通讯软件负责与下位机以及数据采集装置进行数据双向通讯。

图3为本发明M-SVR估计模型某一小时的中心温度估计效果图，与实测值相比可以看出，本发明估计精确，误差较小，变化趋势基本一致。此外，本发明方法速度快、精度高、泛化能力强且有严格的数学推导，相比于其他方法有较高的优越性。

可以理解的是，以上关于本发明的具体描述，仅用于说明本发明而并非受限于本发明实施例所描述的技术方案，本领域的普通技术人员应当理解，仍然可以对本发明进行修改或等同替换，以达到相同的技术效果；只要满足使用需要，都在本发明的保护范围之内。

Claims (7)

1.一种高炉十字测温装置中心温度估计方法，其特征在于包括以下步骤：

(1)模型输入变量的选择

需要估计的高炉十字测温中心温度点为T5,T6,T16,T15,T17，根据工艺机理和变量之间的相关性，共选取以下9个变量作为输入：十字测温点3、十字测温点4、十字测温点8、十字测温点10、十字测温点19、顶温东南、顶温西北、顶温西南、顶温西北；

(2)建模策略的确定

考虑到高炉炼铁过程的动态特性以及输入输出之间的时序和滞后关系，引入输出自反馈，并将上一时刻的输入也纳入模型的输入变量中，构成如下23输入5输出估计模型；

$\{{\tilde{y}}_1\left(t\right),...,{\tilde{y}}_5\left(t\right)\}=f\{x_1\left(t\right),...,x_9\left(t\right),x_1(t-T),...,x_9(t-T),{\tilde{y}}_1(t-T),...,{\tilde{y}}_5(t-T)\}---\left(1\right)$

(3)模型训练样本的确定

从厂方数据库分别读入正常生产情况下一段时间的输入和输出样本数据集，根据(2)建模策略的确定中所提建模策略，将上一时刻的输入输出采样数据也纳入模型的输入中，得训练模型的输入变量为x_i＝{x₁(t),...,x₉(t),x₁(t-T),...,x₉(t-T),y₁(t-T),...,y₅(t-T)},输出变量为y_i＝{y₁(t),...,y₅(t)}；

式中，y₁(t),...,y₅(t)为初始训练模型时中心五点温度(T5,T6,T15,T16,T17)当前时刻的实际输出；

首先采用噪声尖峰滤波算法用于剔除高炉生产的噪声尖峰跳变数据；然后采用移动平均滤波算法用于剔除尖峰跳变滤波后训练数据中较小的高频测量噪声波动干扰；

将输入输出数据归一化，缩至[-1,1]区间；

得初始训练模型所需的输入数据集为X＝{x_i|i＝1,2,...,L},输出数据集Y＝{y_i|i＝1,2,...,L}；

式中，x_i＝{x₁(t),...,x₉(t),x₁(t-T),...,x₉(t-T),y₁(t-T),...,y₅(t-T)}，y_i＝{y₁(t),...,y₅(t)}，L为训练样本组数，y₁(t),...,y₅(t)为初始训练模型时中心五点温度(T5,T6,T15,T16,T17)的实际输出；

(4)M-SVR估计模型的训练和使用

(A)开始：初始化所有变量；

(B)如果选择为十字测温温度估计模型训练，转至(C)，读取估计模型训练所需的样本数据集；如果选择为中心温度估计，转至(I)，读取已完成训练的十字测温中心温度估计模型；

(C)读取模型训练数据集：读入(3)模型训练样本的确定中处理完毕的样本输入数据集X和样本输出数据集Y；

(D)M-SVR模型相关待定参数确定：包括目标函数的惩罚因子C，允许偏差ε，核函数以及核函数的相关参数；根据本发明所提出的估计模型的实际情况，核函数选择Sigmoid函数如式(2)；

式中，d为核函数的伸缩变量；

(E)M-SVR模型的训练过程及模型预设参数的确定：

基于模型训练样本集以及(D)中预先设定的模型相关待定参数，惩罚因子C和允许偏差ε进行估计模型训练；本发明所依据的M-SVR算法通过最小化公式(3)中性能指标L_p(W,b)来进行估计模型的训练，求出与最优分类超平面对应的M-SVR估计模型参数向量从而得到十字测温中心温度估计模型式中x_i为模型输入数据，为模型输出，N为模型输出向量个数；

$\begin{array}{c}{L}_p(W,b)=\frac{1}{2}\underset{j=1}{\overset{N}{\&Sigma;}}\left|\right|w_i|{|}^2+C({\underset{i=1}{\overset{L}{\&Sigma;}}L\left(u_i\right)+\frac{dL\left(u_i\right)}{du}|}_{u_i}\frac{{u_i}^2-{\left(u_i\right)}^2}{2u_i})\\ =\frac{1}{2}\underset{j=1}{\overset{N}{\&Sigma;}}\left|\right|w_i|{|}^2+\frac{1}{2}\underset{i=1}{\overset{L}{\&Sigma;}}{a}_i{u_i}^2+CT\end{array}---\left(3\right)$

式中， $u_i=\left|\right|y_i-{\tilde{y}}_i\left|\right|$ 为估计误差， $a_i={\frac{C}{u_i}\frac{dL\left(u_i\right)}{du}|}_{u_i}=\left\{\begin{array}{cc}0,& u_i<\mathrm{\&epsiv;}\\ \frac{2C(u_i-\mathrm{\&epsiv;})}{u_i},& u_i\&GreaterEqual;\mathrm{\&epsiv;}\end{array}\right.,$ C为目标函数惩罚因子， $L\left(u_i\right)=\left\{\begin{array}{cc}0,& u_i<\mathrm{\&epsiv;}\\ {(u_i-\mathrm{\&epsiv;})}^2,& u_i\&GreaterEqual;\mathrm{\&epsiv;}\end{array}\right.,T=\underset{i=1}{\overset{L}{\&Sigma;}}L\left(u_i\right)-{\frac{dL\left(u_i\right)}{du}|}_{u_i}\frac{u_i^2}{2u_i}$ 为不依赖 $W^{*}={\&lsqb;w_i^{*}\&rsqb;}_{i=1}^N,b^{*}={\&lsqb;b_i^{*}\&rsqb;}_{i=1}^N$ 的常量ε为允许偏差；

采用最小二乘迭代法求解 $W^{*}={\&lsqb;w_i^{*}\&rsqb;}_{i=1}^N,b^{*}={\&lsqb;b_i^{*}\&rsqb;}_{i=1}^N,$ 再依据模型 ${\tilde{y}}_i=\mathrm{\&Phi;}\left(x_i\right)W^{*}+b^{*}$ 得到估计值；

(F)建模效果评估：引入公式(4)相对误差δ(％)对估计模型的建模效果进行评估，如果满足实际工况标准条件，则结束本次M-SVR模型训练过程，转至(G)；如果不符合，则重新训练，转(D)；

$\mathrm{\&delta;}=\frac{1}{L}\underset{i=1}{\overset{L}{\&Sigma;}}\frac{\left|y_i-{\tilde{y}}_i\right|}{y_i}\&times;100\%---\left(4\right)$

式中，为模型的估计值，y_i为模型的实际输出值，L为样本个数；

(G)保存M-SVR模型：估计模型训练结束，符合实际工况标准条件，得到的M-SVR模型保存，用于高炉十字测温中心温度的估计；

(H)读取M-SVR模型：调出训练完成的高炉十字测温中心温度估计模型 ${\tilde{y}}_i=\mathrm{\&Phi;}\left(x_i\right)W^{*}+b^{*};$

(I)读取估计模型输入的数据：

根据(2)建模策略中所述的输入变量，从厂方数据库中读取相应的输入数据；首次使用模型进行温度估计时，上一时刻温度的估计值使用上一时刻实际值代替；

(J)判断数据是否异常或者缺失；判断模型23个输入数据是否有噪声尖峰跳变以及是否有数据缺失情况；若有则转(K)进行数据处理，否则转(L)进行估计运算；

(K)数据处理：若有噪声尖峰跳变数据，则利用数据库的历史数据进行噪声尖峰滤波处理；若有数据缺失情况，则用上一时刻的相应变量数据进行替换；

(L)温度在线估计运算：将输入变量数据进行归一化处理，调用(G)中保存的M-SVR估计模型进行中心温度在线估计；

(M)十字测温中心温度估计结果显示：在估计系统人机界面上显示本次十字测温中心温度的在线估计结果；

公式中涉及的符号表示含义如下：

当前采样时刻测温点3温度x₁(t)(℃)；

当前采样时刻测温点4温度x₂(t)(℃)；

当前采样时刻测温点8温度x₃(t)(℃)；

当前采样时刻测温点10温度x₄(t)(℃)；

当前采样时刻测温点19温度x₅(t)(℃)；

当前采样时刻东南顶温温度x₆(t)(℃)；

当前采样时刻西北顶温温度x₇(t)(℃)；

当前采样时刻东北顶温温度x₈(t)(℃)；

当前采样时刻西南顶温温度x₉(t)(℃)；

上一采样时刻测温点3温度x₁(t-T)(℃)；

上一采样时刻测温点4温度x₂(t-T)(℃)；

上一采样时刻测温点8温度x₃(t-T)(℃)；

上一采样时刻测温点10温度x₄(t-T)(℃)；

上一采样时刻测温点19温度x₅(t-T)(℃)；

上一采样时刻东南顶温温度x₆(t-T)(℃)；

上一采样时刻西北顶温温度x₇(t-T)(℃)；

上一采样时刻东北顶温温度x₈(t-T)(℃)；

上一采样时刻西南顶温温度x₉(t-T)(℃)；

上一采样时刻测温点5温度y₁(t-T)(℃)；

上一采样时刻测温点6温度y₂(t-T)(℃)；

上一采样时刻测温点16温度y₃(t-T)(℃)；

上一采样时刻测温点15温度y₄(t-T)(℃)；

上一采样时刻测温点17温度y₅(t-T)(℃)；

上一采样时刻测温点5估计值

上一采样时刻测温点6估计值

上一采样时刻测温点16估计值

上一采样时刻测温点15估计值

上一采样时刻测温点17估计值

2.根据权利要求1所述一种高炉十字测温装置中心温度估计方法，其特征在于所述步骤(3)一段时间的输入和输出样本数据集是指：一小时的输入和输出样本数据集。

3.根据权利要求1所述一种高炉十字测温装置中心温度估计方法，其特征在于所述步骤(3)剔除尖峰跳变滤波后训练数据中较小的高频测量噪声波动干扰是指：剔除尖峰跳变滤波后训练数据中温度值小于10℃的高频测量噪声波动干扰。

4.根据权利要求1所述一种高炉十字测温装置中心温度估计方法，其特征在于所述式(2)中d＝0.5；式(3)中C＝0.6，ε＝0.1。

5.根据权利要求1所述一种高炉十字测温装置中心温度估计方法，其特征在于所述采用最小二乘迭代法求解 $W^{*}={\&lsqb;w_i^{*}\&rsqb;}_{i=1}^N,b^{*}={\&lsqb;b_i^{*}\&rsqb;}_{i=1}^N,$ 具体步骤如下：

Step1：初始化相关参数，设定k＝0，模型参数值W(k)＝0，b(k)＝0，计算相应的估计误差u_i和支持向量权值a_i；

Step2：根据上一步骤中的初始化条件，由公式(3)求出对应的估计模型参数极值点W_s,b_s，得到线性搜索的下降方向 $Dir\left(k\right)=\left[\begin{array}{c}{W}_s-W\left(k\right)\\ b_s-b\left(k\right)\end{array}\right];$

Step3：确定线性搜索步长λ＝1，并由上一步骤中得到的线性搜索的下降方向Dir(k)，根据 $\left[\begin{array}{c}W(k+1)\\ b(k+1)\end{array}\right]=\left[\begin{array}{c}W\left(k\right)\\ b\left(k\right)\end{array}\right]+\mathrm{\&lambda;}Dir\left(k\right)$ 求出下一步迭代的解W(k+1),b(k+1)；若得到性能指标L_p(W,b)_k+1≥L_p(W,b)_k，则减小搜索步长λ，重新进行 $\left[\begin{array}{c}W(k+1)\\ b(k+1)\end{array}\right]=\left[\begin{array}{c}W\left(k\right)\\ b\left(k\right)\end{array}\right]+\mathrm{\&lambda;}Dir\left(k\right),$ 直到性能指标减小为止，得到下一个迭代点，其中，L_p(W,b)_k+1,L_p(W,b)_k分别表示迭代k+1次和迭代k次的式(3)定义的性能指标取值；

Step4：由迭代解W(k+1),b(k+1)计算模型误差u_i和支持向量权值a_i，令k＝k+1，跳转至Step2，继续执行以上步骤，直到性能指标L_p(W,b)不再减小，跳出循环，此时得到最优的模型参数W^*＝W(k),b^*＝b(k)。

6.根据权利要求1所述一种高炉十字测温装置中心温度估计方法，其特征在于所述满足实际工况标准条件是指：δ≤5％。

7.根据权利要求1所述一种高炉十字测温装置中心温度估计方法，其特征在于还包括以下步骤：

(N)数据保存：将本次估计的相关输入输出数据保存到相应的历史数据库，供之后的估计系统评估、修正以及查询；

(O)在线估计结束？若需继续进行十字测温中心温度的在线估计，则转至(I)；否则转至(P)；

(P)结束。