Scholarship at UWindsor Scholarship at UWindsor Multi-objective optimal reactive power planning under load Multi-objective optimal reactive power planning under load demand and wind power generation uncertainties using ε - demand and wind power generation uncertainties using constraint method constraint method

: This paper presents an improved multi-objective probabilistic Reactive Power Planning (RPP) in power systems considering uncertainties of load demand and wind power generation. The proposed method is capable of simultaneously (1) reducing the reactive power investment cost, (2) minimizing the total active power losses, (3) improving the voltage stability, and (4) enhancing the loadability factor. The generators’ voltage magnitude, the transformer’s tap settings, and the output reactive power of VAR sources are taken into account as the control variables. To solve the probabilistic multi-objective RPP problem, the ε -constraint method is used. To test the e ﬀ ectiveness of the proposed approach, the IEEE 30-bus test system is implemented in the GAMS environment under ﬁve di ﬀ erent conditions. Finally, for a better comprehension of the obtained results, a brief comparison of outcomes is presented.


Introduction
Reactive Power Planning (RPP) in power systems can be considered as one of the most difficult and complicated problems due to its complex variables, constraints, and optimization algorithms [1]. It is related to optimal sizing and allocation of VAR sources in power systems to satisfy prescheduled objectives, such as determining the optimal allocation and minimizing the operation costs [2,3]. The main objective of RPP is to achieve feasible operation with a satisfactory voltage profile with a lack of VAR support conditions. According to the concept of VAR planning in power systems, various objectives functions can be defined for the RPP problem. These objectives may consist of cost-based objective functions or objective functions that maximize or minimize indices, such as voltage stability margin or system loadability [4,5]. Moreover, it is possible to express the RPP as a multi-objective optimization problem, which optimizes several goals simultaneously [1].
Moreover, there is an increasing interest in using Renewable Energy Resources (RESs), such as wind farms and solar power plants, in power systems due to their technical, environmental, and economic advantages [4][5][6]. However, with the high penetration of RESs in power systems, the challenges associated with RPP are dramatically increased. One of the main challenges that can affect the RPP is the uncertainty in the generation availability of RESs. Uncertainty in the sources' parameters leads to difficulties in proper decision-making in the planning of power systems. Furthermore, owing to the stochastic nature of the load demands in electric power systems, additional uncertainties should be considered in RPP.
In RPP research studies, the probabilistic decision-making process based on either source uncertainty or load demand uncertainty is a well-developed research topic. Nevertheless, the probabilistic multi-objective RPP in power systems considering the uncertainty of loads and wind farms at the same time has not been fully investigated. In [7], a novel approach for dynamic VAR planning to improve the short-term voltage stability and transient stability is proposed. The impact of FACTS devices in RPP is analyzed in [8,9]. However, in both studies, an attempt has been made to explain the problem in a deterministic context. A multi-objective RPP that mainly focuses on voltage stability is introduced in [3]. Nonetheless, it is modeled based on a deterministic approach. In [10], a multi-objective approach for RPP with wind generations is presented. In this study, various objectives, such as system loadability, power losses, and cost of reactive power investment are considered. In [11], the RPP is solved using the Genetic Algorithm (GA) to reach coordination in controlling the reactive power in the presence of wind farms and FACTS devices. The loadability factor of the system is optimized by the optimal allocation of wind farms and FACTS devices. This procedure is implemented when loads with constant Power Factor (PF) and wind farms without uncertainty are assumed. Using the Benders decomposition method and considering the high penetration of wind generation, the RPP problem is tackled as two-stage stochastic programming in [12]. Using the Differential Evolutionary Algorithm (DEA), the RPP is solved in a wind integrated system in [13]. A major problem with the suggested model is that it only includes the uncertainty in wind power generation. In [14], a multi-stage stochastic model for RPP is extended, which involves the uncertainty of loads. Nonetheless, the proposed model describes the probabilistic behavior of the system in the absence of wind farms. In [15], a mixed-integer quadratic model for long term VAR planning is proposed. An attempt was made to minimize the operation and investment cost of new VAR sources and the load shedding risk through multi-objective optimization. Though the uncertainty in demand is completely taken into account, the proposed model does not consider uncertainty in wind power generation. In [16], a stochastic model based on chanced constrained programming for RPP is defined. The proposed model is solved using GA. Although the uncertainty is modeled in the power generation, it optimizes only one objective, including operational and investment costs. A chanced constrained model is proposed for probabilistic RPP in [17]. The proposed model is solved through two-stage stochastic programming. The main disadvantage of the proposed model is that it only considers the load as a random parameter. Besides, only the investment cost of new VAR sources is taken into account as the main objective function.
The main drawback of all the mentioned research studies is that the optimal RPP considering load demand and wind power generation uncertainties at the same time are not fully investigated. This paper aims to address RPP as a probabilistic multi-objective problem in order to reduce the total cost of reactive power investment, minimize the active power losses, maximize the voltage stability index, and improve the loadability factor. The generators' voltage magnitude, the transformers tap settings, and the output reactive power of the VAR sources are considered as the main control variables. To cope with the probabilistic multi-objective RPP problem, the ε-constraint technique is employed. To validate the efficiency of the proposed method, the IEEE 30-bus test system is implemented in the GAMS environment under five various conditions. The obtained results show the effectiveness and high accuracy of the proposed method. Table 1 shows a comparison between the proposed probabilistic multi-objective RPP and the previously published research papers.
The rest of this paper is organized as follows. Section 2 deals with the uncertainty modeling. Problem formulation is presented in Section 3. Section 4 describes the optimization method. Simulation results are given in Section 5. Finally, some brief conclusions are summarized in Section 6.

