CN106055520A - 循环流化床生活垃圾焚烧锅炉汽包水位的预测系统及方法 - Google Patents

Abstract

本发明公开了一种循环流化床生活垃圾焚烧锅炉汽包水位的预测系统及方法。本发明先从汽包水位变化特性的机理出发，在合理简化假设的基础上，依据质量守恒方程、能量守恒方程和一些基础的方程，并结合CFB生活垃圾焚烧锅炉水冷壁独特的热量吸收分布规律，建立起适度复杂的基于微分方程形式的汽包水位动态特性机理模型。然后利用ANFIS建模来挖掘运行历史数据中隐含的知识，以补偿汽包水位的机理建模过程当中，由于模型简化假设、模型降阶或线性化处理、对事物变化过程机理的认知不完备、对象特性不同以及内外扰动所带来的误差。充分发挥机理建模和ANFIS建模的优势，提升汽包水位的预测精度。

Description

循环流化床生活垃圾焚烧锅炉汽包水位的预测系统及方法

技术领域

本发明涉及能源工程领域，特别地，涉及一种循环流化床生活垃圾焚烧锅炉汽包水位的预测系统及方法。

背景技术

垃圾焚烧由于能够良好实现垃圾处理技术的减容化、减量化、无害化和资源化，近十几年内，在国家相关产业政策的引导下，国内垃圾焚烧行业取得了蓬勃的发展。循环流化床(Circulating Fluidized Bed,CFB)生活垃圾焚烧技术作为主要的生活垃圾焚烧处理技术之一，已经在国内的多个城市进行了推广应用，截止2014年年底，国内已建成垃圾焚烧锅炉70余台，日处理垃圾量6.4万吨，为我国的垃圾焚烧处理行业做出了重要的贡献。汽包水位作为锅炉重要的运行控制目标之一，是一个反映锅炉内能量和工质平衡的关键参数，对锅炉的安全运行造成了严重的影响。汽包水位过高会引起蒸汽带水，恶化蒸汽的品质，并有可能会损坏管道、汽轮机等设备，汽包水位过低，会导致下降管出口带汽，影响正常水循环，甚至造成汽包干锅和水冷壁烧损。由于汽包液面下存在大量的水蒸气，并且不断地进行膨胀和收缩，使汽包水位的建模和控制变得困难。

目前，国内外并没有针对CFB生活垃圾焚烧锅炉汽包水位的动态特性进行专门的研究，只是针对常规电站的燃煤锅炉进行了研究，并提出了集中参数的机理模型。这些模型通过对汽包、下降管和水冷壁分别进行质量守恒、动量守恒和能量守恒分析，并在简化假设的基础上构建了基于微分方程的动态特性模型，可以用于对汽包水位进行仿真分析，或服务于一些基于模型的控制算法。这些单纯依靠汽包水位建模方法，很好的描述了汽包水位动态变化的本质，解释性好，外推性强，大多数的工况下都能够反映汽包水位的变化趋势。然而这种建模方法却有以下不足之处：

1)机理建模是建立在诸多的简化和假设的基础之上的，与实际的模型和物理过程之间有着一定的误差，必然会导致模型的预测结构和实际值有出入；

2)机理建模方法得到的模型比较复杂，往往都是高阶模型，或者具有非线性，为了更好的求解和实际应用，通常会进行降阶或者线性化处理，从而使模型的精度降低；

3)汽包水位的形成过程牵涉到了复杂的两相流和气泡的生成-膨胀/收缩-破裂过程，人们对这些过程的研究和理解并不是非常的完备，因此会在建模的过程无法完全表达出这些过程的特性，甚至进行错误的数学描述，从而给模型带来误差；

4)这些机理模型都是针对燃煤电站锅炉的，并且是在假设水冷壁均匀吸热的前提下提出的。然而国内的城市生活垃圾多为混合收集，入厂、入炉垃圾成分较为复杂，表现出低热偏值、水分较高和波动性较大的特征，并且垃圾的理化性质和煤相比有截然不同的特性，因此，在CFB锅炉中燃烧时具有独特的燃料热值释放规律。另外，生产工艺的特点也加剧了CFB垃圾焚烧锅炉中水冷壁的吸热的不均匀性。因此，这些基于水冷壁均匀受热的汽包水位变化模型不一定适用于CFB生活垃圾焚烧锅炉。

5)实际的生产过程中，建模对象和物理过程总是面临着各种各样可知或者不可知的扰动，而机理建模的过程中通常会忽略这些扰动，导致创建的模型过于理想化。

在使用机理建模方法创建汽包水位的预测模型时存在的这些不足之处可以通过黑箱建模的方法来进行弥补。随着电子技术、计算机技术和信息技术的发展，集散控制系统(Distributed Control System,DCS)广泛的应用于CFB生活焚烧锅炉的运行过程，包含温度、压力、流量等参数在内的过程数据都被完善得保存下来，这些历史数据中包含丰富过程信息，是人们认识和了解生产过程的重要途径之一，具有很高的挖掘价值。自适应糊神经网络(AdaptiveNeuro-Fuzzy Inference System,ANFIS)融合了神经网络并行计算、分布式信息存储、容错能力强、具备自适应学习功能等优点和模糊算法能够有效表达先验知识的能力，成为构建不确定性和非线性模型的有力工具。ANFIS建模时无需掌握复杂的对象机理，建模过程也相对简单，能够表征那些尚未被人们掌握的过程特性和各种扰动。但是，单纯的ANFIS建模方法存在着模型外推性差、模型的解释性不强、依赖训练样本并且涵盖的工况有限等不足之处。

发明内容

本发明的目的在于针对现有技术的不足，提供一种循环流化床生活垃圾焚烧锅炉汽包水位的预测系统及方法，本发明将机理建模和ANFIS建模联合起来，创建一个基于混合建模方法的汽包水位动态特性表征模型。利用ANFIS建模来补偿汽包水位的机理建模过程当中，由于模型简化假设、模型降阶或线性化处理、对事物变化过程机理的认知不完备、对象特性不同以及内外扰动所带来的误差。充分发挥机理建模和ANFIS建模的优势，提升汽包水位模型的预测精度。

本发明解决其技术问题所采用的技术方案是：一种循环流化床生活垃圾焚烧锅炉汽包水位的预测系统，该系统与循环流化床锅炉的集散控制系统相连，包括数据通讯接口和上位机，所述上位机包括：

第一信号采集模块。利用该模型采集CFB生活垃圾焚烧锅炉在正常运行时的运行工况状态参数和操作变量，并组成ANIFS汽包水位补偿模型输入变量的训练样本矩阵X(m×n)，m表示样本个数，n表示变量的个数。

数据预处理模块。对X(m×n)进行粗大误差处理和随机误差处理，剔除训练样本中的野值，排除异常工况，将训练样本输入变量经归一化处理后映射到[0，1]区间内，得到标准化后的训练样本X^*(m×n)。

专家知识库模块。将水位计测量出来的实际汽包水位值l(m×1)与机理模型计算得到的水位值l_d(m×1)相减，将两者之间的差值作为ANFIS汽包水位补偿模型训练样本的输出向量l^*(m×1)，X^*和l^*共同组成ANFIS模型的训练样本，并进行保存。

混合建模模块。该模块主要包括机理建模和ANFIS建模两个部分。机理建模的步骤如下：

2.1)对水冷壁中的工质进行质量守恒、动量守恒和能量守恒分析，可得：

$\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v\]=0---\left(2\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s{v}^2+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l{v}^2\]+\frac{{\mathrm{dp}}_r}{dz}\\ =-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]g-{\mathrm{\τ}}_{wm}\end{array}---\left(3\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v(h_s+\frac{1}{2}{v}^2)\]+\frac{d}{dz}\[(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v(h_l+\frac{1}{2}{v}^2)\]\\ ={\stackrel{\·}{q}}_{wm}-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]vg\end{array}---\left(4\right)$

其中，z为坐标轴，正向为工质的上升方向；ρ_s、ρ_l分别表示研究对象内饱和水蒸气的密度和饱和水的密度，密度的单位是kg/m³；α_s表示截面容积含汽率；v表示工质的流速，单位是m/s；τ_wm是壁面剪切应力，单位是Pa/m³；p_r是上升管中的压力，单位Pa；表示热流密度，单位是W/m³；g是重力加速度，单位是m/s²。h_l为饱和水的比焓。

饱和温度和工质压力之间存在着以下关系：

T＝T_b(1+a₁ln(P/P_b)+a₂ln²(P/P_b)+a₃ln³(P/P_b)+a₄ln⁴(P/P_b)) (5)

其中P是当前位置压力，T是当前位置的温度，T_b参考点的温度，P_b表示参考点的压力，a₁、a₂、a₃、a₄是常系数。

壁面的剪切应力可以通过下式求取：

${\mathrm{\τ}}_{wm}=\frac{2f}{d_t}{\mathrm{\ρ}}_m{v}^2---\left(6\right)$

f＝0.0014+0.125Re^-0.32 (7)

ρ_m＝α_sρ_s+(1-α_s)ρ_l (8)

式中Re表示雷诺数，f为摩擦系数，d_t为水冷壁上升管的管径，ρ_m为汽水混合物的密度。

水冷壁的热流密度可以通过下式求取：

${\stackrel{\·}{q}}_{wm}=\left(\frac{{\mathrm{\πd}}_t}{A_t}\right)U_{wm}(T_w-T)---\left(9\right)$

式中，A_t表示截面积，U_wm表示传热系数，T_w表示管壁的温度。

对(2)～(4)式进行展开，并将(5)～(9)式代入，可得到如下方程组：

$\left[\begin{array}{ccccc}{A}_{11}& A_{12}& ...& ...& A_{15}\\ A_{21}& A_{22}& ...& ...& A_{25}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ A_{51}& A_{52}& ...& ...& A_{55}\end{array}\right]\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=\left[\begin{array}{c}{B}_1\\ B_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ B_5\end{array}\right]---\left(10\right)$

其中，A₁₁＝(ρ_s-ρ_l)v，A₁₂＝α_sv，A₁₃＝α_sρ_s+(1-α_s)ρ_l，A₁₄＝0，A₁₅＝0；

$A_{21}=\{{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v-\{{\mathrm{\ρ}}_l\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v,$

$A_{22}={\mathrm{\α}}_{s}v\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\],$

$A_{23}={\mathrm{\α}}_s{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{3}{2}{v}^2\]+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l\[c_p^l(T-T_0)+\frac{3}{2}{v}^2\],$

A₂₅＝0；A₃₁＝(ρ_s-ρ_l)v²，A₃₂＝α_sv²，

A₃₃＝2[α_sρ_s+(1-α_s)ρ_l]v，A₃₄＝0，A₃₅＝1；A₄₁＝0，A₄₂＝0，A₄₃＝0，A₄₄＝1，

A₅₁＝0，A₅₂＝1，A₅₃＝0，

y₁＝α_s，y₂＝ρ_s，y₃＝v，y₄＝T，y₅＝P；B₁＝0，

$B_2=\left(\frac{4}{d_t}\right)U_{wm}(T_w-T)-({\mathrm{\α}}_s{\mathrm{\ρ}}_s+(1-{\mathrm{\α}}_s\left){\mathrm{\ρ}}_l\right)vg,$

B₄＝0，B₅＝0。式中，λ_V表示汽化潜热，表示饱和水的定压比热，表示饱和蒸汽的定压比热，T₀表示指定状态下工质的温度。最终可以算得：

$\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=A^{-1}B---\left(11\right)$

2.2)通过对汽包-下降管-水冷壁中的所有工质进行质量守恒分析和能量守恒分析，可得：

$k_1\frac{dp}{dt}=Q-q_f(h_l-h_f)-q_s(h_s-h_l)---\left(12\right)$

$e_1\frac{{\mathrm{dV}}_{lt}}{dt}+e_2\frac{dp}{dt}=q_f-q_s---\left(13\right)$

其中，e₁＝ρ_l-ρ_s；

式中，p表示汽包压力，单位Pa；Q表示整个蒸发区的单位时间内的吸热量，单位是W；V_lt表示汽包-下降管-水冷壁中所有饱和水的体积，单位是m³；V_st表示回路中所有水蒸气的体积，单位是m³；q_f、q_s分别表示给水流量和主汽流量和下降管中循环工质的流量，流量的单位是kg/s；h_f、h_l、h_s分别表示给水比焓、饱和水的比焓和饱和水蒸气的比焓，单位是kJ/kg；m_t表示水冷壁金属和汽包金属的总质量，单位是kg；C_p表示金属的比热容；t_s表示饱和温度。

2.3)对汽包中液面下的工质进行质能守恒分析，忽略液面上饱和蒸汽的凝结，可得：

$\frac{d\left({\mathrm{\ρ}}_s{V}_s^{+}\right)}{dt}={\mathrm{\α}}_r{q}_r-q_s^{+}-q_{cd}---\left(14\right)$

$\frac{{\mathrm{dM}}_{dt}^{+}}{dt}=q_f+q_r-q_s^{+}-q_{dc}---\left(15\right)$

其中，

式中，α_r表示上升管出口的截面质量含汽率；表示汽包内液面下工质的质量，单位是kg；q_r、q_cd、q_dc分别表示上升管中的流量、饱和蒸汽逸出液面的流量、饱和蒸汽凝结的流量和下降管中的流量，单位是kg/s；t_s表示饱和蒸汽的温度，单位是℃；m_d表示汽包金属的质量，单位kg；C_p表示水冷壁和汽包中金属的比热容，单位是kJ/(kg·℃)；A_d是汽包液位的截面面积，单位是m²；是饱和蒸汽逸出的速度，单位是m/s；σ是液面的表面张力，单位是N/m；V_l分别表示液面下饱和水蒸气的体积和饱和水的体积，单位是m³。

