Power grid reliability evaluation method based on optimal important sampling density function successive smoothing estimation

Abstract

The invention relates to the field of large power grid reliability evaluation, in particular to a power grid reliability evaluation method based on optimal important sampling density function successive smoothing estimation. A more efficient IS-PDF is obtained by directly estimating the optimal zero variance IS-PDF step by step, the efficiency of the important sampling method is improved, and the method comprises the following steps: constructing a unified important sampling framework of all random variables in a multivariate standard normal space through a normal transformation method; on the basis of the important sampling framework, constructing an intermediate smooth joint IS-PDF sequence, taking truncation Gaussian distribution as initial distribution of the intermediate smooth joint IS-PDF sequence, carrying out step-by-step estimation on subsequent distribution in the intermediate smooth joint IS-PDF sequence in sequence, and meanwhile, carrying out sampling and power grid reliability index evaluation by using the estimated smooth joint IS-PDF; and taking the current power grid reliability index as a final evaluation index until the variance coefficient of the reliability index obtained during the Kth estimation meets the preset precision.

Description

Technical Field

[0001] The method relates to the field of large power grid reliability assessment, and specifically is a power grid reliability assessment method based on successive smoothing estimation of an optimal important sampling density function. Background Art

[0002] How to effectively speed up the Monte Carlo simulation and improve the efficiency of reliability assessment has long been a focus of research in the field of large power grid reliability assessment. Among the many variance reduction techniques, the cross-entropy based importance sampling method (CE-IS) has attracted much attention due to its excellent acceleration effect. As one of the most effective methods to improve the efficiency of reliability assessment, the traditional cross-entropy importance sampling method has been widely used in large power grid reliability assessment. The CE-IS method uses cross-entropy optimization to iteratively optimize the parameters of the importance sampling probability density function (IS-PDF), thereby providing strict theoretical rules and a universal implementation approach for effectively obtaining the IS-PDF.

[0003] The traditional cross entropy important sampling method also has some defects that are difficult to solve.

[0004] First, the traditional cross entropy important sampling method uses the theoretically optimal zero-variance IS-PDF as the target to perform cross entropy iterative optimization on the parameters of the actual IS-PDF. However, only samples in the fault domain can participate in the iterative optimization calculation of the IS-PDF parameters, so the parameter update of the traditional cross entropy important sampling method does not make full use of the samples. When the reliability of the system is low, it is relatively easy to obtain fault samples, so the traditional cross entropy important sampling method can obtain enough fault samples in a limited sample capacity for the iterative optimization calculation of the IS-PDF parameters. However, when the reliability of the system is high, it is very difficult to obtain fault samples, so the fault samples that the traditional cross entropy important sampling method can obtain in a limited sample capacity will become very scarce, or even there may be no fault samples, which will lead to the phenomenon of reduced efficiency or even failure of the iterative optimization of the IS-PDF parameters.

[0005] Secondly, the traditional cross entropy important sampling method is essentially a parameter estimation method. It is necessary to select the IS-PDF type of various random variables in advance before the IS-PDF parameter estimation can be performed. However, there is no strict theoretical rule as the basis for the selection of the IS-PDF type. Usually, the selection of the IS-PDF type is based on the consideration of the derivation of the parameter update formula and subjective experience. Since the IS-PDF g(x|Q) to be determined appears in the factor ln(g(x|Q)), the exponential family distribution is selected as the form of g(x|Q) to reduce the difficulty of deriving the parameter update analytical formula. In the reliability assessment of the power system, for the discrete variables in the system, that is, the state variables of the system components, their original PDF is a Bernoulli distribution belonging to the exponential distribution family. Therefore, the IS-PDF of such variables usually chooses the same function form, and the component unavailability rate is used as the parameter to be updated for iterative update. For continuous variables, due to their complexity and variability, the selection of the IS-PDF for such variables is mainly based on the feasibility of the derivation of the IS-PDF parameter update formula. For example, for system load, the normal distribution, lognormal distribution, and exponential family distributions such as GMM are usually selected as its IS-PDF. For the normal distribution and lognormal distribution, they are unimodal distributions with fixed shapes, so they may not necessarily highlight the important areas of the system load reasonably. For GMM, it can construct an IS-PDF of any shape, but its parameter update is very troublesome, which may reduce the computational efficiency of the traditional cross-entropy important sampling method. Therefore, it is difficult to select a suitable IS-PDF type, and the traditional cross-entropy important sampling method may suffer from a serious performance degradation due to improper selection of the IS-PDF type, and even cause serious errors in the reliability evaluation results.

[0006] In summary, the defects of the traditional cross entropy important sampling method are that its iterative optimization of IS-PDF parameters can only use fault samples and lacks strict theoretical rules as the basis for selecting the IS-PDF type. These defects will have different degrees of impact on the efficiency and performance of the traditional cross entropy method, leading to more serious problems in some occasions. Summary of the invention

[0007] The present invention proposes an important sampling method based on stepwise estimation of intermediate smooth distribution, and obtains a more efficient IS-PDF by directly and stepwise estimating the optimal zero-variance IS-PDF, thereby further improving the efficiency of the important sampling method.

[0008] The power grid reliability assessment method based on the successive smoothing estimation of the optimal important sampling density function in the present invention comprises:

[0009] Step 1: The original mixed variable space of discrete variables and continuous variables is transformed into a multivariate standard normal space through the normal transformation method, and on this basis, a unified important sampling frame of all random variables in the multivariate standard normal space is constructed;

[0010] Step 2: Based on the important sampling frame, an intermediate smoothed joint IS-PDF sequence is constructed, with the truncated Gaussian distribution as its initial distribution, and the subsequent distributions in the intermediate smoothed joint IS-PDF sequence are estimated step by step in sequence, and the estimated smoothed joint IS-PDF is used for sampling and grid reliability index evaluation; until the coefficient of variance (Cov) of the reliability index obtained at the Kth estimation meets the preset accuracy, the current grid reliability index is used as the final evaluation index.

[0011] In order to avoid the shortcomings of the traditional CE-IS method and the iCEM method in selecting IS-PDF, the present invention proposes an important sampling method based on the stepwise estimation of the intermediate smooth distribution, and obtains an IS-PDF similar to the theoretical optimal IS-PDF and more efficient without relying on the cross entropy optimization method. The method first establishes a unified important sampling framework of all random variables in the multivariate standard normal space, and realizes the continuity of the original distribution f(x|P). Then, under the important sampling framework, a set of smooth intermediate joint IS-PDF sequences are constructed in combination with the smooth function, and the intermediate joint IS-PDF sequence is sampled and estimated step by step by the sample resampling method and the conditional MH sampling method. As the estimation proceeds, the intermediate joint IS-PDF sequence will gradually and smoothly approach the theoretical optimal IS-PDF, and the theoretical optimal IS-PDF is directly approximated, and the joint IS-PDF of all random variables is obtained. Since the method does not need to select the type of IS-PDF and does not rely on the cross entropy optimization method to update the parameters, the shortcomings of the important sampling method based on the cross entropy are avoided, and the efficiency of reliability evaluation is improved. Numerical analysis shows that the importance sampling method based on stepwise estimation of intermediate smooth distribution can achieve an approximate estimation of the theoretical optimal IS-PDF and obtain a more reasonable IS-PDF, which has better performance and efficiency than the cross entropy importance sampling method.