Uncertainty Modeling
In this section, the uncertainties in load demand and wind power generation in the RPP problem are modeled to cope with the stochastic nature of the load demand and wind power generation. In the following subsections, modeling of both the load demand and wind power uncertainties are described. Finally, modeling of the system uncertainty via scenario generation is presented.

Modeling the Load Demand Uncertainty
The uncertainty of the load is usually modeled by the normal distribution with mean (µ) and standard deviation (σ) [18]. In this paper, it is assumed that all the loads have constant PF, the same mean, and standard deviation. Therefore, for simplicity, a normal distribution is applied at the load level (λ) instead of applying in each load independently. The probability of each load level is shown by (π l ), and is calculated using Equation (1). The associated value of each load level is denoted by (λ l ), and can be obtained using Equation (2) [19]. It is worth mentioning that λ Min,l and λ Max,l are known as the minimum and maximum levels of the system loading at the l th load level, respectively.

Modeling the Wind Power Generation Uncertainty
Considering the intermittent nature of the wind speed, the Weibull distribution is often considered as the probability density function that can approximate the behavior of the wind with a reasonable error. Therefore, by defining the Weibull distribution for wind speed, the probability of wind speed at different intervals (scenarios) can easily be calculated. Equation (3) is a general expression for Weibull distribution [20]. Equation (4) can be used to calculate the probability of a wind speed interval (scenario). The corresponding value of each wind speed interval can be achieved using Equation (5).
where v denotes the wind speed, and α and β are the wind speed parameters that vary depending on the region in which the wind blows. Considering v i,w and v f ,w as the initial speed and final speed of the hypothetical scenarios for wind speed, the probability (π w ) of occurrence of any wind speed scenario can simply be obtained. Thereafter, the wind speed (v w ) associated with each scenario is gained using the calculated probabilities. The output power of a wind turbine is highly dependent on wind speed. Therefore, any wind turbine has a characteristic named power curve that exactly shows the capability of a wind turbine in power generation versus existing wind speed. Knowing a specific wind speed (v w ), one can estimate the output power of a wind turbine (P est w ) through its power curve. The power curve is generally defined by a set of equations as it is stated in Equation (6) [21], which in terms, (v c in ), (v rated ), and (v c out ) denote the cut-in wind speed, rated wind speed, and cut-out wind speed for a wind turbine, respectively. The rated power (P r w ) and the estimated output power of the wind turbine are also evident from Equation (6).
In most research studies, the concept of the power curve is extended to a wind farm. Hence, instead of studying a single wind turbine, it is preferable to focus on a group of wind turbines that are in a special area and usually known as wind farms.
Considering several scenarios for a probabilistic problem is generally not an easy procedure. Depending on the problem type, various methods exist for scenario generation [22,23]. However, in this paper, a technique based on [19,24,25] is applied to generate a desirable number of scenarios with reasonable accuracy. In order to have a combination of load and wind scenarios, the following steps are taken:

1.
Several scenarios for the load level are considered.

2.
The probability of each system loading scenario (level of the load) and its corresponding value using Equations (1) and (2) are calculated.

3.
Several scenarios for wind speed are considered. 4.
The probability of each wind speed scenario and its corresponding value using Equations (4) and (5) are calculated.

5.
The output power of the wind farm using the estimated wind speed in each scenario and Equation (6) is generated. 6.
The final number of combined load-wind scenarios is obtained by multiplying the number of load scenarios by the number of wind scenarios. By multiplying the probability of the load scenario by the probability of wind speed scenario, the probability of the combined load-wind scenarios (π s ) can be calculated as follows [19]: π s = π l × π w (7)

Problem Formulation
As mentioned earlier, a wide range of objective functions for the RPP in power systems can be represented. This matter enormously affects the control variables, state variables, and all constraints of the RPP problem. Thus, by a proper formulation, all objectives can be achieved, all constraints can be satisfied, and the feasibility of the problem can be ensured. Due to the fact that the probabilistic nature of the problem has a major impact on its formulation, it is very important to use the probabilistic variables accurately in the problem formulation.

Variables
The same as the other optimization problems in power systems, such as Optimal Power Flow (OPF), two types of variables, named control variables and state variables, are defined for the RPP.
Normally, for a typical RPP, control variables are defined as generators' voltage magnitudes, transformers tap settings, and output reactive power of VAR sources. Considering a scenario-based approach to model the uncertainty of the problem, the control variables set (U) for a probabilistic RPP are expressed as Equation (8) [14][15][16][17]. where V g i,s shows the voltage magnitude of the i th generator for the s th scenario, t k i,s is used to assign the settings of the i th tap-changing transformer for the s th scenario, and Q C i,s shows the output reactive power of the i th VAR compensator device for the s th scenario. Likewise, Ω g , Ω s , Ω TapCh , and Ω Comp symbolize the set of generators, set of scenarios, set of tap-changing transformers, and set of VAR compensator devices, respectively.
The state variables in a typical RPP consist of the generated active power by the slack bus, the generated reactive power by each of the existing generators, the voltage magnitude of the load buses, and the flow of the transmission lines.
Using a scenario-based approach to model the uncertainty of the problem, the state variables set (X) for a probabilistic RPP are expressed as Equation (9) [14][15][16][17].
where P G Slack,s indicates the generated active power by the slack generator (bus) for the s th scenario, Q G i,s is used to denote the generated reactive power by the i th generator for the s th scenario, V L i,s shows the voltage magnitude of the i th load bus for the s th scenario, and S From l,s and S To l,s show the apparent power flow of the sending and receiving ends of the l th line for the s th scenario, respectively. Additionally, Ω PQ and Ω Lines specify the set of the load buses and the set of transmission lines, respectively.

