CN118153311A

CN118153311A - 一种基于重要抽样函数改进非参数估计的电网可靠性评估方法

Info

Publication number: CN118153311A
Application number: CN202410281462.9A
Authority: CN
Inventors: 赵渊; 蒋林志; 谢开贵; 孙昊; 黄吉星
Original assignee: Chongqing University
Current assignee: Chongqing University
Filing date: 2024-03-12
Publication date: 2024-06-07

Abstract

本发明涉及可靠性评估技术领域，具体涉及一种基于重要抽样函数改进非参数估计的电网可靠性评估方法，包括以下步骤：构建电网系统仿真模型和系统状态向量以原始概率密度函数为初始IS‑PDF，以逼近系统故障状态子空间为目的，迭代的计算最终的系统状态IS‑PDF；直至满足收敛条件，得到最终的系统状态IS‑PDF；采用获取的系统状态IS‑PDF，抽取并分析系统状态样本，计算可靠性指标。该方法在避免重要采样概率密度函数参数估计方法局限性的同时，有效提高大电网可靠性评估的效率和准确性。

Description

一种基于重要抽样函数改进非参数估计的电网可靠性评估方法

技术领域

本发明涉及可靠性评估技术领域，具体涉及一种基于重要抽样函数改进非参数估计的电网可靠性评估方法。

背景技术

蒙特卡罗模拟法对系统中的各种随机因素和复杂运行条件进行随机数模拟，将复杂的可靠性评估问题转化为简单的重复实验，并获得能够反映系统可靠性水平的相关指标，是一种有效且常用的大电网可靠性评估工具。作为一种有效且常用的方差削减技术，重要抽样法通过提高随机变量的抽样效率，有效加速蒙特卡罗模拟收敛。而如何为随机变量选择合适的重要采样概率密度函数(IS-PDF)是关系重要抽样法高效性和可靠性的核心内容。通常将重要采样概率密度函数的求解转化为参数估计问题，即首先确定重要采样概率密度函数的函数形式，再对相关参数进行更新以实现重要采样概率密度函数的参数估计。作为一种有效且广受关注的重要抽样法，交叉熵重要抽样法以待求重要采样概率密度函数与理论最优重要采样概率密度函数之间的Kullback-Leibler(KL)距离最小为目标，通过参数迭代更新实现重要采样概率密度函数的高效获取，已在大电网可靠性评估等领域得到广泛应用。

传统交叉熵重要抽样法，以及相关改进方法，在本质上都是通过参数估计的方式获取随机变量重要采样概率密度函数，因此都需要在参数更新之前确定重要采样概率密度函数的函数形式。作为一种重要采样概率密度函数参数估计方法，交叉熵重要抽样法需要事先确定随机变量重要采样概率密度函数的函数形式，以及需要更新的参数。对于离散变量，例如服从伯努利分布的两状态系统元件，其重要采样概率密度函数通常保持相同的函数形式。而对于系统负荷、风速和光照强度等连续变量，从参数迭代更新公式推导的可行性出发，正态分布、高斯混合模型、多元正态分布、对数正态分布和威布尔分布等指数族分布经常被选作连续变量重要采样概率密度函数的函数形式。尽管上述各种函数形式均在现有技术中被选作连续变量重要采样概率密度函数的函数形式，但现有技术并未提供连续变量重要采样概率密度函数函数形式选择的统一标准。针对具体问题中的某个连续变量，是否存在最佳的重要采样概率密度函数函数形式选择，各种函数形式是否都适用于该变量重要采样概率密度函数的参数估计，并未给出明确回答。

目前，交叉熵重要抽样法通常基于主观经验进行连续变量重要采样概率密度函数的函数形式选择，这将不可避免地给问题求解注入主观因素，实验证明，基于主观经验的连续变量IS-PDF函数形式选择无法保证交叉熵重要抽样法对可靠性评估的加速效果，甚至会影响评估结果的准确性。

另一方面，为了在不限制重要采样概率密度函数函数形式的情况下实现随机变量重要采样概率密度函数的有效获取，现有技术中还提出了一些基于非参数方法的重要采样概率密度函数获取方法，这些方法的原理与重要采样概率密度函数参数估计方法存在本质上的区别。例如现有技术中提出了一种非参数重要抽样法，在不指定随机变量重要采样概率密度函数函数形式的情况下，根据事先给定的初始抽样函数生成随机变量样本，并按各样本在原始概率密度函数和初始抽样函数下的概率密度值之比调整样本权重，将基于这些样本所构建的概率密度函数作为随机变量的重要采样概率密度函数。但该方法存在初始抽样函数的选择问题，并且在处理高维问题时计算效率较低。在此基础上，现有技术中还引入了一系列改进措施，在一定程度上提高了非参数重要抽样法在应对复杂问题时的计算效率，并在航空、金融、林业等领域中得到应用，但并未从根本上解决这两个问题，且无法有效获取离散变量的重要采样概率密度函数。此外，现有技术中还提出了一种非参数分层重要抽样法，将高斯混合模型与拉丁超立方抽样相结合，通过对系统状态重要区域的有效定位，降低重复样本出现的概率，有效提高考虑可再生能源接入的电力系统可靠性评估效率。但该方法只是采用非参数高斯混合模型作为重要采样概率密度函数的函数形式，相关参数的更新仍采用交叉熵方法，因此在本质上还是一种重要采样概率密度函数参数估计方法。

综上，由于电力系统中含有众多随机变量，且发电机、输电线路等离散变量的抽样效率对大电网可靠性评估效率的影响不容忽视，现有非参数方法难以实现大电网可靠性评估中系统状态抽样效率的有效提高。因此，如何实现大电网可靠性评估中随机变量重要采样概率密度函数的非参数估计，在避免重要采样概率密度函数参数估计方法局限性的同时，有效提高大电网可靠性评估的效率和准确性，是一个亟待解决的问题。

发明内容

针对上述现有技术的不足，本发明提供了一种基于重要抽样函数改进非参数估计的电网可靠性评估方法，采用非参数核密度估计获取各连续变量在系统故障状态子空间内的边缘分布，作为大电网可靠性评估中各连续变量的IS-PDF，进而在避免重要采样概率密度函数参数估计方法局限性的同时，有效提高大电网可靠性评估的效率和准确性。

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

一种基于重要抽样函数改进非参数估计的电网可靠性评估方法，包括以下步骤：

