Abstract

The increased integration of renewable energies (REs) raised the uncertainties of power systems and has changed the approach to dealing with power system challenges. Hence, the uncertain nature of all the power system variables needs to be considered while dealing with the optimal planning and operation of modern power systems. This paper presents a probabilistic optimal active and reactive power dispatch (POARPD) based on the point estimate method (PEM), considering the uncertainties associated with load variation and wind power generation. In the POARPD, the deterministic optimal active and reactive power dispatch (OARPD) is performed in two stages, which gives a deterministic two-stage OARPD (TSOARPD). The objectives of TSOARPD are the operating cost (OC) minimization in stage 1 and voltage stability (VS) maximization in stage 2, whereas the VS is improved by maximizing the system’s reactive power reserve (RPR). In this paper, instead of using multiobjective optimal power flow, this TSOARPD is used to give more importance to VS when the system is substantially loaded. The POARPD problem is solved using PEM for modified IEEE-9 bus and standard IEEE-30 bus test systems by considering the correlation between the loads. The results are compared with Monte Carlo simulation (MCS). While solving POARPD, the voltage-dependent load model is used to account for the real-time voltage dependency of power system loads. This paper discusses the detailed procedure of solving POARPD by considering correlation and the increased nonlinearities by giving more importance to VS when the system is heavily loaded.

1. Introduction

Conventional techniques to solve power system problems use deterministic approaches that ignore system variables’ uncertainty. But most real-time power system variables and characteristics, such as loads and the parameters of transmission lines, are unpredictable. In addition to these uncertainties, modern power systems incorporated with renewable energies gain more unpredictability, thereby increasing the complexity of planning and operation of modern power systems. Thus, modern power systems with increased uncertainties need approaches that consider different power system states, while conventional methods consider only one [1]. As a result, it is critical to think of the power-flow problem as a probabilistic problem. The OPF is the most commonly used tool for the planning and operation of power systems, and it optimizes the power system’s objective function [2]. Thus, to obtain more realistic ranges of the results of OPF, the uncertainties related to power system variables need to be incorporated into the deterministic OPF, which gives the probabilistic optimal power flow. Thus, POPF utilizes stochastic information of the input-random variables and gives the optimal control variables’ settings and the system’s optimal output variables in the form of random variables. This functionality of POPF makes it simple to identify and investigate the possible ranges and violations of deterministic OPF results.

The first step in solving POPF is to select a suitable and efficient method of handling uncertainty. In the literature for solving the POPF problem, numerical or analytical methods are extensively used [36]. Numerical methods such as MCS and its derivatives have been widely used for dealing with power system problems involving uncertainties [2, 7]. Among all the probabilistic methods, MCS is considered the benchmark method due to its high computational accuracy, but its major drawback is its huge computation burden. To overcome this drawback, approximate methods such as the point estimate method are developed for probabilistic analysis [8]. The PEM approach is less computationally expensive and close to MCS in terms of accuracy. The basic scheme of PEM is 2 PEM which uses two estimated points of inputs, and higher schemes up to seven points (7 PEM) have been proposed in the literature [9]. The PEM, along with Nataf transformation, is used for dealing with uncertainties related to wind turbine generating systems while solving power system problems in [5]. And discrete PEM is used for solving probabilistic load flow with WTGS in [10]. The PEM gives statistical moments of output variables; from these output moments, the output variables’ distributions will be computed using suitable approximation series such as Cornish–Fisher [3] and Gram–Charlier [4, 11].

Once the suitable solution method is selected, then the next step is to supply uncertain input variables to that solution method. In the literature, authors have used various probability distributions to model the uncertainties associated with the power system’s input random variables. The power system’s active and reactive loads are often modeled by normal distribution [1214]. For modeling uncertainties related to the prominent RE sources such as wind speed and solar irradiations, probability distributions other than Gaussian are often used in the literature. The Weibull distribution is used extensively for the modeling of wind speed [1315], and few authors have used other distributions, such as the log-normal distribution also [16]. And the majority of studies employed the beta distribution to model solar irradiation uncertainty [13, 17, 18]. In stochastic analysis, the correlation among the input random variables plays an important role, and the consideration of correlation improves the accuracy of the output variable estimation. In literature, the correlation is handled using the combination of Nataf transformation and Cholesky decomposition [16, 19]. Most papers used wind speed as the input variable to PEM directly, but the authors of [20] discuss the importance of considering the wind generator’s output power as the input variable to PLF instead of the wind speed. As only the wind speeds which are in the range of cut-in and cut-out speeds of the wind turbine are eligible to generate power, the consideration of wind generator output power as the input variable to PEM is more accurate.

As the power system loads are considered random variables, the next step is to include load modeling in PLF. Load modeling is essential in solving power system problems such as OPF, but it has yet to be given priority in PLF or POPF in most papers. In a real-time power system, the ZIP model is used to represent the loads. The foundation of the ZIP model is the division of the load into constant impedance (Z), constant current (I), and constant power (P) loads, all of which are connected in parallel [21]. Therefore, instead of taking the load power as the random variable, the ZIP load needs to be considered as the random variable in PLF or POPF [22]. This consideration of ZIP load as the random variable increases the nonlinearity and affects the computational time and accuracy of PLF (or POPF). Even if this consideration affects the speed and accuracy of PLF, the results will be more realistic.

