基于改进神经网络模型控制高炉炼铁过程-数学建模论文

基于改进的神经网络模型的高炉炼铁预测控制

摘要

高炉炼铁是钢铁工业的重要组成成分，高炉炉温控制是实际生产中的重要程序，建立可以进行炉温控制的炉温预报模型对实际生产具有重要意义。本文用铁水含硅量代表高炉炉温，通过建立多个模型并优化，对高炉炼铁铁水含硅量进行了动态预测。

针对问题一，要求建立一步和二步[Si]预测动态数学模型，首先对附件给出的数据进行预处理，修正了异常值。通过查阅文献得知喷煤量PML和鼓风量FL与铁水含硅量[Si]之间的关系具有滞后性，因而首先建立灰色关联度模型，得出喷煤量PML和鼓风量FL与铁水含硅量[Si]之间的关系皆相差6个炉次即12个小时。然后建立RBF神经网络模型，选取连续的m个样本学习对后一个样本进行一步预测；在二步预测过程中，以步长为二选取m 个样本对之后第二个样本进行预测。

针对问题二，要求验证问题一建立的模型的预测成功率，自主选取数据编写MATLAB 程序对问题一建立的RBF神经网络模型进行求解。分别得到各40组的一步预测和两步预测预测值，将预测值与真实值进行比较计算，得到一步预测和二步预测模型的预测数值成功率分别为84.36%和83.04%。再将铁水含硅量实际升降方向与预测方向比较计算得到一步预测和二步预测的预测方向成功率分别为81.58%和73.68%。本次建立的RBF神经网络模型验证结果较为良好，可用于高炉炼铁铁水硅含量的动态预测。

针对问题三，要求建立质量指标[S]的优化数学模型，并且讨论优化后的[Si]预测控制的预期效果。在RBF模型的基础上，建立粒子群模型对质量指标参数[S]进行优化，从而得到满足期望[S]参数。基于优化后的数据，选取样本编写MATLAB程序对RBF神经网络模型进行求解，将得到的数据与真实值进行比较计算，得到优化后的预测数值成功率达到99.04%，效果较好。证明经过优化后的质量指标对于[Si]的准确预测控制更加准确。

本文建立多个模型并对预测模型进行优化，得到了合理且准确率高的铁水硅含量预测模型，并且对模型的优缺点进行了合理的评价，对控制高炉炼铁炉温操作具有十分重要的参考价值。

关键词：滞后性；动态预测；RBF神经网络；粒子群；MATLAB

1. 问题重述

炼铁过程生产指标产量、能耗、铁水质量等指标都与冶炼过程的一项控制性中间指标——炉温，即铁水含硅量[Si]（铁水含硅质量百分数）密切相关。对2小时后或4小时后高炉炉温上升或下降的预测，即[Si]时间序列的预测关系着当前高炉各项操作参数的调控方向。因此，[Si]的准确预测控制建模成为冶炼过程优化与预测控制的关键技术。本项目仅提供由铁水含硅量[Si]、含硫量[S]、喷煤量PML和鼓风量FL组成的数据库作为数学建模分析和数据挖掘的基础。试求解以下问题：

（1）从给定数据表中[Si]-[S]-FL-PML依序号排列的1000炉生产大数据中，自主选取学习样本和算法，建立[Si]预测动态数学模型，包括一步预测模型和二步预测模型。

（2）自主选取验证样本，验证你所建立的数学模型的预测成功率，并且讨论其动态预测控制的可行性。

（3）以质量指标铁水含硫量[S]为例，含硫量低，铁水质量好，可以生产优质钢，制造优质装备。试建立质量指标[S]的优化数学模型，并且讨论按照优化模型计算结果进行[Si]预测控制的预期效果。

2. 模型假设

1)假设除题中给出的影响因素的其他影响因素影响忽略不计。

2)假设给出数据的高炉运行情况良好，无异常运行。

3. 通用符号说明

4.铁水含硅量动态预测模型的建立

4.1问题分析

本文要求解答预测[Si]动态数学模型，在查阅有关炼铁的文献时发现存在着大滞后的现象，即喷煤量PML和鼓风量FL与铁水含硅量[Si]之间的关系具有滞后性，且滞后性不尽相同。因而建立灰色关联度模型，求解出喷煤量PML和鼓风量FL与铁水含硅量[Si]之间分别的滞后时间。拟建立RBF神经网络模型，选取连续的m个样本学习对后一个样本进行一步预测；在二步预测过程中，打算以步长为二选取m个样本对之后第二个样本进行预测。

