CN112217215B - 一种基于psd-bpa的大型电力系统随机潮流计算方法 - Google Patents

Abstract

本发明公开了一种基于PSD‑BPA的大型电力系统随机潮流计算方法，基于历史数据，对风速和光照强度进行概率建模和相关性建模；基于模型采样生成风速和光照强度数据，并转化为出力场景集合作为输入；分析系统运行方式数据，根据每个采样场景分别修改对应运行参数，调用PSD‑BPA程序完成计算，读取并记录各场景下的潮流计算结果；读取随机潮流计算结果，绘制经验累积分布函数和概率密度函数；计算节点电压和支路潮流越限的概率，分析系统局部或整体潮流风险，若概率大于允许值，则认为节点/支路存在风险。本发明开发对PSD‑BPA软件实现了有效的调用与数据对接，能够对大型电力系统实现准确、高效的蒙特卡洛随机潮流计算。

Description

一种基于PSD-BPA的大型电力系统随机潮流计算方法

技术领域

本发明属于电力系统模拟和潮流计算技术领域，具体涉及一种基于PSD-BPA软件的大型基于PSD-BPA的大型电力系统随机潮流计算方法。

背景技术

潮流计算是电力系统规划和运行中不可或缺的重要组成部分，为系统的规划和运行提供基础参考和校验的依据。传统的电力系统潮流计算主要是指根据给定的运行条件及系统接线方式，确定整个电力系统各部分的运行状态。

近年来，伴随着风能、太阳能等新能源大规模并网，电力系统的电源侧出现了大量的随机性和不确定性，使得电力潮流的时空分布也呈现随机波动的特性。因此，需要引入随机潮流，在进行电力系统分析时考虑系统的各种不确定因素的随机性，分析系统运行参数的概率特性，从而使得计算分析更加贴合实际电网的运行状态。

蒙特卡洛模拟是随机潮流计算的常用方法，主要通过概率采样构建场景库，对每一个场景进行确定性分析，并根据多个场景的计算结果得到潮流的近似概率分布。由于蒙特卡洛模拟需要的计算次数较多，常常不得不牺牲单次潮流计算的精度以提高求解速度，尤其在针对大型系统进行分析时难以全面考虑系统中的各种元件和参数。

目前市场上已有PSD-BPA(Power System Department-Bonneville PowerAdministration，由中国电力科学研究院基于引进的BPA软件开发而成)等成熟的潮流分析软件，能够实现潮流迭代计算和变压器分接头调整、无功补偿器投切等复杂调节过程的全面建模，但这些软件通常基于传统的确定性方法，其计算功能难以扩展应用于随机潮流等较复杂的计算分析。因此，设计合理的蒙特卡洛模拟方法，在考虑随机因素波动特性的同时，有效利用传统潮流分析软件的计算能力，实现高效率、高精度的随机潮流计算，对于新形势下的大型电力系统规划与运行具有重大意义。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于PSD-BPA软件的大型基于PSD-BPA的大型电力系统随机潮流计算方法，综合考虑各种随机因素与各种系统构成，同时有效利用传统潮流分析软件的计算能力，实现高效、准确的大型电力系统随机潮流计算。

本发明采用以下技术方案：

一种基于PSD-BPA的大型电力系统随机潮流计算方法，包括以下步骤：

S1、基于历史数据，对风速和光照强度进行概率建模和相关性建模；基于模型采样生成风速和光照强度数据，并转化为出力场景集合，作为步骤S2的输入；

S2、分析系统运行方式数据，根据每个采样场景分别修改对应运行参数，调用PSD-BPA程序完成计算，读取并记录各场景下的潮流计算结果；

S3、读取步骤S2的随机潮流计算结果，绘制经验累积分布函数和概率密度函数；计算节点电压和支路潮流越限的概率，分析系统局部或整体潮流风险，若概率大于允许值，则认为节点/支路存在风险。

具体的，步骤S1具体为：

S101、获取接入系统的多个风电场的历史风速数据、多个光伏电站的历史光照强度数据，并进行异常数据处理；

S102、采用Weibull分布对每个风电场风速的历史数据进行拟合，然后利用极大似然估计对Weibull分布形状参数k和尺度参数λ进行估计；对光照强度数据进行标幺化处理，采用Beta分布对每个光伏电站光照强度的历史数据进行拟合，通过极大似然估计法对Beta分布参数α和β进行估计；

