传统电网可靠性评估在参数为常数时求取可靠性指标的期望值和概率分布。期望值指标可从概率均值角度揭示电网的长期平均可靠性水平，而指标的概率分布则有助于更加深入地理解电网的随机特性及安全可靠性水平^[1]。在电网可靠性评估中，无论采用何种模型和算法，参数都是模型的基础输入数据，对模型输出结果有着深刻的影响。实际中，参数受元件类型、气候条件、使用寿命、样本资料不足等因素影响具有不确定性^[2-3]，若仅取其期望值作为元件随机规律的具体表征，则可能导致评估结果所表述的电网风险水平缺乏令人信服的准确性和代表性，因此其在风险管控中也缺乏实施决策判断的足够可信度。可见，探索参数的不确定性对电网可靠性评估结果的影响具有重要工程实用意义。

目前，针对参数不确定性对电网可靠性评估的影响研究已取得一定的成果。文献[4]用区间数描述参数，求解出区间数形式的可靠性指标。文献[5]用模糊数表达参数，计算出模糊数形式的可靠性指标。文献[6]将盲数理论引入可靠性评估，计算出盲数形式的可靠性指标。文献[7]利用联系数来描述参数，求解出联系数形式的可靠性指标。文献[8]利用概率密度分布来描述参数，实现了期望值风险指标的概率密度分布计算。纵观现有文献，主要探索参数不确定性对期望值指标的影响，而对指标概率分布特征的影响则鲜有报道。

笔者对电网可靠性概率分布受参数不确定性的影响开展研究，探索了3种计算方法，即双循环蒙特卡洛模拟法、点估计法和改进序贯蒙特卡洛模拟法，详细讨论了它们的原理及优缺点。为对比3种方法的计算结果和计算效率，对RBTS系统进行了评估，进一步利用改进序贯蒙特卡洛模拟法对IEEE-RTS79系统进行了评估，验证了其实用性。

1 参数不确定时电网可靠性概率分布的数学描述

在参数给定时，文献[9]采用序贯蒙特卡洛仿真求得电网可靠性指标Y的逐年样本信息，进一步采用非参数核密度估计^[10]计算Y的条件概率密度分布。若参数是随机变量，则Y的条件概率分布也将具有不确定性，它们之间的函数关系式为

$ {f_{Y\left| U \right.}}\left( {y\left| \mathit{\boldsymbol{u}} \right.} \right) = F\left( \mathit{\boldsymbol{U}} \right) = F\left( {{U_1},{U_2}, \cdots ,{U_m}} \right), $

(1)

式中：U为参数随机向量，m为其的维数，f_Y|U(y|u)为U=u条件下Y的条件概率密度分布。由概率论可知，当U不确定时，对Y可采用边缘密度分布f_Y(y)来完整揭示其概率分布特征，即

$ {f_Y}\left( y \right) = \int_{ - \infty }^{ + \infty } {{f_{Y\left| U \right.}}\left( {y\left| \mathit{\boldsymbol{u}} \right.} \right){f_U}\left( u \right){\rm{d}}\mathit{\boldsymbol{u}}} , $

(2)

式中f_U(u)表示随机向量U的联合概率密度分布。

在式(2)的计算上，笔者探索了3种计算思路，即双循环蒙特卡洛模拟法、点估计法和改进序贯蒙特卡洛模拟法。分别如节2~4所示。

2 双循环蒙特卡洛模拟法的基本原理

双循环蒙特卡洛模拟法的基本原理是外层循环对参数进行随机抽样，获得确切的参数后传入内层循环，内层循环以抽样所得参数为输入数据，进行电网可靠性序贯蒙特卡洛仿真，求取Y的条件概率密度分布。重复以上双循环过程，直至达到满足计算精度的总循环次数N，可获得N条Y的条件概率密度分布。此时式(2)可表示为

$ {f_Y}\left( y \right) \approx \frac{1}{N}\sum\limits_{i = 1}^N {{f_{Y\left| U \right.}}\left( {y\left| {{\mathit{\boldsymbol{u}}_i}} \right.} \right)} , $

