CN103472727B - 群智加权的农药焚烧炉有害物排放达标控制系统及方法 - Google Patents

Abstract

本发明公开了一种群智加权的农药焚烧炉有害物排放达标控制系统及方法。它包括焚烧炉、智能仪表、DCS系统、数据接口以及上位机；DCS系统包括控制站和数据库；用于测量易测变量的智能仪表与DCS系统连接，DCS系统通过数据接口与上位机连接。上位机首先对训练样本进行预处理和模糊化后，得到新的输入矩阵，再采用加权最小二乘支持向量机对新的输入矩阵建立回归模型，获得预测输出，然后对预测输出进行反模糊化，最后引入粒子群算法进行优化，获得优化的系统输出。本发明具有在线测量COD、有效监测COD是否超标、控制COD排放达标、抗噪声和泛化能力强、所需样本数少、计算速度快、在线自动优化、鲁棒性强等优点。

Description

群智加权的农药焚烧炉有害物排放达标控制系统及方法

技术领域

本发明涉及农药生产领域，尤其涉及一种群智加权的农药焚烧炉有害物排放达标控制系统及方法。

背景技术

我国是农药生产和使用大国，农药生产企业已达4100家左右，其中原药生产企业为500多家，国家农业部统计数据显示2008年1～11月农药总产量达171.1万吨。我国农药品种结构的不合理性更加大了环境治理的难度。据不完全统计，全国农药工业每年排放的废水约为15亿吨。其中，处理达标的仅占已处理的1%。焚烧法是目前处理农药残液和废渣最有效、彻底、应用最普遍的方法。焚烧后废水的化学耗氧量(COD)是农药废液焚烧有害物排放的最重要指标，但是其无法在线测量，离线测量一次需要四五个小时，无法及时反映工况变化和指导实际生产。因此，在实际焚烧过程中，COD严重超标。

发明内容

为了克服已有的焚烧炉过程COD无法在线测量、COD严重超标的不足，本发明提供一种群智加权的农药焚烧炉有害物排放达标控制系统及方法，其具有在线测量COD、有效监测COD是否超标、控制COD排放达标、抗噪声和泛化能力强、所需样本数少、计算速度快、在线自动优化、鲁棒性强等优点。

本发明解决其技术问题所采用的技术方案是：

群智加权的农药焚烧炉有害物排放达标控制系统，包括焚烧炉、智能仪表、DCS系统、数据接口以及上位机，所述的DCS系统包括控制站和数据库；所述现场智能仪表与DCS系统连接，所述DCS系统与上位机连接，所述的上位机包括：

数据预处理模块，用于将从DCS数据库输入的模型训练样本进行预处理，对训练样本中心化，即减去样本的平均值，然后对其进行标准化：

计算均值： $\stackrel{\‾}{\mathrm{TX}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{TX}}_i---\left(1\right)$

计算方差： ${\mathrm{\σ}}_x^2=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{TX}}_i-\stackrel{\‾}{\mathrm{TX}})---\left(2\right)$

标准化： $X=\frac{\mathrm{TX}-\stackrel{\‾}{\mathrm{TX}}}{{\mathrm{\σ}}_x}---\left(3\right)$

其中，TX_i为第i个训练样本，是从DCS数据库中采集的生产正常时的关键变量、化学耗氧量（COD）和相应的使COD排放达标时的操作变量的数据，N为训练样本数，为训练样本的均值，X为标准化后的训练样本。σ_x表示训练样本的标准差，σ² _x表示训练样本的方差。

模糊方程模块，对从数据预处理模块传过来的标准化后的训练样本X，进行模糊化。设模糊方程系统中有c^*个模糊群，模糊群k、j的中心分别为v_k、v_j，则第i个标准化后的训练样本X_i对于模糊群k的隶属度μ_ik为：

${\mathrm{\μ}}_{\mathrm{ik}}={\left(\underset{j=1}{\overset{c^{*}}{\mathrm{\Σ}}}{\left(\frac{\left|\right|X_i-v_k\left|\right|}{\left|\right|X_i-v_j\left|\right|}\right)}^{\frac{2}{n-1}}\right)}^{-1}---\left(4\right)$

式中，n为模糊分类过程中需要的分块矩阵指数，通常取作2,||·||为范数表达式。

使用以上隶属度值或者它的变形以获得新的输入矩阵，对于模糊群k，其输入矩阵变形为：

Φ_ik(X_i,μ_ik)=[1 func(μ_ik) X_i] （5）

其中func(μ_ik)为隶属度值μ_ik的变形函数，一般取exp(μ_ik)等，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵。加权最小二乘支持向量机作为模糊方程系统的局部方程，对每个模糊群进行优化拟合。设模型训练样本的第i个目标输出为O_i，加权最小二乘支持向量机通过变换把拟合问题等价于如下二次规划问题：

$\min R(w,\mathrm{\ξ})=\frac{1}{2}{w}^{T}w+\frac{1}{2}\mathrm{\γ}{\mathrm{\Σ}}_{i=1}^N{\mathrm{\ω}}_i{\mathrm{\ξ}}_i^2---\left(6\right)$

其中，R(w,ξ)是优化问题的目标函数，minR(w,ξ)是优化问题的目标函数的最小值,N是训练样本数，ξ={ξ₁,…,ξ_N}是松弛变量，ξ_k是松弛变量的第k个分量，w是加权最小二乘支持向量机超平面的法向量，b是相应的偏移量，而ω_i,i=1,…,N和γ分别是加权最小二乘支持 向量机的权重和惩罚因子，是加权最小二乘支持向量机松弛变量的第k个分量ξ_k标准差的估计，c₁为常数，这里取2.5，c₂为常数，这里取3，上标T表述矩阵的转置。μ_ik表示第i个标准化后的训练样本X_i对于模糊群k的隶属度，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵。由（6）（7）（8）式可推导出模糊群k在训练样本i的输出为：

${\hat{y}}_{\mathrm{ik}}=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{k}K<{\mathrm{\&Phi;}}_k(X,\mathrm{\&mu;}),{\mathrm{\&Phi;}}_{\mathrm{ik}}(X_i,{\mathrm{\&mu;}}_{\mathrm{ik}})>+b_i---\left(9\right)$

