自动化学报  2018, Vol. 44 Issue (3): 552-561   PDF    
基于数据驱动多输出ARMAX建模的高炉十字测温中心温度在线估计
周平1, 刘记平1     
1. 东北大学流程工业综合自动化国家重点实验室 沈阳 110819
摘要: 高炉(Blast furnace,BF)炼铁中,十字测温作为炉顶温度和煤气流分布监测的最主要手段,对高炉的安全、稳定和高效运行起着重要作用.然而,由于高炉炉顶中心部位温度较高,造成十字测温装置中心位置传感器极易损坏,并且更换周期长,因而无法及时判断炉顶煤气流分布.针对这一实际工程问题,本文基于时间序列建模思想,集成采用多输出自回归移动平均(Multi-output autoregressive moving average,M-ARMAX)建模、因子分析、Pearson相关分析、基于赤池信息准则(Akaike information criterion,AIC)与模型拟合优度联合定阶等混合技术,提出一种模型结构简单、精度较高且易于工程实现的十字测温中心温度在线估计方法.首先,提出利用因子分析与Pearson相关分析相结合的稳健特征选择方法选取多输出建模输入变量.然后,采用样本均值消去法预处理采集的高炉样本数据,使其成为离散随机数.基于离散随机数,建立算法简单、易于工程实现的M-ARMAX温度模型:为了克服传统基于AIC阶数确定造成模型阶次高、结构复杂的问题,提出在AIC准则基础上进一步引入模型拟合优度来选取模型最小阶,可保证模型估计精度的同时降低模型阶次;同时,采用可快速收敛的递推最小二乘算法辨识M-ARMAX模型参数,并用残差分析方法检验模型.工业试验和比较分析表明:建立的M-ARMAX模型能够根据实时数据同时对十字测温装置多个中心温度点进行准确和稳定估计,且模型估计误差符合高斯白噪声特性.
关键词: 高炉炼铁     十字测温     多输出自回归移动平均建模     温度估计     赤池信息准则     拟合优度     Pearson相关分析    
Data-driven Multi-output ARMAX Modeling for Online Estimation of Central Temperatures for Cross Temperature Measuring in Blast Furnace Ironmaking
ZHOU Ping1, LIU Ji-Ping1     
1. State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University, Shenyang 110819
Manuscript received : December 26, 2016, accepted: May 6, 2017.
Foundation Item: Supported by National Natural Science Foundation of China (61473064, 61290323, 61333007), the Fundamental Research Funds for the Central Universities (N160805001), and the State (Beijing) Key Laboratory of Process Automation in Mining and Metallurgy (BGRIMM-KZSKL-2017-04)
Author brief: LIU Ji-Ping  Master student at Northeastern University. She received her bachelor degree from Henan University of Science and Technology in 2015. Her research interest covers data-driven modeling and control and machine learning algorithm
Corresponding author. ZHOU Ping  Professor at Northeastern University. He received his bachelor, master and Ph. D. degrees from Northeastern University in 2003, 2006, and 2013, respectively. His research interest covers operation feedback control of industrial process, data-driven modeling and control. Corresponding author of this paper
Recommended by Associate Editor XIE Yong-Fang
Abstract: In a blast furnace (BF) ironmaking process, cross temperature measuring is the most important means for monitoring the temperature and gas flow distribution of furnace top. It plays an important role in safe, stable and efficient operation of the whole BF. However, due to the high temperature in the middle of furnace top, the central position sensors of the cross temperature measurement are easily damaged, and the replacement period always takes a long time, thus the gas flow distribution cannot be monitored in time. To solve such a practical engineering problem, a data-driven model with simple structure and high estimation precision is proposed for online estimation of central temperatures for cross temperature measuring. This model integrates hybrid modeling technologies, such as multi-output autoregressive moving average (M-ARMAX) modeling, factor analysis, Pearson correlation analysis, co-determination of model order by Akaike information criterion (AIC) and goodness-of-fit evaluation, etc. The input variables are determined by using factor analysis combined with Pearson correlation analysis with robustness, and the sample mean elimination method is used to preprocess the BF data to make them become discrete time random data. Discrete random data based M-ARMAX modeling includes determination of model order and identification model parameters. The conventional AIC based model order determination leads to too high order of model structure. Thus, the model goodness-of-fit is introduced in the AIC to select the minimum model order, which can not only guarantee the model accuracy but also reduce the model order. Finally, the parameters of M-ARMAX model are identified by recursive least squares algorithm and the obtained models are verified using residual analysis method. Industrial tests and comparative analysis show that the established M-ARMAX model can accurately estimate the center temperatures of cross temperature measuring devices, whose estimated error conforms with Gauss white noise characteristics.
Key words: Blast furnace (BF) ironmaking     cross temperature measurement     multi-output autoregressive moving average (M-ARMAX) modeling     temperature estimation     Akaike information criterion (AIC)     goodness-of-fit     Pearson correlation analysis    

作为钢铁制造最主要的生产环节之一, 高炉冶炼是把固态铁矿石还原成液态生铁的连续生产过程.如图 1所示, 铁矿石、焦炭和熔剂等固体原料按规定配比由炉顶装料装置分批送入高炉, 并使炉喉料面保持一定的高度.焦炭和矿石在炉内形成交替分层结构[1-2].铁矿石在下降过程中逐步被还原、熔化成铁和渣, 聚集在炉缸中, 并从铁口、渣口分别排出.高炉生产中, 炉顶煤气在高炉内的分布直接反映了煤气热能和化学能的利用情况, 因而影响高炉能耗水平和生产成本.然而, 高炉是大型密闭反应器, 且内部生产条件过于严酷, 例如高温高压, 固、液、气三相混合和多场耦合交错, 使得高炉操作人员无法直接观测炉内煤气流的分布及其状态变化.近年来, 十字测温装置在无钟布料高炉得到了广泛应用, 它能连续准确地测出炉喉径向的煤气流温度分布.由于温度高的地方煤气流旺盛, 因而十字测温装置可以有效监测高炉炉顶煤气流的分布状况.十字测温装置位于高炉炉喉位置, 并在炉喉圆周面上的东北、西北、西南、东南方向安装四个测温臂, 每个测温壁分布有不等的温度传感器(例如热电偶), 共有17 $\sim$ 21个测温传感器.这些测温点能够提供实时温度数据, 可比较全面地反映煤气流在炉喉圆周方向上的分布, 使操作者直观地观察到煤气的变化, 为及时准确地进行高炉上部调剂(布料制度调节)和下部调剂(热风和喷煤调节)提供了可靠依据, 同时还避免因煤气取样分析滞后和不全面给高炉操作带来的影响[3-4].