Objective Functions
For probabilistic multi-objective RPP, the aim is to satisfy three main objectives. These objectives include the minimization of total VAR investment cost, minimization of voltage stability index (L-index), and maximization of loadability factor, which lead to a reduction in total active power losses and improvement of voltage stability.

Minimization of Total VAR Investment Cost
One of the important objectives in the RPP is the total cost of VAR planning. In spite of allocating the optimal location and capacity for VAR sources, optimal VAR planning can handle the RPP problem from economic aspects. For this reason, the first objective function is a cost-based objective function comprising two main parts, as follows: (1). The first part evaluates the expected cost of energy loss (W c ) during the generated scenarios and is expressed as follows [16][17][18][19][20][21][22][23][24][25][26]: where P loss,s shows the active power losses during the s th scenario, t s represents the duration of the s th scenario, h is a constant parameter that is related to the first part cost-based objective function and identifies the per-unit energy cost, and π s denotes the probability of the s th scenario.
To calculate the total active power losses, Equation (11) can be used as follows [27][28][29][30]: where V i,s and V j,s are the sending and receiving ends voltage magnitude of the l th transmission line for the s th scenario, respectively, θ i,s and θ j,s are the sending and receiving ends voltage angles of the l th transmission line for the s th scenario, respectively, and G (l,s) is used to designate the conductance of the l th transmission line for the s th scenario. (2). The second part measures the expected cost of VAR investment (I c ) during the generated scenarios and is derived as follows [14][15][16][17]: where e i and C Ci are the fixed and variable installation costs of VAR sources, respectively.
Accordingly, the first objective function ( f 1 ) can be derived as follows: where F c shows the expected Total VAR Cost (TVC).

Minimization of Voltage Stability Index
In this paper, the L-index is proposed as the voltage stability index that is a well-known static voltage index [31]. In order to estimate the static voltage stability of the power system, the L-index should be calculated for all load buses (PQ buses). All the load buses that have higher values of the L-index than others are considered as the weak buses. Weak buses mostly suffer from a lack of reactive power and are prone to the voltage collapse. Equation (14) can be used to calculate the L-index (L j ) for the j th load bus, as follows [31]: where V i shows the voltage of the i th generator, and V j represents the voltage of the j th load bus. F ji can be derived from Y bus matrix of the system. Thus, by rearranging the current and voltage equations in power systems, as shown in Equation (15), the consecutive Y bus matrix is achieved. Thereafter, using the arrays of the consecutive Y bus matrix, the F ji matrix can be calculated as Equation (16).
where I g and I l show the current of generators and loads, respectively, and V g and V l are the voltage of generators and loads, respectively. In addition, Y gg , Y gl , Y lg , and Y ll are the submatrices of the consecutive Y bus matrix. It should be noted that only Y ll arrays of the consecutive Y bus matrix are related to the PQ nodes. Also, the consecutive Y bus matrix is a symmetric matrix. Therefore, Y lggl = Y lg . By minimizing the values of the L-index at the weak buses, there is a possibility to increase the level of static voltage stability in power systems. The voltage stability of power systems can be determined by the L-index when the maximum value of the L-index (L max ) is assigned to the static voltage stability level in power systems, as follows: To improve the static voltage stability of power systems, it is necessary to minimize L max . It should be noted that the equations proposed for the L-index are related to the deterministic problem. In the case of a probabilistic problem, considering all necessary modifications on the Y bus matrix in each scenario, after re-formulating Equations (14)- (17), new equations can be rewritten as follows: where L j,s indicates L-index value for the j th load bus and s th scenario, V i,s and V j,s are the voltage of the i th generator and j th load bus for the s th scenario, respectively. For each scenario, the F ji,s matrix can be calculated as Equation (20).
where I g,s and I l,s denote the current of generators and loads for the s th scenario, respectively, V g,s and V l,s are the voltage of generators and loads for the s th scenario, respectively. Also, Y gg,s , Y gl,s , Y lg,s , and Y ll,s are the submatrices of the consecutive Y bus matrix for the s th scenario.
According to the aforementioned descriptions, Equations (18)-(20) can be obtained for each scenario. The maximum value of the L-index for each scenario can be derived as follows: Consequently, the second objective function ( f 2 ), which is the expected value of the static voltage stability index during the generated scenarios, can be derived as follows:

Maximization of the Loadability Factor
The injected active power (P i ) and reactive power (Q i ) at the i th bus can be expressed in terms of the voltage (V), the elements of the Y bus matrix of the system, and the loadability factor (Γ) as follows [32]: where P Gi and Q Gi are the active and reactive power generation at the i th bus, respectively, P Di and Q Di represent the base-case active and reactive power consumption at the i th bus, respectively, and N B denotes the total number of buses.
Maximizing the loadability factor is defined as the third objective function, in this paper. However, considering the random nature of the problem, a probabilistic formulation is required. Therefore, by re-formulating Equations (23) and (24), a stochastic formula is derived to obtain the expected loadability factor, as shown in Equations (25) and (26).
where P i,s and Q i,s denote the injected active and reactive power at the i th bus for the s th scenario, respectively, P Gi,s and Q Gi,s represent the active and reactive power generation at the i th bus for the s th scenario, respectively, P Di,s and Q Di,s are the base-case active and reactive power consumption at the i th bus for the s th scenario, respectively, Γ(s) denotes the loadability factor for the s th scenario; V i,s and V j,s indicate the voltage of the i th bus j th bus for the s th scenario, respectively, and lastly, the elements of the Y bus matrix for the s th scenario are shown by Y i,j,s . According to the above-mentioned descriptions, the third objective function ( f 3 ), which is the expected value of the loadability factor, can be derived as follows: Finally, the optimization criteria subjected to equality and inequality constraints are as follows:

Constraints
The role of constraints in creating a feasible space for the problem and satisfying optimality conditions to find optimal solutions is undeniable. For this reason, the correct expression of constraints is one of the major priorities in the problem formulation.

Equality Constraints
The power flow equations are taken as the equality constraints for the RPP. Using the output power of wind farms and also considering the probabilistic nature of the problem, Equations (25) and (26) can be rewritten as follows: where P Wi,s and Q Wi,s show the output active and reactive power of the i th wind farm for the s th scenario, respectively. It should be noted that the output reactive power of the wind farms is neglected in this paper.

Inequality Constraints
To keep both the control and state variables within their specific limits, another set of constraints is added to the problem, named as inequality constraints. Those constraints involve Equations (31)- (38) [24,25].
• Limits on the Control Variables The upper limit (V max g i ) and lower limit (V min g i ) of a generator voltage magnitude for the s th scenario can be applied, as follows: For all tap-changing transformers in each scenario, the following constraint should be satisfied.
where t min k i and t max k i show the minimum and maximum settings of the i th tap-changing transformer, respectively.
The output reactive power of the VAR sources in each scenario is as follows: where Q min show the minimum and maximum output reactive power of the i th VAR compensator device, respectively.

• Limits on the State Variables
In terms of the generation units, for each scenario, two important constraints should be satisfied; (1) the limitation on the generated active power of the slack bus and (2) the limitation on the generated reactive power of each generation unit. Those constraints are given as follows: where P min show the maximum and minimum generated reactive power of the i th generator for the s th scenario, respectively. In order to prevent the voltage collapse or insulating problems, it is required to limit the voltage magnitude of loads for each scenario, as follows: where V min L i and V max L i are considered as the lower and upper limits of the voltage magnitude at the i th load bus for the s th scenario, respectively. To reduce the risk of overload in transmission lines, the apparent flow of the transmission lines should be lower than a specified value. Equations (37) and (38) enforce the apparent flow of transmission lines to be at the secure level, as follows: where S max l indicates the maximum apparent flow of the l th transmission line.

Other Considerations in the Problem Formulation
There are other considerations in the problem formulation, which are listed as follows: • The transformers tap settings and output reactive power of the VAR sources are treated as continuous variables. Therefore, the whole problem is stated as a probabilistic multi-objective nonlinear problem.
• Since the Y bus matrix of power systems is dependent on the transformers tap settings and due to the fact that the transformers tap settings are defined as scenario-dependent variables, the Y bus matrix should be calculated for each scenario separately.

•
The L-index value varies between 0 and 1 for power systems. It should be noted that except for the defined boundaries, the L-index value should be obtained without any further restriction during the optimization procedure.

Optimization Method
Once an optimization problem is formulated carefully, it is required to solve the problem via an optimization method. To choose an optimization method to solve a problem, it is necessary to consider the number of optimization variables, the type of variables, the number of objective functions, the number of constraints, and the convexity or non-convexity of the problem and the other characteristics [27][28][29]. In this regard, the optimization methods can be classified into three major groups; (1) exact methods based on mathematical calculations, (2) heuristic methods, and (3) combination of the exact methods and heuristic methods.

Multi-Objective Optimization Using ε-Constraint Method
According to [33,34], the ε-constraint method is considered as one of the classic methods for multi-objective optimization. This method is in line with the exact methods. In addition to its efficiency and simplicity, this method is applicable to both convex and non-convex problems. The main idea of the ε-constraint method is to reformulate the multi-objective problem as a single-objective problem. Then, by iteratively solving the single-objective problem, a Pareto Front is obtained. In the following, the details of the ε-constraint method are explained [34].
Considering a multi-objective problem (Ψ(X)), as shown in Equation (39), subjected to different constraints that should be optimized, the following steps should be taken.
where f i (X) denotes the i th objective function and n shows the maximum number of existing objective functions.

1.
Each objective function ( f i (X)) is optimized with the existing constraints separately and the results are saved in a table, called the payoff table.

2.
According to the priority of the objective functions, one objective function is selected as the main objective function. Then, the rest of the objective functions are treated as new constraints and added to the main constraints. It should be noted that except for the main objective function, if the goal is to minimize and maximize all the objective functions, then, f i (X) ≤ e i and f i (X) ≥ e i , respectively. Also, e i is a variable parameter.