In this paper, the POPF (here, POARPD in this paper) problem is solved considering the correlated load and wind uncertainties by using PEM and MCS, where wind generators’ output powers are taken as input random variables instead of wind speeds. In the POPF, the deterministic OPF is solved in two stages (TSOARPD) to give more importance to VS than OC when the system is heavily loaded. This TSOARPD is used to avoid the use of multiobjective optimal power flow, which may require proper settings of goals or weights for all the objectives while deciding the objectives’ importance. In the deterministic TSOARPD, stage 1 performs optimal active power dispatch for OC minimization, and stage 2 performs optimal reactive power dispatch for VS maximization. In the literature, there are many ways to indicate the voltage stability margin of power systems [23], such as the smallest singular value of the load-flow Jacobian matrix [24, 25], based on the nose point of power-voltage curves [26], and a static voltage stability indicator denoted by L-index [27, 28], and RPR of the system [2933]. In this paper, maximization of the RPR is taken as an objective function (in stage 2 of TSOARPD) for maximizing VS of the system, whereas the RPR of a system indicates the amount of reactive power generating capacity available at the sources (at the present loading point) that the system can further utilize. This paper uses the interior-point method for solving the deterministic TSOARPD optimization problem. The voltage-dependent load model is considered to incorporate the real-time voltage dependency of power system loads. This paper considers handling correlation and the voltage-dependent loads in the POPF problem using PEM, whereas the existing papers considered either correlation or voltage-dependent loads in POPF. As discussed previously, this paper solves POARPD with proper handling of correlation and dealing with the increased amount of nonlinearity in the system by prioritizing VS (maintaining OC close to its optimal value) when the system is heavily loaded.

Contributions to this work are as follows:(i)PEM-based POARPD considering the correlation(ii)Nataf transformation, along with Cholesky decomposition, is used to handle correlation(iii)Gauss–Hermite quadrature is used to solve the double integral equation used in the Nataf transformation(iv)Interior-point method (for optimization) is applied to solve deterministic TSOARPD(v)Two types of WTGS are considered while solving POARPD(vi)Voltage-dependent (composite) load model is also used

The problem formulation, scenarios considered, details of the test systems used, and observations of the work are presented in the following sections.

2. Problem Formulation

2.1. Deterministic OPF

The deterministic OPF problem for minimizing a power system’s objective function [34] can be mathematically expressed as follows:where f represents the objective function such as operating cost, power loss; x and u are the set of dependent and independent variables of the power system, respectively; and represents the equality constraints (such as power balancing equations), whereas denotes the inequality constraints for all the power system variables [35].

2.2. Deterministic TSOARPD

As the two objectives, i.e., OC minimization and VS maximization, mainly depend on the active and reactive control variables’ settings, respectively, they can be achieved one after the other. For each sample of input random variables (supplied by MCS or PEM), the steps for the TSOARPD are given as follows:Stage 1: Minimize OC by using OAPD and obtain the optimal settings of active power generations of all the generator buses.Stage-2: Fix the active power generations of the generator buses (except slack bus) at their optimal values obtained in stage 1. Now, maximize VS of the system through the maximization of the system’s RPR by using ORPD (using optimal settings for the reactive control variables such as voltage magnitudes of the generator buses, switchable shunt capacitors, and tap-changing transformers).

This successive application of OAPD and ORPD does not need the appropriate settings of weights or goals (those required in MOOPF) to give more importance to VS by keeping the OC close to its optimal value when the system is substantially loaded. With the minimum number of control variables in each stage, and as the complexity involved is minor (as only one objective needs to be achieved in each stage) of TSOARPD, the simulation time of the complete TSOARPD will not differ much from the time taken by MOOPF.

2.3. Objective Functions of TSOARPD

The objectives of the considered TSOARPD are as follows:Objective-1: Operating cost minimization in stage 1where PGi indicates the active power generation at bus i, the cost coefficients are denoted by ai, bi, and ci for the ith generator, and NG denotes the number of generators.Objective-2: Maximization of VS of the system by maximizing the system’s effective reactive power reserve (ERPR) in stage 2The mathematical steps for evaluating ERPR are explained in the section as follows.

2.4. Evaluation of Deterministic ERPR

The RPR refers to the generators’ capacity to maintain bus voltages in the event of increased load or disturbance. The armature and field heating, as well as the current operating point and placement of the generator in the network, determine the quantity of reactive power that can be delivered to the system [32]. Because of the field current limit, the upper bound on the reactive power output of an alternator can be written aswhere IGi, f.max is the maximum field current of the ith generator. Xsi, VGi, PGi, and PRi are the synchronous reactance, terminal voltage, present, and rated active power outputs of the ith generator, respectively.

The maximum reactive power generation due to the armature current limit is given as follows:where IGi, a.max represents the maximum armature current for the ith generator. Now, for ith generator, the available reactive power reserve (ARPR) can be given aswhere QGi.max in equation (6) is the smaller among the values obtained in equations (4) and (5).