图 1 高炉炼铁过程与十字测温装置 Figure 1 Blast furnace ironmaking process and cross temperature measuring device

典型高炉十字测温装置如图 1所示, 其中西北方向测温臂有6个测温传感器, 其余三个测温臂上有5个测温测温传感器.从图 1某时刻十字测温装置21个测温点实时温度分布图可以看出, 中心五点5, 6, 15, 16, 17温度相对较高, 对传感器的耐高温性要求较高, 因而寿命较短且容易损坏.另外, 十字测温装置安装在密闭的高炉筒体顶部, 其中高温高压和多相多场耦合并行, 粉尘与高温煤气密布, 运行环境极为恶劣.因此, 十字测温装置传感器损坏时不能及时进行维护, 必须等到大周期(一般一个季度一次或半年一次)的高炉检修才能进行传感器的更换.显然, 这给高炉操作人员判断煤气分布带来影响, 并且将大大影响操作员基于十字测温装置的各种操作和决策.因此, 建立十字测温中心温度点的估计模型可以有效解决由传感器损坏带来的各种不利影响.

目前对高炉的研究多集中在高炉炉温的建模或在线估计上, 对高炉十字测温却较少研究且重点放在十字测温装置技术及使用上[5-6].近年来, 人们对支持向量机--神经网路等学习算法的研究热度增加, 因此也有些关于机器学习的算法在高炉十字测温上的研究报导[7-8].但这些已有方法模型结构复杂、精度低、实现困难, 难以在高炉实际应用.本文从实际工程应用的角度出发, 致力于建立一种模型结构简单、精度较高且易于工程实现的十字测温中心温度在线估计方法.从面向在线估计或软测量的建模角度来看, 建立十字测温中心点温度估计模型需要解决如下几方面的问题: 1)如何从众多影响因素中选择有效、可靠且简洁的建模输入变量; 2)如何对高炉实际数据进行预处理, 以得到高质量的建模数据; 3)如何采用有效并且易于实际工程实现的建模算法来进行温度估计模型的建立.

目前数据驱动建模集中在建模方法研究, 而忽略问题1)所述模型输入变量选取即特征选择.实际工程建模中, 模型输入选择大多根据经验主观确定, 容易造成特征冗余及关键特征的丢失, 影响模型精度和建模效率.因此, 针对多输入输出建模问题如何选取有效的输入变量是本文建模首要解决的关键问题之一.为此, 提出将数据降维因子分析和稳健Pearson相关分析的Filter特征选择方法相集成, 综合选取模型输入变量.首先, 采用因子分析方法预处理输出变量, 提取最大主因子, 然后分析主因子与高炉关键过程变量的相关关系, 进而选出模型输入变量. Pearson相关分析具有以下优点: 1)计算速度快, 易于处理大规模数据; 2) Pearson相关系数的取值区间为$[-1 \;\;1], $这个特点使得Pearson相关系数能够捕捉更丰富的变量互信息, 即符号表示关系的正负, 绝对值表示相关性的大小.因此, 所提方法选取输入变量是简单且稳健的方法.

对于问题2), 从高炉实际数据出发, 基于时间序列建模思想, 用样本均值消去法去除高炉数据趋势项提取其随机分量并对随机分量建模[9-10]; 对于问题3), 提出算法结构简单的多输出自回归移动平均(Multi-output autoregressive moving average, M-ARMAX)建模算法, 并用赤池信息准则(Akaike information criterion, AIC)结合拟合优度指标函数选取模型最优阶次, 之后运用收敛速度快的在线递推最小二乘(Recursive least-square, RLS)算法辨识模型参数, 同时采用残差分析方法评价所建模型性能.在问题3)中模型阶数的确定对模型估计精度影响较大.关于模型定阶主要有AIC和贝叶斯信息准则(Bayesian information criterion, BIC)两种方法[9-11]. AIC和BIC均利用模型极大似然估计值、待估参数个数以及样本参数构成指标函数来确定模型阶数[10].但两种方法的缺点是计算量大且耗时, 并且仅使用AIC或BIC准则得到的模型阶数通常较高, 增大模型的复杂度.文献[11]提出一种在线辨识的模型阶估计准则, 利用RLS辨识系统阶数和参数, 将最小方差函数作为指标函数, 推导出模型阶次估计的递推形式.但该方法比AIC和BIC准则更为复杂, 并且结果相差不大.为此本文在AIC准则基础上, 进一步引入模型拟合优度指标函数进行联合阶数.所提方法能够在保证模型估计精度的同时降低模型阶次.

考虑到十字测温中心五点的温度不仅与输入变量有依存关系, 还与其自身历史数据存在时序上的相关关系, 同时考虑高炉内部随机干扰的存在, 为此本文选用简单有效的ARMAX时间序列建模技术来建立十字测温中心五点的动态温度模型.由于温度模型输出为中心五点温度, 并且炉喉中心5点温度输出变量之间也具有相关性, 因此建模还需考虑输出变量之间的相关关系.为此, 建立多输出ARMAX (M-ARMAX)十字测温中心五点温度模型.基于实际工业数据的实验和比较分析表明, 相比于多输出支持向量机(Multi-output support vector machine, M-SVM)算法、偏最小二乘回归(Partial least squares regression, PLSR)算法、动态PLSR (DPLSR)算法和BP神经网络算法, 本文建立的基于M-ARMAX的十字测温中心温度模型简单有效, 且具有较高的精度, 可以很好地用于高炉实际生产, 为高炉操作人员判断煤气流在炉喉分布和相关操作决策提供依据.