[0012] In a further embodiment, a sampling method based on a sample resampling operation and a conditional MH sampling technique free from the curse of dimensionality is proposed to efficiently sample the high-dimensional joint IS-PDF.

[0013] In a further embodiment, a comprehensive reliability index is designed to fully utilize the samples extracted from all intermediate joint IS-PDFs to further improve the efficiency of SSAM. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 A schematic diagram of changes in the marginal IS-PDF of the system load in the SSAM method in an embodiment of the present invention.

[0015] Figure 2 Schematic diagram of comparison between the original PDF of system load and IS-PDFs obtained by different methods in an embodiment of the present invention.

[0016] Figure 3 In the SSAM method in the embodiment of the present invention, different The smoothing parameter σ k and the normalization parameter p k Schematic diagram of the changes in .

[0017] Figure 4 The figure is a comparison diagram of the acceptance rates of conditional MH sampling and traditional MH sampling in the SSAM method in the embodiment of the present invention. DETAILED DESCRIPTION

[0018] The power grid reliability assessment method (SSAM) based on the successive smoothing estimation of the optimal importance sampling density function in this embodiment includes the following steps:

[0019] 1. Uniform importance sampling framework for all random variables in multivariate standard normal space

[0020] Directly calculate the optimal zero-variance IS-PDFg op (x|Q op )=H(x)f(x|P) / p F Sampling estimation is not feasible because it is quite difficult to estimate the unknown normalization parameters, especially when the system being evaluated has a high reliability. In addition, even if p is obtained F It is also difficult to use the optimal zero-variance IS-PDF for Monte Carlo sampling because not only does the indicator function H(x) have no available analytical form, but the factor f(x|P) is essentially a complex distribution that includes both discrete and continuous variables. Therefore, even the MH sampling method, which can obtain samples from complex or possibly non-analytical distributions, is difficult to sample from the optimal zero-variance IS-PDF because the MH sampling method is only applicable to samples from continuous distributions.

[0021] In order to solve the above problem, it is necessary to make H(x) and f(x|P) continuous. This example proposes a smoothing function Φ(-S(x) / σ), which replaces the original indicator function H(x) with a function with a smaller smoothing parameter σ K The smoothing function of realizes the smoothing of the indicator function and obtains the smoothed optimal zero-variance IS-PDF:

[0022]

[0023] Therefore, in this example, we will focus on solving the problem of f(x|P) continuity, and transform the original probability space that obeys f(x|P) into a probability space that obeys An independent standard normal space is constructed, and a unified important sampling framework for all random variables in the multivariate standard normal space is established. The specific contents are as follows:

[0024] 1) System component state vector x c Each discrete variable in can be converted into the corresponding standard normal variable by formula (2)

[0025]

[0026] Among them, δ j =Φ -1 (ξ j ) is the threshold between the fault state and the normal state. ξ j is the original unavailability rate of the jth component, Φ -1 (·) is the inverse of the standard normal cumulative distribution function. is an indicator function: if but otherwise,

[0027] 2) System load variable x L Transformed into the corresponding standard normal variable u through equal probability transformation L :

[0028]

[0029] Among them, F L -1 () is the original cumulative distribution function The inverse function of .

[0030] 3) Wind speed vector with correlation Can be converted into corresponding independent standard normal random vectors by Nataf transformation

[0031] First, x is transformed into w Transformed into a zero-mean vector and a correlation matrix C z The normal random vector z is in is the cumulative distribution function The inverse function of .

[0032] Let C x yes The correlation matrix, ρ x (i,j) and ρ z (i,j) are C x and C z The element in row i and column j in ρ z (i,j) can be obtained by solving the following formula:

[0033]

[0034] Among them, μ i and μ j They are and The mean value, σ i and σ j They are and The standard deviation of . is a correlation coefficient C z The bivariate standard normal probability density function of (i,j). After obtaining C z Then, C is calculated by Cholesky decomposition. z The lower triangular matrix L z , that is, C z =L z L z T , then the correlated normal random vector z can be obtained from the independent standard normal random vector u w Indicates that z = L z u w Therefore, the wind speed vector with correlation The corresponding independent standard normal random vector The conversion relationship between them can be expressed by formula (5):

[0035]

[0036] In summary, the above normal transformation converts the original x space containing discrete random variables, continuous variables, and correlated random variables into the u space containing only independent standard normal random variables, and transforms the discrete and complex original distribution f(x|P) into a simple and continuous multivariate standard normal distribution. x=[x c ,x w ,x L ] and u=[u c ,u w ,u L ] can be simply expressed as x=T(u), which means that each sample extracted from the u space can be mapped to a sample in the original x space.

[0037] Step 2: Based on the important sampling frame, an intermediate smoothed joint IS-PDF sequence is constructed, with the truncated Gaussian distribution as its initial distribution, and it is gradually estimated; as the estimation proceeds, the intermediate distribution sequence will gradually approach the optimal zero-variance IS-PDF, thereby achieving an approximate estimation of the optimal zero-variance IS-PDF, and finally obtaining a joint IS-PDF that is highly similar to the optimal zero-variance IS-PDF in functional form, as follows:

[0038] According to the above sampling framework, the optimal IS-PDF g in the discrete x space is op (x|Q op ) is converted to the optimal g in the continuous u space op (u|Q op ):

[0039]

[0040] in Represents the n-dimensional standard normal probability density function, n=D+W+1. and They represent D-dimensional, W-dimensional and one-dimensional standard normal probability density functions respectively. S(u) is the system performance function defined in u-space and can be determined by power flow or optimal power flow analysis of the system state x=T(u) in x-space.

[0041] Combining formula (6) with the smoothing function, we can obtain the smoothed optimal zero-variance IS-PDF:

[0042]

[0043] Among them, σ k is a sufficiently small smoothing parameter, p k is the normalization parameter, Q k ={σ k ,p k}. Below, we focus on the principle and steps of the important sampling method based on the stepwise estimation of the intermediate smooth distribution sequence. Although the smooth optimal zero-variance IS-PDF g(u|Q k ), but the exact estimation of the smoothing parameter σ k and the normalization parameter p k is very difficult, so how to achieve g(u|Q k) is a thorny problem. Therefore, based on the idea of ​​stepwise smoothing estimation, the present invention proposes an important sampling method (SSAM) based on stepwise estimation of intermediate smooth distribution to solve this problem. The basic idea of ​​SSAM is to establish a set of intermediate smoothing joint I SP DF sequences {g0(u|Q0),g1(u|Q1),...,g k (u|Q k )}, where the initial distribution g0(u|Q0) is a distribution whose specific form is known and easy to sample.