Due to network constraints, the ARPR of the alternators cannot be fully utilized by the system. As a result, equation (7) gives the effective reactive power reserve that the network can use.where pfGi signifies the participation or weight factor of generator-i, which will be computed based on the generator’s relative participation (reactive power injection) to the increased reactive demand [31]. Thus, if ΔQGi is the change in the reactive power generated by the ith generator due to an increment in the reactive power load ΔQD in the system, then pfGi is computed as

2.5. Probabilistic OPF

Because the power systems contain numerous unpredictable aspects, using a deterministic OPF that ignores these factors may sometimes result in infeasible solutions. Let Ω be a multidimensional vector of random variables that stands for the uncertain parameters of the power system. Then, the nonlinear probabilistic optimization problem by incorporating Ω into the deterministic OPF [2] can be written as

As the input Ω is random, the output will also become random. Thus, the goal of POPF is to obtain the optimal output variable’s distribution while keeping all dependent variables within their bounds. POPF will produce stochastic information in the form of mean, standard deviation, and higher moments for the optimal output and control variables.

2.6. Methods to Solve POPF

The following two probabilistic approaches are utilized in this paper to solve the POPF problem.

2.6.1. MCS for Solving POPF

The MCS method generates output random variables from a large number of sample vectors of input random variables, using each input sample vector to solve the deterministic OPF [2]. Let the ith output random variable to be considered is yi. The mean (µ), standard deviation (σ), skewness (Sk), and kurtosis (Ku) of yi will then be determined using the MCS technique as follows:where NMCS denotes the number of input sample vectors taken. Usually, NMCS will be a large number (10,000–20,000).

2.6.2. PEM for Solving POPF

By employing the data of input random variables, the PEM uses a few deterministic OPF evaluations to generate the output variables’ stochastic information. The 2m+ 1 scheme of PEM will need 2m + 1 evaluation of deterministic OPF, where the number of input random variables (P, Q loads and PWF, QWF, etc.) is denoted by m. In this scheme, two estimated points and a third point formed with the mean values of input variables will be used to solve POPF. The computational steps for solving the POPF problem by using 2m+ 1 PEM [1] are given as follows.

Calculate the standard locations ξl,k, weights l,k, and estimated points (or locations) pl,k of the lth input random variable pl by using the following equations:where skewness and kurtosis of pl are denoted by λl,3 and λl,4, respectively.

The above setting of ξl,3= 0 gives pl,k=µpl, and this indicates that m of the 3m locations will represent the same vector (µp1, µp2,…, µpl,…, µpm). Thus, at this location, only one OPF evaluation is required. Here, 0 denotes the weight for this location.where µpl and σpl denote the mean and standard deviation of pl, respectively.

After obtaining the locations and weights, the deterministic OPF will be solved for each vector (µp1, µp2, … , pl,k, … , µpm). The output of the OPF problem iswhere Y(l,k) denotes the vector of output variables for the kth concentration of pl.

From Y(l,k), the output variables’ raw moments are computed using the following equation:

From the above calculated moments, the PDF and CDF can be evaluated using the appropriate expansion series; the Gram–Charlier series [4] is used in this paper. The details about the Gram–Charlier type-A series are given in Appendix A.

2.7. Importance of considering Correlation between the Input Random Variables

In the real-time data of input random variables of the power system, there always exists some correlation between the input random variables. Ignoring these correlations leads to a reduction in the accuracy of POPF results. The above-discussed conventional 2m+ 1 scheme of PEM considers the input variables independent of each other. Considering correlation among the input random variables into the conventional 2m+ 1 scheme of PEM [16] involves the following steps:Step 1: Convert all the correlated input random variables (Xi) following different marginal distributions into correlated standard normal variables (Zi) by using the Nataf transformation as given as follows:where F and Ф represent the CDF of marginal and standard normal distributions, respectively. Now, calculate the transformed correlation matrix ρZ of Z from the original correlation matrix ρX of X by solving the following relationship between the elements of ρZ and ρX.where ρx(i,j) indicates the element at ith row and jth column of ρX.where ρ = ρz(i,j). µxi, µxj, σxi, and σxj are the mean and standard deviation values of ith and jth variables, respectively.Step 2: Apply Cholesky decomposition to factorize ρZ.where L in the above equation is a lower triangular matrix.Step 3: Convert the correlated standard normal random vector Z into an uncorrelated random vector G of the standard normal type, byCalculate the standard locations (ξ), weights (), and estimated points ,k (l = 1, …, m, k = 1, 2) for the uncorrelated standard normal vector G by using equations (11)–(14).Step 4: Convert the above calculated uncorrelated estimated points of input variables s,k = [1,k, 2,k, …, m,k]T (where k = 1, 2) back into correlated estimated points byThe above transformation gives a column vector containing correlated estimated points of input variables in standard normal form.Step 5: Apply reverse Nataf transformation to convert us,k into the marginal distribution variables.where the column vector xs,k represents the estimated points of input random variables following marginal distributions. Now, use the elements of xs,k instead of pl,k in equation (16) for solving OPF.

When employing the Nataf transformation, the computation of the coefficients ρz(i,j) may be challenging because it requires solving the integral equation (20), which is not always certain to have a solution, especially if ρx(i,j) is too close to 1 or −1. So in this paper, the Gauss–Hermite quadrature is used to solve the integral equation given by equation (20). The mathematical procedure of this method is provided in the section as follows.

