一种含高比例风电并网的自适应线性化概率潮流计算方法与流程

本发明属于电力系统运行分析和控制技术领域，涉及一种含高比例风电并网的自适应线性化概率潮流计算方法。

背景技术：

高比例风电并网将成为智能电网下的重要场景，其大范围波动导致电力系统不确定性进一步加剧。与此同时，“源-网-荷”互动技术的发展使得柔性负荷能够通过需求响应(demandresponse,dr)主动参与源-荷互动，用户用电自主性增强，负荷侧不确定性显著增加。高比例风电并网场景下源-荷强随机性将给系统安全运行带来极大的挑战，因此，定量评估源-荷双侧不确定性对系统运行的影响具有现实意义。

概率潮流(probabilisticpowerflow,ppf)计算是电力系统不确定性分析的有效手段，准确构建源-荷双侧不确定性模型是ppf计算分析的基础。针对风电出力和负荷概率模型，目前大多采用双参数威布尔分布和正态分布对风速及负荷分布进行拟合，模型简单且实现方便。随着高比例风电并网以及负荷峰谷差进一步拉大，传统分布无法准确拟合多峰、不对称等强随机性特征。此外，柔性负荷可以通过dr主动参与源-荷互动，原有负荷曲线形态发生改变。目前，已有学者在风电消纳和发用电协同调度层面考虑了柔性负荷互动响应，但多数基于响应确定性模型，且并未深入探究互动响应对系统潮流分布的影响，事实上，柔性负荷互动响应量亦具有较强的不确定性。鲜有同时考虑源-荷双侧强随机扰动对潮流分布的影响。

如何有效求解计及源-荷强随机性的ppf是一个重要问题，常规ppf计算方法在处理强随机性时难以兼顾计算效率和精度，无法适用于高比例风电强随机波动场景。此外，尚未计及输入变量间客观相关性影响，ppf模型亦无法有效处理相关性，忽略相关性将导致分析结果产生较大偏差。因此，亟需寻求能够同时有效处理源-荷强随机性和相关性的ppf计算方法。

技术实现要素：

发明目的：本发明的目的在于针对现有技术的不足，提出了一种含高比例风电并网的自适应线性化概率潮流计算方法，解决多区域线性化模型对主观因素依赖性高的问题，同时能够解决计及风电出力和负荷强随机性和相关性影响下的电力系统潮流计算，获取更加准确的概率分布特性，适用于高比例风电并网下概率风险分析与评估。

技术方案：本发明提供了一种含高比例风电并网的自适应线性化概率潮流计算方法，包括以下如下步骤：

步骤1：采用高斯混合模型统一刻画源-荷强随机性，考虑源-荷互动随机响应，准确建立源-荷不确定性模型；

步骤2：在步骤1的基础上随机模拟生成输入变量相关性样本；

步骤3：对潮流方程进行自适应线性化处理，采用迭代算法实现随机波动范围自适应多区域划分，并在划分区域内部线性化以减小潮流全局线性化误差；

步骤4：在步骤3的基础上，利用自适应线性化半不变量法的解析法获取区域整合重组后状态变量的各阶半不变量，最后，运用c型gram-charlier级数拟合状态变量概率分布；

进一步，所述步骤1包括以下步骤：

步骤101：直接根据风功率和负荷样本，采用高斯混合模型(gaussianmixturemodel，gmm)拟合两者的概率分布，建立风功率和负荷的统一概率模型；对于数据缺失的输入变量，则采用相应的经验分布进行拟合，高斯混合模型通常由多个高斯子成分进行线性加权，概率密度函数如下：

式中，x表示随机变量，n表示子成分个数，wt、μt和σt分别表示高斯混合模型的第t个子成分的权重、期望以及标准差，fn(·)表示正态分布的概率密度函数；

其中权重系数wt满足以下归一化条件：

步骤102：采用极大似然估计法确定各子成分的参数ωt、μt和σt。构造对数似然函数l：

式中，nl表示输入变量xi的样本容量；