1 建模算法 1.1 建模输入变量选取

模型输入选择对模型建模性能具有重要影响.考虑到传统单一Pearson相关分析和单一因子分析在多输入多输出建模输入变量选择的不足(例如Pearson相关分析选出的输入变量可能与多输出变量中一些变量具有强相关性, 但是与其他输出变量却不相关), 本文针对研究的多输入多输出时间序列建模问题, 提出基于数据降维因子分析与Pearson相关分析相结合的稳健特征选择新方法.首先采用因子分析从多个输出变量中找出一个主因子, 并与众多过程输入变量做Pearson相关分析, 选出与主因子相关性较强的过程变量作为模型输入变量.

因子分析是一种分析多变量之间的依赖关系并降维的数据处理技术.因子分析提取的主因子能够表示原始众多变量方差变化最大的方向.十字测温中心带的温度由于所在空间位置以及温度本身的连续特性, 使得从十字测温中心带五个输出变量提取一个主因子成为可能.因子分析提取的主因子不可直接观测, 但又客观存在, 是众多原始变量的共同影响因素, 每一个变量都可表示成主因子的线性函数与特殊因子之和, 即

$ x_i=a_{i1}f_1+a_{i2}f_2+\cdots+a_{ip}f_p+\varepsilon_i $ (1)

式中, ${i=1, 2, \cdots, m}$, ${x_i}$表示第${i}$个原始变量, 共有${m}$个原始变量. ${f_i}$表示第${i}$个主因子变量, ${p}$为主因子个数.式(1)的向量表达式为

$ \pmb{x}=A \pmb{f}+\pmb{\varepsilon} $ (2)

式中, $\pmb{x}=[x_1 x_2 \cdots x_m]^{\rm{T}}$, $\pmb{f}=[f_1 f_2 \cdots f_p]^{\rm{T}}$, $A= \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1p}\\ a_{21} & a_{22} & \cdots & a_{2p}\\ \vdots & \vdots & \ddots &\vdots \\ a_{m1} & a_{m2} & \cdots & a_{mp} \end{bmatrix}.$式(2)中的矩阵$A$为因子载荷矩阵, 可用主成分分析(Principal component analysis, PCA)技术进行求解, 具体算法和步骤如下:

算法1. 基于PCA的因子载荷矩阵求取

步骤1.计算原始变量数据标准化之后的协方差阵.

步骤2.计算原始变量协方差阵$\Sigma$的特征值$\lambda_1\geq\lambda_2$ $>$ $\cdots >\lambda_m$, 相应的特征向量为$\pmb{t}_1, \pmb{t}_2, \cdots, \pmb{t}_m$, 选取大于1特征值$\lambda_1\geq\lambda_2>\cdots>\lambda_p$及其对应的特征向量$\pmb{t}_1, \pmb{t}_2, \cdots, \pmb{t}_p$.

步骤3.利用下式计算因子载荷矩阵

$ A=[\lambda_1\pmb{t}_1, \lambda_2\pmb{t}_2, \cdots, \lambda_P\pmb{t}_P] $ (3)

求取因子载荷矩阵$A$后, 利用Thomson回归可得出主因子估计值的计算公式为: $\hat{\pmb{f}}={A^{\rm T}}\sum^{-1}\pmb{x}$, 即可估计出主因子的值.另外, 计算出输出变量的主因子之后, 需利用主因子与众多影响因子做Pearson相关分析, 以剔除与主因子不相关的影响因子, 进一步减少模型的维数. Pearson相关分析又称为皮尔逊积矩相关系数(Pearson moment correlation coefficient, PMCC)用来度量两个变量之间的相互关系, 计算公式如下:

$ r=\frac{\sum\limits_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})} {\sqrt{{\sum\limits_{i=1}^{n}(x_i-\bar{x})^2}} \sqrt{{\sum\limits_{i=1}^{n}(y_i-\bar{y})^2}}} $ (4)

式中, $x_i$, $y_i$表示第$i$时刻两个变量的采样值, $n$为采样数.

1.2 数据预处理

高炉内部干扰因素较多, 数据采集与传送的误差等都会影响数据的质量, 因此做好数据预处理是建模的前提.数据预处理包括差分、归一化、样本均值消去、缺省值及垃圾数的处理等.本文针对研究的时间序列建模问题, 首先采用逐样本均值消减法预处理高炉原始数据, 消去其直流成分, 并对其随机成分进行建模.这样可很大程度保证基于时间序列建模的最优性能.

1.3 M-ARMAX建模算法

考虑到十字测温各测量点温度的时序特性, 本文采用时间序列分析方法建立高炉十字测温装置中心五点的温度估计模型.时间序列分析的基本思想是根据观察到的历史数据, 建立能够比较精确地反映时间序列中所包含动态依存关系的数学模型, 来评价事物的现状和估计事物的未来变化, 并用此模型对系统的未来行为进行预测[12].所以可分析从高炉中获得的输入输出时间序列数据之间以及输出时间序列之间的相互关系来建立十字测温中心点温度的估计模型.时间序列数据的主要特点是存在趋势项, 需要将其从原序列中分解出来, 得到由各种因素影响的随机项, 并对随机项进行建模, 对于高炉数据趋势项即为数据的直流分量(DC competent).离散时间随机数线性动态模型的一般结构为

$ A(z^{-1})y(k)=\frac{B(z^{-1})}{F(z^{-1})}u(k) +\frac{C(z^{-1})}{D(z^{-1})}\varepsilon(k) $ (5)

式中,