2.8. Gauss–Hermite Quadrature

The GHQ is a type of Gaussian quadrature [36, 37] that is used in the numerical analysis to approximate the integrals’ values of the following form:

By using GHQ, the above integral can be approximated aswhere n represents the considered number of quadrature nodes. And the kth root (zero) of the Hermite polynomial (represented by equation (27)) is denoted by tk

The weights ωk for the corresponding roots are

When GHQ is applied for solving equation (19), then

A detailed derivation of equation (29) can be found in [36]. Now, equation (29) can be solved to obtain the value of ρz from the known ρx value by using any numerical method or by fitting a polynomial curve for equation (29). This paper uses polynomial curve fitting to solve this equation, which involves the following steps:(i)Assume uniformly distributed values for ρz in the range [−1, 1].(ii)Solve equation (29) for the assumed values of ρz to obtain the corresponding ρx values.(iii)Using all the [ρz, ρx] pairs, establish a polynomial in ρz using curve fitting. This polynomial will be of the following form:(iv)As ρz = 0 for ρx = 0, set the constant coefficient to zero (c0 = 0) in the above equation.(v)Now, for any known value of ρx, the corresponding ρz can be computed from the above-obtained polynomial.(vi)A valid ρz should satisfy the conditions |ρz| ≤ 1 and ρxρz > 0.

It is worth mentioning that when both the correlated variables follow the normal distribution, then ρx and ρz will be equal. Here, in this paper, as only correlations of loads (following Normal distribution) are considered, ρz can be directly taken as ρz = ρx. But to present the complete procedure of calculating ρz, the above steps are discussed, which will work efficiently even when the correlations between other types of random variables (following other distributions) are considered. In this paper, instead of taking ρz = ρx directly, the ρz is calculated using the above steps by taking the 7th-order polynomial. Here, this procedure is able to compute the ρz value with a maximum error of the order of 10−8, proving its efficiency.

2.9. Uncertainty Modeling

This paper considers load and wind uncertainties. The uncertainty related to load variation is modeled using normal distribution [17, 38]. The PDF of a normal random variable X (having mean µ and standard deviation σ) is shown as follows:

The wind speed uncertainty is modeled using Weibull distribution [7, 39]. The PDF of a Weibull random variable (having k > 0, c > 0 as shape and scale factors) is

Now, the wind generator’s active power output () is computed using the quadratic power curve [5, 40] as given as follows:where Pwr and denote the rated and the actual active power outputs of the wind generator, and represent the actual and the rated wind speeds in m/s, and cin and cout denote the cut-in and cut-out speeds of the wind turbine (in m/s), respectively.

2.10. Solution Tool for Solving TSOARPD

While solving POPF problems, using metaheuristic algorithms for solving the deterministic OPF will increase the simulation time of complete POPF. So the simulation time of POPF can be reduced by using conventional optimization methods for solving deterministic OPF. In this paper, the interior-point algorithm [41] is used for solving the deterministic TSOARPD. As the interior-point method has few advantages, such as ease of handling inequality constraints and speed of convergence, and it does not require a strictly feasible initial point, this method is used for solving large-scale problems such as OPF in the literature [4244]. Interior-point can handle large, sparse problems as well as small dense problems. In this method, as the barrier function (internally calculated) iterates away from the inequality constraint boundaries, the results may be slightly less accurate. This inaccuracy is relatively small for most practical purposes. However, this inaccuracy can further be reduced by taking smaller constraints and optimality tolerances.

The flowchart for solving POARPD using 2m+ 1 PEM considering correlation is shown in Figure 1.

2.11. Voltage Stability Index (L-Index)

The maximum L-index value among the load buses [27, 28] is used in this paper as a proximity indicator for the static voltage stability. Based on the power-flow solution, the L-index gives the distance between the present state and the stability limit of a power system. The numerical computation of the L-index is fast and simple, and its calculation steps can be found in [28]. The value of the L-index varies in the VSM between 0 (no-load state) and 1 (voltage collapse state). Thus, for a power system to improve the VSM, system L-index needs to be moved far away from 1 (towards 0). For a power system with NPQ number of load buses, the system’s L-index can be written as follows:

2.12. Inclusion of Voltage-Dependent Loads in OPF

The real-time electrical loads are voltage-dependent in nature; this dependency affects the OPF calculations. The primary models of electrical loads are constant power load, constant current load, and constant impedance load models [22]. But the real-time electrical loads are not purely of one model; i.e., they are composed of all the three models mentioned above [45, 46]. Thus, a composite load model (ZIP model) can be used to represent the real-time loads in OPF. The ZIP model for active and reactive loads [22, 47] can be given as follows:where P0 and Q0 are the active and reactive load powers in the constant power model of load, and V is the bus voltage magnitude per unit. The coefficients a and b depend on the load composition, and the sum of all three components of a (or b) will always be equal to one.

2.13. Wind Generator Models

Two alternative wind generator models are used in this paper. The quadratic power curve indicates the value of the active power output for a given wind speed . Based on the model of the wind generator, the reactive power will be computed.

2.13.1. Synchronous Generator Model (SG)