3.
In order to assign values to e i , the maximum ( f max i ) and minimum ( f min i ) values of each objective function should be considered, as shown in Equation (40). It should be noted that those values can be obtained from the payoff table.
To generate different values for e i,n i , Equations (41) and (42) are used to minimize and maximize the objective function, respectively. By dividing the domain of the i th objective function into q i equal parts using Equations (41) and (42), q i different values are obtained for e i,n i . It should be noted that n i denotes the number of available generated values for e i,n i .

5.
By using the obtained values from Step 4, it can be derived that f i (X) ≤ e i,n i or f i (X) ≥ e i,n i . For different values of e i , a set of solutions is obtained, which forms the Pareto front of the problem.
According to the above-mentioned descriptions, to solve the probabilistic multi-objective RPP, the following equation is formed.

Fuzzy Decision Maker (FDM)
As already mentioned, after solving a multi-objective optimization problem, a set of optimal solutions is obtained, called the Pareto Front. While only one solution from the Pareto Front can be chosen as the final optimal solution to the problem, which is known as the Best Compromise Solution (BCS). One way to choose the BCS is to use the Fuzzy Decision Maker (FDM). Having used a fuzzy membership function, each of the optimal solutions is mapped between 0 and 1. For the k th objective function, F k , the linear fuzzy membership is defined as Equation (45) [24], and it is supposed that all the objective functions are minimized.
whereF k represents the k th normalized objective function. In addition, F min k and F max k are used to express the minimum and maximum values of the k th objective function, respectively.
After obtaining the fuzzy values of each objective function using Equation (45), there are several ways to find the BCS. In this paper, to obtain the BCS, the min-max method, which is introduced in [35], is used. BCS = max min F 1 ,F 2 , . . . ,F k (46)

Simulation Results and Discussions
To evaluate the performance of the proposed ε-constraint method in the presence of various objectives, containing expected total VAR cost ( f 1 ), expected active power losses, expected voltage stability index ( f 2 ), and expected loadability factor ( f 3 ), two deterministic and three probabilistic cases are studied, as follows: A. Deterministic multi-objective RPP without considering the loadability factor (assessing the proficiency of ε-constraint method) B. Deterministic multi-objective RPP considering the loadability factor C. Probabilistic multi-objective RPP considering the load demand uncertainty D. Probabilistic multi-objective RPP considering the wind power generation uncertainty E. Probabilistic multi-objective RPP considering load demand and wind power generation uncertainties at the same time All the cases are implemented in GAMS environment Ver. 25.1.2 [36][37][38][39], and are solved using the CONOPT 3 Solver [40], in an ASUS laptop, with 8 GB of RAM and 2.4 GHz. The descriptions of the case study are presented in the next subsection.

Case Study Descriptions and Simulation Results
The test system is the IEEE 30-bus test system, which has 6 generation units, 4 transformers, and 41 branches. The initial settings of the generators' voltage magnitude and transformers tap settings are obtained from [30]. Figure 1 shows the single line diagram of the IEEE-30-bus test system. Also, both the output active and reactive power of generators are set according to [41]. The loads' data and line data are available in [42]. It is assumed that there is not any VAR source in the case study.
A. Deterministic multi-objective RPP without considering the loadability factor (assessing the proficiency of ε-constraint method) B. Deterministic multi-objective RPP considering the loadability factor C. Probabilistic multi-objective RPP considering the load demand uncertainty D. Probabilistic multi-objective RPP considering the wind power generation uncertainty E. Probabilistic multi-objective RPP considering load demand and wind power generation uncertainties at the same time All the cases are implemented in GAMS environment Ver. 25.1.2 [36][37][38][39], and are solved using the CONOPT 3 Solver [40], in an ASUS laptop, with 8 GB of RAM and 2.4 GHz. The descriptions of the case study are presented in the next subsection.

Case Study Descriptions and Simulation Results
The test system is the IEEE 30-bus test system, which has 6 generation units, 4 transformers, and 41 branches. The initial settings of the generators' voltage magnitude and transformers tap settings are obtained from [30]. Figure 1 shows the single line diagram of the IEEE-30-bus test system. Also, both the output active and reactive power of generators are set according to [41]. The loads' data and line data are available in [42]. It is assumed that there is not any VAR source in the case study. To allocate appropriately the VAR compensation devices, firstly, the -index should be determined for the load buses. Then, the load buses with high values of -index are taken into account as the candidate buses for the VAR compensation devices installation. It is observed that after implementing the proposed methodology, the load at bus 24, bus 25, bus 26, bus 29, and bus 30 obtain the higher values of -index than the other load buses. As a result, the VAR compensator buses are found. After the allocation of VAR compensation devices, it is supposed that the capacities of the VAR compensators can be set to zero. The initial conditions of the system considering the fullload and 1-year planning are stated in Table 2.  To allocate appropriately the VAR compensation devices, firstly, the L-index should be determined for the load buses. Then, the load buses with high values of L-index are taken into account as the candidate buses for the VAR compensation devices installation. It is observed that after implementing the proposed methodology, the load at bus 24, bus 25, bus 26, bus 29, and bus 30 obtain the higher values of L-index than the other load buses. As a result, the VAR compensator buses are found. After the allocation of VAR compensation devices, it is supposed that the capacities of the VAR compensators can be set to zero. The initial conditions of the system considering the full-load and 1-year planning are stated in Table 2. The control variables limits are expressed in Table 3. The per-unit energy cost is equal to 0.06 ($/h) [3], the fixed installation cost (e i ) for all VAR sources is 1000 ($), and the variable installation cost (C Ci ) for all VAR sources is 3.0 ($/kVAR) [44]. Table 3. Control variables limits.