其中，为第k个群智加权模糊系统输出层的预测输出，μ_mk表示第m个训练样本X_m对于模糊群k的隶属度，Φ_mk(X_m,μ_mk)表示第m个输入变量X_m及其模糊群k的隶属度μ_mk所对应的新的输入矩阵。α_k，(k=1,…,N)是对应的拉格朗日乘子的第k个分量，K<·>是加权最小二乘支持向量机的核函数，这里K<·>取线性核函数。

由反模糊方法中的重心法得到最后的模糊方程系统的输出：

${\hat{y}}_i=\frac{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{\mathrm{ik}}{\hat{y}}_{\mathrm{ik}}}{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{\mathrm{ik}}}---\left(10\right)$

粒子群算法优化模块，用于采用粒子群算法对模糊方程中加权最小二乘支持向量机局部方程的惩罚因子和误差容限值进行优化，具体步骤如下：

①确定粒子群的优化参数为加权最小二乘支持向量机局部方程的惩罚因子和误差容限值、粒子群个体数目popsize、最大循环寻优次数iter_max、第p个粒子的初始位置r_p、初始速度v_p、局部最优值Lbest_p以及整个粒子群的全局最优值Gbest。

②设定优化目标函数，将其转换为适应度，对每个局部模糊方程进行评价；通过相应的误差函数计算适应度函数，并认为误差大的粒子适应度小，粒子p的适应度函数表示为：

f_p=1/(E_p+1) （11）

式中，E_p是模糊方程系统的误差函数，表示为：

$E_p=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-O_i)}^2}---\left(12\right)$

式中，是模糊方程系统的预测输出，O_i为模糊方程系统的目标输出；

③按照如下公式，循环更新每个粒子的速度和位置，

v_p(iter+1)=ω×v_p(iter)+ m₁a₁(Lbest_p-r_p(iter))+m₂a₂(Gbest-r_p(iter))

(13)

r_p(iter+1)=r_p(iter)+v_p(iter+1) （14）

式中，v_p表示更新粒子p的速度，r_p表示更新粒子p的位置，Lbest_p表示更新粒子p的个体最优值，Gbest即为整个粒子群的全局最优值，iter表示循环次数，ω是粒子群算法中的惯性权重，m₁、m₂是对应的加速系数，a₁、a₂是[0,1]之间的随机数；

④对于粒子p，如果新的适应度大于原来的个体最优值，更新粒子的个体最优值：

Lbest_p=f_p （15）

⑤如果粒子p的个体最优值Lbest_p大于原来的粒子群全局最优值Gbest，更新原来的粒子群全局最优值Gbest:

Gbest=Lbest_p （16）

⑥判断是否满足性能要求，若是，结束寻优，得到一组优化的模糊方程的局部方程参数；否则返回步骤③，继续迭代寻优，直至达到最大迭代次数iter_max。

迭代终止时的Gbest即为表示第i个标准化后的训练样本X_i的COD预报值和使COD排放达标的操作变量值。

作为优选的一种方案：所述的上位机还包括：判别模型更新模块，用于按设定的采样时间间隔，采集现场智能仪表信号，将得到的实测化学耗氧量与函数预报值比较，如果相对误差大于10%或实测COD数据不达标，则将DCS数据库中生产正常时的达标的新数据加入训练样本数据，更新模型。

进一步，所述的上位机还包括：结果显示模块，用于将COD预报值和使COD排放达标的操作变量值传给DCS系统，并在DCS的控制站显示过程状态，同时通过DCS系统和现场总线将过程状态信息传递到现场操作站进行显示；同时，DCS系统将所得到的使COD排放达标的操作变量值作为新的操作变量设定值，自动执行COD排放达标操作。再进一步，所述的关键变量包括进入焚烧炉的废液流量、进入焚烧炉的空气流量和进入焚烧炉的燃料流量；所述的操作变量包括进入焚烧炉的空气流量和进入焚烧炉的燃料流量。信号采集模块，用于依照设定的每次采样的时间间隔，从数据库中采集数据。

所述的农药焚烧炉有害物排放达标控制系统实现的有害物排放达标控制方法，所述的控制方法具体实现步骤如下：

1）、对农药焚烧炉化学耗氧量排放过程对象，根据工艺分析和操作分析，确定所用的关键变量，从DCS数据库中采集生产正常时所述变量的数据作为训练样本TX的输入矩阵，采集对应的COD和使COD排放达标的操作变量数据作为输出矩阵Y；

2）、将从DCS数据库输入的模型训练样本进行预处理，对训练样本中心化，即减去样本 的平均值，然后对其进行标准化，使得其均值为0，方差为1。该处理采用以下算式过程来完成：

2.1）计算均值： $\stackrel{\&OverBar;}{\mathrm{TX}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{TX}}_i---\left(1\right)$

2.2）计算方差： ${\mathrm{\&sigma;}}_x^2=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}({\mathrm{TX}}_i-\stackrel{\&OverBar;}{\mathrm{TX}})---\left(2\right)$

2.3）标准化： $X=\frac{\mathrm{TX}-\stackrel{\&OverBar;}{\mathrm{TX}}}{{\mathrm{\&sigma;}}_x}---\left(3\right)$

其中，TX_i为第i个训练样本，是从DCS数据库中采集的生产正常时的关键变量、化学耗氧量（COD）和相应的使COD排放达标时的操作变量的数据，N为训练样本数，为训练样本的均值，X为标准化后的训练样本。σ_x表示训练样本的标准差，σ² _x表示训练样本的方差。

3）对从数据预处理模块传过来的训练样本，进行模糊化。设模糊方程系统中有c^*个模糊群，模糊群k、j的中心分别为v_k、v_j，则第i个标准化后的训练样本X_i对于模糊群k的隶属度μ_ik为：

${\mathrm{\&mu;}}_{\mathrm{ik}}={\left(\underset{j=1}{\overset{c^{*}}{\mathrm{\&Sigma;}}}{\left(\frac{\left|\right|X_i-v_k\left|\right|}{\left|\right|X_i-v_j\left|\right|}\right)}^{\frac{2}{n-1}}\right)}^{-1}---\left(4\right)$

式中，n为模糊分类过程中需要的分块矩阵指数，通常取作2,||·||为范数表达式。

使用以上隶属度值或者它的变形以获得新的输入矩阵，对于模糊群k，其输入矩阵变形为：

Φ_ik(X_i,μ_ik)=[1 func(μ_ik) X_i] （5）

其中func(μ_ik)为隶属度值μ_ik的变形函数，一般取exp(μ_ik)等，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵。