$ A(z^{-1})=1+a_1z^{-1}+a_2z^{-2}+\cdots+a_{np}z^{-np}\\ B(z^{-1})=1+b_1z^{-1}+b_2z^{-2}+\cdots+b_{nq}z^{-nq}\\ F(z^{-1})=1+f_1z^{-1}+f_2z^{-2}+\cdots+f_{nf}z^{-nf}\\ C(z^{-1})=1+c_1z^{-1}+c_2z^{-2}+\cdots+c_{nc}z^{-nc}\\ D(z^{-1})=1+d_1z^{-1}+d_2z^{-2}+\cdots+d_{nd}z^{-nd} $

$z^{-1}$为延迟或后移算子, 即$z^{-1}y(k)=y(k-1)$, $\varepsilon(k)$为零均值的高斯白噪声, $np$, $nq$, $nf$, $nc$, $nd$为滞后的阶次. $k$为采样时刻.

此外, 当$F(z^{-1})=D(z^{-1})=1$时, 式(5)为ARMAX模型; 当$F(z^{-1})=D(z^{-1})=C(z^{-1})=1$时, 式(5)为自回归各态历经(Autoregressive exogenous, ARX)模型; 当$A(z^{-1})=D(z^{-1})=C(z^{-1})=1$时, 式(5)为输出误差(Output error, OE)模型; 当$A(z^{-1})=1$时, 式(5)为Box-Jenkins (BJ)模型.

对于多输入多输出系统, $A(z^{-1})$, $B(z^{-1})$, $C(z^{-1})$分别为维数是$N_y\times N_y$, $N_y\times N_u$, $N_y\times 1$的多项式矩阵, $N_y$为输入的维数, $N_u$为输出的维数.考虑到本文需要建立十字测温中心五点的多输出温度估计模型, 以及中心五点温度与各输入的时序和时滞等动态关系, 本文从实际工程易于实现实现的角度出发, 提出建立多输出ARMAX(M-ARMAX)的十字测温中心温度估计模型.

1.3.1 基于AIC与模型拟合优度的M-ARMAX模型定阶

建模过程增加模型阶数可以提高模型精度, 但是模型复杂度会增大, 并易于导致过拟合问题.为此, 采用信息准则作为确定模型阶数的依据.信息准则通过加入模型复杂度的惩罚项来避免模型过拟合问题.最常用信息准则为赤池信息准则(AIC) [13-15]. AIC是衡量统计模型优良性的一种标准, 由日本统计学家赤池弘次在1974年提出, 它建立在熵的概念上, 提供了权衡估计模型复杂度和拟合数据优良性的标准.通常情况下, AIC定义为

$ {\rm{AIC}}=\log({v})+\frac{2d}{N} $ (6)
$ {v}=\det\left(\frac{1}{N}\sum\limits_{t=1}^N\varepsilon(t, \theta_N)(\varepsilon(t, \theta_N))^ {\rm{T}}\right) $ (7)

式(6)中, $d$为模型参数个数, $v$为损失函数, 表达式如式(7)所示, $N$为训练样本容量.通过改变M-ARMAX模型中$np$$nq$的值, 可选择使AIC最小的模型阶数.但是单纯使用AIC定阶, 随着模型阶次的增加, 模型AIC值降低的幅度变化不大, 会最终导致模型阶次仍然较高.为此, 本文在AIC准则基础上, 进一步引入模型拟合优度目标函数来综合选取模型阶数.这可以在确保模型精度的同时减少模型阶次.基于AIC与模型拟合优度的模型定阶算法和步骤如下:

算法2.基于AIC与模型拟合优度的模型定阶算法

步骤1.采用AIC定阶准则进行阶数初步确定.

步骤2.求取模型AIC输入输出阶次三维图中各等势线上的阶次组合$np$$nq$.

步骤3.利用下式计算不同阶次组合M-ARMAX模型的拟合优度数值

$ fit=100\times \left(1-\frac{\|y-\hat{y}\|}{\|y-{\rm mean}(y)\|}\right) $ (8)

步骤4.将拟合优度最高值确定为最优阶次组合.

1.3.2 基于RLS的M-ARMAX模型参数辨识

需要建立的温度估计模型M-ARMAX表达式为

$ y_i(k)= A_i(z^{-1})y_i(k)+\sum\limits_{j=1, j\neq i}^{N_y}-A_j(z^{-1})y_j (k)~+\\ B_i(z^{-1}){\pmb{u}}(k) +C_i(z^{-1})\varepsilon _i(k) $ (9)

式中, $i$, $j=1, 2, \cdots, 5$, $A_i(z^{-1})$为首一多项式, $B_i(z^{-1})$$1$ $\times$ $N_u$多项式矩阵, ${\pmb{u}}(k)=[u_1(k), u_2(k), \cdots, u_{7}(k)]$为输入变量矩阵, $C_i(z^{-1})$为首一多项式, $\varepsilon _i(k)$为零均值高斯白噪声序列.由于建模数据为可实时采集的高炉数据, 包含各种未知动态干扰, 所以式(9)所示M-ARMAX温度模型可简化为

$ y_i(k)= A_i(z^{-1})y_i(k)+ \sum\limits_{j=1, j\neq i}^{N_y}-A_j(z^{-1})y_j(k)~+\\ B_i(z^{-1}){\pmb{u}}(k) $ (10)

定义数据向量和参数向量:

$ \pmb\varphi(k-1)=[\pmb{y}_1, \pmb{y}_2, \cdots, \pmb{y}_{N_y}] $ (11)
$ \pmb\theta^{\rm{T}}=[\pmb{a}_1, \pmb{a}_2, \cdots, \pmb{a}_{N_y}, \pmb{b}_1, \pmb{b}_2, \cdots, \pmb{b}_{N_u}] $ (12)

式中,

$ \pmb{y}_i=[-y_i(k-1), \cdots, -y_i(k-np)], \quad i=1, 2, \cdots, N_y\\ \pmb{u}_j=[u_j(k), \cdots, u_j(k-nq)], \quad j=1, 2, \cdots, N_u\\ \pmb{a}_i=[a_{i1}, a_{i2}, \cdots, a_{i\times np}], \quad i=1, 2, \cdots, N_y\\ \pmb{b}_j=[b_{i1}, b_{i2}, \cdots, b_{i\times nq}], \quad j=1, 2, \cdots, N_u $