(3)

式中：i为循环次数；u_i为第i次循环参数随机向量U的抽样值；N为序贯仿真循环次数。

单次序贯仿真采用缺供电量指标E_NS的样本集合Ω₁={E_NS1，E_NS2，…，E_NSn}(n为仿真年数)的方差系数β₁作为仿真的收敛判据。由于双循环蒙特卡洛模拟法需进行N次序贯仿真，导致上述收敛判据无法应用。笔者将N次序贯仿真所得的期望缺供电量指标E_ENS构成样本集合Ω₂={E_ENS1，E_ENS2，…，E_ENSN}，以该组样本集合的方差系数β₂作为循环次数N是否足够的判据，即

$ {\beta _2} = \frac{{\sqrt {V\left( {{{\bar I}_2}} \right)} }}{{{{\bar I}_2}}} = \frac{{{\delta _2}/\sqrt N }}{{{{\bar I}_2}}}, $

(4)

式中：I₂和δ₂分别表示样本集合Ω₂的均值和样本标准差。

3 点估计法的基本原理

点估计法^[11]的基本原理为：设参数随机向量U=(U₁, U₂, …，U_m)，通过随机变量U_l的概率密度分布f(u_l)可求得其高阶中心矩，再利用高阶中心矩构造出表征f(u_l)的K个估计点u_{l, k}(k=1, 2, …, K)。将m个随机变量各自的K个估计点组合，可得到U联合概率密度f_U(u)的K×m个估计点，如随机变量U_l取其第k个估计点u_{l, k}，其余随机变量均取均值，此时随机向量U=u(l, k)=(μ₁, μ₂, …, u_{l, k}, …, μ_m)，即为第(l, k)个估计点。此时式(2)可以表示为

$ {f_Y}\left( y \right) \approx \sum\limits_{l = 1}^m {\sum\limits_{k = 1}^K {{\omega _{l,k}}{f_{Y\left| U \right.}}\left( {y\left| u \right.\left( {l,k} \right)} \right)} } , $

(5)

式中：ω_{l, k}为u(l, k)对应的权重；f_Y|U(y|u(l, k))为随机向量U=u(l, k)条件下Y的条件概率密度分布。

随机变量U_l的第k个估计点u_{l, k}的表达式为

$ {u_{l,k}} = {\mu _l} + {\xi _{l,k}}{\sigma _l}, $

(6)

式中：μ_l和σ_l分别为随机变量U_l的均值和标准差，ξ_{l, k}为位置系数。ω_{l, k}和ξ_{l, k}满足关系

$ \left. \begin{array}{l} \sum\limits_{k = 1}^K {{\omega _{l,k}} = \frac{1}{m}} ,\\ \sum\limits_{k = 1}^K {{\omega _{l,k}}{{\left( {{\xi _{l,k}}} \right)}^j}} = {\lambda _{l,j}},j = 1, \cdots ,2K - 1。\end{array} \right\} $

(7)

式中：λ_{l, j}为随机变量U_l的j阶标准中心距，即U_l的j阶中心矩M_j(U_l)与标准方差σ_l的j次方之比：

$ {\lambda _{l,j}} = {M_j}\left( {{U_l}} \right)/{\left( {{\sigma _l}} \right)^j}, $

(8)

$ {M_j}\left( {{X_l}} \right) = \int_{ - \infty }^{ + \infty } {{{\left( {{u_l} - {\mu _l}} \right)}^j}f\left( {{u_l}} \right){\rm{d}}x} 。$

(9)

对于两点估计法(K=2)，ξ_{l, k}和ω_{l, k}如下：

$ \left. \begin{array}{l} {\xi _{l,1}} = \frac{{{\lambda _{l,3}}}}{2} + \sqrt {m + {{\left( {\frac{{{\lambda _{l,3}}}}{2}} \right)}^2}} ,\\ {\xi _{l,2}} = \frac{{{\lambda _{l,3}}}}{2} - \sqrt {m + {{\left( {\frac{{{\lambda _{l,3}}}}{2}} \right)}^2}} ; \end{array} \right\} $

