CN107204618A - 基于数字交错技术的拟蒙特卡洛随机潮流计算方法 - Google Patents

Abstract

本发明公开了一种计及风速相关性的电力系统拟蒙特卡洛潮流计算。主要是将数字交错技术应用在了电力系统潮流计算中，有效地提高了拟蒙特卡洛潮流计算的效率和精度。首先建立了多个风电场出力和负荷的概率分布模型，采用等概率转换和Nataf变换得到在标准正态分布下的等效相关系数矩阵；接着采用基于数字交错技术的Sobol序列产生S维随机数，其中数字交错技术是对原始Sobol序列的生成矩阵进行了行变换，有效优化了高维随机数的分布；最后，采用Cholesky分解技术获得具有相关性的随机输入变量样本进行随机潮流计算，得到节点电压幅值和支路注入有功和无功。该算法能够在多个风电场接入系统的情况下，进行拟蒙特卡洛潮流计算，效率高，计算精度高。

Description

基于数字交错技术的拟蒙特卡洛随机潮流计算方法

技术领域

本发明涉及一种拟蒙特卡洛随机潮流计算方法，尤其涉及一种基于数字交错技术的拟蒙特卡洛随机潮流计算方法，属于供电技术领域。

背景技术

风力发电因分布广泛、使用方便、无污染等优点，受到人们的广泛关注。随着大规模风电场的并网运行，电力系统的运行特性由弱随机性转变为强随机性。并且对于地理距离较近的多个风电场，各风电场在相同时刻下基本处于同一风速带，风功率将具有一定的空间相关性。这种空间相关性将加剧电网中的随机波动性。因此研究含多个风电场的电力系统随机潮流计算具有重要意义。随机潮流计算可以计及各种不确定因素对电力系统的影响，综合评估电网运行的薄弱点及可靠性。目前，随机潮流计算的方法主要包括：卷积法、点估计法和蒙特卡洛仿真法。卷积法利用输入变量的卷积运算得到输出变量的概率分布，该方法思路简单，但计算量大。点估计法根据输入变量的分布特性求取输出变量的数字特征，求出各阶矩后利用级数展开求得概率分布；该方法速度快，但是计算高阶矩时精度较低。蒙特卡洛仿真法通过随机采样得到输入变量的样本，再对每个采样点进行确定性潮流计算，统计获得各状态量的分布情况。整个过程原理简单且适用性广，但是要获得较高的精度必须牺牲计算量和计算时间。因此，为了获得较高精度同时提高计算效率，国内外学者提出了多种改进方法。以拉丁超立方抽样为代表的伪随机抽样提高了采样效率，但是不能保证序列的低偏差性。而采用了低偏差性序列后的拟蒙特卡洛仿真法可以保证随机数的均匀性，但是当序列维数较高时，相邻两维随机数会出现聚集现象，具有一定的规律性，这将影响拟蒙特卡洛仿真效果。

发明内容

本发明要解决的技术问题是提供一种基于数字交错技术的拟蒙特卡洛随机潮流计算方法。

为解决上述技术问题，本发明采用的技术方案是：

一种基于数字交错技术的拟蒙特卡洛随机潮流计算方法，包括以下具体步骤：

步骤1：建立各风电场出力和负荷的概率分布模型；

步骤2：基于等概率转换原则和Nataf变换计算标准正态分布下的相关系数矩阵C_Y；

步骤3：基于数字交错技术的Sobol序列产生分布均匀的S维随机数；由以下具体步骤组成：

本发明对Sobol序列的生成矩阵C_m,j进行处理，定义了相应的函数D_d()，使得高维随机数分布更均匀，有效提高了拟蒙特卡洛法仿真效率。

步骤301：定义最高次项为s_j的多项式：

其中系数为0或1；j代表随机数的维数，j≤d×s；d≥1，s≥1，均为整数，

步骤302：计算直接数{v_1,j,v_2,j,…}：

建立迭代关系式：

其中⊕表示二进制异或运算；初始值任意给定，只需要满足：m_k,j,1≤k≤s_j为奇数，并且小于等于2^k；直接数{v_1,j,v_2,j,…}计算公式为：

步骤303：形成矩阵C_m,j，

其中矩阵的元素定义如下：v_k,j＝(0.v_k,j,1v_k,j,2v_k,j,3…)₂。C_m,1,…,C_m,d×s代表d×s个m×m的矩阵；