In this synchronous generator model (used in WTGS of variable speed type), the power factor cos ϕ or the reactive power will be specified [48]. When cos ϕ is specified, then can be computed as

2.13.2. Pitch-Regulated Induction Generator (IG)

Solving a quadratic equation that includes the induction generator’s slip (s) is required to calculate the reactive power for this model [48], as given as follows:wherewhere V denotes terminal voltage.

Now, from the above-calculated values of a, b, and c, the slip will be computed from the following equation:

From the above-obtained s value, the reactive power will be calculated as follows:

The values used in this paper for the parameters in the above equation are given in Appendix B.

3. Test Systems

In this paper, POARPD is solved for modified IEEE-9 bus and IEEE-30 bus test systems. The original IEEE-9 bus system does not contain switchable shunt capacitors and tap-changing transformers, which makes the optimization problem further easy as only the generators’ voltage magnitudes and active powers will be the control variables. So, to increase the complexity of the optimization problem (to test the effectiveness of the methods used in this paper while solving complex optimization problems), a modified version of the IEEE-9 bus system is used in this paper. To modify the IEEE-9 bus system, a switchable shunt capacitor and a tap-changing transformer are added to the original system. The network, load, and generator data have been taken from MATPOWER [49].

3.1. Modified IEEE-9 Bus Test System

A switchable shunt capacitor of 0.035 p.u is added to the original IEEE-9 bus test system at bus 5, and a tap-changing transformer (with a ratio of 0.95) is considered in the line between the buses 8 and 9 to modify the test system. Along with three generators, the test system has base-case active and reactive loads of 3.15 p.u and 1.15 p.u, respectively. The slack bus is bus 1, and the system has nine branches. All the bus voltages have upper and lower limits of 1.1 and 0.9 p.u. The upper and lower constraints on active power generation by generators are [2.5, 3, 2.7] p.u and [0.1, 0.1, 0.1] p.u, respectively. The tap ratios’ lower and upper bounds are 0.9 and 1.1, respectively, while the switchable shunt capacitors’ minimum and maximum limits are 70% and 130% of the base case values, respectively. For the generators, the maximum field current and the synchronous reactance are assumed as [2.1, 2.15, 2.05] p.u and [1.05, 1.15, 1.0] p.u, respectively.

3.2. IEEE-30 Bus Test System

The standard IEEE-30 bus test system contains 6 generators, 41 branches, 2 switchable shunt capacitors, and 4 tap-changing transformers. The active and reactive loads in the base scenario are 2.834 p.u and 1.262 p.u, respectively. All bus voltages have upper and lower limitations of 1.06 p.u and 0.94 p.u, respectively. The generators’ active power generation upper and lower bounds are [3.602, 1.4, 1, 1, 1, 1] p.u and [0, 0, 0, 0, 0, 0] p.u, respectively. The tap ratios’ lower and upper bounds are 0.9 and 1.1, respectively, while the switchable shunt capacitors’ minimum and maximum limits are 70% and 130% of the base case values, respectively. For the generators, the maximum field current and the synchronous reactance are assumed as [2.2, 2.1, 2, 2, 2.15, 2] p.u and [0.9, 1.05, 0.9, 0.9, 1, 1.1] p.u, respectively.

For all the test cases, the load uncertainty is modeled by taking the base-case loads as the mean values and 7% of the base-case loads as the standard deviations. Tables 1 and 2 provide the wind speed characteristics, the type of WTGS employed, the number of wind generators considered, and the location of wind farms (WFs) used in this paper for the test systems.

The rating taken for SG-type wind generators is 0.55 MW, while the rating taken for IG-type wind generators is 0.5 MW. The considered rated speeds for SG- and IG-based wind farms are 16 and 14 m/s, respectively. For the SG- and IG-based wind farms, the cut-in and cut-out speeds are taken as [4.5, 4.2] m/s and [55, 24] m/s, respectively. The SG-type wind generator’s power factor is taken as 0.85 lag. WTGS-rated capacities for the above 9-bus and 30-bus systems are thus 42 MW and 64 MW, respectively. The IG-type wind generator’s parameter details are given in Appendix B.

For all the test cases, all the simulation works using MATLAB environment have been carried out on an Intel Dual-core 3.6 GHz, 16 GB RAM machine.

4. Results and Discussions

In this paper, the POARPD is solved for the following scenarios:Case 1: The constant power load model is considered for solving the objective functions, whereas the complete POARPD is solved by using the 2m + 1 scheme of PEM neglecting correlation (UCPEM) and the 2m + 1 scheme of PEM considering correlation (CPEM). And the same POARPD problem is solved using the MCS method for comparison.Case 2: The POARPD problem is solved by considering purely voltage-dependent loads using CPEM and MCS.

For all the cases discussed, 10000 samples of input random variables are taken, and the POARPD problem is solved for both the modified IEEE-9 bus and IEEE-30 bus test systems. The correlated standard normal samples are generated using the MATLAB inbuilt function mvnrnd() and later converted into the correlated marginal distribution samples based on the CDF of the corresponding correlated input variables. These correlated marginal distribution samples will have correlation coefficients that deviate (small deviations) from the correlation coefficients supplied to the mvnrnd() function. The new correlation coefficients among these correlated marginal distribution samples are considered (to use in CPEM) instead of the supplied ones to maintain higher accuracy.

