基于SA-FOA的水库防洪优化调度研究

王成民，王 昊，刁艳芳

（1.日照水库管理运行中心，山东 日照 276800； 2.山东农业大学水利土木工程学院，山东 泰安 271018）

0 引 言

洪水灾害是目前发生频率最高的自然灾害，据统计，仅2018年全球洪水灾害事件共发生109起，死亡人口和受灾人口分别为1 995 人和1 262 万人，直接经济损失45亿美元［1］，其中，中国、印度、印度尼西亚、美国等国家是洪水灾害发生次数和遭受损失最多的国家。虽然在近30年间，全球洪水灾害死亡人口及受灾规模都呈现持续下降趋势，但经济损失却呈上升趋势。水库是在河流上筑坝进而起到防洪效果的一种重要水利工程，且随着城市供水规模的不断扩大和农业用水灌溉方式的多样化，水库在防洪、供水、发电及生态环境保护等方面发挥了巨大的社会和经济效益。科学的水库防洪调度设计是安全经济建设水库与实时调度发挥设计效益的关键。目前，常见的两类水库调度方法包括水库常规调度和水库优化调度，后者由于可以通过设定优化目标来制定优化作业方案，因此在水库防洪调度中广泛应用。起初，水库优化调度所采用的方法主要有线性规划［2］、非线性规划、动态规划［3］、逐次优化算法［4］及大系统分解协调原理［5］等，然而，面对由水库群、蓄滞洪区、湖泊及其他防洪工程组成的复杂防洪系统，上述优化算法表现出收敛效率低、维度低等明显局限性。

现代计算技术的进步使得应用生物学、物理学和人工智能原理的启发式智能算法得到了广泛发展，这些算法可以有效解决上述问题，主要包括遗传算法［6］、蚁群算法［7］、人工鱼群算法［8］及粒子群算法［9］等。由于智能启发式算法相比传统的优化算法能够解决非线性、间断、不可微和多维的复杂问题，因此在水库优化调度中得到了广泛的应用。然而，启发式智能算法也存在不足，例如遗传算法易早熟且计算量大；蚁群算法、人工鱼群算法过于复杂且计算量大；粒子群算法易陷入局部最优值。上述智能算法因其各自寻优时存在的弊端，使得应用受到限制。果蝇优化算法（FOA）［10］是近年来提出的一种新型智能优化算法，与其他算法相比，它具有原理简单、参数少、计算量小及易于实现等优点，故受到了国内外众多学科领域学者的广泛关注。通过对该算法的大量研究发现，该算法在遇到高维优化问题时，存在过早收敛和易于陷入局部最优等问题，影响了算法的稳定性［11］。为弥补FOA的不足，采用设置变动的搜索距离L和模拟退火算法（Simulated Annealing algorithm，SA）的Metropolis准则［12］对FOA的迭代过程进行扰动，有效避免上述缺陷。

基于优先考虑水库及下游防护区防洪安全的前提，建立水库双目标防洪优化调度模型；然后，提出耦合模拟退火算法的果蝇优化算法（SA-FOA）对模型进行求解；最后，以日照水库2020年的两场洪水调度过程为例，验证了该算法的可行性和合理性。

1 水库双目标防洪优化调度模型的建立和求解

1.1 水库双目标防洪优化调度模型的建立

水库自身、库区和下游防护区是水库防洪调度的3个主要防护对象，其中，为保证水库及库区安全，希望水库汛期尽可能贴近汛限水位运行，故常以水库调洪最高水位最低或最大防洪库容最小作为目标函数；为保障下游防护区安全，希望水库能够最大限度地削减洪峰流量，故常以最大削峰准则作为目标函数［13］。本文综合考虑3个主要防护对象，建立水库双目标防洪优化调度模型，以实现调洪最高水位最低和最大限度削减洪峰流量。

1.1.1 目标函数

（1）水库及库区防洪安全目标。调洪最高水位越小，水库及库区越安全，因此以坝前最高水位作为目标函数，即：

式中：T为调洪历时；Z(t)为t时刻水库坝上水位，m。

（2）下游防护区安全目标。经证明，水库最大削峰准则的目标是水库出库流量过程平方和最小，即下泄流量最均匀［14］。利用该准则所构建的目标函数的表达式为［14-16］：

式中：q(t)为t时刻水库的泄流量，m3/s；R（t）为t时刻水库到下游防护区之间的区间来水流量，如无区间来水，其取值为0，m3/s。

1.1.2 约束条件

（1）水量平衡方程。

（2）水库下泄流量约束。

（3）水库泄流能力约束。

（4）库容约束。

（5）非负约束：以上变量均为非负。

式中：Q(t)为第t时段水库的入库流量，m3/s；V（t）为第t时段水库的库容，m3；Δt为时段长，h；Φ[V(t)]为水库蓄水量为V(t)时的最大下泄能力，m3/s；q(t)min、q(t)max分别为第t时段水库的泄流量下限、上限，m3/s；V(t)min、V(t)max分别为第t时段水库的库容下限、上限，m3。

