技术领域
[0001] 本
发明涉及一种引力搜索算法及基于该算法的
污水处理控制方法。
背景技术
[0002] 只有在污水处理工艺中引入优质的计算机技术、当代
优化算法和控制理论才能提升污水处理的效率，并减少处理成本。万有引力搜索算法(Gravitational search algorithm，GSA)是由伊朗克曼大学的Esmat Rashedi等人于2009年所分析的一种新的启蒙式优化算法。此种GSA算法的实质是模仿物理学中的万有引力现象，并将其演算成随机寻找最优解的过程。可是GSA目前使用较少，对自适应的研究更是寥寥无几。
[0003] 虽然GSA算法是全局收敛的算法，能够寻找到最优值，但是如果算法在收敛之前没有搜索到全局(或者局部)最优解，将导致过早收敛而找不到最优值，也就是说，GSA跟其他优化算法一样，容易陷入局部最优。
[0004] 现在我国已经建设了很多的污水处理厂，只要目的是减少水污染对环境形成的消极影响，当然处理过程的自动控制方法很前沿，但至今仍没有实现最优控制。只有污水处理已经到达《城镇污水处理厂污染物排放标准》(GB18918-2002)，然后才可以思量控制的目标。从环境保护的
角度研究，出水水质指标越小越好，而从运行成本研究，出水水质越好意味着消耗更多的资源和更多的运营
费用。宁夏某污水处理厂使用先进的污水处理设施，采纳改良型SBR技术。
[0005] 改良型SBR工艺还存在两个主要问题：一是污水处理中的水质参数检测滞后，采用的是离线人工检测方式，一般是人工取水，然后人工在实验室化验，得出出水水质指标COD(
化学需氧量)、
氨氮HN3-N等的化验结果，根据水质情况调整鼓
风曝气盘。这样的检测方式需要耗时半小时，明显的对水质控制不够及时，效率低下。二是只有在水质达标的前提下，较少工艺的运行费用，才能够增强其经济效益。该工艺的运行费用主要包括人工费、
污泥处理费、鼓风曝气费这三部分，其中人工费和污泥处理费一般根据工厂规模可以固定，而鼓风曝气费消耗大，且难控制。鼓风曝气量需要根据污泥
排放量和溶解氧的浓度进行实时快速的调整，才能是曝气量最少，从而费用最优。
发明内容
[0006] 本发明所要解决的技术问题是提供一种有效控制算法收敛速度，同时提高算法
精度的差分自适应万有引力搜索算法。
[0007] 本发明解决上述技术问题的技术方案如下：一种差分自适应万有引力搜索算法，包括如下步骤：
[0008] S1、初始化种群；
[0009] S2、计算种群适应度值；
[0010] S3、计算种群的惯性
质量M及
加速度a；
[0012] S5、比对是否达到
迭代次数，当达到迭代次数时，输出最优值；当未达到迭代次数时，重复所述步骤S2、S3和S4。
[0013] 本发明的有益效果是：在加快收敛速度，收敛速度的跟步长有关，现在最好的方法通过实时更新个体参数G0和加速度α，有效控制算法收敛速度，同时提高算法精度。
[0014] 在上述技术方案的
基础上，本发明还可以做如下改进。
[0015] 进一步，所述步骤S4包括将万有引力搜索算法的公式带入更新公式中；
[0016] 所述万有引力搜索算法的公式包括：
[0017] 粒子间的引力为：
[0018]
[0019] 式中，G表示引力常量，M1、M2分别表示两种粒子的惯性质量，R表示两种粒子中间的欧式距离，F表示引力；
[0020] 依据
牛顿第二定律，当某个粒子受到其他粒子作用力时，随即会产生F大小的合力，同时会产生向这个方向运动的加速度a，即
[0021]
[0022] 引力常数G跟随着
宇宙实际年龄变化的趋势为
[0023]
[0024] 式中，G(t)表示在t时间时的引力常数G的数值，G(t0)表示在第一宇宙子间隔时刻t0的数值；
[0025] 所述更新公式包括：
[0026] 有吸引力就会产生相应的加速度，则物质i在第k维上的加速度为[0027]
[0028] 每次迭代更新的公式为
[0029]
[0030] 本发明所要解决的技术问题是提供一种能预测出水水质的污水处理控制方法。
[0031] 本发明解决上述技术问题的技术方案如下：一种污水处理控制方法，包括如下步骤:
[0032] S100、设定该问题的搜索空间维度为2，设置如下的空间上下限范围：
[0033]
[0034] 其中， 是氨氮含量；COD是化学需氧量；ρDO溶氧量；Qw是污泥排放量。
[0035] S200、种群初始化；
[0036] S300、初始化全局最优值gbest；
[0037] S400、计算适应值和加速度；
[0038] S500、更新速度与位置；
[0039] S600、更新全局最优粒子与个体最优粒子，得到全局最优解pbesti和gbest；
[0040] S700、比对是否达到迭代次数，当达到迭代次数时，输出最优解，当未达到迭代次数时，重复所述步骤S300、S400、S500和S600；
[0041] S800、通过输出的最优值，找到函数值最小时方程所对应的解，即溶解氧浓度ρDO优化设定值和污泥排放量Qw的优化设定值。
[0042] 本发明的有益效果是，通过输出的最优值，找到函数值最小时候，此时方程所对应的解，这就是得到了优AGSA算法寻求得到的溶解氧浓度ρDO优化设定值和污泥排放量Qw的优化设定值。该方法结合差分自适应引力搜索算法，使得鼓风曝气量根据进水口量和溶解氧的浓度进行实时快速的调整，使曝气量最少，从而费用最优。
[0043] 所述步骤S200包括：设置万有引力场的粒子数目100，最大迭代步数M＝200；根据电脑随机初始化种群中各个粒子的位置X、速度Y以及速度的上下限固定范围[Vmin,Vmax]，每个粒子个体历史最优位置pbesti赋值：对于每个粒子，初始时都有相应的赋值pbesti。
[0044] 所述步骤S300包括：在规定的目标要求下，依次算出各个初始粒子X的目标函数值f1,f2，电脑随机取其中某个粒子做为初始的gbest。
[0045] 所述步骤S400包括：
[0046] S401、计算其目标函数值，而后根据线性权重的方法，把原先设定好的目标适应度值统一成一个适应度值：
[0047]
[0048] 式中的t是当前的迭代次数；
[0049] S402、通过下式计算加速度：
[0050]
[0051]
[0052] 所述步骤S401中目标函数为：
[0053]
[0054] 目标函数中T为优化周期，P为鼓风机的功率，a为每度电的费用，V是污水处理厂体积，c1＝0.455为每千克去除COD(化学需氧量)的电耗，c2＝9.1为每千克去除NH3-H(氨氮含量)的电耗，系数c1V,c2V可以记为电费的奖惩因子。
[0055] 所述步骤S401中的目标方程可以简化为：
[0056]
[0057] 所述步骤S500包括，通过下式更新速度与位置：
[0058]
[0059] 其中， 是当前位置； 前一时刻位置； 前一时刻速度；Rij粒子间的欧式距离； 是惯性质量。
附图说明
[0060] 图1为本发明
实施例2中功率P与COD的关系图；
[0061] 图2为本发明实施例2中功率P与氨氮含量的关系图；
[0062] 图3为本发明实施例2中COD与时间的关系图；
[0063] 图4为本发明实施例2中氨氮含量与时间的关系图；
[0064] 图5为本发明实施例2中某次运行费用优化曲线。
具体实施方式
[0065] 以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
[0066] 实施例1
[0067] 一种差分自适应万有引力搜索算法，包括如下步骤：
[0068] S1、初始化种群；
[0069] S2、计算种群适应度值；
[0070] S3、计算种群的惯性质量M及加速度a；
[0071] S4、更新个体位置；
[0072] S5、比对是否达到迭代次数，当达到迭代次数时，输出最优值；当未达到迭代次数时，重复所述步骤S2、S3和S4。
[0073] 所述步骤S4包括将万有引力搜索算法的公式带入更新公式中；
[0074] 所述万有引力搜索算法的公式包括：
[0075] 粒子间的引力为：
[0076]
[0077] 式中，G表示引力常量，M1、M2分别表示两种粒子的惯性质量，R表示两种粒子中间的欧式距离，F表示引力；
[0078] 依据牛顿第二定律，当某个粒子受到其他粒子作用力时，随即会产生F大小的合力，同时会产生向这个方向运动的加速度a，即
[0079]
[0080] 引力常数G跟随着宇宙实际年龄变化的趋势为
[0081]
[0082] 式中，G(t)表示在t时间时的引力常数G的数值，G(t0)表示在第一宇宙子间隔时刻t0的数值。
[0083] 惯性质量计算：
[0084] 假使在一个D维的搜寻空间中，包含了N个物质，第i个物质的具体位置为：
[0085]
[0086] 式中 表示第i个物体在第k维上的位置。引力搜索算法中，每个粒子的惯性质量和它所处的位置，都与它的适度值息息相关，在t时间时粒子Xi的质量可以表达为Mi(t)。
[0087] Mi(t)根据上式进行计算得到
[0088]
[0089] 式中fiti(t)代表粒子Xi适应度。best(t)是时刻t的最优解，worst(t)表示时刻t的最差解，依次根据上面的式子计算得到
[0090]
[0091] 式中，mi(t)将粒子适应值规范化到[0,1]之间，根据其占总质量中的比重当做粒子质量Mi(t)。
[0092] 空间引力计算：
[0093] t时刻物体j在第k维上受到物体i的引力为
[0094]
[0095] 式中，ε是一个非常小的常数，Mat(t)是产生在物质j上的惯性质量，G(t)是万有引力常数，表示为
[0096]
[0097] Rij(t)为两粒子之间的欧式距离
[0098] Rij(t)＝||Xi(t),Xj(t)||2
[0099] 所以t时刻，第k维上作用在Xi的力的总和为式
[0100]
[0101] 所述更新公式包括：
[0102] 有吸引力就会产生相应的加速度，则物质i在第k维上的加速度为[0103]
[0104] 每次迭代更新的公式为
[0105]
[0106] 实施例2
[0107] 一种污水处理控制方法，包括如下步骤:
[0108] S100、设定该问题的搜索空间维度为2，设置如下的空间上下限范围：
[0109]
[0110] 其中， 是氨氮含量；COD是化学需氧量；ρDO溶氧量；Qw是污泥排[0111] 放量；
[0112] S200、种群初始化；
[0113] S300、初始化全局最优值gbest；
[0114] S400、计算适应值和加速度；
[0115] S500、更新速度与位置；
[0116] S600、更新全局最优粒子与个体最优粒子，得到全局最优解pbesti和gbest；
[0117] S700、比对是否达到迭代次数，当达到迭代次数时，输出最优解，当未达到迭代次数时，重复所述步骤S300、S400、S500和S600；
[0118] S800、通过输出的最优值，找到函数值最小时方程所对应的解，即溶解氧浓度ρDO优化设定值和污泥排放量Qw的优化设定值。
[0119] 该工艺的运行费用主要包括人工费、污泥处理费、鼓风曝气费这三部分，其中人工费和污泥处理费一般根据工厂规模可以固定，而鼓风曝气费由污泥排放量和溶解氧浓度影响。
[0120] 该污水处理厂体积为V＝L×B×H＝42×22×6.6≈6.1×103m3＝6.1×106L，去除COD(化学需氧量)和NH3-H(氨氮含量)所需的电耗分别为0.455kWgh/kg和9.1kWgh/kg。
[0121] 设定污泥排放量Qw和溶解氧浓度ρDO为控制变量，COD(化学需氧量)和SNH3(氨氮含量)为状态变量，将鼓风曝气费运最低设置为本文的目标函数，，建立目标方程。
[0122]
[0123] 目标函数中T为优化周期，P为鼓风机的功率，a为每度电的费用，V是污水处理厂体积，c1＝0.455为每千克去除COD(化学需氧量)的电耗，c2＝9.1为每千克去除NH3-H(氨氮含量)的电耗，系数c1V,c2V可以记为电费的奖惩因子。
[0124] 化简目标方程：
[0125]
[0126] 分析上面的目标函数可知：1)进水时COD(化学需氧量)很大，需要大量鼓风曝气，则鼓风机功率可在额定功率下工作，为的是给处理池中好氧菌提供大量氧气，使得好氧菌吸收污染物，从而快速降低COD。2)随着COD量的降低，当降低到60时，达到国家标准。因为考虑政府政策，COD含量达到60mg/l时，政府补偿0.6元/吨；COD含量降到10mg/l时，政府补偿0.8元/吨。3)当COD的量降低到国家以下时，若持续降低COD的量，鼓风机则会因为时间过长，而导致鼓风机曝气时间加长，能耗增加。所以要平衡三者之间的关系，使得在水质达标的情况下，能耗最小，即费用最小。
[0127] 鼓风机的作用是吹氧气，它的的功率P与出水水质化学需氧量COD有关联。通过现场采集的数据P和COD(化学需氧量)的量，进行拟合，拟合图像如图1和图2所示。图1显示P与COD的关系成线性，且斜率k1小于0。COD(化学需氧量)与时间的关系如下图3和图4所示。
[0128] NH3-H(氨氮含量)跟COD(化学需氧量)具有同样的原理。在理想状态下，
状态方程可以写成：
[0129]
[0130] 式中k——底物最大比例速度常数；
[0131] KO——氧的饱和常数；
[0132] Ksub——底物饱和系数；
[0134] Y——产率系数。
[0135] 所述步骤S200包括：设置万有引力场的粒子数目100，最大迭代步数M＝200；根据电脑随机初始化种群中各个粒子的位置X、速度Y以及速度的上下限固定范围[Vmin,Vmax]，每个粒子个体历史最优位置pbesti赋值：对于每个粒子，初始时都有相应的赋值pbesti。
[0136] 所述步骤S300包括：在规定的目标要求下，依次算出各个初始粒子X的目标函数值f1,f2，电脑随机取其中某个粒子做为初始的gbest。
[0137] 所述步骤S400包括：
[0138] S401、计算其目标函数值，而后根据线性权重的方法，把原先设定好的目标适应度值统一成一个适应度值：
[0139]
[0140] 式中的t是当前的迭代次数；
[0141] S402、通过下式计算加速度：
[0142]
[0143]
[0144] 所述步骤S401中目标函数为：
[0145]
[0146] 目标函数中T为优化周期，P为鼓风机的功率，a为每度电的费用，V是污水处理厂体积，c1＝0.455为每千克去除COD(化学需氧量)的电耗，c2＝9.1为每千克去除NH3-H(氨氮含量)的电耗，系数c1V,c2V可以记为电费的奖惩因子。
[0147] 所述步骤S401中的目标方程可以简化为：
[0148]
[0149] 所述步骤S500包括，通过下式更新速度与位置：
[0150]
[0151] 其中， 是当前位置； 前一时刻位置； 前一时刻速度；Rij粒子间的欧式距离； 是惯性质量。
[0152] 该实施例中，依托宁夏某污水处理厂，通过调查、实地考察、实时采集数据，在改良型SBR工艺流程进行仿真实验。该厂一个SBR工艺周期大约4.8h，进水COD(化学需氧量)的平均值为175mg/L，出水口的COD(化学需氧量)约是20mg/L，工厂内某曝气池可利用的池体积为V＝L×B×H＝42×22×6.6≈6.1×103m3。处理周围12km2的地区，地区有居住人口13万人。这个污水处理厂的出水执行《城镇污水处理厂污染物排放标准》(GB18918-2002)中一级B的标准。
[0153] 在MATLAB R2010a测试环境下，利用论文中改进的自适应引力搜索(AGSA)与常规遗传GA和GSA计算、比较和分析。其中参数选取：最大迭代步数200，为了避免随机性的干扰，3种算法独立运算20次，图5为其中一次优化过程。
[0154] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