加权最小二乘支持向量机作为模糊方程系统的局部方程，对每个模糊群进行优化拟合。设模型训练样本的第i个目标输出为O_i，加权最小二乘支持向量机通过变换把拟合问题等价于如下二次规划问题：

$\min R(w,\mathrm{\&xi;})=\frac{1}{2}{w}^{T}w+\frac{1}{2}\mathrm{\&gamma;}{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&omega;}}_i{\mathrm{\&xi;}}_i^2---\left(6\right)$

其中，R(w,ξ)是优化问题的目标函数，minR(w,ξ)是优化问题的目标函数的最小值,N是训练样本数，ξ={ξ₁,…,ξ_N}是松弛变量，ξ_k是松弛变量的第k个分量，w是加权最小二乘支持向量机超平面的法向量，b是相应的偏移量，而ω_i,i=1,…,N和γ分别是加权最小二乘支持向量机的权重和惩罚因子，是加权最小二乘支持向量机松弛变量的第k个分量ξ_k标准差的估计，c₁为常数，这里取2.5，c₂为常数，这里取3，上标T表述矩阵的转置。μ_ik表示第i个标准化后的训练样本X_i对于模糊群k的隶属度，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵。由（6）（7）（8）式可推导出模糊群k在训练样本i的输出为：

${\hat{y}}_{\mathrm{ik}}=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{k}K<{\mathrm{\&Phi;}}_k(X,\mathrm{\&mu;}),{\mathrm{\&Phi;}}_{\mathrm{jk}}(X_i,{\mathrm{\&mu;}}_{\mathrm{ik}})>+b_i---\left(9\right)$

其中，为第k个群智加权模糊系统输出层的预测输出，μ_mk表示第m个训练样本X_m对于模糊群k的隶属度，Φ_mk(X_m,μ_mk) 表示第m个输入变量X_m及其模糊群k的隶属度μ_mk所对应的新的输入矩阵。α_k，(k=1,…,N)是对应的拉格朗日乘子的第k个分量，K<·>是加权最小二乘支持向量机的核函数，这里K<·>取线性核函数。

由反模糊方法中的重心法得到最后的模糊方程系统的输出：

${\hat{y}}_i=\frac{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{\mathrm{ik}}{\hat{y}}_{\mathrm{ik}}}{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{\mathrm{ik}}}---\left(10\right)$

4)采用粒子群算法对模糊方程中加权最小二乘支持向量机局部方程的惩罚因子和误差容限值进行优化，具体步骤如下：

①确定粒子群的优化参数为加权最小二乘支持向量机局部方程的惩罚因子和误差容限值、粒子群个体数目popsize、最大循环寻优次数iter_max、第p个粒子的初始位置r_p、初始速度v_p、局部最优值Lbest_p以及整个粒子群的全局最优值Gbest。

②设定优化目标函数，将其转换为适应度，对每个局部模糊方程进行评价；通过相应的误差函数计算适应度函数，并认为误差大的粒子适应度小，粒子p的适应度函数表示为：

f_p=1/(E_p+1) （11）

式中，E_p是模糊方程系统的误差函数，表示为：

$E_p=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-O_i)}^2}---\left(12\right)$

式中，是模糊方程系统的预测输出，O_i为模糊方程系统的目标输出；

③按照如下公式，循环更新每个粒子的速度和位置，

v_p(iter+1)=ω×v_p(iter)+m₁a₁(Lbest_p-r_p(iter))+m₂a₂(Gbest-r_p(iter))

（13）

r_p(iter+1)=r_p(iter)+v_p(iter+1) （14）

式中，v_p表示更新粒子p的速度，r_p表示更新粒子p的位置，Lbest_p表示更新粒子p的个体最优值，Gbest即为整个粒子群的全局最优值，iter表示循环次数，ω是粒子群算法中的惯性权重，m₁、m₂是对应的加速系数，a₁、a₂是[0,1]之间的随机数；

④对于粒子p，如果新的适应度大于原来的个体最优值，更新粒子的个体最优值：

Lbest_p=f_p （15）

⑤如果粒子p的个体最优值Lbest_p大于原来的粒子群全局最优值Gbest，更新原来的粒子群全局最优值Gbest:

Gbest=Lbest_p （16）

⑥判断是否满足性能要求，若是，结束寻优，得到一组优化的模糊方程的局部方程参数；否则返回步骤③，继续迭代寻优，直至达到最大迭代次数iter_max。

迭代终止时得到Gbest即为表示第i个标准化后的训练样本X_i的COD预报值和使COD排放达标的操作变量值。

作为优选的一种方案：所述的方法还包括：5）按设定的采样时间间隔，采集现场智能仪表信号，将得到的实测化学耗氧量与函数预报值比较，如果相对误差大于10%或实测COD数据不达标，则将DCS数据库中生产正常时的达标的新数据加入训练样本数据。

进一步，在所述的步骤4）中，将COD预报值和使COD排放达标的操作变量值传给DCS系统，并在DCS的控制站显示过程状态，同时通过DCS系统和现场总线将过程状态信息传递到现场操作站进行显示；同时，DCS系统将所得到的使COD排放达标的操作变量值作为新的操作变量设定值，自动执行COD排放达标操作。再进一步，所述的关键变量包括进入焚烧炉的废液流量，进入焚烧炉的空气流量，进入焚烧炉的燃料流量；所述的操作变量包括进入焚烧炉的空气流量和进入焚烧炉的燃料流量。

本发明的技术构思为：本发明提供一种群智加权的农药焚烧炉有害物排放达标控制系统及方法，寻找到使得COD排放达标的操作变量值。

本发明的有益效果主要表现在：1、建立了系统关键变量和COD排放之间定量关系的在线软测量模型；2、迅速找到使得COD排放达标的操作条件。

附图说明

图1是本发明所提出的系统的硬件结构图；

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

具体实施方式

下面结合附图对本发明作进一步描述。本发明实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。

实施例1

参照图1、图2，群智加权的农药焚烧炉有害物排放达标控制系统，包括与焚烧炉1连接的现场智能仪表2、DCS系统以及上位机6，所述DCS系统包括数据接口3、控制站5和数据库4，所述现场智能仪表2与数据接口3连接，所述数据接口与控制站5、数据库4和上位机6连接，所述的上位机6包括：

数据预处理模块7，用于将从DCS数据库输入的模型训练样本进行预处理，对训练样本中心化，即减去样本的平均值，然后对其进行标准化：