1.2 模型的求解

1.2.1 FOA计算步骤

FOA由潘文超于2011年提出［17］，是一种通过模拟果蝇种群觅食行为寻求最优解的方法，其寻优过程如下：

（1）初始参数的设置，包括果蝇群里初始位置X0和Y0，最大觅食代数Gmax，搜索距离L，种群规模size。

（2）果蝇个体通过嗅觉随机搜寻食物，其位置为：

式中：i=1，2，…，size。

（3）计算果蝇个体与原点之间的距离Disti，再计算味道浓度判断值Si，该值可取为距离的倒数，即：

（4）将Si代入味道浓度判断函数中，求解果蝇个体的味道浓度Smelli：

（5）找出果蝇种群中味道浓度最佳值Smellb，如果目标函数是求解最小值，则味道浓度最小值的个体为最佳，反之亦然。

（6）保留最佳味道浓度Smellb及其位置坐标，此时果蝇种群利用视觉向该位置飞行。

（7）迭代寻优，重复步骤（2）～（4），判断当前味道浓度是否优于前一次迭代的最佳味道浓度，并判断当前迭代次数是否小于Gmax，若是执行步骤（6）；否则，结束迭代输出最优值。

1.2.2 SA算法

SA算法是Kirkpatrick等人［18］受到固体或热力学系统的物理退火法的启发而提出的一种概率优化方法。在该方法中，热力系统的当前状态、能量方程和基态类似于优化问题的当前调度解、目标函数和全局最优解。SA算法使用基于概率的Me‐tropolis准则搜索全局最优并跳出局部最优，该准则设置接受劣解的概率［18，19］，其定义如下：

式中：G为迭代次数；pG是第G次迭代中的接受概率；EG是第G次迭代中的目标函数值；Eg是历史最优目标函数值；KbO是Boltzmann系数；TG是第G次迭代中的退火温度，该温度最初设置为较大的值，然后通过以下温度控制函数减小到较小的值：

式中：α为退火系数，其取值区间为［0.8，1.0］。因此，随着TG的逐渐减小，劣解的接受概率将接近0。

1.2.3 SA-FOA

FOA中果蝇移动的搜索距离L会影响其寻优的结果，L大有利于跳出局部极值进而寻求全局最优，并且提高收敛效率，但其局部搜索能力将下降；L过小会导致陷入局部极值，且降低收敛效率。在算法迭代的过程中，对L大小的需求不同。在算法迭代初期，需要较大的L以提高收敛速度；在迭代后期，需要较小的L以避免错过最优值。因此，本文采用随着迭代次数逐渐减小的搜索距离，公式如下：

再者引入SA算法，将目标函数增量∆s=SmellbG+1-SmellbG作为衡量标准，其中，SmellbG为第G次迭代中果蝇种群中味道浓度最佳值，根据Metropolis准则依概率对产生的新个体进行选择，选取进入下一代的群体。具体方法见1.2.4节的第（6）步。该方法可以很好地确定搜索方向，又能依概率舍弃，跳出局部最优，使算法能够较快的搜索到最优位置。

1.2.4 SA-FOA确定水库防洪优化调度方案的步骤

SA-FOA确定水库防洪优化调度方案的步骤如下：

（1）初始参数的设置，包括最大觅食代数Gmax，初始搜索距离L0，当前觅食代数G=1，种群规模size，空间维度dim=T，即洪水时段数，初始温度TG，退火系数α。

（2）随机生成初始种群，即水库初始泄流过程线，即：

（3）将种群代入味道浓度判断函数中，求出果蝇个体的味道浓度值，即：

式中：为第i个种群在第G次觅食时的最高水位，m；Zmax和Zmin分别为调度期水库允许的最高水位和最低水位，m；为第i个种群在第G次觅食时泄流量的平方和分别为第G次觅食时，所有种群泄流量平方和的最大值和最小值。

（4）寻求果蝇种群中味道浓度最佳的个体。本文是求最小值问题，则果蝇种群中味道浓度最小值的个体为最佳，即SmellbG。

（5）采用式（13）计算LG，果蝇群体的位置变为：

再由第（3）和第（4）步骤求出SmellbG+1。

（6）计算味道浓度差值∆s=SmellbG+1-SmellbG。依据Metropolis准则，当Δs≤0时，SmellbG+1被接受，保存该果蝇个体；当Δs>0时，依据式（11）的概率pG=exp(-Δs/TG)收敛。若此时的pG>rand，则新的最小值被接受，将(t)替代作为下次寻优的初始位置；否则，保持不变，开始进行下一次寻优。

（7）令G=G+1，如果G≤Gmax，转入第（3）步继续进行迭代，否则停止迭代，输出最优解。

