近似法概率潮流

       通常潮流计算仅针对单一的系统运行方式。网络拓扑结构、变压器变比、节点注入功率等都是给定的。而在实际的电力系统运行过程中存在着诸多不确定性因素，例如负荷功率和发电机出力随时都在变化，网络结构由于故障或检修也会改变，以及测量和估计误差等。在电力系统运行分析中，为考查各种不确定性因素对系统运行的影响，需要大量地重复性潮流计算。结合概率理论形成的概率潮流计算能够充分考虑这些不确定性因素。由概率潮流计算可以得到节点电压、电流、支路潮流的均值、方差、概率密度函数等统计特性，从而更全面地反映系统运行状况。

       在传统电力系统分析中，负荷的波动、电网运行方式的变化和发电机的停运等因素造成了电力系统一定程度上的不确定性。随着电力工业的发展，以太阳能和风能等为代表的新能源接入电网，给电网带了明显的间歇性和随机性；微网、分布式电源和电动汽车等配电网新概念的发展，大大增强了电源、负荷与电网之间的互动性，其直接结果导致了电力系统的不确定性显著增加，用于电力系统分析的概率潮流算法的研究日益重要。

        1974年，Borkowaka 提出概率潮流计算方法叫，用以解决电力系统中诸多不确定因素。在随后四十年的时间里，概率潮流理论与方法得到了发展。与其几乎同时出现的随机潮流叫和概率潮流相互补充融合，逐渐形成处理电力系统不确定因素的体系：一般认为对于电力系统短期不确定因素采用随机潮流处理，而对于长期的具备规律性的不确定因素采用概率潮流处理，后者更趋向于概率分布的计算。

        概率潮流计算的提出与发展，其最显著的意义是在进行电力系统分析时，考虑了系统各种不确定因素的随机性，从而使得计算分析更加贴合实际电网的运行状态。概率潮流的研究问题，主要集中在3个层而：系统模型、计算模型和计算方法。

        系统元件的不确定性是引入概率潮流的根本原因，主要体现在发电机、负荷、输电线路和变压器的随机性。近些年，随着可再生能源并网规模的日益提高和电力用户的市场行为日趋突出，发电机和负荷的功率模型越发复杂。

        就概率潮流计算模型而言，以四大类模型为主。Borkowaka基于简化的直流模型团提出了概率潮流计算方法。为了提高潮流计算的精度，Allan分别在1976年和1981年提出了线性化交流模型和分段线性化交流模型，Sokierajski在1978年提出保留非线性的交流模型。目前概率潮流的计算方法都是基于这4种模型进行。

        对概率潮流算法的研究是概率潮流分析中的热点，具备广阔的研究空间与研究意义。一种性能良好的概率潮流计算方法应满足以下指标：能够求出输出随机变量的数字特征(包括均值和方差)及概率分布；能够处理多个随机变量间的相关性；满足实用化要求，结果具有足够精度的情况下尽量减少计算时间；满足通用性要求，对输入变量的数学模型不应有太高要求。这4项指标构成概率潮流算法研究的重点和难点，专家学者们从一方而或多方而入手展开研究，形成了当前的多种概率潮流计算方法。就目前己有研究成果而言，概率潮流算法可以大致分为模拟法、近似法和解析法3类。解析法结合电力系统复杂卷积计算的简化特性，通过对卷积计算的特殊处理衍生而来，广义上也可归入近似法。纵观现有概率潮流算法，存在的突出问题是难以将计算准确度和时效性有机统一。

        以模拟法、近似法和解析法为分类基准，综述各类方法中具体算法的原理步骤、优缺点评估以及发展趋势。以算法中现有的问题及不足为依托，结合当前电力系统的研究热点与新兴技术，对概率潮流的发展方向进行展望。

潮流计算是电力系统分析的基础，其典型计算潮流方程如下所示：

        式中：W为系统节点功率注入向量；X为节点电压向量；Y为系统网络参数；Z为系统支路潮流向量。概率潮流的计算正是在上式的基础上，通过考虑输入变量W和Y的概率特性，获得系统状态变量X和Z的分布情况，从而全而地给出系统的运行状况和概率特征。

        随着电力系统随机性的显著提升，概率潮流算法在电力系统分析中获得广泛应用，其主要应用方向可以分为以下几类。

        1)电力系统规划，包括电源规划、电网规划和无功规划等规划问题。规划问题的求解均以潮流计算作为基础，在不确定性增加的电力系统中，概率潮流将成为含随机因素规划问题的求解前提。

        2)静态安全分析，作为电力系统分析的基本问题，采用概率潮流的静态安全分析方法可以更加真实地反映电力系统全而信息。

        3)电力系统运行状态实时在线分析，包括机组组合、在线调度、电力市场机制下的源一网嗬互动。应用某些概率潮流算法，可以在不显著增加计算次数与时间的条件下，更为精确地分析电力系统的运行状态及变化趋势，为电力系统实时分析提供了强有力的工具。

        4)其他基于潮流运算的系统分析，例如最优化潮流计算、电力系统风险评估、互动型配电网潮流计算等。

模拟法概率潮流

        模拟法概率潮流，是将电力系统中的不确定因素作为随机变量建立概率模型，然后抽取概率分布的样本，最后统计输出变量的分布特征。传统的模拟法概率潮流计算方法一般是指随机采样的蒙特卡洛模拟法，后来基于随机模拟法改进衍生出重要抽样法、拉丁超立方采样法和拟蒙特卡洛方法等。

随机采样的蒙特卡洛模拟法

        蒙特卡洛模拟是二战时期美国物理学家Metropoli在执行曼哈顿计划的过程中提出的。蒙特卡洛模拟法以随机模拟和统计实验为手段，是一种从随机变量的概率分布中，通过随机选择数字的方法产生一种符合该随机变量概率分布特性的随机数值序列，作为输入变量序列进行特定分析的求解方法。其计算关键与核心步骤如下：①对潮流方程的输入变量W构造相应的概率模型；②产生随机数序列，作为系统的抽样输入进行大量的数字模拟，每一组采样值通过潮流计算得到相应的模拟实验值；③系统计算，对模拟实验结果进行统计处理，给出所求问题的解。

        蒙特卡洛模拟的优点在于样本数量足够大时，计算结果足够精确;并且计算量一般不受系统规模的影响，该方法的抽样次数与抽样精度的平方成反比。缺点在于为提高计算精度，往往需要提高系统抽样规模，从而导致计算时长过大。考虑其精度优势，随机采样的蒙特卡洛模拟法一般用来作为基准方法进行比较，是衡量其他方法准确性的重要参考。

重要抽样法

        重要抽样法认为期望值附近的采样值对计算结果具有更大的影响力，因此可以重点关注期望值附近的点。基于此，重要抽样法的基本思路是保持原有样本期望值不变，通过改变己知变量概率分布来减小其方差，从而达到减少运算时间的目的。

        如何选取新分布中系统的概率分布使得随机变量在期望不变的情况下减小方差是重要抽样法的关键步骤。有文献采用迭代法搜索重要分布函数，给出了若干重要分布函数的定义方法，并结合分散抽样的技巧提高重要抽样法的收敛速度。也有文献利用蒙特卡洛方法模拟出负荷样本，然后利用核密度估计方法估计出负荷模型的密度函数，将之作为重要抽样密度函数，计算出支路潮流和节点电压的概率密度函数 。

        重要抽样法在电力系统的概率估计中有着广泛应用，该方法可以快速准确地计算出系统运行状态的期望值，为系统分析提供参考。但重要抽样法中仅以期望为研究对象，对于概率变量的方差、概率分布等参数分析存在天然缺陷，计算结果局限性较大。

