基于MVTV Copula方法的多风电场电力系统经济调度分析

李 奇 ，潘俞如 ，邱宜彬 ，陈维荣

（西南交通大学电气工程学院，四川 成都 610031）

就不包含风电场的电力系统而言，系统中引入的不稳定因素较少，此时的经济调度方案制定起来较为容易. 随着国家对风力产业的重视与相关研究的发展，引入电力系统中的风电场数量逐渐增加，系统中涉及的不稳定因素也随之增多，经济调度策略的制定变得越来越困难 然而不当的经济调度策略可能使电力系统运行成本增大，严重时甚至会引起电力事故，因此必须对接入风电场或其他新能源的[1].电力系统进行经济调度研究，这是确保电力系统稳定、经济运行的重要措施[2].

以往的经济调度分析主要侧重于风力发电的随机性，对风电出力相关性考虑较少. 近期有许多专家学者开始进行考虑风电出力相关性的经济调度研究[3]. 鉴于电力资产之间非线性的的动态相关结构，文献[4]采用对上、下尾部相关性敏感的时变SJC Copula （symmetrized joe clayton Copula）函数描绘电价序列间的相关性，以此拟合电价序列；为了最大程度上描述风电场机群的风速序列的相关性，建立由多种Copula函数组成的混合Copula风电场出力模型，相较于只采用某一种Copula函数构成的模型更切合实际[5]；文献[6]通过采用3种时变Copula函数分析两风电场出力数据并获得拟合模型，通过与静态Copula函数生成的出力模型对比，证明时变Copula函数对多风电场的相关性拟合精度更高；文献[7]采用Gumbel Copula函数建立考虑相关性的两风电场联合出力模型，并基于该模型采用柔性负荷调峰方法进行随机经济调度.

需要指出的是，以上文献均认为所描述的相关性模型在研究时段内是不变的. 而事实上，考虑时间波动给相关性带来的影响更符合实际情况，因此Patton提出了可以有效反映相关性内部变化的时变Copula[8]，其准确度在金融、股票市场等相关性研究的问题中均取得了较好的验证. 除此之外，大部分文献仅就两维风电出力数据进行研究，这并不能满足当下多风电场并网的现状. 实际中的多维风电出力分布存在随机性，且不同地区具有不同的特征，如若仍利用两维变量相关性分析方法分析多维风电出力数据的相关性，则模拟数据会具有明显的对称性，这与实际情况不符，偏离实际的数据又会在分析问题与应用时产生更大误差，后果严重者甚至危及电力系统稳定运行与设备安全[9]，因此有必要探索接入多个相关风电场的电力系统经济调度.

从此角度出发，本文提出MVTV Copula （mix vine time-varying Copula）方法对电力系统进行考虑多维风电出力相关性的经济调度研究. MVTV Copula方法是基于藤结构与动态、静态Copula建立的多维变量相关性分析方法，对数据分场景后得到对应的最优Copula函数，并采用混合藤结构建模得到一日内的典型随机场景. 为验证所述方法的有效性，采用IEEE-30节点系统，并结合我国某地区的风电历史出力数据进行研究，分别建立不考虑风电出力相关性的随机场景，以及依据所述MVTV Copula方法的考虑多风电场出力相关性的随机场景. 然后再以预测场景下的煤耗费用和随机场景下的再调度费用最小值为目标，利用回溯搜索算法（backtracking search algorithm，BSA）[10]求解目标函数约束内最优解及其对应调度方案，并通过考虑风电出力相关的经济调度目标函数值、不考虑风电出力相关的经济调度目标函数值与实际风电出力数据经济调度目标函数值的对比分析，对所提出经济调度方法与相关性建模方法的有效性进行验证.

1 Copula理论

1.1 概述

Copula函数由Sklar定理推导而来，Sklar认为，必然存在一个函数能将二维随机变量所对应的联合分布函数结合起来. 已知n维随机变量X=(X1,X2···,Xn)，F1,F2,···,Fn为各变量对应的边缘分布，F为所有变量的联合概率分布函数，则必有一个Copula函数C使得式（1）成立.

Copula函数对变量的边缘分布和相关性先后进行建模，在此基础上构造多元随机变量的联合概率分布，是一种可以简便、清晰地描述数据相关性的计算方法. 常见的Copula函数包含正态Copula函数、Guassian Copula、t-Copula和Gumbel Copula、Clayton Copula、Frank Copula. 由于每种Copula函数的参数性质与计算方式均不同，因此对应不同的相关性结构，这就要求在对二维变量相关结构进行建模时必须选择对应的最优Copula函数.

1.2 时变Copula函数