(10)

$ \left. \begin{array}{l} {\omega _{l,1}} = - \frac{1}{m}\frac{{{\xi _{l,2}}}}{{{\xi _{l,1}} - {\xi _{l,2}}}},\\ {\omega _{l,2}} = \frac{1}{m}\frac{{{\xi _{l,1}}}}{{{\xi _{l,1}} - {\xi _{l,2}}}}。\end{array} \right\} $

(11)

由式(10)可见：随着m的增加，u_{l, 1}和u_{l, 2}将远离其均值μ_l，甚至可能超出U_l的定义域。为避免这种情况，可采用三点估计法(K=3)，ξ_{l, k}和ω_{l, k}如下：

$ \begin{array}{*{20}{c}} {{\xi _{l,k}} = \frac{{{\lambda _{l,3}}}}{2} + \left( { - 1} \right)3 - k\sqrt {{\lambda _{l,4}} - \frac{3}{4}\lambda _{l,3}^2} ,}\\ {k = 1,2,{\xi _{l,3}} = 0;} \end{array} $

(12)

$ \begin{array}{*{20}{c}} {{\omega _{l,k}} = \frac{{\left( { - 1} \right)3 - k}}{{{\xi _{l,k}}\left( {{\xi _{l,1}} - {\xi _{l,2}}} \right)}},}\\ {{\omega _{l,3}} = \frac{1}{m} - \frac{1}{{{\lambda _{l,4}} - \lambda _{l,3}^2}}。} \end{array} $

(13)

对于三点估计法，U的联合概率密度f_U(u)有3m个估计点，但因ξ_{l, 3}=0，估计点u_{l, 3}=μ_l，因此这3m个估计点中会出现m个相同的估计点(μ₁, μ₂, …, μ_l, …, μ_m)，对这m个相同估计点只需进行一次电网可靠性评估，其权重系数为

$ {\omega _0} = \sum\limits_{l = 1}^m {{\omega _{l,3}}} = 1 - \sum\limits_{l = 1}^m {\frac{1}{{{\lambda _{l,4}} - \lambda _{l,3}^2}}} 。$

(14)

基于点估计法的电网可靠性边缘概率密度分布的计算流程如图 1所示。

4 改进序贯蒙特卡洛仿真的基本原理

在电网可靠性评估中，元件的状态持续时间通常服从指数分布，表示为

$ {f_0}\left( t \right) = \lambda {{\rm{e}}^{ - \lambda t}},t \ge 0; $

(15)

$ {\varphi _0}\left( t \right) = \left( {1/r} \right){{\rm{e}}^{ - t/r}},t \ge 0。$

(16)

式中：f₀(t)和φ₀(t)分别为故障前工作时间T_F和故障后修复时间T_R的概率密度函数，参数λ和r分别为故障率和平均修复时间。其概率分布函数为

$ {F_0}\left( t \right) = 1 - {{\rm{e}}^{ - \lambda t}},t \ge 0, $

(17)

$ {\varphi _0}\left( t \right) = 1 - {{\rm{e}}^{ - t/r}},t \ge 0。$

(18)

传统序贯蒙特卡洛仿真需事先设置λ和r为给定常数(一般为期望值)，分别记为λ₀和r₀，此时f₀(t)和φ₀(t)实际为参数已知时T_F和T_R的条件密度函数，根据条件密度函数的定义有

$ f\left( {t\left| {{\lambda _0}} \right.} \right) = f\left( {t,{\lambda _0}} \right)/f\left( {{\lambda _0}} \right), $

(19)

$ \varphi \left( {t\left| {{r_0}} \right.} \right) = \varphi \left( {t,{r_0}} \right)/\varphi \left( {{r_0}} \right)。$

(20)

实际中参数受诸多因素影响而具有不确定性，本质上为随机变量，分别记为R和λ，此时T_F和T_R是参数已知的条件指数分布，但其具体形式却因其参数的不同取值而具有不确定性。此时，对于参数的任意取值λ和r，公式(19)、(20)均成立，故有