拉丁超立方采样法

        为了避免随机采样的蒙特卡洛模拟法的大规模抽样，Mckay等人于1979年提出了拉丁超立方采样法 。它是一种分层采样法，通过改进输入随机变量的样本生成过程，保证其采样值能够有效地反映随机变量的整体分布，算法的出发点就是确保所有的采样区域都能够被采样点覆盖。其基本运算过程分为如下两个步骤:采样和排列。

        拉丁超立方采样法的不足是对输入随机变量的处理较为复杂，一方而要求己知输入随机变量的概率分布函数或累积分布函数，另一方而对不同类型概率分布的随机变量相关性需要特殊变换处理困。但该方法作为一种非常有效的估计输出随机变量期望值的方法，由于采样值能够确保覆盖所有输入随机变量的整个分布区域，无须大规模抽样，并且可以有效处理输入变量之间的相关性和随机性，在准确性、稳健性和时效性上都有较大的优势。

拟蒙特卡洛法

        拟蒙特卡洛法的出发点与拉丁超立方采样法相同，希望通过有效的空间覆盖采样法来规避蒙特卡洛模拟法中的随机抽样。但与拉丁超立方采样法的处理方式不同，拟蒙特卡洛法采用低差异序列实现多维随机变量的空间采样。

        低差异序列，又称伪随机数列，是一系列数值确定的[0,1]区间中的数。在d维变量的空间中，低差异序列中己有n-1个数，生成第n个数的方法是:将这个数插入己有数列中最大的“空白”处，即避免数列在局部空间聚集，从而保证了有限数据的空间全覆盖。

        目前拟蒙特卡洛法己经被应用于概率最优潮流计算和含互动式新能源的电网静态稳定分析中。由于采样过程中一次性生成所需序列，该方法具有比拉丁超立方采样法更高的计算效率。但拟蒙特卡洛法对多变量的高维度问题理论基础薄弱、计算效果差，因此目前多用于小规模电力系统分析计算。

近似法概率潮流

基于半不变量法的概率潮流matlab程序

        近似法是利用输入随机变量的数字特征近似描述系统状态变量统计特性的方法。该方法避开了大规模的重复抽样，因而求解速度较快，又因其能够计及系统输入变量之间的互相关性，因而受到重视。目前研究应用较多的有点估计法、一次二阶矩法和状态变换法。

点估计法

        点估计法是一种概率统计方法，目前所做的应用研究都是基于1998年Hong在己知输入随机变量的连续分布下提出的点估计法。该方法能够根据己知随机变量的概率分布，求得待求随机变量的各阶矩。

        点估计法属于逼近技术的一种，利用输入随机变量的统计信息来逼近输出随机变量的数字特征。其主要运算过程分为以下几步。

        l)用潮流方程中输入随机变量W的各个分布函数求出相应的前2M-1阶中心矩。

        2)通过构造的方式，利用前2M-1阶中心矩独立求出每个输入随机变量的M个离散状态，使得这M个离散状态包含了前2M-1阶中心矩的所有信息。

        3)用所求得的每个输入随机变量的M个离散状态和它们的均值，构造M´K个输入随机变量的离散状态，求出对应输出随机变量的M´K个离散状态。

        4)用求得的潮流方程输出随机变量X和Z的M´K个离散状态逼近相应的期望值与方差等相关数字特征。

        由以上步骤可以分析点估计法的特点如下：

        1)该方法中实际的输入量为输入随机变量前2M-1阶中心矩，此中心矩可以由概率分布函数直接求出，也可以由大量样本逼近拟合方程式展开得到，这样就不必受限于必须己知输入变量概率分布的条件约束。

        2)点估计法不需要知道输入与输出之间的具体函数关系表达式，仅要求每个输入有唯一对应的输出。

        3)输出随机变量有2M-1阶多项式逼近的精度，为了提高估计的精度，可以增加输入变量的高阶矩信息，即增加取点个数。但实际应用中点个数M大于3时不仅急剧增大计算量，而且往往造成解的结果非实数，因此M通常取2或3，即构成常用的两点估计法和三点估计法。

        两点估计法计算简单、容易实现，但其只利用输入变量的前三阶矩信息，计算精度低;三点估计法既能得到较高精度的估计值，又保持了简易性，在点估计法中广为使用。点估计法的缺点在于计算结果中随机变量的高阶矩不够精确，无法准确获得变量的概率分布函数。同时在处理输入变量的时间和空间相关性上具有一定的计算复杂度。

一次二阶矩法

        一次二阶矩法作为一种近似概率仿真方法，己被广泛应用于机械、结构可靠性分析中。该方法通过将状态方程泰勒展开，近似保留一次线性项，形成包含前两阶矩(即均值和方差)的计算方程式。在电力系统概率潮流分析中，其具体步骤如下：

        步骤1：将输入随机变量对输出状态变量的潮流方程按泰勒级数展开为一次项形式。

        步骤2：计算输入变量均值方程式。

        步骤3：计算输入变量协方差方程式。

        步骤4：由步骤2和3联合，通过输入变量的均值和协方差计算输出状态变量的数字特征。

        一次二阶矩法计算简单，效率高;但其计算能力有限，仅能处理输出与输入之间均值和方差的数值计算，算法模型误差较大，并且计算精度受到系统概率潮流模型约束很大，因此研究较少。

状态变换法

        状态变换按照变换方法分为线性变换、多项式变换和无迹变换等。线性变换法基于正态变量线性变换不变性定理，假设节点注入随机变量均服从正态分布，将潮流方程线性化后可得系统状态变量为节点注入变量的线性组合并且仍服从正态分布。多项式变换多用作其他计算方法的辅助手段，用以表征电力系统随机因素的模型转换等问题。无迹变换认为:拟合一个概率分布比求解非线性变换容易得多，基于此，通过较少的样本点和相应的样本权重准确捕获状态分布参数，通过非线性函数传递后输出状态变量的期望与方差。

        状态变换法的优点在于变换过程数学理论清晰，意义明确，计算规模不大，但由于该方法以高斯正态分布为变换基础，使其在新能源(不满足高斯分布)并网问题下的概率分析存在一些不足。

近似法概率潮流展望

        近似在各个科学领域均有应用，在电力系统概率潮流中近似计算的使用也较为普遍。与确定性潮流相比，概率潮流的计算规模显著提升，快速的近似计算显得更为迫切。目前己有的近似法概率潮流算法中，计算结果以均值和方差为主要目标，如何获得状态变量较为准确的整体概率分布是改进的重要方向。另一方而，近似法概率潮流算法计算结果的可信度与误差分析也是研究的内容之一。

        在近似法中，随机变量状态的变换是一种基础计算手段，而泛函分析领域的空间转换与之具有相似的计算逻辑，把发展成熟的泛函空间变换应用于近似法概率潮流，是值得探索的方向。