联立(12)～(15)式并进行展开，可以得到如下方程组：

$\left[\begin{array}{ccccc}{M}_{11}& M_{12}& ...& ...& M_{14}\\ M_{21}& M_{22}& ...& ...& M_{24}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ M_{41}& M_{42}& ...& ...& M_{44}\end{array}\right]\left[\begin{array}{c}{\mathrm{do}}_1/dt\\ {\mathrm{do}}_2/dt\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dt\end{array}\right]=\left[\begin{array}{c}{N}_1\\ N_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ N_4\end{array}\right]$

其中，M₁₂＝0，M₁₃＝0，M₁₄＝0；M₂₂＝ρ_l，M₂₃＝ρ_s，M₂₄＝0；M₃₂＝0，M₃₃＝ρ_s，M₃₄＝0；

M₄₂＝0，M₄₃＝0，M₄₄＝ρ_l-ρ_s；

N₁＝Q-q_f(h_l-h_f)-q_s(h_s-h_l)，

N₄＝q_f-q_s；o₁＝p，o₂＝V_l，o₃＝V_sd，o₄＝V_lt。

最终可得：

$\left[\begin{array}{c}{\mathrm{do}}_1/dz\\ {\mathrm{do}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dz\end{array}\right]=M^{-1}N---\left(16\right)$

采用Runge-Kutta数值解法求解微分方程(11)和(16)。对仿真时间进行离散化处理，将[0，T₀]区间内的仿真时间，离散成n个时间间隔为ΔT的区间，离散后的区间表示为[0,ΔT]，[ΔT，2ΔT]，……，[(n-1)ΔT，nΔT]。ΔT为0.01s～0.1s。

2.4)汽包水位是由汽包中的饱和水和饱和水蒸气的体积共同决定的，忽略汽包体积和汽包水位之间的非线性，可得汽包水位的表达式：

$l_d=\frac{V_l+V_s^{+}}{A_d}-l_0---\left(17\right)$

式中，l_d表示机理模型计算得到的汽包水位，单位mm；l₀是汽包水位的基准水位线。

ANFIS建模的步骤如下：

2.a)初始化粒子群。以聚类半径r_α作为粒子，15个粒子作为一个种群，每个粒子随机赋予[0.20.9]区间内的随机值，其中第i个粒子的位置的向量标示为r_i，i＝1，2，…，15；

2.b)以r_i为聚类半径，进行减法聚类分析，将每个数据点作为可能的聚类中心，并根据各个数据点周围的数据点密度来计算该点作为聚类中心的可能性。

每个数据点X_i作为聚类中心的可能性P_i由式(18)来定义：

$P_i=\underset{j=1}{\overset{m}{\Σ}}\exp (-\frac{\left|\right|X_i-X_j|{|}^2}{{(r_i/2)}^2})---\left(18\right)$

式中m表示n维输入空间中全部的数据点数，X_i＝[X_i1,X_i2,...,X_in]、X_j＝[X_j1,X_j2,...,X_jn]是具体的数据点，r_i是一个正数，定义了该点的邻域半径，||·||符号表示欧式距离。被选为聚类中心的点具有最高的数据点密度，同时该该数据点周围的点被排除作为聚类中心的可能性。第一个聚类中心为X_C1，数据点密度为P_c1。选出第一个聚类中心后，继续采用类似的方法确定下一个聚类中心，但需消除已有聚类中心的影响，修改密度指标的山峰函数如下：

$P_i=P_i-P_{c1}\exp (-\frac{\left|\right|X_i-X_{c1}|{|}^2}{{(r_{\mathrm{\β}}/2)}^2})---\left(19\right)$

其中，r_β＝1.5r_i。循环重复上述过程直到所有剩余数据点作为聚类中心的可能性低于某一阈值δ，即P_ck/P_c1＜δ。

2.c)ANFIS模型训练。根据减法聚类算法得到的聚类中心，按照ANFIS模型结构训练汽包水位补偿模型；对于模糊神经网络模型的所有参数，采用混合最小二乘法的梯度下降算法进行学习。

2.d)计算适应度值。利用训练得到的预测模型计算汽包水位补偿值l^*，并与实际补偿值l_e进行比较，并以误差平方和作为粒子的适应度值MSE，适应度计算公式如下：

$MSE=\underset{i=1}{\overset{m}{\Σ}}{(l_e-l^{*})}^2---\left(28\right)$

2.e)更新极值。以适应度值为评价指标，比较当代粒子与上一代粒子之间的适应度值大小，如果当前粒子的适应度值优于上一代，则将当前粒子的位置设置为个体极值，否则个体极值保持不变。同时获取当代所有粒子适应度值最优的粒子，并与上一代最优粒子进行比较，如果当代最优粒子的适应度值优于上一代最优粒子的适应度值，则将当代粒子的最优适应度值设置为全局最优值，否则全局最优值保持不变。

2.f)更新粒子。根据最新的个体极值和全局极值，按照(29)式和(30)式更新粒子的速度v_id(t)和位置x_id(t)。

v_id(t+1)＝ωv_id(t)+c₁r₁(p_id-x_id(t))+c₂r₂(p_gd-x_id(t)) (29)

x_id(t+1)＝x_id(t)+v_id(t+1) (30)

t是粒子群优化算法的寻优代数。更进一步，为了改善基本粒子群算法容易陷入局部极值和收敛速度慢的缺陷，在PSO算法的基础上引进了动态加速常数c₁、c₂和惯性权重ω：

$\mathrm{\ω}={\mathrm{\ω}}_{max}-\frac{{\mathrm{\ω}}_{max}-{\mathrm{\ω}}_{min}}{T_{\max }}t---\left(31\right)$

$c_1=R_1+\frac{R_2\×t}{T_{max}}---\left(32\right)$

$c_2=R_3-\frac{R_4\×t}{T_{max}}---\left(33\right)$

其中，T_max为最大寻优代数，ω_max为最大惯性权重，ω_min为最小惯性权重，R₁、R₂、R₃、R₄为常数。

2.g)算法停止条件算法判定。判断是否达到最大迭代次数或者到达预测精度的要求，如果没有达到则返回步骤2.b)，利用更新的聚类半径继续搜索，否则退出搜索。

2.h)利用最终寻优得到的聚类半径，对样本进行聚类分析和ANFIS模型训练，得到达到训练终止条件的ANFIS模型，即汽包水位的补偿模型。

第二信号采集模块。用于从数据库中选择需要预测汽包水位的运行工况，或者实时地采集当前锅炉的运行工况。

预测模块。该模块用于对指定的样本进行汽包水位的预测，或者对当前锅炉运行工况下的汽包水位进行实时预测。

结果显示模块。显示汽包水位的预测结果，或者对汽包水位的预测结果进行统计分析。

一种循环流化床生活垃圾焚烧入炉汽包水位的预测方法，该方法包括以下步骤：

1)选择混合模型的输入变量并采集训练样本。混合模型包括机理模型和ANFIS补偿模型。输入变量包括锅炉的设计参数、工质的物性参数和运行参数。设计参数包括水冷壁的当量管径、水冷壁的有效吸热长度、水冷壁的容积、循环水的流量、汽包的容积、汽包的平均截面面积、下降管的容积、汽包-下降管-水冷壁三者总共的金属质量。工质的物性参数包括压力、饱和温度、比焓、比容、汽化潜热、密度。运行参数包括给水流量、主蒸汽流量、垃圾给料量、辅助燃煤量、一次风量、二次风量以及炉膛中所有的温度测点。采集CFB生活垃圾焚烧锅炉在正常运行时的运行工况状态参数和操作变量，并组成混合模型输入变量的训练样本矩阵X(m×n)，m表示样本个数，n表示变量的个数。

2)训练样本预处理。对训练样本X(m×n)进行粗大误差处理和随机误差处理，剔除训练样本中的野值，排除异常工况，将训练样本输入变量经归一化处理后映射到[0，1]区间内，得到归一化后的训练样本X^*(m×n)；所述的异常工况包括锅炉停炉运行工况、炉膛压火运行状况和给料机堵塞工况；

3)利用机理建模的方法构建CFB生活垃圾焚烧锅炉汽包水位的机理模型。

3.1)对水冷壁中的工质进行质量守恒、动量守恒和能量守恒分析，可得：

$\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v\]=0---\left(2\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s{v}^2+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l{v}^2\]+\frac{{\mathrm{dp}}_r}{dz}\\ =-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]g-{\mathrm{\τ}}_{wm}\end{array}---\left(3\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v(h_s+\frac{1}{2}{v}^2)\]+\frac{d}{dz}\[(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v(h_l+\frac{1}{2}{v}^2)\]\\ ={\stackrel{\·}{q}}_{wm}-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]vg\end{array}---\left(4\right)$

其中，z为坐标轴，正向为工质的上升方向；ρ_s、ρ_l分别表示研究对象内饱和水蒸气的密度和饱和水的密度，密度的单位是kg/m³；α_s表示截面容积含汽率；v表示工质的流速，单位是m/s；τ_wm是壁面剪切应力，单位是Pa/m³；p_r是上升管中的压力，单位Pa；表示热流密度，单位是W/m³；g是重力加速度，单位是m/s²。h_l为饱和水的比焓。

饱和温度和工质压力之间存在着以下关系：

T＝T_b(1+a₁ln(P/P_b)+a₂ln²(P/P_b)+a₃ln³(P/P_b)+a₄ln⁴(P/P_b)) (5)

其中P是当前位置压力，T是当前位置的温度，T_b参考点的温度，P_b表示参考点的压力，a₁、a₂、a₃、a₄是常系数。

壁面的剪切应力可以通过下式求取：

${\mathrm{\τ}}_{wm}=\frac{2f}{d_t}{\mathrm{\ρ}}_m{v}^2---\left(6\right)$

f＝0.0014+0.125Re^-0.32 (7)

ρ_m＝α_sρ_s+(1-α_s)ρ_l (8)

式中Re表示雷诺数，f为摩擦系数，d_t为水冷壁上升管的管径，ρ_m为汽水混合物的密度。

水冷壁的热流密度可以通过下式求取：

${\stackrel{\·}{q}}_{wm}=\left(\frac{{\mathrm{\πd}}_t}{A_t}\right)U_{wm}(T_w-T)---\left(9\right)$

式中，A_t表示截面积，U_wm表示传热系数，T_w表示管壁的温度。

对(2)～(4)式进行展开，并将(5)～(9)式代入，可得到如下方程组：

$\left[\begin{array}{ccccc}{A}_{11}& A_{12}& ...& ...& A_{15}\\ A_{21}& A_{22}& ...& ...& A_{25}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ A_{51}& A_{52}& ...& ...& A_{55}\end{array}\right]\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=\left[\begin{array}{c}{B}_1\\ B_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ B_5\end{array}\right]---\left(10\right)$

其中，A₁₁＝(ρ_s-ρ_l)v，A₁₂＝α_sv，A₁₃＝α_sρ_s+(1-α_s)ρ_l，A₁₄＝0，A₁₅＝0；

$A_{21}=\{{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v-\{{\mathrm{\ρ}}_l\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v,$

$A_{22}={\mathrm{\α}}_{s}v\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\],$

$A_{23}={\mathrm{\α}}_s{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{3}{2}{v}^2\]+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l\[c_p^l(T-T_0)+\frac{3}{2}{v}^2\],$

A₂₅＝0；A₃₁＝(ρ_s-ρ_l)v²，A₃₂＝α_sv²，

A₃₃＝2[α_sρ_s+(1-α_s)ρ_l]v，A₃₄＝0，A₃₅＝1；A₄₁＝0，A₄₂＝0，A₄₃＝0，A₄₄＝1，

A₅₁＝0，A₅₂＝1，A₅₃＝0，

y₁＝α_s，y₂＝ρ_s，y₃＝v，y₄＝T，y₅＝P；B₁＝0，

$B_2=\left(\frac{4}{d_t}\right)U_{wm}(T_w-T)-({\mathrm{\α}}_s{\mathrm{\ρ}}_s+(1-{\mathrm{\α}}_s\left){\mathrm{\ρ}}_l\right)vg,$

B₄＝0，B₅＝0。式中，λ_V表示汽化潜热，表示饱和水的定压比热，表示饱和蒸汽的定压比热，T₀表示指定状态下工质的温度。最终可以算得：

$\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=A^{-1}B---\left(11\right)$

3.2)通过对汽包-下降管-水冷壁中的所有工质进行质量守恒分析和能量守恒分析，可得：

$k_1\frac{dp}{dt}=Q-q_f(h_l-h_f)-q_s(h_s-h_l)---\left(12\right)$

$e_1\frac{{\mathrm{dV}}_{lt}}{dt}+e_2\frac{dp}{dt}=q_f-q_s---\left(13\right)$

其中，e₁＝ρ_l-ρ_s；

式中，p表示汽包压力，单位Pa；Q表示整个蒸发区的单位时间内的吸热量，单位是W；V_lt表示汽包-下降管-水冷壁中所有饱和水的体积，单位是m³；V_st表示回路中所有水蒸气的体积，单位是m³；q_f、q_s分别表示给水流量和主汽流量和下降管中循环工质的流量，流量的单位是kg/s；h_f、h_l、h_s分别表示给水比焓、饱和水的比焓和饱和水蒸气的比焓，单位是kJ/kg；m_t表示水冷壁金属和汽包金属的总质量，单位是kg；C_p表示金属的比热容；t_s表示饱和温度。

3.3)对汽包中液面下的工质进行质能守恒分析，忽略液面上饱和蒸汽的凝结，可得：

$\frac{d\left({\mathrm{\ρ}}_s{V}_s^{+}\right)}{dt}={\mathrm{\α}}_r{q}_r-q_s^{+}-q_{cd}---\left(14\right)$