S103、根据各个风电场风速的Weibull累积分布函数及各个光伏电站光照强度的Beta累积分布函数，将出力数据映射为[0，1]区间内的分位数数据；

S104、采用高斯Copula函数对多个风速序列分位数数据和多个光照强度序列分位数数据进行拟合，通过极大似然估计法对高斯Copula函数参数ρ进行估计；

S105、基于建立的模型采样生成风速和光照强度数据，对高斯Copula函数参数ρ进行采样，获得蒙特卡洛样本矩阵W'_m和S'_m，得到满足相关性的[0，1]单位数据；通过各个风电场风速的累积分布函数和各个光伏电站光照强度的累积分布函数，将单位数据转化为风速和光照强度数据；

S106、根据风电场风速与出力的特性曲线，将蒙特卡洛风速数据转化为风电场出力数据；根据光伏电站光照强度与出力的特性曲线，将光照强度数据转化为光伏电站出力数据；将所有出力的场景输出，作为与PSD-BPA交互的输入。

进一步的，步骤S105中，对蒙特卡洛样本矩阵中的每个样本w′_ij或s′_ij，基于步骤S102中拟合得到的累积分布函数F_W(w)和F_B(s)，映射到实际的风速w_ij及光照强度值_ij如下：

w_ij＝F_W(w′_ij)或s_ij＝F_B(s′_ij)

并将映射得到的风速及光照强度值组成新的蒙特卡洛样本矩阵W_m和S_m。

具体的，步骤S2具体为：

S201、读取并分析系统运行方式文件，查找并记录各个相关数据行的位置，包括随机波动的新能源和负荷所在的母线、区域联络线和直流输电线路，存储为内部映射；

S202、从场景集中选取新场景，确定需要修改的系统运行参数，包括新能源机组出力和负荷大小，在内部数据映射中查找参数所在位置，对原有数据进行覆盖或增减，并生成新的运行方式文件；

S203、以PSD-BPA软件潮流计算程序的路径为调用对象，以新生成的此场景运行方式文件路径和指定的输出文件路径为调用参数，生成用于程序调用的Shell命令行语句，在后台通过命令行调用PSD-BPA软件的潮流计算程序，计算对应场景下的电力潮流；

S204、等待调用程序计算完成，读取并解析计算结果文件，检查潮流计算是否收敛：如潮流不收敛，输出警告，并记录场景的信息；如潮流收敛，读取并存储潮流计算结果，包括节点电压的幅值和相位、线路和变压器上流过的有功和无功功率、线路和变压器的有功和无功损耗；

S205、检查是否所有场景均计算完成：如果计算未完成，返回步骤S202；如果计算完成，整理所有场景下的潮流计算结果，包括所有场景下的潮流收敛性和给定运行方式下的系统运行状态，并将数据以文件形式输出。

进一步的，步骤S201具体为：

S2011、以文本形式读取运行方式文件，将全部运行方式数据以行为单位存储为字符串数组，对应数据作为初始运行方式数据的备份，用于在每个场景下修改运行方式时的初始参考值；

S2012、建立运行参数与其对应的数据行位置的内部映射表，初始映射为空；

S2013、选取第一组场景数据，整理出场景中需要修改的运行参数列表，包括随机波动的新能源和负荷所在的母线等，按照运行参数的名称分别查找其所处的具体数据行，记录运行参数到其所在数据行位置的映射，并添加到内部映射表。

进一步的，步骤S202具体为：

S2021、利用原始数据备份重新初始化运行方式数据；

S2022、读取对应采样场景下的风电出力、光伏出力、负荷大小，在内部映射表中查询所在母线对应的数据行，如果查询失败，则在数据中逐行搜索，找到母线对应的数据行位置，并记录到映射表中；如果查询成功，定位对应数据行；

S2023、根据新场景的运行参数采样结果，按照规定的数据格式对所在行的相关数据列进行修改，并保持所在行其他数据列的位置不变；

S2024、将修改后的数据保存为新的运行方式文件，供潮流计算调用。

进一步的，步骤S204具体为：

S2041、等待调用程序运行完成并返回，并读取生成的结果文件；

S2042、在结果文件中搜索潮流收敛信号，判断潮流是否收敛，如果不收敛，记录当前场景信息，并结束场景计算；如果潮流收敛，则进行下一步；