时变Copula函数在静态Copula函数的基础上考虑时间变化给相关性带来的影响，认为涉及参数服从与时间相关的动态方程，再通过相关系数与参数之间的映射关系，从而建立参数的动态计算方程.若以正态函数为例并记C为分布函数，u=(u1,u2,···,un)和v=(v1,v2,···,vn) 为相关变量，φ−1(•)表示正态分布的反函数，s、r为分布函数的相关变量，θ 、ρN,t分别为静态与时变相关系数，后者由对应的时变系数方程求得，其下标N代 表Normal Copula，t为时刻. 静态与时变Copula函数的具体表达式分别如式（2）、（3）所示.

1.3 藤Copula模型

仅单纯采用Copula函数法考虑多维风电出力之间的关系时，其相关性描述仅限于两维变量. 为打破数据维度对Copula函数的限制，有学者提出结合藤结构与Copula函数，采用pair-Copula对多维变量相关性进行分析[11]，该模型先拆分多维变量并通过Copula理论得到两两变量之间的相关性，再结合联合概率分布函数求解. 本文采用C藤（正则藤）和D藤分析数据结构. 当多变量中某个变量与余下变量均具有较强相关性，且整个数据的相关性结构呈现星型放射状时，C藤结构具有较高精度[12]；D藤模型则具有平铺式结构，可以更好地描述两两变量相关程度高的数据相关性[13]. 在进行藤Copula模型建模时，选择适应其相依结构的藤结构有助于提高相关性分析精度.

2 基于MVTV Copula方法的随机场景获取方法

受限于预测技术，风电出力还无法实现精确预测[14]，但实际情况下的风电有功出力数据符合风电出力的随机分布[15]，对于这一过程可以用随机场景来刻画. 本文为获取更贴近实际的多风电场出力随机场景，提出先利用混合藤结构与多种Copula函数构建多维风电出力模型，再在此基础之上利用蒙特卡洛抽样和场景缩减法获取典型随机场景. 所述基于MVTV Copula方法的随机场景获取方法具体如下：

1） 确定最佳聚类参数

统计历史数据中一天内同一时刻的风电出力数据，由此将一段时期内的历史风电出力划分成24组，利用K-means聚类法分析每组数据，并依据中心值进行分场景. 针对传统K-means方法聚类参数需要根据经验值人为设定，不能根据数据的不同自适应选取最优聚类参数的问题，本文利用基于密度的聚类数确定（density based index，DBI）法[16]计算几种聚类数对应的DBI值（D），以此指导K-means聚类方法选取最优聚类参数，如式（4）所示.

式中：m函数总数；B为类边缘密度；C为类中心密度；Dc为第c类聚类数对应的DBI指标.

然后再根据最佳聚类数对待分析数据进行聚类划分.

2） 确定最优Copula函数

以第一组数据（即历史风电出力数据每日时刻01：00的风电出力数据）为例，分别选择各组聚类结果对应的最优Copula，其中静态Copula函数的最优参数均由最大概似估计确定（具体实施可以使用MATLAB自带copulafit函数完成，在此不做赘述），时变最优参数由对应的时变方程确定[17]. 最优函数类型则根据赤池信息（Akaikein formation criterion，AIC）法计算各备选Copula函数和经验Copula函数的AIC值，具体计算方式如式（5）.

式中：k为拟合模型的参数个数，其大小反映了模型的复杂程度；L为似然函数，其大小反映了模型拟合度.

时变Copula函数对应的复杂程度虽然较高，但AIC法会引入惩罚系数提高拟合度. 从多个拟合模型中判断最优模型时，一般较小AIC值对应的模型对应较高的拟合精度.

3） 获取采样点

基于以上理论基础，多维变量的采样点获取方式可以概括如下：

步骤1利用MATLAB自带unifrnd函数模拟生成1组M×N（数据个数 × 维数）均匀变量z，z=(z1,z2, ···,zN)；

步骤2令第1维待求变量u1等于步骤1中的第1列，即u1=z1；

步骤3根据

第2维待求变u2可利用式（7）来计算.

式（6）、（7）中：g、h分别为样本数据中不同维度的采样点向量；hj为h中的第j个变量；h−j为h去除hj后的向量；Xi= (ui,vi).

计算时z2为步骤1生成，C(u1,u2)表示u1、u2在对应Copula函数处的函数值，可以通过copulacdf函数计算. 对等式（7）右边求解，加一个while循环，循环内采用二分法求解u2，直至等式右边计算出来的数据和左边相差很小时结束，此时随机出来的u2即为所求模拟数据.

步骤4令Zi为Xi对应的采样点，则利用

