随着化石能源危机、气候变化、环境污染问题的日益突出，以风、光为代表的新能源发电技术得到迅速发展和大规模的开发利用[1]。然而，由于风力发电、光伏发电具有较强波动性和不确定性特征[2]，高比例的风光发电接入对电力系统运行灵活性提出了新挑战。

常规电力系统运行通过设置一定的备用容量来应对运行中出现的负荷波动或机组故障等不确定性问题，并未考虑新能源发电带来的大幅随机波动问题。实际运行过程中出现的大量弃风弃光问题很大程度上就是因为系统缺乏应对新能源出力快速波动的灵活调节能力。为此，国际学术界提出了电力系统灵活性概念，即在考虑经济性的前提下，系统迅速响应因预见或不可预见事件导致的需求侧和供给侧的大幅功率波动的能力[3]。通过提高电力系统灵活性水平可增加系统应对新能源出力波动的能力，达到增加新能源并网规模、减少碳排放的目的[4]。

要解决高比例新能源接入对系统灵活运行带来的挑战，有2个基本问题需要解决：1）如何准确量化系统的灵活性需求。传统的备用容量确定准则[5]是针对常规发电机和负荷不确定性，并未考虑新能源发电带来的不确定性，尤其是未考虑新能源发电的随机波动性。因此需要研究考虑新能源出力不确定性、波动性和相关性的灵活性需求确定方法。2）灵活性调度。在常规经济调度的基础上，需要考虑新能源出力的不确定性和波动性所引发的潮流随机波动，以及功率波动在电网中的传播，即对灵活性资源的调配提出了新问题。

针对高比例可再生能源电力系统的灵活性需求问题，大多数学者着眼于灵活性需求与供给的平衡机理与灵活性评价研究[6-11]，也有学者对灵活性需求从某一方面进行了定义和量化。文献[12]将灵活性需求定义为在调度时间窗口内为保证电网实时平衡所需的机组爬坡，以净负荷的爬坡时间及功率来衡量，但此方法只从净负荷中提取系统的总爬坡时间及功率，未对可再生能源发电及负荷的不确定性、波动性及相关性进行综合考量。针对高比例可再生能源电力系统的优化调度问题，有学者通过备用的优化方法提升系统灵活性[13-14]，增加新能源的消纳；有学者对风光的出力特征进行建模[15-20]，在调度中以经济性和可再生能源消纳为目标进行优化。文献[21]仅以风电出力预测区间的形式来考虑风电出力的不确定性，采用鲁棒调度模型求解兼顾经济性与风电消纳的调度方案；文献[22]利用历史统计信息将新能源出力由确定性预测序列转化为出力场景集，以此来表征新能源出力的不确定性。文献[23]提出电源灵活性裕度指标，并在调度中对其加以约束以应对风电的随机波动。总体上看，这些研究未对不确定性、波动性及相关性进行综合考量并建模，缺少对灵活性需求的系统量化，也没有考虑风光出力波动在空间上对线路潮流的影响，可能造成灵活性供给与需求的不匹配。

本文针对含高比例可再生能源电力系统的灵活运行问题：1）建立了综合考虑随机源(包括风电、光伏及负荷)出力不确定性、波动性和相关性的灵活性需求模型。以多胞体的形式对多个随机源的联合波动范围进行统一刻画，并利用多胞体投影方法得到随机源出力波动对线路潮流的波动范围，以便准确匹配灵活性供需的空间平衡；2）提出基于灵活性需求约束的灵活性资源滚动调度方法。利用更新的随机源预测值，滚动构造系统的灵活性需求(多胞体)，考虑灵活性需求的分布特征，以经济性和灵活性为目标确定调度计划。最后通过改进的IEEE-39节点系统及我国西北高比例新能源系统进行算例分析，验证所提模型的正确性和有效性。

本文将灵活性需求定义为在一定时间尺度内系统应对新能源出力的不确定和波动及负荷不确定，维持实时供需平衡所需的功率调节能力，包括调节功率、调节速率和调节区域的实时匹配。灵活性需求的产生主要来自2方面：一是随机源的不确定性，包括新能源出力和负荷不确定性，体现在随机源出力的预测值与实际值之间存在偏差；二是新能源出力的波动性，即出力在一段时间内的随机变化。对于风、光这类可再生能源发电，还需要考虑它们之间的相关性，因为不同位置的风资源或风光之间存在的互补性会对灵活性需求产生较大影响。