S2043、在潮流结果文件中定位到节点电压结果和支路潮流结果对应的数据块，逐一记录各节点的电压幅值与相位、线路和变压器上流过的有功无功功率以及线路和变压器产生的损耗；其中，对于系统规模在1000节点以上或采样次数在500次以上的情况，只记录部分主要关注的数据结果。

具体的，步骤S3具体为：

S301、读取随机潮流计算结果，统计并绘制各节点电压、各支路潮流的经验累积分布函数，同时通过核密度估计获得各节点电压、各支路功率的经验概率密度函数；

S302、计算系统中各节点电压越限、各支路潮流越限的概率P_c，分析局部以及整体潮流风险，通过求平均值获得区域或整体潮流风险risk如下：

其中，N_c为对应区域或系统整体的节点或支路的个数，P_ci为第i个节点或支路越限概率。

与现有技术相比，本发明至少具有以下有益效果：

本发明一种基于PSD-BPA的大型电力系统随机潮流计算方法，对大型电力系统实现了准确、高效的蒙特卡洛随机潮流计算，为高比例新能源接入下的电力系统规划和运行提供了重要参考；通过设计交互接口，与主流的潮流计算商业软件实现了有效对接，扩展了现有潮流计算工具的应用范围，也为更高级应用工具的开发提供了途径。

进一步的，通过步骤S1能够实现风速和光照强度不确定性概率建模及相关性建模，并能生成大量的出力场景供后续潮流计算。

进一步的，在步骤S105中基于高斯Copula参数矩阵ρ可生成大量满足相关性的单位数据，进而通过累积分布函数转化为风速和光照强度的样本数据集。

进一步的，通过步骤S2能够有效调用PSD-BPA实现大型电力系统潮流计算，进而得到更加准确的蒙特卡洛概率模拟结果。

进一步的，步骤S201在计算前分析运行方式文件并记录相关数据的位置，可以节省后续每次采样计算时查找数据的时间，提高计算效率。

进一步的，通过步骤202修改运行方式数据并生成新的运行方式文件，能够为PSD-BPA提供新场景下的输入数据，从而能够调用该软件完成对应场景的潮流计算。

进一步的，通过步骤204读取并分析调用程序的输出结果，能够获取该场景下的潮流收敛性、节点电压、支路功率等潮流计算结果，同时完成蒙特卡洛计算的迭代控制。

进一步的，通过步骤S3能够得到各节点电压、支路功率的概率分布，从而计算越限概率，判断系统区域或整体潮流风险。

综上所述，本发明开发对PSD-BPA软件实现了有效的调用与数据对接，能够对大型电力系统实现准确、高效的蒙特卡洛随机潮流计算。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明方法步骤S1的分析流程图；

图2为本发明方法步骤S2的分析流程图；

图3为本发明方法步骤S3的分析流程图。

具体实施方式

本发明提供了一种基于PSD-BPA软件的大型基于PSD-BPA的大型电力系统随机潮流计算方法，基于历史数据对风速和光照强度进行概率建模和相关性建模，基于模型进行蒙特卡洛采样生成大量风速和光照强度数据、并转化为时序出力场景；分析系统运行方式数据(.dat文件)，根据每个采样场景分别修改对应运行参数，利用Windows系统接口后台调用PSD-BPA程序完成计算，读取并记录各场景下的潮流计算结果；读取随机潮流计算结果，绘制经验累积分布函数和概率密度函数，计算节点电压和支路潮流越限的概率，分析系统局部或整体潮流风险。本发明对大型电力系统实现了准确、高效的蒙特卡洛随机潮流计算，为高比例新能源接入下的电力系统规划和运行提供了重要参考；同时设计交互接口，与主流的潮流计算商业软件实现了有效对接，扩展了现有潮流计算工具的应用范围，也为更高级应用工具的开发提供了途径。

本发明一种基于PSD-BPA软件的大型电力系统随机潮流计算方法，包括以下步骤：

S1、基于历史数据对风速和光照强度进行概率建模和相关性建模；基于模型采样生成大量风速和光照强度数据，并转化为时序出力场景；

请参阅图1，包括历史数据的获取及处理、风速及光照强度分布拟合、风速及光照强度分位数映射、多个风速序列和光照强度序列的相关性建模、风速及光照强度蒙特卡洛采样、出力数据转换，具体如下：