$\frac{{\mathrm{dM}}_{dt}^{+}}{dt}=q_f+q_r-q_s^{+}-q_{dc}---\left(15\right)$

其中，

$q_s^{+}=A_d\left(\frac{V_s^{+}}{V_l+V_s^{+}}\right)v_s^{+}{\mathrm{\ρ}}_s;v_s^{+}=1.41{\[\frac{\mathrm{\σ}g({\mathrm{\ρ}}_l-{\mathrm{\ρ}}_s)}{{\mathrm{\ρ}}_l^2}\]}^{1/4};M_{dt}^{+}=V_l{\mathrm{\ρ}}_l+V_s^{+}{\mathrm{\ρ}}_s.$

式中，α_r表示上升管出口的截面质量含汽率；表示汽包内液面下工质的质量，单位是kg；q_r、q_cd、q_dc分别表示上升管中的流量、饱和蒸汽逸出液面的流量、饱和蒸汽凝结的流量和下降管中的流量，单位是kg/s；t_s表示饱和蒸汽的温度，单位是℃；m_d表示汽包金属的质量，单位kg；C_p表示水冷壁和汽包中金属的比热容，单位是kJ/(kg·℃)；A_d是汽包液位的截面面积，单位是m²；是饱和蒸汽逸出的速度，单位是m/s；σ是液面的表面张力，单位是N/m；V_l分别表示液面下饱和水蒸气的体积和饱和水的体积，单位是m³。

联立(12)～(15)式并进行展开，可以得到如下方程组：

$\left[\begin{array}{ccccc}{M}_{11}& M_{12}& ...& ...& M_{14}\\ M_{21}& M_{22}& ...& ...& M_{24}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ M_{41}& M_{42}& ...& ...& M_{44}\end{array}\right]\left[\begin{array}{c}{\mathrm{do}}_1/dt\\ {\mathrm{do}}_2/dt\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dt\end{array}\right]=\left[\begin{array}{c}{N}_1\\ N_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ N_4\end{array}\right]$

其中，M₁₂＝0，M₁₃＝0，M₁₄＝0；M₂₂＝ρ_l，M₂₃＝ρ_s，M₂₄＝0；M₃₂＝0，M₃₃＝ρ_s，M₃₄＝0；

M₄₂＝0，M₄₃＝0，M₄₄＝ρ_l-ρ_s；

N₁＝Q-q_f(h_l-h_f)-q_s(h_s-h_l)，

N₄＝q_f-q_s；o₁＝p，o₂＝V_l，o₃＝V_sd，o₄＝V_lt。

最终可得：

$\left[\begin{array}{c}{\mathrm{do}}_1/dz\\ {\mathrm{do}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dz\end{array}\right]=M^{-1}N---\left(16\right)$

采用Runge-Kutta数值解法求解微分方程(11)和(16)。对仿真时间进行离散化处理，将[0，T₀]区间内的仿真时间，离散成n个时间间隔为ΔT的区间，离散后的区间表示为[0,ΔT]，[ΔT，2ΔT]，……，[(n-1)ΔT，nΔT]。ΔT为0.01s～0.1s。

3.4)汽包水位是由汽包中的饱和水和饱和水蒸气的体积共同决定的，忽略汽包体积和汽包水位之间的非线性，可得汽包水位的表达式：

$l_d=\frac{V_l+V_s^{+}}{A_d}-l_0---\left(17\right)$

式中，l_d表示机理模型计算得到的汽包水位，单位mm；l₀是汽包水位的基准水位线。

(4)建立ANFIS补偿模型，步骤如下：

4.1)初始化粒子群。以聚类半径r_α作为粒子，15个粒子作为一个种群，每个粒子随机赋予[0.20.9]区间内的随机值，其中第i个粒子的位置的向量标示为ri，i＝1，2，…，15；

4.2)以r_i为聚类半径，进行减法聚类分析，将每个数据点作为可能的聚类中心，并根据各个数据点周围的数据点密度来计算该点作为聚类中心的可能性。

每个数据点X_i作为聚类中心的可能性P_i由式(18)来定义：

$P_i=\underset{j=1}{\overset{m}{\Σ}}\exp (-\frac{\left|\right|X_i-X_j|{|}^2}{{(r_i/2)}^2})---\left(18\right)$

式中m表示n维输入空间中全部的数据点数，X_i＝[X_i1,X_i2,...,X_in]、X_j＝[X_j1,X_j2,...,X_jn]是具体的数据点，r_i是一个正数，定义了该点的邻域半径，||·||符号表示欧式距离。被选为聚类中心的点具有最高的数据点密度，同时该该数据点周围的点被排除作为聚类中心的可能性。第一个聚类中心为X_C1，数据点密度为P_c1。选出第一个聚类中心后，继续采用类似的方法确定下一个聚类中心，但需消除已有聚类中心的影响，修改密度指标的山峰函数如下：

$P_i=P_i-P_{c1}\exp (-\frac{\left|\right|X_i-X_{c1}|{|}^2}{{(r_{\mathrm{\β}}/2)}^2})---\left(19\right)$

其中，r_β＝1.5r_i。循环重复上述过程直到所有剩余数据点作为聚类中心的可能性低于某一阈值δ，即P_ck/P_c1＜δ。

4.3)ANFIS模型训练。根据减法聚类算法得到的聚类中心，按照ANFIS模型结构训练汽包水位补偿模型；对于模糊神经网络模型的所有参数，采用混合最小二乘法的梯度下降算法进行学习。

4.4)计算适应度值。利用训练得到的预测模型计算汽包水位补偿值l^*，并与实际补偿值l_e进行比较，并以误差平方和作为粒子的适应度值MSE，适应度计算公式如下：

$MSE=\underset{i=1}{\overset{m}{\Σ}}{(l_e-l^{*})}^2---\left(28\right)$

4.5)更新极值。以适应度值为评价指标，比较当代粒子与上一代粒子之间的适应度值大小，如果当前粒子的适应度值优于上一代，则将当前粒子的位置设置为个体极值，否则个体极值保持不变。同时获取当代所有粒子适应度值最优的粒子，并与上一代最优粒子进行比较，如果当代最优粒子的适应度值优于上一代最优粒子的适应度值，则将当代粒子的最优适应度值设置为全局最优值，否则全局最优值保持不变。

4.6)更新粒子。根据最新的个体极值和全局极值，按照(29)式和(30)式更新粒子的速度v_id(t)和位置x_id(t)。

v_id(t+1)＝ωv_id(t)+c₁r₁(p_id-x_id(t))+c₂r₂(p_gd-x_id(t)) (29)

x_id(t+1)＝x_id(t)+v_id(t+1) (30)

t是粒子群优化算法的寻优代数。更进一步，为了改善基本粒子群算法容易陷入局部极值和收敛速度慢的缺陷，在PSO算法的基础上引进了动态加速常数c₁、c₂和惯性权重ω：

$\mathrm{\ω}={\mathrm{\ω}}_{max}-\frac{{\mathrm{\ω}}_{max}-{\mathrm{\ω}}_{min}}{T_{\max }}t---\left(31\right)$

$c_1=R_1+\frac{R_2\×t}{T_{max}}---\left(32\right)$

$c_2=R_3-\frac{R_4\×t}{T_{max}}---\left(33\right)$

其中，T_max为最大寻优代数，ω_max为最大惯性权重，ω_min为最小惯性权重，R₁、R₂、R₃、R₄为常数。

4.7)算法停止条件算法判定。判断是否达到最大迭代次数或者到达预测精度的要求，如果没有达到则返回步骤2.b)，利用更新的聚类半径继续搜索，否则退出搜索。

4.8)利用最终寻优得到的聚类半径，对样本进行聚类分析和ANFIS模型训练，得到达到训练终止条件的ANFIS模型，即汽包水位的补偿模型。

5)汽包水位的预测。选择需要预测的工况，通过机理模型得到汽包水位的预测值l_d，通过ANFIS模型得到汽包水位的补偿值l^*，则混合模型的输出值即最终预测的汽包水位为

本发明的有益效果主要表现在：本发明针对采用自然循环方式的CFB生活垃圾焚烧锅炉，先从汽包水位变化特性的机理出发，在合理简化假设的基础上，依据质量守恒方程、能量守恒方程和一些基础的方程，并结合CFB生活垃圾焚烧锅炉水冷壁独特的热量吸收分布规律，建立起适度复杂的基于微分方程形式的汽包水位动态特性机理模型。然后利用ANFIS建模来挖掘运行历史数据中隐含的知识，以补偿汽包水位的机理建模过程当中，由于模型简化假设、模型降阶或线性化处理、对事物变化过程机理的认知不完备、对象特性不同以及内外扰动所带来的误差。充分发挥机理建模和ANFIS建模的优势，提升汽包水位的预测精度。所建立的模型能够用于指导CFB生活垃圾焚烧锅炉的设计运行，也能为基于模型的控制算法提供可靠的过程模型。

附图说明

图1是本发明所提出的系统的结构图。

图2是本发明所提出的上位机系统的结构图。

图3本发明所提出的混合建模结构的示意图。

图4是本发明所采用ANFIS模型的系统结构图。

图5是本发明所提出的智能建模方法的流程图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

实施例1

如图1-5所示，本发明提供的一种循环流化床生活垃圾焚烧锅炉汽包水位的预测系统。该系统与循环流化床锅炉的集散控制系统相连，包括数据通讯接口和上位机，所述上位机包括：

第一信号采集模块。在模型输入变量中，锅炉的运行参数是需要从DCS系统中采集的。利用第一信号采集模型采集CFB生活垃圾焚烧锅炉在正常运行时的运行工况状态参数和操作变量，并组成ANIFS汽包水位补偿模型输入变量的训练样本矩阵X(m×n)，m表示样本个数，n表示变量的个数，其中给水流量、给水温度、主汽温度和主汽压力这几个参数需要同时提供给机理建模使用。

数据预处理模块。对X(m×n)进行粗大误差处理和随机误差处理，以摒除那些并不是反映锅炉正常运行工况的虚假信息，将锅炉停炉、压火、给料机堵塞等异常工况排除掉，为了避免预测模型的参数之间量纲和数量级的不同对模型性能造成的不良影响，训练样本输入变量均经过归一化处理后映射到[0，1]区间内，得到标准化后的训练样本X^*(m×n)。预处理过程采用以下步骤进行：

1.1)根据拉伊达准则，剔除训练样本X(m×n)中的野值；

1.2)剔除锅炉停炉运行工况，锅炉停炉时炉膛给煤机和给料机的开度为零，并且炉膛中温度接近常温；

1.3)剔除炉膛压火运行状况，锅炉压火时一次风机、二次风机引风机炉膛给煤机和给料机的开度为零，但是炉膛密相区的温度维持在350℃～450℃；

1.4)剔除给料机堵塞工况，给料机堵塞需要运行人员通过给料口的摄像头拍摄的画面对给料情况进行判断，给料机堵塞时，运行人员会显著地调高给料机的开度，反映在运行数据上，即给料机的开度大于35％；

1.5)数据归一化处理。按照式(1)将数据变量映射到[01]的区间内。

(1)

式中表示第J变量所组成的向量，min()表示最小值，max()表示最大值。

专家知识库模块。将水位计测量出来的实际汽包水位值l(m×1)与机理模型计算得到的水位值l_d(m×1)相减，将两者之间的差值作为ANFIS汽包水位补偿模型训练样本的输出向量l^*(m×1)，X^*和l^*共同组成ANFIS模型的训练样本，并进行保存。

混合建模模块。混合建模模块是CFB生活垃圾焚烧锅炉汽包水位预测系统的核心部分，主要包括机理建模和ANFIS建模两个部分。首先从汽包水位变化特性的机理出发，在合理简化假设的基础上，依据质量守恒方程、能量守恒方程和一些基础的方程，建立起适度复杂的基于微分方程形式的汽包水位动态特性模型即机理模型。然后计算不同工况下汽包水位的预测值，并与锅炉的实际汽包水位值进行比较，将差值作为ANFIS模型训练样本的输出向量，利用ANFIS模型的非线性动力学特性、自学习能力和泛化能力去学习并表征表达机理模型和实际物理过程之间的误差，以弥补汽包水位机理模型的不足之处。

机理建模的步骤如下：

2.1)对水冷壁中的工质进行质量守恒、动量守恒和能量守恒分析，可得：

$\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v\]=0---\left(2\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s{v}^2+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l{v}^2\]+\frac{{\mathrm{dp}}_r}{dz}\\ =-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]g-{\mathrm{\τ}}_{wm}\end{array}---\left(3\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v(h_s+\frac{1}{2}{v}^2)\]+\frac{d}{dz}\[(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v(h_l+\frac{1}{2}{v}^2)\]\\ ={\stackrel{\·}{q}}_{wm}-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]vg\end{array}---\left(4\right)$

其中，z为坐标轴，正向为工质的上升方向；ρ_s、ρ_l分别表示研究对象内饱和水蒸气的密度和饱和水的密度，密度的单位是kg/m³；α_s表示截面容积含汽率；v表示工质的流速，单位是m/s；τ_wm是壁面剪切应力，单位是Pa/m³；p_r是上升管中的压力，单位Pa；表示热流密度，单位是W/m³；g是重力加速度，单位是m/s²。h_l为饱和水的比焓。

饱和温度和工质压力之间存在着以下关系：

T＝T_b(1+a₁ln(P/P_b)+a₂ln²(P/P_b)+a₃ln³(P/P_b)+a₄ln⁴(P/P_b)) (5)

其中P是当前位置压力，T是当前位置的温度，T_b参考点的温度，P_b表示参考点的压力，a₁、a₂、a₃、a₄是常系数。

