基于模拟退火粒子群算法的水资源优化配置

杜佰林，张建丰，高泽海，李涛，黄子奇，张娜

(西安理工大学省部共建西北旱区生态水利国家重点实验室，陕西 西安 710048)

当前，随着中国区域经济的快速发展和城市化进程的加快，水资源开发不平衡、地下水开采量大、水资源利用效率低等问题日趋严峻，造成水资源短缺现象逐渐凸显[1].同时，区域降雨量时空分布不均、用水管理水平不高，导致水资源供需矛盾有所恶化，制约了社会经济快速发展[2].为了缓解上述危机，研究区域水资源优化配置显得尤为必要，而传统的配置方法面对水资源系统的复杂性和不确定性时，往往只进行经验估计或者忽略，致使配置结果存在偏差.因此，研究如何实现区域水资源可持续开发利用和水资源优化配置，具有十分重要的现实意义.

国内外众多学者针对水资源优化配置中模型的建立开展了大量研究工作.国际上利用计算机技术对水资源系统的复杂关系进行了模拟仿真，逐步成为一种被普遍认可的配置思路[3].MASSE[4]针对水资源配置中的单目标问题进行了研究，该成果的特点是将配置理论应用到了实践中.ABDULBAKI等[5]提出并解决了水资源配置的多目标问题.WADA等[6]建立了不同水质下的区域多水源、多用水户水资源配置模型.EL-NAQA等[7]在综合考虑地下水超采、水质恶化的基础上建立了多水源联合调度模型.

中国的水资源分配研究工作以水库优化调度为基点，经历了水资源的供需平衡配置、经济效益最大化配置、面向生态效益合理配置这3个阶段.粟晓玲等[8]提出了干旱地区补给水的同时考虑生态效益的配置供需理念.刘玒玒等[9]应用蚁群和粒子群的混合智能优化算法解决了同时考虑水质和水量的多目标水资源优化配置模型.张倩等[10]应用改进的粒子群算法解决了当前灌区有限农业水资源在不同作物中合理分配的问题.

可以发现，水资源优化配置模型由于自身的多元化呈指数增长，简单的线性规划求解已不能满足配置模型的发展趋势，甚至一些智能优化算法，例如遗传算法、蚁群算法、粒子群算法等往往由于自身的局限，同样导致配置结果不理想.因此，在水资源优化配置过程中选择适宜的模型，克服算法自身寻优的缺陷，提高算法收敛精度及效率，是目前配置过程中亟待解决的问题.

鉴于此，文中建立综合效益相对最优的大荔县水资源优化配置模型，以实现该区域水资源的可持续开发利用；通过对粒子群算法(particle swarm optimization, PSO)和模拟退火算法(simulated annea-ling, SA)深入分析，将模拟退火算法在一定概率控制下，暂时接受一些劣质解的特性，用于改善基本粒子群算法的局部寻优能力，从而获得具有较好全局寻优能力的模拟退火粒子群算法(simulated annealing-particle swarm optimization, SA-PSO)，以期为今后水资源规划模型提供新的求解方法，为区域水资源管理者提供更多的决策依据.

1 模型的建立与求解算法

1.1 目标函数

水资源优化配置是以公平性、持续性和高效性为原则，在一些特定的区域内运用一些技术方法，将有限的水资源科学合理地分配给各用水部门，以满足社会、经济和环境的可持续发展，寻求综合效益的最优值.因此构建以社会、经济、生态效益为目标的综合评价函数，表达式为

Z=opt[Fecon(X),Fsoci(X),Fecol(X)]，

(1)

式中：X为不同赋存形式、不同质量和数量的水资源所构成的决策变量；Fecon(X)，Fsoci(X),Fecol(X)分别为本模型求解的经济效益、社会效益和生态效益.

1) 经济效益目标.以水资源开发利用中不同水平年下各用水户的直接经济效益最大化作为经济效益目标，表达式为

(2)

式中：i为供水水源类别；j为用水户类别；Qij为水源i向用水户j的供水量，万m3；bij为水源i向用水户j的效益系数，元/m3；cij为水源i向用水户j的费用系数，元/m3；αi为水源i的供水次序系数；βj为用水户j的用水公平系数.

2) 社会效益目标.考虑到社会效益很难量化，而社会的发展和稳定受区域缺水程度影响较大，且区域缺水率又可以度量区域的缺水程度，故以水资源开发利用中区域综合缺水率最小化作为社会效益目标，表达式为

(3)

式中：Dj为用水户j的需水量，万m3.

3) 生态效益目标.以各用水户在对应规划水平年中所排放的重要污染因子(chemical oxygen demand, COD)的总量最小化作为生态效益目标，表达式为

(4)

式中：dj为用水户j排放单位体积废水中所含重要污染因子的量，mg/L；pij为用水户j的污水排放系数.