Control Variable Value
The voltage magnitude at the load buses, which are considered as the state variables, must be limited between 0.95 (p.u.) and 1.05 (p.u.). To show the effectiveness of the proposed method, two deterministic and three probabilistic cases are considered in the following subsections.

Case A: Deterministic Multi-Objective RPP without Considering the Loadability Factor
In order to validate the efficiency of the proposed method for multi-objective RPP, the ε-constraint method is applied to a deterministic multi-objective RPP problem. The obtained results are also compared with the approach presented in [3]. The deterministic multi-objective RPP aims to minimize the total VAR cost and voltage stability index. Thus, it is expected to achieve a reduction in active power losses and an improvement in the voltage stability index. Table 4 shows the obtained results from deterministic multi-objective RPP without considering the loadability factor in IEEE 30-bus test system with the initial settings. The duration of the load for deterministic multi-objective RPP is assumed to be 8760 h for full-load condition and without changes in the load level. According to Table 4, 15 Pareto optimal solutions are obtained by the ε-constraint method. After that, the min-max approach chooses the fifth solution (highlighted row) as the BCS. The active power losses for the BCS are 4.9813 MW. Considering the same operating condition, the ε-constraint method is compared with the Multi-Objective Differential Evolution (MODE) algorithm, which is recommended to solve the deterministic multi-objective RPP [3]. For the BCS, the results of the comparison are presented in Table 5. As it can be observed from Table 5, the ε-constraint method shows better performance compared with the MODE algorithm in minimizing total VAR cost and active power losses. The superiority of the ε-constraint method is confirmed by a 6.2578 % reduction in total VAR cost and a 9.3815 % decrease in the active power losses over the Base Case. However, as it can be observed from Table 5, the voltage stability index of the conventional method is better than the proposed approach. The optimal values of the control variables for Case A are represented in Table 6. As it can be observed, only one VAR source has a value of zero. Note that the fixed installation VAR cost is also considered for all VAR sources with the value of zero during the planning studies. It is apparent from Case A that the ε-constraint is an effective method to generate Pareto optimal solutions for the multi-objective RPP. In this part, the ε-constraint method is applied to the multi-objective RPP problem in a complex form. However, the uncertainties of the load demand and wind power generation are not considered in this case. In comparison with Case A, another objective, which is called the loadability factor, is added to the problem. Therefore, the main objectives in this part include minimizing the total VAR cost, reducing the active power losses, improving the voltage stability index, and maximizing the loadability factor. In order to solve a deterministic multi-objective RPP considering the loadability factor, it is assumed that the system is under full-load condition. The duration of the load is assumed to be 8760 h. Table 7 provides the simulation results of the deterministic multi-objective RPP problem considering the loadability factor. In addition, this table illustrates that among the 15 generated Pareto optimal solutions, the eighth solution (highlighted row) is the BCS through the min-max approach. The active power losses for the BCS is 9.3494 MW. Also, it can be observed that the total VAR cost and active power losses are dramatically increased for BCS in comparison with Case A. Moreover, the voltage stability index is not improved compared with Case A. Nevertheless, the loadability factor is improved in Case B. The main reason behind the deterioration of active power losses and voltage stability index is due to the enhancement of the loadability factor. It should be noted that the loadability factor can hugely affect the active power losses and voltage stability of power systems.  Table 8 depicts the optimal values of the control variables for Case B. As it can be observed, three VAR sources have a value of zero. It should be noted that the fixed installation VAR cost is also considered for all VAR sources with the value of zero during the planning studies. In Cases A and B, the multi-objective RPP in power systems is solved using the deterministic approach. However, with the increasing level of uncertainty, probabilistic multi-objective is required for the RPP problem. In Case C, the probabilistic multi-objective RPP considering three different scenarios for the load level is performed. Each scenario for the load level consists of two main parts: (1) probability of the load level and (2) duration of the load. The overall duration of the load is assumed to be 8760 h, which is the expected time horizon for the RPP. The specifications of the system loading are described in Table 9. The simulation results obtained from the probabilistic multi-objective RPP considering the load demand uncertainty are given in Table 10. As it can be observed from this table, 15 Pareto optimal solutions are generated using the ε-constraint method. Using the min-max approach, the eighth solution (highlighted row) is chosen as the BSC. It is worth mentioning that the expected active power losses are 9.5049 MW for the BCS. From Table 10, it is clear that the expected total VAR cost is reduced compared with the Base Case. The expected voltage stability index and the expected loadability factor also show improvement towards the initial conditions. However, with more considerations, it is revealed that the expected active power losses are increased. This fact stems from the evident increase in the loadability factor. As a common incidence in power systems, following the escalation of the loadability factor, the active power losses increase and the system becomes voltage unstable. Generally, from the power systems operators' perspective, monitoring of the voltage magnitude at the load buses as a way of preventing voltage collapse is in high priority. Therefore, the voltage profile of the load buses for each loading scenario is plotted for the BCS, as shown in Figure 2.  The optimal values of the control variables for the BSC among the different load scenarios are given in Table 11. Having checked this table closely, it is noticed that some VAR sources have a value of zero. Consequently, the variable cost of those VAR sources is equal to zero. However, those VAR sources are allocated, and their fixed VAR installation costs are considered for the planning studies in this paper. It is clear from Figure 2 that the voltage magnitude of the load buses remains in the range of 0.95 p.u. and 1.05 p.u. for all three load scenarios. Although the loadability factor is improved, the voltage stability index of the system is ensured from the voltage magnitude point of view. Therefore, by making an allowance for the load demand uncertainty, the obtained total VAR cost seems to be more realistic. In the same way, the voltage stability index and the loadability factor are more reliable due to including more scenarios for the planning horizon. The optimal values of the control variables for the BSC among the different load scenarios are given in Table 11. Having checked this table closely, it is noticed that some VAR sources have a value of zero. Consequently, the variable cost of those VAR sources is equal to zero. However, those VAR sources are allocated, and their fixed VAR installation costs are considered for the planning studies in this paper. It is clear from Figure 2 that the voltage magnitude of the load buses remains in the range of 0.95 p.u. and 1.05 p.u. for all three load scenarios. Although the loadability factor is improved, the voltage stability index of the system is ensured from the voltage magnitude point of view. Therefore, by making an allowance for the load demand uncertainty, the obtained total VAR cost seems to be more realistic. In the same way, the voltage stability index and the loadability factor are more reliable due to including more scenarios for the planning horizon.