4.2模型准备

异常值的存在直接影响到建模准确率和预测精度，因此对异常值进行合理、有效的修正。利用Excel画出原始样本铁水含硅量、铁水含硫量、风量和喷煤量的时间序列曲线图，如图1、2、3、4。

图1 原始样本铁水含硅量的时间序列曲线

图2 原始样本铁水含硫量的时间序列曲线

图3 原始样本集中喷煤量的时间序列曲线

图4 原始样本风量的时间序列曲线

由图1、2、3、4可得每一时刻点的铁水含硅量、铁水含硫量、风量和喷煤量，对波动较大的数据，予以剔除处理。

采用依拉达准则（3σ准则）处理异常值，对本组数据样本{}121000,,...,X x x x =，如果存在偏差大于3σ的数值，则认为它是异常数值。其中

σ=

= （4-1）

其中，x 为数据样本的平均值。

编写MATLAB 程序经3σ准则检测得出：初始样本集1000组数据中，一共有64组异常数据，异常率6.4%。异常数据见表1。程序见附录1。

表1 异常数据

都有缺失值。选用“高次线性插值法”填补缺失值。填补的数据值见表2。

表2 异常数据替换数值

表24.3模型建立 4.3.1 灰色关联度模型 （1）输入变量的选择

引入灰色关联度计算高炉各输入变量对硅含量的影响程度。灰色关联度能够计算因素之间的相互关系，并将对目标影响程度表示为量化值。灰色关联度越高，表示变量对目标变量的相关性越大。灰色关联方法步骤如下：

Step1：整理得到灰色关联序列()()()()00001,2,...,X x x x n =和对比序列

()()()()1,2,...,i i i i X x x x n =，1,2i =，其中0X 为硅含量序列， 12,x x 分别为喷煤，风量的序列。

Step2：计算灰色关联系数。

()()()min max

0max

i i k y k y k ρξρ?+?=

-+? （4-2）

式中

()i k ξ——关联系数

ρ——差分系数，满足[1,0)ρ∈

min ?——()()()min 0min min i t

k

y k y k ?=-

max ?——()()()max

max max i

t

k

y k y k ?=-

()i k ξ的值即为参考序列的每一个值和对比序列对应值的关联程度。

Step3：通过式（4-2）计算得到灰色关联系数之后，关联系数的平均值常被用来作为该参考序列的灰色关联度。

()()01

1,n

i i k x x k n γξ==∑ （4-3）

式中()0,i x x γ——灰色关联度

Step4：设置灰色关联度阈值，选取输入变量。 （2）输入变量滞后时间确定

高炉运行是典型的大时滞过程，输入变量对硅含量的影响存在着较大的时延。考虑到精度和实际生产的需求，时延问题不能被忽略。采用相关系数分析方法计算不同时延的输入变量对硅含量的影响程度，完成之后时间的确定。相关系数定义如下：

()()

n

ij

i

j

i x

x y

y

R --=

∑ （4-4）

式中

i R ——代表第i 个输入的相关系数

i x ——输入变量的平均值 y ——输出变量的平均值

时间序列设置为（0,1,2,3,4,5,6,7,8,9），通过比较相关系数的大小得出喷煤、风量的滞后时间。

4.3.2 RBF 神经网络预测模型

RBF 神经网络有很强的逼近能力、分类能力和学习速度。其工作原理是把网络看成对未知函数的逼近，任何函数都可以表示成一组基函数的加权和，也即选择各隐层神经元的传输函数，使之构成一组基函数来逼近未知函数。RBF 神经网络由一个输入层、一个隐含层和一个输出层组成。

图5 RBF 神经网络示意图

图5即为神经网络示意图。

设输入层的输入为12[,,...,]n X x x x =，实际输出为12[,,...,]n Y y y y =。输入层实现从

()i X R X →的非线性映射，输出层实现从()i k R X y →的线性映射，输出层第k 个神经网络输出为

()1,1,...,m

k ik i i y w R X k p ===∑ （4-5）

式中

