# Probabilistic Load Flow of an Islanded Microgrid with WTGS and PV Uncertainties Containing Electric Vehicle Charging Loads

Authors: Neeraj Gupta; B. Rajanarayan Prusty
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

With the emphasis of utilities on decentralized power systems due to the considerable growth of renewable energy resources (RES), microgrids (MGs) as controllable small grids are need of an hour. For an objective analysis of microgrids with large grid parity of renewable energy sources, probabilistic methods are required along with the already established deterministic methods. A Gauss quadrature-based probabilistic load flow (GQPLF) method for microgrid has been proposed along with intermittent generation and uncertain loads, including electric vehicles (EVs). The distributed generation consists of wind turbine generation systems (WTGS) and photovoltaic systems (PVs) connected to 14 bus microgrid systems with six droop-controlled distributed generators. Extensive simulations have been carried out on a 14-bus microgrid with distributed generation and uncertain loads, and results are compared with the benchmark Monte Carlo simulation (MCS) method, which validates the efficacy of the proposed method. The time taken by the proposed GQPLF for the 14-bus microgrid system simulation is 60.9 seconds, which is significantly less compared to the benchmark Monte Carlo simulation method, in which time taken for a 14-bus microgrid system with 10,000 simulations is 9318.7 seconds. Comparing the proposed GQPLF method with the benchmark method validates its efficacy and efficiency.


---

## Body

## 1. Introduction

The significant grid parity of renewable energy systems with the increase in intermittent generation and aggregation of loads and generating units at low and medium voltage levels leads to the development of microgrids that can operate in self-control mode in a decentralized grid environment [1, 2]. A microgrid is a group of different power sources that work together as a single system that can be controlled to provide reliable power to a local area [3]. Being the promising paradigm of power systems, the microgrid integrates the controllable loads with distributed generation viz., solar and wind, along with the battery energy storage system (BESS) for the maintenance of power supply [4].Microgrids can be operated in grid-connected or islanded mode to provide an uninterruptible power supply during disturbances to critical loads, which can be realized by using an appropriate control scheme. The control scheme automatically provides a signal to detach from the main grid during contingencies, with sufficient capacity to meet the critical loads [4–6]. A central controller governs the bidirectional flow of power at the point of common coupling (PCC), which depends on whether the microgrid supplies power or receives power from the grid [7–9]. In the case of an islanded microgrid, distributed generation (DG) droop control techniques are used to meet load demands while maintaining the system frequency and bus voltages in a decentralized mode. In some cases, in addition to droop-control DG units, pv and pq control strategies have also been implemented in an islanded microgrid [10]. In the grid-connected mode, the frequency regulation of microgrid by the utility is accomplished through PCC and for DG units, and pv and pq control strategies are used [11].Load flow analysis is needed for an autonomous microgrid to find out if the supply from DGs is enough without affecting the voltage profile and to figure out the state of the system [2]. Power flow studies are also required for optimal operation of conventional power systems and expansion planning, which are equally beneficial in the case of microgrid operation. In the literature, different load flow techniques are divided into direct methods, based on the Newton–Raphson (NR) method and methods based on a backward sweep [2]. The direct methods’ effectiveness and convergence criteria are determined by the number of nodes and lines in the impedance matrix. Specific techniques employ the NR method to compute bus voltages and power flows in distribution networks, while the radial method employs the radial topology of distribution systems. A loss allocation strategy based on the decomposition of branch currents has been described in [12], while [13] proposed a rigid power flow method based on a series impedance model. A nondivergent study of load flow based on the NR technique is proposed in [14]. By [15], the NR-based technique was extended to unbalanced systems. Most NR-based techniques exhibit a rapid convergence rate but ineffective in exploiting the system topology. In the sweep-based approaches, the backward, forward sweep-based load flow analysis was proposed in [16], the modified backward and forward sweep method was introduced in [17], and the three-phase radial distribution systems were extended in [18]. Due to the small size of DGs and inability to be used as infinite bus or slack bus, conventional load flow methods like Newton–Raphson, Gauss–Seidel cannot be used in the case of steady-state analysis of microgrid [19, 20].Towards this, various power flow methods for microgrids have been proposed by various researchers. The distributed slack bus model was extended with a Newton iterative approach, and the effect of the participation of generators concerning loads and losses was examined in [21]. The notion of the modified slack bus based on the participation factor of generators is described in [22]. The concept of “domains” and “commons” of individual generators in a microgrid is presented in [23]. To solve the load flow problem of an islanded microgrid, a modified conventional load flow has been proposed in [1]. The slack bus issues and frequency dependency have been taken care of in [24, 25] in which DGs droop characteristics have been applied for an islanded microgrid. Due to the significant penetration of wind energy systems and photovoltaic generation, as well as the rising unpredictability of power systems, there has been a greater emphasis on uncertainty analysis of power systems [26, 27].For achieving these objectives, in the literature, various probabilistic load flow methods have been proposed for distribution and transmission systems which can be numerical or analytical [28, 29]. In the Monte Carlo simulation (MCS)-based numerical method, using an system model, the input parameters as random variables are sampled, and output parameters are estimated from each sample [28, 29]. This method is used as a benchmark for probabilistic analytical-based power flow methods, whose accuracy is guaranteed by probabilistic limit theory [30]. The points and weights are obtained in analytical based methods using the point estimate method, and output moments are obtained with a weighted sum of cumulants of input variables. Various researchers have proposed the extended efficient variants of analytical-based PLF-like PLF-based convolution method in [31], Fourier transform-based convolution in [32], extended PLF based on PEM in [33] and 7 PEM in [34].With the inclusion of uncertain wind generation, extended PEM and discrete PEM have been proposed in [35] and [36], respectively, as different variants of PEM. Nataf transformation-based PEM has been proposed in [32]. With wind and PV uncertainties and to handle the ambiguities in the distribution system, an unsymmetrical two point estimate method has been proposed in [37]. With beta-PV distribution, this method has been extended to 5-point estimate for the transmission system in [38]. For quantifying the wind-based uncertainties, Rayleigh and Weibull distributions have been used in [32, 33]. Queuing theory has been utilised in [39] to characterise EVs’ overall charging and discharging power for probabilistic restricted power flow under unpredictable EV charging loads. Based on queuing theory [39] and the daily recharge energy of electric vehicles, [40] suggests a statistical analysis to find the best way to meet the charging needs of electric vehicles.However, the majority of the methods in the literature concerning microgrid are deterministic, taking into account worst-case scenarios, but for power system analysis, particularly for microgrid, where these unpredictable sources have greater share, uncertainty quantification methods are also necessary in addition to deterministic methods. In addition, the penetration of solar, wind, and electric vehicles in the grid has increased in recent years, and numerous researchers are exploring this area and proposing diverse models [41–44]. Using probabilistic methods, a system operator can get an alternate view of different output parameters and the likelihood that they will happen. In the case of microgrids, significant grid parity of uncertain resources necessitates using these methods, which is the primary motivation behind the current work.Also, the research mentioned above suggests that an uncertainty analysis needs a PLF-based algorithm that is easy to compute and gives a reasonable estimate in a small number of points when the system is unpredictable because of the high number of WTGS, PVs, and EVs on the grid. Since the input parameters, i.e., active and reactive power loads, WTGS power output, and PV power output, are random variables with a known PDF, only a probabilistic technique is adequate for quantifying this type of uncertainty [45].To address the above issues, this paper proposes an analytical Gauss quadrature-based probabilistic flow method (GQPLF) for islanded microgrid probabilistic analysis. The proposed method has been implemented using the modified Gauss–Siedel deterministic load flow along with wind, solar, and EV uncertainties. The contributions of this paper are summarized as follows:(1)
A novel analytical Gauss quadrature-based probabilistic flow method has been proposed for an islanded microgrid which can handle large grid parity of uncertain resources and uncertainty(2)
The proposed Gauss quadrature method has been integrated with different kind of uncertainties, i.e., wind, solar, including EV by taking their distribution functions with droop-controlled generators(3)
The proposed method has been tested on a 14-bus microgrid system and compared with the benchmark numerical Monte Carlo simulation method with large grid parity of uncertain resources under different scenariosThe organization of paper is as follows. The microgrid basic load flow is explained in Section2. In Sections 3 and 4, the Gauss quadrature-based load flow and uncertainty characterization of source and loads has been discussed. Lastly, results and conclusions of the work are presented in Sections 5 and 6, respectively.