为了能够准确描述多个随机源的出力波动情况，本文将每个随机源出力作为一个随机时间序列，采用概率统计特性描述随机源的不确定性、波动性及相关性，利用多胞体构建一定置信度下多随机源的联合波动域。整体思路为：对于每个随机源的出力预测序列，利用其历史预测误差的概率统计信息描述其出力的不确定性；以出力序列的差分描述其波动性；利用主成分分析法考虑多个随机源之间的相关性；利用多随机序列的相关运算得到综合考虑多随机源不确定性、波动性和相关性的多胞体，并将其合理转化为系统的灵活性需求。通过对灵活性需求多胞体在不同维度上的投影，得到系统净负荷(负荷减去风光出力)的波动区间，以及对线路潮流的影响，以描述随机源出力波动在空间上的传播，以便准确量化灵活性需求的空间分布。

本文借鉴文献[24]的方法，在其考虑波动性和相关性的基础上，考虑到新能源出力较大的预测误差，通过对模型的改进，增加对随机源出力不确定性的考虑，形成同时考虑多随机源不确定性、波动性及相关性的联合波动域。模型能够根据需要设定的置信水平，进而将联合波动域转换为系统的灵活性需求。具体过程如下。

1）数据标准化。

模型基于当前时刻t0及预测窗口ΞH(长度为TH)内的预测值，以各随机源出力的历史数据(长度为Ts)及其对应预测数据为输入，最终以多胞体形式得到预测窗口ΞH内随机源的联合波动域。

由于各随机源的大小存在差异，为便于统一计算，需要对各随机数据序列进行标准化处理。标准化处理包括2个方面：一是标幺化；二是去周期化。由于风电没有统一的日周期性特征，因此只进行标幺化处理；而光伏出力具有明显的日周期性，因此通过标准化过程消除其出力趋势与时间的关联，即同时实现标幺化和去周期化，以便进行多随机源波动域的计算。

① 光伏出力。对长度为Ts的历史实际数据$\delta _n^{\rm{s}}(t)$(t=t0-Ts+1, …, t0-1, t0)及长度为Ts+TH的预测序列$\delta _{n, {\rm{pre}}}^{\rm{s}}(t)$(t=t0-Ts+1, …, t0, …, t0+TH)，标准化方法为

式中$\delta _n^{{\rm{CS}}}(t)$为节点n处光伏电站在t时刻晴空条件下的理论最大出力[25]，将其作为基准值。则$\delta _n^{{\rm{s'}}}(t)$和$\delta _{n, {\rm{pre}}}^{{\rm{s'}}}(t)$为标准化后的数据序列。

② 对于风电场或负荷的实际数据ωn(t) (t=t0- Ts+1, …, t0-1, t0)、预测出力数据ωn, pre(t) (t=t0-Ts+ 1, …, t0, …, t0+TH)，标准化方法为

式中μn和σn分别为ωn(t)的均值和标准差。则${\omega '_n}(t)$和${\omega '_{n, {\rm{pre}}}}(t)$为对数据进行标准化后的序列。

2）提取相关性。

由于系统中某些随机源之间具有一定的相关性，比如地理位置相近或一个风带上多风场的出力具有较大的一致性，风光之间可能有互补性。考虑和不考虑这种相关性对确定整体波动域大小有较大影响，进而会影响系统所需灵活性需求的大小，因此，恰当考虑这种相关性很必要。本文采用主成分分析法(principal component analysis，PCA)对多随机源序列进行分析。首先对上述标准化后的序列组成的历史数据矩阵D和预测数据矩阵Dpre进行主成分分析，将其变换为相互独立的序列矩阵Z和Zpre，Z和Zpre即为D和Dpre的主成分(principal component，PC)。

式中：J为随机源的数量；$[\boldsymbol{U}_1^\boldsymbol{D} \ \ \boldsymbol{U}_2^\boldsymbol{D} \cdots \boldsymbol{U}_J^\boldsymbol{D}]$为由D的协方差矩阵的特征向量所组成的矩阵。

上述步骤实现了利用正交变换将可能存在相关性的变量转换成一组线性不相关的变量，即主成分。