S101、获取接入系统的多个风电场的历史风速数据、多个光伏电站的历史光照强度数据，并进行异常数据的处理；

获取系统中所有风电场和光伏电站所处位置的历史风速和光照强度序列，组合成风速样本矩阵W_h和光照强度样本矩阵S_h：

其中，W_h和S_h分别为风速和光照强度的样本矩阵，列为每个风电场(光伏电站)的风速(光照强度)序列，共N_w(N_s)个，行为同一时刻所有风电场(光伏电站)的风速(光照强度)值，共有T_w(T_s)个时刻。

通过3σ原则检验异常数据，即对矩阵W_h和S_h的每一列j，计算其均值μ_j和方差σ_j，再依次检验该列所有时刻的数据1≤i≤T_w或T_s，若存在数据满足：

w_ij(s_ij)＜μ_j-3σ_j或w_ij(s_ij)＞μ_j+3σ_j (2)

则判定为该时刻存在异常数据，删除或该行的全部数据。

S102、采用Weibull分布对每个风电场风速的历史数据进行拟合，然后利用极大似然估计对Weibull分布形状参数k和尺度参数λ进行估计；对光照强度数据进行标幺化处理，采用Beta分布对每个光伏电站光照强度的历史数据进行拟合，通过极大似然估计法对Beta分布参数α和β进行估计；

采用Weibull分布对每个风电场风速的历史数据进行拟合具体为：

采用Beta分布对每个光伏电站光照强度的历史数据进行拟合：

S103、根据各个风电场风速的Weibull累积分布函数及各个光伏电站光照强度的Beta累积分布函数，将出力数据映射为[0，1]区间内的分位数数据；

对样本矩阵中的每个样本w_ij或s_ij，基于步骤S102中拟合得到的累积分布函数F_W(w)和F_B(s)(即f_W(w)和f_B(s)的积分函数)，通过反函数法得到风速的分位数映射w′_ij及光照强度的分位数映射s′_ij如下：

或

并将分位数数据组合成新的样本分位数矩阵W'_h和S'_h。

S104、采用高斯Copula函数对多个风速序列分位数数据和多个光照强度序列分位数数据进行拟合，通过极大似然估计法对高斯Copula函数参数ρ进行估计；

对样本分位数矩阵W'_h和S'_h，采用高斯Copula函数进行拟合如下：

其中，V为样本分位数矩阵或，v₁到v_n为矩阵各列代表的分位数随机变量，ρ为高斯Copula相关性参数矩阵，I为单位矩阵，Φ^-1为标准正态分布的逆分布。

S105、基于建立的模型采样生成大量风速和光照强度数据，首先对高斯Copula函数进行采样，获得蒙特卡洛样本矩阵W'_m和S'_m，得到满足相关性的[0，1]单位数据；然后通过各个风电场风速的累积分布函数和各个光伏电站光照强度的累积分布函数，将单位数据转化为风速和光照强度数据；

对蒙特卡洛样本矩阵中的每个样本w′_ij或s′_ij，基于步骤S102中拟合得到的累积分布函数F_W(w)和F_B(s)，映射到实际的风速w_ij及光照强度值_ij如下：

w_ij＝F_W(w′_ij)或s_ij＝F_B(s′_ij) (7)

并将映射得到的风速及光照强度值组成新的蒙特卡洛样本矩阵W_m和S_m。

S106、根据风电场风速与出力的特性曲线，将蒙特卡洛风速数据转化为风电场出力数据；根据光伏电站光照强度与出力的特性曲线，将光照强度数据转化为光伏电站出力数据；将上述所有出力的场景输出，作为与PSD-BPA交互的输入。

蒙特卡洛风速数据转换为风电出力数据P_w(w)如下：

其中，v_ci、v_r、和v_co分别为风力发电机的切入风速，额定风速和切出风速，P_m为风力发电机的额定出力。

蒙特卡洛光照强度数据转换为出力数据P_s(s)如下：

其中，rmin和rmax分别为可出力的最小光照强度和输出功率而定时的光照强度；μ为并网逆变器的输出效率；A为光伏电池阵列面积；η为光伏组件光电转化效率；P_m为光伏阵列的额定出力。