[0044] Then, using the distribution g obtained from the previous k-1 (u|Q k-1 ) to predict the next unknown distribution g k (u|Q k ) parameter Q k ={σ k ,p k} is estimated. As the estimation proceeds, the smoothing parameter will continue to decrease, that is, {∞=σ0>σ1>…σ K >0}, and the intermediate smoothed joint IS-PDF sequence {g k (u|Q k ),k=0,1,...,K} will gradually approach the optimal zero-variance IS-PDF, that is, g k (u|Q k )≈g op (u|Q op ).

[0045] Note that the implementation of SSAM needs to solve two key problems: 1) How to use the distribution g obtained from the previous k-1 (u|Q k-1 ) to estimate the next unknown distribution g k (u|Q k ) parameter Q k ={σ k ,p k}; 2) In order to obtain Q k After that, how to generate g k (u|Q k ) to ensure that the stepwise estimation process can be continued. The following will introduce the solutions to these two key problems in detail. In addition, in the following, the function symbol g k (u|Q k ) will be abbreviated to g k (u) to make the expression more concise.

[0046] (1) Parameter Q k Estimation method

[0047] G k (u) is rewritten as follows:

[0048]

[0049] in,

[0050] It can be found by observation that the intermediate smoothed joint IS-PDF sequence {g k (u), k=0,1,…,K} has the same smooth distribution sequence as the intermediate one {h k (x|B k ),k=0,1,…,K} of the same form, and g k (u) and h k (x|B k ) contain the same unknown parameters, namely the smoothing parameter σk and the normalization parameter pk. Therefore, g k The parameter estimation method of (u) is similar to h k (x|B k ), the specific process is as follows.

[0051] Normalization parameter p k With p k-1 The ratio is:

[0052]

[0053] Ratio η k An approximate estimate of is:

[0054]

[0055] Among them, u i is the intermediate smoothed joint IS-PDFg from the previous known k-1 (u) is the i-th sample drawn, and N is the number of samples drawn.

[0056] The estimated accuracy of The coefficient of variance (Cov) is used to measure:

[0057]

[0058] The approximate expression is:

[0059]

[0060] In getting after, It can be calculated by the following formula:

[0061]

[0062] As shown in formula (11), An estimate of π requires knowing the exact expression k (u)=Φ(-S(u) / σ k ) / Φ(-S(u) / σ k-1 ). However, π k There is an unknown smoothing parameter σ in (u) k Therefore, when using formula (11) to estimate Before that, we must first get σ k It should be noted that in order to ensure The estimated accuracy and variance coefficient Should satisfy a pre-specified target value. As shown in formulas (10)-(13), if given the previously known intermediate smoothed joint IS-PDF g k-1 The sample u drawn from (u) i ,1≤i≤N, then is a function with only a single variable σ k Taking this into account, the unknown parameters σ can be obtained by solving the following optimization model k :

[0063]

[0064] in is the target value, and ||·||2 represents the 2-norm.

[0065] In this method, in order to avoid a significant impact on the accuracy of the subsequent estimation process, the parameter σ k and p k Therefore, this method is effective for There are relatively strict restrictions on the value of . It is worth noting that The value of should be neither too large nor too small. May lead to The inaccurate estimation leads to the subsequent k and g k (u) estimated inaccuracy. However, the smaller This may result in obtaining an intermediate smoothed joint IS-PDF g k (u) Combined with the previous intermediate smoothing IS-PDF g k-1 (u) are too close, which requires continuous estimation of more intermediate smoothed joint IS-PDFs to finally achieve g k (u) for the theoretical optimal g op(u) gradually approaches, which greatly increases the computational cost. Usually, A balance between accuracy and efficiency can be achieved.

[0066] (2) For intermediate smoothing combined IS-PDF g k Sampling method of (u)

[0067] In addition to the initial intermediate joint IS-PDF g k (u), using traditional sampling methods such as inverse transform sampling, acceptance-rejection sampling, etc., from other intermediate joint IS-PDFg k (u), k ≥ 1, it is extremely difficult to extract samples from the performance function S(u); this is because it is impossible to obtain the intermediate joint IS-PDF g that contains the performance function S(u) k (u), k ≥ 1. To solve this problem, a sampling method based on sample resampling and MH sampling technique that avoids the curse of dimensionality is proposed.

[0068] Assume that we have already smoothed the joint IS-PDF g from the previous intermediate k-1 (u) obtained N iid samples {u1,u2,...,u N}. Then sample u i The weight ω(u i ) can be defined by the following formula:

[0069]

[0070] If u i =u j ,but otherwise, use To express u i It can be proved that when N approaches infinity, (N(u i ) / N)ω(u i ) will approach g k (u i ). Based on this, N can be obtained by resampling. c <N obey g k (u) seed sample: with ω(u i ), 1≤i≤N as u i ,1≤i≤N, from the total sample {u1,u2,...,u N} randomly select N with replacement c samples.

[0071] After obtaining Nc seed samples {u [1] ,u [2] ,...,u[Nc]}After that, MH sampling will be started from these seed samples, while simulating N c The length of the strip is L = N / N c For each Markov chain, from the current known sample u i Get the target distribution g k The next random sample u in (u) i+1 The sampling process is as follows:

[0072] 1) From the proposed distribution function q(θ|u i ) generates a candidate sample θ.

[0073] 2) Calculate the ratio:

[0074]

[0075] 3) P a =min{1,r(u i ,θ)} Accept θ: Generate a uniformly distributed random number r in the interval (0,1). If r is less than P a , then set u i+1 =θ; otherwise, set u i+1 =u i .

[0076] In general, choose the current sample u i is a symmetric probability density function centered on q(θ|u i )=q(u i The probability density function of |θ) is used as a proposed distribution function, such as a multivariate Gaussian probability density function. Therefore, the ratio in equation (17) is simplified to:

[0077]

[0078] As shown in formula (18), if the probability density function If the dimension n is very high, then the ratio may be very small, resulting in r(u i ,θ) is also very small. For this reason, the probability of the MH algorithm obtaining duplicate samples is very high, resulting in low efficiency. In order to solve this problem, the conditional MH sampling method is introduced, and its specific content is described as follows.

[0079] In the conditional MH sampling method, the current state vector u is selected i The multivariate Gaussian probability density function conditional on is used as the proposed distribution function:

[0080]

[0081] Among them, I is the identity matrix and γ∈[0,1] is a coefficient that can be adaptively adjusted.

[0082] Substituting formula (19) into formula (18), r(u i ,θ) can be further simplified as:

[0083]

[0084] It can be seen that there is no longer a ratio in formula (20) And r(u i ,θ) becomes a ratio that depends only on the one-dimensional smooth function. i ,θ) depends only on the current state u i The similarity between the system performance function and the candidate state θ is not the dimension of the random variable space, so the conditional MH sampling method is suitable for smoothing the joint IS-PDFg from high-dimensional intermediate k (u) is sampled.

[0085] It is worth noting that the choice of parameter γ affects the performance of the MH algorithm. If γ is chosen too small, many candidate samples will be rejected; if γ is chosen too large, the obtained samples will have too high correlation. Therefore, an adaptive method is used to adjust γ:

[0086]