3）确定每个主成分的波动区间。

对于步骤2中形成的每一个主成分PC，即可根据其历史预测数据的误差统计信息形成预测窗口ΞH内一定置信度下的不确定性区间$\{ \bar \varepsilon _j^ \uparrow (\tau ), \bar \varepsilon _j^ \downarrow (\tau )\} $ $({\varXi _{\rm{H}}} = \{ 1, 2, \cdots , {T_{\rm{H}}}\} , \tau \in {\varXi _{\rm{H}}})$，方法如下：

首先计算每个主成分的预测误差序列${\bar Z_j}(t)$：

对${\bar Z_j}(t)$进行统计分析，取一定的置信度α(比如95%)，即可得到该随机源指定置信度下的不确定性区间：

式中：P(·)代表事件发生的概率；$\bar \varepsilon _j^ \uparrow $和$\bar \varepsilon _j^ \downarrow $即为第j 个PC不确定性区间的上下界。

同样，利用对序列的差分可以得到不同时间间隔的波动功率序列，根据其统计信息形成一定置信度下的波动区间$\{ \hat \varepsilon _j^ \uparrow (\tau ), \hat \varepsilon _j^ \downarrow (\tau )\} $[19]：

式中：${\hat Z_j}(t, \tau )$为第j个主成分时间尺度为τ的波动功率序列；β为主成分波动区间的置信度；$\hat \varepsilon _j^ \uparrow $和$\hat \varepsilon _j^ \downarrow $为第j(1≤j≤J)个PC波动区间的上下界。

为同时衡量不确定性和波动性的影响，采用如下方法计算综合考虑不确定性和波动性的出力区间$\{ e_j^ \uparrow (\tau ), e_j^ \downarrow (\tau )\} $(τ∈ΞH)：

4）确定多随机源波动域。

上述对每个随机源(主成分)构造了其不确定和波动区间，要确定J个随机源的联合波动域可通过笛卡尔积构造J维多胞体来表示。基于当前时刻t0，针对预测窗口ΞH内的每一个时间节点τ∈ΞH，利用步骤3中构造的主成分出力区间，运用笛卡尔积得到多胞体E(t0, τ)：

5）考虑相关性的波动域。

上述基于主成分构造的多随机源波动域是不含相关性的，且是以标准化后的数据序列构造的，因此，需要对构造的波动域多胞体进行处理：一是重新插入序列的相关性；二是对标准化的数据进行逆变换，还原到各自量级，具体反变换方法见附录A。

6）确定灵活性需求。

由于在每个随机源的不确定区间和波动区间构造时都是基于大量历史数据的统计分析，所得到的波动域多胞体覆盖了随机场景的绝大多数，比如：利用历史预测误差95%的分位数作为不确定区间，意味着所构造的不确定域能覆盖将来95%的不确定性；波动区间也一样。通过调整置信度α和β的大小，可得到不同灵活性水平下的灵活性需求大小。因此，本文将所构造的多随机源联合波动域直接等价为系统所需要的灵活性需求(多胞体)。可以认为上述方法是对传统备用需求确定方法(比如：备用系数对应的区间法)的拓展，即建立了综合考虑多随机源的不确定性、波动性和相关性的灵活性需求确定方法。

图 1展示了以2个随机源为例的灵活性需求多胞体示意图，图中${y_{{n_1}}}$轴和${y_{{n_2}}}$轴为节点n1和n2两个随机源的输出功率，将t0时刻后的多个二维多胞体E(t0, τ)按时间轴方向排列，由于随机源的不确定及波动区间会随时间的推移而增大，所以多个多胞体排列后会呈现近似于圆锥状的波动域。

根据上述方法可以得到系统的灵活性需求多胞体，通过对其在45°方向进行投影，可以得到净负荷的波动区间$\{ e_{\rm{A}}^ \uparrow ({t_0}, \tau ), e_{\rm{A}}^ \downarrow ({t_0}, \tau )\} $。还可利用节点对线路的功率转移分配因子(power transfer distribution factor)，可以将多胞体向线路方向进行投影，得到线路l上的潮流波动区间$\{ e_{{F_l}}^ \uparrow ({t_0}, \tau ), $$e_{{F_l}}^ \downarrow ({t_0}, \tau )\} $。具体投影原理及计算过程参考附录B。净负荷波动区间描述了多个随机源的联合波动范围，潮流波动区间描述了随机源输出功率的波动性在空间上的传播，以便准确调度灵活性需求的空间分布。