步骤304：生成随机数x_i,j：对于第j维随机数，共有2^m个随机数；给定一整数i,满足1≤i≤2^m，对i进行二进制分解：则对于第j维中的第i个随机数x_i,j，其生成的原则如下：

x_i,j＝(b^-1,b^-2,…,b^-m)·C_m,j(i₁,i₂,…,i_m)^T＝ξ_j,1b^-1+ξ_j,2b^-2+…+ξ_j,mb^-m

步骤305：对于步骤304的d×s维随机数x＝(x₁,x₂,…,x_d×s)∈[0,1)^ds数字交错处理:

经过上式的处理后，可以得到s维随机数D_d,1,D_d,2,…,D_d,(s-1),D_d,s，简记为D＝(d₁,d₂,…,d_s)；

步骤4：采用Cholesky分解技术获得具有相关性的随机输入变量样本进行随机潮流计算，得到节点电压幅值和支路注入有功和无功。

进一步，步骤1由以下具体步骤组成：

步骤101：采用正态分布拟合负荷和有功功率的概率密度函数，标准差为期望的5％，无功功率按照给定的功率因数跟随有功功率变化。

步骤102：采用Weibull概率分布模型拟合各个风电场的风功率的概率密度函数：

式中，p_i为第i风电场的风功率i＝1,2,…,S；λ_i>0是对应的比例参数；k_i>0是对应的形状参数。

进一步，步骤2由以下具体步骤组成：

步骤201：P(p₁,p₂,…,p_S)表示维数为S的风功率输入向量；风功率p_i和p_j的相关系数为ρ_ij；f(p_i)和F(p_i)分别为风功率p_i的概率密度函数和累计概率分布函数；根据等概率准则和Nataf变换，建立起p_i和y_i对应的转换关系：φ(y_i)＝F(p_i)；φ表示标准正态分布随机变量Y＝(y₁,y₂,…,y_S)的累积概率分布，其中y_i和y_j间的相关系数为ρ_yij；

步骤202：采用二重积分法可获得随机向量Y的相关系数矩阵C_Y：

式中：u_pi,u_pj表示为风功率p_i,p_j的均值。表示具有相关系数ρ_yij的标准正态随机变量y_i和y_j的联合概率密度函数。

进一步，步骤305也可以通过对生成矩阵C_m,1,…,C_m,d×s的行向量进行变换形成新的生成矩阵C’_m,1,C’_m,2,…,C’_m,d。

进一步，步骤4由以下具体步骤组成：

步骤401：将步骤202所得到的相关系数矩阵C_Y进行Cholesky分解，得到下三角矩阵L：C_Y＝LL^T；

步骤402：基于步骤3产生的样本矩阵其中d_i∈U(0,1)；由式y_i＝φ^-1(d_i)i＝1,2,…,s得到标准正态分布随机向量Y＝(y₁,y₂,…,y_S)；

步骤403：根据Y^corr＝LY求取矩阵Y^corr，矩阵Y^corr的协方差矩阵即为C_Y；根据式依次转换矩阵Y^corr中的各元素，得到相关随机输入样本P；对每组样本分别进行潮流计算，利用统计学方法得到节点电压和支路潮流的数字特征和概率分布。

采用上述技术方案所取得的技术效果在于：

1)本发明适用于解决多个风电场接入对系统节点电压和支路功率的影响；

2)本发明克服了传统蒙特卡洛模拟方法收敛速度慢，效率低的缺点，通过采用了低偏差序列，有效提高了仿真速度和精度；

3)本发明通过引入数字交错函数D_d()，有效解决了低偏差序列在生成高维随机数时存的聚集现象，提高拟蒙特卡洛的仿真效率。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明的流程图；

图2是Sobol序列产生的9维和10维随机数分布图；

图3是本发明实施例中Sobol序列产生的9维和10维随机数分布图；

图4是本发明中基于数字交错技术的生成矩阵形成原理框图。

具体实施方式

实施例1：

以IEEE118节点为例，在节点20、35、38、51、54、67、81、93、97和110处分别接入额定功率为10MW的风电场。风机均以恒定功率因数1运行。假设各个节点的风功率存在相互关系，并且相关系数为0.5。所有负荷服从正态分布，均值为节点负荷的稳态值，标准差为均值的5％，功率因数为0.9。然后进行如下计算：