[0087] in, is the newly generated N c The average acceptance probability of a sample is the ratio of the number of accepted candidate samples to the total number of candidate samples. a The selection of ζ should ensure that the change of γ gradually decreases with each adjustment. For example, we can set ζ a =(i a +1) -1 , where i a Represents the number of adjustments performed so far. γ max ≤1 is the upper limit of γ. For example, we can set γ max =0.99. * =0.44 is the ideal acceptance probability of the MH algorithm.

[0088] (3) Initial intermediate smoothing joint IS-PDF g0(u)

[0089] Starting from the known initial intermediate joint IS-PDF g0(u), the estimated intermediate smoothed joint IS-PDF will gradually move towards the failure region and gradually approach the theoretical optimal IS-PDF g op (u).

[0090] Obviously, the choice of the initial intermediate joint IS-PDFg0(u) has a certain impact on the performance of SSAM. Usually, we can choose As g0(u), but this is only valid when evaluating a system with low reliability. However, if the reliability of the system being evaluated is high, will be away from the failure domain, so that In addition, in order to avoid excessive computational burden in the estimation process, the number of samples N generated from each intermediate smoothed joint IS-PDF is limited, so The N samples generated in the stepwise estimation process may not contain any faulty samples. Since the initial stage of the stepwise estimation process lacks faulty sample information, the intermediate IS-PDFg k (u), 1≤k≤K, the movement to the fault domain may become slow and difficult. Therefore, the final g K (u) may be mainly concentrated in the fault area where the performance function S(u) does not perform too badly, resulting in inaccurate estimated reliability indicators. From this perspective, under the finite sample size N, the initial intermediate smoothed joint IS-PDF g0(u) that can easily generate fault samples is a better choice. Based on this idea, this example constructs a truncated normal distribution that ignores small load levels and uses it as g0(u) to make it easier to obtain fault samples. Specifically, in this example, the truncated normal probability density function is used as the initial intermediate joint IS-PDF, where is the truncated system load distribution, which takes the form shown below.

[0091]

[0092] in, is a function that depends on the truncated load level (for example ) threshold.

[0093] The parameters of g0(u) are Among them, F L (x L ) is the original cumulative distribution function of system load. It is worth noting that the truncated load level Can't be too big. Very large, in the initial stage of SSAM, a large number of fault samples with extremely poor performance function S(u) may be extracted from g0(u). Although the intermediate joint IS-PDF estimate estimated later can quickly approach the fault domain, the final g K (u) may be overly concentrated in the fault areas where S(u) performs extremely poorly, which may lead to deviations in the estimated reliability indicators.

[0094] Therefore, in some embodiments of SSAM, the intermediate joint IS-PDF sequence {g1(u),…,g K (u)} until the variance coefficient (COV) of the reliability index obtained at the Kth estimation meets the preset accuracy, that is, the variance coefficient is less than or equal to β max Based on the final g K The estimated values ​​of the load loss probability (LOLP), expected short supply (EENS) and its COV of (u) are as follows:

[0095]

[0096] Among them, u i is the joint IS-PDFg from the middle K The samples drawn from (u), N K is the number of samples drawn, L c (u i ) represents the system state u determined by optimal load shedding analysis i The load reduction amount, is the likelihood ratio.

[0097] In summary, in some embodiments of the present invention, the SSAM algorithm can be implemented according to the following specific steps:

[0098] (1) Setting Target value and the maximum variance coefficient β of the reliability index max , and set k=0.

[0099] (2) From the initial intermediate joint IS-PDF N samples {u1,…,u N}, where σ0=∞, p0=0.5.

[0100] (3) Let k = k + 1. Using the g known from the previous k-1 (u) is the N samples {u1,…,u N}, use formula (15) to calculate the next intermediate joint IS-PDFg k The smoothing parameter σ of (u) k .

[0101] (4) Estimate the ratio using formula (11)

[0102] (5) Use formula (14) to calculate g k Normalization parameter of (u)

[0103] (6) Using formula (16) for g k-1 (u) samples {u1,…,u N}Calculate the weights {ω(u1),…,ω(u N )}. Perform resampling to obtain distribution g k (u)N c Seeds {u [1] ,u [2] ,…,u [Nc]}. Let j = 1, L = N / N c .

[0104] (7) Run N in parallel c Markov chain simulations are performed, and each Markov chain is pushed forward one state according to the conditional MH sampling algorithm to obtain N c Then, these samples are used to update the parameter γ through formula (21), and N k =j·N c .

[0105] (8) Let K = k, then evaluate the reliability indicators LOLP and EENS by using formulas (23) and (25) respectively, and then obtain by using formulas (24) and (26): and Cov(EENS). If Then set j = j + 1 and go to step (9); if Then the intermediate joint IS-PDF g k (u) is the last distribution of the sequence and outputs the calculation results and ends the algorithm.

[0106] (9) If j≤L, go to step (7); otherwise, obtain N=L·N c samples as the k (u) and go to step (3).

[0107] In order to further improve the convergence speed of SSAM, some embodiments of the present invention also propose a comprehensive reliability index that can fully utilize the samples extracted from all intermediate joint IS-PDFs.

[0108] Assume that, as SSAM proceeds, the IS-PDF g is combined from the middle k (u), 0≤k≤K, N is generated k Sample u i ,(1≤i≤N k ). Based on the k (u), 0≤k≤K, the expected value and variance of various reliability indicators (such as LOLP, EENS, etc.) can be estimated by the following unified expression:

[0109]

[0110] in, is the likelihood ratio, u i From g k (u), F(u) is the reliability index function in u space. If F(u) = I {S(u)＜0} , This is the LOLP indicator. If F(u)=E c (u), It represents the EENS index.

[0111] Assume κ k for The weight of The definition of the comprehensive reliability index is as follows:

[0112]

[0113] Definition The weight κ with the smallest variance k , 0≤k≤K, is the optimal weight, and its specific form is as follows:

[0114]

[0115] In the optimal weight κ k (0≤k≤K), Variance for:

[0116]

[0117] Coefficient of variance It can be estimated by the following formula:

[0118]

[0119] The exemplary steps of the SSAM algorithm using the comprehensive reliability index are as follows:

[0120] Step (1) Set the following parameters:

[0121] ① The number of the next estimated intermediate joint IS-PDF is k=K=1, the maximum number of samples drawn from the intermediate joint IS-PDF each time is N (usually 10000-20000), and the number of seed samples N c (e.g. 10%-20% of N), and the Markov chain length L = N / N c ;

[0122] ② Maximum variance coefficient β max(usually 1%-2%), Target value (usually 1%-2%);

[0123] ③Truncation load level (e.g. 0.8pu), the initial value of the parameter γ (usually 0.5).

[0124] Step (2) transforms the x-space with different types of random variables into the u-space with independent standard normal random variables.

[0125] Step (3) uses truncation As the intermediate joint IS-PDF, Extract N0=N samples {u1,…,u N}, and by transforming the system state x in x space i =T(u i ) to perform power flow or optimal power flow analysis to determine S(u i )(1≤i≤N s ).

[0126] Step (4) calculate g in sequence k The smoothing parameter σ of (u) k ,ratio and normalization parameters