进行迭代求解后续维数模拟数据.

4） 确定最佳藤结构本文选取CvM （cramer von mise）距离（dCvM）来作为模型模拟效果评判指标，通过对比不同藤结构建模结果对应指标值来判断藤结构模拟效果的好坏，从而完成最优藤结构的选取工作，具体实施方式如式（8）所示[17].

式中：Ca(•)与Ce(•)分别为各备选Copula函数与经验Copula函数的联合分布函数值.

最后选取最小的CvM距离所对应的藤结构为该类数据的最优藤结构.

5） 获取典型随机场景

上述第 4）步所确定的最优藤结构对应的随机场景即为所求随机场景，但因其数量一般较多，所以有必要对所获取的随机场景进行相应场景削减以获取典型随机场景，本文采用根据距离进行场景削减的同步回代法完成随机场景的获取[18]：

步骤1对待分析对象进行场景初始化，利用场景生成方法随机抽取S个场景，其中各生成场景s（s∈S）相应的概率为

步骤2获取两两场景的欧氏距离，具体表达为

步骤3利用概率与欧式距离进行场景缩减步骤，获取与场景s相距最近的场景s′，并使用s′代替s，概率距离计算方式如式（9）.

式中：Ds,s′为场景s和s′的欧氏距离；J为所有被削减的场景的集合，其概率修正方式为

式中：ps′为场景s′ 的概率.

步骤4重复步骤2与步骤3，削减至预定的保留数目即可停止.

重复进行以上建模步骤24次，将24次典型随机场景进行组合，即可获取一日内的典型随机场景.

3 考虑风电出力相关性的电力系统经济调度模型

3.1 目标函数

鉴于风电场正常运行期间所需的运行费用较低[19]，本文在建立电力系统经济调度目标函数时忽略风电场的运行成本，仅考虑常规的火电机组（以下简称机组）在运行期间的煤耗费用以及在随机场景下的再调度燃料费用作为目标函数（f），如式（11）.

式中：fy为电力系统在预测风电场出力场景对应的机组燃料成本；fw为随机场景对应的电力系统再调度费用.

对于fy来说，其燃料成本与机组的发电量呈二次函数关系，具体函数形式如式（12）.

式中：T为一个调度周期内总共包含的时间；Q为电力系 统 中包 括 的机 组 的个数；PGq(t) 为 在时刻t机组q的输出功率；aq、bq、cq为机组q的煤耗费用价格系数，由机组本身性质决定.

对于fw，其具体表达式为

式中：hs与 βq分别为各个随机场景对应的概率与机组q的再调度费用；PGq,s(t) 为场景s下时刻t机组q的功率.

3.2 约束条件

在进行经济调度研究时，必须考虑实际情况对经济调度模型的约束限制，只有这样才能获取符合实际情况的调度结果. 除了机组出力限值约束和机组爬坡约束外，还包括系统功率平衡约束. 各约束条件表达如式（14）[20-21].

式中：NW为系统中接入的风电场的总数；PWnW(t)为第nW个风电场在时刻t输出的有功功率；PLoad(t)以及PLoss(t)分别为系统在时刻t的系统负载以及系统网损；PGq,min、PGq,max分别为机组q的有功功率下限、上限；Pq,t和Pq,t−1分别为机组q于时刻t和时刻t−1发出的有功功率；Pq,up和Pq,down分别为最大上升功率和最大下降功率.

4 算例分析

为验证所述基于MVTV Copula方法的电力系统经济调度方式的有效性，本文借助IEEE30节点系统对所述方法进行验证，风电出力原始样本取自我国某地域3个毗邻风电场的历史数据. 其中，风电接入节点分别为8、16、21节点，常规机组接入节点分别为1、2、13、22、23、27节点，调度间隔取1 h. 所述IEEE30节点的拓扑图如图1所示[22-23].

图1 仿真系统拓扑Fig. 1 Topology of simulation system

因实际运行中负荷曲线一般较为固定，所以算例中不考虑负载曲线的预测，而仅对调度周期内的风电场出力数据进行预测. 算例中使用2012年4月某天的风电实际出力及预测出力为数据进行经济调度分析，各风电场实际出力曲线及预测出力曲线如图2所示.

鉴于同一地区气候条件一般不会发生较大改变，所以本文利用该地区历史同期数据来构建2012年4月某天不同时刻的风电出力MVTV Copula模型.以01：00时刻风电出力数据的MVTV Copula模型为例，首先利用DBI指标获取各模型的最佳分类数目，具体操作过程中将K-means聚类参数分别设置为2～14，并选取DBI指标最大的聚类方案确定多维数据的最佳聚类数，所得结果如图3所示.