步骤103：采用期望最大化算法得到原始样本xi(i＝1,2,…,nl)隶属于第t个子成分的概率p(i,t)为：

进一步得到各子成分权重、期望及方差的估计值：

对上述步骤不断进行迭代，直至满足预设精度。

步骤104：考虑源-荷互动随机响应，柔性负荷响应量取决于源-荷互动调度，定义调度系数v代表负荷响应灵活性，此处更侧重于单个时间断面下的潮流分析，因而v直接被指定。负荷互动响应量pdr如下：

pdr＝vpl；

式中，pdr表示柔性负荷互动响应量，pl表示柔性负荷互动响应量；

步骤105：柔性负荷互动响应过程中，实际负荷响应量具有较强的不确定性。对负荷侧可调度资源进行分析，采用正态分布描述其概率特征：

式中，f(pdr)表示柔性负荷随机响应的概率密度函数，分别表示柔性负荷互动响应量期望和标准差；

为进一步量化随机互动响应量，定义随机主动系数s刻画：

进一步：所述步骤2包括以下步骤：

步骤201：实际输入变量之间具有一定的相关性，根据原始样本得到输入变量x和y相关系数：

式中，分别表示输入变量x、y的样本均值；

步骤202：根据输入变量高斯混合模型(gmm)拟合得到各子成分参数(wt,μt,σt)，随机生成各个正态分布子成分样本，采样总体规模为n。依次对n个随机变量进行采样，形成n×n阶采样矩阵。同时，根据顺序矩阵ls调整输入变量样本矩阵的排列顺序，使其相关系数基本和原始样本ρx,y保持一致，得到最终的输入变量样本矩阵xn×n。

进一步：所述步骤3包括以下步骤：

步骤301：采用全线性化潮流模型，在基准运行点处将潮流方程中电压二次项进行泰勒级数展开，忽略2阶及以上高阶项，可得到如下近似关系：

式中，θij、vi分别表示节点i-j电压相角差、节点i电压幅值，下标k-1表示变量第k-1次迭代结果；

进一步代入潮流方程，得到线性化潮流方程，可简化成如下表达式：

式中，x、z、w分别表示节点电压、支路潮流及节点注入功率；s0和t0表示灵敏度矩阵，其中j0为雅可比矩阵，下标0表示基准运行点；

步骤302：考虑源-荷双侧随机扰动，针对波动范围内不同情形选择不同的线性化区域步长，取状态变量期望值μx为线性化精度衡量指标，定义自适应迭代收敛判据为连续两次迭代所得μx误差绝对值的最大值达到精度要求ε。选择随机波动上下限wf,min～wf,max为初始线性化区域，区域注入功率期望值为ef1，在该区域内采用常规半不变量法计算ppf得到

步骤303：将步骤302所得线性化区域对分，波动范围划分为子区域wf,min～wf1、wf1～wf,max，其中wf1满足：

子区域对应注入功率期望值为ef11、ef12，同理在相应子区域内分别进行半不变量法概率潮流计算得到衡量指标

步骤304：判断本次区域划分后结果是否满足下式，如果满足，则迭代停止；否则继续进行下一步；

步骤305：将上一次迭代的波动范围子区域进一步划分为wf,min～wf2、wf2～wf1与wf1～wf3、wf3～wf,max，其中wf2、wf3满足：

子区域对应注入功率期望值为(ef21、ef22)与(ef31、ef32)，并得到

步骤306：针对子区域wf,min～wf1、wf1～wf,max，考察下式是否成立，对满足精度要求的子区域不在细分，不满足的子区域需要重复上述步骤：

步骤307：不断将上一次线性化ppf计算迭代结果与下一次对应区域划分后的结果进行对比，直到所有的子区域划分结果均满足对应的误差精度要求。在自适应区域划分的基础上，将潮流方程在各区域内部基准运行点处线性化展开，得到第m个线性化区域表达式：

δxm＝smδwm；

式中，sm表示在第m个线性化区域进行潮流计算所得灵敏度矩阵。

进一步：所述步骤4包括以下步骤：