[0127] Step (5) calculates the sample {u1,...,u N}~g k-1 The weights of (u) {ω(u1),...,ω(u N )}. Let ω(u1) be u i The probability of occurrence of , and perform a resampling operation, through which {u1,...,u Ns} randomly draw N with replacement c Seeds Assume j=1.

[0128] Step (6) runs N c Markov chain simulations are performed, and each Markov chain is pushed forward one state according to the conditional MH sampling algorithm to obtain N c In the sampling process, the candidate sample θ is extracted from the proposed PDF in formula (19), S(θ) is determined by the power flow or optimal power flow analysis of the system state x = T(θ) in the x space, and the acceptance rate r(u i ,θ) is calculated by (20).

[0129] Step (7) uses formula (21) to use the newly generated N cThe parameters γ are updated by N samples, and N k =j·N c .

[0130] Step (8) uses formulas (27) and (28) to calculate the expected values ​​of LOLP and EENS indicators respectively. and variance

[0131] Step (9) calculates the expected values ​​of LOLP and EENS indicators by using formula (30) The weight κ k , (0≤k≤K).

[0132] Step (10) calculates the comprehensive reliability index of LOLP and EENS indexes by formulas (29), (30) and (31): variance and the coefficient of variance

[0133] Step (11) If the LOLP indicator Greater than β max , then go to step (12); otherwise, go to step (13).

[0134] Step (12) Set j = j + 1. If j ≤ L, go to step (7); otherwise, obtain N = L·N c samples as the k (u), and set k=k+1, K=k, and then go to step (4).

[0135] Step (13) If the EENS index Less than β max , then output the obtained reliability index and stop the simulation; otherwise, set j=j+1 and go to step (6).

[0136] Calculation Example

[0137] In this example, the reliability of the IEEE-RTS79 and IEEE-RTS96 test systems and their modified versions was evaluated. The cut-in wind speed, rated wind speed and cut-out wind speed of the wind turbine are 3m / s, 13m / s and 25m / s respectively. This example comprehensively investigates the efficiency and performance of the original non-sequential Monte Carlo simulation (NCMCS), the SSAM method proposed in this invention and the improved cross entropy method (iCEM) in reliability evaluation. In iCEM, normal, lognormal and GMM are selected as the IS-PDF of the system load, and three different iCEM methods are constructed based on this, namely, iCEM N 、iCEM L and iCEMG , and set the multi-level number parameter ρ a =0.1, target value In SSAM, the target value, seed sample size, and cutoff load level are set as N c =0.2N, In addition, unless otherwise specified, the sample size N is set to 20,000 by default for each parameter update in iCEM or for each intermediate joint IS-PDF parameter estimation in SSAM. In addition, the maximum variance coefficient β is set max = 0.02 was used as the convergence criterion for all methods.

[0138] (1) Reliability evaluation results of the IEEE-RTS79 test system and its modified version

[0139] a. Original IEEE-RTS79 test system

[0140] In order to verify the effectiveness of the proposed SSAM method, the reliability of the original IEEE-RTS79 is first evaluated. Since the reliability of the original IEEE-RTS79 is not very good, the number of samples N required for each iteration is set s =10000 is more appropriate. The evaluation results are shown in Table 1, where K represents the total number of iterations in iCEM or the total number of consecutive estimates in SSAM, and N K N is the number of samples required for the main sampling phase of iCEM or the last estimation phase of SSAM. t =K·N+N K Represents the total sample size during the entire simulation.

[0141] Table 1 IEEE-RTS79 system reliability evaluation results obtained by different methods

[0142]

[0143] As shown in Table 1, all iCEM N 、iCEM L 、iCEM G Both the SSAM and NCMCS methods can obtain reliability indices similar to those of NCMCS, but their total sample size N tand CPU time T are significantly less than NCMCS because IS-PDF significantly improves their ability to capture rare fault events. The results also show that SSAM has higher efficiency than CEM. In addition, the SSAM method proposed in this paper can obtain a "real" joint IS-PDF with a similar functional form to the theoretical zero-variance IS-PDF, while the iCEM method can only obtain a "pseudo" joint IS-PDF that simply combines the IS-PDFs of all random variables. In order to test the gap between the actual IS-PDF and the theoretical zero-variance IS-PDF, this paper uses some commonly used distance measurement methods, such as Kullback-Leibler distance, Hellinger distance, Pearsonχ 2 The metrics and α-divergence (parameter α = 0.5, 0.0 and -0.5) are used to measure the closeness between the IS-PDF obtained in the fault domain and the theoretical optimal IS-PDF.

[0144] As shown in Table 2, SSAM is significantly better than the three iCEM methods because all distance metrics of SSAM are smaller than the corresponding metrics of the iCEM methods, which indicates that SSAM can obtain an IS-PDF that is closer to the theoretical optimal IS-PDF.

[0145] Table 2 Distances between the IS-PDFs obtained by different methods and the optimal IS-PDF

[0146]

[0147] b. IEEE-RTS79 test system to reduce peak load

[0148] In order to demonstrate the superiority of SSAM in a more reliable scenario, the original IEEE-RTS79 was modified to reduce the load level to 80% of the original, which greatly improved the reliability. This modified IEEE-RTS79 is called MRTS79-1, and the evaluation results are shown in Table 3.

[0149] Table 3 Reliability evaluation results of MRTS79-1 system obtained by different methods

[0150]

[0151] As shown in Table 3, the LOLP index is on the order of 10 -5 , which indicates that the reliability of MRTS979-1 is very high, so the total sample size N required using the NCMCS method is tCompared with the NCMCS method, which suffers from heavy computational burden, the iCEM and SSAM methods still show significant performance advantages, even if the evaluated system has extremely high reliability. While ensuring accurate reliability indicators, the total sample size N required for evaluation is greatly reduced. t and CPU time T. However, the results in Table 3 show that iCEM G The computational efficiency is inferior to other iCEM methods. The root cause is that iCEM G It relies on a complex multi-step EM algorithm to update the parameters of the IS-PDF constructed based on GMM. Compared with the iCEM method, SSAM reduces the number of samples and computation time required while obtaining accurate reliability indicators. Unlike the iCEM method, which requires K = 3 parameter iterations, SSAM only requires K = 2 consecutive estimates of the intermediate joint IS-PDF.

[0152] Since the intermediate joint IS-PDF is a high-dimensional PDF, it is difficult to visualize its actual movement. Figure 2 The marginal PDF curves of the system load derived from the high-dimensional intermediate joint IS-PDF are plotted in Figure 2 to further intuitively demonstrate the performance of the SSAM method.

[0153] exist Figure 1 middle, and g0(u L ) are the original and truncated PDFs of the system load in u space, respectively, g1(u L ) and g2(u L ) are the marginal PDFs of system load derived from the estimated intermediate joint IS-PDF, and f(x L ),g0(x L )~g2(x L ) is the corresponding PDF in x space. Figure 1 As shown in Figure 2, as the estimation proceeds, the density peak of the system load edge PDF will gradually move to the risk area (large load area), which shows that the SSAM method can effectively identify and highlight the important areas of the variables.