壁面的剪切应力可以通过下式求取：

${\mathrm{\τ}}_{wm}=\frac{2f}{d_t}{\mathrm{\ρ}}_m{v}^2---\left(6\right)$

f＝0.0014+0.125Re^-0.32 (7)

ρ_m＝α_sρ_s+(1-α_s)ρ_l (8)

式中Re表示雷诺数，f为摩擦系数，d_t为水冷壁上升管的管径，ρ_m为汽水混合物的密度。

水冷壁的热流密度可以通过下式求取：

${\stackrel{\·}{q}}_{wm}=\left(\frac{{\mathrm{\πd}}_t}{A_t}\right)U_{wm}(T_w-T)---\left(9\right)$

式中，A_t表示截面积，U_wm表示传热系数，T_w表示管壁的温度。

对(2)～(4)式进行展开，并将(5)～(9)式代入，可得到如下方程组：

$\left[\begin{array}{ccccc}{A}_{11}& A_{12}& ...& ...& A_{15}\\ A_{21}& A_{22}& ...& ...& A_{25}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ A_{51}& A_{52}& ...& ...& A_{55}\end{array}\right]\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=\left[\begin{array}{c}{B}_1\\ B_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ B_5\end{array}\right]---\left(10\right)$

其中，A₁₁＝(ρ_s-ρ_l)v，A₁₂＝α_sv，A₁₃＝α_sρ_s+(1-α_s)ρ_l，A₁₄＝0，A₁₅＝0；

$A_{21}=\{{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v-\{{\mathrm{\ρ}}_l\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v,$

$A_{22}={\mathrm{\α}}_{s}v\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\],$

$A_{23}={\mathrm{\α}}_s{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{3}{2}{v}^2\]+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l\[c_p^l(T-T_0)+\frac{3}{2}{v}^2\],$

A₂₅＝0；A₃₁＝(ρ_s-ρ_l)v²，A₃₂＝α_sv²，

A₃₃＝2[α_sρ_s+(1-α_s)ρ_l]v，A₃₄＝0，A₃₅＝1；A₄₁＝0，A₄₂＝0，A₄₃＝0，A₄₄＝1，

A₅₁＝0，A₅₂＝1，A₅₃＝0，

y₁＝α_s，y₂＝ρ_s，y₃＝v，y₄＝T，y₅＝P；B₁＝0，

$B_2=\left(\frac{4}{d_t}\right)U_{wm}(T_w-T)-({\mathrm{\α}}_s{\mathrm{\ρ}}_s+(1-{\mathrm{\α}}_s\left){\mathrm{\ρ}}_l\right)vg,$

B₄＝0，B₅＝0。式中，λ_V表示汽化潜热，表示饱和水的定压比热，表示饱和蒸汽的定压比热，T₀表示指定状态下工质的温度。最终可以算得：

$\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=A^{-1}B---\left(11\right)$

2.2)通过对汽包-下降管-水冷壁中的所有工质进行质量守恒分析和能量守恒分析，可得：

$k_1\frac{dp}{dt}=Q-q_f(h_l-h_f)-q_s(h_s-h_l)---\left(12\right)$

$e_1\frac{{\mathrm{dV}}_{lt}}{dt}+e_2\frac{dp}{dt}=q_f-q_s---\left(13\right)$

其中，e₁＝ρ_l-ρ_s；

式中，p表示汽包压力，单位Pa；Q表示整个蒸发区的单位时间内的吸热量，单位是W；V_lt表示汽包-下降管-水冷壁中所有饱和水的体积，单位是m³；V_st表示回路中所有水蒸气的体积，单位是m³；q_f、q_s分别表示给水流量和主汽流量和下降管中循环工质的流量，流量的单位是kg/s；h_f、h_l、h_s分别表示给水比焓、饱和水的比焓和饱和水蒸气的比焓，单位是kJ/kg；m_t表示水冷壁金属和汽包金属的总质量，单位是kg；C_p表示金属的比热容；t_s表示饱和温度。

2.3)对汽包中液面下的工质进行质能守恒分析，忽略液面上饱和蒸汽的凝结，可得：

$\frac{d\left({\mathrm{\ρ}}_s{V}_s^{+}\right)}{dt}={\mathrm{\α}}_r{q}_r-q_s^{+}-q_{cd}---\left(14\right)$

$\frac{{\mathrm{dM}}_{dt}^{+}}{dt}=q_f+q_r-q_s^{+}-q_{dc}---\left(15\right)$

其中，

$q_s^{+}=A_d\left(\frac{V_s^{+}}{V_l+V_s^{+}}\right)v_s^{+}{\mathrm{\ρ}}_s;v_s^{+}=1.41{\[\frac{\mathrm{\σ}g({\mathrm{\ρ}}_l-{\mathrm{\ρ}}_s)}{{\mathrm{\ρ}}_l^2}\]}^{1/4};M_{dt}^{+}=V_l{\mathrm{\ρ}}_l+V_s^{+}{\mathrm{\ρ}}_s.$

式中，α_r表示上升管出口的截面质量含汽率；表示汽包内液面下工质的质量，单位是kg；q_r、q_cd、q_dc分别表示上升管中的流量、饱和蒸汽逸出液面的流量、饱和蒸汽凝结的流量和下降管中的流量，单位是kg/s；t_s表示饱和蒸汽的温度，单位是℃；m_d表示汽包金属的质量，单位kg；C_p表示水冷壁和汽包中金属的比热容，单位是kJ/(kg·℃)；A_d是汽包液位的截面面积，单位是m²；是饱和蒸汽逸出的速度，单位是m/s；σ是液面的表面张力，单位是N/m；V_l分别表示液面下饱和水蒸气的体积和饱和水的体积，单位是m³。

联立(12)～(15)式并进行展开，可以得到如下方程组：

$\left[\begin{array}{ccccc}{M}_{11}& M_{12}& ...& ...& M_{14}\\ M_{21}& M_{22}& ...& ...& M_{24}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ M_{41}& M_{42}& ...& ...& M_{44}\end{array}\right]\left[\begin{array}{c}{\mathrm{do}}_1/dt\\ {\mathrm{do}}_2/dt\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dt\end{array}\right]=\left[\begin{array}{c}{N}_1\\ N_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ N_4\end{array}\right]$

其中，M₁₂＝0，M₁₃＝0，M₁₄＝0；M₂₂＝ρ_l，M₂₃＝ρ_s，M₂₄＝0；M₃₂＝0，M₃₃＝ρ_s，M₃₄＝0；

M₄₂＝0，M₄₃＝0，M₄₄＝ρ_l-ρ_s；

N₁＝Q-q_f(h_l-h_f)-q_s(h_s-h_l)，

N₄＝q_f-q_s；o₁＝p，o₂＝V_l，o₃＝V_sd，o₄＝V_lt。

最终可得：

$\left[\begin{array}{c}{\mathrm{do}}_1/dz\\ {\mathrm{do}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dz\end{array}\right]=M^{-1}N---\left(16\right)$

汽包水位是由汽包中的饱和水和饱和水蒸气的体积共同决定的，忽略汽包体积和汽包水位之间的非线性，可得汽包水位的表达式：

$l_d=\frac{V_l+V_s^{+}}{A_d}-l_0---\left(17\right)$

式中，l_d表示机理模型计算得到的汽包水位，单位mm；l₀是汽包水位的基准水位线。

2.4)模型参数设置。模型中关于汽包、水冷壁结构的参数可以通过查阅锅炉的设计资料获得，饱和水与饱和水蒸气的物性参数可以通过IAPWS-IF97公式来计算，给水流量、水蒸气流量等实时变化参数与DCS系统中的显示的参数同步变化。

2.5)通过求解(11)式可以构建水冷壁中基于空间维度的水蒸气特性模型(称之为水冷壁模型)，求得水冷壁中饱和蒸汽的截面容积含汽率、工质的流速和压力，但该模型是静态的。通过求解(16)和(17)式，可以得到基于时间维度的汽包压力和汽包水位的动态表达式(称之为汽包模型)。汽包模型和水冷壁模型有机结合在一起，并相互提供模型信息支持，采用Runge-Kutta数值解法求解微分方程。为了获取水冷壁模型的动态特性，需要对仿真时间进行离散化处理，将[0，T₀]区间内的仿真时间，离散成n个时间间隔为ΔT的区间，离散后的区间表示为[0,ΔT]，[ΔT，2ΔT]，……，[(n-1)ΔT，nΔT]。ΔT的大小对模型的性能有重要的影响，ΔT太大则会导致仿真的精度较低，甚至出现不收敛的情况；ΔT太小则会耗费过多的计算资源，导致仿真时间过长。通常将ΔT选择在0.01s～0.1s之间较为合适。

在获得汽包水位的机理预测模型之后，利用汽包水位的机理模型，计算不同工况下的汽包水位值，并与实际的汽包水为值进行比较，得到两者之间的误差l_e，作为ANFIS模型的输出训练向量。在建立ANFIS汽包水位补偿模型时，先利用减法聚类算法(SubtractiveClustering,SC)对样本数据进行特征提取，自适应的确定初始模糊规则和模糊神经网络的初始结构参数，再结合最小二乘估计法和误差反向传播算法对模糊神经网络的参数进行学习训练。在该模型中聚类半径是影响建模性能的关键参数，因此以预测精度为目标，利用粒子群优化算法(Particle Swarm Optimization,PSO)寻找聚类半径的最优值。算法步骤如下：

2.a)初始化粒子群。以聚类半径r_α作为粒子，15个粒子作为一个种群，每个粒子随机赋予[0.20.9]区间内的随机值，其中第i个粒子的位置的向量标示为ri，i＝1，2，…，15；

2.b)以r_i为聚类半径，进行减法聚类分析。减法聚类算法用于对建模数据样本的空间进行初始划分以及模糊规则的确定，K-均值聚类算法和模糊C-均值聚类算法均需预设聚类中心的数目，没有充分利用样本空间的蕴含的对象特征信息。而减法聚类算法是一种基于山峰函数的聚类算法，它将每个数据点作为可能的聚类中心，并根据各个数据点周围的数据点密度来计算该点作为聚类中心的可能性。

每个数据点X_i作为聚类中心的可能性P_i由式(18)来定义：

$P_i=\underset{j=1}{\overset{m}{\Σ}}\exp (-\frac{\left|\right|X_i-X_j|{|}^2}{{(r_i/2)}^2})---\left(18\right)$

式中m表示n维输入空间中全部的数据点数，X_i＝[X_i1,X_i2,...,X_in]、X_j＝[X_j1,X_j2,...,X_jn]是具体的数据点，r_i是一个正数，定义了该点的邻域半径，||·||符号表示欧式距离。被选为聚类中心的点具有最高的数据点密度，同时该该数据点周围的点被排除作为聚类中心的可能性。第一个聚类中心为X_C1，数据点密度为P_c1。选出第一个聚类中心后，继续采用类似的方法确定下一个聚类中心，但需消除已有聚类中心的影响，修改密度指标的山峰函数如下：

$P_i=P_i-P_{c1}\exp (-\frac{\left|\right|X_i-X_{c1}|{|}^2}{{(r_{\mathrm{\β}}/2)}^2})---\left(19\right)$

其中，r_β定义了一个密度指标显著减小的邻域，为了避免出现十分接近的聚类中心，r_β＝1.5r_i。循环重复上述过程直到所有剩余数据点作为聚类中心的可能性低于某一阈值δ，即P_ck/P_c1＜δ。

2.c)ANFIS模型训练。不失一般性，假定减法聚类算法得到两个聚类中心(X₁,Y₁)和(X₂,Y₂)，得到两条模糊规则：

Rule 1:IFisandisand…andisTHEN

Rule 2:IFisandisand…andisTHEN

ANFIS系统的第一层为系统的输入层，由n个节点组成，它的作用是将输入向量j＝1,2,…,m按原值传递到下一层。

第二层为模糊化层，由2n个节点组成，它的作用是计算各输入分量属于各语言变量值模糊集合的隶属函数该层采用高斯函数进行模糊化处理，每个节点的输出：

${\mathrm{\μ}}_i^j=\exp (-\frac{(x_i-c_{ij})}{{\mathrm{\σ}}_{ij}}),i=1,2,...,2n;j=1,2,...,2n---\left(20\right)$

式中，c_ij和σ_ij分别表示隶属函数的中心和宽度。

第三层共有两个节点，每个节点代表一条模糊规则，他的作用是完成模糊逻辑的“与”操作，采用product推理计算出每条规则的适应度值，即

$\begin{array}{c}{w}_1={\mathrm{\μ}}_1^1{\mathrm{\μ}}_2^1...{\mathrm{\μ}}_{2n}^1\\ w_2={\mathrm{\μ}}_1^2{\mathrm{\μ}}_2^2...{\mathrm{\μ}}_{2n}^2\end{array}---\left(21\right)$

第四层为归一化层，节点数与第三层一样，它所实现的是归一化计算，即

${\stackrel{\‾}{w}}_j=\frac{w_j}{w_1+w_2},j=1,2---\left(22\right)$

第五层为结论层，该层与第三层的节点数相同。其节点输出为

$O_{4i}=\stackrel{\‾}{w_i}{f}_i=\stackrel{\‾}{w_i}(r_i+p_i^1{X}^1+p_i^2{X}^1+...+p_i^n{X}^n),i=1,2---\left(23\right)$

第六层为输出层及去模糊化层，只有一个节点，使用面积中心法进行解模糊，得到网络的最终输出为

$O_5=\frac{w_1{f}_1+w_2{f}_2}{w_1+w_2}---\left(24\right)$

对于模糊神经网络模型的所有参数，采用混合最小二乘法的梯度下降算法进行学习，步骤如下：