从图3可以看出：根据DBI指标所确定的最佳聚类数为5类，据此对待分析数据各聚类结果分别进行最优藤结构及最优Copula判别，备选静态Copula函数有：Gumbel Copula函数、正态Copula函数、Frank Copula函数、Clayton Copula函数与t-Coopula函数；备选时变Copula函数有：时变Gumbel Copula函数、时变正态Copula函数与时变SJC Copula（symmetrised Joe Clayton Copula）函数. 所得各参数详细情况如表1所示.

图3 DBI指标与场景数关系Fig. 3 Relationship between DBI index and scenario number

表1 各场景下最优藤结构及最优参数Tab. 1 Optimal vine structure and optimal parameters for different scenarios

同时，为了更加直观地反映混合藤结构能在模拟多风电场相关性建模中相较孤立的藤结构获取更高的准确性，本文仍采用AIC值作为模型与原始数据接近程度的衡量标准，采用混合藤与两种孤立藤结构的AIC值对比，如表2所示.

表2 采用混合藤与孤立藤结构的AIC值对比Tab. 2 AIC values of mix vine and isolated vine structures

在获取最佳聚类数和最优Copula函数的基础之上，利用采样方法获取24 h内不同时刻风电出力的720个随机场景，此时按照同步回代削减法可获取若干典型随机出力场景. 本文中典型随机场景的数量定为10个，具体分布如图4所示.

图4 基于MVTV Copula模型的典型随机场景Fig. 4 Typical random scenarios based on MVTV Copula model

由图4可以看出：基于MVTV Copula方法生成的多风电场拟合数据基本分布在坐标轴对角线附近，说明所模拟的风电数据表现出同增同减的正相关特性，该分布较为符合实测相邻风电场实测出力数据情况. 而图5所示为不考虑风电出力相关性的随机场景，风电数据则分布较为零落，未表现出明显分布特征，与实际相邻风电出力数据分布差别较大.根据以上两类典型随机场景进行电力系统经济调度研究，采用BSA算法进行求解，分别计算相应的再调度费用，如表3所示.

由表3可以看出：虽然不考虑风电出力相关性的系统再调度费用估计值小于考虑风电数据相关性的经济调度方式，前者仅为实际调度费用的43.6%，后者为实际调度费用的82.0%，但不考虑相关性所得结果却与实际情况存在较大偏差，而考虑相关性则较为接近实际调度结果. 这是由于实际的风电数据之间存有一定正相关性，某一个风电场的出力数据增大将大大提高其余风电场的出力数据增大的概率，考虑风电出力相关性的经济调度方式一定程度上能够反映这种相关性，得出更贴合实际状况的分析结果. 考虑风电出力相关性时各典型随机场景下不同时段常规机组出力和风电出力总和曲线以及负荷曲线如图6所示.

图5 不考虑相关性的典型随机场景Fig. 5 Typical random scenarios ignoring dependence of wind outputs

表 3 两类典型随机场景情况调度费用对比Tab. 3 Economic dispatch results of two typical random scenarios USD

图6 负载曲线及各场景下机组和风电总出力曲线Fig. 6 Load curve and total output curve of unit and wind power output under different scenarios

从图6中可以看出：对于所述10个典型随机场景，利用BSA算法所获取的考虑风电出力相关性的常规机组有功出力方案，在和风电场有功出力的共同作用下，其总出力曲线可以很好地跟踪负荷曲线，说明系统的有功出力可以满足负荷以及网损需求，所得调度策略在满足经济最优的同时也满足系统约束，策略具备一定实际参考价值.

5 结 论

针对涉及风电出力相关性的经济调度分析忽略多维风电场出力相关性这一问题，通过建立基于MVTV Copula方法的多风电出力模型对电力系统经济调度进行研究. 该经济调度模型先对历史数据采用K-means聚类分析，并利用MVTV Copula方法采样生成随机场景，将预测风电出力场景对应的常规机组运行成本和随机场景对应的电力系统再调度费用最低做为调度目标，利用BSA进行经济调度问题的求解. 为验证所述方法的有效性，采用IEEE30节点系统，并结合我国某地3个相邻风电场的历史数据对所述方法进行仿真验证，所得结果表明，根据MVTV Copula模型模拟获得的风电数据所呈现的分布特性更贴合实测数据情况，并且计及相关性的再调度费用与实际再调度费用的更为接近，进一步说明多风电场相关性在电力系统研究中的重要性. 同时，调度人员可以根据所得随机场景制定合理调度计划，减小风电出力相关性对系统安全经济运行的危害.