步骤1、构建电网系统仿真模型，构建由随机变量构成的系统状态向量x＝[x_d,x_c]，其中x_d是描述D个相互独立的两状态系统元件的离散变量状态向量，状态取值为0和1，分别表示元件故障和元件正常；x_c＝[x_l,x_w]为连续变量状态向量，其中x_l是系统负荷状态变量，x_w是用于描述系统中多个相关能源的多维状态向量，维数分别为1和F；有系统状态x＝[x_d,x_l,x_w]的原始概率密度函数f(x|Ω)的表达式如下：

式中，Ω是系统状态空间，f(x_d|Ω)、f(x_l|Ω)和f(x_w|Ω)＝f(x_w,1,…,x_w,F|Ω)分别是系统元件、系统负荷和多个相关能源的原始概率密度函数；

步骤2、以原始概率密度函数为初始IS-PDF，即以逼近系统故障状态子空间为目的，迭代的计算最终的系统状态IS-PDF；

每次迭代中，基于前一次所得/初始IS-PDF生成多个系统状态样本x_i，再选取相应系统性能度量函数值S(x_i)较小的指定数量个样本，按照多维核密度估计构建系统负荷和多个相关能源的概率密度函数，并利用概率密度函数似然比对所构建的概率密度函数进行权重校正，以得到新的概率密度函数作为当前IS-PDF，其中Ω_k为第k次迭代所得到的系统状态子空间，直至满足收敛条件，将此时的/>作为最终的系统状态IS-PDF；

步骤3、采用获取的系统状态IS-PDF，抽取并分析系统状态样本，计算可靠性指标。

进一步的，步骤2中包括：

步骤21、设定比例系数ρ，初始化迭代计数器k＝1，令

步骤22、基于概率密度函数生成N个系统状态样本x_i，1≤i≤N，并进行系统状态分析，将各系统状态的系统性能度量函数值S(x_i)按从小到大的顺序排列，有：S_[1]≤S_[2]≤…≤S_[N]，并将对应系统状态样本按相同顺序排列为x_[j]，1≤j≤N；

步骤23、若S_[ρN]＞0，令N_u＝ρN；否则，令/>N_u＝N_f，N_f表示N个系统状态样本中所含的故障样本数；

步骤24、令系统状态样本x_[j]，1≤j≤N_u，确定的系统状态子空间为Ω_k，从x_[j]中提取系统负荷样本x_lj和多维能源样本x_wj；

若k＝1，分别采用以下方式估计系统负荷和多维能源在Ω_k内的边缘分布；

式中，h和H分别是核密度估计的带宽系数和带宽矩阵，以渐进积分均方误差最小为目标进行求解，K(·)是高斯核函数，有：

否则，采用以下方式对f(x_l|Ω_k)和f(x_w|Ω_k)进行无偏估计：

其中，W(x_[j])为系统状态x_[j]的似然比，有：

式中，为指示函数，若x_[j]∈Ω_k，/>若/>p(x_[j]∈Ω_k)是系统状态x_[j]∈Ω_k的概率；

步骤25、令若满足收敛条件：/>得到/>作为系统状态IS-PDF，并结束迭代过程；否则，令k＝k+1，返回抽样步骤22进行下一次IS-PDF迭代估计。

进一步的，步骤1中还包括，选定一个定义域为(-∞,+∞)的概率密度函数h(·)，并得到对应的累积分布函数H(·)；改用分布区间为(-∞,+∞)的连续元件状态变量x_d,i原件状态，并将h(·)作为第i，1≤i≤D，个系统元件状态的原始概率密度函数，即f(x_d,i|Ω)＝h(x_d,i)，确定第i个系统元件正常状态与故障状态的分界点ε_i：

ε_i＝H^-1(u_i) (34)

式中，H^-1(·)是累积分布函数H(·)的反函数，u_i为第i个系统元件可用度；

步骤21中，初始化

步骤22中，第i个系统元件的状态抽样方法为：基于f(x_d,i|Ω)抽取(-∞,+∞)上的随机数x_d,i，若x_d,i≤ε_i，则第i个系统元件处于故障状态，否则该元件处于正常状态；

步骤24中，还从x_[j](1≤j≤N_u)中提取第i个系统元件的状态样本x_dj,i，1≤j≤N_u，若k＝1，采用下式估计Ω_k内第i个系统元件的边缘分布f(x_d,i|Ω_k)；

否则，采用下式估计f(x_d,i|Ω_k)：

步骤25中，令还包括，若满足收敛条件：/>采用下式计算系统元件的新不可用率v_i，1≤i≤D：

步骤3中，采用步骤22中的方式对第i个系统元件的状态进行抽样方法，并且，采用新不可用率v_i计算第i个系统元件正常状态与故障状态的分界点ε_i。

优选的，概率密度函数h(·)为标准正态分布或对数正态分布概率密度函数。

优选的，步骤24中，所述收敛条件为：或者从第k-1次迭代到第k次迭代时，/>的变化量小于20％。

优选的，设定比例系数ρ＝0.1。

进一步的，步骤24中，引入归一化因子NF，并采用如下所示的W'(x_[j])代替式(11)和式(12)中的系统状态似然比W(x_[j])：

进一步的，所述步骤3包括:

步骤31、设置初始化参数：方差系数β_max，样本计数器N_MC＝t＝0；

步骤32、令t＝t+1，N_MC＝t，基于系统状态IS-PDF g(x)生成第t个系统状态x_t，计算似然比W(x_t)＝f(x_t)/g(x_t)；

步骤33、通过系统状态分析，计算系统状态x_t的负荷削减量LC(x_t)，若LC(x_t)>0，指示函数H(x_t)＝1，否则H(x_t)＝0；

步骤34、根据式(22)和式(23)计算可靠性指标：失负荷概率LOLP和期望缺供电量EENS；

步骤35、计算可靠性指标EENS的方差系数β_EENS，若β_EENS<β_max，结束迭代过程，并输出可靠性指标LOLP和EENS的计算结果；否则，返回步骤32继续计算可靠性指标。

优选的，设置方差系数β_max＝0.2。

本发明在大电网可靠性评估中提出了一种连续变量IS-PDF非参数估计方法，采用非参数核密度估计获取各连续变量在系统故障状态子空间内的边缘分布，作为大电网可靠性评估中各连续变量的IS-PDF。在此基础上，提出基于似然比的IS-PDF迭代估计，通过对系统故障状态子空间的逐步快速逼近，降低高可靠性系统中连续变量IS-PDF的获取成本。构建基于连续变量IS-PDF非参数估计的大电网可靠性评估方法，有效提高可靠性评估的效率和准确性。