## 2. Microgrid Load Flow

In the case of an islanded microgrid, the system frequency is variable, which leads to the variable load demand, reactance, and bus admittance matrix. The load demand is represented using the exponential function (as this pattern is followed in majority of the loads), which in-turn depends on frequency and voltage which is given by [1, 5](1)PLk=PLkoVkVoα1+kpω−ωo,QLk=QLkoVkVoβ1+kqω−ωo,where Vo is the nominal voltage, deviation in angular frequency is represented by ω−ωo, the active and reactive power corresponding to the nominal operating voltage are represented by PLko and QLko, respectively, the frequency sensitivity parameters are represented by kp and kq, and the exponents for different type of loads are represented by α and β.The bus admittance matrix is represented byYknω, which in-turn depends on the frequency, and as in the case of a microgrid, frequency is variable, so Yknω is given by [1](2)Yknω=Y11ωY12ω…Y1nωY21ωY22ω…Y2nω⋮⋮⋮⋮Yk1ωYk2ω…Yknω.In autonomous microgrids with no help from the electrical distribution grid, droop control is the primary means of distributing demand power across generators. DGs can act as droop-controlled buses in an islanded microgrid due to the absence of infinite bus or slack bus, and is represented as(3)V=Vo−nqQG,ω=ωo−mpPG,where V and ω are the voltage and frequency, respectively; nq and mp are the voltage and frequency droop coefficients, respectively; and PG and QG are the generated active and reactive power, respectively.For the load flow analysis of microgrid, the Gauss–Siedel (GS) method has been modified with initial values of frequency and voltage as 1 per unit and1∠0, respectively. The main procedure involved in power flow analysis is as follows [1]:(1)
The vector with variable parameters (voltage and frequency) is represented asX=VTωT(2)
The voltages at PQ buses can be calculated by solving the conventional GS method (to obtain the initial guess) as(4)Vki+1=1YkkPk−jQkVki∗−∑n=1k−1YknVni+1−∑n=k+1NYknVni.Here,Pk and Qk are the net active and reactive powers at the kth bus, respectively, and Vki+1 is the voltage at the kth bus for iteration i+1th iteration.(3)
The reactive power violation limits are checked at the PV bus (by calculating the reactive power at all the buses) using the following equation:(5)Qki+1=−imagVki∗∑n=1k−1YknVni+1+∑n=kNYknVni.(4)
After retaining the angle, the magnitude of a PV bus is reset to the prespecified value after computing the voltage using equation (4)(5)
Using the droop equations, the power (active and reactive) for frequency dependent voltage buses are calculated using equation (3) as(6)PGki+1=1mpkωo−ωi,QGki+1=1nqkVo−Vki.Limit violations are verified, and if they are, the limits are reset, and lastly, the voltage is determined.(6)
The system’s active power losses are computed.(7)
All the DGs will operate as droop-based DGs, with active power supplied by all the droop buses in the microgrid at same angular frequency. Consequently, the active power generation of microgridPsystem is calculated, which in-turn corresponds to the sum of all the active powers of the DGs. Finally, Psystem is the addition of active power demand and active power loss and is calculated as(7)Psystem=∑k=1dPGk=∑k=1d1mpkωo−ωi,PLoad+PLoss=∑k=1dPGk=∑k=1d1mpkωo−ωi.(8)
Now, updated values of frequency is obtained after calculation of updated power values, so equation (7) is rearranged to calculate the updated frequency as(8)ωi+1=∑k=1d1/mpkωo−PLoadi+1+PLossi+1∑k=1d1/mpk.(9)
SinceYbus is a function of frequency, so finally, Ybus is updated with the updated value of frequency in each iteration.(10)
Consequently, the load flow is solved, voltages are calculated, the error is evaluated, and the line flows are calculated if the convergence criterion has been met as shown in Figure1.Figure 1
Flowchart of microgrid load flow.