将式(11)和式(12)代入式(10), 可得

$ y(k)=\pmb\varphi(k-1)\pmb\theta $ (13)

因此, 上述问题可用最小二乘(Least square, LS)一次算法来辨识参数, 参数辨识表达式为

$ \pmb\theta=(\pmb\varphi^{\rm{T}}(k-1)\pmb\varphi(k-1))^{-1}\pmb\varphi^{\rm{T}}(k-1)y(k) $ (14)

引入记号: $\pmb{y}_N= \begin{bmatrix} y(1)\\ y(2)\\ \vdots\\ y(N) \end{bmatrix} $, $\Phi_N= \begin{bmatrix} \pmb{\varphi}(0)\\ \pmb{\varphi}(1)\\ \vdots\\ \pmb{\varphi}(N-1) \end{bmatrix} $, 则参数表达式可转化为: $\pmb{\theta}=(\Phi^{\rm{T}}_N\Phi_N)^{-1}\Phi^{\rm{T}}_N\pmb{y}_N$, 从上述表达式可以看出, LS一次算法需要计算逆矩阵, 对于多变量高维数的数据向量, LS的计算量往往很大.为此采用递推最小二乘(RLS)算法消除矩阵求逆计算.这不但可大大提高计算效率, 同时还可根据新的观测数据不断修正参数估计值直至参数精度达到满足.实际上, RLS基于LS一次算法, 在每次取得一次新的观测数据后, 就在旧的参数估计值基础上, 利用新引入的观测数据对前一时刻的参数估计值根据递推算法进行修正从而得到当前时刻参数估计值.随着观测数据的不断引入, 参数估计值不断得到修正直到达到满意的精度为止[16-18].

$\pmb{\theta}(k-1)$是基于时刻$k-1$为止的所有观测数据在$k$ $-$ $1$时刻的未知参数$\pmb{\theta}$的LS估计, 为采集到一组新数据后得到的$k$时刻的参数估计, 根据LS一次算法可得

$ \pmb{\theta}(k-1)=(\Phi^{\rm{T}}_{k-1}\Phi_{k-1})^{-1}\Phi^{\rm{T}}_{k-1}\pmb{y}_{k-1} $ (15)
$ \pmb{\theta}(k)=(\Phi^{\rm{T}}_{k}\Phi_{k})^{-1}\Phi^{\rm{T}}_{k}\pmb{y}_{k} $ (16)

考虑到$ \Phi_{k}= \begin{bmatrix} \Phi_{k-1}\\ \pmb\varphi(k-1) \end{bmatrix} $, $ \pmb{y}_{k}= \begin{bmatrix} \pmb{y}_{k-1}\\ y(k) \end{bmatrix} $, 结合式(15)和式(16), 则

$ \pmb{\theta}(k)= (\Phi^{\rm{T}}_{k-1}\Phi_{k-1}+\pmb\varphi^{\rm{T}}(k-1)\pmb\varphi(k-1))^{-1}~\times \notag \\ (\Phi^{\rm{T}}_{k-1}\pmb{y}_{k-1}+\pmb\varphi^{\rm{T}}(k-1)y(k)) $ (17)

$P(k)=(\Phi^{\rm{T}}_{k}\Phi_{k})^{-1}$, 则有$P(k -1)=(\Phi^{\rm{T}}_{k-1}\Phi_{k-1})^{-1}$, 于是有$P(k)=(P^{-1}(k-1)+\pmb\varphi^{\rm{T}}(k-1)\pmb\varphi(k-1))^{-1}$.根据著名的矩阵求逆引理[13], 可得最终RLS算法公式为

$ K(k)=\frac{P(k-1)\pmb\varphi(k-1)}{1+\pmb\varphi^{\rm{T}}(k-1)P(k-1)\pmb\varphi(k-1)} $ (18)
$ \pmb{\theta}(k)=\pmb{\theta}(k-1) +K(k)(y(k)-\pmb\varphi^{\rm{T}}(k-1)\pmb{\theta}(k-1)) $ (19)
$ P(k)=(I+K(k)\pmb\varphi^{\rm{T}}(k-1)P(k-1)) $ (20)

式中, $K(k)$, $\pmb{\theta}(k)$, $P(k)$分别为增益矩阵、参数向量和逆矩阵.注意到上述RLS算法需要事先给定初值$\pmb{\theta}(0)=0$, $P(0)$ $=$ $\alpha I$, 其中$\alpha$为足够大的正数, 本文选用$10^6$. RLS算法具体实施步骤为:

算法3.基于RLS的M-ARMAX模型参数辨识算法

步骤1.设置初值$\pmb{\theta}(0)=0$, $P(0)=10^6I$;

步骤2.构造数据向量$\pmb\varphi(k-1)$;

步骤3.根据式(17)计算增益矩阵$K(k)$;

步骤4.根据式(18)计算参数估计向量$\pmb{\theta}(k)$;

步骤5.根据式(19)计算逆矩阵矩阵$P(k)$;

步骤6.递推一步$k-1\rightarrow k$, 返回步骤2, 构造新的数据向量.

参数估计完成后, 还需要进行模型诊断以检验模型是否过拟合, 所以在检验模型拟合优度的同时还需要检验模型残差是否为白噪声, 残差的白噪声性质可从残差自相关函数和高斯曲线看出[19-20].

2 工业实验及结果分析

