水光互补是促进光伏消纳、提高光伏利用效益的一种重要手段，光伏的接入势必会对水库群的运行方式造成影响，因此，研究水光互补日内优化运行策略对于提高发电主体的收益水平，促进光伏消纳都具有重要意义。目前，有关专家和学者在多能互补优化运行策略方面取得了卓有成效的成果。例如：李守东^[1]将风电、光伏、抽水蓄能视为一个虚拟电源组，提出了考虑风光不确定性的协调调度策略。崔杨等^[2]基于风光火及需求响应资源提出了一种源荷日前–日内多时间尺度协调调度策略，该策略考虑了需求响应资源的响应速度和发电机组的特性。但这些成果的研究对象大多是含风光等多种能源的微电网，微电网的互补往往依靠储能设备^[3-4]、可中断负荷^[5]、交流电网^[6]，以维持功率平衡。单独针对水光互补发电系统优化运行策略的研究十分罕见。魏宏阳等^[7]根据水电调节光伏波动的不同程度、水电与光伏单独或联合并网以及基荷的不同大小提出了3种运行策略，并总结了各策略的优劣，但未曾考虑水电实际的水位优化运行策略。光伏发电本身具有与其他新能源不同的“昼发夜歇”的发电特性，在仅有水光的发电系统中，光伏如何影响水电的运行及水电如何进行蓄泄水都是值得深入探讨的。

为此，作者针对仅含水光两种能源的互补系统重点研究日内协调优化运行策略。首先，建立了兼顾光伏出力局部特性和全局特性的指标体系，并采用聚类方法对光伏实际出力场景进行缩减，选出典型日；其次，提出波动电量比指标用以衡量互补系统的发电量及出力波动双重目标，并以此为目标函数建立了水光互补日内优化运行模型；接着，结合光伏出力特性，对粒子群算法进行改进，以适应水光互补模型的求解；最后，以分布于两条河流的串并联混合开发的3座梯级水电站及其周边的光伏发电组成的复杂互补系统为例提取水库群的蓄水消落策略。

1 光伏发电场景筛选

水光互补优化是一个随机性优化问题，时序模拟法^[8-9]、典型日法^[10]、场景聚类^[11-12]等多场景分析是解决随机性问题的一种有效方法^[13]。其中，聚类作为一种无监督的方法在光伏发电功率场景缩减方面具有广泛的应用。聚类指标（变量）的选取尤为关键，通常有太阳辐照度^[14]、历史功率序列^[15-16]、基于历史数据提取的特征指标^[17]等。其中，基于特征指标的聚类可以降低变量的维度并且保留较多的原始信息，但特征指标的选取需要充分反映出光伏发电曲线的波动情况。本研究以光伏发电日内出力序列数据为基础，建立基于曲线形状的聚类指标体系，具体包含反映整体和局部情况的两类指标。

1.1 局部指标

1）最大波动幅度。其指相邻时段光伏发电功率之差的绝对值的最大值，计算式为：

${\;\;\;\;\;\;\;\;A_{\max }}=\mathop {{\rm{Max}}}\limits_{t=0}^{T - 1} \left({\rm{Abs}}\left(\frac{{{F_{t + 1}} - {F_t}}}{\mu }\right) imes 100{ ext{%}} \right)$

(1)

式中， ${A_{\max }}$ 为光伏发电功率的最大波动幅度， ${F_t}$ 、 ${F_{t + 1}}$ 分别为第t和t+1时段光伏发电总功率， $\;\mu $ 为光伏发电总装机容量，t为时段变量，T为计算周期内总时段数。

2）最大功率。其指光伏功率曲线的峰值，可表征曲线局部高低。

3）峰现时间。其指光伏曲线最大功率对应的时间。

1.2 整体指标

1）平均波动幅度^[18]。其表征的是光伏发电曲线相邻时段波动的平均情况，其值越大，表明平均波动越强。计算式为：

${\;\;\;\;\;\;\;\;\;\;A_{\rm{avg}}}=\frac{1}{{T - 1}}\sum\limits_{t=1}^{T - 1} \left(\frac{{{F_{t + 1}} - {F_t}}}{\mu } imes 100{ ext{%}} \right)$