## 3. PLF Based on Gauss Quadrature

In the present work, the statistical analysis of microgrid has been carried out using the Gauss quadrature-based estimate method, which is generally applied to calculate statistical moments of any RV [29]. The main objective of this method is to calculate the output bus voltage distribution functions from the input probability density functions (PDFs) of bus voltages and active and reactive power load distributions. The estimation of weights and points on the PDF has been successfully obtained using the proposed method, and a complete PDF has been created as a result (from those points).Letμl and σl be the mean and standard deviation of random variable xl, respectively, of PDF fl. The lth random variable of xl can be obtained by varying ll=1,2,…n. The moments are the quantitative measurements associated with the form of a function and by calculating the first moment about the origin and square-root of the mean’s second central moment, and μl and σl can be obtained, respectively. Further, the kth estimated point of xl, as denoted by xl,k, is estimated by the Gauss quadrature points and weights method.


### 3.1. Calculation of Weights and Points with Generalized Gauss Quadrature

The distribution functions associated with various uncertainties (load and generation) are obtained, and the subsequent steps involve calculating the distribution function’s associated points and weights. The Gauss quadrature method provides the abscissaξ and weights ωi in the form of equation (9) for the ∫−∞∞e−ξϕξdξ type normal integrals for some function ϕ [46] as(9)∫−∞∞e−ξϕξdξ≈∑i=1nωiϕξi.The expectation offX for a random variable (RV) X, which is normally distributed with mean μ and variance σ2, can be obtained by applying the transformation ξ=x−μ/2σ,ϕξ=fμ+2σξ/π and dξ=dx/2σ2 as [28](10)∫−∞∞12πσ2e−x−μ2/2σ2fxdx≈∑1nωiπfμ+2σξi.Similarly, weightswol,k and points xl,k for any distribution function can be obtained.The following algorithm has been used to make the weights representative of the whole of the function.(1)
Let us consider active and reactive power loads asn random variables, i.e., xll=1,2,…n(2)
Obtain the summation of the random variables which is represented asSxl(3)
To calculate the weightage for each random variable, the following expression has been usedΠx1=x1/Sxl,Πx2=x2/Sxl,….,Πxn=xn/Sxl(4)
Finally, the weightswl,k for RV are obtained as w1,k=Πx1×wo1,k,w2,k=Πx2×wo2,k,….,wl,k=Πxl×wol,k


### 3.2. Gauss Quadrature-Based PLF for the Microgrid

The following procedure has been adopted for microgrid probabilistic load flow after weights, and points are calculated using the quadrature rules (as shown in Figure2) [29]:(1)
Form the input matrices according to the equation (11), which contains the estimated points of load and generation and the mean values of load and generation as(11)Xk=x1,kμx2…μxnμx1x2,k…μxn⋮⋮⋱⋮μx1μx2…xn,k.Here, loads (active and reactive power) are represented asxll=1,2,…n, μx1 as the mean of lth random variable, and xl,k as the kth point. The value of k=1,..,m, where m=3,5,10,12 for 3,5,10,12 point quadrature.(2)
A deterministic power flow (as explained in Section3.1) is carried out for each row of Xk(3)
While simulating the GQPLF for each point, the system has been checked for any constraint violation, i.e. if the reactive power at the generator bus is violated, that bus has been considered as the load bus, and the updated values of voltage have been calculated, which is used in next iteration and the system will not converge. However, the weight will not be affected by any constraint violations.(4)
As all points needs to be evaluated, by repeating step 2 for all the rows of the matricesX1X2,..,Xk, a total of ‘nm’ load flow computations would be performed(5)
For each output variable of interest,yi,lk, calculate the jth moment (quantitative measure which provides the form of a distribution function) as (12)Eyi,lkj=∑l=1n∑k=1mwl,kEyi,lkj,j=1,…,no.of  moments.Here,yi,lk is the value of ith variable of interest corresponding to lkth load flow with l=1,…n and k=1,…m. In the present work, first eight moments have been considered.(6)
Finally, by using the above obtained parameters (moments) and appropriate expansion series, the distribution functionyi can be obtained [28].Figure 2
Gauss quadrature-based load flow of an islanded microgrid.


## 3.1. Calculation of Weights and Points with Generalized Gauss Quadrature

The distribution functions associated with various uncertainties (load and generation) are obtained, and the subsequent steps involve calculating the distribution function’s associated points and weights. The Gauss quadrature method provides the abscissaξ and weights ωi in the form of equation (9) for the ∫−∞∞e−ξϕξdξ type normal integrals for some function ϕ [46] as(9)∫−∞∞e−ξϕξdξ≈∑i=1nωiϕξi.The expectation offX for a random variable (RV) X, which is normally distributed with mean μ and variance σ2, can be obtained by applying the transformation ξ=x−μ/2σ,ϕξ=fμ+2σξ/π and dξ=dx/2σ2 as [28](10)∫−∞∞12πσ2e−x−μ2/2σ2fxdx≈∑1nωiπfμ+2σξi.Similarly, weightswol,k and points xl,k for any distribution function can be obtained.The following algorithm has been used to make the weights representative of the whole of the function.(1)
Let us consider active and reactive power loads asn random variables, i.e., xll=1,2,…n(2)
Obtain the summation of the random variables which is represented asSxl(3)
To calculate the weightage for each random variable, the following expression has been usedΠx1=x1/Sxl,Πx2=x2/Sxl,….,Πxn=xn/Sxl(4)
Finally, the weightswl,k for RV are obtained as w1,k=Πx1×wo1,k,w2,k=Πx2×wo2,k,….,wl,k=Πxl×wol,k


## 3.2. Gauss Quadrature-Based PLF for the Microgrid

