中国在应对气候变化及能源结构转型的过程中，尽可能地利用水能这一可再生能源，对于减轻环境压力和能源可持续发展意义重大。梯级水电站年总发电量主要取决于年径流量与其分布，以及水库调度管理。为提高水能利用率和供电水平，应充分发挥有调节性能水库的容量效益，实行梯级水电站经济调度^[1–3]。由于梯级水库间复杂而密切的水力、电力联系，梯级水电站经济调度呈多维数、动态性、非线性等特点，涉及水库入流、库容、水头、出力等多个变量，求解复杂。

当前，为解决水电站经济调度问题，学者们引进了较多的求解算法。传统的经典算法有非线性规划（NLP）^[4]、拉格朗日松弛法（LR）^[5]、动态规划（DP）^[6]、逐步优化算法（POA）^[7]等，现代启发式随机搜索方法有遗传算法（GA）^[8]、差分进化算法（DE）^[9]、人工蜂群算法（ABC）^[10]、粒子群算法（PSO）^[11]、蝙蝠算法（BA）^[12]等。其中，蝙蝠算法是2010年提出的新兴启发式随机搜索群智能算法，具有并行性、分布式、收敛速度快等特点，目前在科学领域已得到广泛应用^[13]，但该算法也存在一些不足。

针对BA中缺乏变异机制而导致的种群多样性较低，全局搜索能力较差，易陷入局部最优、收敛精度不高等缺点，将差分算法变异、选择操作引入蝙蝠算法中，改进种群更新策略，提出一种改进的蝙蝠算法（IBA）。将其运用于梯级水电站经济调度问题中。

1 数学模型

通常梯级水电站群中长期优化调度是根据给定的入流过程和综合利用要求，以发电流量为决策变量，制定并实施水电站中长期最优化运行调度方式，以获得调度周期内总发电量最大或总发电效益最大。随着中国社会经济发展及能源结构转型，电网安全稳定运行需要大容量机组参与调峰运行。尤其是对以水电为主的四川电网而言，在枯水期需提供给电网尽可能大而稳定的出力值，以保证电力系统供电的稳定性。因此，为充分发挥有调节性能水库容量效益，保障电网安全稳定运行，同时替代火力发电，减少碳排放，本文建立了兼顾梯级最小出力最大化的梯级总发电量最大模型。

1.1 目标函数

1）梯级总发电量最大

$E = { max}\sum\limits_{i = 1}^n {\sum\limits_{t = 1}^T {{K_{i,t}}} } \cdot {Q_{i,t}} \cdot {H_{i,t}} \cdot {M_t}$

(1)

2）梯级电站年内最小出力最大化

$NP = { max}\;\mathop { min}\limits_{1 \le t \le T} \sum\limits_{i = 1}^n {{K_{i,t}} \cdot {Q_{i,t}} \cdot {H_{i,t}}} $

(2)

式中： $E$ 为调度周期内梯级总发电量； $n$ 为梯级电站数量； $i$ 为电站编号， $i = 1,2, \cdots , n$ ； $T$ 为调度周期时段数（中长期优化调度周期通常为一年，计算时段为月或者旬，本文以月为计算时段，即 $T =12 $ ）； $t$ 为计算时段编号， $t = 1,2, \cdots ,T$ ； ${K_{i,t}}$ 、 ${Q_{i,t}}$ 、 ${H_{i,t}}$ 分别为第 $i$ 个电站 $t$ 时段的综合发电系数、发电流量、净水头； ${M_t}$ 为 $t$ 时段的时长； $NP$ 为整个调度期内梯级最小出力的最大值。

上述梯级水库调度模型是一个多目标复杂优化决策类问题，模型求解较为复杂，通常采用线性加权法、约束法将其转化为单目标优化问题^[1]。本文采用约束法将整个梯级调度周期内最小出力最大化的目标转化为一个约束条件，转化后有：

$\mathop { min}\limits_{1\le t \le T} \sum\limits_{i = 1}^n {{K_{i,t}} \cdot {Q_{i,t}} \cdot {H_{i,t}} \ge \varepsilon } $

(3)

式中， $\varepsilon $ 为假定的最大化后的最小出力，根据整个梯级调度周期内最小出力计算结果动态调整^[14]。