步骤401：引入风功率多区域相关性。假设区域m下风功率相关系数矩阵为cw,m，且为正定矩阵：

式中，ρij,m为第i和j座风电场在区域m下的风功率相关系数；

步骤402：为快速处理相关性，将功率变化量δw进行标准化，在步骤307的基础上重新得到区域m下的线性化表达式：

式中，表示注入功率变化量δw(m)的标准化向量；σ(m)为δw(m)标准差对应的方阵；

步骤403：进一步地，对相关系数矩阵cw,m进行cholesky分解：

式中，gw,m为分解得到的下三角矩阵；

步骤404：将区域m下标准化后随机变量转化为互不相关的变量δy(m)的组合：

步骤405：考虑σ(m)及gw.m的影响，得到新的灵敏度矩阵s1(m)、t1(m)以及线性化表达式：

步骤406：为求得状态变量的概率分布，对各区域进行整合，结合半不变量的定义，通过相应中心矩求得状态变量各阶半不变量，区域划分整合后随机变量δx的k阶中心矩αk表示如下：

式中，r表示最终自适应区域划分数目；

步骤407：将输入变量样本δy离散化，进一步得到：

式中，δym,i表示区域m下注入功率第i个样本；pm,i表示δym,i发生的概率；nm表示区域m样本总量。

步骤408：由半不变量与中心矩的关系，得到状态变量各阶半不变量κ1、…、κr+1：

对于服从正态分布的随机变量，采用常规数值解析方法求取其半不变量；对于服从非正态分布的相关性输入变量，进行自适应线性化处理并处理区域相关性，在此基础上采用自适应线性化半不变量法的解析法获取状态量各阶半不变量，将上述所得状态量半不变量相加得到待求状态变量的各阶半不变量；

步骤409：采用c型gram-charlier级数拟合得到状态量的概率分布。

工作原理：本发明首先采用高斯混合模型统一刻画风电和负荷强随机性，同时考虑柔性负荷随机响应，准确建立源-荷不确定性模型。然后，随机模拟生成输入变量相关性样本。接着，对潮流方程进行自适应线性化处理，采用迭代算法实现随机波动范围自适应多区域划分，并在划分区域内部线性化以减小潮流全局线性化误差，在此基础上，利用自适应线性化半不变量法的解析法获取区域整合重组后状态变量的各阶半不变量。最后，运用c型gram-charlier级数拟合状态变量概率分布。

有益效果：与现有含风电并网的ppf计算相比，本发明具有如下优点和技术效果：

(1)能够有效处理源-荷强随机性和相关性，同时兼顾ppf计算精度和效率，获取更加准确的概率分布特性；

(2)考虑柔性负荷随机互动响应，可以有效促进高比例风电的消纳，降低系统越界概率；

(3)能够通过迭代算法适应不同的波动剧烈程度，合理实现自适应多区域划分，有效解决多区域线性化模型对主观因素依赖性高的问题。

附图说明

图1为改进ieee30节点系统结构图；

图2含高比例风电并网的自适应线性化ppf计算流程图；

图3为实际风功率gmm拟合图；

图4为节点21电压幅值的概率分布曲线；

图5为柔性负荷互动响应对节点21电压幅值的影响。

具体实施方式

以下结合附图和实例对本发明的实施作进一步说明，但本发明的实施和包含不限于此。

一种含高比例风电并网的自适应线性化概率潮流计算方法，包括以下如下步骤：

步骤1：采用高斯混合模型统一刻画源-荷强随机性，考虑源-荷互动随机响应，准确建立源-荷不确定性模型；

步骤2：在步骤1的基础上随机模拟生成输入变量相关性样本；

步骤3：对潮流方程进行自适应线性化处理，采用迭代算法实现随机波动范围自适应多区域划分，并在划分区域内部线性化以减小潮流全局线性化误差；

步骤4：在步骤3的基础上，利用自适应线性化半不变量法的解析法获取区域整合重组后状态变量的各阶半不变量，最后，运用c型gram-charlier级数拟合状态变量概率分布；

进一步，所述步骤1包括以下步骤：