计算均值： $\stackrel{\&OverBar;}{\mathrm{TX}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{TX}}_i---\left(1\right)$

计算方差： ${\mathrm{\&sigma;}}_x^2=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}({\mathrm{TX}}_i-\stackrel{\&OverBar;}{\mathrm{TX}})---\left(2\right)$

标准化： $X=\frac{\mathrm{TX}-\stackrel{\&OverBar;}{\mathrm{TX}}}{{\mathrm{\&sigma;}}_x}---\left(3\right)$

其中，TX_i为第i个训练样本，是从DCS数据库中采集的生产正常时的关键变量、化学耗氧量（COD）和相应的使COD排放达标时的操作变量的数据，N为训练样本数，为训练样本的均值，X为标准化后的训练样本。σ_x表示训练样本的标准差，σ² _x表示训练样本的方差。模糊方程模块，对从数据预处理模块传过来的标准化后的训练样本X，进行模糊化。设模糊方程系统中有c^*个模糊群，模糊群k、j的中心分别为v_k、v_j，则第i个标准化后的训练样本 X_i对于模糊群k的隶属度μ_ik为：

${\mathrm{\&mu;}}_{\mathrm{ik}}={\left(\underset{j=1}{\overset{c^{*}}{\mathrm{\&Sigma;}}}{\left(\frac{\left|\right|X_i-v_k\left|\right|}{\left|\right|X_i-v_j\left|\right|}\right)}^{\frac{2}{n-1}}\right)}^{-1}---\left(4\right)$

式中，n为模糊分类过程中需要的分块矩阵指数，通常取作2,||·||为范数表达式。

使用以上隶属度值或者它的变形以获得新的输入矩阵，对于模糊群k，其输入矩阵变形为：

Φ_ik(X_i,μ_ik)=[1 func(μ_ik) X_i] （5）

其中func(μ_ik)为隶属度值μ_ik的变形函数，一般取exp(μ_ik)等，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵。加权最小二乘支持向量机作为模糊方程系统的局部方程，对每个模糊群进行优化拟合。设模型训练样本的第i个目标输出为O_i，加权最小二乘支持向量机通过变换把拟合问题等价于如下二次规划问题：

$\min R(w,\mathrm{\&xi;})=\frac{1}{2}{w}^{T}w+\frac{1}{2}\mathrm{\&gamma;}{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&omega;}}_i{\mathrm{\&xi;}}_i^2---\left(6\right)$

其中，R(w,ξ)是优化问题的目标函数，minR(w,ξ)是优化问题的目标函数的最小值,N是训练样本数，ξ={ξ₁,…,ξ_N}是松弛变量，ξ_k是松弛变量的第k个分量，w是加权最小二乘支持向量机超平面的法向量，b是相应的偏移量，而ω_i,i=1,…,N和γ分别是加权最小二乘支持向量机的权重和惩罚因子，是加权最小二乘支持向量机松弛变量的第k个分量ξ_k标准差的估计，c₁为常数，这里取2.5，c₂为常数，这里取3，上标T表述矩阵的转置。μ_ik表示第i个标准化后的训练样本X_i对于模糊群k的隶属度，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵。由（6）（7）（8）式可推导出模糊群k在训练样本i的输出为：

${\hat{y}}_{\mathrm{ik}}=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{k}K<{\mathrm{\&Phi;}}_k(X,\mathrm{\&mu;}),{\mathrm{\&Phi;}}_{\mathrm{ik}}(X_i,{\mathrm{\&mu;}}_{\mathrm{ik}})>+b_i---\left(9\right)$

其中，为第k个群智加权模糊系统输出层的预测输出，μ_mk表示第m个训练样本X_m对于模糊群k的隶属度，表示第m个输入变量X_m及其模糊群k的隶属度μ_mk所对应的新的输入矩阵。α_k，(k=1,…,N)是对应的拉格朗日乘子的第k个分量，K<·>是加权最小二乘支持向量机的核函数，这里K<·>取线性核函数。

由反模糊方法中的重心法得到最后的模糊方程系统的输出：

${\hat{y}}_i=\frac{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{\mathrm{ik}}{\hat{y}}_{\mathrm{ik}}}{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{\mathrm{ik}}}---\left(10\right)$

粒子群算法（PSO）模块，用于采用粒子群算法对模糊方程中加权最小二乘支持向量机局部方程的惩罚因子和误差容限值进行优化，具体步骤如下：

①确定粒子群的优化参数为加权最小二乘支持向量机局部方程的惩罚因子和误差容限值、粒子群个体数目popsize、最大循环寻优次数iter_max、第p个粒子的初始位置r_p、初始速度v_p、局部最优值Lbest_p以及整个粒子群的全局最优值Gbest。

②设定优化目标函数，将其转换为适应度，对每个局部模糊方程进行评价；通过相应的误差函数计算适应度函数，并认为误差大的粒子适应度小，粒子p的适应度函数表示为：

f_p=1/(E_p+1) （11）

式中，E_p是模糊方程系统的误差函数，表示为：

$E_p=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-O_i)}^2}---\left(12\right)$

式中，是模糊方程系统的预测输出，O_i为模糊方程系统的目标输出；

③按照如下公式，循环更新每个粒子的速度和位置，

v_p(iter+1)=ω×v_p(iter)+m₁a₁(Lbest_p-r_p(iter))+m₂a₂(Gbest-r_p(iter))

（13）

r_p(iter+1)=r_p(iter)+v_p(iter+1) （14）

式中，v_p表示更新粒子p的速度，r_p表示更新粒子p的位置，Lbest_p表示更新粒子p的个体最优值，Gbest即为整个粒子群的全局最优值，iter表示循环次数，ω是粒子群算法中的惯性权重，m₁、m₂是对应的加速系数，a₁、a₂是[0,1]之间的随机数；

④对于粒子p，如果新的适应度大于原来的个体最优值，更新粒子的个体最优值：

Lbest_p=f_p （15）

⑤如果粒子p的个体最优值Lbest_p大于原来的粒子群全局最优值Gbest，更新原来的粒子群全局最优值Gbest:

Gbest=Lbest_p （16）

⑥判断是否满足性能要求，若是，结束寻优，得到一组优化的模糊方程的局部方程参数；否则返回步骤③，继续迭代寻优，直至达到最大迭代次数iter_max。

迭代终止时的Gbest即为使COD即为表示第i个标准化后的训练样本X_i的COD预报值和使COD排放达标的操作变量值。