The following procedure has been adopted for microgrid probabilistic load flow after weights, and points are calculated using the quadrature rules (as shown in Figure2) [29]:(1)
Form the input matrices according to the equation (11), which contains the estimated points of load and generation and the mean values of load and generation as(11)Xk=x1,kμx2…μxnμx1x2,k…μxn⋮⋮⋱⋮μx1μx2…xn,k.Here, loads (active and reactive power) are represented asxll=1,2,…n, μx1 as the mean of lth random variable, and xl,k as the kth point. The value of k=1,..,m, where m=3,5,10,12 for 3,5,10,12 point quadrature.(2)
A deterministic power flow (as explained in Section3.1) is carried out for each row of Xk(3)
While simulating the GQPLF for each point, the system has been checked for any constraint violation, i.e. if the reactive power at the generator bus is violated, that bus has been considered as the load bus, and the updated values of voltage have been calculated, which is used in next iteration and the system will not converge. However, the weight will not be affected by any constraint violations.(4)
As all points needs to be evaluated, by repeating step 2 for all the rows of the matricesX1X2,..,Xk, a total of ‘nm’ load flow computations would be performed(5)
For each output variable of interest,yi,lk, calculate the jth moment (quantitative measure which provides the form of a distribution function) as (12)Eyi,lkj=∑l=1n∑k=1mwl,kEyi,lkj,j=1,…,no.of  moments.Here,yi,lk is the value of ith variable of interest corresponding to lkth load flow with l=1,…n and k=1,…m. In the present work, first eight moments have been considered.(6)
Finally, by using the above obtained parameters (moments) and appropriate expansion series, the distribution functionyi can be obtained [28].Figure 2
Gauss quadrature-based load flow of an islanded microgrid.


## 4. Uncertainty Characterization of Sources and Loads

### 4.1. Modelling of WTGS

The power from intermittent wind can be extracted using wind turbine generator system (WTGS). For the active power output of the generator, for a given wind speedVw, a quadratic power curve is used, which is usually provided by the manufacturer (equation (13) [28]).(13)P=0Vw≤VcPRVw2−Vc2VR2−Vc2Vc<Vw≤VRPRVR<Vw≤VF0VF<Vw.In equation (13), VR,Vc and VF are turbine’s rated, cut-in and cut-out speeds, respectively, and PR is the turbine’s rated power.To calculate the consumed reactive power by WTGS, a suitable model of induction generator can be solved, and in this work, a doubly fed induction generator (DFIG) model has been used.DFIG consists of wound rotor induction generator that is connected to the converter through slip rings. Furthermore, by manipulating the converter,±30% of the speed variation around the synchronous speed, as well as active and reactive powers, may be regulated. The reactive power absorbed or consumed in this case depends on the power factor and can be calculated as [29](14)Q=P1−cos2θcosθ.


### 4.2. Modelling of the Photovoltaic System

The PV module’s output power is a function of solar irradiance, module characteristics and sites’s ambient temperature. As observed from statistical studies, irradiance of solar has the Beta distribution which is represented as [28](15)PDFpvs=Γα+βΓα+Γβsα−11−sβ−1,if0≤s≤1,0≤α,β,0,else,,where α, β are the beta-PDF parameters (obtained from mean value and standard deviation of beta distribution) and s is the solar irradiance in kW/m2.The output power is determined by calculatingIs and Vs using the following equation from the generated solar irradiance having beta distribution [29]:(16)Ppvs=N×FF×Vs×IsFF=VMPP×IMPPVoc×Isc,Vs=Voc−Kv×TcIs=sa×Isc+KiTc−25Tc=TA+sa×NOT−200.8,where sa is the average solar irradiance, NOT is the nominal operating temperature of PV cell in °C, TA is the ambient temperature in °C, Tc is the cell temperature in °C, Ki and Kv are current and voltage temperature coefficients in A/°C and V/°C respectively, Isc is the short circuit current in Amperes, Voc is the open circuit voltage in Volts, FF is the fill factor, and VMPP and IMPP are the voltage and current at maximum power point in Volts and Amperes, respectively.


### 4.3. Charging Load of Electric Vehicle

The charging load of an electric vehicle is determined by its daily driving range, operating status, and battery capacity.The EV’s recharge power load follows the Weibull distribution. In [28, 29], the KS goodness-of-fit test (Kolmogorov–Smirnov) was used to test and postulate, and the same distribution is utilised in this work:(17)fxev|Φ,φ=φΦxevΦφ−1e−xevΦφ,xev≥0,0,xev<0,where φ,Φ are the parameters (scale and shape) of Weibull distribution, respectively.


## 4.1. Modelling of WTGS

The power from intermittent wind can be extracted using wind turbine generator system (WTGS). For the active power output of the generator, for a given wind speedVw, a quadratic power curve is used, which is usually provided by the manufacturer (equation (13) [28]).(13)P=0Vw≤VcPRVw2−Vc2VR2−Vc2Vc<Vw≤VRPRVR<Vw≤VF0VF<Vw.In equation (13), VR,Vc and VF are turbine’s rated, cut-in and cut-out speeds, respectively, and PR is the turbine’s rated power.To calculate the consumed reactive power by WTGS, a suitable model of induction generator can be solved, and in this work, a doubly fed induction generator (DFIG) model has been used.DFIG consists of wound rotor induction generator that is connected to the converter through slip rings. Furthermore, by manipulating the converter,±30% of the speed variation around the synchronous speed, as well as active and reactive powers, may be regulated. The reactive power absorbed or consumed in this case depends on the power factor and can be calculated as [29](14)Q=P1−cos2θcosθ.


## 4.2. Modelling of the Photovoltaic System