Case D: Probabilistic Multi-Objective RPP Considering the Wind Power Generation Uncertainty
In Case D, it is assumed that a wind farm is located at a PQ node. After generating the wind speed scenarios using the Weibull distribution and power curve, a probabilistic multi-objective RPP is performed. In this case, the IEEE 30-bus test system is modified based on [45]. Hence, a wind farm with a rated power of 40 MW is added to bus 22. The wind farm data is derived from [45] and is presented in Appendix A. To evaluate the impact of the wind farm, six scenarios for the output power of the wind farm are generated. The duration of the load is assumed to be 1460 h and without changes in the load level. The generated wind scenarios and their details are given in Table 12. The obtained results from the probabilistic multi-objective RPP using the generated wind scenarios are represented in Table 13. As it is observed from this table, among the 15 generated Pareto optimal solutions using the ε-constraint technique, the eighth solution (highlighted row) is selected as the BSC after applying the min-max approach. The corresponding value of expected active power losses is 8.5777 MW. Compared with the Base Case and Case A, it is clear that the expected total VAR cost has had a remarkable reduction for the best compromise solution. In addition, the enhancement of expected voltage stability index and expected loadability factor is undeniable towards the Base Case and Case A. In addition, the expected active power losses are elevated in contrast with the Base Case. However, the expected active power losses show a reduction of roughly 1 MW, when it is compared with Case A. The main reason for this reduction is the existence of the wind farm in the case study. Table 13. Obtained Pareto optimal solutions for the probabilistic multi-objective RPP considering the wind generation uncertainty. The optimal values of the control variables for the BCS over the wind scenarios are represented in Table 14. Taking a look at Table 14, it is shown that the VAR sources gain the value of zero in almost all scenarios. Therefore, the variable VAR investment cost for those VAR sources equals to zero. However, the VAR sources are allocated and their fixed VAR investment costs are taken into account during the planning horizon. The voltage profile of the load buses for the BCS is plotted in Figure 3. As it can be observed, the voltage magnitude of the load buses is kept in the interval of 0.95 p.u. and 1.05 p.u. during all wind scenarios. Hence, it can be concluded that based on a proper RPP and having adequate reactive power reserve, the voltage magnitude of the load buses are restricted with specific limits.  To have a more comprehensive overlook of the probabilistic RPP, in Case E, it is preferred to perform RPP in the presence of two stochastic input variables, including the load demand and wind power generation. Hence, considering the level of the load and wind power generation as the stochastic input variables, combined load-wind scenarios are generated via the proposed technique. After determining the load-wind scenarios, a probabilistic multi-objective RPP is performed to evaluate the existing objectives, while the number of random input variables increases. Taking the