(2)

式中， ${A_{\rm{avg}}}$ 为光伏发电功率的平均波动幅度，其他参数意义同上。

2）反向波动次数 $C$ ^[18]。其指当波动趋势反转时的计数值，即当 $ ({F}_{t+1}-{F}_{t}) imes ({F}_{t}-{F}_{t-1})<0$ 时， $C$ 的值增加1。 $C=0$ 表示在时间间隔内光伏发电功率变化趋势是连续增加或连续减小； $ C>0$ 则表示在时间间隔内光伏发电功率变化的方向是不确定的，并且该值越大，不确定性越大。

3）功率分布偏度^[19]。偏度是统计数据分布偏斜方向和程度的度量，是统计数据分布非对称程度的数字特征。光伏发电功率分布曲线的偏度可以表征曲线整体波动剧烈程度，曲线波动越剧烈，其分布偏度值越大。定义式为：

$S=\frac{1}{T}\sum\limits_{t=1}^T \left(\frac{{{{({F_t} - \overline F )} / \mu }}}{\sigma } \right)^3$

(3)

式中， $S$ 为光伏发电功率曲线分布偏度， $\overline F $ 为光伏发电的日平均功率， $\sigma $ 为光伏发电功率曲线的标准差，其他参数意义同上。

4）功率分布峰度。峰度是表征概率密度分布曲线在平均值处峰值高低的特征数。光伏发电功率分布曲线的峰度可以反映曲线峰部的尖度，峰部越尖，曲线变化速率较快，波动越大。定义式为：

$P=\frac{1}{T}\sum\limits_{t=1}^T \left(\frac{{{{({F_t} - \overline F )} /\mu }}}{\sigma } \right)^4$

(4)

式中， $P$ 为光伏发电功率曲线分布峰度，其他参数意义同上。

本研究所建立的光伏功率聚类指标体系包含上述7个指标，兼顾曲线整体与局部，对曲线形状有整体把握，将有利于光伏功率场景的缩减和筛选。

2 水光互补优化运行策略 2.1 模型构建

光伏属于新兴的间歇性能源，其出力过程具有很强的随机性，因此，对水光互补发电系统而言，发电总量和出力波动至关重要。其中，发电总量直接关系着系统的售电收入，而出力波动则关乎电网的安全稳定运行。为了兼顾系统的发电总量和出力波动，本研究提出波动电量比指标，并以此为目标函数进行水光日内互补优化模型构建。波动电量比是指一定时段内系统电量波动平均值占总发电量的比值，按照式（5）计算。

$\varepsilon=\frac{f}{e} imes 100{ ext{%}} $

(5)

式中： $\varepsilon $ 为波动电量比； $f$ 、 $e$ 分别为系统电量波动平均值和发电总量，按照式（6）计算：