步骤101：直接根据风功率和负荷样本，采用高斯混合模型(gaussianmixturemodel,gmm)拟合两者的概率分布，建立风功率和负荷的统一概率模型；对于数据缺失的输入变量，则采用相应的经验分布进行拟合，gmm(高斯混合模型)通常由多个高斯子成分进行线性加权，概率密度函数如下：

式中，x表示随机变量，n表示子成分个数，wt、μt和σt分别表示高斯混合模型的第t个子成分的权重、期望以及标准差，fn(·)表示正态分布的概率密度函数；

其中权重系数wt满足以下归一化条件：

步骤102：采用极大似然估计法确定各子成分的参数ωt、μt和σt。构造对数似然函数l：

式中，nl表示输入变量xi的样本容量；

步骤103：采用期望最大化算法得到原始样本xi(i＝1,2,…,nl)隶属于第t个子成分的概率p(i,t)为：

进一步得到各子成分权重、期望及方差的估计值：

对上述步骤不断进行迭代，直至满足预设精度。

步骤104：考虑源-荷互动随机响应，柔性负荷响应量取决于源-荷互动调度，定义调度系数v代表负荷响应灵活性，此处更侧重于单个时间断面下的潮流分析，因而v直接被指定。负荷互动响应量pdr如下：

pdr＝vpl；

式中，pdr表示柔性负荷互动响应量，pl表示柔性负荷互动响应量；

步骤105：柔性负荷互动响应过程中，实际负荷响应量具有较强的不确定性。对负荷侧可调度资源进行分析，采用正态分布描述其概率特征：

式中，f(pdr)表示柔性负荷随机响应的概率密度函数，分别表示柔性负荷互动响应量期望和标准差；

为进一步量化随机互动响应量，定义随机主动系数s刻画：

进一步：所述步骤2包括以下步骤：

步骤201：实际输入变量之间具有一定的相关性，根据原始样本得到输入变量x和y相关系数：

式中，分别表示输入变量x、y的样本均值；

步骤202：根据输入变量高斯混合模型(gmm)拟合得到各子成分参数(wt,μt,σt)，随机生成各个正态分布子成分样本，采样总体规模为n。依次对n个随机变量进行采样，形成n×n阶采样矩阵。同时，根据顺序矩阵ls调整输入变量样本矩阵的排列顺序，使其相关系数基本和原始样本ρx,y保持一致，得到最终的输入变量样本矩阵xn×n。

进一步：所述步骤3包括以下步骤：

步骤301：采用全线性化潮流模型，在基准运行点处将潮流方程中电压二次项进行泰勒级数展开，忽略2阶及以上高阶项，可得到如下近似关系：

式中，θij、vi分别表示节点i-j电压相角差、节点i电压幅值，下标k-1表示变量第k-1次迭代结果；

进一步代入潮流方程，得到线性化潮流方程，可简化成如下表达式：

式中，x、z、w分别表示节点电压、支路潮流及节点注入功率；s0和t0表示灵敏度矩阵，其中j0为雅可比矩阵，下标0表示基准运行点；

步骤302：考虑源-荷双侧随机扰动，针对波动范围内不同情形选择不同的线性化区域步长，取状态变量期望值μx为线性化精度衡量指标，定义自适应迭代收敛判据为连续两次迭代所得μx误差绝对值的最大值达到精度要求ε。选择随机波动上下限wf,min～wf,max为初始线性化区域，区域注入功率期望值为ef1，在该区域内采用常规半不变量法计算ppf得到

步骤303：将步骤302所得线性化区域对分，波动范围划分为子区域wf,min～wf1、wf1～wf,max，其中wf1满足：

子区域对应注入功率期望值为ef11、ef12，同理在相应子区域内分别进行半不变量法ppf计算得到衡量指标

步骤304：判断本次区域划分后结果是否满足下式，如果满足，则迭代停止；否则继续进行下一步；

步骤305：将上一次迭代的波动范围子区域进一步划分为wf,min～wf2、wf2～wf1与wf1～wf3、wf3～wf,max，其中wf2、wf3满足：