2.c.1)在固定高斯型隶属函数的中心和宽度的前提下，利用最小二乘法计算线性结论参数{r,p}；

2.c.2)固定结论参数，采用误差反向传播算法对高斯函数的中心和宽度进行学习可得：

$S(t+1)=S\left(t\right)-{\mathrm{\η}}_s\left(t\right)\frac{\∂E}{\∂S}+\mathrm{\α}\mathrm{\Δ}S\left(t\right)---\left(25\right)$

ΔS(t)＝S(t+1)-S(t) (26)

$E=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{(f-f_{*})}^2---\left(27\right)$

式中S为参数c和σ，η_s为学习率，α为动量项，f为预测输出，f_*为实际输出，E为平方误差和，n为迭代步数。利用训练样本数据，重复上述步骤，直到满足误差指标或者达到最大训练次数。

2.d)计算适应度值。利用训练得到的预测模型计算汽包水位补偿值l^*，并与实际补偿值l_e进行比较，并以误差平方和作为粒子的适应度值MSE，适应度计算公式如下：

$MSE=\underset{i=1}{\overset{m}{\Σ}}{(l_e-l^{*})}^2---\left(28\right)$

2.e)更新极值。以适应度值为评价指标，比较当代粒子与上一代粒子之间的适应度值大小，如果当前粒子的适应度值优于上一代，则将当前粒子的位置设置为个体极值，否则个体极值保持不变。同时获取当代所有粒子适应度值最优的粒子，并与上一代最优粒子进行比较，如果当代最优粒子的适应度值优于上一代最优粒子的适应度值，则将当代粒子的最优适应度值设置为全局最优值，否则全局最优值保持不变。

2.f)更新粒子。根据最新的个体极值和全局极值，按照(29)式和(30)式更新粒子的速度v_id(t)和位置x_id(t)。

v_id(t+1)＝ωv_id(t)+c₁r₁(p_id-x_id(t))+c₂r₂(p_gd-x_id(t)) (29)

x_id(t+1)＝x_id(t)+v_id(t+1) (30)

t是粒子群优化算法的寻优代数。更进一步，为了改善基本粒子群算法容易陷入局部极值和收敛速度慢的缺陷，在PSO算法的基础上引进了动态加速常数c₁、c₂和惯性权重ω：

$\mathrm{\ω}={\mathrm{\ω}}_{max}-\frac{{\mathrm{\ω}}_{max}-{\mathrm{\ω}}_{min}}{T_{\max }}t---\left(31\right)$

$c_1=R_1+\frac{R_2\×t}{T_{max}}---\left(32\right)$

$c_2=R_3-\frac{R_4\×t}{T_{max}}---\left(33\right)$

其中，T_max为最大寻优代数，ω_max为最大惯性权重，ω_min为最小惯性权重，R₁、R₂、R₃、R₄为常数。

2.g)算法停止条件算法判定。判断是否达到最大迭代次数或者到达预测精度的要求，如果没有达到则返回步骤2.b)，利用更新的聚类半径继续搜索，否则退出搜索。

2.h)利用最终寻优得到的聚类半径，对样本进行聚类分析和ANFIS模型训练，得到达到训练终止条件的ANFIS模型，即汽包水位的补偿模型。

第二信号采集模块。用于从数据库中选择需要预测汽包水位的运行工况，或者实时地采集当前锅炉的运行工况。

预测模块。该模块用于对指定的样本进行汽包水位的预测，或者对当前锅炉运行工况下的汽包水位进行实时预测。

结果显示模块。显示汽包水位的预测结果，或者对汽包水位的预测结果进行统计分析。

实施例2

参照图1、图2、图3、图4，本发明提供的一种循环流化床生活垃圾焚烧锅炉汽包水位预测方法，该方法包括以下步骤：

1)选择混合模型的输入变量并采集训练样本。锅炉汽包水位主要是由工质的状态和汽包中的物质能量平衡所决定的，影响汽包水位变化的因素很多，主要包括锅炉的设计参数、工质的物性参数和运行参数。设计参数主要有水冷壁的当量管径、水冷壁的有效吸热长度、水冷壁的容积、循环水的流量、汽包的容积、汽包的平均截面面积、下降管的容积、汽包-下降管-水冷壁三者总共的金属质量。工质的状态可以通过工质的物性参数来表征，包括压力、饱和温度、比焓、比容和汽化潜热等参数，这些参数可以通过IAPWS-IF97公式来计算。锅炉的设计参数在没有发生锅炉改造的情况下是不变的，而工质的物性参数则可以通过显式的公式进行计算，这两种类型的参数主要用于汽包水位机理建模使用。

汽包中物质平衡通过给水流量和主蒸汽流量之间的差值来决定的，因此给水流量和主蒸汽流量需要作为模型的输入变量；能量平衡主要通过主蒸汽焓值和燃料热值和之间的差值所决定的，主汽焓值由主汽流量和蒸汽比焓所决定的，蒸汽比焓可以通过蒸汽热力参数的计算程序获得。燃料热值的获取比较困难，因为生活垃圾的组成成分十分复杂，并且波动十分频繁，目前尚无可靠的垃圾热值在线测量的仪器，但是可以通过炉膛中温度场的变化情况和排烟温度变化情况间接反映出垃圾热值的变化情况，因此将垃圾给料量、辅助燃煤量、一次风量、二次风量以及炉膛中所有的温度测点作为模型的输入变量。另外一方面，垃圾在循环流化床中的燃烧是一个十分复杂的剧烈物理化学变化过程，垃圾在进入炉膛之后会经历几个过程：干燥加热、挥发分析出及燃烧、焦炭燃烧。垃圾中质轻易碎的组分如纸纸张、塑料以及细颗粒等会在流化风的作用下进入炉膛上部，经历干燥、挥发分的析出及燃烧以及残炭的燃烧等一系列过程；而密度较大、含水率高以及颗粒尺寸较大的组分如木头、厨余垃圾等终端速度大于流化速度的组分会落入密相区，并在密相区中被床料加热、燃烧。由于采取分级送风的燃烧方式，密相区处于缺氧燃烧的状态。与煤不同，垃圾中热量来源主要是挥发分，挥发分含量占有较高通常是固定碳含量的5～7倍，同时，垃圾中水分会降低挥发分的燃烧速率，大量的挥发分在炉膛的中部和上部燃烧，大量的热量在炉膛的中上部释放出来，而传统的燃煤CFB锅炉中燃料的热值集中在密相区和过渡区燃烧，大量的热量在炉膛的中下部释放。于此同时，为了保证床温的稳定性，提高床料的温度水平，CFB生活垃圾的密相区会浇注一层保温材料，受热面半绝热，进一步减少了水冷壁在密相区吸热的比例。所以，CFB生活垃圾焚烧锅炉沿水冷壁方向的热流密度分布规律与常规电站锅炉是不同的。而热流密度的分布与炉膛中温度场的分布有密切的联系，为了表征热流密度分布规律的影响，需要将炉膛中所有的温度测点作为模型的输入变量，不过，在表征燃料热值的影响时已经选择过了。这些运行参数主要用于创建ANFIS汽包水位补偿模型时使用。

在模型输入变量中，锅炉的运行参数是需要从DCS系统中采集的。采集CFB生活垃圾焚烧锅炉在正常运行时的运行工况状态参数和操作变量，并组成ANIFS汽包水位补偿模型输入变量的训练样本矩阵X(m×n)，m表示样本个数，n表示变量的个数，其中给水流量、给水温度、主汽温度和主汽压力这几个参数需要同时提供给机理建模使用。

2)训练样本预处理。对X(m×n)进行粗大误差处理和随机误差处理，以摒除那些并不是反映锅炉正常运行工况的虚假信息，将锅炉停炉、压火、给料机堵塞等异常工况排除掉，为了避免预测模型的参数之间量纲和数量级的不同对模型性能造成的不良影响，训练样本输入变量均经过归一化处理后映射到[0，1]区间内，得到标准化后的训练样本X^*(m×n)。

3)利用机理建模的方法构建CFB生活垃圾焚烧锅炉汽包水位的机理模型。

3.1)对水冷壁中的工质进行质量守恒、动量守恒和能量守恒分析，可得：

$\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v\]=0---\left(2\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s{v}^2+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l{v}^2\]+\frac{{\mathrm{dp}}_r}{dz}\\ =-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]g-{\mathrm{\τ}}_{wm}\end{array}---\left(3\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v(h_s+\frac{1}{2}{v}^2)\]+\frac{d}{dz}\[(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v(h_l+\frac{1}{2}{v}^2)\]\\ ={\stackrel{\·}{q}}_{wm}-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]vg\end{array}---\left(4\right)$

饱和温度和工质压力之间存在着以下关系：

T＝T_b(1+a₁ln(P/P_b)+a₂ln²(P/P_b)+a₃ln³(P/P_b)+a₄ln⁴(P/P_b)) (5)

其中P是当前位置压力，T是当前位置的温度，T_b参考点的温度，P_b表示参考点的压力，a₁、a₂、a₃、a₄是常系数。

壁面的剪切应力可以通过下式求取：

${\mathrm{\τ}}_{wm}=\frac{2f}{d_t}{\mathrm{\ρ}}_m{v}^2---\left(6\right)$

f＝0.0014+0.125Re^-0.32 (7)

ρ_m＝α_sρ_s+(1-α_s)ρ_l (8)

式中Re表示雷诺数，f为摩擦系数，d_t为水冷壁上升管的管径，ρ_m为汽水混合物的密度。

水冷壁的热流密度可以通过下式求取：

${\stackrel{\·}{q}}_{wm}=\left(\frac{{\mathrm{\πd}}_t}{A_t}\right)U_{wm}(T_w-T)---\left(9\right)$

式中，A_t表示截面积，U_wm表示传热系数，T_w表示管壁的温度。

对(2)～(4)式进行展开，并将(5)～(9)式代入，可得到如下方程组：

$\left[\begin{array}{ccccc}{A}_{11}& A_{12}& ...& ...& A_{15}\\ A_{21}& A_{22}& ...& ...& A_{25}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ A_{51}& A_{52}& ...& ...& A_{55}\end{array}\right]\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=\left[\begin{array}{c}{B}_1\\ B_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ B_5\end{array}\right]---\left(10\right)$

其中，A₁₁＝(ρ_s-ρ_l)v，A₁₂＝α_sv，A₁₃＝α_sρ_s+(1-α_s)ρ_l，A₁₄＝0，A₁₅＝0；

$A_{21}=\{{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v-\{{\mathrm{\ρ}}_l\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v,$

$A_{22}={\mathrm{\α}}_{s}v\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\],$

$A_{23}={\mathrm{\α}}_s{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{3}{2}{v}^2\]+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l\[c_p^l(T-T_0)+\frac{3}{2}{v}^2\],$

A₂₅＝0；A₃₁＝(ρ_s-ρ_l)v²，A₃₂＝α_sv²，

A₃₃＝2[α_sρ_s+(1-α_s)ρ_l]v，A₃₄＝0，A₃₅＝1；A₄₁＝0，A₄₂＝0，A₄₃＝0，A₄₄＝1，

A₅₁＝0，A₅₂＝1，A₅₃＝0，

y₁＝α_s，y₂＝ρ_s，y₃＝v，y₄＝T，y₅＝P；B₁＝0，

$B_2=\left(\frac{4}{d_t}\right)U_{wm}(T_w-T)-({\mathrm{\α}}_s{\mathrm{\ρ}}_s+(1-{\mathrm{\α}}_s\left){\mathrm{\ρ}}_l\right)vg,$

B₄＝0，B₅＝0。式中，λ_V表示汽化潜热，表示饱和水的定压比热，表示饱和蒸汽的定压比热，T₀表示指定状态下工质的温度。最终可以算得：

$\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=A^{-1}B---\left(11\right)$

3.2)通过对汽包-下降管-水冷壁中的所有工质进行质量守恒分析和能量守恒分析，可得：

$k_1\frac{dp}{dt}=Q-q_f(h_l-h_f)-q_s(h_s-h_l)---\left(12\right)$

$e_1\frac{{\mathrm{dV}}_{lt}}{dt}+e_2\frac{dp}{dt}=q_f-q_s---\left(13\right)$

其中，e₁＝ρ_l-ρ_s；

3.3)对汽包中液面下的工质进行质能守恒分析，忽略液面上饱和蒸汽的凝结，可得：

$\frac{d\left({\mathrm{\ρ}}_s{V}_s^{+}\right)}{dt}={\mathrm{\α}}_r{q}_r-q_s^{+}-q_{cd}---\left(14\right)$

$\frac{{\mathrm{dM}}_{dt}^{+}}{dt}=q_f+q_r-q_s^{+}-q_{dc}---\left(15\right)$

其中，

$q_s^{+}=A_d\left(\frac{V_s^{+}}{V_l+V_s^{+}}\right)v_s^{+}{\mathrm{\ρ}}_s;v_s^{+}=1.41{\[\frac{\mathrm{\σ}g({\mathrm{\ρ}}_l-{\mathrm{\ρ}}_s)}{{\mathrm{\ρ}}_l^2}\]}^{1/4};M_{dt}^{+}=V_l{\mathrm{\ρ}}_l+V_s^{+}{\mathrm{\ρ}}_s.$

联立(12)～(15)式并进行展开，可以得到如下方程组：

$\left[\begin{array}{ccccc}{M}_{11}& M_{12}& ...& ...& M_{14}\\ M_{21}& M_{22}& ...& ...& M_{24}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ M_{41}& M_{42}& ...& ...& M_{44}\end{array}\right]\left[\begin{array}{c}{\mathrm{do}}_1/dt\\ {\mathrm{do}}_2/dt\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dt\end{array}\right]=\left[\begin{array}{c}{N}_1\\ N_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ N_4\end{array}\right]$