所述上位机6还包括：信号采集模块11，用于依照设定的每次采样的时间间隔，从数据库中采集数据；

所述的上位机6还包括：判别模型更新模块12，按设定的采样时间间隔，采集现场智能仪表信号，将得到的实测COD与函数预报值比较，如果相对误差大于10%或实测COD数据不达标，则将DCS数据库中生产正常时的达标的新数据加入训练样本数据，更新模型。所述的关键变量包括进入焚烧炉的废液流量、进入焚烧炉的空气流量和进入焚烧炉的燃料流量；所述的操作变量包括进入焚烧炉的空气流量和进入焚烧炉的燃料流量。

所述系统还包括DCS系统，所述的DCS系统由数据接口3、控制站5、数据库4构成；智能仪表2、DCS系统、上位机6通过现场总线依次相连；上位机6还包括结果显示模块10，用于将COD预报值和使COD排放达标的操作变量值传给DCS系统，并在DCS的控制站显示过程状态，同时通过DCS系统和现场总线将过程状态信息传递到现场操作站进行显示；同时，DCS系统将所得到的使COD排放达标的操作变量值作为新的操作变量设定值，自动执行COD排放达标操作。

当废液焚烧过程已配有DCS系统时，样本实时动态数据的检测、存储利用DCS系统的实时和历史数据库，得到COD预报值和使COD排放达标的操作变量值的功能主要在上位机上完成。

当废液焚烧过程没有配备DCS系统时，采用数据存储器来替代DCS系统的实时和历史数据库的数据存储功能，并将得到COD预报值和使COD排放达标的操作变量值的功能系统制造成包括I/O元件、数据存储器、程序存储器、运算器、显示模块几大构件的不依赖于DCS系统的一个独立的完整的片上系统，在不管焚烧过程是否配备DCS的情况下，都能够独立使用，更有益于推广使用。

实施例2

参照图1、图2，群智加权的农药焚烧炉有害物排放达标控制方法，所述控制方法具体实 现步骤如下：

1）、对农药焚烧炉化学耗氧量排放过程对象，根据工艺分析和操作分析，确定所用的关键变量，从DCS数据库中采集生产正常时所述变量的数据作为训练样本TX的输入矩阵，采集对应的COD和使COD排放达标的操作变量数据作为输出矩阵Y；

2）、将从DCS数据库输入的模型训练样本进行预处理，对训练样本中心化，即减去样本的平均值，然后对其进行标准化，使得其均值为0，方差为1。该处理采用以下算式过程来完成：

2.1）计算均值： $\stackrel{\&OverBar;}{\mathrm{TX}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{TX}}_i---\left(1\right)$

2.2）计算方差： ${\mathrm{\&sigma;}}_x^2=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}({\mathrm{TX}}_i-\stackrel{\&OverBar;}{\mathrm{TX}})---\left(2\right)$

2.3）标准化： $X=\frac{\mathrm{TX}-\stackrel{\&OverBar;}{\mathrm{TX}}}{{\mathrm{\&sigma;}}_x}---\left(3\right)$

其中，TX_i为第i个训练样本，是从DCS数据库中采集的生产正常时的关键变量、化学耗氧量（COD）和相应的使COD排放达标时的操作变量的数据，N为训练样本数，为训练样本的均值，X为标准化后的训练样本。σ_x表示训练样本的标准差，σ² _x表示训练样本的方差。

3）对从数据预处理模块传过来的经过标准化后的训练样本，进行模糊化。设模糊方程系统中有c^*个模糊群，模糊群k、j的中心分别为v_k、v_j，则第i个标准化后的训练样本X_i对于模糊群k的隶属度μ_ik为：

${\mathrm{\&mu;}}_{\mathrm{ik}}={\left(\underset{j=1}{\overset{c^{*}}{\mathrm{\&Sigma;}}}{\left(\frac{\left|\right|X_i-v_k\left|\right|}{\left|\right|X_i-v_j\left|\right|}\right)}^{\frac{2}{n-1}}\right)}^{-1}---\left(4\right)$

式中，n为模糊分类过程中需要的分块矩阵指数，通常取作2,||·||为范数表达式。

使用以上隶属度值或者它的变形以获得新的输入矩阵，对于模糊群k，其输入矩阵变形为：

Φ_ik(X_i,μ_ik)=[1 func(μ_ik) X_i] （5）

其中func(μ_ik)为隶属度值μ_ik的变形函数，一般取exp(μ_ik)等，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵。

加权最小二乘支持向量机作为模糊方程系统的局部方程，对每个模糊群进行优化拟合。设模型训练样本的第i个目标输出为O_i，加权最小二乘支持向量机通过变换把拟合问题等价 于如下二次规划问题：

$\min R(w,\mathrm{\&xi;})=\frac{1}{2}{w}^{T}w+\frac{1}{2}\mathrm{\&gamma;}{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&omega;}}_i{\mathrm{\&xi;}}_i^2---\left(6\right)$

其中，R(w,ξ)是优化问题的目标函数，minR(w,ξ)是优化问题的目标函数的最小值,N是训练样本数，ξ={ξ₁,…,ξ_N}是松弛变量，ξ_k是松弛变量的第k个分量，w是加权最小二乘支持向量机超平面的法向量，b是相应的偏移量，而ω_i,i=1,…,N和γ分别是加权最小二乘支持向量机的权重和惩罚因子，是加权最小二乘支持向量机松弛变量的第k个分量ξ_k标准差的估计，c₁为常数，这里取2.5，c₂为常数，这里取3，上标T表述矩阵的转置。μ_ik表示第i个标准化后的训练样本X_i对于模糊群k的隶属度，Φ_ik(X_i,μ_ik) 表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵。由（6）（7）（8）式可推导出模糊群k在训练样本i的输出为：

${\hat{y}}_{\mathrm{ik}}=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{k}K<{\mathrm{\&Phi;}}_k(X,\mathrm{\&mu;}),{\mathrm{\&Phi;}}_{\mathrm{ik}}(X_i,{\mathrm{\&mu;}}_{\mathrm{ik}})>+b_i---\left(9\right)$

其中，为第k个群智加权模糊系统输出层的预测输出，μ_mk表示第m个训练样本X_mX_m对于模糊群k的隶属度，Φ_mk(X_m,μ_mk)表示第m个输入变量X_m及其模糊群k的隶属度μ_mk所对应的新的输入矩阵。α_k，(k=1,…,N)是对应的拉格朗日乘子的第k个分量，K<·>是加权最小二乘支持向量机的核函数，这里K<·>取线性核函数。