The PV module’s output power is a function of solar irradiance, module characteristics and sites’s ambient temperature. As observed from statistical studies, irradiance of solar has the Beta distribution which is represented as [28](15)PDFpvs=Γα+βΓα+Γβsα−11−sβ−1,if0≤s≤1,0≤α,β,0,else,,where α, β are the beta-PDF parameters (obtained from mean value and standard deviation of beta distribution) and s is the solar irradiance in kW/m2.The output power is determined by calculatingIs and Vs using the following equation from the generated solar irradiance having beta distribution [29]:(16)Ppvs=N×FF×Vs×IsFF=VMPP×IMPPVoc×Isc,Vs=Voc−Kv×TcIs=sa×Isc+KiTc−25Tc=TA+sa×NOT−200.8,where sa is the average solar irradiance, NOT is the nominal operating temperature of PV cell in °C, TA is the ambient temperature in °C, Tc is the cell temperature in °C, Ki and Kv are current and voltage temperature coefficients in A/°C and V/°C respectively, Isc is the short circuit current in Amperes, Voc is the open circuit voltage in Volts, FF is the fill factor, and VMPP and IMPP are the voltage and current at maximum power point in Volts and Amperes, respectively.


## 4.3. Charging Load of Electric Vehicle

The charging load of an electric vehicle is determined by its daily driving range, operating status, and battery capacity.The EV’s recharge power load follows the Weibull distribution. In [28, 29], the KS goodness-of-fit test (Kolmogorov–Smirnov) was used to test and postulate, and the same distribution is utilised in this work:(17)fxev|Φ,φ=φΦxevΦφ−1e−xevΦφ,xev≥0,0,xev<0,where φ,Φ are the parameters (scale and shape) of Weibull distribution, respectively.


## 5. Results and Discussion

For testing and checking the proposed method’s efficacy, the Gauss quadrature-based microgrid load flow has been tested with WTGS and PV uncertainties containing EV loads. The studies have been conducted on a 14-bus microgrid system containing droop-controlled generators with WTGS, PV, and EV.


### 5.1. Description of the System

The test system of an islanded microgrid under consideration consists of 14 bus systems with two buses having wind turbines and two buses containing photovoltaic systems, as shown in Appendix A. EV load is connected across one bus with six buses containing droop-controlled DGs.In this system, droop controlled DGs are connected at buses5,6,8,12,13,14 with wind generation systems at buses 1,2 and photovoltaic generation systems at buses 3,4. The loads (L1 to L7) are connected at buses 1,2,3,7,9,10,11 with an EV load at bus no. 7.The DGs are of10kVA 3 ϕ, 220 V (L-L), 60 Hrz with base value of power, voltage and frequency as Sbase=1000, Vbase=220 and ωbase=2π60 respectively with droop coefficients (mp and mq) as 0.001583 and 0.005905, respectively. Using parameters Vc=3m/s, VR=12m/s, and VF=20m/s and speed power curve of the turbine, the generated output power of wind farm has been calculated (Section 4.1) with the frequency counting technique [29]. The solar irradiance at the buses where PV modules are connected is having a Beta distribution function with α=12.62 and β=2.21 [29] at 43°C of nominal operating temperature. Weibull distribution has been considered for EV recharge power with parameters (shape and scale, i.e., Φ,φ) as 0.192329 and 1.87103 [40]. The parameters of EV, WTGS and PV are given in [29] and various other technical parameters of the islanded microgrid are given in Appendix B. In the case of an islanded microgrid, the standard deviation of 10% of the mean value of the load has been taken, with no slack bus as in typical load flow. All the simulation studies have been carried out on Intel Core i7-7500U 2.90 GHz using MATLAB environment.


### 5.2. Simulation and Methodology

The proposed methodology involves using statistical models as distribution functions for various uncertain loads and generations, calculating points and weights using GQPLF to calculate the moments of output variables of interest (voltage and power flow), and finally obtaining the distribution functions of output variables using appropriate expansion series. The methodology and procedure involved in the simulation are described (Figure3).(1)
Statistical models of various parameters, viz., WTGS, PV, EV, and loads, have been used to generate their respective distributions. As explained in Sections4.1–4.3, Weibull distribution has been considered for WTGS and EV loads, and beta distribution has been used for handling PV uncertainties. The distributions of WTGS, PV, and EV have been stored in various subroutines. As PV and WTGS generate active power, which is fed to the microgrid, active power is subtracted from the load at the buses where PV and WTGS are connected. Since these sources take reactive power for this operation, reactive power is added to the reactive load at these buses. However, the active power consumed by the EV has been added to the load, as it consumes power.(2)
From the distribution functions of WTGS, PV, and EV (obtained in step 1), weights and points have been calculated (as explained in equations (9) and (10) of Section 3.1 and Figure 2). The weights represent the weightage of points calculated, and no. of points can be varied in GQPLF.(3)
After obtaining weights and points for each distribution, the deterministic load flow of the microgrid has been calculated for every point to obtain points and weights of output variables (as explained in Sections2, 3.1 and Figures 1 and 2). Also, reactive power is required for the operation of WTGS, so this power is updated for a new value of updated voltage in each iteration. However, the active power of WTGS remains the same in every iteration.(4)
The output moments are calculated from the obtained output points and weights (as explained in equation (12) of Section 3.2 and Figure 2).(5)
From the output moments (i.e., moments of voltage and power flow), the distribution functions have been obtained using the appropriate expansion series, as explained in Section3.2 and Figure 2.For checking the efficacy of the methodology proposed above for different scenarios, the following cases have been considered:(i)
GQPLF with normal loads(ii)
GQPLF with different loadings(iii)
GQPLF with variation in system parametersFigure 3
Adopted methodology and simulation details.


### 5.3. GQPLF with Normal Loads