n 为输入层节点数； m 为隐含层节点数； p 为输出层节点数；

ik w 为隐含层第i 个神经元与输出层第k 个神经元的连接权值；

()i R X 为隐含层第i 个神经元的作用函数，即

()(

)

2

2exp /2,1,...,i i

i R X X C i m σ=--= （4-6）

式中

X 为n 维输入向量；

i C 为第i 个基函数的中心，与X 具有相同维数的向量；

i σ为第i 个基函数的宽度；

m 为感知单元的个数（隐含层节点数）；

i X C -为向量i X C -的范数，它通常表示x 与i C 之间的距离；

()i R X 在i C 处有唯一的最大值，随着i X C -的增大，()i R X 迅速衰减到0。

当确定了RBF 神经网络的聚类中心i C 、权值ik w 、及i σ以后，就可以求出给定某一输入时，网络对应输出值。

①一步预测

当时间延时为t 时，对第i 个si 含量预测的输出值为

()()1,1,...,m

k i i k t i y w R X k p -===∑ （4-7）

②二步预测

对第（i+1）个si 含量预测的输出值为

()()11,1,...,m

k i i k t i y w R X k p -+===∑ （4-8）

4.3.3 预测成功率模型

Step1数值预测成功率和成功率相对误差（ε）的计算

100%c s s

x x x ε-=

? （4-9）

其中

s x 为高炉铁水[si]含量实际值； c x 为高炉铁水[si]含量预测值。 Step2炉温升降方向预测

()()11i i i i y x x x x --=--预测 （4-10）

i x 为第i 次硅的实测含量，i x 预测为第i 次硅的预测含量

若0y ≥，则炉温升降方向预测正确 若0y ≤，则炉温升降方向预测不正确 4.4模型求解

通过编写MATLAB 程序，求解灰色关联度模型，得到了不同时延的相关性系数，如表3。程序见附录2。

由表n 个[Si]具有最大关系的是第6n -个喷煤量和风量。

5. 铁水含硅量动态预测模型的求解

5.1问题分析

本题要求自主选取样本对问题一建立的模型进行验证，并求解预测成功率。首先一步

预测针对给出的1000个样本，拟采取任意选取连续50组数据预测一个[Si]的方法，得到40个预测值，对得到的预测值跟实际值进行比较，求得数值预测成功率。之后再将这40个结果的升降方向与实际值进行对比，得到炉温升降方向预测成功率。二步预测以步长为二选取50组数据预测 [Si]，将预测数值的比较变为与前两个数值的比较。 5.2模型建立

采用4.3中建立的模型，在此不多做赘述。 5.3模型求解

Step1神经网络一步预测模型

通过编写MATLAB 程序，得到预测的40个数据，将之与真实值进行对比，画出图1。程序见附录3。

图6 一步预测硅含量预测值与真实值对比

由图6得硅含量预测值与真实值变化趋势相近。

再将第n 个实际值减去第1n -个实际值的结果与第n 个预测值减去第1n -个实际值的结果相乘，若结果为正，则方向预测正确；若结果为负，则方向预测存在偏差。得到的结果如表4。

根据表

（1）由于高炉炼铁存在相当大的影响因素，而本题只考虑两项，因而数据也存在较大的误差。所以将相对误差的最大值0.5043和最小值0.0017去除，计算其余各项的平均值，为0.1564，即数值预测成功率为84.36%。

（2）由（1）去除数据后，得到预测方向一致的有31个，不一致的有7个，即炉温升降方向预测成功率为81.58%。

Step2神经网络二步预测模型

通过编写MATLAB 程序，得到预测的40个数据，将之与真实值进行对比，画出图2。程序见附录4。

图7 二步预测硅含量预测值与真实值对比

由图7得硅含量预测值与真实值变化趋势相近。

再将第n 个实际值减去第1n -个实际值的结果与第n 个预测值减去第1n -个实际值的结果相乘，若结果为正，则方向预测正确；若结果为负，则方向预测存在偏差。得到的结果如表5。

根据表5

（1）由于高炉炼铁存在相当大的影响因素，而本题只考虑两项，因而数据也存在较大的误差。所以将相对误差的最大值0.9003和最小值0.0031去除，计算其余各项的平均值，为0.1696，即数值预测成功率为83.04%。