Case E: Probabilistic Multi-Objective RPP Considering Load Demand and Wind Power Generation Uncertainties
To have a more comprehensive overlook of the probabilistic RPP, in Case E, it is preferred to perform RPP in the presence of two stochastic input variables, including the load demand and wind power generation. Hence, considering the level of the load and wind power generation as the stochastic input variables, combined load-wind scenarios are generated via the proposed technique. After determining the load-wind scenarios, a probabilistic multi-objective RPP is performed to evaluate the existing objectives, while the number of random input variables increases. Taking the IEEE 30 bus-test system with initial settings as the benchmark, 18 combined load-wind scenarios are generated. In general, three load scenarios and six wind scenarios are used to generate 18 combined load-wind scenarios. The descriptions of generated load-wind scenarios are given in Table 15.  Table 16 shows the obtained results for the probabilistic multi-objective RPP using the generated load-wind scenarios. As seen from Table 16, by applying the min-max method among the 15 generated Pareto optimal solutions using the ε-constraint approach, the eighth solution (highlighted row) is selected as the BCS. The associated value to expected active power losses is 8.5575 MW. It can be observed that the expected Total VAR cost is considerably reduced, while the expected voltage stability index and the expected loadability factor are not significantly improved towards Case B. In contrast with Case A and the Base Case, both the expected voltage stability index and the expected loadability factor are improved compared with case B. Considering the expected active power losses, no substantial decrease is observed in Case C when it is compared with Case B. Compared with Cases A and C, a reduction of about 1 MW in expected active power losses can be estimated. Due to enhancing the expected loadability factor, the expected active power losses escalate relative to the Base Case. The optimal values of the control variables among 18 generated load-wind scenarios for the BCS are represented in Table 17. As it can be observed from Table 17, the VAR sources gain the value of zero in most of the scenarios for the BCS. Hence, the fixed installation cost is calculated for the VAR sources that gain the value of zero.
In order to investigate the impact of VAR planning in bus voltage magnitude over the different load-wind scenarios, the voltage profile of load buses is plotted for each load-wind scenario in Figure 4. This Figure shows that the voltage magnitude of the load buses is limited to the range of 0.95 p.u. and 1.05 p.u. for all scenarios. As a result, the voltage magnitude of the load buses is regulated within the predefined limits. In order to investigate the impact of VAR planning in bus voltage magnitude over the different load-wind scenarios, the voltage profile of load buses is plotted for each load-wind scenario in Figure  4. This Figure shows that the voltage magnitude of the load buses is limited to the range of 0.95 p.u. and 1.05 p.u. for all scenarios. As a result, the voltage magnitude of the load buses is regulated within the predefined limits. In order to evaluate the impact of generated reactive power by wind farms on RPP, the reactive power of wind farms is taken into account during the planning process. Assuming the constant PF operation for the proposed wind farms, the generated reactive power is calculated for each scenario based on Table 18. The value of PF for the proposed wind farms is taken to be 0.98. It should be noted that the generated reactive power by the proposed wind farms is calculated using Equation (47), as follows: where , and , indicate the generated reactive and active power by the proposed wind farms, respectively. In order to evaluate the impact of generated reactive power by wind farms on RPP, the reactive power of wind farms is taken into account during the planning process. Assuming the constant PF operation for the proposed wind farms, the generated reactive power is calculated for each scenario based on Table 18. The value of PF for the proposed wind farms is taken to be 0.98. It should be noted that the generated reactive power by the proposed wind farms is calculated using Equation (47), as follows: Q W i,S = P W i,S tan cos −1 (PF) where Q W i,S and P W i,S indicate the generated reactive and active power by the proposed wind farms, respectively.  8.1223 Similar to the former case studies, after performing the probabilistic multi-objective RPP considering the reactive power injection by the proposed wind farms, 15 Pareto optimal solutions are obtained, as shown in Table 19. Table 19. Obtained Pareto optimal solutions for the probabilistic multi-objective RPP considering load demand and wind power generation uncertainties incorporating the generated reactive power by the proposed wind farms. From Table 19, it is clear that for the BCS, all the objectives are slightly improved towards Case E. Although this enhancement does not seem to be significant, it shows the penetration of generated reactive power by the proposed wind farms on planning studies. Moreover, the related active power loss reaches 8.4807 MW, which shows a reduction with respect to Case E. The optimal values of the control variables among 18 generated load-wind scenarios incorporating the generated reactive power by the proposed wind farms for the BCS are represented in Table 20. Table 20. The optimal values of the control variables for probabilistic multi-objective RPP considering load demand and wind power generation uncertainties incorporating the generated reactive power by the proposed wind farms.

Control Variable
It can be inferred from Table 3 that in almost all scenarios, the VAR sources gain the value of zero, except for the installed VAR compensator at bus 30. In addition, it is revealed that the expected value of the required VAR compensator device at bus 30 reduces while the generated reactive power of the hypothetical wind farms is taken into account. As a result, wind farms have the capability to participate in VAR planning. This leads to a reduction in the size and amount of VAR sources. Therefore, practical power systems show less desire to install new VAR support while numerous large-scale wind farms with sufficient generated reactive power are available. Table 21 compares the performance of Case E with other cases. As it can be observed, compared with Case B, the performance of Case E in terms of obtaining better values for f 1 , f 2 , and P loss is improved. Case E also shows better performance rather than Case C. All objectives are improved considerably compared with Case C. Note that, due to the presence of wind farms, f 2 is enhanced, while improving on the loadability index. In addition, Case E is compared with Case D, and it can be observed that all objectives are improved. f 1 is improved significantly. However, f 2 , f 3 , and P loss are not enhanced considerably. This is due to the fact that wind power generation uncertainty has a great impact on all objectives, which in both Case E and Case D are considered. While, the performance of Case F is better than Case E considering f 1 , f 2 , and P loss , and is slightly more than f 3 .

Conclusions
A multi-objective RPP in power systems considering load demand and wind power generation uncertainties to minimize reactive power investment cost, reduce active power losses, improve voltage stability level, and enhance loadability factor is presented in this paper. The ε-Constraint method is used to solve the probabilistic multi-objective RPP. For this purpose, using the L-index, the VAR compensation buses are found at the first stage. Then, to distinguish the exact difference between the deterministic and probabilistic VAR planning studies, five different cases are investigated. In order to test the efficiency of the proposed method, the IEEE 30-bus test system is implemented in GAMS software under five various conditions. The simulation results show that the proposed probabilistic multi-objective RPP considering load demand and wind power generation uncertainties is effective in reducing the VAR installation cost, improving the voltage stability of the system, and enhancing the loadability, simultaneously.
Author Contributions: All authors participated equally in conceptualization, methodology, implementation, and writing-review and editing the paper. Also, All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.