步骤101：对于负荷，有功功率采用正态分布描述，期望值取原始数据值，标准差为期望的某个百分比，无功功率按照给定的功率因数跟随有功功率变化。

步骤102：对于各个风电场的风功率，采用Weilbull概率分布模型进行描述。

式中，p_i为第i风电场的风功率i＝1,2,…,S；λ_i>0是对应的比例参数；k_i>0是对应的形状参数。

步骤201：P(p₁,p₂,…,p_S)表示维数为S的风功率输入向量；风功率p_i和p_j的相关系数为ρ_ij；f(p_i)和F(p_i)分别为风功率p_i的概率密度函数和累计概率分布函数。Y＝(y₁,y₂,…,y_S)为标准正态分布随机向量，其中y_i和y_j间的相关系数为ρ_yij。根据等概率准则和Nataf变换，可以建立起p_i和y_i对应的转换关系：φ(y_i)＝F(p_i)。

式中：φ表示标准正态随机变量的累积概率分布。

步骤202：采用二重积分法可获得随机向量Y的相关系数矩阵C_Y：

式中：u_pi,u_pj表示为风功率p_i,p_j的均值。

步骤301：给出相关初始值定义：基b为一质数，一般取2；并定义d≥1，m≥1，s≥1为整数，其中d×s代表随机数的最高维数，2^m代表随机数的个数。

定义最高次项为s_j的多项式：

其中j代表该随机数的维数，j≤d×s；系数为0或1。

步骤302：计算“直接数”{v_1,j,v_2,j,…}：

给出下列迭代关系式：

其中⊕表示二进制异或运算；初始值任意给定，只需要满足下列关系：m_k,j,1≤k≤s_j为奇数，并且小于等于2^k。“直接数”{v_1,j,v_2,j,…}计算公式如下：

步骤303：形成矩阵C_m,j，

其中矩阵的元素定义如下：v_k,j＝(0.v_k,j,1v_k,j,2v_k,j,3…)₂。C_m,1,…,C_m,d×s代表d×s个m×m的矩阵

步骤304：生成随机数x_i,j：对于第j维随机数，共有2^m个随机数；给定一整数i,满足1≤i≤2^m，对i进行二进制分解：则对于第j维中的第i个随机数x_i,j，其生成的原则如下：

x_i,j＝(b^-1,b^-2,…,b^-m)·C_m,j(i₁,i₂,…,i_m)^T＝ξ_j,1b^-1+ξ_j,2b^-2+…+ξ_j,mb^-m (4)

步骤305：对于步骤304的d×s维随机数x＝(x₁,x₂,…,x_d×s)∈[0,1)^ds数字交错处理:

经过上式的处理后，可以得到s维随机数D_d,1,D_d,2,…,D_d,(s-1),D_d,s，简记为D_s＝(d₁,d₂,…,d_s)。

步骤401：将步骤202所得到的相关系数矩阵C_Y进行Cholesky分解，得到下三角矩阵L：C_Y＝LL^T。

步骤402：基于步骤3产生的样本矩阵其中d_i∈U(0,1)。由式y_i＝φ^-1(d_i)i＝1,2,…,s得到标准正态分布随机向量Y_s＝(y₁,y₂,…,y_S)。

步骤403：根据Y^corr＝LY求取矩阵Y^corr，此时矩阵Y^corr的协方差矩阵即为C_Y。根据式可以依次转换矩阵Y^corr中的各元素，得到相关随机输入样本P。对每组样本分别进行潮流计算，利用统计学方法得到节点电压和支路潮流的数字特征和概率分布。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (5)

1.一种基于数字交错技术的拟蒙特卡洛随机潮流计算方法，其特征在于：包括以下具体步骤：

步骤1：建立各风电场出力和负荷的概率分布模型；

步骤2：基于等概率转换原则和Nataf变换计算标准正态分布下的相关系数矩阵C_Y；

步骤3：基于数字交错技术的Sobol序列产生分布均匀的S维随机数；由以下具体步骤组成：

本发明对Sobol序列的生成矩阵C_m,j进行处理，定义了相应的函数D_d()，使得高维随机数分布更均匀，有效提高了拟蒙特卡洛法仿真效率。

步骤301：定义最高次项为s_j的多项式：