2 实例研究

2.1 日照水库简介

日照水库位于山东省日照市东港区，控制流域面积548 km2，总库容3.18亿m3，兴利库容1.82 亿m3，是一座集防洪、供水、发电、灌溉等多种功能于一体的大（II）型水库，同时也是日照市区重要饮用水水源地。日照水库保护着下游后村镇、高兴镇、南湖镇等7个乡镇、街道办事处90个自然村，9.06万人民群众的生命财产安全，及204国道、兖石铁路、日兰和沈海高速公路和青连高铁等交通干线以及临港大中型企业等防洪重点保护对象。日照水库防洪标准为100年一遇设计、5 000年一遇校核，汛限水位42.50 m、兴利水位43.00 m、防洪高水位43.50 m、设计洪水位43.97 m、校核洪水位46.52 m。汛期洪水调度原则为，当雨前或预计雨后水位低于42.50 m时，不泄洪；当雨前或预计雨后水位42.50～43.50 m时，控泄流量1 000 m3/s；当雨前或预计雨后水位高于43.50 m时，溢洪闸自由敞泄；当水位高于43.97 m时，北放水洞参与泄洪。

2.2 结果对比

收集了日照水库20200813和20200826两场洪水资料，包括水库的水位、库容、入库流量和泄流量，时段长为2 h，同时收集了水库的水位～库容～泄量关系等数据，满足建立水库防洪优化调度模型的要求。参考已有文献［11，17-19］及对两场洪水的计算，确定SA-FOA的初始参数为Gmax=2 000，L0=2，size=100，T1=106，α=0.9。采用SA-FOA对日照水库双目标防洪优化调度模型进行求解，得到非劣解集中各调度方案的目标函数见表1和表2，同时采用FOA进行求解，两个算法求得的调度方案非劣解前端见图1和图2，两场洪水部分方案的泄流过程见图3和图4。

图1 20200813洪水调度非劣方案集Fig.1 Non-inferior scheme set of 20200813 flood

图2 20200826洪水调度非劣方案集Fig.2 Non-inferior scheme set of 20200826 flood

图3 20200813洪水部分调度方案下泄流量过程Fig.3 Discharges of partial operation schemes of 20200813 flood

图4 20200826洪水部分调度方案下泄流量过程Fig.4 Discharges of partial operation schemes of 20200826 flood

表1 20200813洪水调度非劣方案集Tab.1 Non-inferior scheme set of 20200813 flood

表2 20200826洪水调度非劣方案集Tab.2 Non-inferior scheme set of 20200826 flood

（1）由表1、2和图1、2可见，由SA-FOA和FOA求得的调度方案集的最高水位和泄流量平方和近似成反比关系，这是因为两者代表的目标存在冲突导致的。由图1和图2可见，在同一最高水位时，SA-FOA求得的泄流量平方和小于FOA，在同一泄流量平方和时，SA-FOA求得的最高水位小于FOA，由此可见，SA-FOA能够比FOA求得更优的目标值。

（2）20200821洪水的实测最高水位和泄流量平方和分别为42.74 m和3 477 174，由表1可见，均大于SA-FOA求得方案集的目标值，比方案集中最高水位的最大值和最小值分别大0.17%和0.25%，比方案集中的泄流量平方和的最大值和最小值分别大14.87%和22.35%。20200826洪水的实测最高水位和泄流量平方和分别为42.90 m和1 133 554，由表2可见，同样均大于SA-FOA求得方案集的目标值，比方案集中最高水位的最大值和最小值分别大0.27%和0.35%，比方案集中的泄流量平方和的最大值和最小值分别大5.63%和71.54%。综上，SAFOA求得的泄流方案优于实测值。

（3）由图3可见，20200821洪水的4个优化方案的最大泄流量均小于实测值，且比实测泄流过程提前加大了泄流，这是导致优化方案最高水位低于实测值的原因。由图4可见，20200826洪水的4个优化方案的最大泄流量均大于实测值，然而其不超过日照水库的控泄流量1 000 m3/s，故对下游防护区也是安全的。

3 结 论

（1）鉴于FOA存在过早收敛、易于陷入局部最优的不足，本文提出了SA和FOA耦合算法SA-FOA，该算法采用SA的Metropolis准则引入FOA，使其在迭代过程中依概率接受最优解，同时采用逐渐递减的搜索距离L，并将该算法求解水库防洪优化调度模型。

（2）由日照水库两场洪水调度的实例可以看出，SA-FOA求解的调度方案集的最高水位和泄流量平方和均小于实测值和FOA求得的调度方案集，由此表明SA-FOA求得的优化调度方案集既能够保证水库和库区安全，又能提高下游防护区安全。综上所述，本文建立的水库防洪优化调度模型及SA-FOA能够制定出满意的洪水调度决策方案，为决策者提供技术支持。