1.2 约束条件

文中研究目的旨在配置各用水户在不同水源地下的供水量.该供水量在满足配置目标的同时，仍受供水水源、需水能力、环境能力及各变量非负等约束.约束条件表示如下.

1) 供水能力约束.水源i向用水户j的供水总量应小于其最大可供水量，表达式为

(5)

式中：Wimax为要求水源i的最大可供水量，元/m3.

2) 需水能力约束.配水过程中各用水户从水源i所得到的水量应该结合自身需水能力，既不小于最低约束也不得超过额定约束，表达式为

(6)

3) 环境能力约束.各用水户所排放单位体积废水中所含重要污染因子的含量应在国家允许的排放指标内；且排放的重要污染因子总量应不超过该区域的最大允许排放量，表达式为

(7)

式中：m0j为在各行业中国家标准规定所排放的重要污染物COD质量浓度，mg/L；Mmax为区域内最大允许排放重要污染物的总量.

4) 变量非负约束.水源i向用水户j的供水量要满足非负约束，表达式为

Qij≥0,i=1,2,3,…;j=1,2,3,….

(8)

1.3 模拟退火粒子群算法

1.3.1 模拟退火算法

模拟退火算法是由METROPOLIS等[11]根据物理学中固体物质的冷却退火过程和一般组合优化问题求解过程之间的相似性，提出的一种随机全局优化算法.该求解步骤如下.

Step1 初始一个足够大的温度T0>0，降温次数k=0，随初始化x(0)，使得x(i)=x(0)，并计算其能量值E(x(0)).

Step3 根据降温公式d(Tk)，计算Tk+1=d(Tk),k=k+1，判断是否满足Metropolis准则，若满足则计算结束，输出结果；否则返回Step2.

1.3.2 粒子群算法

粒子群算法是由KENNDY等[12]受到鸟群觅食行为的启发后，提出的一种涉及参数少、易实现的高效搜索算法，与模拟退火算法相同属于进化算法，均以迭代的方式探寻最优解.该求解步骤如下.

Step1 初始化含有种群规模N，初始位置popi和速度vi的粒子群体.

Step2 计算每个粒子的适应度值fitness(i).

Step3 对每个粒子，比较个体极值Pbest(i)与适应度值fitness(i)，若fitness(i)>Pbest(i)，则令Pbest(i)=fitness(i).

Step4 对每个粒子，比较全局极值gbest(i)与适应度值fitness(i)，若fitness(i)>gbest(i),则令gbest(i)=fitness(i).

Step5 更新粒子的位置popi和速度vi，即

vid(t+1)=wvid(t)+c1r1[pid-popid(t)]+c2r2[pgd-popid(t)]，

(9)

popid(t+1)=popid(t)+vid(t+1)，

(10)

式中：w为惯性权重；c1，c2为学习因子，分别代表i的个体认知和社会加速常数；r1，r2为[0,1]相互独立的随机数.

Step6 判断是否满足终止条件，若符合其精度要求或达到最大迭代次数，则计算结束，输出结果；否则返回Step2.

1.3.3 模拟退火粒子群算法

粒子群算法虽然能够有效地优化各种函数，但对离散函数、多目标优化等问题，往往容易陷入局部极小值；模拟退火算法则可通过在一定概率下接受一些劣质解，从而跳出原状态并继续搜索.文中在求解过程中将两者融合，构建模拟退火粒子群算法，使算法在搜索过程中有较强的全局寻优能力，有效避免搜索陷入局部极小解.其求解步骤如下.

Step1 初始化含有种群规模N，初始位置popi和速度vi的粒子群体.

Step2 计算每个粒子的适应度fitness(i)，比较个体极值Pbest(i)与适应度值fitness(i)，若fitness(i)>Pbest(i),则令Pbest(i)=fitness(i)；同理比较全局极值gbest(i)与适应度值fitness(i)，若fitness(i)>gbest(i),则令gbest(i)=fitness(i).

Step3 设定足够大的初始温度T0>0并进行相应的退温操作，即

T0=f(pi)/ln 5，

(11)

Tk+1=λTk，

(12)

式中：λ为退火过程中的衰减参数.

Step4 确定当前温度T下每个粒子pi的适应值，即

(13)

式中：TF为粒子适应值.

Step5 从所有pi中确定其全局最优的代替值p′i，并采用式(9)和(10)更新其位置和速度.

Step6 计算确定粒子目标值，并更新Pbest(i)和gbest(i)，随后进行退温操作.

Step7 判断是否满足终止条件，若符合条件，则计算结束并输出结果；否则返回Step4.

该算法流程如图1所示.

图1 模拟退火粒子群算法计算流程

2 实例研究

2.1 研究区域概况