解析法概率潮流

        解析法概率潮流主要关注如何利用随机变量间的关系进行卷积计算得到状态量的概率分布，即核心思想在于有效处理复杂的卷积计算。在卷积计算的过程中，往往会做一些近似处理，因此解析法广义上可归为近似法。传统的卷积技术是一种获得概率潮流的基本方法，假设变量独立，在己知注入概率分布的情况下，可以通过卷积技术得到待求状态变量的概率密度函数。但对于多元线性方程的卷积计算量十分庞大，因此限制了其使用。己有解析法中卷积计算多采用快速傅里叶变换、半不变量法〕和序列运算理论。

        解析法处理卷积计算的应用前提是假设输入随机变量相互独立，导致解析法在处理变量相关性上具有固有缺陷，因此，解析方法中变量相关性处理是研究的热点。

快速傅里叶变换

        在信号处理学科领域，快速傅里叶变换是处理卷积问题的最佳方式，其具备良好的精度和效率，对于处理小规模数据系统具备很大的优越性。在电力系统方而，随着系统规模的增大，系统分析的输入变量急剧增加，使得快速傅里叶变换不再具备精度和效率上的优势，因此，在20世纪80年代经历了短暂的研究后，快速傅里叶变换逐渐退出电力系统概率潮流分析领域。

半不变量法

        目前在电力系统广泛使用的处理卷积的方法是半不变量法。该方法的核心思想是将复杂的卷积运算转换为半不变量之间简单的算术运算，从而大大降低计算过程的复杂度。具体计算步骤可以概括如下:将潮流方程中的随机变量w进行概率分布拟合，经过中心矩计算出输入随机变量的各阶半不变量；然后，通过雅可比矩阵和灵敏度矩阵进行简单的数学运算获得潮流方程输出变量X和Z的各阶半不变量，结合不同级数扩展方式得到相应输出变量的概率密度扩展方程，从而获得输出随机变量的概率分布。有文献深入分析了半不变量法计算随机潮流时各环节的假设条件及可能引起的误差，并提出如何处理节点功率相关性、故障列表和调度策略等问题，使得半不变量法更加实用化。

        半不变量法具备的最大优势在于计算方法简单、计算效率高，虽然计算结果精确度存在一定争议，但满足工程应用要求，因此受到广泛研究。关于半不变量法的算法改进研究主要集中在3个问题上:一是如何处理输入随机变量之间的相关性；二是如何正确分析静态安全稳定问题;三是如何更为精确地描述系统运行状态。有文献提出一种基于Cholesky分解的计及输入变量相关性的半不变量法概率潮流计算方法，并提出基于蒙特卡洛抽样的方法解决一些输入变量的半不变量难以被常规数值方法求解的问题。

部分代码如下：

基于半不变量法的概率潮流matlab程序

clc
clear
%===================考虑分布式电源、发电机和负荷随机波动的概率潮流计算================================
tic                        %半不变量法计算计时开始
%% 基础参数------------------------------------------------------------------
[Nodes,linenum,SB,maxIters,OPdata1,precision,OPdata2,balanceID,balancenotes,...
  lineID,linei,linej,liner,linex,lineb,...
  branchi,branchb,...
  transID,transi,transj,transr,transx,transk,transkMin,transkMax,...
  PQi,PG,QG,PD,QD,...
  PVi,PVV,PVQmin,PVQmax...
  NGi,OP_0,OP_1,OP_2,NGmin,NGmax]=dataIn('IEEE34.txt');  %% 将数据放入各变量后以列向量的格式输出
%% 首先进行基础潮流计算，形成雅克比矩阵
%形成交流系统节点导纳矩阵----------------------------------------------------
[Y,Y0] = formACY(Nodes,branchi,branchb,linei,linej,liner,...
               linex,lineb,transi,transj,transr,transx,transk);
%潮流计算-------------------------------------------------------------------
[V,deta,PQ_loss,S,detaS,Colab,Jacco,Jacco2 ]...
      = NR_main(PVi,PVV,balancenotes,Y,Y0,linei,linej,transi,transj,...
                PG,PD,QG,QD,maxIters,precision,Nodes);      
%%   计算输入的半不变量         
%发电机的随机参数输入--------------------------------------------------------
%%pdfgen(i,1)为发电机序号
%%pdfgen(i,2)为发电机的节点号
%%pdfgen(i,3)为发电机的有功出力
%%pdfgen(i,4)为发电机的无功出力
%%pdfgen(i,5)为发电机的出力的概率
%发电机的八阶半不变量形成-------------------------------------------------------
pdfgen=textread('IEEE34gen.txt');%%普通发电机出力服从二项分布;
ngen=length(pdfgen(:,1));
PgPx=zeros(Nodes,8);
PgQx=zeros(Nodes,8);
PgPx(pdfgen(:,2),:)=NcalGCum(pdfgen(:,3),pdfgen(:,5));
PgQx(pdfgen(:,2),:)=NcalGCum(pdfgen(:,4),pdfgen(:,5));
%负荷的八阶半不变量形成-------------------------------------------------------
%%pdfload(i,1)为负荷序号
%%pdfload(i,2)为负荷的节点号
%%pdfload(i,3)为负荷有功均值
%%pdfload(i,4)为负荷无功均值
%%pdfload(i,5)为负荷有功标准差
%%pdfload(i,6)为负荷无功标准差     %%标准差给定可以参照“3Sita原则” 
%负荷的八阶半不变量---------------------------------------------------------
pdfload=textread('IEEE34load_30%.txt');%%负荷负荷正态分布
nload=length(pdfload(:,1));
PlPx=zeros(Nodes,8);
PlQx=zeros(Nodes,8);
PlPx(pdfload(:,2),:)=NcalPLCum(-pdfload(:,3),-pdfload(:,5));
PlQx(pdfload(:,2),:)=NcalPLCum(-pdfload(:,4),-pdfload(:,6));

