Probabilistic approach to investigate the impact of distributed generation on voltage deviation in distribution system

The increasing penetration of renewable-based generation in the distribution network adds randomness to the bus voltage due to its stochastic characteristics. Consequently, the distribution network faces voltage limit violations at multiple buses. In this scenario, the traditional voltage control method fails to maintain the voltage profile of the network due to its slow response time. A fast and dynamic voltage control technique can be achieved by utilizing the voltage-controlling capability of available DGs in the network with the information about the impact of DG-integrated bus on change in voltage at any bus in the network. Evaluating this information using the traditional probabilistic load flow method is computationally inefficient for real-time applications. Therefore, in this paper, a principal component analysis (PCA)-based novel probabilistic voltage sensitivity index (PVSI) is proposed. It ranks the DG-integrated buses by evaluating the impact of power perturbation at DG-integrated buses on changes in voltage at any bus in the network. The PVSI is analytically derived to reduce the computation time for ranking. Further, the proposed PVSI is validated on the 69-bus and the 141-bus distribution systems by comparing it with the traditional Monte Carlo simulation (MCS) and joint differential entropy (JDE) in terms of accuracy and computation time. From the simulation, it is observed that PVSI gives a similar ranking as both methods in less computation time. Due to the advantage in computation time, the proposed method enables the distribution system operator (DSO) to use this information for fast and dynamic voltage control in a real-time scenario.