在进一步的实施例中，提出一种离散变量IS-PDF非参数估计方法，实现系统元件IS-PDF的非参数估计。分析核密度估计的边界效应对系统元件IS-PDF迭代估计的影响，提出一种离散变量区间映射方法，避免边界效应影响系统元件状态抽样效率。最后将连续变量和离散变量的IS-PDF非参数估计方法相结合，构建基于IS-PDF非参数估计的大电网可靠性评估方法，实现可靠性评估效率的进一步提高。本发明将离散变量IS-PDF非参数估计(NEIS-D)方法和连续变量IS-PDF非参数估计(NEIS-C)方法相结合，构成IS-PDF非参数估计(NEIS)方法，进行大电网可靠性评估中离散变量和连续变量的IS-PDF联合估计。

附图说明

为了使发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的详细描述，其中：

图1为本发明实施例中对系统故障状态子空间逐步逼近的示意图。

图2为本发明实施例中预抽样阶段的流程图。

图3为本发明实施例中第i个系统元件在连续化处理后的概率密度函数。

图4为本发明实施例中示例性的边界效应可能导致的的分布情况的示意图。

图5为本发明实施例中第i个系统元件在区间映射后的概率密度函数。

图6为本发明实施例中基于IS-PDF非参数估计的大电网可靠性评估方法的流程图。

图7为本发明实施例中的方法应用于IEEE-RTS96测试系统所得到的系统元件不可用率的更新结果示意图。

图8为本发明实施例中的方法应用于IEEE-RTS96测试系统所得到的系统负荷原始PDF和不同方法得到的系统负荷IS-PDF的对比图。

图9为本发明实施例中的方法应用于含多个风电场的IEEE-RTS96测试系统所得到的系统元件不可用率的更新结果示意图。

图10为本发明实施例中的方法应用于含多个风电场的IEEE-RTS96测试系统所得到的系统负荷原始PDF和不同方法得到的系统负荷IS-PDF。

图11为本发明实施例中的方法应用于含多个风电场的IEEE-RTS96测试系统所得到的节点13、节点37和节点61的风速原始PDF和IS-PDF的对比图。

具体实施方式

下面通过具体实施方式进一步详细的说明：

假设x是由n_v个随机变量构成的系统状态向量，Ω是系统状态空间，f(x|Ω)是系统状态的原始概率密度函数，系统故障概率p_g可按式(1)进行计算：

式中，x_i(1≤i≤N)是基于原始概率密度函数f(x|Ω)抽取的系统状态样本；H(x_i)为指示函数，若x_i为故障状态，H(x_i)＝1，否则H(x_i)＝0。

重要抽样法采用IS-PDF g(x)代替f(x|Ω)进行系统状态样本的高效抽取，此时可按(2)进行系统故障概率p_g的无偏估计：

式中，x_i(1≤i≤N)是基于IS-PDF g(x)抽取的系统状态样本；W(x_i)＝f(x_i|Ω)/g(x_i)是各系统状态样本对应的似然比，用于确保式(2)是对p_g的无偏估计。

由式(1)和式(2)可知，在所有系统状态样本x_i(1≤i≤N)中，只有故障样本对可靠性指标p_g的计算是有贡献的，因此系统故障状态子空间Ω_f就是系统状态x的重要区域。如果能够获取系统状态x在系统故障状态子空间Ω_f内的联合概率密度函数f(x|Ω_f)，并在可靠性评估中代替原始概率密度函数f(x|Ω)，则系统状态抽样效率将得到有效提高。因此，如式(3)所示的条件概率密度函数f(x|Ω_f)可作为系统状态IS-PDF g(x)：

式中，是指示函数，若系统状态x属于系统故障状态子空间Ω_f，/>否则/>p(x∈Ω_f)是系统状态x属于Ω_f的概率，即系统故障概率p_g。

根据上述思路，本例中提出一种系统状态IS-PDF的非参数估计方法：

基于原始概率密度函数f(x|Ω)生成N'个系统状态样本，并进行系统状态分析，假设其中有N_f个故障样本x_j(1≤j≤N_f)，采用非参数多维核密度估计构建概率密度函数：

式中，H是d×d维对称正定带宽矩阵，以渐进积分均方误差最小为目标进行求解；n_v是系统状态x的维数，即系统中的随机变量数。

当抽取的系统状态样本足够多，即N'→∞时，故障样本x_j(1≤j≤N_f)能够全面且相对均匀地覆盖系统状态x的重要区域，式(4)就成为对系统故障状态子空间Ω_f内系统状态x的实际概率密度函数f(x|Ω_f)的无偏估计。

从连续变量IS-PDF的非参数估计着手，通过分别估计各连续变量在系统故障状态子空间Ω_f内的边缘分布，即可实现大电网可靠性评估中连续变量IS-PDF的有效获取。

假设系统状态x＝[x_d,x_c]，其中x_d是离散变量的状态向量，x_c是连续变量的状态向量。假设x_c＝[x_l,x_w]，其中x_l是系统负荷状态变量，x_w是系统中多个相关风电场的风速状态向量，维数分别为1和F。则系统状态x＝[x_d,x_l,x_w]的原始概率密度函数f(x|Ω)的表达式如下：

式中，x_w＝[x_w,1,…,x_w,F]是F个相关风电场的风速状态向量；f(x_d|Ω)、f(x_l|Ω)和f(x_w|Ω)＝f(x_w,1,…,x_w,F|Ω)分别是系统元件、系统负荷和多维风速的原始概率密度函数。

基于系统状态x的原始概率密度函数f(x|Ω)生成N'个系统状态样本，找出所有故障样本x_j(1≤j≤N_f)，并从中提取系统负荷样本x_lj(1≤j≤N_f)和多维风速样本x_wj＝[x_wj,1,…,x_wj,F]^T(1≤j≤N_f)。分别采用传统核密度估计和如式(4)所示的多维核密度估计构建概率密度函数：

式(6)中的h和式(7)中的H分别是核密度估计的带宽系数和带宽矩阵，以渐进积分均方误差最小为目标进行求解；K(·)是高斯核函数，表达式如下：

当抽取的系统状态样本足够多，即N'→∞时，从故障样本中提取的系统负荷和多维风速样本，将能够对各自的重要区域实现全面且相对均匀的覆盖，式(6)和式(7)就成为对系统故障状态子空间Ω_f内系统负荷边缘分布f(x_l|Ω_f)，和多维风速边缘分布f(x_w|Ω_f)的无偏估计。此时，将和/>分别作为系统负荷IS-PDF g(x_l)和多维风速IS-PDFg(x_w)，可获得如式(9)和式(10)所示的连续变量联合IS-PDF g(x_c)，以及系统状态IS-PDF g(x)：

