浦阳江流域梯级水库防洪调度两阶段启发式搜索算法研究

陈 佳 吴国强

(诸暨市水利局，浙江 诸暨 311800)

动态规划(DP)类算法(包括DP[1]、逐次渐进DP[2]、离散微分DP[3]、折叠DP[4]等)是解决多阶段序贯决策问题的经典方法，在梯级水库优化调度中获得了广泛应用。应用DP类算法求解梯级水库优化调度问题须满足目标函数阶段可分和状态过程无后效性两个基本条件。对于梯级水库调度而言，前一条件通常容易满足，而后一条件受河道水流传播延迟(水库的下泄流至其紧邻下游水库需要花费一定的时间)影响，严格来说不能满足。当调度决策时段为旬、月，甚至更长时，对于大多数梯级水库系统而言，河道水流传播时间远远小于调度决策时段长度，使河道水流传播延迟对调度的影响基本可以忽略不计，此时，DP类算法仍具有适用性。随着调度决策时段变短，这种影响将逐渐增强，当调度决策时段长度变至与河道水流传播时间基本相当或者更短时，这种影响将十分显著，此时，DP类算法已失去适用性。

在梯级水库防洪调度中，河道水流传播延迟对调度的影响通常不能忽略不计，梯级系统最优控制方案不仅与系统当前时刻的蓄水状态有关，而且与当前时刻以前若干时刻的蓄水状态有关(即“有后效性”问题)[5]。对于梯级水库防洪调度问题，通常采用智能算法(如遗传算法[6]、粒子群优化法[7]、水循环算法[8]等)或者逐步优化算法(POA)[9-12]进行求解。智能算法能够处理非线性、不可微、水流传播延迟等各种复杂因素，通用性强，且此类算法大多采用随机搜索机制，理论上可以产生全局解，但在实际应用中往往出现“早熟”收敛和“不稳定解”等问题[13]，应用效果受到影响。POA算法将多阶段问题分解为若干两阶段问题进行逐次求解和迭代逼近，虽然能够处理“有后效性”问题，但受限于“盲目搜索”问题和“维数障碍”[14]，对于梯级水库防洪调度问题的求解效率不高。

因此，本文以浦阳江流域防洪控制性水库联合调度为研究对象，在POA算法的基础上，设计一种两阶段启发式搜索算法，以期为解决梯级水库防洪调度“有后效性”问题和“维数障碍”问题提供一条新的途径。

1 问题描述

1.1 流域简介

浦阳江又名“浣江”，系钱塘江支流，发源于浙江省金华市浦江县，干流自南向北流经浦江县、诸暨市、萧山区，沿途接纳大陈江、开化江、五泄江、枫桥江、凰桐江、永兴河等主要支流，在萧山区境内注入钱塘江[15]。浦阳江主流长150km，流域面积3455km2，流域地处亚热带季风性气候区，受地形、天文、气象等综合因素影响，流域饱受洪、涝、台等灾害困扰，是浙江省防汛防台重点流域之一。浦阳江干、支流大部分河段已修筑堤防，但汛期发生洪水时，水面高程往往大于两岸地面高程，易发生水灾，素有“浙江小黄河”之称，夏季的梅雨和台风雨是造成流域洪涝灾害的主要因素。

1949年后，按照“上蓄、中分、下泄”的治理方针，浦阳江上游兴建了通济桥、安华、陈蔡、石壁等控制性水库，中游兴建了高湖蓄滞洪区，下游拓浚河道、裁弯取直和修建堤防，形成了以通济桥、安华、陈蔡、石壁等水库和高湖蓄滞洪区为骨干控制性工程，辅以中小型水库、干支流堤防、电排泵站等的流域防洪工程体系。流域防洪控制性水库拓扑关系和特征参数见图1和表1。

图1 浦阳江流域防洪控制性水库拓扑关系

表1 浦阳江流域防洪控制性水库特征参数 单位：亿m3

1.2 目标函数

本文以通济桥、安华、陈蔡、石壁4座水库为对象，建立浦阳江流域梯级水库联合防洪调度模型。浦阳江流域防洪调度主要目的是在保证各水库运行安全的基础上，利用水库防洪库容进行拦洪错峰，尽可能地削减下游重要防洪控制断面的洪峰流量。为了便于公式化描述，将通济桥、安华、陈蔡、石壁4座水库依次编号为1、2、3、4。目标函数表示为

(1)

(2)

式中：Rt为t时段通过诸暨太平桥控制断面的流量，m3/s；T为调度期总时段数；max(·)为取括号中集合的最大元素；min{·}为最小化括号中的项；ni为i水库与下游诸暨太平桥控制断面之间河道洪水演进水流流达时间，s；Qi,t-ni为i水库t-ni时段的下泄流量，m3/s。