$\left\{\!\!\!\! \begin{array}{l} f=\dfrac{1}{T}\displaystyle\sum\limits_t {\left| {{e_t} - \overline e } \right|}, \\ e=\displaystyle\sum\limits_t {{e_t}} \end{array} \right.$

(6)

式中， $t$ 为时段变量， $T$ 为时段总数， ${e_t}$ 和 $\overline e $ 分别为 $t$ 时段的系统发电量和计算周期内的时段平均电量，其他参数意义同上。

模型的目标包括系统发电总量最大化和出力波动最小化两项，这是一个典型的多目标问题，本研究将多目标问题转化为单目标进行求解，即以波动电量比最小化为目标函数，如式（7）所示：

$o=\min (\varepsilon )$

(7)

式中： $o$ 为目标函数的缩写，不具实际意义；其他参数意义同上。

构建水电模型如下：

$\left\{\!\!\!\! \begin{array}{l} e_{t}^{\rm{h}}=\displaystyle\sum\limits_{i=1}^n ( N_{i,t}^{} imes {M_t}),\\ {N_{i,t}}={k_i}{Q_{i,t}}({H_{i,t}} - {l_{i,t}}) \end{array} \right.$

(8)

式中，i、n分别为水电站编号和数量，h表示水电站， $e_{_t}^{\rm{h}}$ 为第t时段水电站总电量， ${N_{i,t}}$ 为第i个水电站第t时段出力， ${M_t}$ 为t时段内的小时数， ${k_i}$ 为第i个水电站的出力系数， ${Q_{i,t}}$ 、 ${H_{i,t}}$ 、 ${l_{i,t}}$ 分别为第i个水电站第t时段的发电流量、发电毛水头和水头损失，其他参数意义同上。

水电约束条件包含单站时间上的水量平衡、站间空间上的水量平衡、库容、下泄流量、发电出力约束条件，如式（9）～（13）所示：

${V_{i,t + 1}}={V_{i,t}} + ({R_{i,t}} - {Q'_{i,t}})\Delta t$

(9)

${R_{i,t}}={Q'_{i - 1,t - \Delta {T_{i - 1}}}} + {I_{i,t}}$

(10)

$ {V}_{i,t}^{\rm{min}}\le {V}_{i,t}\le {V}_{i,t}^{\rm{max}}$

(11)

$Q'^{\min }_{i,t} \le {Q'_{i,t}} \le Q'^{\max }_{i,t}$

(12)

$N_{i,t}^{\min } \le {N_{i,t}} \le N_{i,t}^{\max }$

(13)

式中， ${V_{i,t}}$ 、 ${V_{i,t + 1}}$ 分别为第i个水电站第t时段初、末水库蓄水量， ${R_{i,t}}$ 、 ${Q'_{i,t}}$ 分别为第i个水电站第t时段入库、下泄流量， $\Delta t$ 为计算时段所含秒数， ${Q'_{i - 1,t - \Delta {T_{i - 1}}}}$ 为第i–1个水电站（第i个水电站的上游电站） $t - \Delta {T_{i - 1}}$ 时刻的下泄流量， $\Delta {T_{i - 1}}$ 为第i–1个水库到第i个水库的水流滞时对应的时段数， ${I_{i,t}}$ 为第 $t$ 时刻第i–1个水电站到第i个水电站的区间平均入流， $V_{i,t}^{\min }$ 、 $V_{i,t}^{\max }$ 分别为第i个水库第t时刻的允许最小、最大蓄水量， ${Q'}^{\min }_{\!\!\!i,t}$ 、 ${Q'}^{\max }_{\!\!\!i,t}$ 分别为第i个水电站第t时段应保证的最小下泄流量和允许最大下泄流量， $N_{i,t}^{\min }$ 、 $N_{i,t}^{\max }$ 分别为第i个水电站第t时段的允许最小、最大出力，其他参数意义同上。

构建光伏发电模型如下：

$ \left\{\!\!\!\! \begin{array}{l} e_{_t}^{\rm{p}}=\displaystyle\sum\limits_{j=1}^m ( {F\!_{j,t}} imes {M_t});\\ {F\!_{j,t}}={\mu _j}{ heta _{j,t}},{\rm{ }}{ heta _{j,t}} \in [0,1) \end{array} \right.$

(14)

式中，j、m分别为光伏电站编号和数量，p表示光伏电站， $e_{_t}^{\rm{p}}$ 为第t时段光伏电站总电量， ${F\!_{j,t}}$ 为第j个光伏电站第t时段功率， ${\;\mu _j}$ 为第j个光伏电站的装机规模， ${ heta _{j,t}}$ 为第j个光伏电站第t时段的容量系数，其他参数意义同上。

2.2 模型求解

上述模型求解主要考虑逐步优化算法（progressive optimal algorithm，POA）、遗传算法（genetic algorithm，GA）和粒子群算法（particle swarm optimization，PSO）。POA算法受初始解的影响较大，易产生局部收敛解；而与GA算法相比，PSO算法没有交叉和变异的遗传操作，并且具有记忆，在大多数情况下会更快收敛，因此本研究选取PSO算法进行模型求解。

PSO算法是一种依据全体极值和个体极值寻优的基于群体智能的全局随机搜索算法。假设在一个D维的目标搜索空间中有N个粒子组成一个群落，其中，第i个粒子表示为一个D维的向量：

${\;\;\;\;\;\;{{X}}_i}=({x_{i1}},{x_{i2}}, \cdots ,{x_{iD}}),i=1,2, \cdots ,N$

(15)

第i个粒子的“飞行”速度也是一个D维向量，记为：

${\;\;\;\;\;\;{{V}}_i}=({v_{i1}},{v_{i2}}, \cdots ,{v_{iD}}),i=1,2, \cdots ,N$

(16)

第i个粒子迄今为止搜索到的最优位置称为个体极值，记为：

${\;\;\;\;\;\;{P}}_{i}^{\rm{best}}=({p_{i1}},{p_{i2}}, \cdots ,{p_{iD}}),i=1,2, \cdots ,N$

(17)

整个粒子群迄今为止搜索到的最优位置为全局极值，记为：

${{{G}}^{{\rm{best}}}}=({p_{{\rm{g}}1}},{p_{{\rm{g}}2}}, \cdots ,{p_{{\rm{g}}D}})$

(18)

在找到这两个最优值时，粒子根据式（19）更新自己的速度 ${v_{id}}$ 和位置 ${x_{id}}$ ：

${\;\;\;\;\;\;\;\;\left\{\!\!\!\! \begin{array}{l} {v_{id}}=\varphi [{v_{id}} + {c_1}{r_1}({p_{id}} - {x_{id}}) + {c_2}{r_2}({p_{{\rm g}d}} - {x_{id}})] ,\\ {x_{id}}={x_{id}} + {v_{id}} \end{array} \right.\!\!\!\!\!\!\!\!\!\!\!\!\!\!}$

(19)

式中， $ {c}_{1}$ 、 $ {c}_{2}$ 为学习因子， $\varphi $ 为收缩系数， $ {r}_{1}$ 、 ${r}_{2}$ 为[0,1]范围内的均匀随机数，v_id为第i个粒子在d维的速度，p_id为第i个粒子在d维的最优位置，p_gd为整个粒子群在d维的最优位置。

在水光互补优化模型中，光伏本身不具备可控性，在求解过程中将其作为一台不可调度的机组，只叠加电量过程，并且由于PSO算法是一种随机算法，其求解结果受随机的初始种群以及速度和位置的随机变化影响而出现不稳定现象，甚至不满足约束条件。为此，本研究结合光伏发电特性提出了一种基于Sigmoid函数的PSO算法求解该模型，并且将计算周期调整为07：00—次日07：00。

基于Sigmoid函数的PSO算法主要是在生成初始水位粒子群时，采用参照Sigmoid函数的变库容蓄水法。在有出力波动控制要求的水光互补系统中，日内理想出力过程如图1所示。由图1可知，梯级水电出力一日之内的出力变化过程应是先减小后增大。总体而言，在计算周期的后半段，也就是光伏发电的不可发电时段（21：00—次日07：00），水电出力一直保持较大值。为了满足这种变化过程，水库通常在光伏发电的可发电时段蓄水，而在不可发电时段泄水发电，这与刘娟楠等^[20]的研究成果一致。根据式（8）水电出力通常与发电流量Q成正相关关系（忽略水头影响），而一日之内的入库流量变化不大，则在光伏发电的可发电时段，水库蓄水量 $\Delta Q=R - Q$ 的变化趋势刚好与光伏出力变化趋势相同，即先增后减。

Sigmoid函数是一个在生物学中常见的S型函数，也称为S型生长曲线（图2），具有变化速率先增后减的特点，其函数表达如下：

$S\!(x)=\frac{1}{{1 + {{\rm{e}}^{ - x}}}}$

(20)

因此，本研究基于光伏出力特性，参照Sigmoid函数以变库容蓄水初始化水位。采用基于Sigmoid函数的PSO算法求解水光互补优化模型，基本流程如下：

Step1：参照Sigmoid函数对蓄水量曲线进行随机初始化，由蓄水量推求水位粒子群，包括群体规模N、每个粒子的位置 ${{{X}}_i}$ 和速度 ${{{V}}_i}$ 。

Step2：计算每个粒子的适应度 ${\zeta _{i,t}}$ 。

Step3：将每个粒子适应度值 ${\zeta _{i,t}}$ 和个体极值 ${{P}}^{\rm{best}}_{i}$ 的适应度比较，如果 ${\zeta _{i,t}}$ 大于 ${{P}}^{\rm{best}}_{i}$ 的适应度值，则用该粒子 ${{{X}}_i}$ 替换掉 ${{P}}^{\rm{best}}_{i}$ 。

Step4：对每个粒子，用其适应度值 ${\zeta _{i,t}}$ 和全局极值 ${{{G}}^{{\rm{best}}}}$ 的适应度进行比较，如果 ${\zeta _{i,t}}$ 大于 ${{{G}}^{{\rm{best}}}}$ 的适应度值，则用 ${{{X}}_i}$ 替换掉 ${{{G}}^{{\rm{best}}}}$ 。

Step5：根据式（19）更新位置 ${{{X}}_i}$ 和粒子速度 ${{{V}}_i}$ 。

Step6：如果满足精度要求或达到最大循环次数则退出；否则，返回Step2。

PSO算法求解优化问题时，算法是否收敛及其收敛速度与初始种群的生成方法密切相关。基于Sigmoid函数的PSO算法在生成初始种群时采用与光伏发电出力特性相匹配的Sigmoid生长曲线，可以大大提高算法的收敛速度，节约计算成本。

3 案例仿真 3.1 系统简介及计算输入 3.1.1 水电系统

本研究以四川省某流域梯级水电站及当地的光伏为例进行案例研究，计算时段为1 h，计算周期为1 d（07：00—次日07：00）。该流域梯级电站为串并联混合开发，各站均具有日调节能力，流域水电站群拓扑关系、主要参数及约束条件分别见图3、表1。

鉴于大部分水电机组在汛期都处于满发状态，且优先考虑防洪，无暇顾及发电问题，本研究主要侧重于平枯期的互补问题研究。

3.1.2 光伏系统

以2018年该流域所在地某装机50 MW光伏电站的日内可发电时段（07：00—21：00）功率曲线为基础数据，采用快速聚类法进行聚类分析。按照聚类结果，光伏场景分为3类，各类曲线的特点总结见表2。

表2中：Ⅰ类曲线十分光滑，是典型的晴天功率曲线；Ⅱ类曲线在峰值附近多出现锯齿状波动，是由云层暂时遮光所致，是典型的多云天气功率曲线；Ⅲ类曲线具有较低的峰值，且峰度较高在峰值附近变率较大，是典型的阴雨天功率曲线。取每一类中平均功率最大的一日为计算典型日，选取结果如图4所示。

3.2 结果与讨论 3.2.1 联合优化效果

不同运行方案下各天气典型日的各电站发电量及出力过程分别见表3和图5。

由表3可知：与水电单独优化相比，水光联合优化后，梯级各水电的发电量几乎不变，表明各方案下水电的运行结果都是较优的。受光伏发电量的影响，晴天、多云、阴雨天气各典型日水电电量占比逐渐增大，分别为85.07%、88.21%、93.25%。

由图5可知：出力过程有所差异，联合优化方案下，晴天、多云、阴雨典型日梯级水电的极小出力分别为52.19、59.84、67.23 MW，分别出现在13：00—14：00、12：00—13：00、13：00—14：00时间段；极大出力分别为97.94、104.38、89.67 MW，分别出现在次日02：00—03：00、次日01：00—02：00、次日04：00—05：00时间段。在联合优化方案下，晴天、多云和阴雨典型日下互补发电系统的总出力过程均较为平滑，表明水电为适应光伏出力波动利用水库的调蓄能力进行了自我优化和调节，同时，证明采用波动电量比控制互补发电系统的发电量和出力波动是可行的。

本研究采用皮尔逊相关系数度量各水电站与光伏发电出力的互补性，统计结果如表4所示。由表4可知：总体来看，梯级中各水电站发电出力与光伏出力的互补系数均小于0，表明了水光之间的互补性。从互补系数的绝对值看，晴天、多云、阴雨3个典型日下，均是C站、B站、A站逐渐减小，表明光伏出力与C站出力之间的互补性最好，与A站出力之间的互补性最次。

为了探讨联合优化的效果，本研究将水光联合优化（联合方案）与水光单独优化再叠加的结果（分散方案）进行对比，见表5。由表5可知：从晴天、多云和阴雨3种天气来看，联合方案的出力波动均在1 MW以内，波动电量比均小于0.05%，表明本研究提出的参照Sigmoid函数的变库容水位初始化方法是有效的。联合方案的出力极差、出力波动以及波动电量比显著小于分散方案，并且两种方案的电量相差无几，表明联合优化对光伏出力波动具有显著的调节作用。

3.2.2 光伏对水库群蓄泄水的影响

本研究采用库容利用率（实际利用库容与调节库容的比值）探讨光伏的接入对水电运行的影响。表6列出了各方案下各水库的库容利用率。

由表6可知：总体而言，光伏的接入增加了水库的库容利用率，晴天、多云、阴雨天的水电总计库容利用率分别为72.14%、69.05%、44.51%，与水电单独优化时的情况相比，分别提升了44.90、41.81、17.27%。A、B、C站的水电总计库容利用率最高分别增加了31.87、46.96、69.74%，分别出现在晴天、多云和多云典型日。

光伏发电接入系统后，水电为了适应其随机波动特性，会进行蓄放水的自我优化和调节。各方案下梯级水库群的日内末水位过程如图6所示。

由图6可知：相比水电单独优化时的结果，晴天、多云和阴雨典型日下，各水电站均出现了为光伏让道而加速蓄泄水的情况，且提高了水库最高水位。在各方案下，梯级水库群水位总体变化趋势均是泄水—蓄水—泄水。水电单独优化时，A站、B站、C站分别在07：00—16：00、07：00—09：00、07：00—14：00泄水；接着，开始蓄水，并分别于次日04：00、次日03：00、次日04：00蓄至最高水位2 448.48、2 706.34、2 572.45 m。光伏的接入加速了蓄泄水过程，以晴天典型日为例，A站、B站、C站分别在10：00、9：00、9：00降至最低水位，并开始蓄水，分别于17：00、20：00、19：00 蓄至最高水位2 449.32、2 708.10、2 573.60 m，最后再降至末水位。

因此，在调节光伏后，梯级水电在1 d（07：00—次日07：00）之内的运行策略是：为了提高互补发电系统的保证出力，控制出力波动，梯级水库群在光伏出力尚小时，先泄水发电，随着光伏出力的逐渐增大，水库群逐步蓄水，并在进入光伏不可发电时段前后开始泄水，以补足光伏的间歇性，直至在次日07：00到达日末水位。

4 结　论

本研究以1 h为计算时段，1 d（07：00—次日07：00）为计算周期，进行了水光日内互补优化研究，得出以下结果：

1）以光伏发电日内出力序列数据为基础，建立了基于曲线整体趋势和局部变化的聚类指标体系，聚类结果表明该指标体系对曲线形状有整体把握，筛选出的光伏功率场景与晴天、多云、阴雨天气典型功率曲线的特点契合度高，具有代表性。

2）提出了波动电量比指标，用于衡量互补发电系统的发电量和出力波动，并以此为目标函数建立了水光互补优化模型，结合光伏发电特性提出了一种基于Sigmoid函数的PSO算法来求解该模型。

3）案例仿真结果显示，水光互补后3种天气典型日的水光总出力波动均在1 MW以内，波动电量比均小于0.05%，验证了所提算法和模型的可行性和合理性。

4）光伏的接入提升了梯级水库群的库容利用率。在调节光伏后，梯级水电的日内运行策略是：在光伏出力尚小时先泄水发电，随着光伏出力的逐渐增大，水库群逐步蓄水，并在进入光伏不可发电时段前后开始泄水，以补足光伏的间歇性，直至到达日末水位。

除此之外，本研究中仍然存在一些不足，例如未考虑时滞的影响等。另外，本研究提出的基于Sigmoid函数的PSO算法适用于以光伏发电的可发电时段为初始时段的情况，对于其他情况的适用性需要进一步论证。