Introduction
Nowadays, renewable-based distributed generation (DG) is being given priority for power generation due to environmental concerns. Consequently, the integration of renewable energy resources into the distribution system is increasing at a rapid pace [1]. It imposes voltage fluctuation in the distribution system due to its stochastic nature. As a result, overvoltage or undervoltage is observed at several buses in the distribution system (DS). The voltage limit violation in DS not only damages the electrical equipment but also reduces its hosting capacity [2]. Therefore, an effective volt-B Digamber Kumar 2018eez0002@iitrpr.ac.in Bibhu Prasad Padhy bibhu@iitrpr.ac.in 1 Electrical Engineering, IIT, Ropar, Rupnagar, Punjab 140001, India age control technique is required to keep the voltage level at all buses within their permissible limit. The traditional method of voltage control, such as on-load tap-changing transformer (OLTC) [3] and voltage regulator [4], have slow response time and are also not made for bidirectional power flow done by integrated DGs. With advancements in inverter technology, the voltage control using active and reactive power of inverter interfaced DGs in the network gives a fast and dynamic solution [5,6]. However, utilizing all available DGs or random DGs from the distribution system to maintain the voltage profile is not considered as a economical fair solution for the utility grid. In this regard, the information about the impact of power perturbation at the DG-integrated buses on the voltage deviation gives insightful information for maintaining the voltage profile within its permissible limits.
The voltage control technique is mainly classified into centralized and decentralized control [7]. In centralized control, all the voltage-controlling components are mutually coordinated. Thus, it requires an extensive communication infrastructure and sensors. The centralized voltage control utilizes all the voltage-controlling parameters simultaneously to maintain the voltage profile at all buses. Hence, it is a computational burden for a large distribution system [8,9]. In contrast, the decentralized voltage control method locally controls the voltage deviation. The decentralized voltage control technique reduces computational and communication burden in the distribution system [10,11]. However, both control methods do not directly consider the uncertainty present in the load and generation. Therefore, the solution suggested by both these methods is valid for specific scenarios of load and generations.
The stochastic nature of DGs impacts voltage deviation significantly. In [12], the impact of uncertainty in DGs output power on voltage profile has been discussed using probabilistic load flow, i.e., MCS, by P-V and Q-V curves. However, the voltage profile of the distribution system is not only influenced by uncertainty in the output of DGs but also by varying daily load, the structure of the network and the operating characteristics of voltage-compensating devices. A probabilistic approach has been utilized in [13] to determine the impact of uncertainty present in integrated DGs, random network configuration and operating characteristics of voltage-controlling devices on voltage profile. Determining the impact of uncertainty present in load and the generation on bus voltage is a computational process. The analytic method has been presented in [14] to evaluate the change in voltage due to power perturbation at multiple buses. The implementation of this analytic method is complex due to its formulation. It needs to be simplified for implementation in real-time applications. To attenuate the voltage limit violation due to the stochastic behavior of load and generation requires a fast and dynamic method of voltage control.
The method of voltage control, such as reinforcement of DS, gives a fast solution for maintaining the voltage profile [15], but increasing integration of DGs into DS, it becomes very costly. In this scenario, the voltage control using active and reactive power of inverter interfaced DG offers a costeffective solution. This control strategy is broadly classified into active power curtailment (APC) and reactive power compensation (RPC). In RPC, the inverter manages the reactive power demand by operating the DGs into leading or lagging power factor [16]. It makes the inverter similar to reactive power-compensating devices. RPC gives effective voltage control, but it has limited real-time application due to economic reasons. Additionally, RPC does not provide effective voltage control in the case of overvoltage phenomena in lowvoltage DS due to its low X/R ratio. In this case, active power curtailment from the DGs is considered a potential solution.
In [17], the overvoltage is mitigated by curtailing the active power from DGs through the centralized droop control method, and the amount of active power curtailment from DGs is done on the basis of their droop coefficient. However, all the DGs participate together to mitigate the overvoltage, due to which a large amount of active power is curtailed. Centralized voltage control with the objective of minimization of active power curtailment from DGs has been proposed in [18]. It has been implemented in two levels. In the lower level, a droop control has been implemented among the grouped grid-connected inverters. Further, using the available power information from the lower level, minimization of active power curtailment has been performed. In this work, optimal active power curtailment is performed from all the DGs in the network but utilizing all available DGs in DS for active power curtailment to mitigate the overvoltage issues is not an economical and fair solution for utility grids. The size of the distribution network is huge; therefore, it takes ample computation time to implement voltage control techniques using available DGs. A decentralized voltage control method is preferred to reduce the computation time and amount of active power curtailment from the available DGs in the network.
In [19], a double-layer decentralized voltage control strategy has been proposed. To implement this, DS is partitioned into several clusters using a voltage sensitivity-based index. Further, establishing the association between autonomous cluster optimization and inter-cluster coordination optimization, the optimal operating point of DGs in terms of active and reactive power has been determined. In this method, the optimal power flow problem is solved locally. Thus, it takes less computation time to determine the optimal operating point of DGs in comparison with any centralized voltage control method. However, all the DGs are participating together in power curtailment, which is not good for utility grids in terms of economy. For resolving unfair power curtailment from the DGs, in [20,21], the zonal voltage control method has been proposed. To implement this, DS is partitioned into several zones based on the modularity index, which is derived from voltage sensitivity. Further, to mitigate the overvoltage issues in the network, optimal power flow is solved to determine the optimal operating point of DGs for the cluster where overvoltage is observed. In this method to formulate the modularity index, a voltage sensitivity matrix is used, which is difficult to evaluate in highly radial structured DS due to load flow convergence. In [19], [20], and [21], stochastic characteristics have not been considered while partitioning the network. Therefore, the proposed clustering method is scenario specific. Evaluating the impact of DGs on bus voltage with their stochastic characteristics by traditional voltage sensitivity matrix is a computationally cumbersome process and also not applicable for highly radial structured DS. Therefore, the analytic method is required to evaluate the impact of DGs on bus voltage with considering the stochastic characteristics of DGs. In [22], differential entropy, Kullback-Leibler distance, and the Frechet distance index have been proposed to determine the dominant influencer of voltage fluctuation at any bus. It gives the impact of power perturbation at DG-integrated buses on voltage deviation. Evaluating this information by these indexes is computationally high process due to the complex analytic formula of change in voltage due to power perturbation at DG-integrated bus has been used to formulate these logarithmic indexes.
Therefore, to reduce the computation without hampering the accuracy, PCA-based, a novel probabilistic voltage sensitivity index (PVSI) has been proposed in this work to evaluate the impact of power perturbation at DG-integrated bus on bus voltage. The PVSI is analytically derived with simplified formula of change in voltage due to power perturbation at DG-integrated bus, due to which it ranks the impact of DGintegrated bus on bus voltage without using MCS. Thus, it takes very less computation time and also produces a similar result to the traditional MCS method. The highlights of the paper are as follows: • A PCA-based, novel probabilistic voltage sensitivity index (PVSI) is proposed. It ranks the impact of the DGintegrated bus on the voltage deviation at any bus. The DG-integrated bus, which has the highest value of PVSI for any randomly selected bus, is considered as the most impactful bus for the change in voltage over that randomly selected bus. With this information, centralized and decentralized voltage control can be made dynamic and economical. • Eigenvalues of the covariance matrix for linearly corelated variable data points give their variance along the principal component axis. Based on this observation, a novel PVSI is formulated, and in order to reduce computation time, PVSI is analytically derived with a simplified formula of change in voltage due to power perturbation at DG-integrated buses. • Efficiency of the proposed novel PVSI is investigated on 69-bus distribution and 141-bus distribution test systems by comparing with joint differential entropy (JDE) [22], and traditional MCS methods in terms accuracy and computation time.
This paper is organized as follows. In Sect. 2.1, the traditional method for determining the impact of power perturbation at DG-integrated bus on voltage deviation is discussed. A brief explanation of the analytic method for calculating the upper bound of voltage deviation is detailed in Sect. 2.2. The formulation of the PVSI using PCA is presented in Sect. 3. The validation of PVSI on 69-bus and 141-bus distribution systems is discussed in Sect. 4. The paper is concluded in Sect. 5.  Figure 1 shows a schematic of a generic active distribution system coupled with the main AC grid. The a ∈ {a 1 , a 2 ...a n } symbolize the active bus and o symbolize the observer bus.
The bus with distributed generations is taken as an active bus. Similarly, the bus which is selected for observing the impact of power perturbation at DG-integrated buses is termed an observer bus. If the complex power changes at active bus changes S a 1 to S a 1 + S a 1 , then voltage at observer bus becomes V o to V o + V o . A brief description of evaluating voltage sensitivity ( V o ) by traditional method using MCS and analytic methods is explained in subsections (A) and (B).