由反模糊方法中的重心法得到最后的模糊方程系统的输出：

${\hat{y}}_i=\frac{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{\mathrm{ik}}{\hat{y}}_{\mathrm{ik}}}{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{\mathrm{ik}}}---\left(10\right)$

4)采用粒子群算法对模糊方程中加权最小二乘支持向量机局部方程的惩罚因子和误差容限值进行优化，具体步骤如下：

①确定粒子群的优化参数为加权支持向量机局部方程的惩罚因子和误差容限值、粒子群个体数目popsize、最大循环寻优次数iter_max、第p个粒子的初始位置r_p、初始速度v_p、局部最 优值Lbest_p以及整个粒子群的全局最优值Gbest。

②设定优化目标函数，将其转换为适应度，对每个局部模糊方程进行评价；通过相应的误差函数计算适应度函数，并认为误差大的粒子适应度小，粒子p的适应度函数表示为：

f_p=1/(E_p+1) （11）式中，E_p是模糊方程系统的误差函数，表示为：

$E_p=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-O_i)}^2}---\left(12\right)$ 式中，是模糊方程系统的预测输出，O_i为模糊方程系统的目标输出；

③按照如下公式，循环更新每个粒子的速度和位置，

v_p(iter+1)=ω×v_p(iter)+m₁a₁(Lbest_p-r_p(iter))+m₂a₂(Gbest-r_p(iter))

（13）

r_p(iter+1)=r_p(iter)+v_p(iter+1)（14）式中，v_p表示更新粒子p的速度，r_p表示更新粒子p的位置，Lbest_p表示更新粒子p的个体最优值，Gbest即为整个粒子群的全局最优值，iter表示循环次数，ω是粒子群算法中的惯性权重，m₁、m₂是对应的加速系数，a₁、a₂是[0,1]之间的随机数；

④对于粒子p，如果新的适应度大于原来的个体最优值，更新粒子的个体最优值：

Lbest_p=f_p （15）

⑤如果粒子p的个体最优值Lbest_p大于原来的粒子群全局最优值Gbest，更新原来的粒子群全局最优值Gbest:

Gbest=Lbest_p （16）

⑥判断是否满足性能要求，若是，结束寻优，得到一组优化的模糊方程的局部方程参数；否则返回步骤③，继续迭代寻优，直至达到最大迭代次数iter_max。

迭代终止时的Gbest即为表示第i个标准化后的训练样本X_i的COD预报值和使COD排放达标的操作变量值。

所述的方法还包括：5）、按设定的采样时间间隔，采集现场智能仪表信号，将得到的实测COD与函数预报值比较，如果相对误差大于10%或实测COD数据不达标，则将DCS数据库中生产正常时的达标的新数据加入训练样本数据，更新模型。

6）、在所述的步骤4）中，将COD预报值和使COD排放达标的操作变量传给DCS系统并在DCS的控制站显示过程状态，同时通过DCS系统和现场总线将过程状态信息传递到现 场操作站进行显示；同时，DCS系统将所得到的使COD排放达标的操作变量值作为新的操作变量设定值，自动执行COD排放达标操作。

所述的关键变量包括进入焚烧炉的废液流量，进入焚烧炉的空气流量，进入焚烧炉的燃料流量。；所述的操作变量包括进入焚烧炉的空气流量和进入焚烧炉的燃料流量。

Claims (2)

1.一种群智加权的农药焚烧炉有害物排放达标控制系统，包括焚烧炉、现场智能仪表、DCS系统、数据接口以及上位机，所述的DCS系统包括控制站和数据库；所述现场智能仪表与DCS系统连接，所述DCS系统与上位机连接，其特征在于：所述的上位机包括：

数据预处理模块，用于将从DCS数据库输入的模型训练样本进行预处理，对训练样本中心化，即减去样本的平均值，然后对其进行标准化：

计算均值： $\stackrel{\&OverBar;}{TX}=\frac{1}{N}\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{TX}}_i---\left(1\right)$

计算方差： ${\mathrm{\&sigma;}}_x^2=\frac{1}{N-1}\underset{i=1}{\overset{N}{\&Sigma;}}({\mathrm{TX}}_i-\stackrel{\&OverBar;}{TX})---\left(2\right)$

标准化： $X=\frac{TX-\stackrel{\&OverBar;}{TX}}{{\mathrm{\&sigma;}}_x}---\left(3\right)$

其中，TX为训练样本，TX_i为第i个训练样本，是从DCS数据库中采集的生产正常时的关键变量、化学耗氧量COD和相应的使COD排放达标时的操作变量的数据，N为训练样本数，为训练样本的均值，X为标准化后的训练样本；σ_x表示训练样本的标准差，σ² _x表示训练样本的方差；

模糊方程模块，对从数据预处理模块传过来的标准化后的训练样本X，进行模糊化；设模糊方程系统中有c^*个模糊群，模糊群k、j的中心分别为v_k、v_j，则第i个标准化后的训练样本X_i对于模糊群k的隶属度μ_ik为：

${\mathrm{\&mu;}}_{ik}={\left(\underset{j=1}{\overset{c^{*}}{\&Sigma;}}{\left(\frac{\left|\right|X_i-v_k\left|\right|}{\left|\right|X_i-v_j\left|\right|}\right)}^{\frac{2}{n-1}}\right)}^{-1}---\left(4\right)$

式中，n为模糊分类过程中需要的分块矩阵指数，取作2,||·||为范数表达式；

使用以上隶属度值或者它的变形以获得新的输入矩阵，对于模糊群k，其输入矩阵变形为：

Φ_ik(X_i,μ_ik)＝[1 func(μ_ik) X_i] (5)

其中func(μ_ik)为隶属度值μ_ik的变形函数，取或exp(μ_ik)，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵；

加权最小二乘支持向量机作为模糊方程系统的局部方程，对每个模糊群进行优化拟合；设模型训练样本的第i个目标输出为O_i，加权最小二乘支持向量机通过变换把拟合问题等价于如下二次规划问题：

$\min R(w,\mathrm{\&xi;})=\frac{1}{2}{w}^{T}w+\frac{1}{2}{\mathrm{\&gamma;\&Sigma;}}_{i=1}^N{\mathrm{\&omega;}}_i{\mathrm{\&xi;}}_i^2---\left(6\right)$