S2、分析系统运行方式数据(.dat文件)，根据每个采样场景分别修改对应运行参数，利用Windows系统接口后台调用PSD-BPA程序完成计算，读取并记录各场景下的潮流计算结果；

请参阅图2，包括运行方式文件的读取和解析、场景运行参数的修改、PSD-BPA潮流计算程序的调用、潮流结果的解析、计算结果的输出，具体为：

S201、读取并分析系统运行方式文件(.dat文件)，查找并记录各个相关数据行的位置，包括随机波动的新能源和负荷所在的母线、区域联络线、直流输电线路等，并存储为内部映射，以便实现快速仿真修改；

S2011、以文本形式读取运行方式文件，将全部运行方式数据以行为单位存储为字符串数组，该数据作为初始运行方式数据的备份，用于在每个场景下修改运行方式时的初始参考值；

S2012、建立运行参数与其对应的数据行位置的内部映射表，初始映射为空；

S2013、选取第一组场景数据，整理出场景中需要修改的运行参数列表，包括随机波动的新能源和负荷所在的母线等，按照运行参数的名称分别查找其所处的具体数据行，记录运行参数到其所在数据行位置的映射，并添加到内部映射表。

S202、从场景集中选取新场景，确定需要修改的系统运行参数，包括新能源机组出力和负荷大小等，在内部数据映射中查找参数所在位置，对原有数据进行覆盖或增减，并生成新的运行方式文件；

S2021、利用原始数据备份重新初始化运行方式数据，消除此前所有修改；

S2022、读取该采样场景下的风电出力、光伏出力、负荷大小等运行参数，在内部映射表中查询其所在母线对应的数据行，如果查询失败，则在数据中逐行搜索，找到该母线对应的数据行位置，并记录到映射表中，避免后续反复查找；如果查询成功，则直接定位该数据行；

S2023、根据新场景的运行参数采样结果，按照规定的数据格式对该行的相关数据列进行修改，并保持该行其他数据列的位置不变；

S2024、将修改后的数据保存为新的运行方式文件(.dat文件)，以供潮流计算调用。

S203、利用Windows系统的Shell命令行接口，以新生成的运行方式文件路径作为参数，在后台通过命令行调用PSD-BPA软件的潮流计算程序，计算该场景下的电力潮流；

S2031、以PSD-BPA软件潮流计算程序的路径为调用对象，以新生成的此场景运行方式文件路径和指定的输出文件路径为调用参数，生成用于程序调用的Shell命令行语句；

S2032、打开Windows系统命令行调用接口(ShellExecuteEx)，以系统控制台(cmd.exe)为执行程序，执行上一步骤中生成的命令行语句。

S204、等待调用程序计算完成，读取并解析计算结果文件(.pfo文件)，检查潮流计算是否收敛：如潮流不收敛，输出警告，并记录下该场景的信息；如潮流收敛，读取并存储潮流计算结果，包括节点电压的幅值和相位、线路和变压器上流过的有功和无功功率、线路和变压器的有功和无功损耗等；

S2041、等待调用程序运行完成并返回，并读取生成的结果文件(.pfo文件)；

S2042、在结果文件中搜索潮流收敛信号，判断潮流是否收敛，如果不收敛，记录下当前场景信息，并结束此场景的计算；如果潮流收敛，则进行下一步；

S2043、在潮流结果文件中定位到节点电压结果和支路潮流结果对应的数据块，逐一记录各节点的电压幅值与相位、线路和变压器上流过的有功无功功率以及线路和变压器产生的损耗；其中，对于规模较大或采样次数较多的系统，可以只记录部分主要关注的数据结果，以减少计算用时与存储空间占用。

S205、检查是否所有场景均计算完成：如否，返回步骤S202；如是，整理所有场景下的潮流计算结果，包括所有场景下的潮流收敛性和给定运行方式下的系统运行状态，并将数据以文件形式输出。

等待所有场景下的潮流均计算完成，将各场景下节点电压和支路潮流结果以csv表格文件形式分别输出，同时输出各场景的收敛性信息。

输出的节点电压结果文件包含收下数据列：采样场景编号、节点类型(交流节点/直流节点)、节点名称、节点电压等级、电压相位、电压幅值、净有功功率注入、净无功功率注入、对地有功功率、对地无功功率、安排有功功率、安排无功功率。