Traditional method for voltage sensitivity using Monte Carlo simulation (MCS)
The voltage sensitivity gives a direct relationship between power perturbation at DG-integrated buses with the change in bus voltage, as shown in (1) [23].
To obtain V i , dV i d P j , dV i dQ j must be acquired. Traditionally, Newton-Raphson (NR) and perturbation and observation (P&O) methods are used to determine dV i d P j , dV i dQ j .
In NR method, dV i d P j , dV i dQ j is determined by Jacobian matrix as in (2). When the Jacobian matrix is not accessible, the P&O method can be used. dV i d P j and dV i d P j are obtained through the multi-run simulation by means of perturbing the power at jth active bus and observing the change in voltage at ith observer bus, and then formulate it as shown in (3).
For considering the stochastic nature of load and generation, these two methods are applied with MCS to determine the uncertainty in the voltage deviation. It is computationally expensive; hence, it lacks practical implementation. In this paper, it is used to measure the accuracy of the proposed PVSI. Further, the accuracy study is measured in terms of a ranking index. The steps involved in ranking the impact of power perturbation at the active bus on the change in voltage at the observer bus with the MCS method are as follows: • Determine the variance of change in voltage at the observer bus due to power perturbation at active buses using the NR method with MCS. • Consider zero power perturbation for the active bus, whose effect has to be calculated on change in voltage at the observer bus. Repeat step one. • The active bus, which reduces more in the variance of change in voltage at the observer bus, is considered as the most impactful bus for the observer bus.

