来源：《中国电力》2024年第12期
引文：王蕊, 傅质馨, 王健, 等. 基于复杂特征提取和Sinkhorn距离的风光荷多阶段场景树生成方法[J]. 中国电力, 2024, 57(12): 30-40.
随着国家“双碳”目标战略提出和能源转型加速推进，高比例新能源和电动汽车等新型负荷大规模接入电网，其不确定因素增加了电力系统规划运行的复杂度。因此，须研究新能源出力和负荷在长时间尺度上的多维不确定性的降维方法，进一步为电力系统动态规划、电网承载力提升等应用提供支持。《中国电力》2024年第12期刊发了王蕊等撰写的《基于复杂特征提取和Sinkhorn距离的风光荷多阶段场景树生成方法》一文。文章提出了一种考虑风光出力和负荷增长不确定性的多阶段场景树生成方法。首先，将SSAE和基于密度峰值的改进AP算法得到的典型日风光荷场景作为场景树的根节点。其次，考虑负荷增长的不确定性，假定负荷以不同概率表现为高、中和低速增长模式，逐年生成风光荷场景树；然后，采用Sinkhorn距离算法削减场景树，由此生成一个多阶段风光荷场景树。
新能源发电出力和负荷长期增长的不确定性增加了电网规划复杂性，开展新能源出力和负荷长时间尺度上的不确定性分析，对电网的规划与建设具有重要意义。提出了一种基于复杂特征提取和Sinkhorn距离的风光荷多阶段场景树生成方法。首先，为提高风光荷场景的聚类效率，提出基于堆叠稀疏自编码器的风光荷场景特征提取方法，并采用基于密度峰值改进的近邻传播算法对风光荷场景特征集合进行聚类，获得风光荷典型曲线，作为场景树的根节点；然后，考虑负荷不同增长率，逐年生成风光荷场景树，并提出基于Sinkhorn距离的场景树削减方法以降低场景树的规模；最后，算例仿真结果表明，所提方法计算效率高，生成的风光荷多阶段场景树可反映风光出力和负荷增长的不确定性。多阶段风光荷场景树的根节点是确定的场景，本文将全年风光荷典型日作为场景树的根节点。为此，收集全年风光荷历史数据，PWT、PPV和PL分别为风电、光伏和负荷标幺值数据集合，表达式分别为
式中：num为全年总天数；PWT(i)、PPV(i)和PL(i)分别为全年第i天各小时风电、光伏和负荷的标幺值数据集合；PWT,k(i)、PPV,k(i)和PL,k(i)分别为全年内第i天第k小时风电、光伏和负荷的标幺值。PWT、PPV和PL三者构成的全时序场景集合sf为由式（4）可知，sf中有num个样本，每个样本包含一个风光荷场景。记sf中第i个样本矩阵为![]()
m为sf中样本维度。利用SSAE对sf进行特征提取得到该场景集的提取结果。建立SSAE模型，该模型由M个SAE构成。SAE包括输入层、隐含层与输出层，SSAE将前一个SAE的隐含层用作后一个SAE输入层的输入，如图1所示。SSAE输入层、隐含层和输出层神经元个数设置为m，改变输出层的神经元个数可达到降维的目的，其他参数设置参考文献[25]。
![]()
Fig.1 Stacked sparse autoencoder model
SAE的训练过程包括编码与解码。在编码过程中提取输入的风光荷场景集数据特征，即
式中：![]()
为隐含层的输出矩阵；Wa为编码权值矩阵；ba为编码偏置；σh(x)=1/(1+e−x)为激活函数。式中：![]()
为输出层输出；Wb为解码权值矩阵；bb为解码偏置。式中：JSAE为误差函数；ρ为设定的稀疏性系数；hj(sf,i)为输入是sf,i时隐含层神经元j的激活值；![]()
为隐含层第j个神经元的平均激活量；n为隐含层神经元数量。将sf输入到SSAE模型中，对SSAE模型中的SAE逐个按式（5）~（7）进行训练，直至所有SAE均完成训练，得到sf的特征集合yM，第i个特征样本为![]()
通过SAE的堆叠，SSAE能获取风光荷场景更低维度、更深层次的特征。 1.2 基于密度峰值改进的AP聚类的风光荷典型日选取方法以1.1节得到的sf的特征集合yM为聚类指标进行聚类。将生成的风光荷典型时序场景集sc作为风光荷多阶段场景树的根节点。本文采用基于密度峰值改进的AP聚类算法，具体如下。
1）相似度计算。![]()
分别为yM中第i′个和第j′个样本，计算样本之间的相似度![]()
构成相似度矩阵![]()
2）密度峰值计算。由局部密度和相对距离计算样本i′和其他样本之间的密度峰值![]()
3）更新相似度矩阵。计算得到的密度峰值![]()
为偏好值，作为相似度矩阵中对角元素，即4）计算归属度矩阵![]()
和吸引度矩阵![]()
各元素分别为式中：![]()
为第j′个样本对第i′个样本的归属度；![]()
为样本j′从样本i′获得的吸引度值。5）引入阻尼因子，控制迭代速度和防止迭代过程中发生振荡，即式中：d为迭代次数；λ∈(0,1)为阻尼因子，本文取默认值0.5。6）判断聚类中心。根据式（14）判断![]()
是否为聚类中心，如果式（14）大于0，则![]()
对应的全时序场景sf,i可作为聚类中心。7）判断是否终止迭代。当满足：1）达到设置的迭代次数；2）聚类中心不再改变，可终止迭代。若不满足，则返回步骤5）。由此得到包含聚类中心c1,c2,⋯,cK的若干个风光荷场景集合![]()
K为聚类中心的个数。聚类中心为全年风光荷典型日场景，作为场景树的根节点。此外，统计聚类得到的每个场景集合的个数，并计算频率，作为该类场景集合的概率。聚类完成后，选择轮廓系数、CH指标和DB指标3个常用的聚类指标对聚类结果分析，检验聚类的有效性。
考虑负荷增长的不确定性，设定高、中、低3种不同的负荷增长方案。根据场景树根节点K个典型日场景，假定负荷逐年增长，并且以不同概率表现为高、中、低速增长模式，逐步生成场景获得一个场景树，如图2所示。该场景树的每个点代表一个典型日72维的风光荷场景，各个场景的概率由根节点各场景概率和负荷不同增长概率决定。
![]()
由图2可知，第1年有K个可能的典型日风光荷场景，第2年有3K个可能的典型日风光荷场景，以此类推，场景树的规模是随着阶段数成比增长，须削减场景树以降低其规模。场景约简的本质是原始场景集与约简后场景集之间概率距离的优化问题，概率距离越小，两者越接近。因此本文采用基于传统的最优运输方法改进的Sinkhorn距离（Sinkhorn distance，SD）的方法约简各阶段场景集以削减风光荷场景树。另外，引入过程距离（process distance，PD）来评估削减后生成的风光荷场景树与原始风光荷场景树之间的近似度。图2中原始场景树的每个节点记为Tm,τ，每个节点代表一个72维的风光荷场景，其中，τ为原始场景树的时间步长年，m为τ下各个节点的编号。Tm,τ中任意节点对应的场景概率记为pm,τ，每一列的场景概率和为1。以nyear年为单位划分规划阶段，记为t，场景树的根节点记为“阶段1”。接下来的阶段编号为2、3等，各阶段的场景概率和也为1。每个阶段都有多个风光荷场景集待约简，寻找一个最优子集，使得子集中各个场景的新概率能够尽可能接近于原来场景集合的概率，以降低场景树的规模。在任何阶段t，需要约简的风光荷场景集个数等于前一阶段的节点数量。每个阶段待约简的风光荷场景集合用![]()
表示，t为阶段数，t=1,2,⋯,T；w=1,2,⋯,nt−1依次表示前一阶段每个节点对应的待约简场景集的索引，nt−1为前一阶段节点个数，默认情况下n1为第1节生成的典型日风光荷场景个数。设前一阶段各节点产生的叶节点数量为![]()
则下一阶段的叶节点数量为![]()
当t=2时，前一阶段节点个数为n1，则此时待约简的场景集有n1个，各个待约简场景集![]()
为由原始场景树中属于第2阶段的那一列的场景组成，其行号是前一阶段各节点对应的编号，采用基于SD方法将超集![]()
约简为![]()
由此获得各簇中场景![]()
和相应的概率；当t⩾3时，各个待约简场景集![]()
为由原始场景树中属于第t阶段的那一列的场景组成，其行号是前一阶段每簇中各点对应的编号，由此逐阶段输出叶节点![]()
的场景树以及相应的概率![]()
形成多阶段场景树。图3给出了一个简单的3阶段场景树生成示例。
![]()
图3 3阶段场景树约简示例
Fig.3 Example of three-stage scenario tree generation
SD方法是在最优运输方法基础上提出来的。最优传输方法可以推广到最优风光荷场景集的生成，即一个有限的风光荷场景集合I运输到有限的风光荷场景集合S，寻找超集I的最优子集S∗，使得2个风光荷场景集合的概率距离F(I,S∗)最小，由此获得最优的风光荷场景集S∗。
但原始风光荷场景集维度较高，采用传统的最优运输理论在计算时效率较低。因此，提出SD方法，通过将正则化系数γ和熵引入最优传输方法的目标函数，提高寻找最优风光荷场景集的计算效率，模型为
式中：πi,j为原始风光荷场景集I中的第i个场景向约简后的风光荷场景集S中第j个场景的转移概率；ci,j为原始风光荷场景集I中的第i个场景向约简后的风光荷场景集S中第j个场景的转移成本，通常使用2范数进行计算；ai和bj分别为2个风光荷场景集合中各场景的概率；2个风光荷场景集合的概率距离的最优目标值记为SD。采用迭代算法求解式（15），引入拉格朗日乘子将式（15）转化为一个无约束的问题，即式中：![]()
为拉格朗日乘子向量；Π为原始风光荷场景集中场景向约简后的风光荷场景集中场景的转移概率πi,j构成的矩阵，为一个N1×N2决策矩阵，其中N1和N2分别为原始风光荷场景集合和约简后的风光荷场景集合的大小。式中：M为原始风光荷场景集中场景向约简后的风光荷场景集中场景的转移成本ci,j和正则化系数构成的矩阵；exp(−αi/γ)和exp(−βj/γ)构成向量u和v，由此构成矩阵U和V，其相应迭代算法如下。式中：g为迭代次数，当g=1时，记v(g)中所有元素为1。与向量u和v相关的误差为向量u和v继续迭代计算，直至误差项εa和εb低于指定的误差阈值![]()
此时的u和v构成矩阵U和V，结合式（18）计算得到风光荷场景转移概率πi,j，由此得到SD为基于此，最优风光荷场景集生成（原始风光荷场景集I作为输入，将其约简生成一个包含r′个场景的最优子集S∗的算法流程如下。1）待约简的风光荷场景集为I，约简后场景集大小r′，误差阈值β2，初始化迭代次数，d′=1。2）定义风光荷子集，当d′=1时，S(d)是在超集中I随机选中的r′个场景，当迭代次数d⩾2时，S(d′)=S(d′−1)。3）使用式（18）~（20）计算I和S(d′)之间的Sinkhorn距离![]()
和最优运输计划矩阵Π(d′)。4）通过求运输矩阵Π(d′)每行最大值对应的列索引，可确定场景集I中各场景所在簇，由此围绕在S(d′)周围形成簇C1,C2,⋯,Cr′。5）计算相对误差，当d′=1时，继续步骤6），当d′⩾2时，计算相对误差![]()
6）如果相对误差大于误差阈值时，遍历簇内每个风光荷场景作为簇内中心计算运输距离，得到最小的簇内运输距离作为新的簇内中心，以此得到更新后的子集S∗:=S(d′)进行下一次迭代。7）当相对误差小于误差阈值时，得到最终约简子集S∗:=S(d′)，以及相应的聚类C1,C2,⋯,Cr′。8）设定S∗中所有场景的出现概率，由此计算最优子集S∗的概率为式中：Sp为Cr′各簇内场景概率和；Cp为原始风光荷场景集I中场景概率和。9）输出最终输出最优子集S∗包含r′个风光荷场景，各簇C1,C2,⋯,Cr′和对应的概率向量。本文采用PD衡量原始风光荷场景树和约简后场景树的相似性，数值越小表示计算精度越高，即式中：NT和![]()
分别为2个场景树的所有节点集；P1和P2分别为2个场景树的概率集合；符号f→j和e→i为中间节点f和e分别是叶节点j和i的前身；pi、![]()
pe和![]()
分别为到达节点i、j、e和f的概率。式（23）是一个线性规划问题，可采用分解的方式从一个阶段移动到另一个阶段进行求解。
某地2023年8760h风电功率、光伏功率和负荷标幺值组成风光荷全时序场景如图4所示。
![]()
Fig.4 Full temporal scenario set of wind-solar load
使用本文所提出的基于SSAE和改进AP聚类算法所得的聚类结果以及3个典型日的概率分别如图5所示和表1所示。
![]()
Fig.5 Clustering effect based on SSAE and improved AP
Table 1 Typical day probabilities![]()
从图5和表1可以看出，采用本文所提方法聚类出3个风电荷典型日，3个典型日包含的数据范围比较大，构成场景树的根节点。为进一步比较聚类效果的优劣，分别计算原始AP聚类、改进AP聚类、传统k-means聚类、PCA+改进AP聚类和基于SSAE+改进AP聚类的3个聚类有效性检验指标、迭代次数d和求解时间，结果如表2所示。
Table 2 Clustering results![]()
由表2可知，原始AP聚类算法得到23个聚类簇数，数量较多，且聚类中心的代表性较弱。改进AP聚类算法得到的最优聚类簇数为3，大大降低了聚类簇数；3个聚类指标中轮廓系数越大、CH指标越大、DB指标越小，簇自身越紧密、簇间越分散、聚类的效果越佳，结合表2发现，经过特征提取后的聚类效果优于特征提取前的结果，经SSAE特征提取后的聚类效果优于PCA，说明相比于PCA，SSAE具有更有效的提取场景特征的能力，提高改进AP算法的聚类质量；另外，从迭代次数和求解时间来看，相比于其他方法，经过SSAE特征提取后再进行聚类，可以有效降低聚类的迭代次数和求解时间，提高风光荷场景的聚类效率。根据3.1节生成的场景树根节点，假定负荷逐年增长，以0.22、0.55和0.23的概率表现为高速（4.5%）、中速（4%）和低速（2%）增长模式，逐步生成场景，获得一个5年的风光荷场景树，场景树的每个点代表一个典型日72维的风光荷场景。按本文所提方法对上述生成的大规模风光荷场景树进行削减，原始风光荷场景树的时间步长为1~5，将原始场景树构建一个3阶段的风光荷场景树。第3列和第5列对应于阶段2和阶段3。在生成场景树时，不考虑阶段之间的中间值或信息，只使用与阶段t对应的时间步τ上的值。阶段1对应3.1节生成的场景树初始场景，概率分别为0.2329、0.2658、0.5013；在任何阶段，需要约简的超集与前一阶段的节点数量一样，因此阶段2需要约简的集合为3个：![]()
分别对应阶段1中3个场景生成的在第3列的场景集，设置约简后子集的大小为2，基于SD方法将3个待约简集合分别约简为包含2个节点的子集，以及围绕这2个点形成的簇![]()
![]()
再由概率计算公式得到各簇概率。因此，在阶段3，有6个待约简的场景集，再按照上述步骤获得各簇和对应概率，由此获得一个3阶段风光荷场景树，每个节点代表一个72维风光荷场景，场景序号和概率如图6所示。
![]()
Fig.6 Three-stage wind-solar load scenario tree
图6生成的3阶段风光荷场景树的节点1、2、3代表的风光荷场景如图7所示，节点1、4、10代表的风光荷场景如图8所示。
![]()
Fig.7 Wind-solar load scenarios represented by nodes 1, 2 and 3![]()
Fig.8 Wind-solar load scenarios represented by nodes 1, 4 and 10
图7中场景1的风光出力水平最低，场景2的风光出力水平最高，3个场景的负荷水平接近。从各场景的变化趋势来看，光伏出力和负荷具有相近的波动趋势，光伏呈现单峰特性，负荷具有明显的双峰特性。风电的变化趋势各不相同，峰谷值可能出现在不同的时间段内，能够表征一定的差异性。图7中的3个场景作为场景树的根节点，覆盖的波动范围较大。由图8可以发现，负荷水平逐阶段增长，且增幅存在差异，反映出中长期负荷的增长情况。场景树中的每个节点代表的场景皆能反映风光荷的时序信息，较好地表征风电出力和负荷变化的随机性、波动性和差异性。此外，为评估本文所提的基于Sinkhorn削减场景树方法的优越性，采用基于Kantorovich距离和线性规划结合方法、k-means和线性规划结合方法和基于传统最优传输理论的线性规划方法作为对比方法，从削减的计算精度、计算效率2个方面对4种方法进行评价。4种方法计算得到的过程距离和总计算时间如表3所示。
Table 3 Process distance and time calculated by three methods![]()
从计算精度来看，过程距离为原始场景树和削减后场景树之间相似程度的测度，数值越小表示精度越高，本文所提方法的过程距离仅为0.0448，可见所提方法具有较高精度；从计算效率来看，本文所提方法的计算时间小于其他方法，用时0.1947s，计算效率相对更高。这是由于本文方法将其他方法中的线性规划求解问题替换为迭代算法计算过程，缓解其他方法对于风光荷维度较大的场景削减存在的内存瓶颈问题，因而本文方法在计算精度、效率方面有一定的提升。
本文从根节点到叶节点的完整路径定义场景树，提出了一种风光荷多阶段场景树的生成方法，经过算例验证，得到以下结论。1）本文提出的基于SSAE特征提取和改进AP算法选取风光荷典型日场景作为场景树根节点的方法，解决了原始AP聚类算法聚类簇数较大的问题，簇数由23降为3。此外，本文方法计算得到的3个聚类有效性检验指标较优，且经SSAE提取特征后求解时间仅需要17.2524s，优于其他方法，能够更有效提取场景的特征、提高聚类的质量和效率。2）本文采用基于Sinkhorn距离方法削减风光荷场景树，并引入过程距离来表征原始场景树和削减后场景树之间的相似度，数值越小表明计算精度越高。与其他方法相比，本文所提的场景树削减方法计算的过程距离为0.0448，用时0.1947s，证明了本文所提方法在计算的精度、时效方面兼具优势。本文提出的场景树生成方法可反映中长期风光出力和负荷发展的不确定性，具有一定可行性和优越性，适用于电力系统规划运行、调度分析、电网承载力评估和提升、随机模型预测控制等应用场景。
根据国家版权局最新规定，纸媒、网站、微博、微信公众号转载、摘编《中国电力》编辑部的作品，转载时要包含本微信号名称、二维码等关键信息，在文首注明《中国电力》原创。个人请按本微信原文转发、分享。欢迎大家转载分享。