其中，M₁₂＝0，M₁₃＝0，M₁₄＝0；M₂₂＝ρ_l，M₂₃＝ρ_s，M₂₄＝0；M₃₂＝0，M₃₃＝ρ_s，M₃₄＝0；

M₄₂＝0，M₄₃＝0，M₄₄＝ρ_l-ρ_s；

N₁＝Q-q_f(h_l-h_f)-q_s(h_s-h_l)，

N₄＝q_f-q_s；o₁＝p，o₂＝V_l，o₃＝V_sd，o₄＝V_lt。

最终可得：

$\left[\begin{array}{c}{\mathrm{do}}_1/dz\\ {\mathrm{do}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dz\end{array}\right]=M^{-1}N---\left(16\right)$

汽包水位是由汽包中的饱和水和饱和水蒸气的体积共同决定的，忽略汽包体积和汽包水位之间的非线性，可得汽包水位的表达式：

$l_d=\frac{V_l+V_s^{+}}{A_d}-l_0---\left(17\right)$

3.4)模型参数设置。模型中关于汽包、水冷壁结构的参数可以通过查阅锅炉的设计资料获得，饱和水与饱和水蒸气的物性参数可以通过IAPWS-IF97公式来计算，给水流量、水蒸气流量等实时变化参数与DCS系统中的显示的参数同步变化。

3.5)通过求解(11)式可以构建水冷壁中基于空间维度的水蒸气特性模型(称之为水冷壁模型)，求得水冷壁中饱和蒸汽的截面容积含汽率、工质的流速和压力，但该模型是静态的。通过求解(16)和(17)式，可以得到基于时间维度的汽包压力和汽包水位的动态表达式(称之为汽包模型)。汽包模型和水冷壁模型有机结合在一起，并相互提供模型信息支持，采用Runge-Kutta数值解法求解微分方程。为了获取水冷壁模型的动态特性，需要对仿真时间进行离散化处理，将[0，T₀]区间内的仿真时间，离散成n个时间间隔为ΔT的区间，离散后的区间表示为[0,ΔT]，[ΔT，2ΔT]，……，[(n-1)ΔT，nΔT]。ΔT的大小对模型的性能有重要的影响，ΔT太大则会导致仿真的精度较低，甚至出现不收敛的情况；ΔT太小则会耗费过多的计算资源，导致仿真时间过长。通常将ΔT选择在0.01s～0.1s之间较为合适。

4)建立ANFIS汽包水位补偿模型。在获得汽包水位的机理预测模型之后，利用汽包水位的机理模型，计算不同工况下的汽包水位值，并与实际的汽包水为值进行比较，得到两者之间的误差l_e，作为ANFIS模型的输出训练向量。在建立ANFIS汽包水位补偿模型时，先利用减法聚类算法(Subtractive Clustering,SC)对样本数据进行特征提取，自适应的确定初始模糊规则和模糊神经网络的初始结构参数，再结合最小二乘估计法和误差反向传播算法对模糊神经网络的参数进行学习训练。在该模型中聚类半径是影响建模性能的关键参数，因此以预测精度为目标，利用粒子群优化算法(Particle SwarmOptimization,PSO)寻找聚类半径的最优值。算法步骤如下：

4.1)初始化粒子群。以聚类半径r_α作为粒子，15个粒子作为一个种群，每个粒子随机赋予[0.20.9]区间内的随机值，其中第i个粒子的位置的向量标示为r_i，i＝1，2，…，15；

4.2)以r_i为聚类半径，进行减法聚类分析。减法聚类算法用于对建模数据样本的空间进行初始划分以及模糊规则的确定，K-均值聚类算法和模糊C-均值聚类算法均需预设聚类中心的数目，没有充分利用样本空间的蕴含的对象特征信息。而减法聚类算法是一种基于山峰函数的聚类算法，它将每个数据点作为可能的聚类中心，并根据各个数据点周围的数据点密度来计算该点作为聚类中心的可能性。

每个数据点X_i作为聚类中心的可能性P_i由式(18)来定义：

$P_i=\underset{j=1}{\overset{m}{\Σ}}\exp (-\frac{\left|\right|X_i-X_j|{|}^2}{{(r_i/2)}^2})---\left(18\right)$

式中m表示n维输入空间中全部的数据点数，X_i＝[X_i1,X_i2,...,X_in]、X_j＝[X_j1,X_j2,...,X_jn]是具体的数据点，r_i是一个正数，定义了该点的邻域半径，||·||符号表示欧式距离。被选为聚类中心的点具有最高的数据点密度，同时该该数据点周围的点被排除作为聚类中心的可能性。第一个聚类中心为X_C1，数据点密度为P_c1。选出第一个聚类中心后，继续采用类似的方法确定下一个聚类中心，但需消除已有聚类中心的影响，修改密度指标的山峰函数如下：

$P_i=P_i-P_{c1}\exp (-\frac{\left|\right|X_i-X_{c1}|{|}^2}{{(r_{\mathrm{\β}}/2)}^2})---\left(19\right)$

其中，r_β定义了一个密度指标显著减小的邻域，为了避免出现十分接近的聚类中心，r_β＝1.5r_i。循环重复上述过程直到所有剩余数据点作为聚类中心的可能性低于某一阈值δ，即P_ck/P_c1＜δ。

4.3)ANFIS模型训练。根据减法聚类算法得到的聚类中心，按照ANFIS模型结构训练汽包水位补偿模型；对于模糊神经网络模型的所有参数，采用混合最小二乘法的梯度下降算法进行学习。

4.4)计算适应度值。利用训练得到的预测模型计算汽包水位补偿值l^*，并与实际补偿值l_e进行比较，并以误差平方和作为粒子的适应度值MSE，适应度计算公式如下：

$MSE=\underset{i=1}{\overset{m}{\Σ}}{(l_e-l^{*})}^2---\left(28\right)$

4.5)更新极值。以适应度值为评价指标，比较当代粒子与上一代粒子之间的适应度值大小，如果当前粒子的适应度值优于上一代，则将当前粒子的位置设置为个体极值，否则个体极值保持不变。同时获取当代所有粒子适应度值最优的粒子，并与上一代最优粒子进行比较，如果当代最优粒子的适应度值优于上一代最优粒子的适应度值，则将当代粒子的最优适应度值设置为全局最优值，否则全局最优值保持不变。

4.6)更新粒子。根据最新的个体极值和全局极值，按照(29)式和(30)式更新粒子的速度v_id(t)和位置x_id(t)。

v_id(t+1)＝ωv_id(t)+c₁r₁(p_id-x_id(t))+c₂r₂(p_gd-x_id(t)) (29)

x_id(t+1)＝x_id(t)+v_id(t+1) (30)

t是粒子群优化算法的寻优代数。更进一步，为了改善基本粒子群算法容易陷入局部极值和收敛速度慢的缺陷，在PSO算法的基础上引进了动态加速常数c₁、c₂和惯性权重ω：

$\mathrm{\ω}={\mathrm{\ω}}_{max}-\frac{{\mathrm{\ω}}_{max}-{\mathrm{\ω}}_{min}}{T_{\max }}t---\left(31\right)$

$c_1=R_1+\frac{R_2\×t}{T_{max}}---\left(32\right)$

$c_2=R_3-\frac{R_4\×t}{T_{max}}---\left(33\right)$

其中，T_max为最大寻优代数，ω_max为最大惯性权重，ω_min为最小惯性权重，R₁、R₂、R₃、R₄为常数。

4.7)算法停止条件算法判定。判断是否达到最大迭代次数或者到达预测精度的要求，如果没有达到则返回步骤3.2)，利用更新的聚类半径继续搜索，否则退出搜索。

4.8)利用最终寻优得到的聚类半径，对样本进行聚类分析和ANFIS模型训练，得到达到训练终止条件的ANFIS模型，即汽包水位的补偿模型。

5)汽包水位的预测。选择需要预测的工况，通过机理模型得到汽包水位的预测值l_d，通过ANFIS模型得到汽包水位的补偿值l^*，则混合模型的输出值即最终预测的汽包水位为

本发明所提出的循环流化床锅汽包水位预测系统及方法，已通过上述具体实施步骤进行了描述，相关技术人员明显能在不脱离本发明内容、精神和范围内对本文所述的装置和操作方法进行改动或适当变更与组合，来实现本发明技术。特别需要指出的是，所有相类似的替换和改动对本领域的技术人员是显而易见的，它们都会被视为包括在本发明精神、范围和内容中。

Claims (2)

1.一种循环流化床生活垃圾焚烧锅炉汽包水位的预测系统，该系统与循环流化床锅炉的集散控制系统相连，包括数据通讯接口和上位机，所述上位机包括：

第一信号采集模块。利用该模型采集CFB生活垃圾焚烧锅炉在正常运行时的运行工况状态参数和操作变量，并组成ANIFS汽包水位补偿模型输入变量的训练样本矩阵X(m×n)，m表示样本个数，n表示变量的个数。

数据预处理模块。对X(m×n)进行粗大误差处理和随机误差处理，剔除训练样本中的野值，排除异常工况，将训练样本输入变量经归一化处理后映射到[0，1]区间内，得到标准化后的训练样本X^*(m×n)。

专家知识库模块。将水位计测量出来的实际汽包水位值l(m×1)与机理模型计算得到的水位值l_d(m×1)相减，将两者之间的差值作为ANFIS汽包水位补偿模型训练样本的输出向量l^*(m×1)，X^*和l^*共同组成ANFIS模型的训练样本，并进行保存。

混合建模模块。该模块主要包括机理建模和ANFIS建模两个部分。机理建模的步骤如下：

2.1)对水冷壁中的工质进行质量守恒、动量守恒和能量守恒分析，可得：

$\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v\]=0---\left(2\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s{v}^2+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l{v}^2\]+\frac{{\mathrm{dp}}_r}{dz}\\ =-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]g-{\mathrm{\τ}}_{wm}\end{array}---\left(3\right)$

$\begin{array}{c}\frac{d}{dz}\[{\mathrm{\α}}_s{\mathrm{\ρ}}_{s}v(h_s+\frac{1}{2}{v}^2)\]+\frac{d}{dz}\[(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_{l}v(h_l+\frac{1}{2}{v}^2)\]\\ ={\stackrel{\·}{q}}_{wm}-\[{\mathrm{\α}}_s{\mathrm{\ρ}}_s+\left(1-{\mathrm{\α}}_s\right){\mathrm{\ρ}}_l\]vg\end{array}---\left(4\right)$

其中，z为坐标轴，正向为工质的上升方向；ρ_s、ρ_l分别表示研究对象内饱和水蒸气的密度和饱和水的密度，密度的单位是kg/m³；α_s表示截面容积含汽率；v表示工质的流速，单位是m/s；τ_wm是壁面剪切应力，单位是Pa/m³；p_r是上升管中的压力，单位Pa；表示热流密度，单位是W/m³；g是重力加速度，单位是m/s²。h_l为饱和水的比焓。

饱和温度和工质压力之间存在着以下关系：

T＝T_b(1+a₁ln(P/P_b)+a₂ln²(P/P_b)+a₃ln³(P/P_b)+a₄ln⁴(P/P_b)) (5)

其中P是当前位置压力，T是当前位置的温度，T_b参考点的温度，P_b表示参考点的压力，a₁、a₂、a₃、a₄是常系数。

壁面的剪切应力可以通过下式求取：

${\mathrm{\τ}}_{wm}=\frac{2f}{d_t}{\mathrm{\ρ}}_m{v}^2---\left(6\right)$

f＝0.0014+0.125Re^-0.32 (7)

ρ_m＝α_sρ_s+(1-α_s)ρ_l (8)

式中Re表示雷诺数，f为摩擦系数，d_t为水冷壁上升管的管径，ρ_m为汽水混合物的密度。

水冷壁的热流密度可以通过下式求取：

${\stackrel{\·}{q}}_{wm}=\left(\frac{{\mathrm{\πd}}_t}{A_t}\right)U_{wm}(T_w-T)---\left(9\right)$

式中，A_t表示截面积，U_wm表示传热系数，T_w表示管壁的温度。

对(2)～(4)式进行展开，并将(5)～(9)式代入，可得到如下方程组：

$\left[\begin{array}{ccccc}{A}_{11}& A_{12}& ...& ...& A_{15}\\ A_{21}& A_{22}& ...& ...& A_{25}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ A_{51}& A_{52}& ...& ...& A_{55}\end{array}\right]\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=\left[\begin{array}{c}{B}_1\\ B_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ B_5\end{array}\right]---\left(10\right)$

其中，A₁₁＝(ρ_s-ρ_l)v，A₁₂＝α_sv，A₁₃＝α_sρ_s+(1-α_s)ρ_l，A₁₄＝0，A₁₅＝0；

$A_{21}=\{{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v-\{{\mathrm{\ρ}}_l\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\]\}v,$

$A_{22}={\mathrm{\α}}_{s}v\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{1}{2}{v}^2\],$

$A_{23}={\mathrm{\α}}_s{\mathrm{\ρ}}_s\[c_p^l(T-T_0)+{\mathrm{\λ}}_V+\frac{3}{2}{v}^2\]+(1-{\mathrm{\α}}_s){\mathrm{\ρ}}_l\[c_p^l(T-T_0)+\frac{3}{2}{v}^2\],$

A₂₅＝0；A₃₁＝(ρ_s-ρ_l)v²，A₃₂＝α_sv²，A₃₃＝2[α_sρ_s+(1-α_s)ρ_l]v，A₃₄＝0，A₃₅＝1；A₄₁＝0，A₄₂＝0，A₄₃＝0，A₄₄＝1，A₅₁＝0，A₅₂＝1，A₅₃＝0，y₁＝α_s，y₂＝ρ_s，y₃＝v，y₄＝T，y₅＝P；B₁＝0，

$B_2=\left(\frac{4}{d_t}\right)U_{wm}(T_w-T)-({\mathrm{\α}}_s{\mathrm{\ρ}}_s+(1-{\mathrm{\α}}_s\left){\mathrm{\ρ}}_l\right)vg,$

B₄＝0，B₅＝0。式中，λ_V表示汽化潜热，表示饱和水的定压比热，表示饱和蒸汽的定压比热，T₀表示指定状态下工质的温度。最终可以算得：

$\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=A^{-1}B---\left(11\right)$

2.2)通过对汽包-下降管-水冷壁中的所有工质进行质量守恒分析和能量守恒分析，可得：