Analytic method for voltage sensitivity
The analytic method gives an upper bound of voltage deviation due to power change at active bus [14]. The upper bound of voltage deviation ( V o ) at observer has been shown in (4).
Z oa 1 is the common impedance shared between active bus a 1 and observer bus O. In Fig. 1, it is represented by the green line. It is determined from the impedance matrix (Z bus ), which is shown in (5). The off-diagonal element of (Z bus ) is considered as common impedance. To formulate the (Z bus ), impedance building algorithm is used [24]. V * a 1 is a complex conjugate of pre-perturb voltage at the active bus a 1 . If complex power changes at multiple buses, then voltage deviation is calculated at the observation bus by superposition theorem [14]. The expression for voltage deviation is shown in (6).
V o = − a n S a n Z oa n V * a n = a n V oa n (6) after simplifying (4) V r oa n = − 1 |V a n | ( P a n (R oa n cos θ a n − X oa n sin θ a n ) − Q a n (R oa n sin θ a n + X oa n cos θ a n )) ( Q a n (R oa n cos θ a n − X oa n sin θ a n ) + P a n (R oa n sin θ a n + X oa n cos θ a n )) where V r oa n and V i oa n are real and imaginary parts of change in voltage at observer bus O due to power perturbation at active bus a n . P a n and Q a n are active and reactive power perturbation at nth active bus. V a n , θ a n is pre-perturb voltage magnitude and angle of nth active bus. R oa n , X oa n is common resistance and reactance between observer bus and nth active bus. It is real and imaginary part of the common impedance Z oa n . Using (8) and (9), the real and imaginary part of the change in voltage due to power perturbation at any buses can be calculated. In general, the pre-perturb voltage at the active bus of the distribution system always remains within the permissible limit [25]. Therefore, |V a n | θ a n = 1 0 then, (8) and (9) becomes V r oa n = − P a n R oa n − Q a n X oa n (10) V i oa n = − Q a n R oa n + P a n X oa n (11) Equations (10) and (11) give a straightforward relationship between the active and reactive power perturbation with voltage deviation. Due to the stochastic nature of injected complex power by the renewable-base generation, it is highly important to define voltage deviation in a stochastic sense. Therefore, the mean, variance, and probability density function (PDF) of voltage deviation are the essential components to be obtained. For validating the stochastic properties of analytic formulation of change in voltage as shown in (10) and (11), 69-bus distribution system is considered.
Complex power perturbation is considered at twelve buses, from bus 2 to bus 69, in the interval of six buses. It is modeled by normal random vector N (0,σ 2 ), as shown in (12). In (12), the first and fourth elements of covariance matrix (σ 2 ) represent the variance of active (σ 2 p ) and reactive (σ 2 q ) power, respectively. The third and fourth elements represent covariance between active and reactive power (ρ pq σ 2 p σ 2 q ). DGs provide reactive power by curtailing active power. Hence, the active-reactive power correlation coefficient (ρ pq ) is taken as −0.1. The change in voltage at bus 27 is calculated using the analytic method. The mean, variance, and PDF of voltage deviation are determined using (A3) to (A6), given in "Appendix A". For evaluating the accuracy of the analytic method, PDF of the change in voltage at bus 27 is also determined by the NR method using MCS for 10,000 scenarios. In this case, mean and variance are estimated using the normal distribution mean, and variance formula [26]. After comparing both the PDF of V r 27 and V i 27 , it is observed that both methods give a similar PDF for change in voltage at bus 27, as shown in Fig. 2. But the analytic method reduces huge computation time.
In this paper, analytic method of calculating the change in voltage at the observer bus, as shown in (10) and (11), is used to formulate PV S I for ranking the impact of power perturbation at the active bus on change in voltage at observer bus.

Formulation of PVSI using PCA
The traditional method of finding out the impact of power perturbation at active buses on observer buses is computationally inefficient. Thus, in this paper, novel PVSI is proposed. The proposed method is based on principal component analysis (PCA).
PCA gives the variance of linearly correlated variable data points by reducing their dimensionality [27]. Linear correlation between the data points is determined by the correlation coefficient as shown in (13).
The maximum eigenvalue along the principal component axis of the covariance matrix for linearly correlated variable data points is considered as their variance. Considering this point, in this paper, PCA is used for calculating the vari-  (14) to (22).
In Fig. 3, transformation of data point ( v r oa 1 , v r o 1 ) along Z 1 and Z 2 is shown. Similarly, all the data points present on axis V r oa and V r o are transformed onto the new axis Z 1 and Z 2 and it is expressed in matrix form as follows: where, The Z 1 and Z 2 are orthogonal to each other. Therefore, it can be written as a T n a n = 1 n ∈ {1, 2} where A is the transformation matrix. Now, variance of data point along the Z n is calculated using (14) and (15): In (17) and (18), is covariance matrix for calculating the impact on real part of change in voltage at observer bus due to power perturbation at active bus. Its diagonal element represents variance and off-diagonal element represents covariance of variable V r oa and V r o , respectively. Similarly, covariance matrices and ( V oa , V o ) are formulated to determine the impact of power perturbation at active bus on imaginary and absolute part of change in voltage at observer bus. Elements of covariance is analytically derived as give in "Appendix A". Orthogonality of Z 1 and Z 2 is show in (15). Using (15) and (16), Lagrange multiplier is formulated for determining the maximum variance of variable data points ( V r oa , V r o ) along the new axis Z n , as shown in (20).
For determining the maximum variance of variable data points along the new axis Z n , Lagrange multiplier is differentiated in term nth row element (a n ) of transformation matrix A because the position of variable data points on new transformed axis is decided by (a n ). After simplifying (20), it is found that, (21) represents characteristic of covariance matrix ( V r oa , V r o ) . Where, λ represents the eigenvalue of (17) and (20), it is seen that the trace of the covariance matrix and summation of their eigenvalues are equal due to Z 1 and Z 2 , as shown in (22).
From (22), it is concluded that the maximum variance of data point along the principal component axis is the eigenvalue of the covariance matrix. Now a selection of principal component axis for determining the PVSI is done on the basis η h as given in (23). η h gives a fraction of the total variance shared by the principal component axis. The higher value of η h means that higher value data point lies along that axis.
On the basis of variance along that axis, the variance of variable data point can be estimated. In this work, the axis which has the highest value of η h that principal component axis is considered for evaluating PVSI.
where k represents active bus number, i.e., k ∈ {2, 4, 6, .......} To evaluate the impact of power perturbation on the real part of the change in voltage PVSI k r is formulated as shown in (24). Similarly, PVSI k i and PVSI k are formulated to determine the impact of power perturbation on imaginary and absolute parts of change in voltage, respectively. Algorithm 1 gives the steps to calculate the PVSI k r , PVSI k i and PVSI k for ranking the impact of power perturbation at active buses on V r o , V i o and V o , respectively.