上述连续变量IS-PDF非参数估计方法与交叉熵重要抽样法存在本质区别。交叉熵重要抽样法需要基于主观经验选择连续变量IS-PDF的函数形式，并将所选函数形式作为故障样本信息的载体，以参数迭代更新的方式逼近理论最优IS-PDF。本例提出的连续变量IS-PDF非参数估计方法无须限制IS-PDF的函数形式，而是从故障样本中提取各连续变量样本，对系统故障状态子空间Ω_f内各连续变量的边缘分布进行直接估计，作为各连续变量的IS-PDF，在避免IS-PDF函数形式选择带来主观因素影响的同时，本例所提方法能够获得各连续变量在其重要区域上更加真实的分布情况，有利于变量抽样效率的有效提高。

由式(6)和式(7)可知，只有故障样本参与连续变量IS-PDF的非参数估计，当抽取的系统状态样本数量有限时，故障样本x_j(1≤j≤N_f)可能无法全面且相对均匀地覆盖系统故障状态子空间Ω_f，导致由式(6)和式(7)得到的概率密度函数和/>无法反映Ω_f内系统负荷和多维风速的实际分布情况。为准确估计系统故障状态子空间Ω_f内各连续变量的边缘分布，需要抽取充足的系统状态样本，但在系统故障状态子空间稀疏且分散的高可靠性系统中，故障样本通常难以抽取。为避免抽取并分析大量系统状态样本带来的高计算成本，本例采用一种对系统故障状态子空间Ω_f逐步逼近的思路，提高连续变量IS-PDF的非参数估计效率。

图1中，Ω₀＝Ω是原始系统状态空间，Ω₁是包含系统故障状态子空间Ω_f的系统状态子空间。对Ω₁进一步缩小，可得到包含Ω_f的系统状态子空间Ω₂。这种将Ω₀不断缩小至Ω_f的过程，也是对各连续变量重要区域不断逼近的过程。本发明将这种思路运用于连续变量IS-PDF的非参数估计，提出了一种基于似然比的IS-PDF迭代估计算法：每次迭代中抽取一定数量的系统状态样本，允许故障和接近故障的系统状态样本共同参与连续变量IS-PDF的非参数估计；每次迭代获取的连续变量IS-PDF被用于下一次迭代中系统状态样本的高效抽取；通过多次迭代，对系统故障状态子空间Ω_f进行逐步快速逼近。

假设系统状态x＝[x_d,x_c]，其中x_d是系统元件的状态向量(离散变量)，x_c是连续变量的状态向量。假设x_c＝[x_l,x_w]，其中x_l表示系统负荷状态变量，x_w是多个相关风电场的风速状态向量，维数分别为1和F。引入系统性能度量函数S(x)(Geng L,Zhao Y,LiW.Enhanced cross entropy method for composite power system reliabilityevaluation[J].IEEE Transactions on Power Systems,2019,34(4):3129-3139.),比例系数ρ和非负阈值系统负荷和多维风速的IS-PDF迭代估计步骤如下：

①设置初始化参数：每次迭代需要抽取的系统状态样本数N，比例系数ρ，迭代计数器k＝1。系统元件概率密度函数f(x_d)，系统负荷初始概率密度函数f(x_l|Ω)，多维风速初始概率密度函数f(x_w|Ω)，构成系统状态初始概率密度函数

②基于抽取N个系统状态样本x_i(1≤i≤N)，并进行系统状态分析，将各系统状态的系统性能度量函数值S(x_i)(1≤i≤N)按从小到大的顺序排列：S_[1]≤S_[2]≤…≤S_[N]，并将对应系统状态样本按相同顺序排列：x_[j](1≤j≤N)；

③若S_[ρN]＞0，令N_u＝ρN；

否则，令N_u＝N_f(故障样本数)；

④令系统状态样本x_[j](1≤j≤N_u)确定的系统状态子空间为Ω_k，从x_[j](1≤j≤N_u)中提取系统负荷样本x_lj(1≤j≤N_u)和多维风速样本x_wj(1≤j≤N_u)，若k＝1，分别采用式(6)和式(7)估计系统负荷和多维风速在Ω_k内的边缘分布；

否则，采用式(11)和式(12)对f(x_l|Ω_k)和f(x_w|Ω_k)进行无偏估计：

式中，K(·)是高斯核函数，W(x_[j])是如式(13)所示的系统状态x_[j]的似然比：

式中，为指示函数，若x_[j]∈Ω_k，/>若/>

⑤令若/>分别将/>和/>作为系统负荷和多维风速的IS-PDF，即/>得到系统状态IS-PDF/>并结束迭代过程；

否则，令k＝k+1，返回步骤②进行下一次IS-PDF迭代估计。

在第k(k≥2)次迭代中，式(11)和式(12)采用系统状态似然比W(x_[j])对系统负荷样本x_lj(1≤j≤N_u)和多维风速样本x_wj(1≤j≤N_u)进行权重校正，这是因为：

(2)N个系统状态样本x_i(1≤i≤N)中，哪些样本参与连续变量IS-PDF的非参数估计是由系统性能度量函数值S(x_i)(1≤i≤N)决定的。以系统负荷IS-PDF的迭代估计为例，参与估计f(x_l|Ω_k)的系统负荷样本x_lj(1≤j≤N_u)的出现概率不仅与边缘分布有关，同时也受到系统中其他变量分布情况的影响。

假设离散变量x有K种取值{X₁,…X_K}，离散变量y有M种取值{Y₁,…Y_M}，向量(x,y)服从联合分布P₁(x,y)，且离散变量x的边缘分布可以表示为：

式中，δ(·)是狄拉克函数，即单位脉冲函数；P(X_i)是x＝X_i(1≤i≤K)的概率。

需要解决的问题是：假设f(x)未知，如何根据由另一个联合分布P₂(x,y)生成的N_xy个样本(x_j,y_j)(1≤j≤N_xy)，获取对随机变量x的原始边缘分布f(x)的无偏估计。

假设样本(x_j,y_j)(1≤j≤N_xy)中，x＝X_i(1≤i≤K)出现的次数为N_i，(x_i,y_m)(1≤m≤M)出现的次数为则有：

从样本(x_j,y_j)(1≤j≤N_xy)中提取变量x的样本x_j(1≤j≤N_xy)，采用非参数核密度估计构建概率密度函数：