1.3 约束条件

浦阳江流域梯级水库防洪调度应满足以下约束条件：

a.水库水量平衡方程：

Vi,t=Vi,t-1+(Ii,t-Qi,t)Δt

(3)

b.梯级水库水力联系：

I2,t=At+Q1,t-n

(4)

c.水库蓄水量上、下限：

(5)

d.水库泄流能力限制：

(6)

e.下游河道安全泄量限制：

(7)

f.水库下泄流量变幅限制：

(8)

g.边界条件：

(9)

h.流域洪水调度原则要求：各水库控制运用应遵循浦阳江干流洪水调度原则，详见文献[16]。

2 问题求解

2.1 启发式搜索策略

POA采用两阶段优化策略在问题的解空间内进行迭代搜索，因无须求解贝尔曼方程，且能够适应河道水流传播延迟引起的“有后效性”问题，因此在梯级水库短期调度中获得了广泛应用。然而，采用POA求解目标函数式(1)所示问题时发现，当且仅当所求解的两阶段问题涵盖洪峰流量所在时段时，才能达到削减洪峰的目的。但是，这样的两阶段问题在POA所有两阶段问题中的占比很小，仅为2/(T-1)(考虑洪峰流量为单值的情况)。换句话说，POA大部分两阶段问题的求解达不到削减洪峰的目的。因此，POA求解目标函数式(1)所示问题的效率并不高。

为了提高POA的求解效率，对POA两阶段优化策略进行改进。研究发现，目标函数式(1)的预期效果是洪水总量在不同时段之间分配尽可能均匀、洪峰流量尽可能小。因此，在峰、谷时段之间进行水量优化分配，可以达到快速改善目标函数式(1)的效果。因此，固定其他T-2个时段的流量值，仅对峰、谷两时段的流量值进行优化。不妨假定峰、谷流量对应时段号分别为p和q，由水量平衡原理可知，p、q两时段的流量之和应当为定值。因此，两阶段启发式搜索算法(以下简称“改进算法”)的优化问题可以表示为

s.t.Rp+Rq=C

(10)

式中：Rp、Rq分别为p、q时段(峰、谷流量时段)诸暨太平桥控制断面的流量，m3/s；C为常数。

本文的改进算法正是通过反复构建并求解式(10)所示两阶段问题，而达到实现洪峰流量最小化的目的的。

2.2 两阶段问题求解

对于不同形式的POA两阶段问题，因无通用的方法可供使用，通常采用穷举法进行求解。即先将连续的水库状态区间离散为由若干状态值构成的有限集合，然后对所有可能的离散状态组合方案进行比选，获得最优状态组合方案。采用穷举法求解POA两阶段问题，最大的障碍为所需计算量和存储量会随水库数量呈指数增长，即“维数障碍”。因此，采用增量法对式(10)所示问题进行求解，主要步骤如下：

a.确定各水库待优化的变量。将式(2)代入式(10)可得，安华、陈蔡、石壁3座水库待优化的变量为Qi,p-ni和Qi,q-ni(i=2,3,4)；通济桥水库待优化的变量为Q1,p-n2-n和Q1,q-n2-n。

b.固定其他变量的值，对各水库待优化的变量的值进行优化。根据水量平衡原理可知，i(i=1,2,3,4)水库2个待优化的变量之和为定值，因此，只需对1个量进行寻优，不妨假定为Qi,m,在Qi,m上、下各变动一个增量ΔQi(ΔQi取较小的正数，且对于不同的水库，ΔQi取值可以不相同)，形成3个流量值，即Qi,m+ΔQi、Qi,m和Qi,m-ΔQi。对4座水库81种流量组合方案进行比选，获得max(Rp，Rq)值最小的方案，并按此方案更新4座水库待优化的变量的值和Rp、Rq的值。

c.反复执行b步骤，直到|Rp-Rq|满足指定的精度要求为止。

2.3 主要步骤和流程图

本文的改进算法的主要步骤如下：

a.产生一条可行的水库系统初始下泄轨线，并计算诸暨太平桥控制断面的流量过程{Rt|t=1,2,…,T}。水库系统初始下泄轨线通常可以自上而下依次对各水库采用等流量调节方法产生。

b.构建式(10)所示两阶段优化问题。值得注意的是，若式(10)中的Rq已无法上调(到达约束边界上)，则应当用{Rt|t=1,2,…,T}序列中第2小流量对应时段替换q，若替换后依然无法上调，再用序列中第3小流量对应时段进行替换，依次类推，直到某个时段的流量可以上调为止。

c.采用增量法对b步骤构建的两阶段优化问题进行求解。