[0154] c. IEEE-RTS79 test system with multiple wind farms

[0155] In order to verify the performance of SSAM in scenarios with correlated random variables, IEEE-RTS79 was modified as follows:

[0156] ①Increase the load level to 110% of the original level.

[0157] ② Six 200MW wind farms were added on nodes 13, 15, 16, 18, 21 and 23.

[0158] The modified version is named MRTS79-2, and the number of samples required for each iteration is set to N = 15000. The reliability evaluation results of MRTS79- are shown in Table 4.

[0159] Table 4 Reliability evaluation results of MRTS79-2 system obtained by different methods

[0160]

[0161] As shown in Table 4, the SSAM method can obtain the reliability index as accurate as the CEM method, but the required sample size N is t The CPU time T used is also shorter. In addition, the number of consecutive estimations required in SSAM is K = 2, while the number of iterative parameter updates required in iCEM is K = 3, which shows that SSAM is more efficient than iCEM in finding IS-PDF.

[0162] (2) Reliability evaluation results of IEEE-RTS96 test system and its modified version

[0163] a. Original IEEE-RTS96 test system

[0164] In order to further demonstrate the effectiveness of the SSAM method in reliability analysis of power systems with higher reliability and larger scale, this example further conducts a reliability evaluation on the IEEE-RTS96 test system. The evaluation results are shown in Table 5.

[0165] Table 5 IEEE-RTS96 system reliability evaluation results obtained by different methods

[0166]

[0167] As can be seen from Table 5, when conducting reliability assessment on power systems with higher reliability and larger scale, the NCMCS method needs to extract a large number of system status samples, which increases the computational burden and causes it to take a considerable amount of time to complete the reliability assessment. Thanks to the acceleration effect of the important sampling method on reliability assessment, the iCEM method and the proposed SSAM method converge quickly, and the CPU time used is approximately 1‰ of that of NCMCS. Compared with the iCEM method, SSAM not only obtains a reliability index that is very close to NCMCS, but also uses the least CPU time T among all methods, so it has better performance and efficiency. In addition, iCEM N and iCEM LThe obtained reliability indexes have large errors, especially the LOLP index. The reason for this phenomenon is that the normal distribution and lognormal distribution are symmetrical distributions and right-skewed distributions respectively, and only a fixed-shape IS-PDF can be constructed, which makes iCEM N and iCEM L The obtained IS-PDF cannot fully highlight the important areas of the system, so that some samples of small load levels that contribute to the system reliability index cannot be obtained, resulting in a low reliability index.

[0168] In order to more intuitively demonstrate the difference between SSAM and iCEM and reflect the superiority of SSAM, Figure 2 The IS-PDF curves of system load obtained by different methods are plotted in .

[0169] Figure 2 (a) shows the IS-PDF of the system load initially obtained by the iCEM method and the SSAM method. In order to highlight the load level that has a significant contribution to the system failure, the iCEM N and iCEM L The peak of the IS-PDF curve obtained by the method is located at a higher load level (about 0.95 pu) and a large part of its right tail exceeds the annual peak load (1.00 pu), and the system load sample corresponding to this part cannot be extracted in practice. Therefore, iCEM N and iCEM L The method actually uses the truncated IS-PDF in the main sampling stage, such as Figure 3 (b) As shown in Figure 2, the normal distribution and lognormal distribution are symmetric and right-skewed distributions, respectively, and are only suitable for describing important areas with similar shapes to them. Therefore, iCEM N and iCEM L The left tail of the IS-PDF obtained by the method only extends to the area with higher load level (about 0.85 pu). This means that iCEM N and iCEM L The methods can hardly obtain samples with small load levels, especially those with small load levels (such as 0.80-0.85 pu) that have a non-negligible contribution to the system reliability index, so the reliability index obtained by these methods is too small. GMM can construct IS-PDF of any shape, making iCEM G The obtained IS-PDF has a right tail that does not exceed the annual peak load, while its left tail can extend to the area with a smaller load level (about 0.80 pu). G The reliability index error is small. However, the GMM parameter update is complex and involves a multi-step EM algorithm, which leads to iCEMG Compared with iCEM, SSAM can obtain an edge IS-PDF with a reasonable shape without selecting the IS-PDF type of system load and without complicated parameter updating process. Therefore, SSAM can obtain accurate reliability indicators in a shorter time.

[0170] b. IEEE-RTS96 test system with multiple wind farms

[0171] In order to further verify the performance of SSAM in scenarios with correlated random variables, the following modifications are made to IEEE-RTS96:

[0172] ①Increase the load level to 110% of the original level.

[0173] ② Six 200MW wind farms were added on nodes 13, 15, 16, 18, 21 and 23.

[0174] The modified version is named MRTS79-2, and its reliability evaluation results are shown in Table 6.

[0175] Table 6 Reliability evaluation results of MRTS96 system obtained by different methods

[0176]

[0177] As shown in Table 6, as expected, SSAM can obtain accurate reliability indicators with minimal computation. However, since the IS-PDF constructed based on normal distribution and lognormal distribution cannot fully describe the important areas of the system, some samples of small load levels that contribute to the system reliability indicators cannot be obtained, so iCEM is used to construct the IS-PDF based on the normal distribution and lognormal distribution. N and iCEM L The obtained reliability index has a large negative error.

[0178] (3) Analysis results of factors affecting SSAM

[0179] Target value The value of has an important influence on the smoothing parameter σ of the intermediate joint IS-PDF k and the normalization parameter p k The cut-off load level has a direct impact on the estimate of The choice of plays a key role in constructing the initial intermediate joint IS-PDF. In addition, the choice of MH sampling method has a great impact on the effectiveness of sampling from the intermediate joint IS-PDF. The influence of all the above factors will be explored below.

[0180] a. Target value The results of the impact analysis

[0181] To analyze the target value Impact on SSAM performance, using different The SSAM under the value of IEEE-RTS96 is evaluated. As the continuous estimation proceeds, the smoothing parameter σ k and the normalization parameter p k The change curve of Figure 3 The reliability evaluation results are shown in Table 7.

[0182] Table 7 Different The evaluation results of IEEE-RTS96

[0183]

[0184] like Figure 3 As shown, in each given Under this condition, the parameter σ k (Note that σ0 = +∞) and p k As expected, they decrease with the increase of the number of stepwise estimations k, which indicates that as the stepwise estimation proceeds, the estimated intermediate joint IS-PDF g k (u)(0≤k≤K) will be as close as possible to the theoretical optimal IS-PDF. However, as shown in Table 9, The choice of does have a significant impact on the performance of SSAM. Essentially, it is to estimate the parameter σ k and p k Therefore, a very small target value, such as This means very high accuracy requirements. In order to meet such strict requirements, the intermediate joint IS-PDFg k (u) will not be much different from the previous gk-1(u), so SSAM needs to perform K = 5 intermediate joint IS-PDF estimates to achieve an approximate estimate of the theoretical optimal IS-PDF, which greatly increases the computational burden and makes the evaluation take longer. On the contrary, a larger target value, such as This means that the accuracy requirement is lower, so the next estimated g k (u) will be combined with the previous g k-1 (u) is quite different, which means that SSAM only needs to perform K = 2 intermediate joint IS-PDF estimations to converge quickly, thus shortening the evaluation time. However, due to the low accuracy requirement, the estimation of the intermediate joint IS-PDF is not accurate, so the obtained reliability index has a large error. Therefore, as shown in Table 7, in order to achieve a balance between accuracy and efficiency, The value should be kept between 1% and 3%.