（2）由（1）去除数据后，得到预测方向一致的有28个，不一致的有9个，即炉温升降方向预测成功率为73.68%。

Step3可行性讨论

本问题对问题一建立的模型进行了验证，发现一步预测和二步预测模型的数值预测成功率均高于80%，一步预测模型的炉温升降方向预测成功率也高于80%，成功率都比较高，这两个模型都具有较好的可行性。

6. 铁水含硅量动态预测模型的优化、求解

6.1问题分析

本问题要求建立质量指标[S]的优化数学模型，并且讨论按照优化模型计算结果进行[Si]预测控制的预期效果。在众多优化方法中由于粒子群算法更加方便简洁，因而拟采用粒子群算法对神经网络模型进行优化。之后选取连续的50个样本进行训练，对之后的40个数据进行预测，对预测的结果与真实值作比较，讨论预计效果。 6.2模型建立 6.2.1 粒子群算法

粒子群算法是一种基于信息共享机制的优化算法。在其算法中，每个粒子均被视为是D 维解空间的一个点。()12,,...,i i i iD X x x x =表示第i 个粒子，()12,,...,i i i iD P P P P =为其经过的最好位置。G 为群体中最好粒子的位置索引号。()12,,...,i i i iD V v v v =表示第i 个粒子的速度。依照（1）（2）来调整粒子的速度和位置。

()()12**()**()*id id id id gd gd V w V c rand P X c rand P X =+-+- （6-1）

id id id X X V =+ （6-2）

其中

1,2,...,;1,2,...,d D i N ==（N 为群体规模）

； 1c 、2c 为常数，称为学习因子；

()rand 是[0,1]上的随机数；

w 是惯性权重。

粒子速度与粒子位置的调整公式可以分为三个部分来解读；*id w V 表示惯性权重和粒子前一时刻速度的乘积，说明下一时刻时，该粒子的速度与前一时刻的关系，它体现粒子的记忆性，决定了算法的全局寻优和局部寻优的能力；()()*id id rand P X -表示粒子本身的历史属性，说明了粒子的历史最优位置对当前时刻速度调整的影响力，表明了算法具有全局寻优能力；()()*gd gd rand P X -是群体影响部分，说明子的群体最优位置对当前速度调整的影响程度。因此，这三个部分必须给予恰当、合理的设置，才能保证算法的性能最大化。

图8 粒子群算法过程图

6.2.2 高炉冶炼过程参数寻优

Step1操作参数寻优的约束条件和目标函数

各操作变量的期望值范围均由高炉专家经验确定。[]Si 含量、喷煤量PML 和风量FL 的最佳期望值分别为0.57，13.79，1799.53。综上所述，通过粒子群算法求解所示的优化目标函数，也即求解（6-3）的最小化问题，最终实现对操作参数的寻优。

2

2

2

min 0.310.310.410.5713.791799.530.30.7..101616001750z

Si PML FL z Si s t PML FL ?????

?=-+-+- ? ? ???????

≤≤??≤≤??≤≤?

（6-3）

其中[]Si ，PML ，FL 分别为[]Si 含量、风量和喷煤量的实际值。 Step2粒子群算法参数测定 （1）惯性权重w

惯性权重代表着粒子前一时刻速度对当前时刻速度的影响程度。若w 相对比较大则粒子群算法在寻优过程中的局部搜索能力强；若0w =，则此时粒子群算法仅仅只能完成局部寻优。经过考虑我们将w 设置为0.7 （2）学习因子c

学习因子1c ，2c c2分别表示飞向个体局部最优位置best P 和全局最优值best g 得加速权重。若c 较小，则粒子会在目标E 域附近浮动；若c 较大，则粒子可能直接飞向目标区域。我们通过实验认为122c c ==时为最好的选择。 （3）最大速度max V

max V 的设置不能超过粒子的搜索范围，max max 0.8V X =。 群体规模取500个。 6.3模型求解

编写MATLAB 程序，对粒子群优化后的模型进行求解。分别画出未优化模型与优化后模型的预测数据与实际值的对比图，得到图8。数据对比见表4。

图8 未优化模型与优化后模型的预测数据与实际值的对比

由图8可得，优化后的模型预测准确率较好，基本与原始数据相同。

根据表699.04%。