其中，R(w,ξ)是优化问题的目标函数，minR(w,ξ)是优化问题的目标函数的最小值，ξ＝{ξ₁,…,ξ_N}是松弛变量，ξ_k是松弛变量的第k个分量，ξ_i是松弛变量的第i个分量，w是加权最小二乘支持向量机超平面的法向量，b是相应的偏移量，而ω_i,i＝1,…,N和γ分别是加权最小二乘支持向量机的权重和惩罚因子，是加权最小二乘支持向量机松弛变量的第k个分量ξ_k标准差的估计，c₁为常数，这里取2.5，c₂为常数，这里取3，上标T表述矩阵的转置；由(6)(7)(8)式可推导出模糊群k在第i个训练样本TX_i的输出为：

${\hat{y}}_{ik}=\underset{k=1}{\overset{N}{\&Sigma;}}{\mathrm{\&alpha;}}_{k}K<{\mathrm{\&Phi;}}_k(X,\mathrm{\&mu;}),{\mathrm{\&Phi;}}_{ik}(X_i,{\mathrm{\&mu;}}_{ik})>+b_i---\left(9\right)$

其中，为第k个群智加权模糊系统输出层的预测输出，α_k，k＝1,…,N是对应的拉格朗日乘子的第k个分量，K<·>是加权最小二乘支持向量机的核函数，这里K<·>取线性核函数；

由反模糊方法中的重心法得到最后的模糊方程系统的输出：

${\hat{y}}_i=\frac{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{ik}{\hat{y}}_{ik}}{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{ik}}---\left(10\right)$

其中，是模糊方程系统的预测输出；

粒子群算法优化模块，用于采用粒子群算法对模糊方程中加权最小二乘支持向量机局部方程的惩罚因子和误差容限值进行优化，具体步骤如下：

①确定粒子群的优化参数为加权支持向量机局部方程的惩罚因子和误差容限值、粒子群个体数目popsize、最大循环寻优次数iter_max、第p个粒子的初始位置r_p0、初始速度v_p0、局部最优值Lbest_p0以及整个粒子群的全局最优值Gbest；

②设定优化目标函数，将其转换为适应度，对每个局部模糊方程进行评价；通过相应的误差函数计算适应度函数，并认为误差大的粒子适应度小，粒子p的适应度函数表示为：

f_p＝1/(E_p+1) (11)

式中，E_p是模糊方程系统的误差函数，表示为：

$E_p=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\&Sigma;}}{({\hat{y}}_i-O_i)}^2}---\left(12\right)$

③按照如下公式，循环更新每个粒子的速度和位置，

v_p(iter+1)＝ω×v_p(iter)+m₁a₁(Lbest_p-r_p(iter))+m₂a₂(Gbest-r_p(iter)) (13)

r_p(iter+1)＝r_p(iter)+v_p(iter+1) (14)

式中，v_p表示更新粒子p的速度，r_p表示更新粒子p的位置，Lbest_p表示更新粒子p的个体最优值，iter表示循环次数，ω是粒子群算法中的惯性权重，m₁、m₂是对应的加速系数，a₁、a₂是[0,1]之间的随机数；

④对于粒子p，如果新的适应度大于原来的个体最优值，更新粒子的个体最优值：

Lbest_p＝f_p (15)

⑤如果粒子p的个体最优值Lbest_p大于原来的粒子群全局最优值Gbest，更新原来的粒子群全局最优值Gbest:

Gbest＝Lbest_p (16)

⑥判断是否满足性能要求，若是，结束寻优，得到一组优化的模糊方程的局部方程参数；否则返回步骤③，继续迭代寻优，直至达到最大循环寻优次数iter_max；

迭代终止时的Gbest即为表示第i个标准化后的训练样本X_i的COD预报值和使COD排放达标的操作变量值；

所述上位机还包括：

判别模型更新模块，用于按设定的采样时间间隔，采集现场智能仪表信号，将得到的实测化学耗氧量与函数预报值比较，如果相对误差大于10％或实测COD数据不达标，则将DCS数据库中生产正常时的达标的新数据加入训练样本数据，更新模型；结果显示模块，用于将COD预报值和使COD排放达标的操作变量值传给DCS系统，在DCS的控制站显示，并通过DCS系统和现场总线传递到现场操作站进行显示；同时，DCS系统将所得到的使COD排放达标的操作变量值作为新的操作变量设定值，自动执行COD排放达标操作；信号采集模块，用于依照设定的每次采样的时间间隔，从数据库中采集数据；所述的关键变量包括进入焚烧炉的废液流量、进入焚烧炉的空气流量和进入焚烧炉的燃料流量；所述的操作变量包括进入焚烧炉的空气流量和进入焚烧炉的燃料流量。

2.一种用如权利要求1所述的群智加权的农药焚烧炉有害物排放达标控制系统实现的方法，其特征在于：所述控制方法具体实现步骤如下：

1)、对农药焚烧炉有害物排放过程对象，根据工艺分析和操作分析，确定所用的关键变量，从DCS数据库中采集生产正常时所述变量的数据作为训练样本TX的输入矩阵，采集对应的COD和使COD排放达标的操作变量数据作为输出矩阵Y；

2)、将从DCS数据库输入的模型训练样本进行预处理，对训练样本中心化，即减去样本的平均值，然后对其进行标准化，使得其均值为0，方差为1；该处理采用以下算式过程来完成：

2.1)计算均值： $\stackrel{\&OverBar;}{TX}=\frac{1}{N}\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{TX}}_i---\left(1\right)$

2.2)计算方差： ${\mathrm{\&sigma;}}_x^2=\frac{1}{N-1}\underset{i=1}{\overset{N}{\&Sigma;}}({\mathrm{TX}}_i-\stackrel{\&OverBar;}{TX})---\left(2\right)$

2.3)标准化： $X=\frac{TX-\stackrel{\&OverBar;}{TX}}{{\mathrm{\&sigma;}}_x}---\left(3\right)$

其中，TX为训练样本，TX_i为第i个训练样本，是从DCS数据库中采集的生产正常时的关键变量、化学耗氧量COD和相应的使COD排放达标时的操作变量的数据，N为训练样本数，为训练样本的均值，X为标准化后的训练样本；σ_x表示训练样本的标准差，σ² _x表示训练样本的方差；