为了验证所提方法, 选取某钢厂2 650 ${\rm{m}}^3$高炉实际生产数据.此数据为采样频率为10 s的离散时间序列数据, 通过样本均值消去法预处理这些数据使其成为离散时间随机数.建模输出变量为需要估计的十字测温装置中心测温点5, 6, 15, 16, 17的温度.建模数据确定和预处理后, 首先采用本文方法提取十字测温中心五点温度建模的主因子. 表 1为数据的Kaiser-Meyer-Olkin (KMO)与Bartlett检验表, 从表 1可以看出, ${\rm{KMO}}>0.6$, Bartlett显著性检验为0, 因此十字测温中心五点适合做因子分析.图 2为中心五点协方差阵的碎石图, 可以看出协方差阵$\Sigma$特征值大于1的因子有1个, 因此可以从5个输出变量中提取一个具有代表性的主因子, 与输入变量进行Pearson相关分析, 选出与5个输出变量相关性较高的过程变量. 表 2为通过PCA计算得到的该因子载荷矩阵${A}$值.计算出输出变量的主因子之后, 利用这个主因子与众多高炉过程变量做Pearson相关分析, 结果如图 3所示.为了进一步减少模型的维数, 选出与主因子相关性较强的过程变量作为模型的输入变量, 本文选取相关性大于$r_0=0.5$的过程变量为模型输入.得到建模输入变量为:十字测温点3、十字测温点4、十字测温点8、十字测温点10、十字测温点20、顶温东南、顶温西北、顶温西南、顶温西北.因子分析所选的主因子是代表的5个输出变量方差变化最大的方向, 没有实际物理意义, 并且这一主因子并不能解释五个输出变量的所有方差信息.为此将主因子选出的各个影响因子再与十字测温中心五个温度点做Pearson相关性分析, 以剔除各变量与中心温度点不相关的影响因子, 结果如表 3所示.从表 3可以看出, 十字测温点4和十字测温点20与中心测温点6不相关, 十字测点20与中心测温点16不相关, 所以将十字测温点4和十字测温点20剔除.最后所选建模输入变量为十字测温点3、十字测温点8、十字测温点10、顶温东南、顶温西北、顶温西南、顶温西北.

图 2 碎石图 Figure 2 Scree plot
图 3 所有输入与主因子的Pearson相关系数 Figure 3 Pearson correlation coefficients between all inputs and the main factor
表 1 KMO和Bartlett分析结果 Table 1 The results of KMO and Bartlett analysis
表 2 因子载荷矩阵 Table 2 Factor loading matrix
表 3 输入输出变量的Pearson相关系数 Table 3 Pearson correlation coefficients between inputs and outputs

然后, 采用本文方法对模型阶数进行确定.模型AIC值与不同输入输出阶次的三维立体图如图 4所示, 其中等势线上不同阶次组合对应的模型拟合优度如表 4所示.从表 4可以看出, 模型拟合优度函数最大值为96.7, 对应的M-ARMAX模型阶数为.最后, 选用RLS算法进行M-ARMAX模型参数估计, 参数收敛曲线如图 5所示.从图 5可以看出, 参数经过200次迭代后, 均基本能够收敛.

图 4 不同阶次所对应的模型AIC值 Figure 4 The AIC value corresponding to different order
图 5 参数收敛曲线 Figure 5 Parameter convergence curves
表 4 不同阶次对应的模型拟合优度 Table 4 The goodness of model fit value corresponding to different order

M-ARMAX温度估计模型训练效果如图 6所示, 而图 7 (a)为模型训练时输出估计值与实际值回归图, 横坐标为高炉实际温度输出数据, 纵坐标为模型输出估计值. 图 7 (a)上的点基本分布在过原点斜率为1的曲线上及其周围, 说明训练的估计值可以很好地跟随实际值变化, 模型精度较高. 图 7 (b)图 7 (c)分别为建模误差自相关函数(Auto correlation function, ACF)曲线和概率密度函数(Probability density function, PDF)曲线.从建模误差自相关函数曲线分布可以看出, 当$t=0 $时, ACF = 1, 其他时刻自相关函数值在零线上下波动, 从建模误差PDF分布曲线可以看出, 残差均值基本为0, 曲线只有一个尖峰且光滑, 符合高斯曲线特性, 说明测试残差符合白噪声性质.综上, 采用M-ARMAX方法建立的十字测温中心点温度模型具有较好的建模精度, 对原建模数据可以已很高的精度进行拟合.

图 6 所提方法十字测温中心温度建模效果 Figure 6 Modeling results of center temperature estimation model for cross temperature measuring
图 7 建模散点图(左)及建模误差自相关函数(中)和PDF曲线(右) Figure 7 Scatter diagram of modeling and autocorrelation function and PDF curve of modeling error

采用新的工业数据对模型进行测试, 为此将本文基于M-ARMAX的建模方法与常见M-SVM模型、BP神经网络模型、PLSR模型以及PLSR的改进方法即动态偏最小二乘(Dynamic partial least squares, DPLS)方法进行多方面对比, 相关模型参数选取如下:

1) PLSR和DPLS主元个数根据累积方差百分比(Cumulative percent variance, CPV)大于0.8的原则, 输出变量主元个数为1和输入变量主元个数为3.其中DPLS输入输出得分主元之间滞后阶次根据交叉验证分别选取为3和2.

2) BP-NN的隐层节点和输出层节点的传输函数分别采用双曲正切S形函数tansig和线性的purelin函数.根据输入个数7, 隐层节点数设置为19.学习率初值设为0.05, 网络训练采用梯度自适应学习率训练算法traingdx, 即在训练过程中通过检查权重的修正值是否真正降低了误差函数来自动调整学习率.

3) M-SVM激活函数选为Sigmod函数.另外, M-SVM中的惩罚因子$C$表示对误差的容忍度, 而核函数参数$\sigma$表示所选的支持向量的影响范围.这里采用交叉验证法将其分别确定为$C=100$$\sigma=2$.

图 8为不同模型对十字测温装置中心温度点的实际估计效果图.可以看出, PLSR线性模型拟合值不能很好地跟随实际值得变化趋势, 模型估计精度最差.而M-SVM和BP模型得到的输出估计值基本上可以跟随实际值的变化趋势, 但是模型精度不高, 温度估计具有较大的误差.相比而言, 本文所提基于M-ARMAX的温度模型估计精度最好, 可以准确地跟随实际值的趋势变化并且模型估计效果最好.

图 8 不同建模方法下的十字测温中心温度估计效果对比 Figure 8 Estimation results of center temperature estimation model for cross temperature measuring with different method

