1、.基于遗传算法的温度预报模型参数优化1 问题描述.近年来,随着纯净钢生产技术的进步和连铸技术的发展,炉外精炼工艺与设备迅速普及。其中,LF以其优异的综合性能,在实际生产中得到了广泛应用。而点测钢水温度的高成本,使精炼炉温度预报成为了极具实际意义的研究。因此,钢水温度预报模型的建立显得至关重要。目前,温度预报模型的建立主要采用 3 种方法,即:机理模型、 “黑箱模型”和“灰箱模型” 。机理模型是指用尽可能准确的数学方程来描述过程机理而建立的模型,而“黑箱模型”则采用一些数学方法(智能算法、回归算法等) ,结合实验数据或实际生产中的数据,生成一种只有输入和输出的模型,而完全不考虑过程机理。然而,机
2、理模型需要考虑的因素太多,且这些因素具有相当的不确定性,因此模型中的许多参数很难得到,严重影响了温度预测的精度;“黑箱模型”则完全建立在数据的基础上,如果生产环境和条件改变导致了数据改变,原先模型的准确性就得不到保障,即模型的可移植性差。因此,本文主要研究“灰箱模型”的建立,即采用机理建模与数据建模相结合的方式,应用智能优化算法对机理模型中较难获得的参数进行辨识与确定,其流程如下图所示:图 1 温度预报模型建立的方案流程图而参数的辨识过程,即是基于智能优化算法,根据输入输出条件,对机理模型不易获得的参数进行寻优的过程。2 理论基础2.1 机理模型.本研究所基于的机理模型,是根据 LF 炉精炼过
3、程中的传热基本方程、能量守恒方程和质量守恒方程等来建立的。并运用有限差分法和有限元法等,通过控制初始条件和边界条件来对模型进行求解,得到钢包内的温度情况。2.1.1 热平衡分析3.1.1 吸热与散热精炼过程中,钢水热量的来源与去向大致如图 2 所示。图 2 热量的来源与去向将钢水和炉渣作为一个系统,来推导吸热和散热与其温度变化的关系。由于系统在加热与散热的过程中,其内能变化都体现在温度的变化上,所以系统实际吸收的热量为通过电弧加热所吸收的热量与冶炼过程中散去的热量之差,即:()stslesachlscmTQ(1)式中, 、 为钢水和渣的比热容, 、 为钢水和渣的重量, 为stcl stlmT系
4、统总的温度变化, 、 、 、 分别代表通过电弧加热所吸收的热量,eQsalch和通过渣面及炉盖散热、炉壁及炉底散热、加渣料及加合金影响所损失的热量。为最终求出温度随时间变化的曲线,可对式(1)两端对时间取微分。由此,对每个温度的影响因素进行分析与建模,分别求解它们的热量变化速率,再将模型综合后求积分,即可得出温度随时间变化规律的曲线,温度的计算公式为:(2)()()esachlstslsdQdQTmcdttt2.1.2 机理模型描述.(1)电弧加热模型系统所吸收的热量应与总的耗电量相关,二者间是一个乘系数的关系,即:earclQ(3)其中 为该次冶炼过程中的总耗电量,而 是系统的能量吸收系数,
5、即alQarc电弧效率。(2)渣面及炉盖散热模型LF 精炼过程中,通过渣面损失的热量 包括通过渣面的辐射和对流的热saQ损失 ,还包括吹氩带走的热量 ,即:srQg(4)sarg渣面及炉盖散热的近似计算方法,是将渣面散热和炉盖散热视为一个整体,通过计算水冷炉盖中冷却水带走的热量,来确定渣面及炉盖的散热量,即有:(5)()srwaloutwinQcGtT式中, 、 分别为冷却水的比热容和密度, 为冷却水流量,其单位wc为 ; 为冶炼时间,单位为 ; 、 分别为冷却水的出口、入口3/mhalt hwoutin温度。而计算高温气体带走热量的经验计算方法为:对小于 120t 吨位的钢包炉,气体带走热量
6、约占有功功率的 5%-6%;而对于大于 120t 吨位的钢包炉,气体带走的热量大约占有功功率的 7%-9%。因此有:(6)gQP其中 为气体带走热量占有功功率的百分数, 为有功功率。将 和g srQ综合考虑,并考虑精炼时间和吹氩时间,以冷却水流量、冷却水出入口温度,gQ以及有功功率作为输入,以 为模型待确定参数,即可搭建该部分模型。g.(3)炉壁及炉底散热模型炉壁散热可当作一维非稳态导热问题建模,利用其初始条件和边界条件来求解。为简化问题,将炉壁的倾斜度忽略不计,将其视为圆柱体;并且将炉壁材料的导热系数和比热容取平均值,视为常量。计算时,可将炉壁视为无限长圆筒壁,在柱坐标下求解。就炉壁而言,其
7、导热微分方程式为:2()ccppTTtr(7)上式中, 、 、 分别为钢包材料的密度、比热容和导热系数; 为pcp cT炉壁径向上某点的温度, 为该点距钢包中轴线的距离(即半径) ; 为时间,r t单位是 。式(7)的初始条件( )为:s0t201ln(/)clstTr(8)式中, 、 分别为首次测温时钢水温度和钢包外壁的温度; 、 分0stTl 1r2别为钢包的内径和外径。式(8)的边界条件是,在 时有:1rcstT(9)在 有:2r2|()cpracThT(10)以上两式中, 、 为钢水和环境的实时温度; 为钢包侧壁的综合对stTa ac流换热系数,其单位为 。运用有限差分的方法,可对该微
8、分方程进2/()WmC行求解。同理,对于炉底散热问题,也可用同样的方法进行讨论,这里不再赘述。.(4)加渣料与合金的影响加料是LF 的精炼流程中的一个重要环节,不同渣料的熔化热不同,因此加渣料对于钢水温度的影响取决于不同渣料的热物性参数和加入量的多少;与加渣料类似,不同合金的熔化热、熔解热以及化学反应吸收或放出的热量都不同。根据大量的实际数据统计,以下给出一些常见渣料、合金的温降系数值,如表1所示。其中,加渣料、合金对钢水的平均温降系数 ,单位为 。ik/10Ckg表1 部分常见渣料与合金的温降系数需要注意的是,上表中的负值表示降温。于是,加渣料与合金的热效应可由式(11)计算。10chime
9、ltdQkt(11)其中 为加某种渣料或合金的重量( ) , 为熔化时间。imkgelt2.1.3 机理模型仿真依据各部分机理模型,以其所有输出与钢水、渣的重量及比热容作为总模型输入,最终通过仿真输出钢水温度。该模型的Simulink框图及仿真结果如下所示。.图3 总机理模型Simulink仿真框图图4 总机理模型Simulink仿真结果上图中,横坐标表示时间(s) ,纵坐标表示钢水温度( ) ,该图显示了一个精炼周期内(约45min)钢水温度的变化情况。显然,在电弧加热的时间段内,钢水温度显著上升;在其它时段,由于散热,钢水温度有缓慢下降的趋势。而且在一个精炼周期内,钢水大约升温70 ,基本
10、符合实际情况。然而,由于电C弧效率、钢包材料的导热系数与比热容等参数很难准确得知,机理模型中只能大致在一定范围内取值,因此机理模型仅能准确反映温度变化趋势,难以达到精确预测温度的效果。因此,本研究将利用遗传算法来对这些难以准确获得的参数进行辨识与优化。.2.2 遗传算法2.2.1 算法简介遗传算法(GA)是基于 Darwin 进化论和 Mendal 遗传学说的一种优化搜索方法,从 20 世纪 60 年代开始兴起。它是基于自然选择和基因遗传学原理的搜索方法,将“优胜劣汰,适者生存”的生物进化原理引入待优化参数形成的编码串群体中,按照一定的适配值函数及一系列遗传操作对各个个体进行筛选,从而使适配值
11、高的个体被保留下来,组成新的群体。新群体包含上一代的大量信息,并且引入了新的优于上一代的个体。这样周而复始,群体中各个体适应度不断提高,直至满足一定条件。此时,群体中适配值最高的个体即为待优化参数的最优解。遗传算法主要有以下特点:(1)遗传算法的处理对象不是参变量本身,而是参变量编码后的称为人工染色体的位串,使得遗传算法可直接对集合、队列、矩阵、图表等结构对象进行操作。这一点使得遗传算法的应用范围很广。(2)遗传算法是多点搜索,而不是单点,从而避免陷入局部最优,逐步逼近全局最优解。也正是其固有的并行性,是遗传算法优于其他算法的关键。(3)与自然界类似,遗传算法对求解问题的本身可以一无所知,如对
12、搜索空间没有特殊要求(如连续、凸性等) ,对目标函数也不要求诸如连续性、导数存在和单峰等假设,它所需要的仅是对算法所产生的每个染色体进行评价,并基于适应值排列等级来选择染色体,使适应值好的染色体比适应值差的染色体有更多的繁殖机会。遗传算法利用简单的编码技术和繁殖机制来表现复杂的现象,从而解决非常困难的问题。(4)遗传算法是一种自适应随机搜索方法。遗传算法是以概率原则指导搜索,适应值高的个体,在进化过程中将被赋予更高的选择概率,在下一代中有更大繁殖机会。遗传算法可以更有效地利用已有的信息来搜寻那些有希望改善解质量的串,比一般的随机搜索有更高的搜索效率。因此,利用遗传算法的上述优点与广泛适用性,对
13、机理模型中难以获得的.参数加以辨识后,再应用到机理模型中,是一种能进一步发展和完善机理模型的切实可靠的方法。2.2.2 编码方法编码是应用遗传算法时要解决的首要问题,也是设计遗传算法时的一个关键步骤。编码除了决定个体的染色体排列形式之外,还决定了个体从搜索空间的基因型变化到解空间的表现型时的解码方法,编码方法也影响到交叉算子、变异算子等遗传算子的运算方法。迄今为止,最常用的是二进制编码方法和浮点数编码方法。本例中采用浮点数编码方法,即将个体的每一个基因值用某一范围内的一个浮点数来表示,个体编码长度等于其决策变量的个数。因为这种编码方法使用的是决策变量的真实值,所以浮点数编码方法也称真值编码方法
14、或实参编码方法。浮点数编码方法大致具有以下几个优点: (1)适合于在遗传算法中表示范围较大的数。(2)适合于精度要求较高的遗传算法。(3)便于较大空间的遗传搜索。(4)改善了遗传算法的计算复杂性,提高了运算效率。(5)便于设计针对问题的专门知识的知识型遗传算子。根据之前建立起的机理模型,选择系统中不易获得的参数进行辨识,它们是:电弧效率 、钢包材料导热系数 、钢包材料比热容 。将这些参数作arcppc为遗传算法中组成种群中个体的基因,并以向量 来表示。采用浮(,)Tarc点数编码方法,即用这些参数的真值作为其编码,并且限定这些参数的范围。根据经验,给出电弧效率、钢包材料导热系数以及比热容的参考
15、范围,即:, , 0.5,8arc1.3,65p10,38pc(12)2.2.3 适应度函数的建立.个体适应度函数的值不仅代表了该个体本身的好坏,还决定了在下一代遗传中,该个体存留下来的概率。本研究利用模型计算的终点温度值与实际测量的终点温度值间的误差来构造适应度函数,即:21()injfT(13)式中, 表示某个体的适应度, 为输入输出数据的组数,而 为第if njT组数据的计算温度值与实际温度值的误差值。为方便计算,将图 3 所示的jSimulink 机理模型编写为一个 m 函数,供遗传算法所调用。该函数以总耗电量、加热时间、有功功率、钢包尺寸、加料量、钢水初始温度等信息为输入,用作为待优
16、化参数,以钢水终点温度为函数计算的输出,以钢水的(,)Tarcp终点测温值作为比较。而输入信息和钢水的终点测温值,均取自多炉次的现场操作记录与数据。2.2.4 基本操作(1)选择操作选择是从一个旧种群中选择生命力强的个体位串,产生新种群的过程。本例中选择操作采用适度比例法,它根据某个体的适应度与该代全部个体适应度之和的比值,来决定其子孙遗留的可能性,即在第 代中,某个体 被选取的概ni率为:isiNfP(14)其中, 为某个体的适应度, 为该代全部个体的适应度之和,可根据if f式(13)来进行计算。确定概率后,可用轮盘赌的选择方法来实现选择操作。.例如,用计算机生成 01 之间均匀分布的随机
17、数,若 ,则当产生的随0.5siP机数在 00.5 之间时,该串被复制,否则被淘汰。(2)交叉操作交叉运算是指两个相互配对的染色体按某种方式相互交换某部分基因,从而形成两个新的个体。在遗传算法中,将种群中的 个个体以随机的方式组成N对配对个体组,交叉操作在这些配对个体组中的两个个体间进行的。/2N考察某配对个体组中的两个个体 、 ,交叉操作采用一定方式将它们变ab为两个新的个体 、 。在遗传算法中,交叉操作过程需要满足:ab(15)基于上式,浮点数编码的交叉操作采用如下方式来实现:(1)abba(16)其中, 、 为区间 上均匀分布的随机数,且有 ,调整 的大小(0,)r 1rr即可控制交叉操
18、作的变化范围。(3)变异操作使用变异算子来调整个体编码串中的部分基因值,可以从局部的角度出发使个体更加逼近最优解,提高遗传算法的局部搜索能力。变异操作将某个个体的参数 ,操作变为域内的另一个值 。浮点数编码的遗传算法采用下式进行:cc(,)N(17)其中, 表示平均值为 ,方差为 的正态分布随机数。可以看出,(,)Ncc变异操作即是以当前值为中心,主要在一个小范围内进行随机扰动的变化。.3 算法设计根据以上介绍,就遗传算法本身而言,其步骤主要如下:(1)确定编码方式,选取一定大小的初始种群;(2)计算所选种群中每一个个体的适应度函数值及复制概率;(3)用轮盘赌方法,淘汰相应个数的函数值较小的个
19、体,替换为相应个数的函数值较大的个体;(4)对个体随机两两配对,按指定概率 进行交叉操作;cP(5)对每一个体中的每一参数,按指定概率 进行变异操作;m(6)若满足收敛条件则输出最优解并退出,否则返回执行步骤(2),如此往复。而根据本机理模型具体实例而言,现将具体的算法流程作以介绍。步骤一:产生待辨识参数的初始种群根据式(12)中所示的待辨识参数上下限范围,采用随机的方式,生成3个参数的初始种群,种群大小为 ,并设置最大进化代数 。若10N10G为种群中某个体的向量表示,则用MATLAB 随机生成的种群可用(,)Tarcp一个 的矩阵来表示。310步骤二:计算个体的适应度选取50炉钢的操作记录
20、(即数据组数 )作为输入数据及温度比较标50n准,对于选中的每个个体,通过 组输入数据的计算结果,根据式(13)计算适应度函数的值,并判断当前的进化代数是否达到设定值。若达到设定值,则结束遗传算法寻优过程,此时最大适应度值对应的个体即为机理模型参数的辨识结果;若进化代数未达到设定值,则继续执行步骤三。步骤三:进行遗传算法操作.根据步骤二中计算出的适应度值与式(14)计算出的选择概率以及预先设定的交叉和变异概率,进行遗传算法的选择、交叉和变异操作,产生下一代群体,而后将遗传代数 增加1再返回步骤二。G根据以上的三个步骤,绘制算法的流程图如下所示:图 5 算法流程图4 仿真实验与结果分析综合前文所
21、述,参与参数优化的输入数据有:总耗电量,冷却水流量,冷却水入口、出口温度,有功功率,精炼周期,加渣、加合金重量,钢包内外径,钢包高度,包底厚度,钢水初温,钢水进站与出站重量等。利用 MATLAB 中英国谢菲尔德大学遗传算法工具箱,很容易利用工具箱内的函数进行生成种群、选择、交叉、变异操作。经 MATLAB 运行仿真,最终遗传代数到达 100 代时得到种群中的最优个体如下:(,)(0.7352,4.13)TTarcp(18)而随着遗传代数的增加,种群中的个体逐步优化,计算温度与实际温度的误差减小,即个体适应度增大,图 6 为每代最优个体适应度值的变化趋势。.图 6 误差趋势图可见,该趋势符合实际
22、情况。得到参数后,将最终得到优化的参数代入到钢水温度预报的机理模型中,即可得到更加精确的温度预报曲线。为检验模型的精确性,另取 20 组现场数据,以式(18)所给出的参数计算终点温度,并与实测温度比较,结果如下表所示。表 2 模型验证结果.为使结果更加直观,用折线图的方式绘制温度计算值与实测值的偏差如下:图 7 预报温度与实际温度误差折线图可见,在20组数据中,有11组的预报结果误差在5以内,有16组的预报结果误差在10以内,而全部预报结果的误差均在15以内,平均预报温度误差为6.2。可以看出,精炼过程中钢水温度的预报值与实测值比较吻合,说明该模型建立即参数优化比较合理。.5 总结与展望本文在
23、钢水温度预报机理模型的基础上,应用遗传算法对机理模型中不易获得的参数进行辨识与寻优,取得了初步的成功。然而,在温度预报方面,虽然本文中的机理模型可以正确预测温度变化趋势,混合模型可以利用智能方法辨识出部分机理模型的参数,使模型更加准确,但由于冶炼现场的复杂性,任何确定性或非确定性因素都有可能影响温度变化的走势。因此,若想建立更加精确的温度预报模型,还需要对现场的实际情况做更加深入的考察与分析。此外,应用遗传算法时,待辨识参数的选择也尤为关键。若选择较多的参数进行辨识,有可能使辨识后模型的结果更贴近于实际情况,但计算相对复杂;而选择较少的参数辨识,可使计算较为简化,但模型精度变差。因此,结合现场
24、经验,合理把握辨识参数的程度亦十分关键。而在执行遗传算法的过程中,若种群中个体的个数为 ,数据的组数为 ,则每代种群需进行 次计算,NnnN效率相对较低,所以研究并开发其它相对高效的智能算法做钢水温度预报研究,将成为我们今后研究方向之一。.附录 MATLAB 源程序代码clcclear alldata=xlsread(C:wrydata.xls);P=data(:,1:end-1);T=data(:,end);N=3;threshold=-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;%*%N
25、IND=100;%MAXGEN=100;PRECI=3;%*GGAP=0.95;px=0.7;pm=0.01;trace=zeros(N+1,MAXGEN);Chrom=crtrp(NIND,0.5 1.3 1000;0.8 6.5 1380);%gen=0;for k=1:MAXGENObj=zeros(size(P,1),NIND);tin=25;R1=1.85;R2=2;Tevn=30;Height=4;Z1=0.25;for j=1:NINDfor i=1:size(P,1)qall=P(i,1);iter=Chrom(j,1);flow=P(i,2);tout=P(i,3);pact
26、ive=P(i,4);RealTime=P(i,5);Cp=Chrom(j,3);lamad=Chrom(j,2);zha=P(i,6);cao=P(i,7);casi=P(i,8);sife=P(i,9);.mnfe=P(i,10);Tst0=P(i,11);Tls0=P(i,12);mstin=P(i,13);mstout=P(i,14);Obj(i,j)=mechanicalmodel(qall,iter,flow,tout,tin,pactive,RealTime,Cp,lamad,zha,cao,casi,sife,mnfe,R1,R2,Tst0,Tls0,Tevn,Height,Z
27、1,mstin,mstout);endendObjV=zeros(1,NIND);fprintf(%dn,gen);for j=1:NINDObjV(j)=sum(Obj(:,j)-T).2);endObjV=ObjV;FitnV=ranking(ObjV);SelCh=select(sus,Chrom,FitnV,GGAP);SelCh=recombin(reclin,SelCh,px);ObjSel=zeros(size(P,1),size(SelCh,1);for j=1:size(SelCh,1)for i=1:size(P,1)qall=P(i,1);iter=SelCh(j,1);
28、flow=P(i,2);tout=P(i,3);pactive=P(i,4);RealTime=P(i,5);Cp=SelCh(j,3);lamad=SelCh(j,2);zha=P(i,6);cao=P(i,7);casi=P(i,8);sife=P(i,9);mnfe=P(i,10);Tst0=P(i,11);Tls0=P(i,12);mstin=P(i,13);.mstout=P(i,14);ObjSel(i,j)=mechanicalmodel(qall,iter,flow,tout,tin,pactive,RealTime,Cp,lamad,zha,cao,casi,sife,mnf
29、e,R1,R2,Tst0,Tls0,Tevn,Height,Z1,mstin,mstout);endendObjSelV=zeros(1,size(SelCh,1);for j=1:size(SelCh,1)ObjSelV(j)=sum(ObjSel(:,j)-T).2);endObjSelV=ObjSelV;Chrom,ObjV=reins(Chrom,SelCh,1,1,ObjV,ObjSelV);gen=gen+1;Y,I=min(ObjV);trace(1:N,gen)=Chrom(I,:);trace(end,gen)=Y;end%outcome=trace;figure(1);pl
30、ot(1:MAXGEN,trace(end,:);grid on;xlabel()ylabel()title()bestX=trace(1:end-1,end);bestErr=trace(end,end);function y=mechanicalmodel(qall,iter,flow,tout,tin,pactive,RealTime,Cp,lamad,zha,cao,casi,sife,mnfe,R1,R2,Tst0,Tls0,Tevn,Height,Z1,mstin,mstout)y1=absorbedheat(qall,iter);y2=surfacedissipation(flo
31、w,tout,tin,pactive,RealTime);y3=materials(zha,cao,casi,sife,mnfe);y4=sidedissipation(R1,R2,Tst0,Tls0,Tevn,Height,RealTime,Cp,lamad);y5=bottomdissipation(Z1,Tst0,Tls0,Tevn,R1,RealTime,Cp,lamad);y=(y1-y2-y4-y5)/(460*mstin+1050*(mstout-mstin)-y3+Tst0;.function y=absorbedheat(qall,iter)y=qall*iter*3.6*1
32、000000;function y=surfacedissipation(flow,tout,tin,pactive,RealTime)y=RealTime/3600*flow*1000*4200*(tout-tin)+pactive*0.065*RealTime/3600*3.6*1000000;function y=materials(zha,cao,casi,sife,mnfe)y=(zha+cao+casi*1.05-sife*0.9+mnfe*0.9)*0.01;function y=sidedissipation(R1,R2,Tst0,Tls0,Tevn,Height,RealTi
33、me,Cp,lamad)deltaR=(R2-R1)/9;for i=1:1:10R(i)=R1+deltaR*(i-1);endT0=zeros(1,10);for i=1:1:10T0(i)=Tls0+(Tst0-Tls0)*(log(R2/R(i)/(log(R2/R1);endlamad1=lamad;deltaTau=10;hac=10;rho=2000;B=hac*deltaR/lamad1;Fo=(lamad1*deltaTau)/(rho*Cp*deltaR*deltaR);T=zeros(270,10);for j=1:1:270for i=1:1:10switch icas
34、e 1T(j,i)=Tst0;case 10T(j,i)=(1-2*Fo*(1-deltaR/(2*R(i)-2*Fo*B)*T0(i)+2*Fo*(1-deltaR/(2*R(i)*T0(i-1)+2*Fo*B*Tevn;otherwiseT(j,i)=Fo*(1-deltaR/(2*R(i)*T0(i-1)+(1-2*Fo)*T0(i)+Fo*(1+deltaR/(2*R(i)*T0(i+1);endendT0=T(j,:);endfor j=1:1:270Qac(j)=lamad1*(T(j,1)-T(j,2)/deltaR*2*pi*R1*Height;endx=RealTime/10
35、;if (x=0).Result=0;elseResult=Qac(fix(x);endy=Result*RealTime;function y=bottomdissipation(Z1,Tst0,Tls0,Tevn,R1,RealTime,Cp,lamad)deltaZ=Z1/9;for i=1:1:10Z(i)=deltaZ*(i-1);endT0=zeros(1,10);for i=1:1:10T0(i)=Tst0-(Tst0-Tls0)*Z(i)/Z1;endlamad2=lamad;deltaTau=10;hab=5;rho=2000;B=hab*deltaZ/lamad2;Fo=(
36、lamad2*deltaTau)/(rho*Cp*deltaZ*deltaZ);T=zeros(270,10);for j=1:1:270for i=1:1:10switch icase 1T(j,i)=Tst0;case 10T(j,i)=(1-2*Fo*(1+B)*T0(i)+2*Fo*T0(i-1)+2*Fo*B*Tevn;otherwiseT(j,i)=Fo*T0(i-1)+(1-2*Fo)*T0(i)+Fo*T0(i+1);endendT0=T(j,:);endfor j=1:1:270Qab(j)=lamad2*(T(j,1)-T(j,2)/deltaZ*2*pi*R1*R1;endx=RealTime/10;if (x=0)Result=0;elseResult=Qab(fix(x);endy=Result*RealTime;