输出的支路潮流结果包含以下数据列：采样场景编号、支路类型(输电线路/变压器)、起始节点名称与电压等级、终止节点名称与电压等级、平行线路编号、正向流过的有功和无功功率、反向流过的有功和无功功率、线路对地充电功率、变压器分接头位置。

S3、读取随机潮流计算结果，绘制经验累积分布函数和概率密度函数；计算节点电压和支路潮流越限的概率，分析系统局部或整体潮流风险。

请参阅图3，包括经验累积分布函数和概率密度函数的绘制、电压潮流越界概率及风险分析，具体为：

S301、读取随机潮流计算结果，统计并绘制各节点电压、各支路潮流的经验累积分布函数，同时通过核密度估计获得各节点电压、各支路功率的经验概率密度函数；

读取PSD-BPA随机潮流计算结果，对每个节点电压或支路潮流的结果数据集，统计得到各分位数α_i,i＝1,…,A_n(A_n为绘制点个数)对应的电压或潮流值y_αi：

P{y＜y_αi}＝α_i (10)

以二维数组{y_αi,α_i}作图即得到节点电压或支路潮流的经验累积分布函数。

此外，对每个节点电压或支路潮流的结果数据集，采用核密度估计获得经验概率密度函数

如下：

其中，x_i为第i个数据，K是核函数，L是数据集长度，h是平滑参数。

S302、计算系统中各节点电压越限、各支路潮流越限的概率，分析局部以及整体潮流风险。

对于某一节点电压或支路潮流Y，计算电压或潮流的越界概率P_c如下：

P_c＝F_Y(y_min)+[1-F_Y(y_max)] (12)

其中，y_min和y_max分别为随机变量Y可允许的最小值和最大值，F_Y为步骤S301中所得经验累积分布函数。

对于某一区域或者系统整体，通过求平均值获得区域或整体潮流风险risk如下：

其中，N_c为该区域或系统整体的节点或支路的个数。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明应用于大规模海上风电接入下的南方电网随机潮流计算。本发明调用PSD-BPA的随机潮流计算与基于MATLAB的牛顿-拉夫逊法蒙特卡洛计算随机潮流相比，两者在计算时间和精度上的对比如表1所示。

表1本发明与常用MATLB蒙特卡洛模拟法对比(10000次)

随机潮流方法	本发明	MATLAB蒙特卡洛
			计算时间	29076.42秒	44669.45秒
计算精度(误差)	参考0％	3.636％

由表1可以看出，采用本发明后，相对于通常基于MATLAB的蒙特卡洛随机潮流，计算时间大幅度减少。这是因为PSD-BPA为成熟的潮流计算商业平台，其程序是经过优化后的，因此在计算时间上具有一定优势。与此同时，两种方法的随机潮流结果非常接近，误差在3.6％左右，因此本发明的计算精度也能够得到保证。

综上所述，本发明一种基于PSD-BPA的大型电力系统随机潮流计算方法，与主流的潮流计算商业软件实现了有效对接，能够针对大型电力系统实现准确、高效的随机潮流计算，为高比例新能源接入下的电力系统规划和运行提供了重要参考；同时扩展了现有潮流计算工具的应用范围，也为更高级应用工具的开发提供了途径。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (4)

1.一种基于PSD-BPA的大型电力系统随机潮流计算方法，其特征在于，包括以下步骤：

S1、基于历史数据，对风速和光照强度进行概率建模和相关性建模；基于模型采样生成风速和光照强度数据，并转化为出力场景集合，作为步骤S2的输入，具体为：

S101、获取接入系统的多个风电场的历史风速数据、多个光伏电站的历史光照强度数据，并进行异常数据处理，组合成风速样本矩阵W_h和光照强度样本矩阵S_h：

其中，W_h和S_h分别为风速和光照强度的样本矩阵，列为每个风电场的风速序列、每个光伏电站的光照强度序列，共N_w、N_s个，行为同一时刻所有风电场的风速值、同一时刻所有光伏电站的光照强度值，共有T_w、T_s个时刻；

通过3σ原则检验异常数据，即对矩阵W_h和S_h的每一列j，计算其均值μ_j和方差σ_j，再依次检验该列所有时刻的数据1≤i≤T_w或T_s，若存在数据满足：

w_ij(s_ij)＜μ_j-3σ_j或w_ij(s_ij)＞μ_j+3σ_j