图 9为不同方法对不同测温点的估计效果散点图, 图 10为相应估计误差的PDF分布曲线图.可以看出, 相对于其他几种方法, 本文所提M-ARMAX方法的温度估计值与实际值基本在对角线附件波动, 其估计误差PDF形状的中心更接近中轴, PDF形状更接近零均值高斯分布的白噪声.图 9图 10从不同方面验证了所提M-ARMAX温度估计模型的有效性和优越性.同时, 图 8 $\sim$ 10也说明所提基于M-ARMAX的温度估计模型能够解释输出变量的绝大部分信息.

图 9 不同建模方法温度估计值与实际值散点图 Figure 9 Scatter diagram of estimated temperature and actual temperature by different models
图 10 不同建模方法温度估计误差PDF曲线 Figure 10 PDF curve of temperature estimation error by different models

从算法结构上分析, SVM通过Sigmiod核函数将输出数据映射为高维Hilbert空间里的数据[21-23], 而BP-NN通过隐含层的激活函数映射输入数据, 两者都是在完成输入数据的映射或变换之后建立与输出数据的回归关系[24-25].虽然这两种算法都表现出较强的非线性逼近能力, 但应用在本文十字测温中心点温度建模和估计时效果不佳, 泛化性能较差.为此, 通过多次试验和比较分析, 我们转向于采用线性建模方法建立十字测温温度估计模型.本文所提M-ARMAX方法以及常规PLSR及其改进的DPLS都是基于线性回归的数据驱动建模方法. PLSR模型是综合了典型相关性分析、主成分分析以及回归分析的综合性数据驱动建模方法, 但此种方法在处理输入数据和输出数据的同时会丢失部分原始数据信息.此外, 常规PLSR模型关注的是当前时刻变量之间的联系, 不能表示包含历史信息的动态特性[26-27].而DPLS模型是在PLSR的基础上,通过引入过程变量的时滞信息来描述过程的动态特性, DPLS一方面解决了PLSR不能描述过程动态特征的问题, 另一方面继承了PLSR在通过提取输入输出数据特征向量时丢失信息的缺点.本文M-ARMAX时间序列模型能够很好地保留输入输出变量的全部方差变化信息, 通过自回归与滑动平均过程来描述变量的动态特性, 并运用离散时间序列分析数据的变化特点.因此, 采用M-ARMAX算法能够很好地说明十字测温温度系统数据之间的动态依存关系, 拟合精度高、建模效果好、算法简单且易于工程实现和工业应用.

3 结论

本文旨在建立一种模型结构简单、精度较高且易于工程实现的高炉十字测温中心温度在线估计方法及系统, 主要工作包括: 1)考虑到传统单一Pearson相关分析和单一因子分析在多输入多输出建模输入变量选择的不足, 提出基于因子分析与Pearson相关分析相结合的稳健特征选择新方法, 即首先采用因子分析从多个输出变量中找出一个主因子, 并与众多过程输入变量做Pearson相关分析的Filter特征选择, 选出与主因子相关性较强的过程变量作为最终建模输入变量; 2)对传统基于AIC模型阶数确定方法易造成模型阶次高、结构复杂的问题, 提出在AIC准则基础上进一步引入模型拟合优度目标函数来选取模型最小阶的方法, 可保证模型估计精度的同时降低模型阶次; 3)考虑到所研究问题的时序动态特性, 也为了易于工程实现, 采用M-ARMAX时间序列建模技术建立十字测温中心温度模型, 另外采用收敛速度较快的RLS技术用于辨识模型参数, 并基于残差分析方法检验建模效果.工业试验及比较分析表明, 相比于M-SVM和BP-NN方法, 所提M-ARMAX模型简单有效, 并且估计精度更高; 相比于PLS和DPLS回归方法, M-ARMAX模型能够准确地很好地说明十字测温温度系统数据之间的动态依存关系, 同时模型的拟合精度得到很大提高.另外, 统计分析以及工业应用表明, 所提方法能够满足实际生产对十字测温中心温度的测量精度和测量稳定性的生产操作要求.