For the modified IEEE-9 bus system, all the active and reactive loads are considered correlated while generating the input random sample vectors for all the cases. And for the IEEE-30 bus system, the active and reactive loads at the buses 19, 21, 30 are considered correlated for generating input random sample vectors for all the cases. For both the test systems, the input sample vectors are generated such that the wind speed samples are independent of each other; furthermore, there is no correlation between wind speeds and loads for all the test cases (i.e., only the loads are correlated). The test systems’ correlation coefficient matrices are given in Appendix B. The PDF of the wind speed at bus 7 and the PDF of the total active power of the considered IEEE-9 bus system are shown in Figures 2(a) and 2(b).

The deterministic results for the base case and other types of OPF are compared in Table 3. The results for the base-case settings of control variables are denoted by “Base.” The results of multiobjective OPF for OC minimization and ERPR maximization considering all the active and reactive control variables are denoted by MOOPF. The results of the two-stage OARPD used in this paper are denoted by TSOARPD. For ERPR and L-index, the better performance of TSOARPD compared to MOOPF can be observed from the results in Table 3.

In Table 3, the L-index is the maximum among all the system load bus L-index values. From the deterministic results (Table 3) for the 100% and 120% loading conditions, it can be observed that TSOARPD gives an optimal OC, which is a little bit higher than the OC provided by MOOPF, but at the same time, the optimized ERPR is higher and gives a better L-index value (smaller than the L-index obtained by MOOPF). But this increment in OC will be a very small percentage. Thus, the TSOARPD gives better VS by keeping the OC close to its minimum value when the system is substantially loaded (or in contingencies when it is required to provide more importance to VS than OC). As improving the VS leads to minimizing the active power losses of the system, it will again help maintain the OC close to its optimal value. Thus, optimizing OC and VS in two stages will achieve the desired active and reactive objectives one after the other. Table 3 shows that the MOOPF and TSOARPD are optimizing the reactive power generations (QG) of the conventional generators using the reactive control variables (as the active control variables are fixed in stage 2). This process of optimizing QG is based on the generator’s participation factor (pfG). The QG of the generator having high pfG will be moved away (below the maximum value of QG) so that the ARPR of that generator is increased, maximizing ERPR. Now, as the generator with high pfG has high ERPR, it will respond (more) to the load increment or contingency in the system and provide more reactive power to the system. Thus, the reactive power that the system can utilize (at that particular loading state) increases, increasing the system’s capability to accommodate more load and improving the system’s VS. The convergence characteristics of deterministic TSOARPD are given in Figure 3. Figure 3(a) presents the optimization of OC (stage 1). In stage 2 of TSOARPD; the ERPR is maximized by minimizing the negative value of ERPR. Thus, Figure 3(b) shows the minimization of “-ERPR”. The time taken by MOOPF and TSOARPD for the 100% loading case is 1.21 and 1.34 seconds, respectively.

The load bus voltages of the considered IEEE-30 bus system are shown in Figure 4. For constant power load case (Figure 4(a)), TSOARPD is able to maintain higher load bus voltages which will increase the amount of load that can be further served by the system. For voltage-dependent load (Figure 4(b)) (completely voltage-dependent) case, the load bus voltages provided by MOOPF and TSOARPD are almost the same.

The POARPD results for the IEEE-9 bus system for case 1 are given in Table 4, in which VLAvg represents the average load bus voltage. For the POARPD results obtained by PEM methods, the errors that occurred in the calculation of mean and standard deviations are calculated with respect to MCS as shown in equation (41). In equation (41), U denotes the parameter (µ or σ, etc.) for which the error is to be calculated. For the IEEE-9 bus system, the corresponding errors in the calculation of POARPD results by UCPEM and CPEM are given in Figure 5.

From the numerical results in Table 4, it can be observed that the mean values of the optimal output variables are not the same as the optimal values of output variables obtained from the deterministic TSOARPD (given in Table 3). This difference is due to the nonlinear nature of the power-flow equations (the mean values of output variables will not precisely match with the output variables’ values obtained at the mean values of input variables). Table 4 shows that with the reduction in L-index value, the active and reactive power losses are minimized. The average active power losses for the base case are 0.049 p.u which are reduced to 0.0370 p.u after POARPD, and the average reactive power losses for the base case are 0.474 p.u which are reduced to 0.3732 p.u, after POARPD. The average load bus voltage is increased from 1.0196 p.u (base case) to 1.0876 p.u (after POARPD), and the weak bus voltage (Vbus9) is improved from 1.0196 p.u (base case) to 1.0615 p.u (after POARPD). As the load bus voltages are improved, now the system can accommodate more loads without losing stability.