基于上文建立的灵活性需求模型，进行日内滚动灵活性调度。调度模型每隔一个调度间隔执行一次，同时在求取调度计划时设置一定长度的预测窗口，利用随机源的超短期预测数据，考虑预测窗口时间内系统的灵活性需求，以经济性和灵活性为目标，滚动确定当前时刻的最优调度计划。

目标函数是预测窗口内的总运行成本，由2部分组成：基于负荷预测和风光出力预测的发电成本；满足灵活性需求的调度成本，并设置权重因子γ调整灵活性调度成本在总目标中的大小。同时在每一部分成本中考虑到了可能的弃风、弃光、切负荷费用。

式中：Ci(gi(τ))为火电机组i的燃料成本；gi(τ)为机组i的出力；ui(t)为机组i的运行状态变量；Ξi和Ξn分别为火电机组及随机源所在节点的集合；πd、πw和πs分别为切负荷、弃风、弃光的惩罚成本，$x_n^{\rm{d}}(\tau )$、$x_n^{\rm{w}}(\tau )$、$x_n^{\rm{s}}(\tau )$为切负荷、弃风及弃光功率。式中下标为pre的变量为利用预测数据得到的调度计划，下标为q的变量为利用灵活性需求多胞体数据得到的调度计划。预测窗口ΞH={0, 1, 2, …, TH}为τ的集合，方向集合Ξq={↑, ↓}，包含上下两个方向，q代表Ξq中的某一方向。

上述目标函数的涵义是：每个随机源的出力有一条预测序列，多个随机源形成一个灵活性需求多胞体。由于多胞体的投影为区间，存在区间的上下边界。因此，本文实际上是采用一种区间优化方法：以预测值序列、多胞体投影区间的上下两条边界分别求取调度计划(发电机出力、弃风、弃光、切负荷)，并设置初值约束使每一种调度方案在t0时刻的调度结果相同。即目标函数包含了3套调度计划的成本，最终只取t0时刻的调度结果为实际的调度方案，但在制定该调度方案时考虑了未来时段的灵活性需求，在实现调度方案的经济性最优的同时，保证系统运行的灵活性要求。

1）系统功率平衡约束。

式中：t0为当前时刻，$e_{\rm{A}}^q({t_0}, \tau )$为由系统灵活性需求多胞体的投影得到的净负荷波动区间的上界或下界，dn, pre(τ)、wn, pre(τ)和sn, pre(τ)分别为负荷、风电及光伏在τ时刻的出力预测值。式(19)包含2个边界的功率平衡，式(20)对应预测值功率平衡。

2）线路传输功率约束。

考虑到随机源不确定性与波动性在空间上的传播，计及其所带来的潮流波动，采用直流潮流计算方法对线路传输功率加以约束。

式中：$f_l^{\max }$为线路l的最大传输功率；fl, pre(τ)为利用随机源预测值得到的线路潮流；fl, q(τ)为考虑灵活性需求空间特性得到的线路潮流；hli为节点i对线路l的功率转移分布因子；$e_{{F_l}}^q({t_{\rm{0}}}, \tau )$为系统灵活性需求多胞体投影到线路l上造成的潮流波动区间的上、下界。

3）弃风弃光切负荷约束。

4）初值约束。

上述由预测值及灵活性需求包络线分别得到的调度方案的初值进行约束，使各调度方案在t0时刻的调度结果相同，从而将gi, pre(0)、$x_{n, {\rm{pre}}}^{\rm{d}}(0)$、$x_{n, {\rm{pre}}}^{\rm{w}}(0)$和$x_{n, {\rm{pre}}}^{\rm{s}}(0)$作为最终的调度计划。

5）火电机组运行约束。

火电机组的出力约束为

式中$P_i^{\min }$和$P_i^{\max }$为机组i的最小/最大技术出力。

爬坡约束为

式中$r_m^{{\rm{up}}}$和$r_m^{{\rm{down}}}$为机组i的上爬坡和下爬坡功率。