$ f\left( {t\left| \lambda \right.} \right) = f\left( {t,\lambda } \right)/f\left( \lambda \right), $

(21)

$ \varphi \left( {t\left| r \right.} \right) = \varphi \left( {t,r} \right)/\varphi \left( r \right)。$

(22)

式中：f(t, λ)为T_F和λ的联合密度函数，f(t|λ)为λ=λ时T_F服从的条件指数分布，f(λ)为λ服从的概率密度分布；φ(t, r)为T_R和R的联合密度函数，φ(t|r)为R=r时T_R服从的条件指数分布，φ(r)为R服从的概率密度分布。

将式(21)、(22)进行变换，可得

$ f\left( {t,\lambda } \right) = f\left( {t\left| \lambda \right.} \right) \cdot f\left( \lambda \right), $

(23)

$ \varphi \left( {t,r} \right) = \varphi \left( {t\left| r \right.} \right) \cdot \varphi \left( r \right)。$

(24)

将式(23)、(24)分别对R和λ进行积分，可求得T_F和T_R的边缘密度函数f(t)和φ(t)，进一步将f(t)和φ(t)分别对t积分可求得T_F和T_R的边缘分布函数F(t)和φ(t)，其具体表达式为

$ \begin{array}{*{20}{c}} {F\left( t \right) = \int_0^t {{\rm{d}}t} \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {f\left( {t,\lambda } \right){\rm{d}}\lambda } = }\\ {\int_0^t {{\rm{d}}t} \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {f\left( {t\left| \lambda \right.} \right) \cdot f\left( \lambda \right){\rm{d}}\lambda } ,} \end{array} $

(25)

$ \begin{array}{*{20}{c}} {\varphi \left( t \right) = \int_0^t {{\rm{d}}t} \int_{{r_{\min }}}^{{r_{\max }}} {\varphi \left( {t,r} \right){\rm{d}}r} = }\\ {\int_0^t {{\rm{d}}t} \int_{{r_{\min }}}^{{r_{\max }}} {\varphi \left( {t\left| r \right.} \right) \cdot \varphi \left( r \right){\rm{d}}r} ,} \end{array} $

(26)

附录阐述了F(t)的推导过程。

由式(25)、(26)可知，F(t)和φ(t)分别计及了λ和R的不确定性，将其纳入序贯仿真中，改进传统序贯仿真的状态持续时间抽样模型，从而实现参数不确定时电网可靠性指标边缘概率分布的计算。

F(t)和φ(t)的表达式与f(λ)和φ(r)的具体形式相关，当f(λ)和φ(r)的表达式较为复杂时，F(t)和φ(t)的解析表达式将非常复杂甚至根本无法给出具体的解析表达式。以T_F为例，当λ服从三角形分布时，F(t)的表达式见附录，对于T_F的概率抽样值F(t_i)=A(A为[0, 1]之间服从均匀分布的随机数)，难以通过F(t)的反函数t_i=F^-1(A)求出t_i。鉴于此，笔者采取数值积分将连续分布F(t)进行离散化表征，即令t在数值序列{t₁, t₂, …, t_n}上依次取值，计算相应时刻的概率值{F(t₁), F(t₂), …, F(t_n)}，直到F(t_n)的计算值达到1。当{t₁, t₂, …, t_n}之间的间隔足够小时，F(t_n)可近似代替故障前工作时间的边缘分布。图 2给出了T_F随机抽样的原理。

首先得到一个服从[0, 1]均匀分布的随机数A，若A满足

$ F\left( {{t_i}} \right) \le A \le F\left( {{t_{i + 1}}} \right),i \in \left[ {1,n - 1} \right], $

(27)

则相应的故障前工作时间抽样值为

$ t = \left( {{t_i} + {t_{i + 1}}} \right)/2, $

(28)

φ(t)的随机抽样与F(t)类似，不再赘述。

5 算例分析

