技术领域
[0001] 本
发明属于综合能源系统可靠性技术领域，具体设计一种综合能源系统概率能流的计算方法。
背景技术
[0002] 近年来，全球能源形势越来越严峻，与能源消耗相关的环境问题如大气变暖、天气污染问题也逐渐恶化，能源问题已经影响全人类生存与发展。探索可行、高效、经济的能源格局，减少与能源消耗对环境的影响已经成为世界的共识。综合能源系统能够对各个能源进行优势互补，协调运行，提高
能量利用效率和经济效益，减少能源消耗对环境的影响，因而得到了世界各国的广泛关注。
[0003] 综合能源系统能够整合各种形式的能源，因而不可避免地存在可靠性问题：新能源的接入给电
力系统带来了很大的不确定性问题；气网负荷也具有一定的
波动性和不确定性，同时各种形式的能量耦合环节也会给综合能源系统带来不确定性问题。因此，实际上综合能源系统的能量流呈现一定概率的形式，传统的稳态能流问题已经不能完全反应综合能源系统的状态信息，因此有必要进行综合能源系统的概率能流计算。
[0004] 要对综合能源系统的概率能流进行计算，必先计算该系统的稳态能量流。目前，
牛顿法是稳态能量流计算的通用方法。此外，根据雅可比矩阵的结构，通过
迭代各种能量系统使用该方法获得的解可以分为两种形式：独立牛顿迭代法和统一牛顿迭代法。独立牛顿迭代法中每个能量系统独立形成自己的雅可比矩阵，并分别使用牛顿法进行计算，在能量耦合单元对
不平衡量迭代直到整个系统的能量流收敛。这种方法比较简单，但计算速度稍慢。统一牛顿迭代法是构造一个包含诸多能源系统的雅可比矩阵并进行统一的牛顿迭代，以此获得整个网络的综合能量流。这种方法雅克比矩阵的构造比较复杂，但由于不需要进行能量耦合单元的不平衡量迭代，因此收敛速度比较快。
[0005] 目前，综合能源系统的概率能流进行计算研究并不多，主要方法类似于传统电力系统中的概率潮流计算方法，包括蒙特卡罗模拟法，
解析法和近似方法。根据大数定理，蒙特卡罗模拟需要许多样本(3000个以上)才能获得良好的
精度，因此这种方法需要耗费相当长的时间，但是蒙特卡罗模拟法可以获得最准确的结果，因此，它是衡量其他方法精确度的参考。解析法通过线性化交流模型和卷积方法解决概率能流问题，在系统相对较大且需要考虑许多个组件的情况下，使用此方法将使得计算量呈现指数增长的趋势。因此，该方法不像其他两种方法那样广泛使用。近似方法是解决概率能流问题的重要方法，该方法主要包括一阶二阶矩法和点估计法。一阶二阶矩法具有快速的解决方案速度和良好的精度。由于采用泰勒展开法，只考虑线性项，忽略高阶项，因此，当系统高度非线性或概率变量非常不对称时，将导致比较大的误差。点估计法现在被广泛用于解决概率能流的问题，通过高斯积分技术，可以轻松获得估计点，该方法的求解速度非常高。但是原始的点估计法并不能处理相关变量的问题，并且估计点可能超过约束的边界，这将导致结果失去准确性。因此，点估计法仅限于特定应用。
[0006] 综上所述，目前的大多数方法存在计算速度和计算准确性的问题，本
专利提出一种结合Nataf变换并考虑边界约束的快速点估计法，对综合能源系统的概率能流问题进行研究，有效地提高了求解精度和计算效率。
发明内容
[0007] 本发明的目的是提供一种考虑相关性和边界约束的综合能源系统概率能流的计算方法。
[0008] 为实现上述目的，本发明采用的技术方案如下：
[0009] 一种综合能源系统概率能流的计算方法，根据综合能源系统的电力系统潮流和
天然气系统能流，得到综合能源系统稳态能流；
[0010] 对
风力发电单元、
光伏发电单元和气网负荷进行概率建模，得到
风力发电概率模型、光伏出力模型以及气网负荷概率模型，按照风力发电概率模型、光伏出力模型和气网负荷概率模型进行
采样，得到概率采样样本，利用Nataf变换对概率采样样本添加相关性，然后采用多点估计的方法对采样样本进行点估计，根据所得的估计点求解综合能源系统稳态能流，得到概率能流。
[0011] 本发明进一步的改进在于，电力系统潮流通过以下过程得到：
[0012] 令Pig＝0，Pic＝0， 对下面式子采用牛顿法进行求解，得到电力系统潮流；
[0013]
[0014]
[0015]
[0016]
[0017] 其中，Pis和 是电源注入
节点i的有功功率和
无功功率，Pil和 是第i节点负荷的有功功率和无功功率需求，Pic和 是
压缩机消耗的有功功率和无功功率，Pig和 则是燃气发
电机注入i节点的有功功率和无功功率，Pij与Qij是i节点与j之间的传输有功概率和无功概率，n是与i节点直接相连的节点数量； Pij和Qij是节点i与节点j之间输电线的有功和无功功率；
[0018] 其中，节点i与节点j之间输电线的有功功率Pij和节点i与节点j之间输电线的无功功率Qij如下：
[0019]
[0020]
[0021] 其中，Ui和Uj是节点i与节点j
电压，θij＝θi-θj，θi是
母线i的电压
角，θj是母线j的电压角；Gij和Bij分别为传输线电导和电纳。
[0022] 本发明进一步的改进在于，天然气系统能流通过以下过程得到：
[0023] 利用牛顿法对下式进行迭代求解；
[0024]
[0025] 其中，Qi为节点处i天然气的净注入量；fim、fin分别为向下游节点注入流量和自上游节点输出流量； Fj为压缩机的气体消耗；Gij为关联系数；
[0026] 在第k次迭代求解时，修正方程如下所示：
[0027]
[0028] 其中，
[0029]
[0030]
[0031] 式中， 为所求函数的误差向量；J(k)为第k次迭代的雅可比矩阵； 为第k次迭代的修正量向量，其中g＝1…i；
[0032] 当第k次迭代的修正量向量 中最大的一个小于一个给定值时，得到天然气系统能流；
[0033] 其中，当压缩机从节点i取气时，关联系数Gij取1，否则为0,若压缩机所需的能量取自于
电网，则关联系数Gij也取0。
[0034] 本发明进一步的改进在于，当为输气网时，向下游节点注入流量fim为：
[0035]
[0036] 式中，Wim为系数，pi与pm分别为i节点与m节点的气压；
[0037] 当为配气网时，向下游节点注入流量fim为：
[0038]
[0039] 其中，D与L为管道的直径和长度，SG为相对度f为
摩擦系数，pi与pm分别为i节点与m节点的气压；
[0040] 对于含压缩机支路，如果是燃气给压缩机供能，则压缩机的气体消耗为：
[0041]
[0042] 式中：η为压缩机效率；a、m、k1、k2为有关的常系数，pi与pm分别为i节点与m节点的气压。
[0043] 本发明进一步的改进在于，根据综合能源系统的电力系统潮流和天然气系统能流，得到综合能源系统稳态能流的具体过程如下：
[0044] 令节点i
燃气轮机的输出功率PG,i和压缩机消耗
电能 分别为燃气发电机注入i节点的有功功率Pig和压缩机消耗的有功功率Pic，燃气发电机注入i节点的无功功率压缩机消耗的无功功率 再次计算电力系统潮流，得到新的Pig和Pic，然后计算天然气系统能流，并不断重复上述计算过程，直到前后两次计算得到的Pig的差值和Pic的差值都小于0.001，得到综合能源系统稳态能流；
[0045] 其中，压缩机消耗电能 表达式如下：
[0046]
[0047] 其中，Qc,k为流过压缩机的气流量，pi与pm分别为i节点与m节点的气压，Bk为相应的系数，Qc,k为流过压缩机的气流量。
[0048] 本发明进一步的改进在于，风力发电概率模型为：
[0049]
[0050] 式中：kp为威布尔分布的形状参数；cp为尺度参数；
[0051] 光伏出力模型为：
[0052]
[0053] 式中，Gt,j、 分别为t时刻能源局域网j中光伏组件的光照强度、输出电功率；Gstc、 分别为光伏组件在标准条件下的光照强度、最大电功率输出；Tt,j、Tτ，分别为光伏设备在t时刻的
温度、参考温度；kT为功率温度系数；
[0054] 气网负荷概率模型如下:
[0055]
[0056] 其中，μ是气网负荷的均值，δ是气网负荷的方差，x是气网负荷消耗的气量。
[0057] 本发明进一步的改进在于，利用Nataf变换对概率采样样本添加相关性，具体过程如下：
[0058] 设输入随机向量X的相关系数矩阵为ρ,根据Nataf分布模型推导出相关系数矩阵ρ各分量ρij的计算表达式：
[0059]
[0060] 式中，ρ0ij为标准正态随机向量Y相关系数ρ的分量，σi和σj为正态空间的标准差，μi和μj为正态空间的均值，xi和xj为X中的元素，通过相关系数矩阵ρ各分量ρij的计算表达式得到xi＝F-1(Φ(yi))、 xj＝F-1(Φ(yj))，yi和yj为标准正态空间Y中的元素，ρ是已知原非正态空间的变量相关系数；
[0061] 通过反求相关系数矩阵ρ各分量ρij的计算表达式，得到所有的标准正态随机向量Y相关系数ρ的分量ρ0ij，并组成矩阵ρ0，ρ0为标准正态空间变量的相关系数；将标准正态空间变量的相关系数ρ0进行 Cholesky分解得：
[0062]
[0063] 式中，L0为相关系数矩阵经Cholesky分解得到的下三角阵；
[0064] 利用相关系数矩阵经Cholesky分解得到的下三角阵L0将相关的标准正态随机向量Y转换为随机变量向量X：
[0065]
[0066] 利用Nataf分布模型产生若干含相关性的风机、光伏采样样本X，X＝{x1,x2...xn}，[0067] 其中，Pise为未考虑风机、光伏出力时综合能源系统在平衡点处的电源有功功率，为未考虑风机、光伏出力时综合能源系统在平衡点处的电源无功功率，fi,me为未考虑天然气气网具有波动性负荷时综合能源系统在平衡点处的向下游节点注入流量，fi,ne为未考虑天然气气网具有波动性负荷时综合能源系统在平衡点处的向上游节点输出流量；
[0068] 将采样样本与平衡点处值相加得到sa
其中，Pi 为考虑风机、光伏
出力时综合能源系统在平衡点处的电源有功功率， 为考虑风机、光伏出力时综合能源系统在平衡点处的电源无功功率，fi,ma为考虑天然气气网具有波动性负荷时综合能源系统在平衡点处的向下游节点注入流量， fi,na为考虑天然气气网具有波动性负荷时综合能源系统在平衡点处的向上游节点输出流量；
[0069] 将 代替 得到若干综合能源系统的概率能流结果，将若干综合能源系统的概率能流结果进行方差、均值、概率
密度和累计概率密度的求解，得到综合能源系统概率能流。
[0070] 本发明进一步的改进在于，对于具有多个风电和光伏的综合能源系统，假设m是点估计的数量，n是响应函数的数量，m点估计的表达式为：
[0071] xij＝Gi-1(Φ(yij))i＝1:n j＝1:m
[0072]
[0073]
[0074] 根据m点估计的表达式进行计算，得到综合能源系统的概率统计量：
[0075]
[0076]
[0077] 其中，h(·)为响应函数，μh是综合能源系统能流的均值，D(h)是综合能源系统能流的方差，n为响应函数的数目，m为估计的点的数目，μ(·)为函数的均值，σh为函数h(·)的方差。
[0078] 本发明进一步的改进在于，根据综合能源系统的概率统计量，得到综合能源系统方差、均值的具体过程如下：
[0079] 将计算量从m×n降低为(m-1)×n+1：
[0080]
[0081]
[0082] 这里，
[0083]
[0084] μGn＝μ(Gn-1(Φ(y)))
[0085]
[0086]
[0087] 当k＝1时有：
[0088]
[0089] 其中，μhi为函数hi(·)的均值， 为函数Gi(·)的均值，pn为权重系数，Φ(yn为Nataf变换函数。
[0090] 本发明进一步的改进在于，当m＝2或3时，若估计点超过约束，则对超过限制的点进行约束，具体过程如下：
[0091] 假设x是任意空间的样本，μx为x的均值，δx是x的方差，λx,i是x的i阶中心矩，f(x)是x的概率密度函数，Hong的2点估计和3点估计法表述为下式：
[0092]
[0093]
[0094]
[0095] 式中，μx为x的均值，ξ2，i,i＝1,2，ξ3,i,i＝1,2,3代表
位置参数，p2，j,j＝1,2，p3，j,j＝1,2,3是权重系数， 和 为所估计的点：
[0096]
[0097]
[0098] 式中，n是变量的数量，ξ2，i,i＝1,2，ξ3,i,i＝1,2,3代表位置参数，p2，j,j＝1,2，p3，j,j＝1,2,3 是权重系数， 和 为所估计的点；
[0099] 2点估计和3点估计的均值和方差用下式进行计算：
[0100] μz＝p2，1h(Z(x2，1))+p2，2h(Z(x2，2))
[0101] μz＝p3，1h(Z(x3，1))+p3，2h(Z(x3，2))+p3，3h(Z(x3，3))
[0102]
[0103]
[0104] 其中，μz为响应函数h(·)的均值，D(x)为响应函数h(·)的方差，ξ2，i,i＝1,2，ξ3,i,i＝1,2,3代表位置参数，p2，j是权重系数；j＝1,2；p3，j是权重系数,j＝1,2,3；和 为所估计的点；
[0105] 假设F(x)是x原来风机概率模型中的约束函数，满足下式：
[0106]
[0107] 其中，xmin为约束最小值，xmax为约束最大值，F(·)为变换函数，δ、Lmax为形状参数，将x转换为与原来约束函数F(x)等价的形式：
[0108]
[0109] 然后对x′使用Hong的2点估计法或者3点估计法，得到估计点x′21,x′22和x′31,x′32,x′33，将这些点变换到原来空间为：
[0110]
[0111] 当k＝1时，得到实际估计点x21,x22和x31,x32,x33；
[0112] 将实际估计点x21,x22和x31,x32,x33应用下式进行计算，得到最终概率能流的均值、方差：
[0113] μz＝p2，1h(Z(x2，1))+p2，2h(Z(x2，2))
[0114] μz＝p3，1h(Z(x3，1))+p3，2h(Z(x3，2))+p3，3h(Z(x3，3))
[0115]
[0116]
[0117] 其中，μz为响应函数h(·)的均值，D(x)为响应函数h(·)的方差，ξ2，i,i＝1,2，ξ3,i,i＝1,2,3代表位置参数，p2，j,j＝1,2，p3，j,j＝1,2,3是权重系数， 和为所估计的点。
[0118] 与
现有技术相比，本发明具有的有益效果：本发明首先根据综合能源系统的电力系统潮流和天然气系统能流，得到综合能源系统稳态能流；对风力发电单元、光伏发电单元和气网负荷进行概率建模，得到风力发电概率模型、光伏出力模型以及气网负荷概率模型，按照风力发电概率模型、光伏出力模型和气网负荷概率模型进行采样，得到概率采样样本，利用Nataf变换对概率采样样本添加相关性，然后采用多点估计法对采样样本进行点估计，根据所得的估计点求解综合能源系统稳态能流，得到概率能流。本发明综合考虑了电力系统中的
风能、
太阳能以及天然气系统中气体负荷的不确定性的，建立了计算概率能流的实用化模型。考虑到现实中不同风力发电和不同太阳能发电之间的相关性，利用Nataf变换准确有效地解决相关问题，达到解决综合能源系统中的综合能流问题，并且有效地提高了求解精度和计算效率。考虑到蒙特卡罗方法在解决概率能流问题的效率低下，本发明使用多点估计方法提高了解决的效率，本发明利用快速的多点估计法，降低了计算复杂度和计算量。
[0119] 进一步的，对于估计点超限的情况，结合Hong的点估计法，本发明利用幂变换方法和等变约束变换法，改进了2点和3点估计并减少计算误差，达到解决综合能源系统中的综合能流问题，并且有效地提高了求解精度和计算效率。
附图说明
[0120] 图1为综合能源系统的总体图。
[0121] 图2为综合能源系统稳态能流计算
流程图。
[0122] 图3为Nataf变换以后风电样本之间和光伏样本之间相关性的控制。其中，(a)为风机1与风机2的采样点，(b)为风机1与风机3的采样点，(c)为风机2与风机3的采样点，(d)为光伏1与光伏2的采样点。
[0123] 图4为多点估计法与蒙特卡洛方法的对比。其中，(a)为气网气压均值，(b)为气压方差，(c)为电压均值，(d)为电压方差。
[0124] 图5为不经过Nataf变换风电样本之间和光伏样本之间的相关性。其中，(a)为风机1与风机2的采样点，(b)为风机1与风机3的采样点，(c)为风机2与风机3的采样点，(d)为光伏1与光伏2的采样点。
[0125] 图6为Nataf变换以后与不经过Nataf变换的蒙特卡洛模拟结果对比。其中，(a)为气网气压的均值， (b)为电网气压的均值。
[0126] 图7为有超过约束边界的点估计方法与蒙特卡洛模拟方法的结果对比。其中，(a)为天然气系统气压均值，(b)为天然气系统气压方差，(c)为电压均值，(d)为电压方差。
[0127] 图8为本发明的流程图。
具体实施方式
[0128] 为使本发明
实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0129] 参见图8，具体过程如下：
[0130] S1，采用独立构造电网、气网的雅克比矩阵的方法对电-气互相耦合的综合能源系统中的与风机出力、光伏出力、气网负荷有关的状态量拓展为概率量，并将概率量进行概率建模，再将概率量的概率空间运用Nataf变换的方法变换至正态空间，将正态空间中具有相关性概率样本通过查阅Gauss–Hermite积分表得到多点估计的估计点和相对应的权重，并将这些得到的点运用Nataf反变换变换到原来的概率空间。并且通过Choleskey分解公式对正态空间的概率样本加入相关性，从而得到综合能源系统的概率样本，运用Nataf 变换对综合能源系统的概率样本加入相关性，并结合多点估计法进行快速计算，得到综合能源系统的概率能流。
[0131] S2，在结合多点估计法对综合能源系统的概率能流进行快速计算的过程中，若其中的两点和三点估计方法出现超过约束边界的现象，则对两点和三点估计方法中出现的超界估计点进行约束，达到更加精确的计算概率能流的目的。
[0132] 采用牛顿法分别对电力系统潮流和天然气系统的能流进行求解，求解的不平衡量通过耦合系统进行迭代直到整体收敛，得到某一状态下的综合能源系统的稳态能流。
[0133] 下面为具体说明。
[0134] 综合能源系统的总体结构如图1所示，该系统包括电力系统和天然气系统。将综合能源系统稳态能流计算分为电力系统潮流计算部分和天然气系统能流计算部分，本发明采用牛顿迭代法对电力系统潮流计算部分和天然气系统能流计算部分分别进行计算，得到综合能源系统稳态能流。参见图2，具体计算过程如下：
[0135] 1、计算电力系统潮流
[0136] 电力系统潮流计算通过AC潮流建模，电力系统AC潮流方程的状态变量主要是节点电压U和相角θ，这两个状态量已知的时候便可以得到电力系统的潮流，电网中的基本潮流计算格式是：
[0137]
[0138]
[0139] 其中，Pij和Qij是节点i与节点j之间输电线的有功和无功功率。Ui和Uj是节点i与节点j电压，而θij＝θi-θj，θi是母线i的电压角，θj是母线j的电压角。Gij和Bij分别为传输线电导和电纳。
[0140] 节点的有功无功平衡则满足方程：
[0141]
[0142]
[0143]
[0144]
[0145] 其中，Pis和 是电源注入节点i的有功功率和无功功率，Pil和 是第i节点负荷的有功功率和无功功率需求，Pic和 是压缩机消耗的有功功率和无功功率，Pig和 则是燃气发电机注入i节点的有功功率和无功功率。Pij与Qij是i节点与j之间的传输有功概率和无功概率，n是与i节点直接相连的节点数量。
[0146] 令Pig＝0，Pic＝0， 对上述式子应用牛顿法进行求解，得到电力系统潮流。
[0147] 2、计算天然气系统能流
[0148] 天然气系统能流计算与能流计算与电力系统类似，气网各个量的基本方程满足一个类似的方程：
[0149]
[0150] 式中：Qi为节点处i天然气的净注入量；fim、fin分别为向下游节点注入和自上游节点输出流量；Fj为压缩机的气体消耗；当压缩机从节点i取气时，关联系数Gij取1，否则为0,若压缩机所需的能量取自于电网，则关联系数Gij也取0。
[0151] 当为输气网的时候，气压比较大，管道气流量通常按照Weymouth方程计算：
[0152]
[0153] 式中，系数Wim取决于管道长度、直径和摩阻系数及天然气压缩率等因素，pi与pm分别为i节点与m节点的气压。
[0154] 当为配气网的时候，气压比较低，管道中气体流量通常可以描述为：
[0155]
[0156] 其中，D与L为管道的直径和长度，SG为相对度f为摩擦系数。pi与pm分别为i节点与m节点的气压。
[0157] 对于含压缩机支路，如果是燃气给压缩机供能，则其压缩机的气体消耗为：
[0158]
[0159] 式中：η为压缩机效率；a、m、k1、k2为有关的常系数，关联系数Gij取决于气体流向，pi与pm分别为i节点与m节点的气压。
[0160] 天然气网潮流计算与电力系统潮流计算类似，以下做一个对比：
[0161] 表1天然气网潮流计算与电力系统潮流计算对比
[0162]
[0163] 利用牛顿法进行迭代求解，在第k次迭代求解时，其修正方程如下所示：
[0164]
[0165] 其中，
[0166]
[0167]
[0168] 式中， 为所求函数的误差向量；J(k)为第k次迭代的雅可比矩阵； 为第k次迭代的修正量向量， 等于前后两次迭代xg的差值。通过反复迭代上述式
子，当第k次迭代的修正向量 中最大的一个小于一个给定值时候，一般为
0.01或0.001，满足收敛条件，输出结果，得到天然气系统能流，即得到Qi，fim，fin和Fj。
[0169] 3、计算综合能源系统稳态能流
[0170] 计算综合能源系统稳态能流是通过耦合系统进行迭代计算的，其中耦合系统包括燃气轮机和压缩机。
[0171] 电力系统与天然气系统间耦合后电力系统通过压缩机向天然气系统供应能量，天然气系统通过燃气轮机为电力系统提供电能。燃气轮机在电力系统相当于电源，在天然气网里相当于负荷，其输入天然气流量与输出电功率之间的转化关系如下：
[0172]
[0173]
[0174] 其中，Hg,i为燃气轮机输入热量值，PG,i节点i燃气轮机的输出功率； 为天然气中节点m等效气负荷；GHV为固定的高热值；ag,i、bg,i、cg,i由燃气轮机的耗热率曲线决定。
[0175] 压缩机对于电力系统相当于负荷，其消耗
电能表达式如下：
[0176]
[0177] 其中，Qc,k为流过压缩机的气流量， 为消耗电能，pi与pm分别为i节点与m节点的气压。Bk为相应的系数，Qc,k为流过压缩机的气流量。
[0178] 令节点i燃气轮机的输出功率PG,i和压缩机消耗电能 分别为燃气发电机注入i节点的有功功率 Pig和压缩机消耗的有功功率Pic，燃气发电机注入i节点的无功功率压缩机消耗的无功功率 再次计算电力系统潮流，得到新的Pig和Pic，然后计算天然气系统能流，并不断重复上述计算过程，直到前后两次计算得到的Pig的差值和Pic的差值都小于0.001，得到综合能源系统稳态能流；
[0179] 由于风力发电单元、光伏发电单元和气网负荷呈现一定的概率特性，使得综合能源系统的能流表现为概率能流，本发明继而对风力发电单元、光伏发电单元和气网负荷进行概率建模，得到各自的概率模型，即风力发电概率模型、光伏出力模型以及气网负荷概率模型，具体过程如下：
[0180] 风力发电概率模型满足两参数威布尔分布，风速为v的风力发电概率模型f(v)为：
[0181]
[0182] 式中：kp为威布尔分布的形状参数，取值一般在1-3之间；cp为尺度参数，取值通常略高于平均风速。
[0183] 风机出力与风速有以下关系：
[0184]
[0185]
[0186] 其中，P和δ是相关的形状参数，PX是风机出力，xci为
切入风速，xco为
切出风速。
[0187] 光伏出力模型为：
[0188]
[0189] 式中，Gt,j、 分别为t时刻能源局域网j中光伏组件的光照强度、输出电功率；Gstc、 分别为光伏组件在标准条件下的光照强度、最大电功率输出；Tt,j、Tτ，分别为光伏设备在t时刻的温度、参考温度；kT为功率温度系数。其中，光伏组件在标准条件下的光照强度Gstc在一段时间内可近似看成Beta 分布，其概率密度函数如下：
[0190]
[0191] 式中：Gstc(t)、Gstcmax分别为t时段的实际光强和某一时段内最大光强；α和β分别为Bate分布的形状参数。
[0192] 气网负荷概率模型如下:
[0193]
[0194] 其中，μ是气网负荷的均值，δ是气网负荷的方差，x是气网负荷消耗的气量。
[0195] 按照上述风力发电概率模型、光伏出力模型和气网负荷概率模型进行采样，得到概率采样样本，由于风机之间和光伏之间存在相关性，所以本发明利用Nataf变换对概率采样样本添加相关性，具体过程如下：
[0196] 1)设输入随机变量向量
[0197] X＝{x1,x2...xn}
[0198] 其中，x1,x2...xn为采样样本，即为符合上述概率模型的随机数列。随机变量xt(t＝1,...n)的概率密度函数 fi(x1)和
累积分布函数Fi(xi)已知，通过等概率转换原则进行对应：
[0199]
[0200] 定义Y＝{y1,y2...yn}为标准正态随机向量，式中Φ(.)和Φ-1(.)分别为标准正态累积分布函数和逆累积分布函数。
[0201] 利用隐函数求导法则，可以推导出随机向量X的联合概率密度函数为
[0202]
[0203] 式中，fi(xi)，i＝1…n为概率密度函数，φ(.)为标准正态分布的概率密度函数，同时可以得到Nataf 分布模型：
[0204]
[0205] 这个分布模型被称为Nataf分布模型，可以应用于均值为0、方差为1及相关系数为ρ0的n维标准正态分布。
[0206] 2)控制随机向量变换和相关性
[0207] 设输入随机向量X的相关系数矩阵为ρ,根据上面的Nataf分布模型可推导出相关系数矩阵ρ各分量ρij的计算表达式：
[0208]
[0209] 式中，ρ0ij为标准正态随机向量Y相关系数ρ的分量，σi和σj为正态空间的标准差，μi-1和μj为正态空间的均值，xi和xj为X中的元素，通过上述式子可以知道xi＝F (Φ(yi))、xj＝F-1(Φ(yj))，yi和yj为标准正态空间Y中的元素，ρ是已知原非正态空间的变量相关系数。
[0210] 通过反求上式解非线性方程，可以得到所有的标准正态随机向量Y相关系数ρ的分量ρ0ij，将这些组成矩阵ρ0，ρ0为标准正态空间变量的相关系数。显然ρ0是一对称矩阵,将标准正态空间变量的相关系数ρ0进行Cholesky分解可得：
[0211]
[0212] 式中，L0为相关系数矩阵经Cholesky分解得到的下三角阵。利用相关系数矩阵经Cholesky分解得到的下三角阵L0可将相关的标准正态随机向量Y转换为随机变量向量X：
[0213]
[0214] 利用上述Nataf变换的方法(即Nataf分布模型)可以产生大量含相关性的风机、光伏采样样本X， X＝{x1,x2...xn}， Pise为未考虑风机、光伏出力时综合能源系统在平衡点处的电源有功功率， 为未考虑风机、光伏出力时综合能源系统在平衡点处的电源无功功率， fi,me为未考虑天然气气网具有波动性负荷时综合能源系统在平衡点处的向下游节点注入流量，fi,ne为未考虑天然气气网具有波动性负荷时综合能源系统在平衡点处的向上游节点输出流量；fim、fin分别为向下游节点注入和自上游节点输出流量；将这些将采样样本与平衡点处值相加，将相加得到的值加入电力系统潮流计算中，得到 其中，Pisa为考虑风
机、光伏出力时综合能源系统在平衡点处的电源有功功率， 为考虑风机、光伏出力时综合能源系统在平衡点处的电源无功功率，fi,ma为考虑天然气气网具有波动性负荷时综合能源系统在平衡点处的向下游节点注入流量，fina为考虑天然气气网具有波动性负荷时综合能源系统在平衡点处的向上游节点输出流量，将这些 代替之前的
并利用蒙特卡洛模拟的方法加入综合能源系统能稳态流的求解中，就可
以得到许多综合能源系统的概率能流结果，将这些结果运用统计的方法进行方差、均值、概率密度和累计概率密度的求解，便得到综合能源系统概率能流的计算结果。
[0215] 由于蒙特卡洛需要大量样本进行计算，效率低下，本发明利用Nataf变换与多点估计结合的方法进行求解，使得求解综合能源系统稳态能流的次数由几千次降为几十次，大大提高了求解效率。
[0216] 3)采用多点估计法计算概率能流
[0217] 由上述Nataf变换的结果 可以得到X＝{x1,x2...xn}，其中：
[0218] xj＝G-1(Φ(yj))j＝1:k
[0219] yj为上述相关的标准正态随机向量Y＝{y1,y2...yn}中的样本，假设m是点估计的数量，m-点估计的表达式为：
[0220] yi＝μy+ξiσy,i＝1,2…m
[0221] μx＝p1G-1(Φ(y1))+p2G-1(Φ(y2))+…pmG-1(Φ(ym))
[0222] 其中，Φ代表正态空间的样本，μy是未进行Nataf变换前正态空间样本的均值，σy是未进行Nataf 变换前正态空间样本的方差，ξi是第i点估计的位置参数，pi是第i点估计的权重。ξi和pi可以查阅下面的Gauss–Hermite积分公式得到：
[0223] 表1 Gauss–Hermite积分公式表
[0224]
[0225] 然后，将估计的m(m＝2、3、5或7)个点进行确定性稳态能流计算，得到综合能源系统的概率统计量：
[0226] μh＝p1h(G-1(Φ(y1)))+…+pmh((G-1(Φ(ym)))
[0227]
[0228] 其中，h(·)为响应函数，m为估计的点的数目，μ(·)为函数的均值，μh是综合能源系统能流的均值，D(h) 是综合能源系统能流的方差，即代表上述综合能源系统的稳态能流计算。
[0229] 对于综合能源系统这个具有多个风电和光伏的系统，实际上相当于多个响应函数，假设m是点估计的数量，n是响应函数的数量，m点估计的表达式为：
[0230] xij＝Gi-1(Φ(yij))i＝1:n j＝1:m
[0231]
[0232]
[0233] 然后，进行确定性稳态能流计算，得到综合能源系统的概率统计量:
[0234]
[0235]
[0236] 其中，h(·)为响应函数，即代表上述综合能源系统的稳态能流计算。μh是综合能源系统能流的均值，D(h)是综合能源系统能流的方差，n为响应函数的数目，m为估计的点的数目，μ(·)为函数的均值，σh为函数h(·)的方差；上述得到综合能源系统的概率统计量可以通过快速的点估计法降低计算次数，提高计算效率，得到综合能源系统方差、均值这些统计量，具体过程如下：
[0237] 将计算量从m×n降低为(m-1)×n+1：
[0238]
[0239]
[0240] 这里，
[0241]
[0242] μGn＝μ(Gn-1(Φ(y)))
[0243]
[0244]
[0245] 当k＝1时候有：
[0246]
[0247] 其中，μhi为函数hi(·)的均值， 为函数Gi(·)的均值，pn为多点估计法得到的权重系数，Φ(yn为 Nataf变换函数。
[0248] 由于上述方法中m＝2、3的时候，即用2点或者3点估计可能会出现估计点超过约束的情况，因此必须对超过限制的点进行约束，本发明结合Hong的点估计法实现该目标。
[0249] 假设x是任意空间的样本，μx为x的均值，δx是x的方差，λx,i是x的i阶中心矩，f(x)是x的概率密度函数，Hong的2点估计和3点估计法可表述为以下的式子：
[0250]
[0251]
[0252]
[0253] 式中，μx为x的均值，ξ2，i,i＝1,2，ξ3,i,i＝1,2,3代表位置参数，p2，j,j＝1,2，p3，j,j＝1,2,3是权重系数， 和 为所估计的点：
[0254]
[0255]
[0256] 式中，n是变量的数量，ξ2，i,i＝1,2，ξ3,i,i＝1,2,3代表位置参数，p2，j,j＝1,2，p3，j,j＝1,2,3 是权重系数， 和 为所估计的点：
[0257] 2点估计和3点估计的均值和方差可以用如下的式子进行计算：
[0258] μz＝p2，1h(Z(x2，1))+p2，2h(Z(x2，2))
[0259] μz＝p3，1h(Z(x3，1))+p3，2h(Z(x3，2))+p3，3h(Z(x3，3))
[0260]
[0261]
[0262] 其中，μz为响应函数h(·)的均值，D(x)为响应函数h(·)的方差，ξ2，i,i＝1,2，ξ3,i,i＝1,2,3代表位置参数，p2，j,j＝1,2，p3，j,j＝1,2,3是权重系数， 和为所估计的点。
[0263] 这些点估计法适用于任意空间的概率变量，当出现估计点超过边界时候，本发明提出一种变换方法保证一定能够找到被约束的xi，这种方法称作等约束变换法，这种方法的优点是估计点肯定在约束范围内，因此它是一种简单且更通用的方法。
[0264] 假设F(x)是x原来风机概率模型中的约束函数，满足如下的关系：
[0265]
[0266] 其中，xmin为约束最小值，xmax为约束最大值，F(·)为变换函数，δ、Lmax为形状参数，将x转换为与原来约束函数F(x)等价的形式：
[0267]
[0268] 然后对x′使用Hong的2点估计法或者3点估计法，得到估计点x′21,x′22和x′31,x′32,x′33，将这些点变换到原来空间为：
[0269]
[0270] 很容易看出当k＝1的时候，该约束就是原来的约束，得到实际估计点x21,x22和x31,x32,x33。
[0271] 将这些得到的点应用下面式子进行计算：
[0272] μz＝p2，1h(Z(x2，1))+p2，2h(Z(x2，2))
[0273] μz＝p3，1h(Z(x3，1))+p3，2h(Z(x3，2))+p3，3h(Z(x3，3))
[0274]
[0275]
[0276] 其中，μz为响应函数h(·)的均值，D(x)为响应函数h(·)的方差，ξ2，i,i＝1,2，ξ3,i,i＝1,2,3代表位置参数，p2，j,j＝1,2，p3，j,j＝1,2,3是权重系数， 和为所估计的点，得到最终概率能流的均值、方差这些统计量。
[0277] 下面通过两个实施例对本发明实施例进行进一步说明，以下仅为本发明实施例的实施例，并不以此为限。
[0278] 实施例1
[0279] 本发明算例的综合能源系统如图1所示，其中风电系统的风速满足威布尔分布，太阳能系统的光强满足Beta分布，假设风速为正常风速(8.5m/s)，气网负荷3为波动负荷，满足正态分布，风电系统的风速分布和发电函数、光伏光强满足如下关系：
[0280]
[0281]
[0282] f(x)＝20(x)2(1-x)0≤x≤1
[0283] 运用Nataf变换给样本加入相关性，图3为Nataf变换后样本之间的分布图。从图3中(a)可以看出，风机1和风机2的相关性为0.32，即ρ12＝0.32，从图(b)中可以看出，风机1和风机3的相关性为0.45，即ρ13＝0.45，从图(c)中可以看出，风机2和风机3的相关性为0.37，即ρ23＝0.37，从图(d)中可以看出，光伏1和光伏2的相关性为0.37，即ρ＝0.37，从图3可以看出来，样本之间的分布图呈现椭圆关系，说明Nataf变换对样本相关性的控制结果很好。
[0284] 图4为各种点估计法与蒙特卡洛模拟结果的对比，从图4中(a)、(b)、(c)、(d)的对比结果可以看出来，点估计法的精度比较高，特别是7点估计法具有更加好的效果。
[0285] 表2和表3计算了各种点估计法的相对误差平均值，可以看出电网的求解精度比较高。
[0286] 表2气网的平均误差
[0287]
[0288] 表3电网的平均误差
[0289]
[0290] 将通过Nataf变换加入相关性的样本与不经过Nataf变换的样本进行对比，不经过Nataf变换的样本分布如图5所示，由图5中(a)可以看出，风机1与风机2的相关性为0.012，即ρ12＝0.012，从图(b)中可以看出，风机1与风机3的相关性为0.005，即ρ23＝0.005，从图(c)中可以看出，风机2和风机3的相关性为0.043，即ρ13＝0.043，从图(d)中可以看出，光伏1与光伏2的相关性为0.003，即ρ＝0.003。由图5的样本分布可以看出来，不同风机之间和不同光伏之间的样本分布呈现圆形或者杂乱的方形，因此样本之间相关性很低，可以看成不具有相关性，其仿真结果与加入Nataf变换的仿真结果对比如图6中(a) 和(b)所示，从图6可以看出，不加入Nataf变换对电网影响比气网大一些。
[0291] 实施例2
[0292] 本发明算例的综合能源系统仍然如图1所示，其中风电系统的风速满足威布尔分布，太阳能系统的光强满足Beta分布，此时天气处于无风状态，风速为非正常风速(7.14m/s)，气网负荷3为波动负荷，满足正态分布，风电系统的风速分布和发电函数、光伏光强满足如下关系：
[0293]
[0294]
[0295] f(x)＝20(x)2(1-x)0≤x≤1
[0296] 通过仿真计算可以得到多点估计的结果如图7所示，从图7中(a)、(b)、(c)、(d)可以看出来，2 点估计和3点估计由于存在超过约束的估计点，导致精度比较低，表4和表5进一步说明了该问题。
[0297] 表4气网的平均误差
[0298]
[0299] 表5电网的平均误差
[0300]
[0301] 运用幂变换法和等约束变换法对估计点进行约束，得到表7、表8和表9的结果：
[0302] 表6不同方法估计点的比较
[0303]
[0304] 表7等约束变换法约束估计点后气网的平均相对误差
[0305]
[0306] 表8等约束变换法约束估计点后电网的平均相对误差
[0307]
[0308] 从表6、表7和表8可以看出来，等约束变换法后限制估计点在约束以内，精度得到了很大的提高
[0309] 本发明实施例提供的一种综合能源系统考虑约束和相关性的概率能流计算方法，通过将Nataf变换与多点估计相结合，给概率样本加入相关性，并同时大大提高综合能源系统概率能流的计算效率，进一步，通过等约束变换法，将估计点限制在约束内，提高概率能流计算的准确性。
[0310] 本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的
硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
[0311] 通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到各实施方式可借助
软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件。基于这样的理解，上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在计算机可读存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，
服务器，或者网络设备等)执行各个实施例或者实施例的某些部分所述的方法。
[0312] 最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行
修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