则判定为该时刻存在异常数据，删除或该行的全部数据；

S102、采用Weibull分布对每个风电场风速的历史数据进行拟合，然后利用极大似然估计对Weibull分布形状参数k和尺度参数λ进行估计；对光照强度数据进行标幺化处理，采用Beta分布对每个光伏电站光照强度的历史数据进行拟合，通过极大似然估计法对Beta分布参数α和β进行估计；采用Weibull分布对每个风电场风速的历史数据进行拟合具体为：

采用Beta分布对每个光伏电站光照强度的历史数据进行拟合：

S103、根据各个风电场风速的Weibull累积分布函数及各个光伏电站光照强度的Beta累积分布函数，将出力数据映射为[0，1]区间内的分位数数据，对样本矩阵中的每个样本w_ij或s_ij，基于步骤S102中拟合得到的累积分布函数F_W(w)和F_B(s)，通过反函数法得到风速的分位数映射w′_ij及光照强度的分位数映射s′_ij如下：

或

并将分位数数据组合成新的样本分位数矩阵W'_h和S'_h；

S104、采用高斯Copula函数对多个风速序列分位数数据和多个光照强度序列分位数数据进行拟合，通过极大似然估计法对高斯Copula函数参数ρ进行估计，对样本分位数矩阵W'_h和S'_h，采用高斯Copula函数进行拟合如下：

其中，V为样本分位数矩阵或，v₁到v_n为矩阵各列代表的分位数随机变量，ρ为高斯Copula相关性参数矩阵，I为单位矩阵，Φ^-1为标准正态分布的逆分布；

S105、基于建立的模型采样生成风速和光照强度数据，对高斯Copula函数参数ρ进行采样，获得蒙特卡洛样本矩阵W'_m和S'_m，得到满足相关性的[0，1]单位数据；通过各个风电场风速的累积分布函数和各个光伏电站光照强度的累积分布函数，将单位数据转化为风速和光照强度数据，对蒙特卡洛样本矩阵中的每个样本w′_ij或s′_ij，基于步骤S102中拟合得到的累积分布函数F_W(w)和F_B(s)，映射到实际的风速w_ij及光照强度值_ij如下：

w_ij＝F_W(w′_ij)或s_ij＝F_B(s′_ij)

并将映射得到的风速及光照强度值组成新的蒙特卡洛样本矩阵W_m和S_m；

S106、根据风电场风速与出力的特性曲线，将蒙特卡洛风速数据转化为风电场出力数据；根据光伏电站光照强度与出力的特性曲线，将光照强度数据转化为光伏电站出力数据；将所有出力的场景输出，作为与PSD-BPA交互的输入，蒙特卡洛风速数据转换为风电出力数据P_w(w)如下：

其中，v_in、v_rated、和v_out分别为风力发电机的切入风速，额定风速和切出风速，P_m为风力发电机的额定出力；

蒙特卡洛光照强度数据转换为出力数据P_s(s)如下：

其中，r_min和r_max分别为可出力的最小光照强度和输出功率而定时的光照强度；μ为并网逆变器的输出效率；A为光伏电池阵列面积；η为光伏组件光电转化效率；P_m为光伏阵列的额定出力；

S2、分析系统运行方式数据，根据每个采样场景分别修改对应运行参数，调用PSD-BPA程序完成计算，读取并记录各场景下的潮流计算结果，具体为：

S201、读取并分析系统运行方式文件，查找并记录各个相关数据行的位置，包括随机波动的新能源和负荷所在的母线、区域联络线和直流输电线路，存储为内部映射；

S202、从场景集中选取新场景，确定需要修改的系统运行参数，包括新能源机组出力和负荷大小，在内部数据映射中查找参数所在位置，对原有数据进行覆盖或增减，并生成新的运行方式文件；

S203、以PSD-BPA软件潮流计算程序的路径为调用对象，以新生成的此场景运行方式文件路径和指定的输出文件路径为调用参数，生成用于程序调用的Shell命令行语句，在后台通过命令行调用PSD-BPA软件的潮流计算程序，计算对应场景下的电力潮流；

S204、等待调用程序计算完成，读取并解析计算结果文件，检查潮流计算是否收敛：如潮流不收敛，输出警告，并记录场景的信息；如潮流收敛，读取并存储潮流计算结果，包括节点电压的幅值和相位、线路和变压器上流过的有功和无功功率、线路和变压器的有功和无功损耗；