$k_1\frac{dp}{dt}=Q-q_f(h_l-h_f)-q_s(h_s-h_l)---\left(12\right)$

$e_1\frac{{\mathrm{dV}}_{lt}}{dt}+e_2\frac{dp}{dt}=q_f-q_s---\left(13\right)$

其中，e₁＝ρ_l-ρ_s；

式中，p表示汽包压力，单位Pa；Q表示整个蒸发区的单位时间内的吸热量，单位是W；V_lt表示汽包-下降管-水冷壁中所有饱和水的体积，单位是m³；V_st表示回路中所有水蒸气的体积，单位是m³；q_f、q_s分别表示给水流量和主汽流量和下降管中循环工质的流量，流量的单位是kg/s；h_f、h_l、h_s分别表示给水比焓、饱和水的比焓和饱和水蒸气的比焓，单位是kJ/kg；m_t表示水冷壁金属和汽包金属的总质量，单位是kg；C_p表示金属的比热容；t_s表示饱和温度。

2.3)对汽包中液面下的工质进行质能守恒分析，忽略液面上饱和蒸汽的凝结，可得：

$\frac{d\left({\mathrm{\ρ}}_s{V}_s^{+}\right)}{dt}={\mathrm{\α}}_r{q}_r-q_s^{+}-q_{cd}---\left(14\right)$

$\frac{{\mathrm{dM}}_{dt}^{+}}{dt}=q_f+q_r-q_s^{+}-q_{dc}---\left(15\right)$

其中，

$q_s^{+}=A_d\left(\frac{V_s^{+}}{V_l+V_s^{+}}\right)v_s^{+}{\mathrm{\ρ}}_s;v_s^{+}=1.41{\[\frac{\mathrm{\σ}g({\mathrm{\ρ}}_l-{\mathrm{\ρ}}_s)}{{\mathrm{\ρ}}_l^2}\]}^{1/4};M_{dt}^{+}=V_l{\mathrm{\ρ}}_l+V_s^{+}{\mathrm{\ρ}}_s.$

式中，α_r表示上升管出口的截面质量含汽率；表示汽包内液面下工质的质量，单位是kg；q_r、q_cd、q_dc分别表示上升管中的流量、饱和蒸汽逸出液面的流量、饱和蒸汽凝结的流量和下降管中的流量，单位是kg/s；t_s表示饱和蒸汽的温度，单位是℃；m_d表示汽包金属的质量，单位kg；C_p表示水冷壁和汽包中金属的比热容，单位是kJ/(kg·℃)；A_d是汽包液位的截面面积，单位是m²；是饱和蒸汽逸出的速度，单位是m/s；σ是液面的表面张力，单位是N/m；V_l分别表示液面下饱和水蒸气的体积和饱和水的体积，单位是m³。

联立(12)～(15)式并进行展开，可以得到如下方程组：

$\left[\begin{array}{ccccc}{M}_{11}& M_{12}& ...& ...& M_{14}\\ M_{21}& M_{22}& ...& ...& M_{24}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ M_{41}& M_{42}& ...& ...& M_{44}\end{array}\right]\left[\begin{array}{c}{\mathrm{do}}_1/dt\\ {\mathrm{do}}_2/dt\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dt\end{array}\right]=\left[\begin{array}{c}{N}_1\\ N_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ N_4\end{array}\right]$

其中，M₁₂＝0，M₁₃＝0，M₁₄＝0；M₂₂＝ρ_l，M₂₃＝ρ_s，M₂₄＝0；M₃₂＝0，M₃₃＝ρ_s，M₃₄＝0；M₄₂＝0，M₄₃＝0，M₄₄＝ρ_l-ρ_s；

N₁＝Q-q_f(h_l-h_f)-q_s(h_s-h_l)，N₄＝q_f-q_s；o₁＝p，o₂＝V_l，o₃＝V_sd，o₄＝V_lt。

最终可得：

$\left[\begin{array}{c}{\mathrm{do}}_1/dz\\ {\mathrm{do}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dz\end{array}\right]=M^{-1}N---\left(16\right)$

采用Runge-Kutta数值解法求解微分方程(11)和(16)。对仿真时间进行离散化处理，将[0，T₀]区间内的仿真时间，离散成n个时间间隔为ΔT的区间，离散后的区间表示为[0,ΔT]，[ΔT，2ΔT]，……，[(n-1)ΔT，nΔT]。ΔT为0.01s～0.1s。

2.4)汽包水位是由汽包中的饱和水和饱和水蒸气的体积共同决定的，忽略汽包体积和汽包水位之间的非线性，可得汽包水位的表达式：

$l_d=\frac{V_l+V_s^{+}}{A_d}-l_0---\left(17\right)$

式中，l_d表示机理模型计算得到的汽包水位，单位mm；l₀是汽包水位的基准水位线。

ANFIS建模的步骤如下：

2.a)初始化粒子群。以聚类半径r_α作为粒子，15个粒子作为一个种群，每个粒子随机赋予[0.20.9]区间内的随机值，其中第i个粒子的位置的向量标示为r_i，i＝1，2，…，15；

2.b)以r_i为聚类半径，进行减法聚类分析，将每个数据点作为可能的聚类中心，并根据各个数据点周围的数据点密度来计算该点作为聚类中心的可能性。

每个数据点X_i作为聚类中心的可能性P_i由式(18)来定义：

$P_i=\underset{j=1}{\overset{m}{\Σ}}\exp (-\frac{\left|\right|X_i-X_j|{|}^2}{{(r_i/2)}^2})---\left(18\right)$

式中m表示n维输入空间中全部的数据点数，X_i＝[X_i1,X_i2,...,X_in]、X_j＝[X_j1,X_j2,...,X_jn]是具体的数据点，r_i是一个正数，定义了该点的邻域半径，||·||符号表示欧式距离。被选为聚类中心的点具有最高的数据点密度，同时该该数据点周围的点被排除作为聚类中心的可能性。第一个聚类中心为X_C1，数据点密度为P_c1。选出第一个聚类中心后，继续采用类似的方法确定下一个聚类中心，但需消除已有聚类中心的影响，修改密度指标的山峰函数如下：

$P_i=P_i-P_{c1}\exp (-\frac{\left|\right|X_i-X_{c1}|{|}^2}{{(r_{\mathrm{\β}}/2)}^2})---\left(19\right)$

其中，r_β＝1.5r_i。循环重复上述过程直到所有剩余数据点作为聚类中心的可能性低于某一阈值δ，即P_ck/P_c1<δ。

2.c)ANFIS模型训练。根据减法聚类算法得到的聚类中心，按照ANFIS模型结构训练汽包水位补偿模型；对于模糊神经网络模型的所有参数，采用混合最小二乘法的梯度下降算法进行学习。

2.d)计算适应度值。利用训练得到的预测模型计算汽包水位补偿值l^*，并与实际补偿值l_e进行比较，并以误差平方和作为粒子的适应度值MSE，适应度计算公式如下：

$MSE=\underset{i=1}{\overset{m}{\&Sigma;}}{(l_e-l^{*})}^2---\left(28\right)$

2.e)更新极值。以适应度值为评价指标，比较当代粒子与上一代粒子之间的适应度值大小，如果当前粒子的适应度值优于上一代，则将当前粒子的位置设置为个体极值，否则个体极值保持不变。同时获取当代所有粒子适应度值最优的粒子，并与上一代最优粒子进行比较，如果当代最优粒子的适应度值优于上一代最优粒子的适应度值，则将当代粒子的最优适应度值设置为全局最优值，否则全局最优值保持不变。

2.f)更新粒子。根据最新的个体极值和全局极值，按照(29)式和(30)式更新粒子的速度v_id(t)和位置x_id(t)。

v_id(t+1)＝ωv_id(t)+c₁r₁(p_id-x_id(t))+c₂r₂(p_gd-x_id(t)) (29)

x_id(t+1)＝x_id(t)+v_id(t+1) (30)

t是粒子群优化算法的寻优代数。更进一步，为了改善基本粒子群算法容易陷入局部极值和收敛速度慢的缺陷，在PSO算法的基础上引进了动态加速常数c₁、c₂和惯性权重ω：

$\mathrm{\&omega;}={\mathrm{\&omega;}}_{max}-\frac{{\mathrm{\&omega;}}_{max}-{\mathrm{\&omega;}}_{min}}{T_{\max }}t---\left(31\right)$

$c_1=R_1+\frac{R_2\&times;t}{T_{max}}---\left(32\right)$

$c_2=R_3-\frac{R_4\&times;t}{T_{max}}---\left(33\right)$

其中，T_max为最大寻优代数，ω_max为最大惯性权重，ω_min为最小惯性权重，R₁、R₂、R₃、R₄为常数。

2.g)算法停止条件算法判定。判断是否达到最大迭代次数或者到达预测精度的要求，如果没有达到则返回步骤2.b)，利用更新的聚类半径继续搜索，否则退出搜索。

2.h)利用最终寻优得到的聚类半径，对样本进行聚类分析和ANFIS模型训练，得到达到训练终止条件的ANFIS模型，即汽包水位的补偿模型。

第二信号采集模块。用于从数据库中选择需要预测汽包水位的运行工况，或者实时地采集当前锅炉的运行工况。

预测模块。该模块用于对指定的样本进行汽包水位的预测，或者对当前锅炉运行工况下的汽包水位进行实时预测。

结果显示模块。显示汽包水位的预测结果，或者对汽包水位的预测结果进行统计分析。

2.一种循环流化床生活垃圾焚烧入炉汽包水位的预测方法，其特征在于，该方法包括以下步骤：

1)选择混合模型的输入变量并采集训练样本。混合模型包括机理模型和ANFIS补偿模型。输入变量包括锅炉的设计参数、工质的物性参数和运行参数。设计参数包括水冷壁的当量管径、水冷壁的有效吸热长度、水冷壁的容积、循环水的流量、汽包的容积、汽包的平均截面面积、下降管的容积、汽包-下降管-水冷壁三者总共的金属质量。工质的物性参数包括压力、饱和温度、比焓、比容、汽化潜热、密度。运行参数包括给水流量、主蒸汽流量、垃圾给料量、辅助燃煤量、一次风量、二次风量以及炉膛中所有的温度测点。采集CFB生活垃圾焚烧锅炉在正常运行时的运行工况状态参数和操作变量，并组成混合模型输入变量的训练样本矩阵X(m×n)，m表示样本个数，n表示变量的个数。

2)训练样本预处理。对训练样本X(m×n)进行粗大误差处理和随机误差处理，剔除训练样本中的野值，排除异常工况，将训练样本输入变量经归一化处理后映射到[0，1]区间内，得到归一化后的训练样本X^*(m×n)；所述的异常工况包括锅炉停炉运行工况、炉膛压火运行状况和给料机堵塞工况；

3)利用机理建模的方法构建CFB生活垃圾焚烧锅炉汽包水位的机理模型。

3.1)对水冷壁中的工质进行质量守恒、动量守恒和能量守恒分析，可得：

$\frac{d}{dz}\&lsqb;{\mathrm{\&alpha;}}_s{\mathrm{\&rho;}}_{s}v+(1-{\mathrm{\&alpha;}}_s){\mathrm{\&rho;}}_{l}v\&rsqb;=0---\left(2\right)$

$\begin{array}{c}\frac{d}{dz}\&lsqb;{\mathrm{\&alpha;}}_s{\mathrm{\&rho;}}_s{v}^2+(1-{\mathrm{\&alpha;}}_s){\mathrm{\&rho;}}_l{v}^2\&rsqb;+\frac{{\mathrm{dp}}_r}{dz}\\ =-\&lsqb;{\mathrm{\&alpha;}}_s{\mathrm{\&rho;}}_s+\left(1-{\mathrm{\&alpha;}}_s\right){\mathrm{\&rho;}}_l\&rsqb;g-{\mathrm{\&tau;}}_{wm}\end{array}---\left(3\right)$

$\begin{array}{c}\frac{d}{dz}\&lsqb;{\mathrm{\&alpha;}}_s{\mathrm{\&rho;}}_{s}v(h_s+\frac{1}{2}{v}^2)\&rsqb;+\frac{d}{dz}\&lsqb;(1-{\mathrm{\&alpha;}}_s){\mathrm{\&rho;}}_{l}v(h_l+\frac{1}{2}{v}^2)\&rsqb;\\ ={\stackrel{\&CenterDot;}{q}}_{wm}-\&lsqb;{\mathrm{\&alpha;}}_s{\mathrm{\&rho;}}_s+\left(1-{\mathrm{\&alpha;}}_s\right){\mathrm{\&rho;}}_l\&rsqb;vg\end{array}---\left(4\right)$

其中，z为坐标轴，正向为工质的上升方向；ρ_s、ρ_l分别表示研究对象内饱和水蒸气的密度和饱和水的密度，密度的单位是kg/m³；α_s表示截面容积含汽率；v表示工质的流速，单位是m/s；τ_wm是壁面剪切应力，单位是Pa/m³；p_r是上升管中的压力，单位Pa；表示热流密度，单位是W/m³；g是重力加速度，单位是m/s²。h_l为饱和水的比焓。