3)对从数据预处理模块传过来的训练样本，进行模糊化；设模糊方程系统中有c^*个模糊群，模糊群k、j的中心分别为v_k、v_j，则第i个标准化后的训练样本X_i对于模糊群k的隶属度μ_ik为：

${\mathrm{\&mu;}}_{ik}={\left(\underset{j=1}{\overset{c^{*}}{\&Sigma;}}{\left(\frac{\left|\right|X_i-v_k\left|\right|}{\left|\right|X_i-v_j\left|\right|}\right)}^{\frac{2}{n-1}}\right)}^{-1}---\left(4\right)$

式中，n为模糊分类过程中需要的分块矩阵指数，取作2,||·||为范数表达式；

使用以上隶属度值或者它的变形以获得新的输入矩阵，对于模糊群k，其输入矩阵变形为：

Φ_ik(X_i,μ_ik)＝[1 func(μ_ik) X_i] (5)

其中func(μ_ik)为隶属度值μ_ik的变形函数，取或exp(μ_ik)，Φ_ik(X_i,μ_ik)表示第i个输入变量X_i及其模糊群k的隶属度μ_ik所对应的新的输入矩阵；

加权最小二乘支持向量机作为模糊方程系统的局部方程，对每个模糊群进行优化拟合；设模型训练样本的第i个目标输出为O_i，加权最小二乘支持向量机通过变换把拟合问题等价于如下二次规划问题：

$\min R(w,\mathrm{\&xi;})=\frac{1}{2}{w}^{T}w+\frac{1}{2}{\mathrm{\&gamma;\&Sigma;}}_{i=1}^N{\mathrm{\&omega;}}_i{\mathrm{\&xi;}}_i^2---\left(6\right)$

其中，R(w,ξ)是优化问题的目标函数，minR(w,ξ)是优化问题的目标函数的最小值，ξ＝{ξ₁,…,ξ_N}是松弛变量，ξ_k是松弛变量的第k个分量，ξ_i是松弛变量的第i个分量，w是加权最小二乘支持向量机超平面的法向量，b是相应的偏移量，而ω_i,i＝1,…,N和γ分别是加权最小二乘支持向量机的权重和惩罚因子，是加权最小二乘支持向量机松弛变量的第k个分量ξ_k标准差的估计，c₁为常数，这里取2.5，c₂为常数，这里取3，上标T表述矩阵的转置；由(6)(7)(8)式可推导出模糊群k在第i个训练样本TX_i的输出为：

${\hat{y}}_{ik}=\underset{k=1}{\overset{N}{\&Sigma;}}{\mathrm{\&alpha;}}_{k}K<{\mathrm{\&Phi;}}_k(X,\mathrm{\&mu;}),{\mathrm{\&Phi;}}_{ik}(X_i,{\mathrm{\&mu;}}_{ik})>+b_i---\left(9\right)$

其中，为第k个群智加权模糊系统输出层的预测输出，α_k，k＝1,…,N是对应的拉格朗日乘子的第k个分量，K<·>是加权最小二乘支持向量机的核函数，这里K<·>取线性核函数；

由反模糊方法中的重心法得到最后的模糊方程系统的输出：

${\hat{y}}_i=\frac{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{ik}{\hat{y}}_{ik}}{{\mathrm{\&Sigma;}}_{k=1}^{c^{*}}{\mathrm{\&mu;}}_{ik}}---\left(10\right)$

其中，是模糊方程系统的预测输出；

4)采用粒子群算法对模糊方程中加权最小二乘支持向量机局部方程的惩罚因子和误差容限值进行优化，具体步骤如下：

①确定粒子群的优化参数为加权最小二乘支持向量机局部方程的惩罚因子和误差容限值、粒子群个体数目popsize、最大循环寻优次数iter_max、第p个粒子的初始位置r_p0、初始速度v_p0、局部最优值Lbest_p0以及整个粒子群的全局最优值Gbest；

②设定优化目标函数，将其转换为适应度，对每个局部模糊方程进行评价；通过相应的误差函数计算适应度函数，并认为误差大的粒子适应度小，粒子p的适应度函数表示为：

f_p＝1/(E_p+1) (11)

式中，E_p是模糊方程系统的误差函数，表示为：

$E_p=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\&Sigma;}}{({\hat{y}}_i-O_i)}^2}---\left(12\right)$

③按照如下公式，循环更新每个粒子的速度和位置，

v_p(iter+1)＝ω×v_p(iter)+m₁a₁(Lbest_p-r_p(iter))+m₂a₂(Gbest-r_p(iter)) (13)

r_p(iter+1)＝r_p(iter)+v_p(iter+1) (14)

式中，v_p表示更新粒子p的速度，r_p表示更新粒子p的位置，Lbest_p表示更新粒子p的个体最优值，iter表示循环次数，ω是粒子群算法中的惯性权重，m₁、m₂是对应的加速系数，a₁、a₂是[0,1]之间的随机数；

④对于粒子p，如果新的适应度大于原来的个体最优值，更新粒子的个体最优值：

Lbest_p＝f_p (15)

⑤如果粒子p的个体最优值Lbest_p大于原来的粒子群全局最优值Gbest，更新原来的粒子群全局最优值Gbest:

Gbest＝Lbest_p (16)

⑥判断是否满足性能要求，若是，结束寻优，得到一组优化的模糊方程的局部方程参数；否则返回步骤③，继续迭代寻优，直至达到最大循环寻优次数iter_max；

迭代终止时的Gbest即为表示第i个标准化后的训练样本X_i的COD预报值和使COD排放达标的操作变量值；

所述方法还包括以下步骤：5)、判别模型更新模块，用于按设定的采样时间间隔，采集现场智能仪表信号，将得到的实测化学耗氧量与函数预报值比较，如果相对误差大于10％或实测COD数据不达标，则将DCS数据库中生产正常时的达标的新数据加入训练样本数据，更新模型；

6)、在所述的步骤4)中得到的COD预报值和使COD排放达标的操作变量值，将结果传给DCS系统，在DCS的控制站显示，并通过DCS系统和现场总线传递到现场操作站进行显示；同时，DCS系统将所得到的使COD排放达标的操作变量值作为新的操作变量设定值，自动执行COD排放达标操作；

所述的关键变量包括进入焚烧炉的废液流量、进入焚烧炉的空气流量和进入焚烧炉的燃料流量；所述的操作变量包括进入焚烧炉的空气流量和进入焚烧炉的燃料流量。