Case study and discussion
To verify PVSI for ranking the impact of power perturbation at active buses on change in voltage at observer buses 69-bus and 141-bus distribution system is taken into account. The system data are given in [28,29]. The test system's nominal voltage and power ratings are taken as 12.47 kV and 100 MVA, respectively. The PVSI is capable of determining the impact of power perturbation at the active bus on voltage deviation at the observer bus with the arbitrary number of DGs at a random position. However, for simplicity, in both  the test system, the even number bus is selected as the active bus, and the far-end odd number bus from the grid has been selected as the observer bus. For finding out PVSI, Algorithm 1 is implemented in MATLAB R2019b with the system configuration of Intel-Core i7-8700, CPU @ 3.20GHz. It is verified in three scenarios on 69-bus and 141-bus distribution systems by comparing the results with the traditional MCS method as explained in section II-A and joint differential entropy (JDE) method [22].

Change in active power at multiple buses
The output power from DGs depends on weather conditions, due to which it is very important to gather information about the impact of active power perturbation on voltage deviation for the fair and dynamic solution of voltage control. The variance of power perturbation at any active bus depends on the size of DG at that bus. In this case, three different sizes of normal random variables are taken to model the active power perturbation at multiple locations in distribution systems, as shown in (25), (26) and (27). The power perturbation on available DGs in the distribution system is correlated. The correlation exists among them due to several factors, such as environmental factors like solar irradiance, wind speed or change in the price of electricity. For incorporating the correlation among the DGs for active power perturbation, the correlation coefficient (ρ p ) for active power perturbation has been taken as 0.2 [14].

69-bus distribution system
Active power perturbation ( P k ) is considered on thirtyfour active buses as shown in (25), (26) and (27), where k represents the active buses k ∈ {2, 4, 6.....68}, and it is shown in Fig. 4 by red color. Observing their effect on change in voltage, bus 27, bus 45 and bus 51 are selected as observer buses.