为了获取随机变量x的原始边缘分布f(x)，采用R(x_j,y_j)＝P₁(x_j,y_j)/P₂(x_j,y_j)，即联合分布P₁(x,y)和P₂(x,y)下样本(x_j,y_j)的概率之比，对样本x_j(1≤j≤N_xy)的权重进行校正，则式(16)变为：

由式(17)得到的就是对式(14)所示边缘分布f(x)的无偏估计，证明如下：

式(17)可进行如下变形：

当基于联合分布P₂(x,y)抽取的向量(x,y)样本足够多，即N_xy→∞时，式(18)中的将收敛于P₂(X_i,Y_m)。因此，式(18)可以改写为：

/>

此时，式(19)中的K_h(x-X_i)和将分别收敛于狄拉克函数δ(x-X_i)和P(X_i)。因此，式(19)就是对式(14)所示边缘分布f(x)的无偏估计，证明完毕。

当M→∞时，离散变量x转变为连续变量。将上述边缘分布的无偏估计方法用于大电网可靠性评估中连续变量IS-PDF的非参数迭代估计时，式(17)中的概率之比R(x_j,y_j)将变为系统状态的似然比。因此，式(11)和式(12)采用系统状态似然比W(x_[j])对系统负荷样本x_lj(1≤j≤N_u)和多维风速样本x_wj(1≤j≤N_u)进行权重校正，保证所得和/>是对目标边缘分布f(x_l|Ω_k)和f(x_w|Ω_k)的无偏估计。

考虑实际仿真过程中抽取的系统状态样本数量有限，引入如式(20)所示的归一化因子NF，并采用如式(21)所示的W'(x_[j])代替式(11)和式(12)中的系统状态似然比W(x_[j])，确保式(11)和式(12)中所有核函数的权重之和为1。得益于NF的引入，式(13)中的p(x_[j]∈Ω_k)无须计算。

上述基于似然比的IS-PDF迭代估计算法能够对系统故障状态子空间Ω_f进行逐步快速逼近，并根据实际样本获取各连续变量在其重要区域上的分布情况，实现大电网可靠性评估中连续变量IS-PDF的高效获取。

基于连续变量IS-PDF非参数估计的大电网可靠性评估由两阶段构成：(1)预抽样阶段：采用基于似然比的IS-PDF迭代估计算法，进行系统中各连续变量IS-PDF的非参数估计，并构成系统状态IS-PDF；(2)主抽样阶段，采用预抽样阶段获取的系统状态IS-PDF，抽取并分析系统状态样本，计算可靠性指标。

预抽样阶段：

预抽样阶段流程如图2所示，其中，进行系统负荷和多维风速的IS-PDF非参数估计，考虑实际仿真过程中迭代终止条件可能无法满足，新增迭代终止条件：若从第k-1次迭代到第k次迭代时，/>的变化量小于20％，迭代结束。

主抽样阶段：

①设置初始化参数：方差系数β_max，样本计数器N_MC＝t＝0；

②令t＝t+1，N_MC＝t，基于系统状态IS-PDF g(x)生成第t个系统状态x_t，计算似然比W(x_t)＝f(x_t)/g(x_t)；

③通过系统状态分析，计算系统状态x_t的负荷削减量LC(x_t)，若LC(x_t)>0，指示函数H(x_t)＝1，否则H(x_t)＝0；

④根据式(22)和式(23)计算可靠性指标：失负荷概率(LOLP)和期望缺供电量(EENS)；

/>

⑤计算可靠性指标EENS的方差系数β_EENS，若β_EENS<β_max，结束迭代过程，并输出可靠性指标LOLP和EENS的计算结果；否则，返回步骤②继续计算可靠性指标。

在大电网可靠性评估中，系统状态分析需要对发电机、输电线路等相互独立的两状态系统元件进行状态抽样。假设系统中有D个两状态系统元件，则第i(1≤i≤D)个系统元件的状态x_d,i可采用系统元件状态抽样方法进行确定：

式中，Randn_i是一个随机生成的(0,1)上的随机数，服从均匀分布U(0,1)；u_i是第i个系统元件的不可用率；x_d,i＝0代表第i个系统元件处于正常状态，x_d,i＝1表示该元件处于故障状态。

由于第i个系统元件的状态变量x_d,i服从如式(25)所示参数为u_i的伯努利分布，且只有x_d,i＝0和x_d,i＝1两种状态，前面提出的连续变量IS-PDF非参数估计方法无法用于获取系统元件的IS-PDF。

因此，本发明首先提出了一种离散变量的连续化方法，将系统元件转变为连续变量，为系统元件IS-PDF的非参数估计创造条件。以第i个系统元件为例，离散变量连续化方法的基本原理如下：

将第i个系统元件的原始概率密度函数转变成区间(0,1)上的均匀分布：

此时，x_d,i转变为服从均匀分布U(0,1)的连续变量，第i个系统元件的概率密度函数f(x_d,i|Ω)如图3所示，元件不可用率u_i将区间(0,1)划分为(0,u_i]和(u_i,1)两个区间，f(x_d,i|Ω)的函数图像也被x_d,i＝u_i分为区域A和区域B两个部分。

连续化处理之后，第i个系统元件的概率密度函数f(x_d,i|Ω)满足：

区域A的面积就是第i个系统元件的不可用率，则元件状态抽样方法变为：基于f(x_d,i|Ω)抽取(0,1)上的随机数x_d,i，若x_d,i≤u_i，则第i个系统元件处于故障状态，否则第i个系统元件处于正常状态。根据新的状态抽样方法确定系统元件状态，并未改变系统元件的故障概率。

假设第i个系统元件由正常状态转入故障状态时，整个系统的安全稳定运行会受到较大冲击，则将第i个系统元件称为关系系统可靠性水平的关键元件。根据连续化处理之后的系统元件状态抽样方法，当随机数x_d,i落入区间(0,u_i]时，第i个系统元件处于故障状态。因此，若第i个系统元件是关键元件，则区间(0,u_i]就是该元件的重要区域。如果能在保持x_d,i＝u_i位置不变的同时，将图3中区域A的面积进行合理的扩大，并同步缩小区域B的面积以满足概率密度函数在定义域上积分为1的性质，即：

式中，s(A)和s(B)分别表示区域A和区域B的面积，即f(x_d,i|Ω)在区间(0,u_i]和区间(u_i,1)上的积分。