The Gauss quadrature-based load flow was performed using the approach described in Section3, and the resulting voltage distribution function (CDF and PDF) at bus no. 9 is depicted in Figure 4. Further, the distribution function of power (i.e., active and reactive) between the line connected to bus no. 10 and 11 is also shown in Figure 5. These distributions were obtained by following the methodology discussed in Section 5.2.Figure 4
Distributions of voltage at bus no. 9.Figure 5
Distribution function of power flows in line connected between bus no. 10-11.From these figures, it is observed that distributions of various output parameters can be obtained using the GQPLF method for a given input uncertain parameters, i.e., WTGS, PV, and EV.However, to check the effectiveness of the proposed method, microgrid load flow has also been carried out by the numerical-based Monte Carlo simulation (MCS) method, which is the most accurate for 10000 simulations [28, 47]. It entails those as follows:(1)
DG output distribution sampling, active and reactive power load distribution sampling, and EV loads sampling utilize the most appropriate sampling approach. In this case, the acceptance-rejection sampling method has been used.(2)
For each sample of input parameters, i.e., active and reactive power loads and generation, a microgrid deterministic load flow is calculated (as described in Section2)(3)
Calculation of output variables, i.e., voltage, angle, and active and reactive power loss for 10000 samples(4)
Calculate the frequency distribution and probability distribution from the output samples by converting the frequency values to probability values for each variable of interestThe resulting distributions of various parameters using MCS is also shown in Figures4 and 5 with statistical parameters in Table 1. The statistical parameters, in the form of moments, are the quantitative measure of the distribution function. The first raw moment μ1 is mean, the second central moment is variance, the third standardized moment is skewness, and the fourth is kurtosis of the distribution function.Table 1
Statistical moments of output variables.
Voltage moments at bus no. 9Methodμ1μ2μ3μ4GQPLF0.972040.945050.918720.89312MCS0.972010.945040.918710.89311Active power flow moments in the line between bus nos. 10-11Methodμ1μ2μ3μ4GQPLF0.318590.101520.032350.01031MCS0.318580.101500.032350.01030Reactive power flow moments in the line between bus nos. 10-11Methodμ1μ2μ3μ4GQPLF0.442860.196210.086970.03857MCS0.442840.196200.086970.03857From Figure4 and Table 1, it has been observed that the distribution function of voltage in the case of GQPLF lies between 0.925 p.u. and 1.015 p.u. with a mean of 0.97204 p.u. and for MCS distribution function, it lies between 0.928 p.u. and 1.013 p.u., with a mean of 0.97201 p.u. From this, it can be inferred that if the voltage limit has been defined as ±5% of normal voltage (1 p.u), then in this case, there is a probability of voltage violation on the lower side. The probable cause of this violation is the lack of reactive power in the system. Even though the probability of the violation is less, still reactive power compensation devices are required to keep the voltage within the defined limits. Similarly, if the active and reactive power flow limits are defined for Figure 5 having statistical parameters in Table 1, similar inferences can be made. Though the proposed method is versatile and can be used up to 15-point estimation, in this study, only three points have been considered, and it has been observed that with three points, reasonable efficacy has been obtained.From these figures and tables, it is observed that the Guass quadrature-based PLF method can be successfully applied to the probabilistic analysis of islanded microgrid with inclusion of EV loads.


### 5.4. GQPLF with Different Loadings

Furthermore, to check the effectiveness of the proposed method, analysis has also been conducted for different loadings. In this case, all the loads of microgrid (active and reactive) are incremented by5% and decremented by 25%, i.e., three cases have been considered—case 1: normal loadings, case 2: 0.75 the normal loadings, and case 3: 1.05 normal loadings. The resulting distributions of voltages for all three cases at bus no. 13 are shown in Figure 6.Figure 6
CDF and PDF of voltage at bus no. 13 for different values of loadings.From this figure, it is observed that with the increase in loadings, there is a decrease in voltage at bus no. 13; that is, when the loadings are increased to 1.05 p.u., the voltage CDF is well below the defined lower limit of 0.95 p.u. The cause of this is the increase in the reactive power requirements due to the increment of reactive power load, which reduces voltage. Furthermore, if the loading is reduced by25% of the normal load, the voltage CDF and PDF shift towards the right and is well within the defined upper and lower limits of voltage (±5% of normal voltage, i.e., 1 p.u.) due to reduction in reactive power requirements.Similar observations have been made for the power loss in the line connecting bus no. 2 and 3, as shown in Figure7. In this case, normal loading, 0.8 of normal loading, and 0.95 of normal loading have been considered. From this figure, it has been observed that with the decrease in loadings, the power loss distribution shifts towards the left, that is, there is a decrease in power loss due to decrease in load in the system.Figure 7
Distribution function of power loss in the line connecting the bus nos. 2 and 3.


### 5.5. GQPLF with Variation in System Parameters

Furthermore, analysis has also been carried out with the variation in line resistance from 0.9 to 1.05 the normal resistance. In this case, three cases have been considered, case 1: normal resistance, case 2: 0.9 of normal resistance, and case 3: 1.5 times normal resistance. For these scenarios, the resulting CDF of voltage at bus no. 12 and power loss in the line connecting between bus nos. 4 and 2 is shown in Figures8 and 9, respectively.Figure 8
Distribution function of voltage at bus no. 12 for different values of line resistances.Figure 9
Distribution function of power loss in the line connecting the bus nos. 4 and 2.From Figure8, it has been observed that with the decrease in line resistance, the bus voltage increases with the shifting of CDF of voltage towards right and vice-versa. In this case, the CDF which lies between 0.963 p.u. and 0.979 p.u. for a normal case has been shifted to 0.968 p.u. and 0.979 p.u. when the resistance is decreased to 0.9 of normal resistance, and there is shifting of CDF towards the left when line resistance is increased to 1.5 times normal resistance. The decrease in voltage due to the increase in resistance is due to the increase in losses.Observations similar to this have been made in power loss across the line connected between the bus nos. 4 and 2, where with the decrease in line resistance, there is a shift of power loss CDF towards the left and vice-versa, as shown in Figure9. This is due to decrease in power loss due to the increase in line resistance and vice-versa.The time taken by the GQPLF for the simulation of the 14-bus microgrid system is 60.9 seconds, which is approximately equal to one minute and is significantly less in comparison to the benchmark MCS method in which time taken for the same case with 10,000 simulations is 9318.7 seconds (155 minutes, approximately). Based on the findings of the preceding analysis, it has been observed that the proposed GQPLF can be successfully, effectively, and efficiently used for the probabilistic analysis of an islanded microgrid with WTGS and PV containing EV loads.The practical application of this method is in the planning and operational analysis of microgrids in the face of large-scale uncertainties due to various renewable energy resources like wind, solar, and electric vehicles. As the proposed method will provide the output variables of the microgrid in the form of probability distributions, i.e., values and likelihood, the severity of the system variable can be accessed by this method. It will benefit the power system planner and operator and provide an alternative perspective on system performance in the face of uncertainty.