子区域对应注入功率期望值为(ef21、ef22)与(ef31、ef32)，并得到

步骤306：针对子区域wf,min～wf1、wf1～wf,max，考察下式是否成立，对满足精度要求的子区域不在细分，不满足的子区域需要重复上述步骤：

步骤307：不断将上一次线性化ppf计算迭代结果与下一次对应区域划分后的结果进行对比，直到所有的子区域划分结果均满足对应的误差精度要求。在自适应区域划分的基础上，将潮流方程在各区域内部基准运行点处线性化展开，得到第m个线性化区域表达式：

δxm＝smδwm；

式中，sm表示在第m个线性化区域进行潮流计算所得灵敏度矩阵。

进一步：所述步骤4包括以下步骤：

步骤401：引入风功率多区域相关性。假设区域m下风功率相关系数矩阵为cw,m，且为正定矩阵：

式中，ρij,m为第i和j座风电场在区域m下的风功率相关系数；

步骤402：为方便处理相关性，将功率变化量δw进行标准化，在步骤307的基础上重新得到区域m下的线性化表达式：

式中，表示注入功率变化量δw(m)的标准化向量；σ(m)为δw(m)标准差对应的方阵；

步骤403：进一步地，对相关系数矩阵cw,m进行cholesky分解：

式中，gw，m为分解得到的下三角矩阵；

步骤404：将区域m下标准化后随机变量转化为互不相关的变量δy(m)的组合：

步骤405：考虑σ(m)及gw.m的影响，得到新的灵敏度矩阵s1(m)、t1(m)以及线性化表达式：

步骤406：为求得状态变量的概率分布，对各区域进行整合，结合半不变量的定义，通过相应中心矩求得状态变量各阶半不变量，区域划分整合后随机变量δx的k阶中心矩αk表示如下：

式中，r表示最终自适应区域划分数目；

步骤407：将输入变量样本δy离散化，进一步得到：

式中，δym,i表示区域m下注入功率第i个样本；pm,i表示δym,i发生的概率；nm表示区域m样本总量。

步骤408：由半不变量与中心矩的关系，得到状态变量各阶半不变量κ1、…、κr+1：

对于服从正态分布的随机变量，采用常规数值解析方法求取其半不变量；对于服从非正态分布的相关性输入变量，进行自适应线性化处理并处理区域相关性，在此基础上采用自适应线性化半不变量法的解析法获取状态量各阶半不变量，将上述所得状态量半不变量相加得到待求状态变量的各阶半不变量；

步骤409：采用c型gram-charlier级数拟合得到状态量的概率分布。

本发明首先采用高斯混合模型统一刻画风电和负荷强随机性，同时考虑柔性负荷随机响应，准确建立源-荷不确定性模型，然后随机模拟生成输入变量相关性样本，接着对潮流方程进行自适应线性化处理，采用迭代算法实现随机波动范围自适应多区域划分，并在划分区域内部线性化以减小潮流全局线性化误差，在此基础上，利用自适应线性化半不变量法的解析法获取区域整合重组后状态变量的各阶半不变量，最后运用c型gram-charlier级数拟合状态变量概率分布。

算例一：

选取ieee30节点系统进行算例测试，在节点21、24上接入两个额定容量为60mw的风电场，风电渗透率水平约为40％，功率因数恒定为0.98，并将节点19、30处的负荷替换为实际负荷，系统结构示意图如图1所示；同时假定节点7、12、21为柔性负荷，可参与源-荷互动场景，负荷响应调度系数v取0.5，节点7、12、21的随机主动系数分别设置为0.20、0.12以及0.16。

采用gmm(高斯混合模型)拟合实际风功率及负荷概率分布，权衡模型精度和效率，选取合适高斯子成分个数并得到各子成分参数，如表1所示。

表1为风电场参数

