发明名称 循环流化床生活垃圾焚烧锅炉汽包水位的预测系统及方法
摘要 本发明公开了一种循环流化床生活垃圾焚烧锅炉汽包水位的预测系统及方法。本发明先从汽包水位变化特性的机理出发,在合理简化假设的基础上,依据质量守恒方程、能量守恒方程和一些基础的方程,并结合CFB生活垃圾焚烧锅炉水冷壁独特的热量吸收分布规律,建立起适度复杂的基于微分方程形式的汽包水位动态特性机理模型。然后利用ANFIS建模来挖掘运行历史数据中隐含的知识,以补偿汽包水位的机理建模过程当中,由于模型简化假设、模型降阶或线性化处理、对事物变化过程机理的认知不完备、对象特性不同以及内外扰动所带来的误差。充分发挥机理建模和ANFIS建模的优势,提升汽包水位的预测精度。
申请公布号 CN106055520B 申请公布日期 2017.05.03
申请号 CN201610398791.7 申请日期 2016.06.06
申请人 浙江大学 发明人 尤海辉;马增益;唐义军;王月兰;倪明江;严建华
分类号 G06F17/11(2006.01)I;G06F17/50(2006.01)I 主分类号 G06F17/11(2006.01)I
代理机构 杭州求是专利事务所有限公司 33200 代理人 邱启旺
主权项 一种循环流化床生活垃圾焚烧锅炉汽包水位的预测系统,该系统与循环流化床锅炉的集散控制系统相连,包括数据通讯接口和上位机,所述上位机包括:第一信号采集模块;利用该模型采集CFB生活垃圾焚烧锅炉在正常运行时的运行工况状态参数和操作变量,并组成ANIFS汽包水位补偿模型输入变量的训练样本矩阵X(m×n),m表示样本个数,n表示变量的个数;数据预处理模块;对X(m×n)进行粗大误差处理和随机误差处理,剔除训练样本中的野值,排除异常工况,将训练样本输入变量经归一化处理后映射到[0,1]区间内,得到标准化后的训练样本X<sup>*</sup>(m×n);专家知识库模块;将水位计测量出来的实际汽包水位值l(m×1)与机理模型计算得到的水位值l<sub>d</sub>(m×1)相减,将两者之间的差值作为ANFIS汽包水位补偿模型训练样本的输出向量l<sup>*</sup>(m×1),X<sup>*</sup>和l<sup>*</sup>共同组成ANFIS模型的训练样本,并进行保存;混合建模模块;该模块主要包括1)机理建模和2)ANFIS建模两个部分;1)机理建模的步骤如下:1.1)对水冷壁中的工质进行质量守恒、动量守恒和能量守恒分析,可得:<img file="FDA0001221101610000011.GIF" wi="1660" he="118" /><img file="FDA0001221101610000012.GIF" wi="1662" he="191" /><img file="FDA0001221101610000013.GIF" wi="1662" he="222" />其中,z为坐标轴,正向为工质的上升方向;ρ<sub>s</sub>、ρ<sub>l</sub>分别表示研究对象内饱和水蒸气的密度和饱和水的密度,密度的单位是kg/m<sup>3</sup>;α<sub>s</sub>表示截面容积含汽率;v表示工质的流速,单位是m/s;τ<sub>wm</sub>是壁面剪切应力,单位是Pa/m<sup>3</sup>;p<sub>r</sub>是上升管中的压力,单位Pa;<img file="FDA0001221101610000014.GIF" wi="70" he="79" />表示热流密度,单位是W/m<sup>3</sup>;g是重力加速度,单位是m/s<sup>2</sup>;h<sub>l</sub>为饱和水的比焓;饱和温度和工质压力之间存在着以下关系:T=T<sub>b</sub>(1+a<sub>1</sub>ln(P/P<sub>b</sub>)+a<sub>2</sub>ln<sup>2</sup>(P/P<sub>b</sub>)+a<sub>3</sub>ln<sup>3</sup>(P/P<sub>b</sub>)+a<sub>4</sub>ln<sup>4</sup>(P/P<sub>b</sub>))   (5)其中P是当前位置压力,T是当前位置的温度,T<sub>b</sub>参考点的温度,P<sub>b</sub>表示参考点的压力,a<sub>1</sub>、a<sub>2</sub>、a<sub>3</sub>、a<sub>4</sub>是常系数;壁面的剪切应力可以通过下式求取:<img file="FDA0001221101610000021.GIF" wi="1661" he="127" />f=0.0014+0.125Re<sup>‑0.32</sup>          (7)ρ<sub>m</sub>=α<sub>s</sub>ρ<sub>s</sub>+(1‑α<sub>s</sub>)ρ<sub>l</sub>       (8)式中Re表示雷诺数,f为摩擦系数,d<sub>t</sub>为水冷壁上升管的管径,ρ<sub>m</sub>为汽水混合物的密度;水冷壁的热流密度可以通过下式求取:<img file="FDA0001221101610000022.GIF" wi="1660" he="127" />式中,A<sub>t</sub>表示截面积,U<sub>wm</sub>表示传热系数,T<sub>w</sub>表示管壁的温度;对(2)~(4)式进行展开,并将(5)~(9)式代入,可得到如下方程组:<img file="FDA0001221101610000023.GIF" wi="862" he="295" />其中,A<sub>11</sub>=(ρ<sub>s</sub>‑ρ<sub>l</sub>)v,A<sub>12</sub>=α<sub>s</sub>v,A<sub>13</sub>=α<sub>s</sub>ρ<sub>s</sub>+(1‑α<sub>s</sub>)ρ<sub>l</sub>,A<sub>14</sub>=0,A<sub>15</sub>=0;<img file="FDA0001221101610000024.GIF" wi="1230" he="117" /><img file="FDA0001221101610000025.GIF" wi="629" he="133" /><img file="FDA0001221101610000026.GIF" wi="1252" he="130" /><img file="FDA0001221101610000027.GIF" wi="611" he="79" />A<sub>25</sub>=0;A<sub>31</sub>=(ρ<sub>s</sub>‑ρ<sub>l</sub>)v<sup>2</sup>,A<sub>32</sub>=α<sub>s</sub>v<sup>2</sup>,A<sub>33</sub>=2[α<sub>s</sub>ρ<sub>s</sub>+(1‑α<sub>s</sub>)ρ<sub>l</sub>]v,A<sub>34</sub>=0,A<sub>35</sub>=1;A<sub>41</sub>=0,A<sub>42</sub>=0,A<sub>43</sub>=0,A<sub>44</sub>=1,<img file="FDA0001221101610000031.GIF" wi="1265" he="111" />A<sub>51</sub>=0,A<sub>52</sub>=1,A<sub>53</sub>=0,<img file="FDA0001221101610000032.GIF" wi="406" he="110" />y<sub>1</sub>=α<sub>s</sub>,y<sub>2</sub>=ρ<sub>s</sub>,y<sub>3</sub>=v,y<sub>4</sub>=T,y<sub>5</sub>=P;B<sub>1</sub>=0,<img file="FDA0001221101610000033.GIF" wi="940" he="125" /><img file="FDA0001221101610000034.GIF" wi="1210" he="134" />B<sub>4</sub>=0,B<sub>5</sub>=0;式中,λ<sub>V</sub>表示汽化潜热,<img file="FDA0001221101610000035.GIF" wi="46" he="70" />表示饱和水的定压比热,<img file="FDA0001221101610000036.GIF" wi="46" he="62" />表示饱和蒸汽的定压比热,T<sub>0</sub>表示指定状态下工质的温度;最终可以算得:<img file="FDA0001221101610000037.GIF" wi="1662" he="295" />1.2)通过对汽包‑下降管‑水冷壁中的所有工质进行质量守恒分析和能量守恒分析,可得:<img file="FDA0001221101610000038.GIF" wi="1662" he="116" /><img file="FDA0001221101610000039.GIF" wi="1659" he="119" />其中,<img file="FDA00012211016100000310.GIF" wi="502" he="131" />e<sub>1</sub>=ρ<sub>l</sub>‑ρ<sub>s</sub>;<img file="FDA00012211016100000311.GIF" wi="429" he="128" />式中,p表示汽包压力,单位Pa;Q表示整个蒸发区的单位时间内的吸热量,单位是W;V<sub>lt</sub>表示汽包‑下降管‑水冷壁中所有饱和水的体积,单位是m<sup>3</sup>;V<sub>st</sub>表示回路中所有水蒸气的体积,单位是m<sup>3</sup>;q<sub>f</sub>、q<sub>s</sub>分别表示给水流量和主汽流量和下降管中循环工质的流量,流量的单位是kg/s;h<sub>f</sub>、h<sub>l</sub>、h<sub>s</sub>分别表示给水比焓、饱和水的比焓和饱和水蒸气的比焓,单位是kJ/kg;m<sub>t</sub>表示水冷壁金属和汽包金属的总质量,单位是kg;C<sub>p</sub>表示金属的比热容;t<sub>s</sub>表示饱和温度;1.3)对汽包中液面下的工质进行质能守恒分析,忽略液面上饱和蒸汽的凝结,可得:<img file="FDA00012211016100000312.GIF" wi="1661" he="118" /><img file="FDA0001221101610000041.GIF" wi="1660" he="119" />其中,<img file="FDA0001221101610000042.GIF" wi="1422" he="134" /><img file="FDA0001221101610000043.GIF" wi="1397" he="133" />式中,α<sub>r</sub>表示上升管出口的截面质量含汽率;<img file="FDA0001221101610000044.GIF" wi="80" he="63" />表示汽包内液面下工质的质量,单位是kg;q<sub>r</sub>、<img file="FDA0001221101610000045.GIF" wi="79" he="63" />q<sub>cd</sub>、q<sub>dc</sub>分别表示上升管中的流量、饱和蒸汽逸出液面的流量、饱和蒸汽凝结的流量和下降管中的流量,单位是kg/s;t<sub>s</sub>表示饱和蒸汽的温度,单位是℃;m<sub>d</sub>表示汽包金属的质量,单位kg;C<sub>p</sub>表示水冷壁和汽包中金属的比热容,单位是kJ/(kg·℃);A<sub>d</sub>是汽包液位的截面面积,单位是m<sup>2</sup>;<img file="FDA0001221101610000046.GIF" wi="46" he="63" />是饱和蒸汽逸出的速度,单位是m/s;σ是液面的表面张力,单位是N/m;<img file="FDA0001221101610000047.GIF" wi="86" he="61" />V<sub>l</sub>分别表示液面下饱和水蒸气的体积和饱和水的体积,单位是m<sup>3</sup>;联立(12)~(15)式并进行展开,可以得到如下方程组:<img file="FDA0001221101610000048.GIF" wi="806" he="295" />其中,<img file="FDA0001221101610000049.GIF" wi="519" he="127" />M<sub>12</sub>=0,M<sub>13</sub>=0,M<sub>14</sub>=0;<img file="FDA00012211016100000410.GIF" wi="434" he="127" />M<sub>22</sub>=ρ<sub>l</sub>,M<sub>23</sub>=ρ<sub>s</sub>,M<sub>24</sub>=0;<img file="FDA00012211016100000411.GIF" wi="286" he="127" />M<sub>32</sub>=0,M<sub>33</sub>=ρ<sub>s</sub>,M<sub>34</sub>=0;<img file="FDA00012211016100000412.GIF" wi="455" he="127" />M<sub>42</sub>=0,M<sub>43</sub>=0,M<sub>44</sub>=ρ<sub>l</sub>‑ρ<sub>s</sub>;N<sub>1</sub>=Q‑q<sub>f</sub>(h<sub>l</sub>‑h<sub>f</sub>)‑q<sub>s</sub>(h<sub>s</sub>‑h<sub>l</sub>),<img file="FDA00012211016100000413.GIF" wi="921" he="69" />N<sub>4</sub>=q<sub>f</sub>‑q<sub>s</sub>;o<sub>1</sub>=p,o<sub>2</sub>=V<sub>l</sub>,o<sub>3</sub>=V<sub>sd</sub>,o<sub>4</sub>=V<sub>lt</sub>;最终可得:<img file="FDA0001221101610000051.GIF" wi="1660" he="299" />采用Runge‑Kutta数值解法求解微分方程(11)和(16);对仿真时间进行离散化处理,将[0,T<sub>0</sub>]区间内的仿真时间,离散成n个时间间隔为ΔT的区间,离散后的区间表示为[0,ΔT],[ΔT,2ΔT],……,[(n‑1)ΔT,nΔT];ΔT为0.01s~0.1s;1.4)汽包水位是由汽包中的饱和水和饱和水蒸气的体积共同决定的,忽略汽包体积和汽包水位之间的非线性,可得汽包水位的表达式:<img file="FDA0001221101610000052.GIF" wi="1421" he="132" />式中,l<sub>d</sub>表示机理模型计算得到的汽包水位,单位mm;l<sub>0</sub>是汽包水位的基准水位线;2)ANFIS建模的步骤如下:2.1)初始化粒子群;以聚类半径r<sub>α</sub>作为粒子,15个粒子作为一个种群,每个粒子随机赋予[0.20.9]区间内的随机值,其中第i个粒子的位置的向量标示为ri,i=1,2,…,15;2.2)以r<sub>i</sub>为聚类半径,进行减法聚类分析,将每个数据点作为可能的聚类中心,并根据各个数据点周围的数据点密度来计算该点作为聚类中心的可能性;每个数据点X<sub>i</sub>作为聚类中心的可能性P<sub>i</sub>由式(18)来定义:<img file="FDA0001221101610000053.GIF" wi="1667" he="171" />式中m表示n维输入空间中全部的数据点数,X<sub>i</sub>=[X<sub>i1</sub>,X<sub>i2</sub>,...,X<sub>in</sub>]、X<sub>j</sub>=[X<sub>j1</sub>,X<sub>j2</sub>,...,X<sub>jn</sub>]是具体的数据点,r<sub>i</sub>是一个正数,定义了该点的邻域半径,||·||符号表示欧式距离;被选为聚类中心的点具有最高的数据点密度,同时该该数据点周围的点被排除作为聚类中心的可能性;第一个聚类中心为X<sub>C1</sub>,数据点密度为P<sub>c1</sub>;选出第一个聚类中心后,继续采用类似的方法确定下一个聚类中心,但需消除已有聚类中心的影响,修改密度指标的山峰函数如下:<img file="FDA0001221101610000054.GIF" wi="1686" he="163" />其中,r<sub>β</sub>=1.5r<sub>i</sub>;循环重复上述过程直到所有剩余数据点作为聚类中心的可能性低于 某一阈值δ,即P<sub>ck</sub>/P<sub>c1</sub><δ;2.3)ANFIS模型训练;根据减法聚类算法得到的聚类中心,按照ANFIS模型结构训练汽包水位补偿模型;对于模糊神经网络模型的所有参数,采用混合最小二乘法的梯度下降算法进行学习;2.4)计算适应度值;利用训练得到的预测模型计算汽包水位补偿值l<sup>*</sup>,并与实际补偿值l<sub>e</sub>进行比较,并以误差平方和作为粒子的适应度值MSE,适应度计算公式如下:<img file="FDA0001221101610000061.GIF" wi="1686" he="126" />2.5)更新极值;以适应度值为评价指标,比较当代粒子与上一代粒子之间的适应度值大小,如果当前粒子的适应度值优于上一代,则将当前粒子的位置设置为个体极值,否则个体极值保持不变;同时获取当代所有粒子适应度值最优的粒子,并与上一代最优粒子进行比较,如果当代最优粒子的适应度值优于上一代最优粒子的适应度值,则将当代粒子的最优适应度值设置为全局最优值,否则全局最优值保持不变;2.6)更新粒子;根据最新的个体极值和全局极值,按照(29)式和(30)式更新粒子的速度v<sub>id</sub>(t)和位置x<sub>id</sub>(t);v<sub>id</sub>(t+1)=ωv<sub>id</sub>(t)+c<sub>1</sub>r<sub>1</sub>(p<sub>id</sub>‑x<sub>id</sub>(t))+c<sub>2</sub>r<sub>2</sub>(p<sub>gd</sub>‑x<sub>id</sub>(t))      (29)x<sub>id</sub>(t+1)=x<sub>id</sub>(t)+v<sub>id</sub>(t+1)          (30)t是粒子群优化算法的寻优代数;更进一步,为了改善基本粒子群算法容易陷入局部极值和收敛速度慢的缺陷,在PSO算法的基础上引进了动态加速常数c<sub>1</sub>、c<sub>2</sub>和惯性权重ω:<img file="FDA0001221101610000062.GIF" wi="1662" he="119" /><img file="FDA0001221101610000063.GIF" wi="1662" he="119" /><img file="FDA0001221101610000064.GIF" wi="1662" he="117" />其中,T<sub>max</sub>为最大寻优代数,ω<sub>max</sub>为最大惯性权重,ω<sub>min</sub>为最小惯性权重,R<sub>1</sub>、R<sub>2</sub>、R<sub>3</sub>、R<sub>4</sub>为常数;2.7)算法停止条件算法判定;判断是否达到最大迭代次数或者到达预测精度的要求,如果没有达到则返回步骤2.2),利用更新的聚类半径继续搜索,否则退出搜索;2.8)利用最终寻优得到的聚类半径,对样本进行聚类分析和ANFIS模型训练,得到 达到训练终止条件的ANFIS模型,即汽包水位的补偿模型;第二信号采集模块;用于从数据库中选择需要预测汽包水位的运行工况,或者实时地采集当前锅炉的运行工况;预测模块;该模块用于对指定的样本进行汽包水位的预测,或者对当前锅炉运行工况下的汽包水位进行实时预测;结果显示模块;显示汽包水位的预测结果,或者对汽包水位的预测结果进行统计分析。
地址 310027 浙江省杭州市西湖区浙大路38号