为进行参数不确定性的量化分析，需事先假设参数服从的概率密度分布。文献[12]在计算概率潮流分布时假设线路参数服从均匀分布和正态分布；文献[13]利用随机集合理论处理参数不确定性时假设参数服从三角形分布文献；文献[14]在进行电子设备可靠性成本效益优化分析时假设元件故障率服从三角形分布和正态分布。参考上述文献，本文算例分析中采用了均匀、三角和正态分布描述参数的不确定性。若为实际工程系统，则可以依据专家经验和历史统计数据对参数的实际变化情况作出判断，并给出符合工程实际的概率分布类型。

为对比3种方法的计算结果和计算效率，对RBTS系统^[15]进行了可靠性评估。由于系统可靠性指标对线路L9的参数灵敏度最高，故先假设线路L9的参数服从均匀分布[a, b]，a，b分别取参数期望值的0.1倍和1.9倍。同时采用双循环蒙特卡洛模拟法(方法1)、三点估计法(方法2)和改进序贯蒙特卡洛模拟法(方法3)进行评估。单次序贯仿真的收敛条件为E_NS的方差系数小于3%，方法1总的收敛判据为E_ENS的方差系数小于1%。评估结果见图 3。

可见，3种方法的计算结果非常接近，但方法1需进行230次序贯仿真，方法2需5次，方法3只需1次，故方法2和方法3的效率远高于方法1。

上述评估仅计入了线路L9的参数不确定性，当计及多个元件参数的不确定性时，若采用方法1进行评估，为达到预定的计算精度其计算量将急剧增加。此时，可采用方法2和方法3进行评估。在参数分别取期望值、服从均匀、三角形和正态分布时对RBTS系统进行评估。均匀分布[a, b]和三角形分布[a, c, b]的a和b分别取参数期望值的0.1倍和1.9倍，三角形分布的众数c取参数期望值；正态分布N(μ, σ²)的期望值μ取参数期望值，标准差σ取参数期望值的0.3倍，序贯仿真的收敛条件为缺供电量E_NS指标的方差系数小于3%。评估结果见图 4~6。

可见，当参数服从不同分布时，两种方法的计算结果接近，但方法2需进行81次序贯仿真，而方法3只需一次，故其计算效率远高于方法2。当参数服从不同分布时，可靠性指标的概率密度分布与参数取期望值时相差较大，在均匀分布下尤为明显。根本原因在于参数服从均匀分布时，参数取分布范围两端数值的概率较其他分布类型有所增大，使得T_F取较大数值和T_R取较小数值的概率较大，导致可靠性指标为较小数值的概率也随之增大。

进一步对IEEE-RTS 79系统^[16]进行评估，该系统包含70个元件，若计及所有参数的不确定性，方法2需进行281次序贯仿真，计算量过于庞大，故采用方法3进行评估。假设参数取期望值、服从均匀分布[a, b]、三角形分布[a, c, b]和正态分布N(μ, σ²)，且参数a, b, c, μ和σ的取值与评估RBTS系统时相同，序贯仿真的收敛条件为缺供电量指标E_NS的方差系数小于3%。评估结果见图 7。可见，对于系统规模较大的IEEE-RTS79系统，采用方法3只需进行一次序贯仿真便能实现参数不确定时电网可靠性指标边缘概率密度分布的计算，证明了方法3的实用性。

6 结论

笔者探索了参数不确定时电网可靠性概率分布的3种计算方法，并比较了它们的计算结果和计算效率。评估结果表明，双循环蒙特卡洛模拟法虽然简单直观，但计算量过大，只能作为一种基准参考方法，验证其他方法的正确性和有效性。点估计法的计算效率较双循环蒙特卡洛模拟法有所提高，但当系统规模较大时，其计算效率仍然较低。改进序贯蒙特卡洛模拟法可高效计及多个元件可靠性参数的不确定性对电网可靠性概率分布的影响，具有较强的工程实用价值。

附录：元件故障前工作时间边缘分布函数F(t)的推导

1) 均匀分布。假设λ服从[a, b]上的均匀分布

$f\left( \lambda \right) = \frac{1}{{b - a}},\lambda \in \left[ {a,b} \right],$

代入式(25)得到