大荔县(109°43′～110°19′E，34°36′～35°02′N)位于陕西省关中东部平原区，地处黄、洛、渭河的交汇带，并以“三秦通衢，三辅重镇”而出名，是一座典型的陕晋豫承接产业转移示范区，总县域面积为1 776 km2，属暖温带半干旱大陆性季风气候.全县年平均气温为14.4 ℃，年均降水量为514 mm，年均蒸发量为968.3 mm，降水量远小于蒸发量，很难满足蒸发量需求，时有水资源短缺现象发生.大荔县水资源供给系统依据《全国水资源综合规划》概要[13]的有关要求，具体描述如图2所示.

图2 大荔县水资源供给系统

2.2 数据来源

考虑到配置的合理性，在2016现状年的基础上选择在规划水平年2020和2030年来水保证率P={50%，75%，95%} 的3种不同来水水平情况下进行水量优化配置.统计大荔县多年外调水量资料和渭南市“十三五”水利发展规划得出县域目前水资源总量为67 938.67万m3.结合大荔县现状库坝工程、调水工程及规划新建工程的设计用途、供水对象、工程布局和调整方案等实际情况，将可供水量按地表水、地下水、其他水等3个供水水源进行统计，即得到不同水平年不同来水保证率下的可供水量Qor；根据1996—2016年的《大荔县统计年鉴》、《大荔县志》、《洛惠渠志》等相关资料，在获得区域社会经济发展指标的基础上，采用水资源规划中定额法、趋势法进行预测，从而得到在2020和2030年中P={50%，75%，95%}保证率时，不同用水行业的需水量Qre具体数值见表1.

表1 大荔县不同水平年不同保证率下需水量预测

2.3 模型参数确定

2.3.1 模型参数

依据上述水资源优化配置模型并结合大荔县实际用水情况，确定本次水资源配置的供水次序：地表水、地下水、其他水；参考式(14)将供水次序转化为0～1的系数，即供水次序系数αi分别为0.50，0.33，0.17.

(14)

式中：I为供水水源的总数；hi为水源i的供水次序序号；hmax为其供水序号的最大值.

配置过程以“先生活、后生产”为前提条件，在坚持以生活和生态用水为优先原则的基础上，设定需水量部门得到用水的顺序依次为生活、生态、工业、第三产业和农业用水.参数设置参考式(14)，即各部门用水公平系数βj分别为0.33，0.27，0.20，0.13，0.07.用水效益系数bij参考《陕西省行业用水定额》(DB 61/T 943—2013)得出，即生活、生态、工业、第三产业和农业的bij分别为500，500，580，620，2元/m3.据大荔县现状年2016年的水费征收标准确定各部门的费用系数cij，即生活、生态、工业、第三产业和农业的cij分别为2.58，2.58，3.98，5.76，0.45元/m3.需水系数rj的设置参照文献[14]，采用概率分布函数和置信区间，设置置信水平为0.90时，rj～N(0.85,0.382)，即生活和生态的rj为1，工业、第三产业和农业的rj均服从rj～N(0.85,0.382).大荔县污水排放主要表现为城镇生活和工业废水排放，故污水排放系数pij统一参照《城市排水工程规划规范》(DB 50318—2000)为0.75.目标函数权重的设置，结合水资源系统中不确定性的特点，在采用概率分布的基础上参照文献[15]中熵权法的相关知识，计算得社会、经济、生态效益的目标函数权重依次为0.5，0.3，0.2.

2.3.2 算法参数

以Matlab R2018b为开发工具编写模拟退火粒子群算法程序，依据前人经验[16]并经过多次试算求解，设定粒子种群数目N为1 000，惯性权重ω为0.729 8，学习因子c1=c2=2，r1和r2均为[0，1]的随机数；设定模拟退火初始温度T0为1 000，终止温度Tf为1，衰减参数λ为0.6，最大迭代次数maxTit为1 000.

3 结果与分析

3.1 最优配水量

利用所建立的基于模拟退火粒子群算法的优化配置模型对大荔县进行水资源优化配置，得出大荔县在不同规划水平年、不同来水保证率下的水资源优化配置结果Qop，结合水资源供需平衡特征，对水资源配置结果进行分析，以P=75%为例，分析结果见表2，表中Qor，Qre，rws分别为供水量、需水量、缺水率.

由表2可知，到2020年和2030年，在“优先保证生活和生态用水”的前提下，大荔县生活和生态用水均得到了完全满足，工业用水缺水情况极少，缺水情况主要表现在农业用水和第三产业用水，符合县域现状.优化配置结果表明，2020年大荔县农业缺水率为3.24%，第三产业缺水率为0.07%，整体缺水率为2.00%；2030年大荔县的农业缺水率为3.57%，第三产业缺水率为0.13%，整体缺水率为1.67%，缺水情况仍然存在.但以上缺水率均未超过4%，可通过采用一定农艺措施或调整农业种植结构，促使该行业用水供需平衡.配置结果符合规划模型中约束条件的设置，满足5类不同用水户的用水需求，配置结果合理.同时，研究区整体缺水率降低，一定程度上提高了各行业用水需求，这也表明基于模拟退火粒子群算法的优化配置模型具有可行性.