对以上结果分析可得，由粒子群模型改进后的预测模型准确率高达99.04%，证明经粒子群模型优化后的质量指标对于[Si]的预测控制的预期效果更好。

7. 模型评价

优点：

（1）问题一建立灰色关联度模型，该模型简单方便，符合题目要求并求解出铁水含硅量相关系数的滞后性，较好的分析了铁水含硅量与相关因素的具体影响效果，结果较好。（2）问题一和二建立神经网络模型并求解，该模型对一步预测和二步预测进行了分析，并且结果良好，神经网络模型适用于高炉炼铁铁水硅含量预测。

（3）本题第三问利用粒子群算法对神经网络模型进行了优化，方法简单，数据分析结果较好。

缺点：

在建立模型时由于模型的适用性，对数据进行了简化，导致结果的不准确性。

参考文献

[1] 吴金花, 高炉冶炼过程分析及其铁水硅含量预测模型研究, 燕山大学，2016年5月。

[2]胡登飞，基于神经网络预测模型的高炉喷煤操作模式优化，内蒙古科技大学，2014年6

月。

[3]曲红玉，沈婷，曾璐，基于灰色神经网络的高炉冶炼优质铁水预测，哈尔滨商业大学金

融学院。

[4]石琳,李志玲,崔桂梅,基于偏最小二乘回归的高炉铁水硅含量模型,内蒙古大学学报：自然

科学版,2010,41（4）：427-430。

[5]李启会，高炉冶炼过程的模糊辨识、预测及控制，浙江大学，2005年l月。

[6]赵敏，高炉冶炼过程的复杂性机理及其预测研究，浙江大学，2008年6月。

附录

附录1 寻找异常数据MATLAB程序

A =[12.14 11.62 10.32 12.05 12.92 12.01 14.43 13.85 14.48 12.58 12

13.81 11.59 8.02 10.21 10.59 11.63 11.92 13.86 17.18 16.67 14.33

13.57 15.51 13.03 13.46 13.71 13.13 15.72 13.39 12.18 13.85 12.93

12.86 15.24 14.13 13.57 14 13.78 15.05 14.87 15.78 13.82 13.28

15.97 15 12.66 6.57 15.03 13.44 13.8 11.66 11.61 12.68 15.38

15.04 12.79 13.94 15.18 16.1 15.18 12 11.71 14.87 13.74 13.55

12.66 14.45 15.09 15.76 13.45 15.11 13.95 12.61 13.73 15.26 13.64

14.88 15.92 14.49 13.98 16.24 13.86 14.1 7.89 13.08 ];

[mu,sigma,muci,sigmaci]=normfit(A,0.05)

附录2 灰色关联度MATLAB程序

clc,clear

load x.txt %把原始数据存放在纯文本文件x.txt 中，其中把数据的"替换替换成. for i=1:3

x(i,:)=x(i,:)/x(i,1); %标准化数据

end

data=x;

n=size(data,2); %求矩阵的列数，即观测时刻的个数

ck=data(1,:); %提出参考数列

bj=data(2:end,:); %提出比较数列

m2=size(bj,1); %求比较数列的个数

for j=1:m2

t(j,:)=bj(j,:)-ck;

end