## 5.1. Description of the System

The test system of an islanded microgrid under consideration consists of 14 bus systems with two buses having wind turbines and two buses containing photovoltaic systems, as shown in Appendix A. EV load is connected across one bus with six buses containing droop-controlled DGs.In this system, droop controlled DGs are connected at buses5,6,8,12,13,14 with wind generation systems at buses 1,2 and photovoltaic generation systems at buses 3,4. The loads (L1 to L7) are connected at buses 1,2,3,7,9,10,11 with an EV load at bus no. 7.The DGs are of10kVA 3 ϕ, 220 V (L-L), 60 Hrz with base value of power, voltage and frequency as Sbase=1000, Vbase=220 and ωbase=2π60 respectively with droop coefficients (mp and mq) as 0.001583 and 0.005905, respectively. Using parameters Vc=3m/s, VR=12m/s, and VF=20m/s and speed power curve of the turbine, the generated output power of wind farm has been calculated (Section 4.1) with the frequency counting technique [29]. The solar irradiance at the buses where PV modules are connected is having a Beta distribution function with α=12.62 and β=2.21 [29] at 43°C of nominal operating temperature. Weibull distribution has been considered for EV recharge power with parameters (shape and scale, i.e., Φ,φ) as 0.192329 and 1.87103 [40]. The parameters of EV, WTGS and PV are given in [29] and various other technical parameters of the islanded microgrid are given in Appendix B. In the case of an islanded microgrid, the standard deviation of 10% of the mean value of the load has been taken, with no slack bus as in typical load flow. All the simulation studies have been carried out on Intel Core i7-7500U 2.90 GHz using MATLAB environment.


## 5.2. Simulation and Methodology

The proposed methodology involves using statistical models as distribution functions for various uncertain loads and generations, calculating points and weights using GQPLF to calculate the moments of output variables of interest (voltage and power flow), and finally obtaining the distribution functions of output variables using appropriate expansion series. The methodology and procedure involved in the simulation are described (Figure3).(1)
Statistical models of various parameters, viz., WTGS, PV, EV, and loads, have been used to generate their respective distributions. As explained in Sections4.1–4.3, Weibull distribution has been considered for WTGS and EV loads, and beta distribution has been used for handling PV uncertainties. The distributions of WTGS, PV, and EV have been stored in various subroutines. As PV and WTGS generate active power, which is fed to the microgrid, active power is subtracted from the load at the buses where PV and WTGS are connected. Since these sources take reactive power for this operation, reactive power is added to the reactive load at these buses. However, the active power consumed by the EV has been added to the load, as it consumes power.(2)
From the distribution functions of WTGS, PV, and EV (obtained in step 1), weights and points have been calculated (as explained in equations (9) and (10) of Section 3.1 and Figure 2). The weights represent the weightage of points calculated, and no. of points can be varied in GQPLF.(3)
After obtaining weights and points for each distribution, the deterministic load flow of the microgrid has been calculated for every point to obtain points and weights of output variables (as explained in Sections2, 3.1 and Figures 1 and 2). Also, reactive power is required for the operation of WTGS, so this power is updated for a new value of updated voltage in each iteration. However, the active power of WTGS remains the same in every iteration.(4)
The output moments are calculated from the obtained output points and weights (as explained in equation (12) of Section 3.2 and Figure 2).(5)
From the output moments (i.e., moments of voltage and power flow), the distribution functions have been obtained using the appropriate expansion series, as explained in Section3.2 and Figure 2.For checking the efficacy of the methodology proposed above for different scenarios, the following cases have been considered:(i)
GQPLF with normal loads(ii)
GQPLF with different loadings(iii)
GQPLF with variation in system parametersFigure 3
Adopted methodology and simulation details.


## 5.3. GQPLF with Normal Loads

The Gauss quadrature-based load flow was performed using the approach described in Section3, and the resulting voltage distribution function (CDF and PDF) at bus no. 9 is depicted in Figure 4. Further, the distribution function of power (i.e., active and reactive) between the line connected to bus no. 10 and 11 is also shown in Figure 5. These distributions were obtained by following the methodology discussed in Section 5.2.Figure 4
Distributions of voltage at bus no. 9.Figure 5
Distribution function of power flows in line connected between bus no. 10-11.From these figures, it is observed that distributions of various output parameters can be obtained using the GQPLF method for a given input uncertain parameters, i.e., WTGS, PV, and EV.However, to check the effectiveness of the proposed method, microgrid load flow has also been carried out by the numerical-based Monte Carlo simulation (MCS) method, which is the most accurate for 10000 simulations [28, 47]. It entails those as follows:(1)
DG output distribution sampling, active and reactive power load distribution sampling, and EV loads sampling utilize the most appropriate sampling approach. In this case, the acceptance-rejection sampling method has been used.(2)
For each sample of input parameters, i.e., active and reactive power loads and generation, a microgrid deterministic load flow is calculated (as described in Section2)(3)
Calculation of output variables, i.e., voltage, angle, and active and reactive power loss for 10000 samples(4)
Calculate the frequency distribution and probability distribution from the output samples by converting the frequency values to probability values for each variable of interestThe resulting distributions of various parameters using MCS is also shown in Figures4 and 5 with statistical parameters in Table 1. The statistical parameters, in the form of moments, are the quantitative measure of the distribution function. The first raw moment μ1 is mean, the second central moment is variance, the third standardized moment is skewness, and the fourth is kurtosis of the distribution function.Table 1
Statistical moments of output variables.
Voltage moments at bus no. 9Methodμ1μ2μ3μ4GQPLF0.972040.945050.918720.89312MCS0.972010.945040.918710.89311Active power flow moments in the line between bus nos. 10-11Methodμ1μ2μ3μ4GQPLF0.318590.101520.032350.01031MCS0.318580.101500.032350.01030Reactive power flow moments in the line between bus nos. 10-11Methodμ1μ2μ3μ4GQPLF0.442860.196210.086970.03857MCS0.442840.196200.086970.03857From Figure4 and Table 1, it has been observed that the distribution function of voltage in the case of GQPLF lies between 0.925 p.u. and 1.015 p.u. with a mean of 0.97204 p.u. and for MCS distribution function, it lies between 0.928 p.u. and 1.013 p.u., with a mean of 0.97201 p.u. From this, it can be inferred that if the voltage limit has been defined as ±5% of normal voltage (1 p.u), then in this case, there is a probability of voltage violation on the lower side. The probable cause of this violation is the lack of reactive power in the system. Even though the probability of the violation is less, still reactive power compensation devices are required to keep the voltage within the defined limits. Similarly, if the active and reactive power flow limits are defined for Figure 5 having statistical parameters in Table 1, similar inferences can be made. Though the proposed method is versatile and can be used up to 15-point estimation, in this study, only three points have been considered, and it has been observed that with three points, reasonable efficacy has been obtained.From these figures and tables, it is observed that the Guass quadrature-based PLF method can be successfully applied to the probabilistic analysis of islanded microgrid with inclusion of EV loads.