参考文献
1
Song He-Da, Zhou Ping, Wang Hong, Chai Tian-You. Nonlinear subspace modeling of multivariate molten iron quality in blast furnace ironmaking and its application. Acta Automatica Sinica, 2015, 42(21): 1664-1679.
( 宋贺达, 周平, 王宏, 柴天佑. 高炉炼铁过程多元铁水质量非线性子空间建模及应用. 自动化学报, 2015, 42(11): 1664-1679.)
2
Zhou P, Yuan M, Wang H, Wang Z, Chai T Y. Multivariable dynamic modeling for molten iron quality using online sequential random vector functional-link networks with self-feedback connections. Information Sciences, 2015, 325: 237-255. DOI:10.1016/j.ins.2015.07.002
3
Jian L, Gao C H, Xia Z H. Constructing multiple kernel learning framework for blast furnace automation. IEEE Transactions on Automation Science and Engineering, 2012, 9(4): 763-777. DOI:10.1109/TASE.2012.2211100
4
Birk M, Marklund O, Medvedev A. Video monitoring of pulverized coal injection in the blast furnace. IEEE Transactions on Industry Applications, 2002, 38(2): 571-576. DOI:10.1109/28.993181
5
Wang Xian. Discussion on optimum design scheme of cross temperature measuring and cooling device for blast furnace. Iron and Steel Technology, 2010(1): 1-3.
( 王贤. 高炉十字测温冷却装置优化设计方案的探讨. 钢铁技术, 2010(1): 1-3.)
6
Tu Zheng-Huan. Research on quick replacement of cross temperature measuring gun for blast furnace. Machine China, 2015(22): 185-186.
( 涂正环. 高炉十字测温枪快速更换的研究. 中国机械, 2015(22): 185-186.)
7
Tang Zhen-Hao, Tang Li-Xin, Yang Yang. Blast furnace cross temperature prediction based on data-driven and intelligent optimization. Information and Control, 2014, 43(3): 355-360.
( 唐振浩, 唐立新, 杨阳. 基于数据驱动和智能优化的高炉十字测温温度预报. 信息与控制, 2014, 43(3): 355-360.)
8
Xu Hua-Bing, Xu Hua-Yan. Blast fuenace throat cross temperature forecast by time-sequence neural network based on TD algorithm. Industrial Furnace, 2013, 35(3): 44-48.
( 徐化冰, 徐华岩. 基于TD算法的时序神经网络高炉炉喉十字测温预报. 工业炉, 2013, 35(3): 44-48.)
9
Ueda S, Natsui S, Nogami H, Yagi J I, Ariyama T. Recent progress and future perspective on mathematical modeling of blast furnace. ISIJ International, 2010, 50(7): 914-923. DOI:10.2355/isijinternational.50.914
10
Nurkkala A, Pettersson F, Saxén H. Blast furnace dynamics using multiple autoregressive models with exogenous inputs. ISIJ International, 2012, 52(10): 1764-1771. DOI:10.2355/isijinternational.52.1764
11
Cai M L, Cai F, Shi A G, Zhou B, Zhang Y S. Chaotic time series prediction based on local-region multi-steps forecasting model. Advances in Neural Networks-ISNN 2004. Lecture Notes in Computer Science. Berlin, Heidelberg, Germany:Springer, 2004. 418-423
12
Xu Mei-Ling, Han Min. Factor echo state network for multivariate chaotic time series prediction. Acta Automatica Sinica, 2015, 41(5): 1042-1046.
( 许美玲, 韩敏. 多元混沌时间序列的因子回声状态网络预测模型. 自动化学报, 2015, 41(5): 1042-1046.)
13
Burnham K P, Anderson D R. Multimodel inference:understanding AIC and BIC in model selection. Sociological Methods and Research, 2004, 33(2): 261-304. DOI:10.1177/0049124104268644
14
Chaurasia A, Harel O. Using AIC in multiple linear regression framework with multiply imputed data. Health Services and Outcomes Research Methodology, 2012, 12(2-3): 219-233. DOI:10.1007/s10742-012-0088-8
15
Resende P A A, Dorea C C Y. Model identification using the efficient determination criterion. Journal of Multivariate Analysis, 2016, 150: 229-244. DOI:10.1016/j.jmva.2016.06.002
16
Xin Bin, Bai Yong-Qiang, Chen Jie. Two-stage ARMAX parameter identification based on bias-eliminated least squares estimation and Durbin's method. Acta Automatica Sinica, 2012, 38(3): 491-496.
( 辛斌, 白永强, 陈杰. 基于偏差消除最小二乘估计和Durbin方法的两阶段ARMAX参数辨识. 自动化学报, 2012, 38(3): 491-496.)
17
Wang C, Tang T. Recursive least squares estimation algorithm applied to a class of linear-in-parameters output error moving average systems. Applied Mathematics Letters, 2014, 29: 36-41.
18
Du Da-Jun, Shang Li-Li, Qi Bo, Fei Min-Rui. Convergence analysis of an online recursive identification method with uncomplete communication constraints. Acta Automatica Sinica, 2015, 41(8): 1502-1515.
( 杜大军, 商立立, 漆波, 费敏锐. 一种不完全信息下递推辨识方法及收敛性分析. 自动化学报, 2015, 41(8): 1502-1515.)
19
Karimi M. On the residual variance and the prediction error for the LSF estimation method and new modified finite sample criteria for autoregressive model order selection. IEEE Transactions on Signal Processing, 2005, 53(7): 2432-2441. DOI:10.1109/TSP.2005.849182
20
Frandsen S. Relevant criteria for testing the quality of models for turbulent wind speed fluctuations. Journal of Solar Energy Engineering Transactions of the Asme, 2008, 130(3): 1047-1057.
21
Gao Yan, Zhou Cheng-Hu, Su Fen-Zhen. Study on SVM classifications with multi-features of OLI images. Engineering of Surveying and Mapping, 2014, 23(6): 1-5, 10.
( 高燕, 周成虎, 苏奋振. 基于OLI影像多参数设置的SVM分类研究. 测绘工程, 2014, 23(6): 1-5, 10.)
22
Ahmad A S, Hassan M Y, Abdullah M P, Rahman H A, Hussin F, Abdullah H, Saidur R. A review on applications of ANN and SVM for building electrical energy consumption forecasting. Renewable and Sustainable Energy Reviews, 2014, 33: 102-109. DOI:10.1016/j.rser.2014.01.069
23
Kuo B C, Ho H H, Li C H, Hung C C, Taur J S. A kernel-based feature selection method for SVM with RBF kernel for hyperspectral image classification. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(1): 317-326. DOI:10.1109/JSTARS.2013.2262926
24
Ren T, Liu S, Yan G C, Wu H P. Temperature prediction of the molten salt collector tube using BP neural network. IET Renewable Power Generation, 2016, 10(2): 212-220. DOI:10.1049/iet-rpg.2015.0065
25
Zhang Y X, Gao X D, Katayama S. Weld appearance prediction with BP neural network improved by genetic algorithm during disk laser welding. Journal of Manufacturing Systems, 2015, 34: 53-59. DOI:10.1016/j.jmsy.2014.10.005
26
Farifteh J, van der Meer F, Atzberger C, Carranza E J M. Quantitative analysis of salt-affected soil reflectance spectra:a comparison of two adaptive methods (PLSR and ANN). Remote Sensing of Environment, 2007, 110(1): 59-78. DOI:10.1016/j.rse.2007.02.005
27
Sarkar A, Karki V, Aggarwal S K, Maurya G S, Kumar R, Rai A K, Mao X L, Russo R E. Evaluation of the prediction precision capability of partial least squares regression approach for analysis of high alloy steel by laser induced breakdown spectroscopy. Spectrochimica Acta Part B:Atomic Spectroscopy, 2015, 108: 8-14. DOI:10.1016/j.sab.2015.04.002