mn=min(min(abs(t'))); %求最小差

mx=max(max(abs(t'))); %求最大差

rho=0.5; %分辨系数设置

ksi=(mn+rho*mx)./(abs(t)+rho*mx); %求关联系数

r=sum(ksi')/n %求关联度

[rs,rind]=sort(r,'descend') %对关联度进行排序

附录3 神经网络一步预测模型MATLAB程序

clc,clear

a=load('si.txt');

a=a'

for i=1:40

P=a([3*i-2:3*i-1],[1:end-1]);[PN,PS1]=mapminmax(P);

T=a(3*i,[1:end-1]);[TN,PS2]=mapminmax(T);

net1=newrb(PN,TN);

x=a([3*i-2:3*i-1],end);xn=mapminmax('apply',x,PS1);

yn1(i)=sim(net1,xn);y1(i)=mapminmax('reverse',yn1(i),PS2);

deltal=abs(a(3*i,50)-y1(i))/a(3*i,50);

net2=feedforwardnet(4);

net2=train(net2,PN,TN);

yn2(i)=net2(xn);y2(i)=mapminmax('reverse',yn2(i),PS2)

epsilon(i)=a(3*i,50)-y2(i) %计算残差

delta(i)=abs(epsilon(i)./a(3*i,50)) %计算相对误差

sj(i)=a(3*i,50);

end

y=1:40;

figure(1);%分别做新数据和旧数据与实际数据的对比值

subplot(1,1,1);

plot(y,y2,'r-o',y,sj,'b--+');

legend('网络输出','实际');

xlabel('炉数');ylabel('含硅量');

title('运用功能工具箱含硅量学习和测试对比图')

title('运用功能工具箱含硅量学习和测试对比图')

附录4 神经网络二步预测模型MATLAB程序

clc,clear

a=load('si2.txt');

a=a'

for i=1:40

P=a([3*i-2:3*i-1],[1:end-1]);[PN,PS1]=mapminmax(P);

T=a(3*i,[1:end-1]);[TN,PS2]=mapminmax(T);

net1=newrb(PN,TN);

x=a([3*i-2:3*i-1],end);xn=mapminmax('apply',x,PS1);

yn1(i)=sim(net1,xn);y1(i)=mapminmax('reverse',yn1(i),PS2);

deltal=abs(a(3*i,25)-y1(i))/a(3*i,25);

net2=feedforwardnet(4);

net2=train(net2,PN,TN);

yn2(i)=net2(xn);y2(i)=mapminmax('reverse',yn2(i),PS2)

epsilon(i)=a(3*i,25)-y2(i) %计算残差

delta(i)=abs(epsilon(i)./a(3*i,25)) %计算相对误差

sj(i)=a(3*i,25);

end

y=1:40;

figure(1);%分别做新数据和旧数据与实际数据的对比值

subplot(1,1,1);

plot(y,y2,'r-o',y,sj,'b--+');

legend('网络输出','实际');

xlabel('炉数');ylabel('含硅量');

title('运用功能工具箱含硅量学习和测试对比图')

title('运用功能工具箱含硅量学习和测试对比图')

附录5 粒子群优化模型程序

function result=fitness(x,D)

sum=0;

z=[0.57 13.79 1799.53 0.006];

for i=1:D

Z=z(i);

sum=sum+(1-x(i)/Z)^2;

end

result=sum;

clear all;

clc;

format long;

%------给定初始化条件----------------------------------------------

c1=2; %学习因子1

c2=2; %学习因子2

w=0.7; %惯性权重

MaxDT=1000; %最大迭代次数

D=3; %搜索空间维数（未知数个数）

M=500; %初始化群体个体数目

eps=10^(-6); %设置精度(在已知最小值时候用)

%------初始化种群的个体(可以在这里限定位置和速度的范围)------------ for i=1:M

for j=1:D

x(i,j)=randn; %随机初始化位置

v(i,j)=randn; %随机初始化速度

end

end

%------先计算各个粒子的适应度，并初始化p(i)和gbest--------------------

for i=1:M

p(i)=fitness(x(i,:),D);

y(i,:)=x(i,:);

end

gbest=x(1,:); %gbest为全局最优

for i=2:M

if fitness(x(i,:),D) < fitness(gbest,D)

gbest=x(i,:);

end

end

%------进入主要循环，按照公式依次迭代，直到满足精度要求------------

for t=1:MaxDT

for i=1:M

v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(gbest-x(i,:));

x(i,:)=x(i,:)+v(i,:);

if fitness(x(i,:),D)

p(i)=fitness(x(i,:),D);

y(i,:)=x(i,:);

end

if p(i)

gbest=y(i,:);

end

end

end

%------显示计算结果

disp('*************************************************************') disp('函数的全局最优位置为：')

Solution=gbest'

disp('最后得到的优化极值为：')

Result=fitness(gbest,D)

disp('*************************************************************') 附录6 优化预测实际对比程序

clear,clc

p=[2.62E-04 1.37E+01 1.80E+03

6.21E-04 1.41E+01 1.82E+03

1.29E-05 1.38E+01 1.80E+03

3.42E-05 1.39E+01 1.81E+03

2.95E-06 1.38E+01 1.80E+03

7.14E-09 1.38E+01 1.80E+03

1.32E-06 1.38E+01 1.80E+03

4.43E-10 1.38E+01 1.80E+03