From the graphs in Figure 5, it is clear that considering correlations among the input random variables while calculating locations (or estimated points) in PEM will increase the accuracy of the results. When the correlation is neglected in UCPEM, even though the errors in the mean values of output variables are small, the corresponding errors in the standard deviation are much higher. When the correlation is considered (in CPEM), the errors related to the mean and standard deviation calculation are significantly reduced. In the above accuracy comparison (Figure 5), the maximum errors by UCPEM are 0.54% for the mean and 33.3% for the standard deviation, and those by CPEM are 0.27% for the mean and 16.7% for the standard deviation. In the case of control variables, the highest errors by UCPEM and CPEM for mean computations are 0.21% and 0.066%, respectively. The highest errors for standard deviation computations using UCPEM and CPEM are 61.5% and 37.5%, respectively. As a result, the CPEM can provide a more accurate standard deviation and other higher moments of output variables, such as skewness and kurtosis, which results in a more precise probability distribution of the output variables. Figure 5(a) presents the errors in the mean value calculation. Figure 5(b) presents the errors in the standard calculation for the output variables of IEEE-9 bus system for case 1.

The most relevant outcomes of the POPF are the mean and the standard deviation values of output variables; nevertheless, the accuracy of the estimated higher moments of the results can be compared using CDF or PDF. This comparison is given in Figure 6.

The CDF graphs in Figure 6 show that not only the values of mean and standard deviation but also the higher moments (of the output variables) estimated by CPEM are similar to those obtained by MCS. Figure 6(a) shows the CDF of OC, Figure 6(b) shows the CDF of ERPR, and Figure 6(c) shows the CDF of maximum L-index for the modified IEEE-9 bus system for case 1.

For case 2, all the active and reactive loads are considered purely voltage-dependent, whereas each load is composed of a 50% constant current and 50% constant impedance nature. This ratio of the composite loads is applied to consider the increased amount of nonlinearities in the system. The numerical results of POARPD for the considered IEEE-9 bus system for case 2 are given in Table 5.

For the considered IEEE-9 bus system, probability distributions of the reactive demand on the system for case 1 and case 2 are given in Figure 7.

When the system loads are considered purely voltage-dependent, the TSOARPD minimizes the load bus voltages, thereby minimizing the load magnitudes in order to optimize the objective values. The OC has decreased due to the reduced active load on the system in case 2. In this scenario, even though the reactive load on the system is reduced (Figure 7), there is a reduction in the ERPR. This reduction in ERPR is due to the decrease in bus voltages. Reduced bus voltages lower branch charging, which causes generators to generate more reactive power. As a result, the ERPR of the system decreases as the generators’ reactive power generation increases.

Table 6 shows the numerical results of POARPD for the IEEE-30 bus system for case 1. As the generated samples of independent variables (apart from the variables considered to be correlated) also contain correlations of small magnitudes, ignoring these small magnitudes of correlations still affects the accuracy. Here, only the correlations between the loads at three buses are considered to demonstrate the significant improvement in the accuracy with the consideration of even a small portion of the correlated behavior of the input variables. As only the correlation between six input variables (P and Q loads) is considered for the IEEE-30 bus system, the results generated by UCPEM and CPEM are nearly identical. If the correlation between all the input variables is assessed in CPEM, the difference will be more significant (and more accurate). The following numerical results show that the accuracy has improved even with a minor consideration of correlation. Even for the IEEE-30 bus system, the accuracy of the 2m + 1 PEM scheme is good because it is independent of the number of input random variables. Table 6 shows that the average active power losses are reduced from 0.176 p.u (base case) to 0.145 p.u after POARPD, and the reactive power losses are reduced from 0.677 p.u (base case) to 0.535 p.u after POARPD. The bar graphs in Figure 8 show the errors in the computation of POARPD output variables of the IEEE-30 bus system for case 1. Figure 8(a) shows the errors in the mean value calculation, and Figure 8(b) shows the errors in the standard calculation for the output variables of IEEE-30 bus system for case 1.

In the comparison in Figure 8, the maximum errors by UCPEM and CPEM for the output variables of POARPD for the IEEE-30 bus system are 1.69%, 1.69% (for mean values), and 62.4%, 46.5% (for standard deviation values), respectively. The further errors are since the correlations between most of the input variables of the IEEE-30 bus system are ignored in this case. The consistency in the accuracy of the 2m + 1 scheme (irrespective of the number of input random variables) can be observed from the above results. Furthermore, the precision with which the higher moments are calculated can be seen by comparing the probability distributions. This comparison is shown in Figure 9.

The numerical results of POARPD for the IEEE-30 bus system for case 2 are given in Table 7. The probability distribution of ERPR using CPEM for the IEEE-30 bus system for case 2 is given in Figure 10. The variation of Ploss for case 1 and case 2 is also given in Figure 10. Figure 10(a) presents the CDF of ERPR for case 2, and Figure 10(b) presents the CDF of Ploss for cases 1, 2 for IEEE-30 bus system. Due to the reduced load (active and reactive) in case 2, there will be a reduced power transfer in the system. Due to this reduction in the line flows, the magnitudes of Ploss and Qloss are decreased in case 2. But due to the reduced bus voltage magnitudes in case 2, the Ploss and Qloss will not be too far below those values obtained in case 1.

From the numerical results of case 2, the effect of the TSOARPD’s increased nonlinearity on the accuracy of CPEM can be observed. As the active and reactive loads are purely voltage-dependent in case 2, the net load values on the system will vary nonlinearly (with an order of 2) in the deterministic TSOARPD. Due to this nonlinear variation of loads, the approximated values (based on a few evaluations) of output and control variables of POARPD will not exactly be the same as those obtained by evaluating POARPD for all the possible vectors of loads (using MCS). Due to this increased nonlinearity in TSOARPD, the accuracy of CPEM will decrease, but the accuracy will still be in an acceptable range.