其中系数为0或1；j代表随机数的维数，j≤d×s；d≥1，s≥1，均为整数，

步骤302：计算直接数{v_1,j,v_2,j,…}：

建立迭代关系式：其中表示二进制异或运算；初始值任意给定，只需要满足：m_k,j,1≤k≤s_j为奇数，并且小于等于2^k；直接数{v_1,j,v_2,j,…}计算公式为：

步骤303：形成矩阵C_m,j，

其中矩阵的元素定义如下：v_k,j＝(0.v_k,j,1v_k,j,2v_k,j,3…)₂。C_m,1,…,C_m,d×s代表d×s个m×m的矩阵；

步骤304：生成随机数x_i,j：对于第j维随机数，共有2^m个随机数；给定一整数i,满足1≤i≤2^m，对i进行二进制分解：则对于第j维中的第i个随机数x_i,j，其生成的原则如下：

x_i,j＝(b^-1,b^-2,…,b^-m)·C_m,j(i₁,i₂,…,i_m)^T＝ξ_j,1b^-1+ξ_j,2b^-2+…+ξ_j,mb^-m

步骤305：对于步骤304的d×s维随机数x＝(x₁,x₂,…,x_d×s)∈[0,1)^ds数字交错处理:

经过上式的处理后，可以得到s维随机数D_d,1,D_d,2,…,D_d,(s-1),D_d,s，简记为D＝(d₁,d₂,…,d_s)；

步骤4：采用Cholesky分解技术获得具有相关性的随机输入变量样本进行随机潮流计算，得到节点电压幅值和支路注入有功和无功。

2.根据权利要求1所述的基于数字交错技术的拟蒙特卡洛随机潮流计算方法，其特征在于：步骤1由以下具体步骤组成：

步骤101：采用正态分布拟合负荷和有功功率的概率密度函数，标准差为期望的5％，无功功率按照给定的功率因数跟随有功功率变化。

步骤102：采用Weibull概率分布模型拟合各个风电场的风功率的概率密度函数：

式中，p_i为第i风电场的风功率i＝1,2,…,S；λ_i>0是对应的比例参数；k_i>0是对应的形状参数。

3.根据权利要求1所述的基于数字交错技术的拟蒙特卡洛随机潮流计算方 法，其特征在于：步骤2由以下具体步骤组成：

步骤201：P(p₁,p₂,…,p_S)表示维数为S的风功率输入向量；风功率p_i和p_j的相关系数为ρ_ij；f(p_i)和F(p_i)分别为风功率p_i的概率密度函数和累计概率分布函数；根据等概率准则和Nataf变换，建立起p_i和y_i对应的转换关系：φ(y_i)＝F(p_i)；φ表示标准正态分布随机变量Y＝(y₁,y₂,…,y_S)的累积概率分布，其中y_i和y_j间的相关系数为ρ_yij；

步骤202：采用二重积分法可获得随机向量Y的相关系数矩阵C_Y：

式中：u_pi,u_pj表示为风功率p_i,p_j的均值。表示具有相关系数ρ_yij的标准正态随机变量y_i和y_j的联合概率密度函数。

4.根据权利要求1所述的基于数字交错技术的拟蒙特卡洛随机潮流计算方法，其特征在于：步骤305也可以通过对生成矩阵C_m,1,…,C_m,d×s的行向量进行变换形成新的生成矩阵C’_m,1,C’_m,2,…,C’_m,d。

5.根据权利要求3所述的基于数字交错技术的拟蒙特卡洛随机潮流计算方法，其特征在于：步骤4由以下具体步骤组成：

步骤401：将步骤202所得到的相关系数矩阵C_Y进行Cholesky分解，得到下三角矩阵L：C_Y＝LL^T；

步骤402：基于步骤3产生的样本矩阵其中d_i∈U(0,1)；由式y_i＝φ^-1(d_i)i＝1,2,…,s得到标准正态分布随机向量Y＝(y₁,y₂,…,y_S)；

步骤403：根据Y^corr＝LY求取矩阵Y^corr，矩阵Y^corr的协方差矩阵即为C_Y；根据式依次转换矩阵Y^corr中的各元素，得到相关随机输入样本P；对每组样本分别进行潮流计算，利用统计学方法得到节点电压和支路潮流的数字特征和概率分布。