Renewable energy sources (RES) have strong uncertainties, which significantly increase the risks of power imbalance and load shedding in composite power systems. It is thus necessary to evaluate the operational reliability for guiding economic dispatch and reducing the risks. Current methods cannot meet the requirement for the operational timeliness of reliability evaluations due to the high computational complexity of the optimal power flow (OPF) calculations of massive contingencies. This paper proposes a fully analytical approach to construct fast-to-run analytical functions of reliability indices and avoid reassessments when the load and RES change. The approach consists of uniform design (UD)-based contingency screening and a modified stochastic response surface method (mSRSM). The contingency screening method is used to select critical contingencies while considering the uncertainties. The mSRSM is used to construct the analytical functions of the load shedding to the load and RES generation for the selected contingencies. An analytical function of a smooth virtual variable that maps to the load shedding is established in such a way that, when the load and RES vary, the reliability can be assessed within a very short time rather than using laborious OPF calculations. Case studies illustrate the excellent performance of the proposed method for real-time reliability evaluation.
Longxun Xu, Bo Hu, Changzheng Shao, Kaigui Xie, Congcong Pan, Heng-Ming Tai, Wenyuan Li.
A Fully Analytical Approach for the Real-Time Dynamic Reliability Evaluation of Composite Power Systems with Renewable Energy Sources.
Engineering, 2025, 51(8): 154-168 DOI:10.1016/j.eng.2024.09.023
Renewable energy resources (RES) have been widely utilized in power systems worldwide due to their excellent environmental and economic benefits [1]. The Renewable Energy Market Update[2] by the International Energy Agency (IEA) reveals a continuous growth in global renewable capacity, with an unprecedented surge in 2023. In 2022, China accounted for nearly half of all new renewable power capacity worldwide. As of 2024, the country’s share is set at 55% of the global annual deployment of renewable capacity. However, the increasing RES uncertainty intensifies the difficulty of achieving a real-time balance of power generation and load during operation. For example, in the 2016 South Australia blackout, failures caused the grid to lose around 52% of wind power, resulting in insufficient power supply and ultimately a blackout [3]. In 2021 in Texas, USA, extreme weather forced the shutdown of one-third of the installed capacity, including RES units, and the operators urgently implemented alternating outages, leading to power interruption for five million users [4]. Many blackouts have revealed that it is necessary to focus on the real-time balance between generation and load during operation under the circumstances of RES uncertainty and random failures of power equipment.
The fundamental requirement for power system reliability is to provide sufficient and secure energy. The operators must make dispatching decisions with high reliability to reduce the likelihood and consequences of outages [5]. A reliability evaluation is a complex procedure involving many calculations. The methods for reliability evaluation can be roughly divided into two categories [6]: Monte Carlo simulation (MCS) and state enumeration (SE). MCS evaluates reliability through sampling techniques; its efficiency decreases greatly in high-reliability systems due to the long convergence time. SE acquires reliability indices precisely, since it establishes explicit analytical functions between reliability indices and the availability of components. It usually takes many hours to evaluate the reliability under a single operating state of a power system characterized by its grid topology, generation, and load, which makes operational reliability results outdated. Many methods have been proposed to reduce the computational burden from the perspective of lowering the number of contingencies and simplifying the contingency analysis.
Contingency screening is an effective method to reduce the number of contingencies. An efficient sorting algorithm has been proposed to select contingencies with a high occurrence probability [7]. However, this method may overlook contingencies with low probability but high impact. Algorithms based on sensitivity analysis and overload risk were proposed to determine contingencies with high impact and overload [8], while ignoring the occurrence probability of contingencies. A bilevel optimization model was proposed to identify the highest-risk contingencies and to screen contingencies according to their degree of risk by repeated optimization [9]. In Ref. [10], the impact of a high-order contingency was represented by the sum of the impacts of several low-order contingencies. In Ref. [11], contingencies that caused the same impact were searched and merged into one. A uniform design (UD) approach has been applied to select contingencies, in which reliability indices are obtained through statistical principles [12], [13]. Machine learning methods have also been used to identify critical contingencies [14], [15]. With the extreme faults of power systems as critical contingencies, a reliable spatiotemporal restoration strategy was effectively determined for important loads [16]. A neural-network-based method was proposed in Ref. [17] to accurately extract operation features from historical data.
The challenge in a contingency analysis is how to simplify the optimization process. A redundant-constraint elimination method was proposed to simplify the optimal power flow (OPF) based on the theory of steady-state security regions [18], but it was still necessary to solve an optimization model. Multistate models of wind speed and generation units were built using a Markov chain, and the load shedding was then calculated using universal generating functions [19]. Nevertheless, a reliability analysis for the continuous variation of the operational states could not be performed due to the discrete nature of multistate models. A load-feasible region was proposed to calculate the optimal load shedding in a reliability evaluation as the minimum distance between an operating point and the region [20]. The load and RES uncertainties were considered to be parameters in the OPF [21], [22]. According to multi-parameter linear optimization and the Lagrange multiplier theory, the mapping between the load shedding and these parameters is assumed to be the same when the load and RES generation vary in a small range. However, this mapping is invalid once these values are beyond the range. A similar relationship between reliability and uncertainty was established using a convolutional neural-network-based regression approach [23]. However, this black-box model with poor interpretability suffered from excessive data scale and dimensionality disasters during training.
The above studies have room for improvement:
(1) Most of the contingency screening methods described above did not consider the uncertainty of the load and RES. With a different load and RES, the impact of contingencies differs. Thus, the contingency analyses obtained by the above methods may be invalid. Although machine learning methods can address uncertainties effectively, laborious OPF calculations are still required to evaluate the reliability when the load and RES generation change.
(2) Although the methods focusing on simplifying contingency analyses reduce the computational complexity of reliability evaluations while considering the uncertainties, they still may not satisfy the requirement for timely operation. Except for that in Ref. [23], the existing methods still need to repeatedly perform OPF calculations for reliability evaluation when the load and RES generation change, mainly because they do not provide a direct analytical relationship between load shedding and the uncertainties of the load and RES.
When the generation units and basic network topology remain the same during operation, there should be an explicit function between the reliability indices and load and RES uncertainties under a stochastic contingency. Such a function has not been addressed. This paper establishes a fully analytical model to reveal this analytical function between the reliability indices and uncertainties based on polynomial chaos expansion (PCE).
PCE is often used as an uncertainty quantification model in stochastic equations [24], [25]. The uncertain output is expressed as the weighted sum of orthogonal polynomial basis functions for the random variables [26]. To the best of our knowledge, no study has indicated that the variables in optimization models can also be expressed as a PCE model. This may be due to the non-smoothness of variables in optimization models [27], [28]. With varying parameters in the OPF, the variables restricted by inequality constraints are continuous but may not be differentiable—that is, unsmooth. For unsmooth models, a continuous differentiable PCE model converges slowly and leads to errors at the non-differentiable points. In order to represent the reliability indices as an analytical function with respect to the uncertainties of the load and RES under stochastic contingencies, this paper modifies the stochastic response surface method (SRSM) to establish a PCE-based analytical function for variables in optimization.
The main contributions of this paper are as follows:
(1) An analytical model of reliability indices with respect to the load and RES uncertainties is established to replace the traditional OPF model. Once the load and RES change, the reliability indices are updated by the fast-to-run analytical model instead of laborious OPF calculations. This is practical for evaluating the reliability of an economic dispatch to provide guidance.
(2) A UD-based contingency screening method is proposed to reduce the number of contingencies and the computational complexity. The contingencies selected are critical for different loads and RES. Therefore, the accuracy of the reliability evaluation can be ensured by only analyzing these contingencies, regardless of how the load and RES generation change.
(3) For every selected contingency, an analytical function is revealed for load shedding with respect to the load and RES. Due to its non-smooth nature caused by the constraints of the OPF, the modified SRSM (mSRSM) is proposed to represent a PCE model for a smooth virtual variable. Then, the analytical function of the load shedding is obtained through its map with the virtual variable.
The rest of this paper is organized as follows. The reliability evaluation framework of the proposed method is introduced in Section 2. In Section 3, the proposed analytical approach using the UD-based contingency screening scheme is described. In Section 4, the analytical function of the reliability indices derived from the mSRSM is presented. In Section 5, case studies are performed in which the accuracy and efficiency of the proposed approach are verified and the impact of RES uncertainties on the reliability is analyzed. Conclusions are drawn in Section 6.
2. Problem statement and solution framework
2.1. Overview of composite power system reliability evaluations
Power system reliability is the ability of a power system to continuously provide electricity and energy to loads, based on the given quality standards and required quantities. In this paper, the evaluation object is a composite generating and transmission system without the topology of the distribution system. Fig. 1 shows the procedure for the reliability evaluation of the composite power system. It includes the following three steps:
Step 1: contingency generation. Obtain a contingency and calculate its probability of occurrence. A contingency is a combination of different component outages. For a system of N components, there are 2N contingencies if every component has both an operation state and an outage state. The probability of the contingency s can be calculated as follows:
where Uj is the unavailability of operating component j; J is the set of all operating components; Al is the availability of outage component l; and L is the set of all outage components.
Step 2: contingency analysis. Calculate how much load needs to be curtailed under the contingency obtained in step 1. This is usually achieved through a direct current (DC)-OPF model with the objective function of minimizing load shedding [20], [21], [22], while the generation cost and RES curtailment cost are not considered. When considering the uncertainties of the load and RES generation, the model is reformulated into a stochastic form:
Subject to
where fLS(s) is the system load shedding under the contingency s; cT is the weight of the load at each node, set as an all-ones vector in this paper; Pd is the load demand; T and θ are the power flow of branches and the voltage phase angle of nodes, respectively; θmin and θmax are the lower and upper boundaries of the voltage phase angle, respectively; Tmax is the maximum transmission capacity of branches; X is the branch impedance matrix; A is the node-branch incidence matrix; D, CG, and CRES are the node-branch, node-thermal unit, and node-RES unit incidence matrices, respectively; PG and PRES are the generation of thermal units and RES units, respectively; and are the lower and upper boundaries of the thermal unit generation, respectively; and are the forecast value and forecast errors of RES generation, respectively; and are the forecast value and forecast errors of the load, respectively. Perror is the forecast error vector and uncertain, and is assumed to follow a Gaussian distribution [29], [30]:
where Σ is the correlation matrix of the forecast errors.
The objective function, Eq. (2), is to minimize the load shedding and supply all loads as much as possible. Eqs. (3a), (3b), (3c), (3d), (3e), (3f), (3g) are the constraints of the DC power flow, power balance of the nodes, voltage phase angle, transmission capacity, generation limits (ramping, reserve power, etc.), and load shedding, respectively. In particular, the generation of every RES unit is limited by its actual maximum, which is represented by the forecast value and error in the constraint Eq. (3f). The load shedding at each node is also limited by the actual load represented by the forecast value and error in the constraint Eq. (3g).
Step 3: reliability calculation. Summarize the analysis results and calculate the reliability indices after all the contingencies have been analyzed. The following reliability indices are commonly used: the loss of load probability (LOLP) and the expected demand not supplied (EDNS).
where S is the set of all critical contingencies and I(s) is the following characteristic function.
2.2. Challenges and solution framework
The computational cost of numerous contingencies and their optimization analyses of a composite system reliability evaluation is huge. It is difficult to apply such a reliability evaluation to operations to guide economic dispatch. The reasons are as follows: ① There are too many contingencies to analyze using the OPF model (Eqs. (2), (3a), (3b), (3c), (3d), (3e), (3f), (3g)) with the determined forecast error Perror. ② Once the forecast error Perror changes, the reliability model must be performed again. Therefore, using the traditional reliability evaluation method during operation is time-consuming and wastes a great deal of computing resources.
Fig. 2 illustrates a power system reliability evaluation performed using the traditional method and the proposed method. It can be seen from Fig. 2(a) that the traditional method repeats the reliability evaluation model by means of the OPF to calculate the reliability indices. Once the loads and generations change, the operational reliability results are outdated, which poses a challenge to applying reliability evaluations for system operation. It is necessary to explore the explicit relationship between fLS(s) and Perror. As a result, reliability indices can be acquired in real time via this relationship, instead of repeating the OPF.
This paper proposes a fully analytical approach to establish explicit analytical functions of reliability indices for the load and RES uncertainties. Fig. 2(b) illustrates this reliability evaluation process. A contingency screening method is proposed to reduce the computational complexity. Then, for every contingency, the analytical function of the load shedding fLS with respect to the forecast errors is revealed. Finally, an analytical model of the reliability indices is established at time t to replace the repeating OPF calculations. As a result, the updated reliability indices at any subsequent time tn can be quickly obtained by feeding the new real-time load and RES generation into the analytical model. In this way, the drawbacks of the traditional method are addressed.
In the proposed fully analytical approach, fLS(s) under the contingency s is represented as an explicit analytical function Fs(·) according to Eq. (7) instead of the OPF calculations in Eqs. (2), (3a), (3b), (3c), (3d), (3e), (3f), (3g).
There are two benefits to replacing the OPF with Eq. (7). First, doing so can significantly reduce the computational time. Second, when the load and RES generation change, the reliability indices can be calculated by simply inputting new inputs into the function Eq. (7) for every contingency.
This fully analytical approach analyzes the entire process of reliability evaluation, which can be divided into two parts. First, in contingency generation, UD-based contingency screening is used to screen the critical contingencies with the forecast uncertainty. The contingency set S and the probability p(s) of every contingency are determined. Second, in contingency analysis, the analytical function Fs(·) is established. Finally, the reliability indices are calculated using Eq. (4).
3. UD-based contingency screening
In this section, the load and RES variables are described. Then, UD-based contingency screening is proposed to screen the critical contingency set S, which lowers the computational burden of the reliability evaluation.
3.1. Correlation transformation of the load and RES variables
For each contingency, the OPF model is surrogated by an analytical function that only relates to the loads and RES generation. Here, a stochastic distributionally robust optimization approach for energy systems has been proposed to reduce the load-shedding risks and operation cost [31]. As shown in the upper boundaries of Eqs. (3f), (3g), the actual maximum values of the uncertain loads and RES generation in operation are modeled by the forecast values and the errors.
Under circumstances in which all the forecast errors are independent of each other, matrix Σ is a diagonal matrix. Then, the vector ξ with an independent standard Gaussian distribution can be used to represent Perror:
On the other hand, if correlations between the forecast errors exist, the Cholesky decomposition [24] is used to transform them into independent ones. Then, Perror can also be represented by ξ:
where L is the lower triangular matrix and
Overall, a stochastic variable ξ, which follows an independent standard Gaussian distribution, can be used to represent the forecast error:
3.2. Procedure of UD-based contingency screening
As the OPF model shows, different loads and RES generations will lead to different fLS. The probability of a contingency is constant, but the load shedding with an uncertain load and uncertain RES generation is not. Therefore, whether a contingency is critical depends on its load shedding, considering the entire stochastic space of the load and RES generation. In this section, a scheme for multi-level reliability evaluation through UD-based contingency screening is developed to judge the criticality of a contingency.
Feature levels for the load and RES generation should be used to represent the entire stochastic space. Due to the assumption of a Gaussian distribution for the forecast errors, based on the three-sigma rule of thumb, a seven-level model can be used to discrete these uncertainties—that is, as ±3, ±2, ±1, and 0 standard deviations of the mean (for the numbered experiments 1–7 in ascending order). Regarding the load and RES generation as two factors, a test of two factors with seven levels can be acquired.
Different combinations of the load and RES levels will lead to different fLS and will screen out distinct contingencies. In a sample test of two influencing factors with seven different levels, 27 experiments should be performed to obtain a complete evaluation, according to the theory of permutation and combination, which poses a significant computational burden. However, only seven experiments are required for the UD method. The main idea of the UD method [32] is to scatter its design points uniformly on the experimental domain under a uniformity measure. This is usually performed using an Ur(rm) or Ur*(rm) UD table, where r and m are the number of factors and levels, respectively.
The UD tables of U7(74) and U7*(74) are listed in Table 1, Table 2, Table 3, while their attached tables are in Table 2, Table 3. In Table 1, Table 2, each row represents one experimental arrangement, where the columns are the factors in the test. For a two-factor experiment, only two columns need to be selected from the four. The results of different experimental arrangements are evaluated by the discrepancy [33]. In the attached tables, the best arrangements with the smallest discrepancy are shown. For example, for the UD of two factors with seven levels, the experiments will be arranged in columns I and III of Table 4. Other columns are not used, because the combination of columns I and III has the smallest discrepancy. Moreover, since the discrepancy of U7*(74) is smaller than that of U7(74), U7*(74) is used.
In the UD of two factors with seven levels, seven experiments will be performed. In the first experiment, the elements for the load forecast errors in ξ are equal to level one (i.e., 3 standard deviation), while the elements for the RES generation forecast errors in ξ are set at level five (i.e., –2 standard deviation). Then, the upper boundaries of the constraints in Eqs. (3f), (3g) are determined, and fLS,1(s) for contingency s in the first experiment is calculated. Similarly, fLS,i(s) for contingency s in the ith experiment is calculated, with the load and RES generation forecast error equal to the values of columns I and III in the ith row, respectively.
In this paper, contingency screening under a level of load and RES generation is regarded as an experiment. The UD is used to arrange the experiments to minimize the computational complexity in comparison with that required to complete 27 evaluations. The load and RES generation are calculated using Eq. (11), with columns I and III in U7*(74) selected as the stochastic variable ξ. In every experiment, the risk of a contingency can be determined:
where Ri(s) is the risk of contingency s in the ith experiment.
A threshold ∊ is given to select critical contingencies with a larger risk R than ∊ in every experiment. Then, we combine these selected contingencies from every experiment. To improve the computational efficiency, contingencies that have been selected as critical events in the previous experiments will be skipped to avoid duplicate calculations in subsequent experiments. The procedure for UD-based contingency screening is given below.
Step 1: determine the experimental schemes with the seven-level loads and RES generation according to the UD table. For the ith experiment, the elements representing the loads and RES generation in ξ are equal to the values of columns I and III in the ith row, respectively. Then calculate Perror using Eq. (11).
Step 2: let the set S be an empty set.
Step 3: in every experiment, calculate the OPF model (Eqs. (2), (3a), (3b), (3c), (3d), (3e), (3f), (3g)) and risk (Eq. (12)) for all the contingencies except those in the set S.
Step 4: select the contingencies with a risk larger than the threshold ∊ and incorporate them in the set S.
Step 5: if all the experiments are performed, output the contingencies in the set S as the critical contingencies. Otherwise, return to step 3.
4. Analytical model of the reliability indices via mSRSM
In this section, the analytical function (Eq. (7)) for all the critical contingencies is established based on a PCE model. In an optimization model with varying stochastic parameters, changes in the decision variables are continuous but not differentiable due to the boundaries of the inequality constraints, also known as unsmooth variables. For example, in the OPF model, with uncertain forecast errors, the decision variable Pd, restricted by the constraint in Eq. (3g), is smooth when it varies inside the boundaries but unsmooth when it reaches the boundaries of constraints. As the sum of Pd, the load shedding fLS is also limited between 0 and the sum of the loads . Analytical functions of such unsmooth variables with respect to an uncertain load and RES cannot be directly created by a smooth PCE [27], [28].
First, this section proposes a general mSRSM method to establish the analytical functions of load shedding with respect to the stochastic parameters in an optimization model. It is applied to a contingency analysis in reliability evaluation. The analytical functions of the reliability indices are then developed based on the critical contingencies and their analytical functions (Eq. (7)).
4.1. Smooth process of decision variables
To obtain an analytical expression of the unsmooth decision variable somewhere in the OPF model, a smooth virtual variable unrestricted by the constraints is defined. When the decision variable fluctuates within the boundaries, the virtual variable is equal to it. When the decision variable lies on the boundaries, the virtual variable differs from it. Without loss of generality, we suppose that the load shedding fLS and its virtual variable are the research objects. Fig. 3 illustrates the relationship between them, which can be expressed as follows:
In other words, the load shedding fLS can be regarded as a truncated form of the virtual variable . Hence, once the function of with respect to the stochastic variable ξ is created, the analytical function of fLS can be established using Eq. (13).
4.2. PCE model of the virtual variable
The relationship between the virtual variable and the stochastic parameters is determined via PCE. In the PCE model, the stochastic outputs are represented through random inputs in the form of a weighted sum of orthogonal polynomial basis functions. An analytical function of the virtual variable can be established by the following PCE with finite terms:
where ai and a(i) are PCE coefficients; i1 and i2 are the labels of variables; is the orthogonal polynomial basis function with as input; and Nd is the number of random input variables.
The number of terms in Eq. (14) can be derived as follows:
where NP is the number of terms and coefficients in the PCE model and P is the maximum order of the orthogonal polynomial basis function of PCE.
Then, the coefficients a(i), a(i,i) can be grouped as Coef = [a0, a1, …, ]T, and all the orthogonal polynomial basis functions can be marked from to .
The PCE has two key components: the orthogonal polynomial basis functions and the coefficients. The functions can be divided into one-dimensional (1D) and multi-dimensional input basis functions. The 1D orthogonal polynomial basis functions in this paper are the Hermite polynomials, since all variables follow a standard Gaussian distribution [26]. The orthogonal polynomial basis functions with multi-dimensional input variables can be acquired by the tensor product based on the 1D functions:
where ξ1, ξ2, and ξNd are the Nd elements of the vector ξ.
According to the properties of polynomial functions, each term in Eq. (14) is continuous and differentiable. Consequently, the PCE model is smooth. Regardless of how large an order P of the PCE model is used to fit the load shedding directly, a systematic error always exists where the load shedding is unsmooth. Therefore, it is better to use the virtual variable in the intermediate process.
Traditionally, the coefficients in PCE can be derived by the SRSM [24]. The core of the SRSM is to obtain M fixed samples ξ1, ξ2, …, ξM, also known as collocation points, and calculate their load shedding samples , , …, using Eqs. (2), (3a), (3b), (3c), (3d), (3e), (3f), (3g). Based on these pairs between the collocation points and the load shedding samples, the coefficients can be obtained via least squares regression.
However, in this paper, only a portion of the collocation points can be used to solve the PCE coefficients of the virtual variable. The virtual variable is determined by the load shedding between the lower and upper boundaries, and it is undefined when the load shedding lies on the boundaries. Therefore, only some valid collocation points, which make the load shedding vary between the boundaries, are useful for virtual variables.
Fig. 4 provides a better illustration. In Fig. 4(a), according to the SRSM, M collocation points are given, and the load-shedding samples of the optimization model are calculated. It can be seen that some samples fall on the lower and upper boundaries, which leads to a decrease in the valid samples. Moreover, since at least one of the parameters X, D, CG, and CRES in Eqs. (2), (3a), (3b), (3c), (3d), (3e), (3f), (3g) is different in two contingencies, the load shedding may be different. The valid collocation points obtained by using fixed collocation points are also different. Consequently, it is difficult to solve the PCE coefficients in all contingencies with the same fixed collocation points.
An iterative Latin hypercube sampling (LHS) method is proposed to obtain valid collocation points. LHS is applied to get collocation points—instead of fixed collocation points—that can cover the entire stochastic space regardless of the different OPF of the contingencies. The undetermined coefficients require at least NP valid collocation points due to the NP terms in Eq. (14). Therefore, the number of intervals, the only parameter in LHS, should be no less than NP. Similarly, to the SRSM, the valid collocation points obtained from a single sampling are not enough. An iterative sampling strategy is proposed to keep sampling the collocation points until the cumulative number NV of valid collocation points is greater than NP. As demonstrated in Fig. 4(b), after sampling the collocation points ξ1,1, ξ1,2, …, ξ1,M by means of LHS in the first iteration, we calculate all the load-shedding samples and count the number of valid collocation points. If it is less than NP, we sample M collocation points again. At the uth sampling, when the cumulative number NV is larger than NP, we stop the sampling and establish the PCE model of the virtual variable for this contingency.
After NV collocation points are obtained, the PCE coefficients of the virtual variable are solved using Eq. (17).
where
and where , …, are the NV load shedding samples corresponding to the NV valid collocation points; , …, are the orthogonal polynomial basis functions with the valid collocation point as input, k = 1, 2, …, NV.
4.3. Analytical function of reliability indices
In this paper, the relationship between fLS and its virtual variable is described as Eq. (13), where is formulated as the PCE model (14). Once the load shedding fLS for all the critical contingencies is calculated, reliability indices are obtained using Eq. (4). However, the OPF model varies for different contingencies. In some contingencies, it is easy to obtain valid collocation points, while the opposite is true in others. Therefore, the number of sampling iterations is also different. It is unrealistic to sample unlimitedly in order to obtain sufficient valid collocation points in contingencies where it is difficult to obtain them.
To improve the computational efficiency, two upper iteration numbers, n1 and n2, are set to stop the sampling. All the contingencies can be partitioned into three groups, sets E, F, and G, according to the two thresholds. A contingency belongs to set E if all collocation points sampled via LHS in this contingency remain invalid in n1 sampling iterations. If NV(s) is smaller than NP after n2 iterations in the contingency s, then the contingency s belongs to set F. Other contingencies for which NV(s) is larger than NP belong in the set G.
As the NV of the contingencies in sets E or F is insufficient, such that the coefficients of the PCE model cannot be obtained, the load-shedding events that occur in these contingencies are ignored during modeling. As a result, the PCE outputs in these contingencies are roughly estimated as constant at 0. On the other hand, enough collocation points can be sampled for contingencies in set G. Therefore, their load shedding is described as the PCE functions.
The procedure of the mSRSM to establish an analytical function of the load shedding fLS(s) is described as follows:
Step 1: suppose a virtual response according to Eq. (13), and let NV(s) = 0.
Step 2: in one sampling, select M collocation points via LHS.
Step 3: determine the load shedding fLS corresponding to the collocation points generated in step 2. Calculate the forecast error of the load and RES generation using Eq. (11) and the load-shedding samples fLS(s) using Eqs. (2), (3a), (3b), (3c), (3d), (3e), (3f), (3g).
Step 4: determine the number of these load-shedding samples that vary within 0 and . Update NV(s) with this number and save these collocation points.
Step 5: stop sampling if either one of the following conditions is satisfied; otherwise, return to step 2: ① After iterating n1 times, all the true responses so far are on the boundaries—that is, NV(s) = 0. ② NV(s) < NP after iterating n2 times, where n1 < n2.
Step 6: solve the coefficients of the virtual response using Eq. (17), with all the valid collocation points and their load-shedding samples.
Step 7: obtain the analytical function, Eq. (7), by means of Eq. (13).
Finally, based on the analytical expressions in Eqs. (1), (4), (13), and (14), the fully analytical model of the reliability indices is established after completing an analytical function of fLS for all contingencies:
where sgn(·) is the signum function. p(s) is calculated in Eq. (1), fLS(s) is calculated in Eq. (13), and is calculated in Eqs. (14), (17).
A flowchart of the fully analytical approach is provided in Fig. 5. In the contingency generation, the critical contingency set S is selected, and then the probability p(s) of every contingency is determined. In the contingency analysis, the analytical functions (Eq. (7)) for critical contingencies are formed by the mSRSM. Then, the fully analytical model (Eq. (14), (17)) of the reliability indices is obtained.
Since the analyses of each contingency are independent from each other, parallel computing methods can be applied to improve the computational efficiency during modeling.
5. Case study
A performance evaluation of the proposed fully analytical approach was carried out on three systems: the Roy Billinton Test System (RBTS) [34], Reliability Test System (RTS)-79 [35], and RTS-96 [36]. The algorithm was tested on a desktop with a 3.0-GHz Intel Core i5-12500H processor and 16 GB memory by means of MATLAB R2019a.
We suppose that all the loads follow a Gaussian distribution, with the means equal to the peak loads and the standard deviations equal to 5% of the means. Some units are replaced by RES units, whose generation also follows a Gaussian distribution, with the means equal to the original maximum generation and the standard deviations equal to 15% of the means. Parallel computation with six cores was adopted in the tests on RTS-79 and RTS-96. Many samples with different random loads and RES generation inputs by means of SE and the proposed method were evaluated.
The performance metrics are the mean absolute percentage errors (MAPEs) for the indices LOLP and EDNS, defined as follows:
where NS is the number of samples in the case studies and SE is state enumeration.
5.1. Accuracy of the mSRSM
The RBTS was modified to test the accuracy of the mSRSM. The RBTS is composed of six buses, nine transmission lines, and 11 units. The total installed capacity is 240 MW, and the peak load is 185 MW. In the modified RBTS, two 20-MW units on buses #1 and #2 are replaced by RES units. The contingencies are enumerated to the fourth contingency order, where the reliability evaluation by means of SE takes 15.8 s, on average.
Table 5 shows the results of applying the SRSM and mSRSM to establish the analytical functions of the reliability indices. The maximum orders of the PCE in the SRSM are P = 1 and 3. The number of collocation points used in the test is twice the dimension of the PCE [27]. The parameters in the mSRSM are M = NP, n1 = 1, and n2 = 10. It can be seen from Table 5 that the SRSM cannot provide accurate LOLP and EDNS, because the MAPE is very large, even more for a large P. This may be caused by the overfitting of the linear optimization model with the nonlinear PCE. The MAPEs for the mSRSM are much smaller. Thus, the mSRSM is used to construct the analytical models for the reliability evaluation.
Table 6 shows the evaluation results from the fully analytical approach using the mSRSM, without considering the UD-based contingency screening. The parameters used include M, the number of sampling cases; n1, the upper bound of the iteration for ① in step 5 of the mSRSM procedure; and n2 for ② in step 5. Table 6 illustrates the effects of these three parameters on the performance of the MAPE and computation time. Card(E) denotes the number of contingencies in set E, and Card(F) and Card(G) are defined similarly. The larger Card(G) is, the more accurate the model is and the smaller the MAPE of the LOLP and EDNS.
Fig. 6 shows the quantile–quantile (Q–Q) diagrams of LOLP and EDNS for the proposed method and the SE method. The parameters M, n1, and n2 are set to NP, 2, and 30, respectively. The Q–Q diagram plots the quantiles of the probability density functions (PDFs) by the two methods and tests the similarity by comparing their similarities to the line y = x. The two PDFs are more similar when the Q–Q diagram is close to the line y = x. It can be inferred from Fig. 6 that the mSRSM exhibits great accuracy, since the analytical models for most contingencies are accurate. Fig. 7 shows the histograms of LOLP and EDNS. Based on the histograms, the operators can obtain the varying ranges and extreme values of two indices with the given load and RES uncertainties during dispatch.
The main error of the mSRSM is caused by the constant 0 for contingencies in sets E and F. As shown in Table 6, the error of the mSRSM decreases with the increase in the cardinality of set G. It can also be observed from Table 6 that MAPELOLP is always larger than MAPEEDNS. According to Eq. (4), the reliability index EDNS is a continuous function of fLS(s), while LOLP is a piecewise function. As a result, any small error in the vicinity of segment places may make the LOLP leap.
During the modeling, the computation time is influenced by three parameters: M, n1, and n2. However, this is not the case for the time of the reliability evaluation. Moreover, since the evaluation is performed using the fully analytical functions, the time required to evaluate the reliability of the composite power system is much less than the modeling time.
A sensitivity test of the parameters in the mSRSM was also conducted. Fig. 8 depicts two MAPEs and the cardinality of set G with different M, n1, and n2 in Table 6. It can be seen from Fig. 8 that MAPELOLP and MAPEEDNS decrease as M, n1, and n2 increase independently, while the modeling time varies contrary to the two MAPEs. The larger n1 and n2 mean that more iterations are needed to obtain collocation points for approximating the virtual variables in some contingencies. The number of collocation points attained in every iteration grows as M increases. As a result, the modeling time rises. Fewer contingencies break the loop in step 5 of the mSRSM, and their analytical models of load shedding are established rather than using the constant 0, which downscales the sets E and F and improves the cardinality of the set G. Therefore, the accuracy of the mSRSM improves as the MAPEs are lowered.
5.2. Benefit of UD-based contingency screening
We now investigate the benefit of the UD-based contingency screening on the reliability evaluation. The test system is a modified RTS-79, which has 24 buses, 38 transmission lines, and 32 units. The units on bus #18 (400 MW), bus #13 ((2 × 197) MW), bus #15 (12 MW), bus #2 (20 MW), and bus #7 (100 MW) are replaced by RES units. The total installed generation capacity is 3405 MW and the peak load is 2850 MW. The contingencies in this test are enumerated to the third order of the contingencies. Each sample evaluated by SE takes 187 s, on average.
A reliability evaluation using the proposed method was conducted with and without the UD-based contingency screening. We let parameters M = NP, n1 = 2, and n2 = 30. The threshold ∊ is determined by the product of a ratio ρ, the median probability of the highest-order contingencies and the total system load [11]. By adjusting the ratio ρ, the screening results of the UD-based contingency screening can be adaptively controlled. Table 7 displays the computation time and evaluation errors. The “—” symbol denotes the reliability evaluation without the UD-based contingency screening. It can be seen from Table 7 that the proposed method provides more accurate results with a smaller ρ. When the UD-based contingency screening is not applied, the proposed method offers the best accuracy but the worst computational efficiency. For the case where ρ = 0.5, the maximum MAPE is still smaller than that in Ref. [37].
Table 8 summarizes the number of contingencies with different ρ. It shows that the UD-based contingency screening helps to lower the number of selected contingencies and the cardinality. It can also be seen from Table 8 that the contingencies in set G, which affect the load shedding, only take up a small proportion of the entire contingency space. As ρ becomes smaller, more contingencies in set G are chosen to participate in the modeling procedure. The scale of the contingencies constructing the analytical models of the load shedding grows. Meanwhile, the number of contingencies whose load shedding is marked as constant at 0 decreases. Consequently, the proposed method enables the reduction of MAPELOLP and MAPEEDNS.
The Q–Q diagrams of LOLP and EDNS using the proposed method and the SE are depicted in Fig. 9. For the same sample, the LOLP and EDNS values calculated by the proposed method increase as ρ decreases, which is near the line y = x. As ρ declines further, the LOLP and EDNS using the proposed method are closer to those obtained using SE. The histograms of the reliability depicted in Fig. 10 also show that, when ρ gets smaller, the height of the bar when using the proposed method is closer to that when using SE. Since the contingencies that are not selected for establishing models possess small risks and contribute less to the reliability indices, their effects on the accuracy are insignificant.
With a small ρ, the fully analytical approach takes more time to model because the load sheddings of more contingencies are established as analytical functions by the proposed method. In contrast, unlike the parameters in the mSRSM, a larger ρ also results in an increase in the evaluation time because of the increasing number of analytical models for contingencies. The results show that the UD-based contingency screening method significantly reduces the modeling time and offers acceptable accuracy.
5.3. Impact of RES uncertainties on reliability
The modified RTS-96 was used to evaluate the impact of RES uncertainties on reliability. It is a system made up of three modified RTS-79s, where the units in the same location in the three modified RTS-79s are replaced by RES units. RTS-96 is a large system containing 73 buses, 120 transmission lines, and 96 units. This system usually has a total installed generation capacity of 10 215 MW and a peak load of 8 550 MW. In the test, the peak load of the modified RTS-96 was up to 9405 MW. The comprehensive outages and unit outages were enumerated to the third contingency order, while the branch outages were only enumerated to the first contingency order. There were a total of 706 337 contingencies, and each sample took about 1.5 h.
Four cases were tested regarding the RES forecast errors and the standard deviations of the RES generation. Table 9 shows the results with the mSRSM parameters ρ = 0.1, M = NP, n1 = 2, and n2 = 30. It is notable that each sample evaluated by SE would take 1.34 h on average. It can be seen from Table 9 that the proposed method exhibits very good performance in terms of both accuracy and evaluation time. The contingencies selected by the UD-based contingency screening account for less than 13% of all contingencies. As the RES forecast errors increase, the errors of the proposed method only increase slightly. This increase is mainly caused by a few extreme samples, as shown in Fig. 11. It can be said that the indices calculated by the mSRSM match well with those calculated by SE under most conditions.
Due to the acceptable accuracy of the fully analytical approach, all the reliability evaluation results were calculated using the proposed fully analytical model. In the test, 10 000 loads and RES generations were sampled. The results for the statistical features (i.e., mean, variance, and skewness) of the reliability indices with different RES forecast errors of LOLP are shown in Table 10 and those of EDNS in Table 11. These statistics reflect the features of the system reliability under different uncertainties. It can be seen from Tables 9 and 10 that the larger the RES forecast errors are, the larger the statistical features are. Details are described below.
(1) Mean. A large mean indicates that larger RES forecast errors will lead to a greater risk of load shedding. If more RES generation occurs in the power system, the generation of thermal units will be decreased, but there will be little impact on reliability. However, less RES generation will lead to more extreme scenarios in which the system’s generation is far less than the load, lowering the system reliability.
(2) Variance. The strong uncertainty of the input RES generation causes large fluctuation (i.e., variance) of the sample outputs.
(3) Skewness. A larger skewness indicates that the distribution of samples inclines more toward the right tail, which is caused by the increase in contingencies with a high risk of load shedding.
Consequently, since the accuracy is verified in Table 9 with different RES forecast errors, it can be concluded that larger RES forecast errors lead to a greater variation in reliability indices. On the other hand, although the forecast errors increase over time, the proposed method still has high accuracy. Therefore, for large-scale systems with a long modeling time, establishing a fully analytical model based on earlier forecast results is feasible.
6. Conclusions
This paper established fully analytical functions of the reliability indices of the load and RES uncertainties to satisfy the requirement of fast operational reliability evaluation in a composite power system. UD-based contingency screening was applied to select critical contingencies that were applicable for various load and RES scenarios. An mSRSM was developed to construct the analytical functions of load shedding for contingencies. In the mSRSM, a smooth virtual variable is used to construct the load shedding so that the PCE method can be applied to obtain the analytical functions of the reliability indices. The fully analytical model enables the reliability evaluation to be performed by means of simple functions when the load and RES change, rather than repeated OPF calculations.
The case study results revealed that the proposed method performed well in terms of accuracy and computational efficiency and met the requirements for the real-time calculation of operational reliability evaluations. Moreover, the fully analytical model obtained the statistical features of the reliability indices for different RES scenarios.
In order to reduce the modeling time, a classifier can be applied to the PCE modeling of virtual variables to determine whether the collocation point is valid, while OPF is only used to calculate the amount of load shedding for the valid collocation points. Parallel computing can also be used to further improve the model’s efficiency. In addition, since the proposed method places no limitations on the OPF model, this method can be extended to the analytical expression of other optimization models.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
This work was supported by the Joint Research Fund in Smart Grid under cooperative agreement between the National Natural Science Foundation of China (NSFC) and the State Grid Cooperation of China (SGCC; U1966601), and the Fundamental Research Funds for the Central Universities of China (2023CDJYXTD-004).
AhsanF, DanaNH, SarkerSK, LiL, MuyeenSM, AliMF.Data-driven next-generation smart grid towards sustainable energy evolution: techniques and technology review.Protect Control Mod Power Syst2023; 8(3):1-42.
[2]
Renewable energy market update—June 2023.Report. Paris: International Energy Agency; 2023.
[3]
YanR, MasoodNA, SahaTK, BaiF, GuH.The anatomy of the 2016 South Australia blackout: a catastrophic event in a high renewable network.IEEE Trans Power Syst2018; 33(5):5374-5388.
[4]
ZhangG, ZhongH, TanZ, ChengT, XiaQ, KangQ.Texas electric power crisis of 2021 warns of a new blackout mechanism.CSEE J Power Energy Syst2022; 8(1):1-9.
[5]
PanC, HuB, ShaoC, XuL, XieK, WangY.Reliability-constrained economic dispatch with analytical formulation of operational risk evaluation.IEEE Trans Power Syst2024; 39(2):4422-4436.
[6]
LiW.Risk assessment of power systems, models, methods, and applications. Wiley, Hoboken (2014)
[7]
LiuH, SunY, WangP, ChengL, GoelL.A novel state selection technique for power system reliability evaluation.Electr Power Syst Res2008; 78(6):1019-1027.
[8]
DavisCM, OverbyeTJ.Multiple element contingency screening.IEEE Trans Power Syst2011; 26(3):1294-1301.
[9]
DingT, LiC, YanC, LiF, BieZ.A bilevel optimization model for risk assessment and contingency ranking in transmission system reliability evaluation.IEEE Trans Power Syst2017; 32(5):3803-3813.
[10]
LeiY, SunY, HouK, ZhangP, ZhuL, YangX.Impact increment based hybrid reliability assessment method for transmission systems.CSEE J Power Energy Syst2022; 8(1):317-328.
[11]
WangL, HuB, XieK, NiuT, LiW, LiaoQ, et al.Screening model of incremental risk events for reliability analysis of transmission system.Int J Electr Power Energy Syst2020; 120:105995.
[12]
XieK, HuangY, HuB, TaiHM, WangL, LiaoQ.Reliability evaluation of bulk power systems using the uniform design technique.IET Gener Transm Distrib2020; 14(3):400-407.
[13]
ShaoC, ZhaoS, FengQ, HuB, XieK.Reliability evaluation of multi-area integrated electricity–gas systems based on the improved uniform design.IEEE Trans Power Syst2024; 39(4):5932-5945.
[14]
PindoriyaNM, JirutitijaroenP, SrinivasanD, SinghC.Composite reliability evaluation using Monte Carlo simulation and least squares support vector classifier.IEEE Trans Power Syst2011; 26(4):2483-2490.
[15]
UrgunD, SinghC.A hybrid Monte Carlo simulation and multi label classification method for composite system reliability evaluation.IEEE Trans Power Syst2019; 34(2):908-917.
[16]
JianJ, ZhaoJ, JiH, BaiL, XuJ, LiP.Supply restoration of data centers in flexible distribution networks with spatial–temporal regulation.IEEE Trans Smart Grid2024; 15(1):340-354.
[17]
ZhaoJ, ZhangZ, YuH, JiH, LiP, XiW.Cloud-edge collaboration-based local voltage control for DGs with privacy preservation.IEEE Trans Ind Inf2023; 19(1):98-108.
[18]
HuaB, BieZ, LiuC, LiG, WangX.Eliminating redundant line flow constraints in composite system reliability evaluation.IEEE Trans Power Syst2013; 28(3):3490-3498.
[19]
DingY, SinghC, GoelL, WangP, et al.Short-term and medium-term reliability evaluation for power systems with high penetration of wind power.IEEE Trans Sustain Energy2014; 5(3):896-906.
[20]
LiX, XieK, ShaoC, HuB.A region-based approach for the operational reliability evaluation of power systems with renewable energy integration.IEEE Trans Power Syst2024; 39(2):3389-3400.
[21]
YongP, ZhangN, KangC, XiaQ, LuD.MPLP-based fast power system reliability evaluation using transmission line status dictionary.IEEE Trans Power Syst2019; 34(2):1630-1640.
[22]
LiuZ, HouK, JiaH, ZhaoJ, WangD, MuY.A Lagrange multiplier based state enumeration reliability assessment for power systems with multiple types of loads and renewable generations.IEEE Trans Power Syst2021; 36(4):3260-3270.
[23]
KamruzzamanM, BhusalN, BenidrisM.A convolutional neural network-based approach to composite power system reliability evaluation.Int J Electr Power Energy Syst2021; 135:107468.
[24]
RenZ, LiW, BillintonR, YanW.Probabilistic power flow analysis based on the stochastic response surface method.IEEE Trans Power Syst2016; 31(3):2307-2315.
[25]
XuY, KorkaliM, MiliL, ChenX, MinL.Risk assessment of rare events in probabilistic power flow via hybrid multi-surrogate method.IEEE Trans Smart Grid2020; 11(2):1593-1603.
ZhuJ, LiG, GuoY, ChenJ, LiuH, LuoY.Real-time risk-averse dispatch of an integrated electricity and natural gas system via conditional value-at-risk-based lookup-table approximate dynamic programming.Protect Control Mod Power Syst2024; 9(2):47-60.
[30]
NajibiF, ApostolopoulouD, AlonsoE.Enhanced performance Gaussian process regression for probabilistic short-term solar output forecast.Int J Electr Power Energy Syst2022; 15(6):106916.
[31]
ZhouY, LiX, HanH, WeiZ, ZangH, SunG, et al.Resilience-oriented planning of integrated electricity and heat systems: a stochastic distributionally robust optimization approach.Appl Energy2024; 353(A):122053.
[32]
FangKT, LiuMQ, QinH, ZhouYD.Theory and application of uniform experimental designs. Springer, Singapore (2018)
[33]
HickernellFJ.Goodness-of-fit statistics, discrepancies and robust designs.Stat Probab Let1999; 44(1):73-78.
[34]
BillintonR, KumarS, ChowdhuryN, ChuK, DebnathK, GoelL.A reliability test system for educational purposes-basic data.IEEE Trans Power Syst1989; 4(3):1238-1244.
[35]
Probability Methods Subcommittee. IEEE reliability test system. IEEE Trans Power Appar Syst 1979; PAS-98(6):2047–54.
[36]
GriggC, WongP, AlbrechtP, AllanR, BhavarajuM, BillintonR.The IEEE reliability test system—1996. A report prepared by the reliability test system task force of the application of probability methods subcommittee.IEEE Trans Power Syst1996; 14(3):1010-1018.
[37]
ZhuY, SinghC.Assessing bulk power system reliability by end-to-end line maintenance-aware learning.IEEE Access2023; 11:49639-49649.