在目标函数1）中加入梯级电站最小出力约束作为惩罚项，最终的目标函数为：

$E \!=\! { max} \!\!\sum\limits_{t ={{1}}}^T {\left[\!\! {\sum\limits_{i = 1}^n {\!\!{K_{i,t}} \cdot {Q_{i,t}}} \cdot {H_{i,t}} \!-\! {\sigma _t} \cdot \lambda \cdot \!\left| {\varepsilon \!-\!\! \sum\limits_{i = 1}^n {{K_{i,t}} \cdot {Q_{i,t}} \cdot {H_{i,t}}} } \right|} \!\right]} \cdot \!{M_t}$

(4)

式中， $\lambda $ 为惩罚因子，为常数。 ${\sigma _t}$ 定义为：

${\sigma _t} = \left\{ {\begin{aligned} & {0,\;\sum\limits_{i = 1}^n {{K_{i,t}} \cdot {Q_{i,t}} \cdot {H_{i,t}} \ge \varepsilon } }; \\ & {1,\;\sum\limits_{i = 1}^n {{K_{i,t}} \cdot {Q_{i,t}} \cdot {H_{i,t}} < \varepsilon } \;} \end{aligned}} \right.$

(5)

1.2 约束条件

1）水库水量平衡

${V_{i,t + 1}} = {V_{i,t}} + ({I_{i,t}} - {Q_{i,t}} - {S_{i,t}}) \cdot {M_t}$

(6)

2）梯级水库水力联系

${I_{i + 1,t}} = {Q_{i,t}} + {S_{i,t}} + {R_{i,i + 1,t}}$

(7)

3）水库水位约束

${\textit{Z}}_{i,t}^{\min } \le {{\textit{Z}}_{i,t}} \le {\textit{Z}}_{i,t}^{\max }$

(8)

4）过机流量约束

$Q_{i,t}^{u,\min } \le {Q_{i,t}} \le Q_{i,t}^{u,{ max}}$

(9)

5）下泄流量约束

$Q_{i,t}^{\min } \le {Q_{i,t}} + {S_{i,t}} \le Q_{i,t}^{\max }$

(10)

6）单站出力约束

$N_{i,t}^{\min } \le {K_{i,t}} \cdot {Q_{i,t}} \cdot {H_{i,t}} \le N_{i,t}^{\max }$

(11)

7）调度期末水位限制

${{\textit{Z}}_{i,{{start}}}} = {{\textit{Z}}_{i,{{end}}}}$

(12)

8）所有变量非负约束（ $ \ge 0$ ）。

其中， ${V_{i,t}}$ 、 ${V_{i,t + 1}}$ 分别为 $t$ 时段初、末第 $i$ 个水库库容， ${I_{i,t}}$ 、 ${S_{i,t}}$ 分别为第 $i$ 个水库 $t$ 时段的入库流量、弃水流量， ${R_{i,i + 1,t}}$ 为第 $i$ 个水库与 $i+1$ 个水库在 $t$ 时段的区间流量， ${{\textit{Z}}_{i,t}}$ 为第 $i$ 个水库 $t$ 时段的水位， ${\textit{Z}}_{i,t}^{\min }$ 、 ${\textit{Z}}_{i,t}^{\max }$ 分别为第 $i$ 个水库 $t$ 时段允许的最低、最高水位， $Q_{i,t}^{u,\min }$ 、 $Q_{i,t}^{u,{ max}}$ 分别为第 $i$ 个电站 $t$ 时段水轮机组允许通过的最小、最大流量， $Q_{i,t}^{\min }$ 、 $Q_{i,t}^{\max }$ 分别为第 $i$ 个水库 $t$ 时段允许的最小、最大下泄流量， $N_{i,t}^{\min }$ 、 $N_{i,t}^{\max }$ 分别为第 $i$ 个电站 $t$ 时段允许的最低、最高出力， ${{\textit{Z}}_{i,{{start}}}}$ 、 ${{\textit{Z}}_{i,{{end}}}}$ 分别为水库 $i$ 调度周期始末水位。

2 改进蝙蝠算法 2.1 蝙蝠算法

蝙蝠算法^[15–17]是Yang根据微型蝙蝠回声定位能力于2010年提出的新型元启发式算法。微型蝙蝠通过超声波进行回声定位，探测猎物，且随着自身与目标距离的减小，逐渐增大脉冲发射频率，减小脉冲音量。Yang将其与优化目标功能相联系，提出一种启发式随机搜索的蝙蝠算法。

对于目标函数为 $\max f\left( { X} \right)$ ，目标变量为 ${ X}= ({x_1},$ $ {x_2}, \cdots ,{x_d}{)^{{T}}}$ 的优化问题，BA算法的实施过程如下：

Step1：参数初始化。模型参数具体包括种群个体数 $TN$ 、脉冲音量最大值 ${A_0}$ 、脉冲发射率最大值 ${R_0}$ 、脉冲频率范围 $\left[ {{f_{\min }},{f_{\max }}} \right]$ 、音量衰减系数 $\alpha $ 、脉冲发射率增强系数 $\gamma $ 、搜索精度 $\varepsilon $ 或最大迭代次数 ${g_{\max }}$ 。

Step2：蝙蝠位置初始化。蝙蝠在 $d$ 维空间中随机扩散分布一组可行解，即随机初始化每只蝙蝠的位置 ${{ X}_j} = {\left( {{x_j}_1,{x_j}_2, \cdots ,{x_{jd}}} \right)^{{T}}}$ ， $j = 1 , 2 , \cdots , TN$ 。计算个体适应度值 $f({{ X}_j})$ ，根据其优劣确定当前种群最优蝙蝠个体位置 ${X_ * }$ 。

Step3：蝙蝠的搜索脉冲频率、速度、位置更新。种群在每一次迭代过程中，更新搜索脉冲频率、速度、位置，称之为全局搜索，公式如下：

${f_j} = {f_{\min }} + ({f_{\max }} - {f_{\min }}) \cdot \beta $

(13)

${ V}_j^{g + 1} = { V}_j^g + ({ X}_j^g - {{ X}_ * }) \cdot {f_j}$

(14)

${ X}_j^{g + 1} = { X}_j^g + { V}_j^{g + 1}$

(15)

式中： $\beta \in [0,1]$ 为均匀分布的随机向量； ${f_j}$ 为第 $j$ 只蝙蝠的搜索脉冲频率， ${f_j} \in \left[ {{f_{\min }},{f_{\max }}} \right]$ ； ${ V}_j^g$ 、 ${ V}_j^{g + 1}$ 分别为第 $j$ 只蝙蝠在第 $g$ 、 $g+1$ 次迭代时的速度； ${ X}_j^g$ 、 ${ X}_j^{g + 1}$ 分别为第 $j$ 只蝙蝠在第 $g$ 、 $g+1$ 次迭代时的位置； ${{ X}_ * }$ 为当前种群中最优蝙蝠个体的位置。

Step4：生成均匀分布随机数 $rand$ ，如果 $rand > {r_j}$ （第 $j$ 只蝙蝠脉冲发射率），则对当前种群最优解进行随机扰动，产生一个新的个体，对其进行越界处理，替代Step3产生的第 $j$ 只蝙蝠的位置。这一步称为局部搜索，随机游走产生新个体的公式为：

${{ X}_{{{new}}}} = {{ X}_{{{old}}}} + [\theta \cdot {{A}^g}]$

(16)

式中， $\theta \in [ - 1,1]$ 为均匀随机数， ${{ X}_{{{old}}}}$ 为当前种群最优解， ${{ X}_{{{new}}}}$ 为局部搜索产生新个体， ${A^g}$ 为第 $g$ 次迭代种群所有蝙蝠脉冲音量的平均值。

Step5：计算种群所有个体适应度函数值 $f({{ X}_j})$ ，生成均匀分布随机数 $rand$ ；若 $rand < {A_j}$ （第 $j$ 只蝙蝠脉冲音量）且 $f({{ X}_j}) > f({{ X}_ * })$ ，接受本次迭代更新产生的蝙蝠个体位置，按照式（17）～（18）对 ${r_j}$ 、 ${A_j}$ 进行更新：

$r_j^{g + 1} = {R_0} \cdot \left[ {1 - \exp ( - \gamma \cdot g)} \right]$

(17)

$A_j^{g + 1} = \alpha A_j^g$

(18)

式中， $\alpha $ 、 $\gamma $ 均为常数。对于任意 $\alpha \in (0,1)$ ， $\gamma > 0$ ，应有 $A_j^g \to 0$ ， $r_j^g \to {R_0}$ ，当 $g \to \infty $ 。

Step6：计算所有蝙蝠的适应度值，确定当前种群最优解和最优值。

Step7：重复Step3～6，直至满足设定的终止条件。

Step8：输出全局最优解和最优值。

2.2 算法的改进

从2.1节可知，在蝙蝠算法中，每只蝙蝠个体脉冲频率随迭代而更新，操作较复杂、费时。并且，蝙蝠算法无条件接受局部搜索产生的新解代替全局搜索产生的解，有条件接受本次迭代更新产生的新解。也即：当前迭代更新产生的所有个体适应度函数值要与当前种群最优值比较，只有优于当前种群最优的蝙蝠个体的位置更新才被接受。然而，对于有多个峰值的目标函数而言，这一更新策略很容易使函数值陷入局部最优而难以跳出。在蝙蝠算法中，只有蝙蝠个体的新解被接受，才更新其脉冲发射率 $r_j$ 和脉冲音量 $A_j$ ，这一更新策略使得蝙蝠个体脉冲发射率 $r_j$ 更新速度缓慢，降低了局部搜索概率。此外，蝙蝠算法中脉冲频率通常为正值，通过式（14）～（15）产生的新个体将与当前种群最优个体产生更大偏离，容易出现不符合约束的无效个体，且这一无效个体很可能将一直存在于种群中，难以被取代，导致重复计算。因此，针对以上缺陷，对蝙蝠算法做了以下改进：

1）在作者提出的改进蝙蝠算法中，每只蝙蝠个体的脉冲频率不随个体 迭代而更新，减少计算时间。蝙蝠个体脉冲频率初始化公式为：

${f_{j,i}} = {f_{\min }} + ({f_{\max }} - {f_{\min }}) \cdot {\beta _{j,i}}$

(19)

式中， ${f_{j,i}}$ 为第 $j$ 只蝙蝠第 $i$ 维度的搜索脉冲频率分量， ${\beta _{j,i}} \in [0,1]$ 为均匀分布随机数。

2）更改新解接受条件，即无条件接受全局搜索产生的新解 ${{ X}_j}(g + 1)$ ；对于局部搜索产生新解 ${{ X}_{1,j}}(g + 1)$ ，只有当 $f({{ X}_{1,j}}(g + 1)) > f({{ X}_j}(g + 1))$ 时，才由 ${{ X}_{1,j}}(g + 1)$ 取代 ${{ X}_j}(g + 1)$ 作为第 $j$ 只蝙蝠的新位置。

3）蝙蝠个体脉冲发射率 $r_j$ 和脉冲音量 $A_j$ 随种群迭代而更新。

4）更改速度更新公式，如式（20）所示，缩小新个体与当前种群最优个体的偏离值，使新个体靠近当前种群最优个体。经实验证明，这一举措可以明显减少计算时间，提高计算精度。

${ V}_j^g = { V}_j^{g - 1} - ({ X}_j^{g - 1} - {{ X}_ * }) \cdot {{ f}_j}$

(20)

5）针对蝙蝠算法中缺乏变异机制，种群多样性较低，算法易陷入局部最优的不足，引入DE算法中的变异、选择操作。

2.3 求解步骤

经以上改进，本文的改进蝙蝠算法为：

Step1：参数初始化。模型参数具体包括种群个体数 $TN$ 、脉冲音量最大值 ${A_0}$ 、脉冲发射率最大值 ${R_0}$ 、脉冲频率范围 $\left[ {{f_{\min }},{f_{\max }}} \right]$ 及按照式（19）初始化的蝙蝠个体脉冲频率矢量、音量衰减系数 $\alpha $ 、脉冲发射率增强系数 $\gamma $ 、缩放因子 $F$ 、搜索精度 $\varepsilon $ 或最大迭代次数 ${g_{\max }}$ 。

Step2：蝙蝠位置初始化。蝙蝠在 $d$ 维空间中随机扩散分布一组可行解，即随机初始化每只蝙蝠的位置 ${X_j} = {({x_j}_1,{x_j}_2, \cdots ,{x_{jd}})^{{T}}}$ ， $j = 1 , 2 , \cdots , TN$ 。

Step3：确定初始化种群最优个体。计算所有个体适应度函数值 $f({{ X}_j})$ ，根据其优劣确定当前种群最优蝙蝠个体位置 ${{ X}_ * }$ 。

Step4：种群全局搜索、更新。种群在每一次迭代过程中，分别按照式（20）、（15）、（17）、（18）更新搜索速度、位置、脉冲发射率、脉冲音量，对蝙蝠位置进行越界处理。

Step5：种群局部搜索。生成均匀分布随机数 $rand1$ 、 $rand2$ ，如果 $rand1 > {r_j}$ ，则按照式（16）对当前种群最优解进行随机扰动，产生一个新的个体 ${{ X}_{1,j}}$ ，对其进行越界处理，计算其适应度函数值。如果 $f({{ X}_{1,j}}) >$ $ f({{ X}_j})$ 且 $rand2 < {A_j}$ ，局部搜索新解优于全局搜索新解，则由局部搜索新解 ${{ X}_{1,j}}$ 取代全局搜索新解 ${{ X}_j}$ 。

Step6：种群差分变异、选择。生成均匀分布随机数 $rand3$ ，如果 $rand3 \le CR$ （变异概率，式（21）），则按照式（22）生成一个新的解 ${{ X}_{2,j}}$ ，进行越界处理，计算适应度函数值。如果 $f({{ X}_{2,j}}) > f({{ X}_j})$ ，则由 ${{ X}_{2,j}}$ 取代 ${{ X}_j}$ 成为第 $j$ 只蝙蝠第 $g$ 次更新的位置。

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!CR = 0.6 + g/(2 \cdot {g_{\max }})$

(21)

${{ X}_{2,j}}(g) = {{ X}_ * } + F \cdot \left( {{{ X}_b}(g) - {{ X}_c}(g)} \right)$

(22)

式中， ${{ X}_b}$ 、 ${{ X}_c}$ 为从种群中随机抽取的两个不同个体。

Step7：计算所有蝙蝠的适应度值，确定当前种群最优解和最优值。

Step8：重复Step4～7，直至满足设定的终止条件。

Step9：输出全局最优解和最优值。

在本文的改进蝙蝠算法中，按式（21）动态更新种群变异概率，使得随着算法迭代次数增加，变异概率增加；由于迭代后期个体间差异较小，也增强了算法局部搜索能力。按式（22）产生的新个体继承了当前种群最优个体的特征，同时，吸收了种群其他个体信息，增加了种群多样性，加快收敛速度。

3 实　例

将改进蝙蝠算法用于求解瀑布沟、深溪沟、枕头坝一级3个梯级水电站经济调度问题。其中，瀑布沟水电站具不完全年调节性能，深溪沟、枕头坝一级水电站均为径流式电站。3个梯级水电站基本参数如表1所示，其中，瀑布沟水电站月水位变幅不超过30 m。

3.1 算法参数测试与分析

改进蝙蝠算法（IBA）的主要控制参数有种群个体数 $TN$ 、缩放因子 $F$ 和最大迭代次数 ${g_{\max }}$ 。根据文献[15]，对于常规的优化问题，15～50个个体即可满足需求。对于控制参数 $F$ 和 ${g_{\max }}$ ，通过选取典型年进行不同参数测试。可获得该梯级水电站1945—2015年共70个水文年历史径流资料，选取 $P = 50\% $ 水文年作为典型年，径流资料见表2。瀑布沟水电站调度期初水位设为790 m。

实验测试平台：操作系统Windows 10中文版，CPU为Intel（R）Core（TM）i5–5200U，主频为2.20 GHz，内存为4.00 GB；基于VS2010平台，采用C#语言编程。

选取种群个体数 $TN = 30$ ，脉冲音量 ${A_0} = 1.0$ ，脉冲发射率最大值 ${R_0} = 1.0$ ，脉冲频率范围[0，1]，音量衰减系数 $\alpha = 0.99 $ ，脉冲发射率增强系数 $\gamma = 0.99 $ 。设置搜索精度 $\varepsilon = 0.001 $ 或最大迭代次数 ${g_{\max }}$ 两种终止条件，满足任一个终止条件算法都将结束。参数测试结果如图1所示。

由图1可看出：在其他参数相同的情况下，最大迭代次数 ${g_{\max }}$ 越大，算法收敛到全局最优点的可能性越大，解的精度越高，但算法运行时间也越长；如果 ${g_{\max }}$ 过小，由于算法搜索时间过短，易收敛至局部最优点。变异因子 $F$ 主要影响当前种群最优个体变异程度， $F$ 越小，如 $F = 0.25$ 时，新个体吸收种群其他个体信息越少，新个体与当前种群最优个体差异不大，难以跳出当前局部最优点；随着 $F$ 增大，新个体吸收种群其他个体信息逐渐增大，新个体与当前最优个体差异也逐渐增大，跳出当前局部最优的可能性也逐渐增大；但 $F$ 过大，如 $F=0.75$ 时，产生的新个体不再服从水力发电最优运行水位内在规律，尽管发电量也较大，却难以保证梯级最小出力和枯期出力的均匀性，如表3所示，因而 $F$ 取值不能过大。总的来看， $F$ 取值范围为[0，1]时，随着 $F$ 取值增大，解的精度由低变高再变低。

3.2 模拟计算结果及分析

采用IBA对典型年径流资料进行模拟计算，其中， $F=0.5$ ， ${g_{\max }} = 1 \; 000$ ，其余参数同上，模拟优化调度结果如表4所示。

由表4可知，梯级水电站在丰水期（6—10月）蓄水，在枯水期（12月—次年4月）供水，没有发生弃水现象，水位变化及梯级电站出力过程符合调度要求，改进蝙蝠算法模拟优化调度结果合理。

为分析比较算法的性能，分别用标准BA算法和容易实现且可快速收敛的POA算法进行模拟计算。实验中，蝙蝠算法的种群个体数、脉冲音量最大值、脉冲发射率最大值、脉冲频率范围、搜索精度、最大迭代次数、算法终止条件、初始种群生成方式等均与改进蝙蝠算法相同，POA算法搜索步长为0.01 m，最大迭代次数 ${g_{\max }} = 10 \; 000$ 。利用蝙蝠算法、POA算法求解问题时，设定的梯级最小出力值与改进蝙蝠算法调节后最大化的最小出力相同。考虑水力发电丰枯特性，枯期出力尽可能大而均匀才能保障电网安全稳定运行，本文将着重分析梯级电站枯期出力情况。表5为不同算法模拟计算的结果比较。

从表5中实际梯级最小出力值来看，IBA、BA、POA算法分别为1 481、1 051、708 MW，设定梯级最小出力值为1 481 MW，BA、POA算法不满足梯级最小出力约束，未能实现梯级最小出力最大化目标。此外，IBA算法得到的梯级最小出力满足梯级电站设计保证出力要求，超过梯级设计保证出力96 MW。从梯级电站枯期出力平均值及标准差来看，IBA算法优化后的枯期出力大而均匀，明显优于另两种算法。这表明采用IBA算法进行梯级电站优化调度，更能保障电网安全稳定运行，这对于接入大量风电、光电等出力波动率大的新能源电力系统而言意义重大。从梯级年发电量和算法运行时间来看，IBA算法获得的发电量较BA算法高出10⁷ kW·h，运行时间略短于后者；尽管POA算法运行时间与年发电量都优于IBA算法，但其出力值却偏小而波动率大。因此，对于梯级水电站经济调度问题，本文的IBA算法明显优于BA算法，也优于传统的POA算法。

4 结　论

为给电网提供尽可能大而稳定的枯期出力，建立了兼顾梯级最小出力最大化的总发电量最大数学模型，研究了基于差分变异、选择因子的改进蝙蝠算法及其在梯级水电站中长期经济调度中的应用，并对算法主要控制参数进行了测试分析。通过工程实例检验，表明本文的改进蝙蝠算法能在给电网提供尽可能大而稳定的枯期出力，同时缩短计算时间，获得更高质量的解。但该算法仍然存在参数选择问题，下一阶段将引入动态控制机制，避免人工试算。