%-------光伏随机特性建模-----------------------------
%选择上海31°8’N、121°35’E作为光照强度分布的考量位置，在HOMERE软件上获取光强分布的期望值和方差。
%miu=0.150314263;
%sita=0.049758487;
%利用HOMER软件获取广州（113°15′E，23°7′N）的光照强度数据样本作为后续应用的模型
M=textread('Guangzhao.txt');
maxM=max(M(:,1))
miuM=mean(M(:,1)./maxM)
sitaM=sqrt(var(M(:,1)./maxM))
a=miuM^2*(1-miuM)/sitaM^2-miuM
b=miuM*(1-miuM)^2/sitaM^2-(1-miuM)
nbins=30;                                        %设置直方图需要绘制的区间数
Max=max(M(:,1));
Min=min(M(:,1));
Length=(Max-Min)/nbins;
[nhist,c_points]=hist(M(:,1),nbins);      %获取所绘制的频数直方图各区间的中心以及各区间对应的元素的个数
figure
x=0:0.001:1;
y1=betapdf(x,a,b);
y2=betacdf(x,a,b);
subplot(1,2,1);
bar(c_points,nhist/size(M,1)/Length,'grouped','white');
hold on;
plot(x,y1,'LineWidth',2);
hold on;
title('Beta分布的概率密度函数','FontSize',16);
grid on;
box on;
subplot(1,2,2);
plot(x,y2,'LineWidth',2);
title('Beta分布的累积分布函数','FontSize',16);
grid on;
box on;
%在节点33上加入2个额定功率为100kW的太阳能电池方阵
PdsPx=zeros(Nodes,8);
%计算太阳能光伏随机输出的八阶原点矩
a1=a/(a+b);
a2=a*(a+1)/((a+b)*(a+b+1));
a3=a*(a+1)*(a+2)/((a+b)*(a+b+1)*(a+b+2));
a4=a*(a+1)*(a+2)*(a+3)/((a+b)*(a+b+1)*(a+b+2)*(a+b+3));
a5=a*(a+1)*(a+2)*(a+3)*(a+4)/((a+b)*(a+b+1)*(a+b+2)*(a+b+3)*(a+b+4));
a6=a*(a+1)*(a+2)*(a+3)*(a+4)*(a+5)/((a+b)*(a+b+1)*(a+b+2)*(a+b+3)*(a+b+4)*(a+b+5));
a7=a*(a+1)*(a+2)*(a+3)*(a+4)*(a+5)*(a+6)/((a+b)*(a+b+1)*(a+b+2)*(a+b+3)*(a+b+4)*(a+b+5)*(a+b+6));
a8=a*(a+1)*(a+2)*(a+3)*(a+4)*(a+5)*(a+6)*(a+7)/((a+b)*(a+b+1)*(a+b+2)*(a+b+3)*(a+b+4)*(a+b+5)*(a+b+6)*(a+b+7));   
%计算太阳能光伏随机输出的八阶半不变量
PdsPx(34,1)=a1;
PdsPx(34,2)=a2-a1^2;
PdsPx(34,3)=a3-3*a1*a2+2*a1^3;
PdsPx(34,4)=a4-3*a2^2-4*a1*a3+12*a1^2*a2-6*a1^4;
PdsPx(34,5)=a5-5*a1*a4-10*a2*a3+20*a1^2*a3+30*a2^2*a1-60*a1^3*a2+24*a1^5;
PdsPx(34,6)=a6-6*a1*a5-15*a2*a4+30*a1^2*a4-10*a3^2+120*a1*a2*a3-120*a1^3*a3+30*a2^3-270*a2^2*a1^2+360*a1^4*a2-120*a1^6;
PdsPx(34,7)=a7-7*a1*a6-21*a2*a5+42*a1^2*a5-35*a3*a4+210*a1*a2*a4-210*a1^3*a4+140*a3^2*a1+210*a2^2*a3-1260*a3*a2*a1^2+840*a3*a1^4-630*a1*a2^3+2520*a2^2*a1^3-2520*a1^5*a2+720*a1^7;
PdsPx(34,8)=a8-7*a1*PdsPx(34,7)-21*a2*PdsPx(34,6)-35*a3*PdsPx(34,5)-35*a4*PdsPx(34,4)-21*a5*PdsPx(34,3)-7*a6*PdsPx(34,2)-a7*PdsPx(34,1);