则第i个系统元件的等效不可用率，即区域A的面积将得到提高，在仿真过程中将以更大的概率抽到该元件的故障状态样本。如果能够对系统中所有关键元件的重要程度进行识别，并根据各元件对系统可靠性的重要程度合理扩大其区域A的面积，就能够通过提高系统元件状态抽样效率，对大电网可靠性评估进行有效加速。

假设系统状态x＝[x_d,x_c]，其中x_d是D个相互独立系统元件的状态向量，x_c是系统中所有连续变量的状态向量，假设Ω是原始系统状态空间，f(x|Ω)＝f(x_d|Ω)f(x_c|Ω)是系统状态x的原始概率密度函数，其中f(x_d|Ω)的表达式如下：

式中，f(x_d,i|Ω)(1≤i≤D)是第i个系统元件的原始概率密度函数。在进行离散变量连续化处理之后，f(x_d,i|Ω)成为(0,1)上的均匀分布。此时，向量x_d中的元素不再是对应元件正常或故障状态的数值0或1，而是与元件状态相关的(0,1)上的随机数。以第i个系统元件为例，离散变量IS-PDF的非参数估计原理如下：

基于原始概率密度函数f(x|Ω)生成N'个系统状态样本，并进行系统状态分析。假设其中有N_f个故障样本x_j(1≤j≤N_f)，从故障样本中提取用于确定第i个系统元件状态的随机数x_dj,i(1≤j≤N_f)，采用非参数核密度估计构建概率密度函数：

式中，h是核密度估计的带宽系数，以渐进积分均方误差最小为目标进行求解；Ω_f表示系统故障状态子空间；K(·)是如式(8)所示的高斯核函数。

假设第i个系统元件是关系系统可靠性水平的关键元件，则在故障样本x_j(1≤j≤N_f)中该元件处于故障状态的概率较高，即元件状态样本x_dj,i(1≤j≤N_f)将以更大的概率落入区间(0,u_i]。因此，在保持x_d,i＝u_i位置不变的情况下，中区域A的面积相比原始概率密度函数f(x_d,i|Ω)将得到一定程度的扩大。当抽取的系统状态样本足够多，即N'→∞时，由式(30)构建的概率密度函数/>将以区域A面积的变化情况反映第i个系统元件的重要程度，/>可作为可靠性评估中第i个系统元件的IS-PDF g(x_d,i)，并得到如式(31)所示第i个系统元件的新不可用率v_i：

然而，采用式(31)得到的新不可用率v_i可能无法准确反映第i个系统元件的重要程度。这是因为当采用非参数核密度估计进行有界区间概率密度估计时，所得概率密度函数的左、右尾部通常会超出变量的原始分布区间，将这种现象称为非参数核密度估计的边界效应。尽管参与估计第i个系统元件IS-PDF的样本x_dj,i(1≤j≤N_f)均为(0,1)上的随机数，但在边界效应的作用下，由式(30)得到的概率密度函数的左、右尾部将分别穿过x_d,i＝0和x_d,i＝1，导致/>在区间(0,1)上的积分小于1，且有如下关系式成立：/>

如果仍以区域A的面积s(A)作为第i个系统元件的新不可用率v_i，将造成对该元件重要程度的低估，并导致元件状态抽样效率偏低。由于x_d,i的分布区间为(0,1)，如式(33)所示的不可用率更符合实际仿真中抽到该元件故障状态样本的概率。尽管采用式(33)得到的新不可用率能够体现出元件状态抽样效率的提高，但无法体现实际样本所反映的该元件对系统可靠性的重要程度。

为提升高可靠性系统中系统元件IS-PDF的非参数估计效率，可引入上面提出的基于似然比的IS-PDF迭代估计算法，实现系统元件IS-PDF的高效获取，避免抽取并分析大量系统状态样本造成的高计算成本。

采用基于似然比的IS-PDF迭代估计算法时，核密度估计的边界效应对系统元件状态抽样效率的影响不容忽视，具体分析如下：

假设第i(1≤i≤D)个系统元件是关系系统可靠性水平的关键元件，其原始不可用率为u_i。在基于似然比的IS-PDF迭代估计中，假设第k次迭代获取的该元件状态变量的边缘分布为由于区间(0,u_i]是第i个系统元件的重要区域，参与估计概率密度函数的样本将大概率落入区间(0,u_i]。若不可用率u_i较小，则区间(0,u_i]在整个分布区间(0,1)中的占比较小，当多数样本落入狭窄区间(0,u_i]时，非参数核密度估计的边界效应可能使/>呈现出类似图4的分布情况。

如图4所示，的函数曲线有部分穿过x_d,i＝0和x_d,i＝1，将/>在x_d,i＝0左侧的部分称为区域C，其面积用s(C)表示。尽管s(A)/[s(A)+s(B)]>u_i，但落入区间(0,u_i]的样本对/>的贡献更多地体现为s(C)的增大，导致s(A)的变化情况无法体现实际样本所反映的第i个系统元件的重要程度。由于x_d,i的分布区间为(0,1)，s(C)的增长对提高下一次迭代中第i个系统元件状态抽样效率的贡献为0。

区域C的面积反映的就是非参数核密度估计边界效应的严重程度，且边界效应对系统元件状态抽样效率的影响将在迭代过程中不断叠加，导致有限迭代次数下得到的第i个系统元件的新不可用率v_i偏小。如果不能有效遏制核密度估计边界效应对系统元件IS-PDF迭代估计的影响，将无法有效提高系统元件状态抽样效率，不利于大电网可靠性评估效率的进一步提高。

为避免非参数核密度估计的边界效应对系统元件状态抽样效率的影响，本发明提出了一种离散变量区间映射方法：在对系统元件进行连续化处理之后，以某个定义域为(-∞,+∞)的概率密度函数为媒介，将系统元件状态变量的分布区间由(0,1)映射至(-∞,+∞)。以第i(1≤i≤D)个系统元件为例，介绍离散变量的区间映射原理：

假设h(·)是一个定义域为(-∞,+∞)的概率密度函数，H(·)是对应的累积分布函数。将h(·)作为第i个系统元件的原始概率密度函数，即f(x_d,i|Ω)＝h(x_d,i)，则元件状态变量x_d,i的分布区间将由(0,1)变为(-∞,+∞)。此时，需要确定如式(34)所示第i个系统元件正常状态与故障状态的分界点ε_i：

ε_i＝H^-1(u_i) (34)

式中，H^-1(·)是累积分布函数H(·)的反函数，由式(34)得到的分界点ε_i满足：