饱和温度和工质压力之间存在着以下关系：

T＝T_b(1+a₁ln(P/P_b)+a₂ln²(P/P_b)+a₃ln³(P/P_b)+a₄ln⁴(P/P_b)) (5)

其中P是当前位置压力，T是当前位置的温度，T_b参考点的温度，P_b表示参考点的压力，a₁、a₂、a₃、a₄是常系数。

壁面的剪切应力可以通过下式求取：

${\mathrm{\&tau;}}_{wm}=\frac{2f}{d_t}{\mathrm{\&rho;}}_m{v}^2---\left(6\right)$

f＝0.0014+0.125Re^-0.32 (7)

ρ_m＝α_sρ_s+(1-α_s)ρ_l (8)

式中Re表示雷诺数，f为摩擦系数，d_t为水冷壁上升管的管径，ρ_m为汽水混合物的密度。

水冷壁的热流密度可以通过下式求取：

${\stackrel{\&CenterDot;}{q}}_{wm}=\left(\frac{{\mathrm{\&pi;d}}_t}{A_t}\right)U_{wm}(T_w-T)---\left(9\right)$

式中，A_t表示截面积，U_wm表示传热系数，T_w表示管壁的温度。

对(2)～(4)式进行展开，并将(5)～(9)式代入，可得到如下方程组：

$\left[\begin{array}{ccccc}{A}_{11}& A_{12}& ...& ...& A_{15}\\ A_{21}& A_{22}& ...& ...& A_{25}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ A_{51}& A_{52}& ...& ...& A_{55}\end{array}\right]\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=\left[\begin{array}{c}{B}_1\\ B_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ B_5\end{array}\right]---\left(10\right)$

其中，A₁₁＝(ρ_s-ρ_l)v，A₁₂＝α_sv，A₁₃＝α_sρ_s+(1-α_s)ρ_l，A₁₄＝0，A₁₅＝0；

$A_{21}=\{{\mathrm{\&rho;}}_s\&lsqb;c_p^l(T-T_0)+{\mathrm{\&lambda;}}_V+\frac{1}{2}{v}^2\&rsqb;\}v-\{{\mathrm{\&rho;}}_l\&lsqb;c_p^l(T-T_0)+{\mathrm{\&lambda;}}_V+\frac{1}{2}{v}^2\&rsqb;\}v,$

$A_{22}={\mathrm{\&alpha;}}_{s}v\&lsqb;c_p^l(T-T_0)+{\mathrm{\&lambda;}}_V+\frac{1}{2}{v}^2\&rsqb;,$

$A_{23}={\mathrm{\&alpha;}}_s{\mathrm{\&rho;}}_s\&lsqb;c_p^l(T-T_0)+{\mathrm{\&lambda;}}_V+\frac{3}{2}{v}^2\&rsqb;+(1-{\mathrm{\&alpha;}}_s){\mathrm{\&rho;}}_l\&lsqb;c_p^l(T-T_0)+\frac{3}{2}{v}^2\&rsqb;,$

A₂₅＝0；A₃₁＝(ρ_s-ρ_l)v²，A₃₂＝α_sv²，A₃₃＝2[α_sρ_s+(1-α_s)ρ_l]v，A₃₄＝0，A₃₅＝1；A₄₁＝0，A₄₂＝0，A₄₃＝0，A₄₄＝1，A₅₁＝0，A₅₂＝1，A₅₃＝0，y₁＝α_s，y₂＝ρ_s，y₃＝v，y₄＝T，y₅＝P；B₁＝0，

$B_2=\left(\frac{4}{d_t}\right)U_{wm}(T_w-T)-({\mathrm{\&alpha;}}_s{\mathrm{\&rho;}}_s+(1-{\mathrm{\&alpha;}}_s\left){\mathrm{\&rho;}}_l\right)vg,$

B₄＝0，B₅＝0。式中，λ_V表示汽化潜热，表示饱和水的定压比热，表示饱和蒸汽的定压比热，T₀表示指定状态下工质的温度。最终可以算得：

$\left[\begin{array}{c}{\mathrm{dy}}_1/dz\\ {\mathrm{dy}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{dy}}_5/dz\end{array}\right]=A^{-1}B---\left(11\right)$

3.2)通过对汽包-下降管-水冷壁中的所有工质进行质量守恒分析和能量守恒分析，可得：

$k_1\frac{dp}{dt}=Q-q_f(h_l-h_f)-q_s(h_s-h_l)---\left(12\right)$

$e_1\frac{{\mathrm{dV}}_{lt}}{dt}+e_2\frac{dp}{dt}=q_f-q_s---\left(13\right)$

其中，e₁＝ρ_l-ρ_s；

式中，p表示汽包压力，单位Pa；Q表示整个蒸发区的单位时间内的吸热量，单位是W；V_lt表示汽包-下降管-水冷壁中所有饱和水的体积，单位是m³；V_st表示回路中所有水蒸气的体积，单位是m³；q_f、q_s分别表示给水流量和主汽流量和下降管中循环工质的流量，流量的单位是kg/s；h_f、h_l、h_s分别表示给水比焓、饱和水的比焓和饱和水蒸气的比焓，单位是kJ/kg；m_t表示水冷壁金属和汽包金属的总质量，单位是kg；C_p表示金属的比热容；t_s表示饱和温度。

3.3)对汽包中液面下的工质进行质能守恒分析，忽略液面上饱和蒸汽的凝结，可得：

$\frac{d\left({\mathrm{\&rho;}}_s{V}_s^{+}\right)}{dt}={\mathrm{\&alpha;}}_r{q}_r-q_s^{+}-q_{cd}---\left(14\right)$

$\frac{{\mathrm{dM}}_{dt}^{+}}{dt}=q_f+q_r-q_s^{+}-q_{dc}---\left(15\right)$

其中，

$q_s^{+}=A_d\left(\frac{V_s^{+}}{V_l+V_s^{+}}\right)v_s^{+}{\mathrm{\&rho;}}_s;v_s^{+}=1.41{\&lsqb;\frac{\mathrm{\&sigma;}g({\mathrm{\&rho;}}_l-{\mathrm{\&rho;}}_s)}{{\mathrm{\&rho;}}_l^2}\&rsqb;}^{1/4};M_{dt}^{+}=V_l{\mathrm{\&rho;}}_l+V_s^{+}{\mathrm{\&rho;}}_s.$

式中，α_r表示上升管出口的截面质量含汽率；表示汽包内液面下工质的质量，单位是kg；q_r、q_cd、q_dc分别表示上升管中的流量、饱和蒸汽逸出液面的流量、饱和蒸汽凝结的流量和下降管中的流量，单位是kg/s；t_s表示饱和蒸汽的温度，单位是℃；m_d表示汽包金属的质量，单位kg；C_p表示水冷壁和汽包中金属的比热容，单位是kJ/(kg·℃)；A_d是汽包液位的截面面积，单位是m²；是饱和蒸汽逸出的速度，单位是m/s；σ是液面的表面张力，单位是N/m；V_l分别表示液面下饱和水蒸气的体积和饱和水的体积，单位是m³。

联立(12)～(15)式并进行展开，可以得到如下方程组：

$\left[\begin{array}{ccccc}{M}_{11}& M_{12}& ...& ...& M_{14}\\ M_{21}& M_{22}& ...& ...& M_{24}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ M_{41}& M_{42}& ...& ...& M_{44}\end{array}\right]\left[\begin{array}{c}{\mathrm{do}}_1/dt\\ {\mathrm{do}}_2/dt\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dt\end{array}\right]=\left[\begin{array}{c}{N}_1\\ N_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ N_4\end{array}\right]$

其中，M₁₂＝0，M₁₃＝0，M₁₄＝0；M₂₂＝ρ_l，M₂₃＝ρ_s，M₂₄＝0；M₃₂＝0，M₃₃＝ρ_s，M₃₄＝0；M₄₂＝0，M₄₃＝0，M₄₄＝ρ_l-ρ_s；

N₁＝Q-q_f(h_l-h_f)-q_s(h_s-h_l)，N₄＝q_f-q_s；o₁＝p，o₂＝V_l，o₃＝V_sd，o₄＝V_lt。

最终可得：

$\left[\begin{array}{c}{\mathrm{do}}_1/dz\\ {\mathrm{do}}_2/dz\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{do}}_4/dz\end{array}\right]=M^{-1}N---\left(16\right)$

采用Runge-Kutta数值解法求解微分方程(11)和(16)。对仿真时间进行离散化处理，将[0，T₀]区间内的仿真时间，离散成n个时间间隔为ΔT的区间，离散后的区间表示为[0,ΔT]，[ΔT，2ΔT]，……，[(n-1)ΔT，nΔT]。ΔT为0.01s～0.1s。

3.4)汽包水位是由汽包中的饱和水和饱和水蒸气的体积共同决定的，忽略汽包体积和汽包水位之间的非线性，可得汽包水位的表达式：

$l_d=\frac{V_l+V_s^{+}}{A_d}-l_0---\left(17\right)$

式中，l_d表示机理模型计算得到的汽包水位，单位mm；l₀是汽包水位的基准水位线。

(4)建立ANFIS补偿模型，步骤如下：

4.1)初始化粒子群。以聚类半径r_α作为粒子，15个粒子作为一个种群，每个粒子随机赋予[0.20.9]区间内的随机值，其中第i个粒子的位置的向量标示为r_i，i＝1，2，…，15；

4.2)以r_i为聚类半径，进行减法聚类分析，将每个数据点作为可能的聚类中心，并根据各个数据点周围的数据点密度来计算该点作为聚类中心的可能性。

每个数据点X_i作为聚类中心的可能性P_i由式(18)来定义：

$P_i=\underset{j=1}{\overset{m}{\&Sigma;}}\exp (-\frac{\left|\right|X_i-X_j|{|}^2}{{(r_i/2)}^2})---\left(18\right)$

式中m表示n维输入空间中全部的数据点数，X_i＝[X_i1,X_i2,...,X_in]、X_j＝[X_j1,X_j2,...,X_jn]是具体的数据点，r_i是一个正数，定义了该点的邻域半径，||·||符号表示欧式距离。被选为聚类中心的点具有最高的数据点密度，同时该该数据点周围的点被排除作为聚类中心的可能性。第一个聚类中心为X_C1，数据点密度为P_c1。选出第一个聚类中心后，继续采用类似的方法确定下一个聚类中心，但需消除已有聚类中心的影响，修改密度指标的山峰函数如下：

$P_i=P_i-P_{c1}\exp (-\frac{\left|\right|X_i-X_{c1}|{|}^2}{{(r_{\mathrm{\&beta;}}/2)}^2})---\left(19\right)$

其中，r_β＝1.5r_i。循环重复上述过程直到所有剩余数据点作为聚类中心的可能性低于某一阈值δ，即P_ck/P_c1<δ。

4.3)ANFIS模型训练。根据减法聚类算法得到的聚类中心，按照ANFIS模型结构训练汽包水位补偿模型；对于模糊神经网络模型的所有参数，采用混合最小二乘法的梯度下降算法进行学习。

4.4)计算适应度值。利用训练得到的预测模型计算汽包水位补偿值l^*，并与实际补偿值l_e进行比较，并以误差平方和作为粒子的适应度值MSE，适应度计算公式如下：

$MSE=\underset{i=1}{\overset{m}{\&Sigma;}}{(l_e-l^{*})}^2---\left(28\right)$

4.5)更新极值。以适应度值为评价指标，比较当代粒子与上一代粒子之间的适应度值大小，如果当前粒子的适应度值优于上一代，则将当前粒子的位置设置为个体极值，否则个体极值保持不变。同时获取当代所有粒子适应度值最优的粒子，并与上一代最优粒子进行比较，如果当代最优粒子的适应度值优于上一代最优粒子的适应度值，则将当代粒子的最优适应度值设置为全局最优值，否则全局最优值保持不变。

4.6)更新粒子。根据最新的个体极值和全局极值，按照(29)式和(30)式更新粒子的速度v_id(t)和位置x_id(t)。

v_id(t+1)＝ωv_id(t)+c₁r₁(p_id-x_id(t))+c₂r₂(p_gd-x_id(t)) (29)

x_id(t+1)＝x_id(t)+v_id(t+1) (30)

t是粒子群优化算法的寻优代数。更进一步，为了改善基本粒子群算法容易陷入局部极值和收敛速度慢的缺陷，在PSO算法的基础上引进了动态加速常数c₁、c₂和惯性权重ω：

$\mathrm{\&omega;}={\mathrm{\&omega;}}_{max}-\frac{{\mathrm{\&omega;}}_{max}-{\mathrm{\&omega;}}_{min}}{T_{\max }}t---\left(31\right)$

$c_1=R_1+\frac{R_2\&times;t}{T_{max}}---\left(32\right)$

$c_2=R_3-\frac{R_4\&times;t}{T_{max}}---\left(33\right)$

其中，T_max为最大寻优代数，ω_max为最大惯性权重，ω_min为最小惯性权重，R₁、R₂、R₃、R₄为常数。

4.7)算法停止条件算法判定。判断是否达到最大迭代次数或者到达预测精度的要求，如果没有达到则返回步骤2.b)，利用更新的聚类半径继续搜索，否则退出搜索。

4.8)利用最终寻优得到的聚类半径，对样本进行聚类分析和ANFIS模型训练，得到达到训练终止条件的ANFIS模型，即汽包水位的补偿模型。

5)汽包水位的预测。选择需要预测的工况，通过机理模型得到汽包水位的预测值l_d，通过ANFIS模型得到汽包水位的补偿值l^*，则混合模型的输出值即最终预测的汽包水位为