在此基础上，生成样本容量为30000的风功率样本和负荷样本，实际风功率之间的相关系数为0.8472，暂不考虑负荷间的相关性；同时将得到的模拟样本与实际样本进行对比，两者期望值和标准差的相对误差均在10^-4以下，pdf曲线总体拟合良好，如图3所示，表明gmm能够准确刻画非正态输入变量强随机性。此外，假定数据缺失的负荷均服从正态分布，其余节点负荷以标准算例中给定负荷数据作为期望值，标准差取期望值的10％。

采用迭代算法对风功率及非正态负荷样本波动范围进行自适应多区域划分，选取电压幅值期望值为衡量指标，迭代误差要求ε为0.02％；风功率样本自适应划分为5个线性化区域(pu)，分别为[0，0.125]、[0.125，0.25]、[0.25，0.50]、[0.50，0.75]、[0.75，1.00]，同理，非正态负荷样本划分为2个区域(pu)，分别为[0.00，0.50]、[0.50，1.00]，采用自适应线性化半不变量法进行ppf计算，为了验证算法的准确程度，以30000次蒙特卡罗模拟法(montecarlosimulationmethod,mcsm)计算结果为基准，同时将常规半不变量法与本计算方法所得ppf计算结果进行对比。

采用状态变量数字特征的相对误差指标度量本计算方法的准确程度，可以得到采用常规半不变量法所得状态变量数字特征的相对误差指标普遍较大，其相对误差均值和最大值最大分别为29.4027％和56.2648％，而采用本计算方法相应结果为1.9775％和7.0339％，均在8％以内，计算精度得到明显提高。另外，采用mcsm、常规半不变量法及本计算方法进行ppf计算，平均耗时分别为1072.30s、7.87s和16.18s，由于本计算方法进行自适应线性化处理，计算平均耗时比常规半不变量法稍多，但计算效率明显优于mcsm，最终结果表明，本计算方法适用于含高比例风电并网下源-荷不确定性分析，能兼顾计算精度和效率，具有一定的准确性和快速性。

通过c型gram-charlier级数拟合得到状态变量概率分布，以节点21电压幅值为例进行说明，如图4所示。

由图4可知，采用本计算方法与mcsm所得累积概率分布函数(cumulativedistributionfunction,cdf)更为接近，节点电压幅值越界(>1.05pu)概率分别为22.24％和22.20％，能够有效保留ppf尾部效应特性；而采用常规半不变量法所得cdf与之相差较大，且电压越界概率仅为12.51％，低估了系统潜在运行风险，导致ppf分析结果偏离客观实际，因此，进一步说明本计算方法能够获取更加准确的概率特性，为运行调度人员的安全分析提供有效信息，从而提前在系统运行薄弱环节做好预防措施。

为验证本计算方法能够更加有效地适应高比例风电并网场景，在未来可预见的风电占比大小前提下，通过改变风电并网渗透率，分别利用mcsm、常规cm以及本自适应线性化方法进行ppf计算，首先假定其余参数与上述3.2小节保持一致，然后采用方差和根均值(averagerootmeansquare,arms)度量ppf计算线性化误差的大小，得到各自区域划分数、计算时间、节点21电压幅值的arms以及电压越界概率值，结果如表2所示：

表2为不同风电渗透率下的三种方法所得结果对比

由表2可得，随着风电渗透率的不断提高，常规cm的arms值由0.110上升至0.550，准确性下降明显，而采用本计算方法能够根据风功率样本波动程度合理进行自适应区域划分，对应区域划分数目分别为1、2、5、10个，所得arms值普遍较小，有效弥补了高比例风电并网下ppf计算线性化精度较差的明显不足，避免了区域多划分造成的时间负担(或少划分时引起的精度不足)。

高比例风电并网后的大范围波动可能导致电力系统的网络安全约束不在允许范围内，节点21电压幅值的cdf随调度系数v变化曲线如图5所示。

由图5可知，当v逐渐提高时，表示将负荷从其余时段转移到该时段或者用户主动增加负荷量，促进了风电的消纳，电压越界概率大大降低；当v提高到0.8时，电压越界的概率降至约10％，由此可见，柔性负荷参与互动响应，可有效降低系统越界概率，减小系统运行风险。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明保护的范围之内。