S205、检查是否所有场景均计算完成：如果计算未完成，返回步骤S202；如果计算完成，整理所有场景下的潮流计算结果，包括所有场景下的潮流收敛性和给定运行方式下的系统运行状态，并将数据以文件形式输出；

S3、读取步骤S2的随机潮流计算结果，绘制经验累积分布函数和概率密度函数；计算节点电压和支路潮流越限的概率，分析系统局部或整体潮流风险，若概率大于允许值，则认为节点/支路存在风险，具体为：

S301、读取随机潮流计算结果，统计并绘制各节点电压、各支路潮流的经验累积分布函数，同时通过核密度估计获得各节点电压、各支路功率的经验概率密度函数；

读取PSD-BPA随机潮流计算结果，对每个节点电压或支路潮流的结果数据集，统计得到各分位数α_i,i＝1,…,A_n 对应的电压或潮流值y_αi ， A_n 为绘制点个数 ：

P{y＜y_αi}＝α_i

以二维数组{y_αi,α_i}作图即得到节点电压或支路潮流的经验累积分布函数；

对每个节点电压或支路潮流的结果数据集，采用核密度估计获得经验概率密度函数

如下：

其中，x_i为第i个数据，K是核函数，L是数据集长度，h是平滑参数

S302、计算系统中各节点电压越限、各支路潮流越限的概率P_c，分析局部以及整体潮流风险，通过求平均值获得区域或整体潮流风险risk如下：

其中，N_c为对应区域或系统整体的节点或支路的个数，P_ci为第i个节点或支路越限概率；

对于某一节点电压或支路潮流Y，计算电压或潮流的越界概率P_c如下：

P_c＝F_Y(y_min)+[1-F_Y(y_max)]

其中，y_min和y_max分别为随机变量Y可允许的最小值和最大值，F_Y为步骤S301中所得经验累积分布函数；

对于某一区域或者系统整体，通过求平均值获得区域或整体潮流风险risk如下：

其中，N_c为该区域或系统整体的节点或支路的个数。

2.根据权利要求1所述的基于PSD-BPA的大型电力系统随机潮流计算方法，其特征在于，步骤S201具体为：

S2011、以文本形式读取运行方式文件，将全部运行方式数据以行为单位存储为字符串数组，对应数据作为初始运行方式数据的备份，用于在每个场景下修改运行方式时的初始参考值；

S2012、建立运行参数与其对应的数据行位置的内部映射表，初始映射为空；

S2013、选取第一组场景数据，整理出场景中需要修改的运行参数列表，包括随机波动的新能源和负荷所在的母线等，按照运行参数的名称分别查找其所处的具体数据行，记录运行参数到其所在数据行位置的映射，并添加到内部映射表。

3.根据权利要求2所述的基于PSD-BPA的大型电力系统随机潮流计算方法，其特征在于，步骤S202具体为：

S2021、利用原始数据备份重新初始化运行方式数据；

S2022、读取对应采样场景下的风电出力、光伏出力、负荷大小，在内部映射表中查询所在母线对应的数据行，如果查询失败，则在数据中逐行搜索，找到母线对应的数据行位置，并记录到映射表中；如果查询成功，定位对应数据行；

S2023、根据新场景的运行参数采样结果，按照规定的数据格式对所在行的相关数据列进行修改，并保持所在行其他数据列的位置不变；

S2024、将修改后的数据保存为新的运行方式文件，供潮流计算调用。

4.根据权利要求2所述的基于PSD-BPA的大型电力系统随机潮流计算方法，其特征在于，步骤S204具体为：

S2041、等待调用程序运行完成并返回，并读取生成的结果文件；

S2042、在结果文件中搜索潮流收敛信号，判断潮流是否收敛，如果不收敛，记录当前场景信息，并结束场景计算；如果潮流收敛，则进行下一步；

S2043、在潮流结果文件中定位到节点电压结果和支路潮流结果对应的数据块，逐一记录各节点的电压幅值与相位、线路和变压器上流过的有功无功功率以及线路和变压器产生的损耗；其中，对于系统规模在1000节点以上或采样次数在500次以上的情况，只记录部分主要关注的数据结果。

Images