## 5.4. GQPLF with Different Loadings

Furthermore, to check the effectiveness of the proposed method, analysis has also been conducted for different loadings. In this case, all the loads of microgrid (active and reactive) are incremented by5% and decremented by 25%, i.e., three cases have been considered—case 1: normal loadings, case 2: 0.75 the normal loadings, and case 3: 1.05 normal loadings. The resulting distributions of voltages for all three cases at bus no. 13 are shown in Figure 6.Figure 6
CDF and PDF of voltage at bus no. 13 for different values of loadings.From this figure, it is observed that with the increase in loadings, there is a decrease in voltage at bus no. 13; that is, when the loadings are increased to 1.05 p.u., the voltage CDF is well below the defined lower limit of 0.95 p.u. The cause of this is the increase in the reactive power requirements due to the increment of reactive power load, which reduces voltage. Furthermore, if the loading is reduced by25% of the normal load, the voltage CDF and PDF shift towards the right and is well within the defined upper and lower limits of voltage (±5% of normal voltage, i.e., 1 p.u.) due to reduction in reactive power requirements.Similar observations have been made for the power loss in the line connecting bus no. 2 and 3, as shown in Figure7. In this case, normal loading, 0.8 of normal loading, and 0.95 of normal loading have been considered. From this figure, it has been observed that with the decrease in loadings, the power loss distribution shifts towards the left, that is, there is a decrease in power loss due to decrease in load in the system.Figure 7
Distribution function of power loss in the line connecting the bus nos. 2 and 3.


## 5.5. GQPLF with Variation in System Parameters

Furthermore, analysis has also been carried out with the variation in line resistance from 0.9 to 1.05 the normal resistance. In this case, three cases have been considered, case 1: normal resistance, case 2: 0.9 of normal resistance, and case 3: 1.5 times normal resistance. For these scenarios, the resulting CDF of voltage at bus no. 12 and power loss in the line connecting between bus nos. 4 and 2 is shown in Figures8 and 9, respectively.Figure 8
Distribution function of voltage at bus no. 12 for different values of line resistances.Figure 9
Distribution function of power loss in the line connecting the bus nos. 4 and 2.From Figure8, it has been observed that with the decrease in line resistance, the bus voltage increases with the shifting of CDF of voltage towards right and vice-versa. In this case, the CDF which lies between 0.963 p.u. and 0.979 p.u. for a normal case has been shifted to 0.968 p.u. and 0.979 p.u. when the resistance is decreased to 0.9 of normal resistance, and there is shifting of CDF towards the left when line resistance is increased to 1.5 times normal resistance. The decrease in voltage due to the increase in resistance is due to the increase in losses.Observations similar to this have been made in power loss across the line connected between the bus nos. 4 and 2, where with the decrease in line resistance, there is a shift of power loss CDF towards the left and vice-versa, as shown in Figure9. This is due to decrease in power loss due to the increase in line resistance and vice-versa.The time taken by the GQPLF for the simulation of the 14-bus microgrid system is 60.9 seconds, which is approximately equal to one minute and is significantly less in comparison to the benchmark MCS method in which time taken for the same case with 10,000 simulations is 9318.7 seconds (155 minutes, approximately). Based on the findings of the preceding analysis, it has been observed that the proposed GQPLF can be successfully, effectively, and efficiently used for the probabilistic analysis of an islanded microgrid with WTGS and PV containing EV loads.The practical application of this method is in the planning and operational analysis of microgrids in the face of large-scale uncertainties due to various renewable energy resources like wind, solar, and electric vehicles. As the proposed method will provide the output variables of the microgrid in the form of probability distributions, i.e., values and likelihood, the severity of the system variable can be accessed by this method. It will benefit the power system planner and operator and provide an alternative perspective on system performance in the face of uncertainty.


## 6. Conclusion

This work proposes a Gauss quadrature-based probabilistic power flow method for an islanded microgrid with wind, solar, and load uncertainties, including electric vehicles. The proposed method has been verified on a 14-bus test microgrid system and compared with the numerical Monte Carlo simulation method to test its efficacy. The findings of simulation studies are as follows:(1)
The Gauss quadrature method can be effectively applied to an islanded microgrid under different scenarios(2)
With variable demand and generation, including electric vehicles, the distribution function of microgrid parameters such as voltage and power flow is accessible(3)
The microgrid operator is provided with an objective perspective on the operation of the system as a result of the distribution of various parameters, which reveals the possibility of voltage and power flow violations(4)
The proposed method can be extended by including multimodal load and generation distributions with the modelling of interval-based renewable energy sources(5)
The proposed method can be extended by including the correlation among the loads, generations, renewable energy resources, and a combination of these



---

*Source: 9569224-2022-11-25.xml*