[0185] (2) Truncation load level The impact analysis results

[0186] For analysis Impact on SSAM performance, using different The reliability evaluation results of IEEE-RTS96 are shown in Table 8. It is worth noting that the minimum load of the system is 0.34, so This means choosing the original distribution As the initial intermediate joint IS-PDF g0(u).

[0187] Table 8 Different The evaluation results of IEEE-RTS96

[0188]

[0189] As shown in Table 8, when the censored load level is very low, e.g. or SSAM not only requires a longer CPU time T, but also obtains a larger reliability index error. This means that the initial intermediate joint IS-PDFg0(u) is far from the theoretical optimal IS-PDF. It should be noted that IEEE-RTS96 is a highly reliable system because its LOLP is very low, reaching 10 -5 The order of magnitude, therefore, in a sample size of N s = 20000, it is difficult to capture the fault event using g0(u), which causes the subsequent estimated intermediate joint IS-PDF to slowly approach the fault domain, making the estimation time longer. In addition, due to the lack of sufficient fault samples to estimate the subsequent intermediate joint IS-PDF, the final intermediate joint IS-PDFg K (u) may overemphasize the fault area with a smaller fault degree, that is, the fault area where the performance function S(u) value is closer to 0, resulting in a smaller estimated value of the reliability index. On the contrary, if the truncated load level is very high, such as The initial intermediate joint IS-PDF g0(u) will be very close to the theoretical optimal IS-PDF, and due to the high load level, it is easy to obtain samples with very serious faults from g0(u), so the subsequent estimated intermediate joint IS-PDF will quickly approach the fault domain, making the estimation time shorter. However, since too many samples with very serious faults are used to estimate the intermediate joint IS-PDF, the final intermediate joint IS-PDF g0(u) is K(u) may overemphasize the fault area with very serious fault degree, that is, the fault area with smaller performance function S(u) value, which will also lead to smaller estimated value of reliability index. Therefore, it can be seen from Table 8 that It is very reasonable to take the value between 0.6 and 0.8.

[0190] (3) Results of the impact analysis of the MH sampling method

[0191] In order to verify the advantages of conditional MH sampling, the reliability evaluation results of IEEE-RTS96 based on SSAM based on traditional MH sampling method and SSAM based on conditional MH sampling method are compared in Table 9. In order to distinguish, SSAM based on traditional MH sampling method is called SSAM T .also, Figure 4 We also compared the intermediate joint IS-PDF g using the two MH sampling methods. k (u),1≤k≤2, the change of acceptance rate when sampling. k (u), 1≤k≤2, N are retained during sampling c = 0.2N = 4000 seed samples, so Figure 4 The number of samples drawn by medium-condition MH sampling and traditional MH sampling is the number after excluding these seed samples.

[0192] Table 9 is different The evaluation results of IEEE-RTS96

[0193]

[0194]

[0195] like Figure 4 As shown in Figure 2, due to the influence of the curse of dimensionality, in the SSAM based on traditional MH sampling T In this method, the sampling acceptance rate is significantly reduced and almost approaches zero. This phenomenon indicates that the newly generated candidate samples are frequently rejected in the traditional MH sampling process, and thus a large number of duplicate samples are generated. This low acceptance rate may cause the samples generated by the traditional MH sampling method to fail to effectively follow the target distribution g. k (u),1≤k≤2, so that through SSAM T The obtained reliability index is significantly biased compared with the index obtained by the NCMCS method. However, in SSAM using conditional MH sampling, the sampling acceptance rate can be maintained at a reasonable level, thus ensuring the accuracy of the obtained reliability index.

[0196] The above is only an embodiment of the present invention, and the common knowledge such as the known specific technical solutions and / or characteristics in the solution is not described in detail here. It should be pointed out that for those skilled in the art, without departing from the technical solution of the present invention, several modifications and improvements can be made, which should also be regarded as the protection scope of the present invention, and these will not affect the effect of the implementation of the present invention and the practicality of the patent. The scope of protection required by the present invention shall be based on the content of its claims, and the specific implementation methods and other records in the specification can be used to interpret the content of the claims.

Claims

1. A power grid reliability assessment method based on successive smoothing estimation of optimal importance sampling density function, characterized in that: The following steps are involved: Step 1: transform the original mixed variable space of discrete variables and continuous variables into a multivariate standard normal space through the normal transformation method, and on this basis, construct a unified important sampling frame of all random variables in the multivariate standard normal space; Step 2: Based on the important sampling frame, an intermediate smoothed joint IS-PDF sequence is constructed, with the truncated Gaussian distribution as its initial distribution, and the subsequent distributions in the intermediate smoothed joint IS-PDF sequence are estimated step by step in sequence, and the estimated smoothed joint IS-PDF is used for sampling and power grid reliability index evaluation; until the variance coefficient of the reliability index obtained at the Kth estimation meets the preset accuracy, the current power grid reliability index is used as the final evaluation index.

2. The method according to claim 1, characterized in that: Step 1 includes: 1.1) System element state vector x c Each discrete variable in Transformed from formula (2) into the corresponding standard normal variable Among them, δ j =Φ -1 (ξ j ) is the threshold between the fault state and the normal state, ξ j is the original unavailability rate of the jth component, Φ -1 (·) is the inverse function of the standard normal cumulative distribution function, is an indicator function: if but otherwise, 1.2) System load variable x L Transformed into the corresponding standard normal variable u through equal probability transformation L : in, is the raw cumulative distribution function The inverse function of 1.3) Wind speed vector with correlation Transformed into corresponding independent standard normal random vectors by Nataf transformation 3. The method according to claim 2, characterized in that In step 1.3), the process of converting to the corresponding independent standard normal random vector using Nataf transformation includes: First, x is transformed into w Transformed into a zero-mean vector and a correlation matrix C z The normal random vector z is in is the cumulative distribution function The inverse function of Let C x yes The correlation matrix, ρ x (i,j) and ρ z (i,j) are C x and C z The element in the i-th row and j-th column, ρ z (i,j) can be obtained by solving the following formula: Among them, μ i and μ j They are and The mean value, σ i and σ j They are and The standard deviation of is a correlation coefficient C z The bivariate standard normal probability density function of (i,j); After obtaining C z Then, C is calculated by Cholesky decomposition. z The lower triangular matrix L z , that is, C z =L z L z T , then the correlated normal random vector z is composed of independent standard normal random vector u w Indicates that z = L z u w ; Wind speed vector with correlation The corresponding independent standard normal random vector The conversion relationship between them can be expressed by formula (5):