Calculation of PVSI in 69-bus distribution system for active bus 68
1. For finding out the impact of active power perturbation at bus 68 on real, imaginary and absolute part of change in voltage at observer bus 27, firstly covariance matrices ( V r Similarly, PVSI k r , PVSI k i and PVSI k are evaluated for thirtyfour active buses in 69-bus distribution system and for seventy active buses in 141-bus distribution test systems. On the basis of PVSI, the impact of power perturbation at active buses on the change in voltage at observer bus is ranked.

Discussion
Ranking of active buses for perturbation in active power are done using proposed PVSI, JDE [22] and the traditional MCS methods, as explained in section II-A for 10,000 scenarios. For measuring the accuracy of PVSI, JDE and the traditional MCS method are used. In Table 1, the ranking of an active bus for the 69-bus distribution system is presented. Similarly, the ranking of active buses for the 141-bus distribution system is presented in Table 3. In Table 1 and Table 3, the top ten active buses are arranged in descending order on the basis of their effect on the change in voltage at the observer bus. In the first column, ranking by the traditional MCS method; in the second column, ranking by the JDE method; and in the third column, ranking by PVSI have presented. The computation time for all three methods of ranking is given in Table 3. The traditional MCS method runs (k + 1) × 10, 000 times load flow for considering 10, 000 scenarios of power perturbation at active buses to rank the active buses on the basis of their impacts on the observer bus voltage. Here, k is the number of active buses. Therefore, it takes a very high computation time for ranking, which makes this method inefficient to use in a real-time scenario. The JDE takes relatively very less computation time than traditional MCS because it is analytically derived. However, the analytic formulation of change in voltage at the observer bus due to power perturbation at active buses used to formulate JDE requires an initial condition voltage profile at the active buses. And it is determined by load flow which increases the computation time. To formulate the PVSI, the analytic formula of change in voltage at the observer bus due to power perturbation at active buses does not require the initial condition of voltage profile at active buses. Therefore, it takes less computation time than the traditional MCS and JDE. After comparing PVSI with JDE and traditional MCS methods in terms of accuracy and computation time, it is observed that PVSI gives a similar ranking as the traditional MCS method in less computation time than JDE, as shown in Tables 1, 3 and 5. Therefore, it can be applied in a real-time scenario to determine the impact of active power perturbation on voltage deviation at any bus. In this case, PVSI is the maximum variance of a twodimensional data point  Fig. 5 shows PVSI values for all the active buses due to active power perturbation at active buses. Considering bus 27 as the observer bus of the 69 bus system, PVSI values are shown in Fig. 5a and considering bus 141 as the observer bus of the 141-bus system, the PVSI values are shown in Fig. 5b. In Fig. 5, at the Y-axis PVSI values of all the active buses and at the X-axis, active-buses k = 2, 3, 6... are represented. The brown line represents the value of PVSI k r for all the active buses, and it shows the impact of active power perturbation on the real part of the change in voltage. The yellow color represents the value of PVSI k i for all active buses, and it shows the impact of active power perturbation on the imaginary part of the change in voltage, and the blue line represents the value of PVSI k and it shows the impact of active power perturbation on the absolute change in voltage.
In Fig. 5, the highest value of PVSI k is observed at active bus 24. Therefore, in the 69-bus distribution system, bus 24 is selected as the most impactful bus for change in voltage at bus 27. From all the active buses, bus 24 has the highest impact on the voltage profile of bus 27 due to active power perturbation at active buses. Similarly, for the observer bus 45 and 51, the highest value of PVSI k is observed. Thus, bus 66 is selected as the most impactful bus for change in voltage at both the observer bus 45 and 51. From all the active buses, 66 has the highest impact on the voltage profile of both the observer buses 25 and 51 due to active power perturbation at active buses. In Fig. 5, the highest value of PVSI k is observed at bus 138. Thus, in the 141-bus distribution system, bus 138 is selected as the most impactful bus for change in voltage at bus 141. Similarly, for 59 and 129 the highest value of PVSI k observed at buses 54 and 126. Hence, buses 54 and 126 have the highest impact on change in voltage at observer buses 59 and 129 due to active power perturbation at active buses. PVSI k r gives the impact of power perturbation on the real part of the change in voltage while PVSI k i gives the impact of power perturbation on the imaginary part of the change in voltage. In Fig. 5, it is seen that the value of PVSI k i is more than PVSI k i for all active buses. Therefore, it is concluded that perturbation in active power has more influence on PVSI k r than PVSI k i at the observer bus.

Change in reactive power at multiple buses
With the development of electronic power inverters, inverterinterfaced DG itself offers reactive power compensation for maintaining the voltage profile. Study the impact of reactive power perturbation on voltage deviation; reactive power perturbation is considered at multiple buses. Their effect on change in voltage at the observer bus is determined based on the basis PVSI, JDE and the traditional MCS method. Perturbation in reactive power is modeled by a zero-mean normal random variable, as shown in (28), (29) and (30). For considering the correlation among reactive power perturbation among the available DGs, the reactive power correlation coefficient (ρ q ) is taken as 0.2.

Discussion
Ranking of active buses for perturbation in reactive power at active buses is done using the proposed PVSI, JDE [22] and traditional MCS methods for 10,000 scenarios. After comparing PVSI with JDE and traditional MCS methods, it is observed that the ranking due to perturbation in reactive power at active buses gives a similar ranking as active power perturbation at active buses, as shown in Table 1 and Table  3. Their computation time is given in Table 2. After comparing the ranking method PVSI with JDE and the traditional method, it is observed that PVSI gives a similar ranking as  Fig. 6 shows PVSI values for all the active buses due to reactive power perturbation at active buses. Considering bus 27 as the observer bus of the 69 bus system, PVSI values are shown in Fig. 6a and considering bus 141 as the observer bus of the 141 bus system, the PVSI values are shown in Fig. 6b.
In Fig. 6a, the highest value of PVSI is observed at bus 24. Therefore, in the case of the 69-bus distribution system, reactive power perturbation at bus 24 is selected as the most impactful bus for change in voltage at bus 27 due to reactive power perturbation at active buses. Similarly, The highest value of PVSI is observed at bus 66 for observer buses 45 and 51. Therefore, active bus 66 is selected as the most impactful bus for observer buses 45 and 51 due to reactive power perturbation at active buses.
In Fig. 6b, the maximum value PVSI is observed at bus 132. Therefore, in the case of the 141-bus distribution system, bus 132 is considered as the most impactful bus for change in voltage at bus 141 due to reactive power pertur- bation at active buses. Similarly, for the observer bus 59 and 129, the higher value of PVSI is observed at 54 and 126, respectively. Therefore, active bus 54 and 126 is selected as the most impactful bus for observer buses 59 and 129 due to reactive power perturbation at active buses. In Fig. 6, it is seen that the value of PVSI k i is more than PVSI k i for all active buses. Therefore, it is concluded that perturbation in reactive power has more influence on V i o than V r o at the observer bus.

Change in active and reactive power both at multiple buses
In this case, complex power perturbation is considered at multiple buses, and their effect on change in voltage at the observer bus is determined on the basis of PV S I , JDE and traditional MCS. Perturbation in complex power is modeled as bivariate normal random variable as shown in (31), (32) and (33). The active and reactive power of inverter-based distributed generation is negatively correlated. So, the correlation coefficient of active and reactive power perturbation

69-bus distribution system
Complex power perturbation ( S k ) is considered at thirtyfour active buses as shown in (31)

Discussion
Ranking of active buses for complex power perturbation at active buses is done using the proposed PVSI, JDE [22] and traditional MCS method. The ranking of active buses for complex power perturbation is the same as active power perturbation and reactive power perturbation at active buses, as shown in Table 1 and Table 3. The computation time is given in Table 4. After comparing all three methods in terms of accuracy and computation time, it is observed that the PVSI gives a similar ranking as traditional MCS in less computation time than JDE due to their simplistic and analytic formulation.
In this case, PVSI is the maximum variance of a twodimensional data point and ( V oak , V o ) along the principal component axis due to perturbation in complex power at active buses. Here, V r oak , V r oak and V oak are the real part of change in voltage, imaginary part of change in voltage and absolute change in voltage at observer bus O due to complex power perturbation at kth active buses and V r o , V i o and V o are the cumulative real part of change in voltage, imaginary part of change in voltage and absolute change in voltage at the observer bus due to complex power perturbation at the active bus together. Fig. 7 shows PVSI values for all the active buses due to complex power perturbation at active buses. Considering bus 27 as the observer bus of the 69 bus system, PVSI values are shown in Fig. 7a and considering bus 141 as the observer bus of the 141 bus system, the PVSI values are shown in Fig. 7b. Fig. 7a shows the maximum value of PVSI observed at bus 24. Therefore, in the case of the 69-bus distribution system,  Similarly, in Fig. 7b, the maximum variation in PVSI occurs at bus 132. Therefore, in the case of the 141-bus distribution system, complex power perturbation at bus 132 is considered as the most impactful bus for change in voltage bus 141. Similarly, for the observer bus 59 and 129, the higher value of PVSI is observed at bus 54 and 126, respectively. Therefore, active bus 54 and 126 is selected as the most impactful bus for observer buses 59 and 129 due to complex power perturbation at active buses. From all three cases, it is observed that the active bus near to observer bus, where the largest power perturbation is considered, has the highest impact on the change in voltage at that observer bus. It happens due to the impact of power perturbation on voltage deviation depending on the size and location of the perturbation. In Fig. 7, it is seen that the value of PVSI k r is more than PVSI k i , because perturbation of active power is dominant over reactive power perturbation due to its size.

Conclusion
Extracting information about the impact of power perturbation at the active buses on voltage deviation at observer buses by the traditional MCS method is a computationally cumbersome process. Hence, a principal component analysis-based novel PVSI is proposed, and the elements of the covariance matrix are analytically derived to reduce the computation time. The effectiveness of the proposed method is validated on the 69-bus and 141-bus distribution systems by comparing it with established JDE and traditional MCS methods. To evaluate the impact of power perturbation at active buses on voltage deviation at observer buses, active, reactive and complex power perturbations are considered at multiple active buses in three different cases. From all these cases, it is observed that PVSI gives a similar ranking as the traditional MCS and JDE methods in 85.75% and 85.8% less computation time than the JDE for 69-bus and 141-bus systems, respectively. Further, PVSI k r and PVSI k i have been used to determine the impact of active, reactive and complex power perturbations on the real part of the change in voltage ( V r O ),   ). However, the impact of a complex power depends on the variance of active and reactive power perturbation. In this work, the variance of active power is taken higher; therefore, complex power perturbation has a higher impact on the ( V r O ). Due to the advantage in computation time and accuracy, PVSI can be utilized in many applications, such as identifying the DG responsible for voltage fluctuation in the distribution system to enhance grid reliability and clustering the distribution network to implement the decentralized voltage control. One of the strengths of the proposed PVSI is that in the simulation, not only unidirectional but also bidirectional power flows are taken into account at several active buses. It is done by taking the variance of power perturbation more than the load at that active bus. Therefore, the ranking proposed by PVSI is valid for the distribution system where DGs are feeding power back to the grid.
Author contributions DK contributed to conceptualization, formal analysis, methodology, investigation, validation, visualization, and writing original draft. BPP contributed to supervision, resources, software, writing-review and editing.
Funding No funds, grants, or other support was received.