此时，第i个系统元件的状态抽样方法变为：基于f(x_d,i|Ω)抽取(-∞,+∞)上的随机数x_d,i，若x_d,i≤ε_i，则第i个系统元件处于故障状态，否则该元件处于正常状态。

函数h(·)的选择是多样的，例如正态分布、对数正态分布等，本例选择标准正态分布作为各系统元件的原始概率密度函数。如图5所示，以x_d,i＝ε_i为分界线，标准正态分布形式的概率密度函数f(x_d,i|Ω)被分成区域A和区域B两个部分，区域A的面积就是第i个系统元件的原始不可用率u_i。

假设在系统元件IS-PDF迭代估计中，由第k次迭代得到第i个系统元件的IS-PDF则可以采用式(36)计算第i个系统元件的新不可用率v_i：

以上就是离散变量区间映射方法的基本原理。经过离散变量连续化和区间映射处理，系统元件状态变量的分布区间将由(0,1)转变为(-∞,+∞)。每次迭代中的系统元件状态样本抽取，以及如式(36)所示的系统元件新不可用率的计算，将不再受到样本分布区间的限制。因此，非参数核密度估计的边界效应将不再影响系统元件IS-PDF的迭代估计，系统元件状态抽样效率得到保证。

将所提出的连续变量IS-PDF非参数估计方法和上述离散变量IS-PDF非参数估计方法相结合，进行系统状态IS-PDF的非参数估计：每次迭代从同一组系统状态样本中提取各随机变量样本，进行离散变量和连续变量的IS-PDF联合估计，将各变量边缘分布相乘，构成系统状态IS-PDF。这种系统状态IS-PDF的非参数估计方法不仅可以避免系统状态样本的重复抽取，并且通过同时对离散变量和连续变量进行重要抽样，进一步提高了系统状态抽样效率，有效减少IS-PDF迭代估计过程中需要抽取并分析的系统状态样本数。在此基础上，构建了基于IS-PDF非参数估计的大电网可靠性评估方法。

基于IS-PDF非参数估计的大电网可靠性评估由两阶段构成：(1)预抽样阶段：采用基于似然比的IS-PDF迭代估计算法，进行离散变量和连续变量的IS-PDF联合估计，并最终构成系统状态IS-PDF；(2)主抽样阶段：基于预抽样阶段获取的系统状态IS-PDF生成系统状态样本，进行系统状态分析和可靠性指标计算，算法流程如图6所示。

假设系统状态x＝[x_d,x_l,x_w]，其中x_d是系统元件状态向量，x_l是系统负荷状态变量，x_w是多个相关风电场的风速状态向量，维数分别为D，1和F，可靠性评估流程如下：

预抽样阶段：

①设置初始化参数：每次迭代抽取的系统状态样本数N，比例系数ρ，迭代计数器k＝1。选择标准正态分布作为系统元件的初始概率密度函数f(x_d,i|Ω)(1≤i≤D)，并采用式(34)确定系统元件正常状态与故障状态的分界点ε_i(1≤i≤D)。系统负荷初始概率密度函数f(x_l|Ω)，多维风速初始概率密度函数f(x_w|Ω)，构成系统状态初始概率密度函数

②基于概率密度函数抽取N个系统状态样本x_n(1≤n≤N)，并进行系统状态分析。将各系统状态的系统性能度量函数值S(x_i)(1≤i≤N)按从小到大的顺序排列：S_[1]≤S_[2]≤…≤S_[N]，并将对应的系统状态样本按相同顺序排列：x_[j](1≤j≤N)；

③若S_[ρN]＞0，令N_u＝ρN；否则，令/>N_u＝N_f(故障样本数)；

⑤令若/>将/>作为系统状态IS-PDF，即/>采用式(36)计算系统元件的新不可用率v_i(1≤i≤D)，并结束迭代过程；否则，令k＝k+1，返回步骤②进行下一次IS-PDF迭代估计。

主抽样阶段：

①设置初始化参数：方差系数β_max，样本计数器N_MC＝t＝0；

③通过系统状态分析，计算负荷削减量LC(x_t)，若LC(x_t)>0，指示函数H(x_t)＝1，否则H(x_t)＝0；

④计算可靠性指标LOLP和EENS；

所提出的离散变量IS-PDF非参数估计(NEIS-D)方法和连续变量IS-PDF非参数估计(NEIS-C)方法相结合，构成IS-PDF非参数估计(NEIS)方法，进行大电网可靠性评估中离散变量和连续变量的IS-PDF联合估计。为验证NEIS法，以及所提出的离散变量区间映射方法的有效性，选择IEEE-RTS96测试系统和含有多个相关风电场的IEEE-RTS96测试系统进行可靠性评估。

为全面对比传统交叉熵重要抽样法(CEIS)、NEIS-C法，以及NEIS法对大电网可靠性评估的加速效果，NEIS法的参数设置为：每次迭代抽取的系统状态样本数N＝20000，比例系数ρ＝10％。与可靠性指标LOLP相比，可靠性指标EENS的收敛速度较慢，因此以EENS的方差系数β_EENS小于给定值β_max＝0.02作为计算可靠性指标时非序贯蒙特卡罗模拟的收敛判据。选择文献Subcommittee P M.IEEEreliability test system[J].IEEE Transactionson power apparatus and systems,1979(6):2047-2054.中的年时序负荷曲线，采用传统核密度估计构建系统负荷的原始概率密度函数。选择文献National MeteorologicalInformation Center.The Hourly Wind Speed Data for Each Site in the Period of10Years.[Online].Available:http://data.cma.cn/,Accessed on:Aug.2012.中的风速数据，采用多维核密度估计构建多维风速的原始概率密度函数。

IEEE-RTS96测试系统：

为检验所提出的NEIS法在应对大规模电力系统分析时的有效性，采用NEIS法和不进行离散变量区间映射的NEIS法对IEEE-RTS96测试系统进行可靠性评估，并与原始非序贯蒙特卡罗模拟法(CNMCS)、CEIS法和NEIS-C法的可靠性评估结果进行对比。采用两种NEIS法进行系统元件和系统负荷的IS-PDF联合估计，并采用式(33)计算系统元件的新不可用率。图7选择不可用率变化较大的22、23、32、54、55、64、86、87和96号系统元件，以各元件的原始不可用率为基准，比较两种NEIS法对系统元件重要程度的识别和凸显结果。从图7可以看出，尽管不进行离散变量区间映射的NEIS法可以对关键元件进行有效识别，但在系统元件不可用率的更新幅度上与NEIS法存在差距。NEIS法中关键元件的新不可用率更高，体现出离散变量区间映射方法对核密度估计边界效应的有效抑制作用。系统元件不可用率更新结果上的差异将导致仿真过程中不同的系统元件状态抽样效率，并对可靠性评估效率产生影响。