采用改进的IEEE 10机39节点系统及我国西北高比例新能源电力系统对本文模型进行测试。模型采用MATLAB R2016a编程实现，采用CPLEX进行求解。

测试系统总装机容量为6327MW，负荷峰值3350MW，谷值为2878MW。为了构造高比例新能源系统，在节点16和30分别接入了容量为900MW的大型风电场1和风电场2，节点6接入了容量为900MW的大型光伏电站。系统装有10台常规火电机组，燃料成本采用二次函数分段线性化方法计算，机组参数见附录C。不考虑机组的强迫停运率。选取15天的数据进行调度模拟，即包含15种场景。其中风光电站及负荷的出力数据取自实际运行数据，预测数据选取对应的超短期预测数据，风电、光伏及负荷的平均预测误差分别为：15.26%、14.27%、6.89%。机组组合计划在日前根据预测值以经济性为目标制定，本文中考虑灵活性需求的滚动调度模型每15min执行一次，设置长度为2h的预测窗口，利用灵活性需求模型计算预测窗口内的灵活性需求多胞体，置信度α、β均取95%。调度模型中灵活性需求考虑权重γ取1。采用长度为1年的历史数据进行区间确定。图 2给出了15个典型天内的风电场及光伏电站出力曲线，图 3给出了新能源总出力，负荷及净负荷曲线。为了得到更好的显示效果，图中将15天的数据纵向连接显示。

图 4是风电场1和2的出力散点图，从中可以看出风电场1和2具有明显的相关性。

在第一天的调度时段内选取波动性较强的时间节点，观察灵活性需求量化模型对功率波动的包络效果。图 5为利用灵活性需求模型构造的多胞体投影得到的波动区间，其中风电场1出力在12:15—13:30时间内向上波动167.3MW，风电场2出力在16:45—18:15时间内向上波动310.2MW，光伏电站出力在12:15—13:30内向下波动273.2MW。在各个场景下，即使功率波动幅度较大或实际功率曲线大幅偏移预测曲线，所给出的灵活性需求多胞体均能很好地包络实际功率波动，可见由此制定的调度方案能更好地应对风光出力的不确定性及波动性。

为了验证所提调度模型在经济性及灵活性方面的优势，进行对比分析，设置4种场景：Case 1仅以经济性为目标，不设置预测窗口，不考虑系统的灵活性需求；Case 2在本文模型中去掉随机源不确定性考虑；Case 3不考虑随机源的波动性；Case 4通过去除PCA变换及反变换的步骤，不考虑随机源之间的相关性。表 1展示了5种调度策略弃风、弃光、切负荷结果及运行成本的对比。

通过与Case 1对比可以发现，本文模型通过对灵活性需求的合理量化，能够在保证经济性的同时减少弃风、弃光及切负荷现象的发生；通过与Case 2、3、4对比可以发现，本文模型通过对随机源不确定性及相关性的考虑，能够更加准确地量化系统的灵活性需求，精确刻画随机源的波动范围，使多胞体更加覆盖系统的实际运行情况，从而能在调度中得到更优的决策结果，降低系统的运行成本，得到更好的灵活性及经济性效果。

为验证相关性对灵活性需求的影响，以第一天16:00时刻为例，以两个风电场和光伏电站为3个维度，做出了预测窗口内4个时间点考虑/不考虑相关性情况下的三维多胞体，如图 6所示。对比分析可以看出，考虑随机源之间相关性，波动域多胞体明显小于不考虑相关性的多胞体大小，该波动域作为灵活性需求进一步会影响调度策略，即考虑相关性能够给出更灵活、经济的调度方案。

为分析系统灵活性调度成本对调度方案的影响，对灵活性成本权重γ进行灵敏度分析。图 7展示了不同权重下的调度结果。通过分析可以得出，随着系统灵活性需求成本在目标函数中的比重逐渐加大，系统的弃风、弃光、切负荷电量，及运行成本都在下降，说明恰当考虑随机源出力波动，能给出经济性和灵活性更佳的调度方案。

利用课题组开发的西北电网算例系统进行验证。该算例系统包含65个大型风电场和35个大型光伏电站、10个等效负荷、20台常规机组及两个直流外送线路。常规机组总装机容量19900MW，风电场总装机容量20500MW，光伏电站总装机容量12500MW，图 8为系统的网络结构图。