表2 大荔县不同水平年P=75%保证率下的水资源优化配置结果

3.2 配置结果

为了说明优化配置后的水量Qop与原始供水量Qor之间的差异，在模型求解得到地表水、地下水和其他水源在不同规划水平年、不同来水保证率下配置水量的基础上，将大荔县不同水平年、不同保证率下供水量对比结果进行分析，具体分析结果如表3所示，表中Q为原始供水量与优化配置后水量的差值.

表3 大荔县不同水平年不同保证率下供水量对比结果

由表3可知，优化配置后的水量较原始供水量均有递减的趋势，其中2020年在P=50%，75%和95%的3种不同来水水平情况下，优化配置后水量较原始供水量的总节水量分别为372.66，333.11，255.14万m3；2030年在P=50%，75%和95%的3种不同来水水平情况下，优化配置后水量较原始供水量的总节水量分别为426.63，404.76和352.85万m3.总节水量均超过250万m3，节水效果明显，达到了节约用水的目的，可满足渭南市水资源综合规划要求.

对大荔县进行水资源优化配置后的3个目标函数效益值进行计算和分析，以P=75%为例，分析结果见表4.与2016现状年相比，大荔县在2020和2030年2个规划水平年中的经济效益fecon(指各用水户的直接经济效益)、社会效益fsoci(以区域综合缺水率表示，即各用水户的缺水量/需水量的百分数)和生态效益fecol(指各用水户所排放的重要污染因子COD的总量)改善均比较明显，配置结果合理，验证了基于模拟退火粒子群算法优化配置模型具有可行性.

表4 大荔县不同水平年P=75%保证率下的效益值

3.3 模型精度

为了检验模拟退火粒子群算法的可靠性，引入经典测试函数——Ackley函数[17]：

-8≤xi≤8.

(15)

以Matlab R2018b为开发工具，编写Ackley函数程序并绘制三维图像，获取Ackley函数的空间几何特性.图3为Ackley函数图.该函数较为复杂，是一个n维函数并有多个局部极小值点，且fmin(0,0)=0，可用于检测一个算法的全局收敛速度.分别采用SA，PSO和SA-PSO 的3种进化算法依次迭代计算测试函数，从而评价比较3种进化算法的性能，整理对比后结果见表5，表中Titf和Tits分别为最快、最慢迭代次数.由表可知3种进化算法都能在有限的Tit内找到其最优解，其中SA的收敛速度最慢，为327次；SA-PSO的收敛速度最快，仅为15次，优化效率显著提升.

图3 Ackley函数图

表5 3种进化算法的性能对比

为了进一步表明模拟退火粒子群算法全局寻优的优越性，在Matlab R2018b中对性能差异较小的基本PSO和SA-PSO的个体最优适应值Fva和迭代次数Tit进行了仿真, 仿真结果如图4所示.

根据图4可以得出：SA-PSO的性能较基本PSO有了一定程度的提高，在初始条件基本相同的情况下，SA-PSO由于结合了SA思想中一定概率控制下暂时接受一些劣质解的特性，克服了基本PSO容易陷入局部最优解的缺点，使其全局寻优能力增强，进而在精度和进化速度两方面都要优于基本PSO，且SA-PSO原理通俗易懂、容易实现，可为今后水资源规划模型提供一个新的求解思路.

图4 大荔县2020，2030年优化配置的进化过程

4 结 论

构建以社会、经济、生态效益为目标的综合评价函数，建立基于模拟退火粒子群算法的水资源优化配置模型；结合大荔县进行实例验证，证明算法的可行性，实现大荔县的不同行业用水量的合理优化配置.

1) 通过将模拟退火算法在一定概率控制下，暂时接受一些劣质解的特性，用于改善基本粒子群算法的局部寻优能力，有效克服了单一粒子群算法和模拟退火算法相应的缺点，避免粒子陷入局部极小值，提高了算法的精度，加快了算法的进化速度.

2) 构建大荔县水资源优化配置模型，并采用模拟退火粒子群算法进行求解，通过粒子编码、设置约束条件等，求解模型得到大荔县不同行业用水量的优化配置方案.优化后节水效果明显，经济效益增长显著，可满足大荔县的可持续发展要求.

3) 通过对比3种进化算法计算测试函数的结果，比较个体最优适应值的求解结果，发现模拟退火粒子群算法能够获得更好的个体适应度值，收敛速度更快，整体运行更稳定，为区域水资源优化配置提供了一种新的解决途径.