Conflict of interest
The authors have no conflicts of interest to declare that are relevant to the content of this article.

Formation of covariance matrix
From (10) var( V r oak ) = var( Q ak X oak − P ak R oak ) (A10) = X 2 oak var( Q ak ) + R 2 oak var( P ak ) − 2R oak X oak cov( P ak Q ak ) (A11) E( V r oak , V r oa j ) = E(( Q ak X oak − P ak R oak ) ( Q aj X oa j − P aj R oa j )) (A12) = E( Q aj Q ak X oak X oa j ) − E( Q aj P ak R oak X oa j ) −E( Q ak P aj R oa j X oak ) + E( P aj P ak R oak X oak ) = X oak X oa j E( Q aj Q ak ) − R oak X oa j E( Q aj P ak ) −R oa j X oai E( Q ak P aj ) + R oai X oak E( P aj P ak ) = X oak X oa j cov( Q aj Q ak ) − R oak X oa j cov( Q aj P ak ) −R oa j X oai cov( Q ak P aj ) + R oai X oak cov( P aj P ak ) (A13) From (11) var( V i oak ) = var(− Q ak R oak − P ak X oak ) (A14) = R oak var( Q ak ) + X oak var( P ak ) + 2X oak R oak cov( Q ak P ak ) (A15) E( V i oak , V i oa j ) = E((− Q ak R oak − P ak X oak ) (− Q aj R oa j − P aj X oa j )) (A16) R oai R oa j cov( Q ak Q aj ) + R oa j X oak cov( P ak Q aj ) X oa j R oak cov( P aj Q ak ) + X oa j X oak cov( P ak P aj ) (A17) E( V r oak , V i o ) = E(( Q ak R oak − P ak X oak ) (− Q aj R oa j − P aj X oa j )) (A18) −X oak R oa j cov( Q ak Q aj ) + R oa j R oak cov( P ak Q aj ) −X oak X oak cov( P aj Q ak ) + X oak R oa j cov( P ak P aj ) E( V i oak , V r o ) = −X oa j R oak cov( Q ak Q aj ) +R oa j R oak cov( P aj Q ak + X oak R oa j cov( P ak P aj ) −X oak X oa j cov( P ak Q aj ) (A19)