$F\left( t \right) = 1 - \frac{{{{\rm{e}}^{ - at}} - {{\rm{e}}^{ - bt}}}}{{t\left( {b - a} \right)}}。$

2) 三角形分布。假设λ服从三角形分布[a, c, b]：

$f\left( \lambda \right) = \left\{ \begin{array}{l} \frac{{2(\lambda - a)}}{{\left( {b - a} \right)\left( {c - a} \right)}},\lambda \in \left[ {a,c} \right],\\ \frac{{2\left( {b - \lambda } \right)}}{{\left( {b - a} \right)\left( {b - c} \right)}},\lambda \in \left[ {c,b} \right]。\end{array} \right.$

代入式(25)得到

$\begin{array}{c} F\left( t \right) = \frac{{2\left( {c - b} \right){{\rm{e}}^{ - at}} + 2\left( {a - c} \right){{\rm{e}}^{ - bt}} + 2\left( {b - a} \right){{\rm{e}}^{ - ct}}}}{{{t^2}\left( {a - b} \right)\left( {a - c} \right)\left( {b - c} \right)}} - \\ \frac{{a\left( {{b^2} - {c^2}} \right) + b\left( {{c^2} - {a^2}} \right) + c\left( {{a^2} - {b^2}} \right)}}{{\left( {a - b} \right)\left( {a - c} \right)\left( {b - c} \right)}} \end{array}$

3) 其他分布(如正态分布、Gamma分布等)。由于f(λ)的表达式过于复杂，无法用初等函数表示积分式$ \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {\lambda {{\text{e}}^{-\lambda t}}} f\left( \lambda \right){\text{d}}\lambda $，可采取数值积分计算。

将区间[λ_min, λ_max]n等分，步长$ h = \frac{{{\lambda _{\max }}-{\lambda _{\min }}}}{n} $，节点λ_i=λ_min+ih，i=0, 1, …, n，根据Simpson公式：

$\begin{array}{c} \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {\lambda {{\rm{e}}^{ - \lambda t}}} \cdot f\left( \lambda \right){\rm{d}}\lambda \approx \\ \sum\limits_{i = 0}^{n - 1} {\int_{{\lambda _i}}^{{\lambda _{i + 1}}} {{\lambda _i}{{\rm{e}}^{ - {\lambda _i}t}}} } \cdot f\left( {{\lambda _i}} \right) = \\ \frac{h}{6}\sum\limits_{i = 0}^{n - 1} {\left[ {{a_i}{{\rm{e}}^{ - {\lambda _i}t}} + {b_i}{{\rm{e}}^{ - {\lambda _i} + \frac{1}{2}t}} + {c_i}{{\rm{e}}^{ - {\lambda _{i + 1}}t}}} \right]} , \end{array}$

式中：$ {\lambda _{i + \frac{1}{2}}} = \frac{{{\lambda _i} + {\lambda _{i + 1}}}}{2} $，a_i=λ_if(λ_i)，$ {b_i} = 4{\lambda _{i + \frac{1}{2}}}f\left( {{\lambda _{i + \frac{1}{2}}}} \right) $，c_i=λ_i+1f(λ_i+1)。

$\begin{array}{c} F\left( t \right) = \int_0^t {{\rm{d}}t} \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {\lambda {{\rm{e}}^{ - \lambda t}}} \cdot f\left( \lambda \right){\rm{d}}\lambda \approx \\ \frac{h}{6}\sum\limits_{i = 0}^{n - 1} {[\frac{{{a_i}}}{{{\lambda _i}}}\left( {1 - {{\rm{e}}^{ - {\lambda _i}t}}} \right) + \frac{{{b_i}}}{{{\lambda _{i + \frac{1}{2}}}}}\left( {1 - {{\rm{e}}^{ - {\lambda _{i + \frac{1}{2}t}}}}} \right) + } \\ \frac{{{c_i}}}{{{\lambda _{i + 1}}}}\left( {1 - {{\rm{e}}^{ - {\lambda _{i + 1}}t}}} \right)]。\end{array}$

φ(t)的推导方法与F(t)相近，可参见F(t)的推导过程，此处不再赘述。