以15天的实际数据为基础，可再生能源出力、负荷曲线及直流外送功率如图 9、10中所示。由于本文不研究直流的控制策略，因此直流外送功率事先给定。15天内负荷总量为2.04×106MW·h，可再生能源发电量2.60×106MW·h，直流外送电量3.01×106MW·h。同时设置对照组，对照组设置方法与3.1节中的Case1相同，表 2展示了两组的调度结果。通过结果可以看出，本文所提出的调度策略在保证系统的灵活可靠运行下，能够有效减少弃风、弃光及切负荷，提高风光消纳水平和运行经济性。

常规主要针对负荷不确定性和机组故障的备用容量确定方法，难以保证大规模波动性电源接入下电力系统的灵活安全运行，为此，本文建立了考虑风光出力不确定性、波动性和相关性的灵活性需求量化模型；进而考虑随机源出力波动对线路潮流的影响(即不确定性及波动性的空间分布)，构建了日内灵活性调度模型，以经济性和灵活性为目标，采用区间优化方法，滚动确定灵活性资源的调度计划。通过算例验证了模型的正确性和有效性，结果表明所提模型能够准确刻画风电、光伏及负荷的不确定性、波动性及相关性对系统灵活运行的影响，能够给出更加经济灵活的调度方案。研究成果为高比例可再生能源电力系统灵活运行提供了方法支撑。

附录见本刊网络版(http://www.dwjs.com.cn/CN/1000-3673/current.shtml)。

对于步骤4中得到的多胞体E(t0, τ)，其本质上为一个正多胞体，多胞体每一个面为平面。首先对多胞体的外包络进行插值，并以J维坐标的形式将多胞体的外包络表示为矩阵Ema：

式中：坐标的右上角标为所在维度，右下角标为坐标个数的编号，N的大小取决于插值的数量，插值越精细，N越大，多胞体外包络也越精密。

随后利用如下公式对多胞体进行主成分分析的反变换：

式中：$[\boldsymbol{U}_1^\boldsymbol{D} \ \ \boldsymbol{U}_2^\boldsymbol{D} \cdots \boldsymbol{U}_J^\boldsymbol{D}]$的含义与正文步骤2中相同，为由D的协方差矩阵的特征向量所组成的矩阵；${\boldsymbol{E}'_{{\rm{ma}}}}$为PCA反变化后多胞体外包络的矩阵表示。

图B1展示了以两个随机源(风电、光伏或负荷)为例的灵活性需求多胞体及其投影，图中${y_{{n_1}}}$轴和${y_{{n_2}}}$轴为节点n1和n2两个随机源的输出功率，以t0时刻为基础，t0+τ时刻的灵活性需求多胞体E(t0, τ)在45°方向上的投影可以得到${y_{{n_1}}}$和${y_{{n_2}}}$在t0+τ时刻的联合波动区间$\{ e_{\rm{A}}^ \uparrow ({t_0}, \tau ), e_{\rm{A}}^ \downarrow ({t_0}, \tau )\} $，其计算方式如下：

式中：Y为随机源输出功率组成的列向量；1为与Y长度相同的由1组成的行向量。

当随机源为系统内所有风电、光伏及负荷时，其联合输出功率即为系统净负荷，所得到的联合输出功率区间即为系统净负荷的波动区间，即有：

式中dn(t0)、wn(t0)和sn(t0)分别为负荷、风电及光伏在t0时刻的功率值。

同样，利用E(t0, τ)在${h_{{l_{{n_1}}}}}{y_{{n_1}}} + {h_{{l_{{n_2}}}}}{y_{{n_2}}}$方向上的投影能够量化随机源在线路l上可能造成的潮流影响。其中${h_{{l_{{n_1}}}}}$和${h_{{l_{{n_2}}}}}$为节点n1和n2对线路l的功率转移分配因子(power transfer distribution factor)。根据随机源在τ时刻输出功率范围，通过投影的方式在线路l上得到潮流波动区间$\{ e_{{F_l}}^ \uparrow ({t_0}, \tau ), e_{{F_l}}^ \downarrow ({t_0}, \tau )\} $，其计算方式如下：

式中Hl为随机源所在节点对线路l的功率转移分配因子组成的行向量。