rmax=maxM*10^-3;                         %最大辐射度
A=800;                                         %光伏组件总面积  400 为100kw
prey=0.13;                                    %光电转换效率
m=A*rmax*prey;
PdsPx(34,1)=m*PdsPx(34,1);
PdsPx(34,2)=(m^2)*PdsPx(34,2);
PdsPx(34,3)=(m^3)*PdsPx(34,3);
PdsPx(34,4)=(m^4)*PdsPx(34,4);
PdsPx(34,5)=(m^5)*PdsPx(34,5);
PdsPx(34,6)=(m^6)*PdsPx(34,6);
PdsPx(34,7)=(m^7)*PdsPx(34,7);
PdsPx(34,8)=(m^8)*PdsPx(34,8);
PlPx(:,1)=0;
PdsPx(:,1)=0;
PlQx(:,1)=0;
%加入光伏后的八阶半不变量矩阵修改
B_s=[PlPx+PdsPx;PlQx];                %矩阵B相当于deltaW的半不变量
S1=Jacco\eye(size(Jacco));          %函数eye创建单位矩阵，并求雅可比矩阵的逆
S2=S1.^2;                    
S3=S1.^3;
S4=S1.^4;
S5=S1.^5;
S6=S1.^6;
S7=S1.^7;
S8=S1.^8;
%电压幅值及相角的八阶半不变量------------------------------------------------
k_s(:,1)=S1*B_s(:,1);
k_s(:,2)=S2*B_s(:,2);
k_s(:,3)=S3*B_s(:,3);
k_s(:,4)=S4*B_s(:,4);
k_s(:,5)=S5*B_s(:,5);
k_s(:,6)=S6*B_s(:,6);
k_s(:,7)=S7*B_s(:,7);
k_s(:,8)=S8*B_s(:,8);
%随机变量电压幅值(DetaV/V)的八阶半不变量----------------------------------------------------------
Vk_s=k_s(Nodes+1:2*Nodes,:);                 %% 前34个量是电压相角的八阶半不变量
%DetaV的八阶半不变量
V1=V;
V2=V.^2;                    
V3=V.^3;
V4=V.^4;
V5=V.^5;
V6=V.^6;
V7=V.^7;
V8=V.^8;
DetaVk_s(:,1)=V1.*Vk_s(:,1);
DetaVk_s(:,2)=V2.*Vk_s(:,2);
DetaVk_s(:,3)=V3.*Vk_s(:,3);
DetaVk_s(:,4)=V4.*Vk_s(:,4);
DetaVk_s(:,5)=V5.*Vk_s(:,5);
DetaVk_s(:,6)=V6.*Vk_s(:,6);
DetaVk_s(:,7)=V7.*Vk_s(:,7);
DetaVk_s(:,8)=V8.*Vk_s(:,8);
%PQ节点的位置---------------------------------------------------------------
Vpqi=[2:1:34];
%随机变量电压(PQ节点)的八阶半不变量--------------------------------------------------
DetaVpqik_s=DetaVk_s(Vpqi',:);                   %% 取出33个PQ节点的电压幅值的八阶半不变量
%支路潮流的灵敏度矩阵--------------------------------------------------------
G0=formG0(Jacco2,lineID,linei,linej,transi,transj,S,Nodes); 
T1=G0*S1;
T2=T1.^2;
T3=T1.^3;
T4=T1.^4;
T5=T1.^5;
T6=T1.^6;
T7=T1.^7;
T8=T1.^8;
%支路潮流的八阶半不变量------------------------------------------------------
kb_s(:,1)=T1*B_s(:,1);
kb_s(:,2)=T2*B_s(:,2);
kb_s(:,3)=T3*B_s(:,3);
kb_s(:,4)=T4*B_s(:,4);
kb_s(:,5)=T5*B_s(:,5);
kb_s(:,6)=T6*B_s(:,6);
kb_s(:,7)=T7*B_s(:,7);
kb_s(:,8)=T8*B_s(:,8);
%支路有功功率的八阶半不变量--------------------------------------------------
b=length(lineID);
kp_s=kb_s(1:b,:);
%支路无功功率的八阶半不变量--------------------------------------------------
kq_s=kb_s(b+1:2*b,:);
%% 用Gram-charlie展开
for i=1:8
    gv_s(:,i)=DetaVpqik_s(:,i)./DetaVpqik_s(:,2).^(i/2);
    gp_s(:,i)=kp_s(:,i)./kp_s(:,2).^(i/2);
    gq_s(:,i)=kq_s(:,i)./kq_s(:,2).^(i/2);
 end
%电压Gram-charlie展开系数---------------------------------------------------
cv_s=zeros(size(gv_s));
cv_s(:,3)=-1*gv_s(:,3);
cv_s(:,4)=gv_s(:,4);
cv_s(:,5)=(-1)*gv_s(:,5);
cv_s(:,6)=gv_s(:,6)+10*gv_s(:,3).^2;
cv_s(:,7)=(-1)*(gv_s(:,7)+35*gv_s(:,3).*gv_s(:,4));
cv_s(:,8)=gv_s(:,8)+56*gv_s(:,3).*gv_s(:,5)+35*gv_s(:,4).^2;

%线路有功Gram-charlie展开系数-----------------------------------------------
cp_s=zeros(size(gp_s));
cp_s(:,3)=-1*gp_s(:,3);
cp_s(:,4)=gp_s(:,4);
cp_s(:,5)=(-1)*gp_s(:,5);
cp_s(:,6)=gp_s(:,6)+10*gp_s(:,3).^2;
cp_s(:,7)=(-1)*(gp_s(:,7)+35*gp_s(:,3).*gp_s(:,4));
cp_s(:,8)=gp_s(:,8)+56*gp_s(:,3).*gp_s(:,5)+35*gp_s(:,4).^2;

%线路无功Gram-charlie展开系数-----------------------------------------------
cq_s=zeros(size(gq_s));
cq_s(:,3)=-1*gq_s(:,3);
cq_s(:,4)=gq_s(:,4);
cq_s(:,5)=(-1)*gq_s(:,5);
cq_s(:,6)=gq_s(:,6)+10*gq_s(:,3).^2;
cq_s(:,7)=(-1)*(gq_s(:,7)+35*gq_s(:,3).*gq_s(:,4));
cq_s(:,8)=gq_s(:,8)+56*gq_s(:,3).*gq_s(:,5)+35*gq_s(:,4).^2;

%电压幅值均方差-------------------------------------------------------------
Vsgm=sqrt(DetaVpqik_s(:,2));      %% 求的是每个PQ节点的电压幅值标准差
%由基准值计算出来的电压期望值------------------------------------------------
Vmui=V(Vpqi);       %% 求的是每个PQ节点的电压幅值期望值，在NR_main.m潮流计算中求出各节点电压幅值初始值
%线路有功\无功均方差--------------------------------------------------------
li = [linei;transi];
lj = [linej;transj];
S_lin=S(sub2ind([Nodes,Nodes],li,lj));     %% 按索引号去获取对应线路的有功无功
psgm=sqrt(kp_s(:,2)); 
qsgm=sqrt(kq_s(:,2));
%由基准值计算出来的线路有功、无功期望值---------------------------------------
pmui=real(S_lin);
%pmui=pmui(1:b,:)+kp(:,1);
qmui=imag(S_lin);
%qmui=qmui(1:b,:)+kq(:,1);
for y=1:length(Vpqi)                     %% 求解的是每个PQ节点的电压幅值的分布情况
    x=0.80:0.0001:1.20;                 
    x1=(x-Vmui(y))/Vsgm(y);             %% 对每个PQ节点的电压随机变量做标准化
    fai_V=normpdf(x1,0,1);       %% 标准化正态分布的概率密度函数
    f1_s(y,:)=fai_V+cv_s(y,3).*(3*x1-x1.^3).*fai_V/factorial(3)+cv_s(y,4).*(x1.^4-6*x1.^2+3).*fai_V/factorial(4)+...
          cv_s(y,5).*(-x1.^5+10*x1.^3-15*x1).*fai_V/factorial(5)+cv_s(y,6).*(x1.^6-15*x1.^4+45*x1.^2-15).*fai_V/factorial(6)+...
          cv_s(y,7).*(-x1.^7+21*x1.^5-105*x1.^3-105*x1).*fai_V/factorial(7)+cv_s(y,8).*(x1.^8-28*x1.^6+210*x1.^4-420*x1.^2+105).*fai_V/factorial(8);
    f3_s(y,:)=f1_s(y,:)/Vsgm(y);
    f2_s(y,:)=normcdf(x1,0,1)+cv_s(y,3).*(1-x1.^2).*fai_V/factorial(3)+cv_s(y,4).*(3*x1-x1.^3).*fai_V/factorial(4)+...
            cv_s(y,5).*(x1.^4-6*x1.^2+3).*fai_V/factorial(5)+cv_s(y,6).*(-x1.^5+10*x1.^3-15*x1).*fai_V/factorial(6)+...
            cv_s(y,7).*(x1.^6-15*x1.^4+45*x1.^2-15).*fai_V/factorial(7)+cv_s(y,8).*(-x1.^7+21*x1.^5-105*x1.^3+105*x1).*fai_V/factorial(8);
end

for m=1:b
    Pi=-1.0:0.001:1.5;
    Px1=(Pi-pmui(m))/psgm(m); %% 对有功功率随机变量做标准化
    fai_P=1./sqrt(2*pi).*exp((-0.5)*(Px1.^2));
    fp1_s(m,:)=fai_P+cp_s(m,3).*(3*Px1-Px1.^3).*fai_P/factorial(3)+cp_s(m,4).*(Px1.^4-6*Px1.^2+3).*fai_P/factorial(4)+...
          cp_s(m,5).*(-Px1.^5+10*Px1.^3-15*Px1).*fai_P/factorial(5)+cp_s(m,6).*(Px1.^6-15*Px1.^4+45*Px1.^2-15).*fai_P/factorial(6)+...
          cp_s(m,7).*(-Px1.^7+21*Px1.^5-105*Px1.^3+105*Px1).*fai_P/factorial(7)+cp_s(m,8).*(Px1.^8-28*Px1.^6+210*Px1.^4-420*Px1.^2+105).*fai_P/factorial(8);
    fp3_s(m,:)=fp1_s(m,:)/psgm(m);
    fp2_s(m,:)=normcdf(Px1,0,1)+cp_s(m,3).*(1-Px1.^2).*fai_P/factorial(3)+cp_s(m,4).*(3*Px1-Px1.^3).*fai_P/factorial(4)+...
            cp_s(m,5).*(Px1.^4-6*Px1.^2+3).*fai_P/factorial(5)+cp_s(m,6).*(-Px1.^5+10*Px1.^3-15*Px1).*fai_P/factorial(6)+...
            cp_s(m,7).*(Px1.^6-15*Px1.^4+45*Px1.^2-15).*fai_P/factorial(7)+cp_s(m,8).*(-Px1.^7+21*Px1.^5-105*Px1.^3+105*Px1).*fai_P/factorial(8);
end

for m=1:b
    Qi=-0.5:0.001:0.5;
    Qx1=(Qi-qmui(m))/qsgm(m); %% 对无功功率随机变量做标准化
    fai_Q=1./sqrt(2*pi).*exp((-0.5)*(Qx1.^2));
    fq1_s(m,:)=fai_Q+cq_s(m,3).*(3*Qx1-Qx1.^3).*fai_Q/factorial(3)+cq_s(m,4).*(Qx1.^4-6*Qx1.^2+3).*fai_Q/factorial(4)+...
          cq_s(m,5).*(-Qx1.^5+10*Qx1.^3-15*Qx1).*fai_Q/factorial(5)+cq_s(m,6).*(Qx1.^6-15*Qx1.^4+45*Qx1.^2-15).*fai_Q/factorial(6)+...
          cq_s(m,7).*(-Qx1.^7+21*Qx1.^5-105*Qx1.^3+105*Qx1).*fai_Q/factorial(7)+cq_s(m,8).*(Qx1.^8-28*Qx1.^6+210*Qx1.^4-420*Qx1.^2+105).*fai_Q/factorial(8);
    fq3_s(m,:)=fq1_s(m,:)/qsgm(m);
    fq2_s(m,:)=normcdf(Qx1,0,1)+cq_s(m,3).*(1-Qx1.^2).*fai_Q/factorial(3)+cq_s(m,4).*(3*Qx1-Qx1.^3).*fai_Q/factorial(4)+...
            cq_s(m,5).*(Qx1.^4-6*Qx1.^2+3).*fai_Q/factorial(5)+cq_s(m,6).*(-Qx1.^5+10*Qx1.^3-15*Qx1).*fai_Q/factorial(6)+...
            cq_s(m,7).*(Qx1.^6-15*Qx1.^4+45*Qx1.^2-15).*fai_Q/factorial(7)+cq_s(m,8).*(-Qx1.^7+21*Qx1.^5-105*Qx1.^3+105*Qx1).*fai_Q/factorial(8);
end

toc                                               %半不变量法计算计时结束
t1=toc

%求解各节点电压越限概率，依据欧洲标准EN 50160的规定
P_Vofflimits_CMGC=ProbCMGC(f2_s,length(Vpqi),x)

%====================================================================================================================================================================================================
%% Cornish_Fisher级数展开
%由于级数展开式中需用到规格化后的各阶半不变量数据，而gv_s、gp_s、gq_s已完成除了第一阶的标准化
%电压分布情况
alpha=0:0.00001:1;
z=norminv(alpha,0,1);                %标准正态分布的上alpha分位数
%z=-4:0.001:4;
%alpha=normcdf(z,0,1);
for i=1:size(gv_s,1)
    xbV(i,:)=z+(z.^2-1)*gv_s(i,3)/6+(z.^3-3*z)*gv_s(i,4)/24-...
    (2*z.^3-5*z)*(gv_s(i,3)^2)/36+(z.^4-6*z.^2+3)*gv_s(i,5)/120-...
    (z.^4-5*z.^2+2)*gv_s(i,3)*gv_s(i,4)/24+(12*z.^4-53*z.^2+17)*(gv_s(i,3)^3)/324;
    ybV(i,:)=xbV(i,:).*sqrt(DetaVpqik_s(i,2))+Vmui(i,1);
end
%有功分布情况
alpha=0:0.00001:1;
z=norminv(alpha,0,1);                %标准正态分布的上alpha分位数
%z=-4:0.001:4;
%alpha=normcdf(z,0,1);
for i=1:size(gp_s,1)
    xbP(i,:)=z+(z.^2-1)*gp_s(i,3)/6+(z.^3-3*z)*gp_s(i,4)/24-...
    (2*z.^3-5*z)*(gp_s(i,3)^2)/36+(z.^4-6*z.^2+3)*gp_s(i,5)/120-...
    (z.^4-5*z.^2)*gp_s(i,2)*gp_s(i,3)/24+(12*z.^4-53*z.^2)*(gp_s(i,3)^3)/324;
    ybP(i,:)=xbP(i,:).*sqrt(kp_s(i,2))+pmui(i,1);
end
%无功分布情况
alpha=0:0.00001:1;
z=norminv(alpha,0,1);                %标准正态分布的上alpha分位数
%z=-4:0.001:4;
%alpha=normcdf(z,0,1);
for i=1:size(gq_s,1)
    xbQ(i,:)=z+(z.^2-1)*gq_s(i,3)/6+(z.^3-3*z)*gq_s(i,4)/24-...
    (2*z.^3-5*z)*(gq_s(i,3)^2)/36+(z.^4-6*z.^2+3)*gq_s(i,5)/120-...
    (z.^4-5*z.^2)*gq_s(i,2)*gq_s(i,3)/24+(12*z.^4-53*z.^2)*(gq_s(i,3)^3)/324;
    ybQ(i,:)=xbQ(i,:).*sqrt(kq_s(i,2))+qmui(i,1);
end
%====================================================================================================================================================================================================
%% MonteCarlo模拟----
daishu=6000;%抽样次数
%tic                                      %MonteCarlo模拟法计算计时开始
%基础参数------------------------------------------------------------------
[Nodes,linenum,SB,maxIters,OPdata1,precision,OPdata2,balanceID,balancenotes,...
  lineID,linei,linej,liner,linex,lineb,...
  branchi,branchb,...
  transID,transi,transj,transr,transx,transk,transkMin,transkMax,...
  PQi,PG,QG,PD,QD,...
  PVi,PVV,PVQmin,PVQmax...
  NGi,OP_0,OP_1,OP_2,NGmin,NGmax]=dataIn('IEEE34.txt');  %% 将数据放入各变量后以列向量的格式输出
pdfload=textread('IEEE34load_30%.txt');     %%负荷负荷正态分布
%利用HOMER软件获取广州（113°15′E，23°7′N）的光照强度数据样本作为后续应用的模型
M=textread('Guangzhao.txt');
maxM=max(M(:,1))
miuM=mean(M(:,1)./maxM)
sitaM=sqrt(var(M(:,1)./maxM))
a=miuM^2*(1-miuM)/sitaM^2-miuM
b=miuM*(1-miuM)^2/sitaM^2-(1-miuM)
Vfenbu_s=zeros(Nodes,daishu);
for iii=1:daishu
  [nPQi,mPQi]=size(PQi);
  for i=1:nPQi
    PD(PQi(i))=normrnd(pdfload(i,3),pdfload(i,5));       %%函数normrand生成正态分布随机数      
    QD(PQi(i))=normrnd(pdfload(i,4),pdfload(i,6));
  end
  nPQi=20;
  rmax=maxM*10^-3;                         %最大辐射度
  A=800;                                     %光伏组件400，每个2m^2
  prey=0.13;                               %光电转换效率13%
  r=rmax*random('beta',a,b);
  PG(PQi(20))=A*r*prey;
%形成交流系统节点导纳矩阵----------------------------------------------------
[Y,Y0] = formACY(Nodes,branchi,branchb,linei,linej,liner,...
               linex,lineb,transi,transj,transr,transx,transk);
%潮流计算-------------------------------------------------------------------
[V,deta,PQ_loss,S,detaS,Colab,Jacco,Jacco2 ]...
      = NR_main(PVi,PVV,balancenotes,Y,Y0,linei,linej,transi,transj,...
                PG,PD,QG,QD,maxIters,precision,Nodes);
li = [linei;transi];
lj = [linej;transj];
S_linefenbu_s(:,iii)=S(ind2sub(size(S),sub2ind(size(S),li,lj)));
P_linefenbu_s(:,iii)=real(S_linefenbu_s(:,iii));
Q_linefenbu_s(:,iii)=imag(S_linefenbu_s(:,iii));
Vfenbu_s(:,iii)=V;
end

%MonteCarlo计算得累积概率分布情况
%节点电压累积概率
nrnd=100;                        %区间生成随机点个数
VA_s=tabulate(Vfenbu_s(34,:));        %给出随机数生成后的信息矩阵，包括value、count、percent
VB_s=linspace(0.85,1.15,nrnd);      %在区间[0.93,1.07]等间距生成nrnd个点
VC_s=zeros(1,nrnd);                   %生成一个空矩阵，将个点概率放入
for i=1:size(VB_s,2)
    for j=1:size(VA_s,1)
        if VB_s(i)>=VA_s(j,1)
            VC_s(1,i)=VC_s(1,i)+VA_s(j,3);
        end
    end    
end
VC_s=VC_s./100;
%支路有功累积概率
pA_s=tabulate(P_linefenbu_s(33,:));        %给出随机数生成后的信息矩阵，包括value、count、percent
pB_s=linspace(-0.5,0.5,nrnd);      %在区间[0.93,1.07]等间距生成nrnd个点
pC_s=zeros(1,nrnd);                   %生成一个空矩阵，将个点概率放入
for i=1:size(pB_s,2)
    for j=1:size(pA_s,1)
        if pB_s(i)>=pA_s(j,1)
            pC_s(1,i)=pC_s(1,i)+pA_s(j,3);
        end
    end    
end
pC_s=pC_s./100;
%支路无功累积概率
qA_s=tabulate(Q_linefenbu_s(33,:));        %给出随机数生成后的信息矩阵，包括value、count、percent
qB_s=linspace(-0.5,0.5,nrnd);      %在区间[0.93,1.07]等间距生成nrnd个点
qC_s=zeros(1,nrnd);                   %生成一个空矩阵，将个点概率放入
for i=1:size(qB_s,2)
    for j=1:size(qA_s,1)
        if qB_s(i)>=qA_s(j,1)
            qC_s(1,i)=qC_s(1,i)+qA_s(j,3);
        end
    end    
end
qC_s=qC_s./100;

%样条插值
VBi_s=linspace(0.85,1.15,nrnd);
VCi_s=spline(VB_s,VC_s,VBi_s);
pBi_s=linspace(-0.5,0.5,nrnd);
pCi_s=spline(pB_s,pC_s,pBi_s);
qBi_s=linspace(-0.5,0.5,nrnd);
qCi_s=spline(qB_s,qC_s,qBi_s);
toc                                                      %MonteCarlo模拟法计算计时结束
t2=toc

nbins=50;                                        %设置直方图需要绘制的区间数
Max=max(Vfenbu_s(34,:));
Min=min(Vfenbu_s(34,:));
Length=(Max-Min)/nbins;
[nhist,c_points]=hist(Vfenbu_s(34,:),nbins);      %获取所绘制的频数直方图各区间的中心以及各区间对应的元素的个数

%求解各节点电压越限概率，依据欧洲标准EN 50160的规定
P_Vofflimits_MC=ProbMC(Vfenbu_s,Nodes,daishu);
P_Vofflimits_MC=P_Vofflimits_MC(2:34,1)           %33个PQ节点的电压越限概率



%求电压样本的八阶原点矩
for i=1:size(Vfenbu_s,1)
    for k=1:8
        Vk=Vfenbu_s(i,:).^k;
        Vai(i,k)=mean(Vk);
    end
end
%求电压样本的八阶半不变量
Vr=zeros(size(Vfenbu_s,1),8);
Vr(:,1)=Vai(:,1);
for k=2:8
    sigma=0;
    for j=1:k-1
        sigma=nchoosek(k-1,j)*Vai(:,j).*Vr(:,k-j)+sigma;
    end
    Vr(:,k)=Vai(:,k)-sigma;
end
for i=2:8
    Vr_st(:,i)=Vr(:,i)./Vr(:,2).^(i/2);
end

alpha=0:0.00001:1;
z=norminv(alpha,0,1);                %标准正态分布的上alpha分位数

for i=1:length(Vpqi)
    xxV(i,:)=z+(z.^2-1)*gv_s(i,3)/6+(z.^3-3*z)*gv_s(i,4)/24-...
    (2*z.^3-5*z)*(gv_s(i,3)^2)/36+(z.^4-6*z.^2+3)*gv_s(i,5)/120-...
    (z.^4-5*z.^2)*gv_s(i,2)*gv_s(i,3)/24+(12*z.^4-53*z.^2)*(gv_s(i,3)^3)/324;
    yyV(i,:)=xxV(i,:).*sqrt(DetaVpqik_s(i,2))+Vmui(i,1);
end

%求有功样本的八阶原点矩
for i=1:size(P_linefenbu_s,1)
    for k=1:8
        Pk=P_linefenbu_s(i,:).^k;
        Pai(i,k)=mean(Pk);
    end
end
%求有功样本的八阶半不变量
Pr=zeros(size(P_linefenbu_s,1),8);
Pr(:,1)=Pai(:,1);
for k=2:8
    sigma=0;
    for j=1:k-1
        sigma=nchoosek(k-1,j)*Pai(:,j).*Pr(:,k-j)+sigma;
    end
    Pr(:,k)=Pai(:,k)-sigma;
end
alpha=0:0.00001:1;
z=norminv(alpha,0,1);                %标准正态分布的上alpha分位数

for i=1:length(kp_s)
    xxP(i,:)=z+(z.^2-1)*gp_s(i,3)/6+(z.^3-3*z)*gp_s(i,4)/24-...
    (2*z.^3-5*z)*(gp_s(i,3)^2)/36+(z.^4-6*z.^2+3)*gp_s(i,5)/120-...
    (z.^4-5*z.^2)*gp_s(i,2)*gp_s(i,3)/24+(12*z.^4-53*z.^2)*(gp_s(i,3)^3)/324;
    yyP(i,:)=xxP(i,:).*sqrt(kp_s(i,2))+pmui(i,1);
end

%求无功样本的八阶原点矩
for i=1:size(Q_linefenbu_s,1)
    for k=1:8
        Qk=Q_linefenbu_s(i,:).^k;
        Qai(i,k)=mean(Qk);
    end
end
%求无功样本的八阶半不变量
Qr=zeros(size(Q_linefenbu_s,1),8);
Qr(:,1)=Qai(:,1);
for k=2:8
    sigma=0;
    for j=1:k-1
        sigma=nchoosek(k-1,j)*Qai(:,j).*Qr(:,k-j)+sigma;
    end
    Qr(:,k)=Qai(:,k)-sigma;
end
alpha=0:0.00001:1;
z=norminv(alpha,0,1);                %标准正态分布的上alpha分位数

for i=1:length(kq_s)
    xxQ(i,:)=z+(z.^2-1)*gq_s(i,3)/6+(z.^3-3*z)*gq_s(i,4)/24-...
    (2*z.^3-5*z)*(gq_s(i,3)^2)/36+(z.^4-6*z.^2+3)*gq_s(i,5)/120-...
    (z.^4-5*z.^2)*gq_s(i,2)*gq_s(i,3)/24+(12*z.^4-53*z.^2)*(gq_s(i,3)^3)/324;
    yyQ(i,:)=xxQ(i,:).*sqrt(kq_s(i,2))+qmui(i,1);
end
P_Vofflimits_CMCF=ProbCMCF(yyV,alpha,length(DetaVpqik_s))


%求解各节点电压越限概率，依据欧洲标准EN 50160的规定
P_Vofflimits_PV_CF=ProbCMCF(ybV,alpha,length(Vpqi));
difference1=P_Vofflimits_CMGC-P_Vofflimits_MC        %比较两种方法求出的电压越限概率差距
difference2=P_Vofflimits_CMCF-P_Vofflimits_MC

%% 画图对比分析-------------------
%绘制加入DG节点33的概率密度曲线
figure
bar(c_points,nhist/daishu/Length,'grouped','white');
hold on;
[f,xi]=ksdensity(Vfenbu_s(34,:));
plot(xi,f,'r-.','LineWidth',2);
hold on;
plot(x,f3_s(33,:),'b-','LineWidth',2);
hold on;
[f,xi]=ksdensity(ybV(33,:));
plot(xi,f,'g--','LineWidth',2);
hold on;
title('节点33电压 概率密度曲线','FontSize',16);
axis auto normal;
xlabel('电压(p.u.)','FontSize',16);
ylabel('概率密度','FontSize',16);
legend('MonteCarlo模拟频率直方图','MonteCarlo模拟概率密度估计','半不变量+Gram-Charlier级数展开','半不变量+Cornish-Fisher级数展开');
grid on;
box on;
hold off;


%绘制加入DG节点33的累积分布曲线
figure
plot(x,f2_s(33,:),'b-','LineWidth',2);
hold on;
plot(VBi_s,VCi_s,'r-.','LineWidth',2);
hold on;
plot(ybV(33,:),alpha,'g--','LineWidth',2);
hold on;
axis auto normal;
title('节点33电压 累积概率曲线','FontSize',16);
xlabel('节点电压(p.u.)','FontSize',16);
ylabel('累积概率','FontSize',16);
legend('半不变量+Gram-Charlier级数展开','MonteCarlo模拟结果','半不变量+Cornish-Fisher级数展开');
grid on;
box on;
hold off;


% 绘制线路有功概率分布情况
figure
plot(Pi,fp3_s(33,:),'b-','LineWidth',2);
hold on;
[f_pline_s,xi_pline_s]=ksdensity(P_linefenbu_s(33,:));
plot(xi_pline_s,f_pline_s,'r-.','LineWidth',2);
hold on;
[f,xi]=ksdensity(ybP(33,:));
plot(xi,f,'g--','LineWidth',2);
hold on;
title('线路31-33 有功概率密度曲线','FontSize',16);
axis auto normal;
xlabel('有功功率(p.u.)','FontSize',16);
ylabel('概率密度','FontSize',16);
legend('半不变量+Gram-Charlier级数展开','MonteCarlo模拟结果','半不变量+Cornish-Fisher级数展开');
grid on;
box on;
hold off;

figure
plot(Pi,fp2_s(33,:),'b-','LineWidth',2);
hold on;
plot(pBi_s,pCi_s,'r:','LineWidth',2);
hold on;
plot(ybP(33,:),alpha,'g--','LineWidth',2);
hold on;
title('线路31-33 有功累积分布曲线','FontSize',16);
axis auto normal;
xlabel('有功功率(p.u.)','FontSize',16);
ylabel('累积概率','FontSize',16);
legend('半不变量+Gram-Charlier级数展开','MonteCarlo模拟结果','半不变量+Cornish-Fisher级数展开');
grid on;
box on;


%绘制线路无功概率分布情况
figure
plot(Qi,fq3_s(33,:),'b-','LineWidth',2);
hold on;
[f_qline_s,xi_qline_s]=ksdensity(Q_linefenbu_s(33,:));
plot(xi_qline_s,f_qline_s,'r-.','LineWidth',2);
hold on;
[f,xi]=ksdensity(ybQ(33,:));
plot(xi,f,'g--','LineWidth',2);
hold on;
title('线路31-33 无功概率密度曲线','FontSize',16);
axis auto normal;
xlabel('无功功率(p.u.)','FontSize',16);
ylabel('概率密度','FontSize',16);
legend('半不变量+Gram-Charlier级数展开','MonteCarlo模拟结果','半不变量+Cornish-Fisher级数展开');
grid on;
box on;

figure
plot(Qi,fq2_s(33,:),'b-','LineWidth',2);
hold on;
plot(qBi_s,qCi_s,'r:','LineWidth',2);
hold on;
plot(ybQ(33,:),alpha,'g--','LineWidth',2)
hold on;
title('线路31-33 无功累积分布曲线','FontSize',16);
axis auto normal;
xlabel('无功功率(p.u.)','FontSize',16);
ylabel('累积概率','FontSize',16);
legend('半不变量+Gram-Charlier级数展开','MonteCarlo模拟结果','半不变量+Cornish-Fisher级数展开');
grid on;
box on;

 

序列运算理论

        为了将原本复杂的卷积计算转换为简单的算术运算，康重庆教授创建了序列运算理论，在此基础上扩展了概率性序列运算方法，从而形成了全新的电力系统不确定性分析框架。

        基于序列运算的概率直流潮流是以简化序列卷积为出发点，其自定义的卷和、序乘等运算方法计算简单，在效率上具有很大的优势。但由于序列的建立和运算的定义都要满足全新的规则和要求，从而限制了该方法的大规模推广应用。

解析法概率潮流展望

        解析法概率潮流的研究重点在卷积运算的高效处理上，快速傅里叶运算难以处理大规模电力系统的多变量计算，序列运算理论体系架构的特殊性使之难以在短时间内大规模推广，因此半不变量法是解析法中研究的热点算法。半不变量法在处理正态分布的变量相关性上己有较多的研究，该算法的研究重点是非正态分布的变量相关性处理，改善各阶矩和半不变量的计算精度与计算效率，获得更为准确的概率密度函数。

总结

       电力系统中不确定因素的增加使得概率潮流成为研究热点之一。从概率潮流算法原理的角度出发，将电力系统概率潮流算法分为模拟法、近似法和解析法3类，对各算法的研究现状和发展进行了综述和总结，并对各类算法的发展趋势提供了分析思路。

        概率潮流与确定性潮流本质上的区别在于输入变量的不确定性，计算目标是合理反映随机分布本身的概率特性对于系统潮流的影响，因此，概率潮流的研究目的可概括为随机变量的概率输入及其对潮流输出的概率贡献。模拟法分析输入随机变量的抽样方法，近似法关注输出与输入之间的转化关系，解析法重点讨论系统方程卷积的计算与简化，而3类方法中随机变量概率分布的表征和相关性分析都是关键因素。

        而对新环境下电力系统的全新特征，未来概率潮流算法研究的总体发展趋势应该包含以下几个主要方向：可再生能源功率输出的概率密度函数并不一定精确满足特定的典型分布，如何生成适用于概率潮流计算的准确的可再生能源功率输出的概率分布非常重要，目前己有的核密度估计法。具有很大的研究空间；在大规模新能源集中并网的趋势下，如何处理区域内新能源功率输出的复杂相关性变得日益重要；电动汽车作为电网互动性负荷，与电网的功率交互是双向的，如何对双向功率负荷进行概率建模也是研究难点之一；电力市场机制下的电源和负荷具备开放性特点，基于市场特征和需求响应的最优概率潮流分析具备很高的研究意义。