图8以CEIS法获取的系统负荷IS-PDF为参照，检验NEIS法和不进行离散变量区间映射的NEIS法对系统负荷重要区域的识别和凸显结果。可以看出，两种NEIS法获取的系统负荷IS-PDF(红色和绿色曲线)对系统负荷重要区域的范围，以及不同区域重要程度的判断结果近似，且均为不对称多峰分布。CEIS法(威布尔分布，τ＝30)与两种NEIS法对系统负荷重要区域范围的判断结果近似。

由表1可知，以CNMCS法得到的可靠性评估结果为基准，NEIS法和不进行离散变量区间映射的NEIS法均得到了较为精确的可靠性评估结果，其中可靠性指标LOLP的误差不超过5％，EENS的误差不超过1％。与NEIS-C法相比，两种NEIS法通过同时对系统元件和系统负荷进行重要抽样，进一步提高了大电网可靠性评估效率。但在离散变量区间映射方法的帮助下，NEIS法凭借更高的系统元件状态抽样效率，与不进行离散变量区间映射的NEIS法相比进一步减少了9％左右的计算时间，仅从计算时间来看已接近CEIS法对可靠性评估的加速效果。同时，凭借对系统负荷重要区域的合理凸显，NEIS法有效避免了CEIS法中系统负荷IS-PDF函数形式选择不当对可靠性评估精度的影响。以上结果证明了离散变量区间映射方法的有效性，以及NEIS法对IEEE-RTS96测试系统可靠性评估的加速效果。

表1不同方法得到的IEEE-RTS96测试系统可靠性评估结果

含多个风电场的IEEE-RTS96测试系统：

为进一步检验NEIS法在应对含有多维相关连续变量的大电网可靠性评估问题时的有效性，采用NEIS法进行系统元件、系统负荷和多维风速的IS-PDF联合估计，并采用式(33)计算系统元件的新不可用率。

图9以系统元件的原始不可用率为参照，选择不可用率变化较大的22、23、32、54、55、64、86、87和96号系统元件，检验NEIS法对系统元件重要程度的识别和凸显结果。可以看出，在离散变量区间映射方法的帮助下，NEIS法对系统元件不可用率进行了有效更新，这将显著提高仿真过程中系统故障状态样本被抽取的概率。为检验NEIS法对系统负荷重要区域和各节点风度重要区域的凸显结果，图10和图11分别以系统负荷原始概率密度函数，以及13、37和61号节点风速的原始概率密度函数为参照，将NEIS-C法和NEIS法获取的系统负荷IS-PDF和各节点风速IS-PDF进行函数图像对比。

从图10和图11可以看出，NEIS法与NEIS-C法得到的系统负荷IS-PDF，以及两种方法得到的各节点风速IS-PDF在一定程度上相互拟合，并且均为不对称多峰分布。两种方法对系统负荷和各节点风速重要区域的范围，以及不同区域重要程度的判断结果近似。为进一步检验NEIS法对各变量重要区域识别和凸显结果的合理性，以及与NEIS-C法相比，在可靠性评估效率上的提升效果，表2将CNMCS法、NEIS-C法和NEIS法的可靠性评估结果进行了对比。

表2不同方法得到的含3个风电场的IEEE-RTS96系统可靠性评估结果

由表2可知，以CNMCS法得到的可靠性评估结果为基准，NEIS-C法和NEIS法均得到了较为精确的可靠性评估结果，其中NEIS法得到的可靠性指标LOLP和EENS的误差均为4.3％左右。在可靠性评估效率方面，NEIS-C法的计算时间仅为CNMCS法的11％左右，而NEIS法通过对系统元件状态抽样效率的有效提高，大幅减少了仿真过程中需要抽取并分析的系统状态样本数，计算时间仅为NEIS-C法的9.1％左右。以上结果证明了NEIS法通过对离散变量和连续变量的IS-PDF联合估计，实现了大电网可靠性评估中系统状态IS-PDF的高效合理获取，有效弥补了NEIS-C法无法获取系统元件等离散变量IS-PDF的不足，在确保结果准确性的同时，进一步提高了大电网可靠性评估效率。

Claims

1.一种基于重要抽样函数改进非参数估计的电网可靠性评估方法，其特征在于，包括以下步骤：

2.根据权利要求1所述的方法，其特征在于，步骤2中包括：

步骤21、设定比例系数ρ，初始化迭代计数器k＝1，令

否则，采用以下方式对f(x_l|Ω_k)和f(x_w|Ω_k)进行无偏估计：

其中，W(x_[j])为系统状态x_[j]的似然比，有：

式中，为指示函数，若x_[j]∈Ω_k，/>若/> p(x_[j]∈Ω_k)是系统状态x_[j]∈Ω_k的概率；

3.根据权利要求2所述的方法，其特征在于，步骤1中还包括，选定一个定义域为(-∞,+∞)的概率密度函数h(·)，并得到对应的累积分布函数H(·)；改用分布区间为(-∞,+∞)的连续元件状态变量x_d,i原件状态，并将h(·)作为第i，1≤i≤D，个系统元件状态的原始概率密度函数，即f(x_d,i|Ω)＝h(x_d,i)，确定第i个系统元件正常状态与故障状态的分界点ε_i：

ε_i＝H^-1(u_i) (34)

步骤21中，初始化

否则，采用下式估计f(x_d,i|Ω_k)：

4.根据权利要求3所述的方法，其特征在于，概率密度函数h(·)为标准正态分布或对数正态分布概率密度函数。

5.根据权利要求2或3所述的方法，其特征在于，步骤24中，所述收敛条件为：或者从第k-1次迭代到第k次迭代时，/>的变化量小于20％。

6.根据权利要求2或3所述的方法，其特征在于，设定比例系数ρ＝10％。

7.根据权利要求2或3所述的方法，其特征在于，步骤24中，引入归一化因子NF，并采用如下所示的W'(x_[j])代替式(11)和式(12)中的系统状态似然比W(x_[j])：

8.根据权利要求1、2或3所述的方法，其特征在于，所述步骤3包括:

步骤34、根据下式计算失负荷概率LOLP和期望缺供电量EENS作为可靠性指标：：

9.根据权利要求8所述的方法，其特征在于，设置方差系数β_max＝0.2。