4. The method according to claim 1, characterized in that In step 2, the process of constructing the intermediate smoothed joint IS-PDF sequence and sequentially estimating the subsequent distributions in the intermediate smoothed joint IS-PDF sequence includes: The optimal IS-PDF g in the discrete x space op (x|Q op ) is transformed into the optimal g in the continuous u space op (u|Q op ): in represents the n-dimensional standard normal probability density function, n = D + W + 1, and Represent the D-dimensional, W-dimensional and one-dimensional standard normal probability density functions respectively, is an indicator function: if S(u i )<0, then otherwise, S(u) is the system performance function defined in u-space and can be determined by power flow or optimal power flow analysis of the system state x = T(u) in x-space; Combine formula (6) with the smoothing function Φ(-S(u) / σ K ) to obtain the smoothed optimal zero-variance IS-PDF: Among them, σ k is a sufficiently small smoothing parameter, p k is the normalization parameter, Q k ={σ k ,p k }; Establish the intermediate smoothed joint IS-PDF sequence {g0(u|Q0),g1(u|Q1),...,g k (u|Q k )}, where the truncated Gaussian distribution is used as the initial distribution g0(u|Q0), and the subsequent parameters Q are estimated in turn k ; Parameter Q k The estimation methods include: G k (u) is rewritten as follows: in, The ratio η k The approximate estimate is: Among them, u i is the intermediate smoothed joint IS-PDF g from the previous known k-1 The i-th sample drawn in (u), N is the number of samples drawn; and the unknown parameter σ is obtained by solving the following optimization model k , and then get in is the target value, ||·||2 represents the 2-norm, where yes The coefficient of variance of In getting after, Calculated by the following formula:

5. The method according to claim 4, characterized in that In step 2, the intermediate smoothing joint IS-PDFg k The sampling of (u) adopts a sampling method based on sample resampling and MH sampling technology that avoids the curse of dimensionality, including: Let the IS-PDF g that has been smoothed from the previous intermediate joint k-1 (u) we obtain N iid samples {u1,u2,...,u N }, sample u i The weight ω(u i ) is defined by the following formula: set up is such an indicator function: if u i =u j ,but otherwise, use To express u i Number of occurrences; By performing a resampling operation, N c <N obey g k (u) seed sample: with ω(u i ),1≤i≤N as ω(u i ),1≤i≤N, the probability of occurrence of the total sample {u1,u2,...,u N } randomly select N with replacement c samples; In obtaining N c seed samples {u [1] ,u [2] ,...,u [Nc] }After that, MH sampling is started from these seed samples, while simulating N c The length of the strip is L = N / N c For each Markov chain, from the current known sample u i Get the target distribution g k The next random sample u in (u) i+1 The sampling process is as follows: 1) From the proposed distribution function q(θ|u i ) generates a candidate sample θ; 2) Calculate the ratio: 3) P a =min{1,r(u i ,θ)} Accept θ: Generate a uniformly distributed random number r in the interval (0,1). If r is less than P a , then set u i+1 =θ; otherwise, set u i+1 =u i .

6. The method according to claim 5, characterized in that The multivariate Gaussian probability density function is selected as the proposed distribution function, and the conditional MH sampling method is adopted, the specific contents of which are as follows: Choose the current state vector u i The multivariate Gaussian probability density function conditional on is used as the proposed distribution function: Among them, I is the identity matrix, γ∈[0,1] is a coefficient that can be adaptively adjusted; r(u i ,θ) is simplified to: Use an adaptive method to adjust γ: in, is the newly generated N c The average acceptance probability of samples is the ratio of the number of accepted candidate samples to the total number of candidate samples, γ max ≤1 is the upper limit of γ, a* is the ideal acceptance probability of the MH algorithm, ζ a is an optional parameter, a The selection of should ensure that the change of γ gradually decreases with each adjustment.

7. The method according to claim 6, characterized in that Setting a =(i a +1) -1 , where i a Represents the number of adjustments performed so far.

8. The method according to claim 4, characterized in that In step 2, the truncated normal probability density function is used as the initial intermediate joint IS-PDF, where is the truncated system load distribution, which takes the form shown below. in, is a function that depends on the truncated load level The threshold value of Yes, the parameters of g0(u) are where F L (x L ) is the raw cumulative distribution function of system load.

9. The method according to claim 6, characterized in that In step 2, based on the obtained g k The grid reliability index evaluation of (u) includes: load loss probability LOLP, expected power shortage EENS and its variance coefficient Cov, which are calculated as follows: With each estimation, the maximum number of samples N is preset, and the maximum Markov chain length L = N / N c , set N K is the number of samples drawn, K = k, N K =j·N c , j = 1, we have: Among them, u i is the joint IS-PDFg from the middle k The sample drawn from (u), L c (u i ) represents the system state u determined by optimal load shedding analysis i The load reduction amount, is the likelihood ratio; If the variance coefficient of the LOLP indicator is greater than the maximum variance coefficient β max , then set j = j + 1. If j ≤ L at this time, return to the seed sample extraction step and use the newly generated N c The seed samples update the parameter γ, and then re-sample N k =j·N c The sample is used to calculate the index; if j>L at this time, then the obtained N=L·N c samples as the k (u) and continue with the next round of estimation; If the variance coefficient of the LOLP indicator is less than or equal to the maximum variance coefficient β max , then the intermediate joint IS-PDF g k (u) is taken as the last distribution of the sequence and the calculation result is output.

10. The method according to claim 6, characterized in that The power grid reliability evaluation index is a comprehensive reliability index that makes full use of samples extracted from all intermediate combined IS-PDFs, and its specific calculation includes: Based on g k (u), 0≤k≤K, the expected value and variance of various reliability indicators are estimated by the following unified expressions: Assume κ k for The weight of The definition of the comprehensive reliability index is as follows: Definition The weight κ with the smallest variance k , 0≤k≤K, is the optimal weight, and its specific form is as follows: In the optimal weight κ k (0≤k≤K), Variance for: Coefficient of variance It is estimated by the following formula: Among them, F(u) is the reliability index function in u space. If F(u)=I {S(u)＜0} , It is the LOLP comprehensive reliability index. If F(u)=E c (u)=I {S(u)＜0} L c (u), It represents the comprehensive reliability index of EENS; If the variance coefficient of the LOLP comprehensive reliability index is greater than the maximum variance coefficient β max , then set j = j + 1. If j ≤ L at this time, return to the seed sample extraction step and use the newly generated N c The seed samples update the parameter γ, and then re-sample N k =j·N c The sample is used to calculate the index; if j>L at this time, then the obtained N=L·N c samples as the k (u) and continue with the next round of estimation; If the variance coefficient of the LOLP comprehensive reliability index is less than or equal to the maximum variance coefficient β max , then determine whether the variance coefficient of the EENS index is less than β max If yes, the intermediate joint IS-PDF g k (u) is the last distribution of the sequence and the calculation result is output. If not, set j = j + 1 and re-run N c Markov chain simulation, get N c new samples to update the parameter γ, and then re-sample N k =j·N c Calculate the indicators for the samples.

Images