d.反复执行b～c步骤，直到达到最大迭代次数或者相邻两代最优下泄轨线的偏差满足指定的精度要求为止。

将上述步骤在计算机中进行编程实现。程序运行之前，需要对相关参数和变量进行初始化。参数初始化是指给搜索步长ΔQi、最大迭代次数Num和收敛精度Lim赋初值；变量初始化是指为水库系统指定一条可行的初始下泄轨线。初始化完成之后，算法将按图2所示流程运行。

图2 两阶段启发式搜索算法流程

3 案例研究

将本文的改进算法应用于浦阳江流域2021年“烟花”台风洪水优化调度中。“烟花”台风期间，浦阳江流域自7月22日开始出现降雨，至7月29日中午降雨基本结束，降雨主要集中在23—25日，流域最大3日降雨量为258.0mm，最大5日降雨量为294.6mm。该案例调度期取7月23日9时至28日7时，以1h为时段长度，总计118个时段，洪水总量为2.40亿m3。在不考虑水库调蓄的情况下，诸暨太平桥控制断面天然洪水洪峰流量为1330.50m3/s，峰现时间为26日6时，峰谷差为1302.88m3/s。为了更好地展示本文的改进算法的应用效果，将其优化结果与POA的优化结果进行对比分析。各水库调度期初蓄水量取实测值，调度期末蓄水量取各自汛限水位对应的蓄水值。POA两阶段问题采用穷举法进行求解。为简单起见，各水库搜索步长ΔQi统一取10m3/s，且与POA的离散精度相同。两种算法优化后，诸暨太平桥控制断面洪水过程对比见图3。

图3 不同算法优化后诸暨太平桥控制断面洪水过程对比

从峰值看，本文的改进算法与POA优化后，诸暨太平桥控制断面洪峰流量分别为593.47m3/s和923.47m3/s(峰现时间均为25日23时)，与天然洪水相比，削峰率分别为55.39%和30.59%，“削峰”效果显著，但前者效果更好，更大程度地保障了水库系统下游城区的防洪安全。

从洪水过程看，本文的改进算法与POA优化后，诸暨太平桥控制断面洪水过程峰谷差分别为369.58m3/s和889.37m3/s，与天然洪水过程相比，分别下降了71.63%和31.74%，“削峰填谷”效果明显，但前者效果更好，洪水过程更平稳、波动更小，方案实用性更强。

从计算耗时看，本文的改进算法耗时为3.4s，远远小于POA(POA耗时为14808.8s)，说明与后者相比，本文的改进算法效率更高(两种算法在同一台PC上进行测试，PC配置Intel Core i5-7500 3.40GHz CPU和8GB RAM)。研究发现，POA耗时长主要有两个方面的原因：ⓐ采用穷举法求解两阶段问题，计算工作量随水库数量呈指数增长，存在“维数障碍”；ⓑ大部分两阶段问题不涵盖洪峰流量所在时段，对其求解达不到削减洪峰的目的。

4 结 语

为了更好地解决梯级水库防洪调度“有后效性”问题和“维数障碍”问题，本文对POA进行改进，设计了一种POA改良版本。本文的改进算法首先采用逐次逼近策略在解空间内进行迭代搜索，解决水库系统状态过程有后效性导致难以使用DP类算法进行求解的难题；其次，对POA两阶段优化策略进行改进，采用启发式搜索策略对峰、谷流量时段进行两阶段优化，提高POA求解效率和解的质量；最后，采用增量法对两阶段问题进行求解，缓解POA的“维数障碍”。根据浦阳江流域防洪控制性水库联合调度应用结果，本文得出以下结论：

a.本文的改进算法通过迭代搜索，逼近问题的解，因无须求解贝尔曼方程，能够适应梯级系统状态过程有后效性的特点，且通过应用启发式搜索策略，与POA相比，削峰率提高了35.73%，峰谷差降低了58.44%，说明本文改进算法能够更好地解决梯级水库防洪调度“有后效性”问题。

b.本文的改进算法采用启发式搜索策略构建两阶段问题，使用增量法对两阶段问题进行求解，与POA相比，在不牺牲解的质量的情况下，算法耗时缩短了99.98%，说明本文改进算法能够更好地解决梯级水库防洪调度“维数障碍”问题。

此外，本文的改进算法概念简单，易于实现，其适用范围不局限于梯级水库防洪调度，亦可推广应用于求解其他具有min{max(·)}或者max{min(·)}形式目标函数的水资源系统优化问题(如水电调峰问题)。因本文的关注点在于改进POA的性能，模型中河道洪水演进采用简易平移法，虽然会引进一定的偏差，但不影响本文的主要结论。下一步研究将致力于完善河道洪水演进模型，进一步提高调度方案的可靠性和实用性。