In all of the above CPEM accuracy comparisons, the maximum errors for the IEEE-30 bus system are slightly greater than those for the considered IEEE-9 bus system (for both mean and standard deviations) in all the cases. And in case 2, due to the increased amount of nonlinearity in the deterministic TSOARPD, the maximum errors for case 2 are slightly higher for both the test systems. However, the errors related to the standard deviation are significantly reduced by using CPEM. A similar improvement exists in the accuracy of the calculation of higher moments also, leading to obtaining accurate data of CDF and PDF of the output variables.

From the numerical results provided by both the deterministic TSOARPD and the POARPD, it can be observed that the absolute error in the optimal mean values of all the output variables obtained by POARPD concerning the deterministic optimal values obtained by TSOARPD is of very small magnitude. Thus, the optimal values given by the deterministic TSOARPD are acceptable as the mean values of system variables. But the primary purpose of POARPD is to provide information about the uncertainty range of output variables, which cannot be obtained by deterministic OARPD. The significant reduction in the simulation time provided by CPEM compared to MCS can be observed in Table 8. The computational time for CPEM is a little higher than UCPEM due to its algorithm’s extra steps, but this increased computational time is negligible compared to MCS. By using CPEM, calculating the stochastic information of the output variables within the practically acceptable time (with a little bit of sacrifice in the accuracy) will enable us to make fast decisions while dealing with a real-time power system, whereas using MCS will sometimes take hours to provide the results.

5. Conclusions

A POARPD method including the load and wind uncertainties considering the correlation between loads is presented in this paper. In this POARPD, the deterministic OARPD is performed in two stages to achieve better optimization of ERPR, thereby better voltage stability compared to MOOPF. This two-stage OARPD can provide better optimal voltage stability (in terms of both ERPR and L-index) with a bit of sacrifice in the optimal operating cost during normal and heavily loaded conditions of the system. To demonstrate this, TSOARPD is solved for 100% and 120% loading conditions and compared with the corresponding results obtained by MOOPF.

To demonstrate the accuracy of CPEM, the POARPD problem is solved for the modified IEEE-9 bus and for the IEEE-30 bus systems, and the performance is compared against MCS. The results indicate that the correlations among the input variables (such as loads and wind speeds) significantly impact the outputs of POARPD. By considering the correlations among the input variables, the CPEM can give more accurate uncertainty ranges of output variables than the traditional uncorrelated PEM (UCPEM). When a power system contains WTGS, the ERPR of the system varies over a wide range with a decreased magnitude (and the L-index with an increased magnitude) compared to the system without WTGS. This wide range of variation is due to the increased amount of system uncertainty (due to the additional uncertainty of WTGS along with load uncertainty), and the change in the magnitude is due to the consumption of reactive power by the WTGS. When the power system loads are considered voltage-dependent (even if the system contains WTGS or not), the system’s ERPR got reduced, and the system’s L-index got increased. Both of these changes indicate the system’s voltage stability reduction when the loads are voltage-dependent. So, considering the voltage-dependent loads gives a more realistic state of voltage stability of the system.

Acronyms

ARPR:Available reactive power reserve
CDF:Cumulative distribution function
CPEM:Correlated PEM (PEM considering correlation)
ERPR:Effective reactive power reserve
GHQ:Gauss–Hermite quadrature
IG:Pitch-regulated induction generator
MCS:Monte Carlo simulation
MOOPF:Multiobjective optimal power flow
OAPD:Optimal active power dispatch
OARPD:Optimal active and reactive power dispatch
OC:Operating cost
OPF:Optimal power flow
ORPD:Optimal reactive power dispatch
PDF:Probability density function
PEM:Point estimate method
PLF:Probabilistic load flow
POARPD:Probabilistic optimal active and reactive power dispatch
POPF:Probabilistic optimal power flow
RE:Renewable energy
RPR:Reactive power reserve
SG:Synchronous generator
TSOARPD:Two-stage optimal active and reactive power dispatch
UCPEM:Uncorrelated PEM (PEM neglecting correlation)
VS:Voltage stability
VSM:Voltage stability margin
WTGS:Wind turbine generating systems.

Data Availability

The data supporting the current study are available from the corresponding author upon request. For data-related queries, kindly contact Baseem Khan [email protected].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Supplementary Materials

In the supplementary file, Appendix A presented the Gram–Charlier type-A series description which is utilized in this work. The expansion of this series is utilized for the computation of the probability distribution function (PDF) and cumulative distribution function (CDF). In Appendix B, Table B1 is utilized for providing the parameters of the induction generator type wind generator. Furthermore, in Appendix B, Table B2 is utilized for showing the diagonal and off-diagonal elements of the load correlation matrix for the IEEE 9 bus system, utilized in this work. Moreover, in Appendix B, Table B3 is utilized for showing the diagonal and off-diagonal elements of the load correlation matrix for the IEEE 30 bus system, utilized in this work. (Supplementary Materials)