Cumulant-based correlated probabilistic load flow considering photovoltaic generation and electric vehicle charging demand

Abstract

This paper applies a cumulant-based analytical method for probabilistic load flow (PLF) assessment in transmission and distribution systems. The uncertainties pertaining to photovoltaic generations and aggregate bus load powers are probabilistically modeled in the case of transmission systems. In the case of distribution systems, the uncertainties pertaining to plug-in hybrid electric vehicle and battery electric vehicle charging demands in residential community as well as charging stations are probabilistically modeled. The probability distributions of the result variables (bus voltages and branch power flows) pertaining to these inputs are accurately established. The multiple input correlation cases are incorporated. Simultaneously, the performance of the proposed method is demonstrated on a modified Ward-Hale 6-bus system and an IEEE 14-bus transmission system as well as on a modified IEEE 69-bus radial and an IEEE 33-bus mesh distribution system. The results of the proposed method are compared with that of Monte-Carlo simulation.

Keywords

battery electric vehicle / extended cumulant method / photovoltaic generation / plug-in hybrid electric vehicle / probabilistic load flow

Introduction

In recent years, increasing penetration of photovoltaic (PV) generations and replacement of internal combustion engine (ICE) by electric vehicles (EVs) have resulted in uncertainty factors in power systems. Probabilistic load flow (PLF) evaluates the steady state system performance of power system in a broader sense by accommodating aforementioned uncertainties []. In PLF studies, the uncertainties associated with input random variables (RVs) are probabilistically modeled. The fundamental task of PLF is to characterize probability distributions of the result variables (bus voltages and branch power flows) from given statistical information of input RVs and their multiple correlation cases. A critical review on PLF research is made in Ref. []. The various methods for PLF assessment include numerical, approximate and analytical methods. Monte-Carlo simulation (MCS) is a numerical method for PLF. It is considered as an accurate PLF method but is computationally burdensome. Because of its accuracy, MCS acts as reference for other PLF methods. The efforts in PLF research are to propose new methods which can accurately estimate the probability distributions of result variables in less computational effort. The cumulant-based analytical method [–] is computationally efficacious but inclusion of dependency among input RVs is a challenging factor.

Nowadays, major interest in PLF has been devoted to accommodating and analyzing the impact of PV generation uncertainty on power system [,–] and effect of uncertainty caused by EV charging on distribution system performance [–]. The authors in Ref. [] have analyzed the impact of high penetration of PV generation on transmission system performance at a certain instant of time using the joint cumulant-based formulation. PLF is also applied to distribution systems considering PV generations [–]. On the other hand, PLF research considering EV charging is also greatly appreciated. In Ref. [], all EV charging processes have been modeled by using a single queuing model. The authors in Ref. [] have proposed an EV charging model by considering the daily arrival time, departure time, and travelling distance of EVs. The EV charging model using real time data collected from a transportation authority has been detailed in Ref. []. The expected number of EVs in the future is predicted from the data. In order to avoid charging during peak load time, a model based on time of use has been developed in Ref. []. In Refs. [] and [], the EV charging process has been modeled considering residential community charging. In Ref. [], the EV charging process is modeled by taking different types of EV classes in both residential communities and charging stations.

From the literature review, it is noticed that none of the PLF research considering PV generation and EV charging demand incorporates multiple correlation cases. Therefore, this paper uses EV demand models at the residential community charging as well as at the charging stations. It adopts two separate models for the plug-in hybrid electric vehicle (PHEV) and the battery electric vehicle (BEV) respectively. Different queuing models are used for residential community charging and station charging. Further, multiple input correlation cases are incorporated. An analytical method such as, extended cumulant method (ECM) is employed for PLF. Since ECM is an analytical method, a sensitivity matrix-based PLF model is developed for both transmission and distribution systems. The results of ECM are compared with those of MCS for validation purpose.

Probabilistic load flow model

Transmission system

In a power system, comprising the ‘n’ number of buses (‘m’P\left|V\right| buses, m≤n) and the ‘1’ number of branches, the linearized PLF model [] is given as

(1)x_i={\displaystyle \sum _{j=\mathrm{1}}^{2n-m-2}{a}_{ij}(y_j^{\mathrm{D}}+y_j^{\mathrm{C}})+x_{i\mathrm{0}}\hspace{0.17em},}

(2)z_i={\displaystyle \sum _{j=\mathrm{1}}^{2n-m-2}{b}_{ij}(y_j^{\mathrm{D}}+y_j^{\mathrm{C}})+z_{i\mathrm{0}}\hspace{0.17em},}

where x_{i\mathrm{0}}=x_i^{\mathrm{0}}-{\displaystyle \sum _{j=\mathrm{1}}^{2n-m-2}{a}_{ij}{y}_j^{\mathrm{0}}}, z_{i\mathrm{0}}=z_i^{\mathrm{0}}-{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{2n}-m-2}{b}_{ij}{y}_j^{\mathrm{0}}}.

In Eqs. (1) and (2), x , y and z are the vectors of bus voltages, bus power injections and line power flows respectively; x⁰, y⁰ and z⁰ are expected values of x, y and z respectively; a_ij and b_ij are sensitivity coefficients; the superscripts “D” and “C” stand for “discrete” and “continuous” respectively, and scalarsx_i0 and z_i0 can be obtained from the deterministic load flow [].

Radial distribution system

Sensitivity matrix for bus voltage magnitudes

Consider a radial distribution system consisting of ‘n’ number of buses, the feeding bus (root bus) being numbered as ‘0’ and the remaining as 1,2,...,n-1, the magnitude of voltage at any kth bus is expressed as Ref. [].

(3)\left|V_k\right|=\left|V_{\mathrm{0}}\right|-\Delta \left|V_k\right|=\left|V_{\mathrm{0}}\right|-{\displaystyle \frac{\mathrm{1}}{\left|V_{\mathrm{0}}\right|}}{\displaystyle \sum _{i=\mathrm{1}}^{n\mathrm{-1}}(P_i{R}_{ik}+Q_i{X}_{ik})},

where \left|V_{\mathrm{0}}\right| denotes the magnitude of voltage at the feeding bus; \Delta \left|V_k\right| is the total voltage drop from feeding point to the kth bus; P_i and Q_i are the real and reactive bus load powers at the ith bus; R_ik and X_ik represents respectively the total resistance and reactance of the common element on the supplying paths of theith and kth bus starting from the root bus.

The compact form of Eq. (3) can be expressed as

(4){[\mathit{V}]}_{n\times 1}={[{\mathit{K}}_{\mathrm{1}}]}_{n\times 2n}\hspace{0.17em}{\left[\begin{array}{l}\mathit{P}\\ \mathit{Q}\end{array}\right]}_{2n\times 1},

where K₁ is the bus voltage magnitude sensitivity matrix, and load powers P and Q are RVs.

Sensitivity matrix for branch power flows

Let the radial system have l number of branches. The development of branch power flow sensitivity matrix is a two-step process. First, the system is assumed as lossless, and lossless branch power flow sensitivity matrixK₂ as a function of bus load powers is calculated. Then, branch loss is expressed as a function of bus load powers, and branch power loss sensitivity matrixK₃ is calculated. Finally, the sum of K₂ and K₃ yields the branch power flow sensitivity matrix K₄.

Now, considering the lossless system, the real and reactive branch power flows in any branch i-j (branch connecting bus i and bus j ) can be obtained as

(5){P^{\prime }}_{l,\hspace{0.17em}i-j}={\displaystyle \sum _{k\in {\Lambda }_j}{P}_k},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{Q^{\prime }}_{l,\hspace{0.17em}i-j}={\displaystyle \sum _{k\in {\Lambda }_j}{Q}_k\hspace{0.17em}},

where {\Lambda }_j denotes the set of all buses supplied via the jth bus, including itself.

In compact form, Eq. (5) can be expressed as

(6){\left[\begin{array}{l}{{\mathit{P}}^{\prime }}_l\\ {Q^{\prime }}_l\end{array}\right]}_{2l\times 1}={[{\mathit{K}}_{\mathrm{2}}]}_{2l\times 2n}\hspace{0.17em}{\left[\begin{array}{l}\mathit{P}\\ \mathit{Q}\end{array}\right]}_{2n\times 1}\hspace{0.17em}.

The real and reactive power loss in branch i-j can be expressed as

(7)\begin{array}{l}{P}_{\mathrm{Loss},i-j}=R_{i-j}(I_{p\hspace{0.17em}(i-j)}^{\mathrm{2}}+I_{q\hspace{0.17em}(i-j)}^{\mathrm{2}})\hspace{0.17em},\\ Q_{\mathrm{Loss},i-j}=X_{i-j}(I_{p\hspace{0.17em}(i-j)}^{\mathrm{2}}+I_{q\hspace{0.17em}(i-j)}^{\mathrm{2}})\hspace{0.17em},\end{array}

where R_i-j and X_i-j are the resistance and reactance of the branch i-j, and I_p(i-j) and I_q(i-j) are the real and imaginary parts of the current in branch i-j i.e., I_(i-j).

The linearization of the non-linear terms I_{p(i-j)}^{\mathrm{2}} and I_{q(i-j)}^{\mathrm{2}} in Eq. (7) are expressed as in Ref. [],

(8)\begin{array}{l}{I}_{p(i-j)}^{\mathrm{2}}=\mathrm{2}{I}_{p\mathrm{0}(i-j)}{I}_{p(i-j)}-I_{p\mathrm{0}(i-j)}^{\mathrm{2}}\hspace{0.17em},\\ I_{q(i-j)}^{\mathrm{2}}=\mathrm{2}{I}_{q\mathrm{0}(i-j)}{I}_{q(i-j)}-I_{q\mathrm{0}(i-j)}^{\mathrm{2}},\end{array}

where I_p0(i-j) and I_q0(i-j) are the expected values of I_p(i-j) and I_q(i-j) respectively.

Referring to Fig. 1, I_i-j can be expressed as

(9)I_{i-j}=I_{j-j}+I_{\mathrm{br}(j)} ,

where I_j-j can be calculated using

(10)I_{j-j}={\displaystyle \frac{P_j-Q_j}{{(V_j)}^{\ast }}}.

Using Eqs. (8)–(10), the compact form of Eq. (7) can be expressed as

(11){\left[\begin{array}{l}{\mathit{P}}_{\mathrm{Loss}}\\ {\mathit{Q}}_{\mathrm{Loss}}\end{array}\right]}_{2l\times 1}={[{\mathit{K}}_{\mathrm{3}}]}_{2l\times 2n}{\left[\begin{array}{l}\mathit{P}\\ \mathit{Q}\end{array}\right]}_{2n\times 1}\hspace{0.17em}.

Finally, K₄ is obtained as K₄=K₂+K₃.

Mesh distribution system

Unlike the radial distribution system, the process of obtaining the sensitivity matrices K₁ and K₄ in case of mesh distribution system is not straight forward. First, the mesh distribution system is converted into radial distribution systems by breaking loops at one of their local buses. This results in loop breaking points (LBPs) as indicated in Fig. 2. The current injections at the LBPsj₁ and j₂ are I_tj and –I_tj respectively where, I_tj is the current flowing through the jth bus before loop breaking. Now the matrices K₁ and K₄ can be obtained using the steps as followed for the radial distribution system.

Calculation of LBP current injections

For the calculation of LBP current injections, the Thevenin equivalent matrix [Z_Th] is developed. The following iterative process is used in order to obtain an elementI_tj of the current injection vector [I_t].

1)Develop [Z_Th].

2)Apply the forward-backward sweep iteration method. Calculate the voltage drops across the ends of LBPs. Note: Initial value of currents for the first iteration is considered as zero.

3)Incremental changes in current injections are calculated using, [DI_t]=[Z_Th]^–1[V_Th], where [V_Th] is the vector including voltage drops between both ends of LBPs.

4)At any kth iteration, update the LBP current injections using, [I_t]^k=[I_t]^k–1+[DI_t]^k.

5)Repeat steps 2) to 4) until convergence is reached. The convergence criterion is the minimum of the voltage drop between LBPs, which should be less than 0.01 V.

Power injections at LBPs

The LBP power injections can be obtained by referring to Fig. 3.

Let S_j be the complex power demand at the jth bus. After creation of LBPs j₁ and j₂, the complex power demand at j₁ is given as S_j1=S_j+(I_tjV_j1)* and that at j₂ is given as S_j2= (–I_tjV_j2)*, V_j1 and V_j2 are the voltages at j₁ and j₂ respectively.

Uncertainty modeling

Photovoltaic generation

In PV systems, the real power generation P_PV strongly depends on solar irradiance which is unpredictable due to climatic conditions (cloud and fog).

According to Ref. [], P_PV is given as

(12)P_{\mathrm{PV}}={\eta }_{g}rA(\mathrm{1}-K_T\Delta T_F)\hspace{0.17em},

where r is the irradiance, DT_F is the forecasting error of PV cell temperature, and K_T is the temperature coefficient.

Introducing random variables R and T as functions of r and DT respectively, Eq. (12) can be rewritten as

(13)P_{\mathrm{PV}}=RT,

where R=h_grA, T=1–K_TDT.

The linearization of Eq. (13) is given as

(14)P_{\mathrm{PV}}={\mu }_R{T}^{\prime }+R\hspace{0.17em},

where m_R is the mean value of R, \hspace{0.17em}{T}^{\prime }=T-\mathrm{1}.

The linearization of Eq. (14) can be inferred from Ref. [].

For a solar park consisting of ‘\mathit{nr}’ number of PV units, aggregate generation is given by

(15)P_{\mathrm{PVA}}=P_{{\mathrm{PV}}_1}+P_{{\mathrm{PV}}_2}+\mathrm{...}+P_{{\mathrm{PV}}_{\mathit{nr}}}.

Electric vehicle power demand

There are two types of EVs, the first type is PHEV powered by battery as well as ICE while the second type is BEV powered only by battery. EVs need electric power for charging process. Since each EV has different battery makes and capacities, their charging process is different and varies from time to time. Probabilistic modeling of a single PHEV charging demand follows the formulations developed in Ref. []. The operating status is given as

(16)P_{\mathrm{EV}}={\displaystyle \frac{E_{\mathrm{BAT}}}{E_{\mathrm{BAT}}+E_{\mathrm{ENG}}}},

where E_BAT is the energy supplied to the vehicle by battery and E_ENG is the energy supplied to the vehicle by ICE.

The operating status P_EV and the battery capacity B_CAP are correlated with each other in the case of PHEV and are modeled using Eq. (17). However, they are independent in the case of BEVs sinceP_EV using Eq. (16) is always unity as E_ENG= 0.

(17)\left(\begin{array}{l}{P}_{\mathrm{EV}}\\ B_{\mathrm{CAP}}\end{array}\right)=\left(\begin{array}{l}{\mu }_{P_{\mathrm{EV}}}\\ {\mu }_{B_{\mathrm{CAP}}}\end{array}\right)+\mathit{L}\hspace{0.17em}\left(\begin{array}{l}{Z}_{\mathrm{1}}\\ Z_{\mathrm{2}}\end{array}\right)\hspace{0.17em}.

In Eq. (17), Z₁ and Z₂ are independent standard Gaussian RVs, and L is the Cholesky decomposition of

\left[\begin{array}{ll}{\sigma }_{P_{\mathrm{EV}}}^{\mathrm{2}}& \mathrm{0.8}{\sigma }_{P_{\mathrm{EV}}}{\sigma }_{B_{\mathrm{CAP}}}\\ \mathrm{0.8}{\sigma }_{P_{\mathrm{EV}}}{\sigma }_{B_{\mathrm{CAP}}}& {\sigma }_{B_{\mathrm{CAP}}}^{\mathrm{2}}\end{array}\right]\hspace{0.17em},

where

\begin{array}{l}{\mu }_{P_{\mathrm{EV}}}={\displaystyle \frac{P_{\mathrm{EV,min}}+P_{\mathrm{EV,max}}}{\mathrm{2}}}\hspace{0.17em},\hspace{0.17em}{\mu }_{B_{\mathrm{CAP}}}={\displaystyle \frac{B_{\mathrm{CAP,min}}+B_{\mathrm{CAP,max}}}{\mathrm{2}}}\hspace{0.17em},\\ {\sigma }_{P_{\mathrm{EV}}}={\displaystyle \frac{P_{\mathrm{EV,max}}-P_{\mathrm{EV,min}}}{\mathrm{4}}}\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}{\sigma }_{B_{\mathrm{CAP}}}={\displaystyle \frac{B_{\mathrm{CAP,max}}-B_{\mathrm{CAP,min}}}{\mathrm{4}}}.\end{array}

The energy consumption per mile driven by the EV is given as

(18)M_{\mathrm{E}}=E_{\mathrm{A}}{(P_{\mathrm{EV}})}^{E_{\mathrm{B}}},

where E_A and E_B are constant coefficients dependent on EV class.

The daily driven miles M_D follows the lognormal distribution and is given as

(19)M_{\mathrm{D}}=e^{{\mu }_{\mathrm{m}}+{\sigma }_{\mathrm{m}}Z},

where Z is standard Gaussian RV, distribution parameters m_m and s_m are the mean and standard deviation values of natural logarithm of M_D, and m_m and s_m are related to the non-logarithmic mean (m_md) and standard deviation (s_md) as

(20){\mu }_{\mathrm{m}}=\ln \hspace{0.17em}\left({\displaystyle \frac{{\mu }_{\mathrm{md}}^2}{\sqrt{{\mu }_{\mathrm{md}}^2+{\sigma }_{\mathrm{md}}^2}}}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\sigma }_{\mathrm{m}}=\sqrt{\ln \hspace{0.17em}\left(\mathrm{1}+{\displaystyle \frac{{\sigma }_{\mathrm{md}}^2}{{\mu }_{\mathrm{md}}^2}}\right)}.

In this paper, the values of m_md and s_md are assumed as 40 miles and 20 miles respectively.

The daily recharge energy of each EV is given as

(21)

where E_M is the maximum driving distance and is given as

(22)E_{\mathrm{M}}={\displaystyle \frac{B_{\mathrm{CAP}}}{M_{\mathrm{E}}}}.

The aggregate charging demand for BEV and PHEV is obtained by using the queuing theory. The queuing modelM/M/C is followed for charging station, in which the first M denotes the inter arrival time of customers following exponential queue with mean time T_l. The second M denotes the service time of customer with mean time T_m, and C denotes the maximum number of customers being served at a time. The number of customers waiting for the service is considered to be infinite. The probability of ‘n₁’ number of customers being charged at any instant is given as

(23)P_{n_1}=\{\begin{array}{l}\begin{array}{ll}{\left\{{\displaystyle \sum _{i=\mathrm{0}}^{C-\mathrm{1}}{\displaystyle \frac{{(\rho C)}^i}{i\hspace{0.17em}!}}+{\displaystyle \frac{{(\rho C)}^C}{C!}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\rho }}\right)}\right\}}^{-\mathrm{1}}& \mathrm{ }{n}_1=\mathrm{0},\end{array}\\ \begin{array}{ll}{\displaystyle \frac{{(\rho C)}^{n_1}}{n_1\hspace{0.17em}!}}{P}_{\mathrm{0}}& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{ }{n}_1=\mathrm{1},\hspace{0.17em}\mathrm{2},\hspace{0.17em}\hspace{0.17em}\mathrm{...},\hspace{0.17em}C,\end{array}\end{array}

where r is the occupation rate given as r= T_m/CT_l; P₀ is the probability of zero number of EV charging.

The value of r should be less than 1 so that the queue does not explode. The queuing model M/M/C/K/N_max (C≤K≤N_max) is followed in residential community charging. Here, K is the maximum number of customers being served or waiting in a queue and N_max is the maximum number of possible customers to be served in a queue. The number of customers to be served in residential community charging is limited as compared to that in a charging station. Hence, this is the difference between charging at residential community and charging station.

The probability of ‘n₁’ number of customers being charged simultaneously in residential community is given as

(24)P_{n_1}=\{\begin{array}{l}\begin{array}{ll}{\left({\displaystyle \sum _{i=\mathrm{0}}^C\left(\begin{array}{l}{N}_{\max }\\ i\end{array}\right)\hspace{0.17em}{(\rho C)}^i+{\displaystyle \sum _{i=C+\mathrm{1}}^K\frac{N_{\max }!}{(N_{\max }-i)\hspace{0.17em}!}}{\displaystyle \frac{{(\rho C)}^i}{C!\hspace{0.17em}{C}^{i-C}}}}\right)}^{-\mathrm{1}}& n_1=\mathrm{0},\end{array}\\ \begin{array}{ll}{(\rho C)}^{n_1}\left({\displaystyle \frac{N_{\max }}{n_1}}\right)P_{\mathrm{0}}& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{n}_1=\mathrm{1},\hspace{0.17em}\mathrm{2},\hspace{0.17em}\cdots ,\hspace{0.17em}C.\end{array}\end{array}

The time taken by each EV to charge follows exponential distribution i.e.,

(25)T=-T_{\mu }\ln U,

where U is the uniformly distributed RV.

PHEVs have to be charged for a minimum time T_min. Moreover, the charging time should not exceed the maximum time T_max due to its battery capacity. The maximum time is not included in BEV as it is allowed to be fully charged. Since charging stations require fast charging, level 3 charging (400V/63A) is used. However, in the case of residential communities, slow charging is acceptable. Therefore, level 1 charging (230V/16A) is considered. The charging voltageV and maximum current I_max are determined from the charging levels. The overall charging demand of EV can be calculated as

(26)P={\displaystyle \sum _{i=\mathrm{1}}^{n_{\mathrm{1}}}V{I}_i},

where I_i is the charging current of the ith EV.

Extended cumulant method

Estimation of cumulants of result variables

ECM is the extension of conventional cumulant method to include input correlations []. Starting with a two input case, a generalization of ‘nr’ correlated inputs is explained. For two correlated input RVs X₁ and X₂, the cumulants of Eq. (27) can be obtained using Eqs. (28)‒(30).

(27)Y=X_{\mathrm{1}}\pm X_{\mathrm{2}},

(28)C_{Y,\hspace{0.17em}k}=A_{\mathrm{1}}(k)\hspace{0.17em}{C}_{X_{\mathrm{1}},\hspace{0.17em}k}+A_{\mathrm{2}}(k)\hspace{0.17em}{C}_{X_{\mathrm{2}},\hspace{0.17em}k},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\sigma }_{X_{\mathrm{2}}}\ge {\sigma }_{X_{\mathrm{1}}},

(29)A_{\mathrm{1}}(k)={(\mathrm{1}+s)}^k-s^k,\hspace{0.17em}{A}_{\mathrm{2}}(k)={(\pm \mathrm{1})}^k,

(30)s=\pm {\rho }_{X_{\mathrm{1}}{X}_{\mathrm{2}}}\left({\displaystyle \frac{{\sigma }_{X_{\mathrm{2}}}}{{\sigma }_{X_{\mathrm{1}}}}}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\sigma }_{X_{\mathrm{2}}}\ge {\sigma }_{X_{\mathrm{1}}\hspace{0.17em}}.

As a generalization, for ‘nr’ correlated input RVs, the cumulants of Eq. (31) can be obtained using Eqs. (32)‒(35).

(31)Y=X_{\mathrm{1}}\pm X_{\mathrm{2}}\pm X_{\mathrm{3}}\pm \cdots \pm X_{nr}.

Introducing W_i=X_i±X_i+1 for i=1,2,...,nr–1, a compact representation of Eq. (31) is given as

(32)Y=W_{nr-\mathrm{2}}\pm X_{nr}\mathrm{.}

For the W_i given by Eq. (33), the standard deviation of W_i and the correlation coefficient between W_i and X_i+2 are given by Eqs. (34) and (35).

Note: W₀=X₁,

(33)W_i=W_{i-\mathrm{1}}\pm X_{i+\mathrm{1}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\le i\le nr-2,

(34){\sigma }_{W_i}=\sqrt{{\sigma }_{X_i}^{\mathrm{2}}\pm \mathrm{2}{\rho }_{X_i{X}_{i+\mathrm{1}}}{\sigma }_{X_i}{\sigma }_{X_{i+\mathrm{1}}}+{\sigma }_{X_{i+\mathrm{1}}}^{\mathrm{2}}},

(35){\rho }_{W_i{X}_{i+\mathrm{2}}}={\displaystyle \frac{{\rho }_{X_i{X}_{i+\mathrm{2}}}{\sigma }_{X_i}\pm {\rho }_{X_{i+\mathrm{1}}{X}_{i+\mathrm{2}}}{\sigma }_{X_{i+\mathrm{1}}}}{{\sigma }_{W_i}}}.

Approximation of probability distributions of result variables

The probability distributions of the result variables are approximated from their first five cumulants using the Cornish-Fisher’s series expansion method (CFM). CFM approximates cumulative probability plot of a result variableY using its q-quantile with the help of quantile obtained from the standard Gaussian distribution i.e. Q(q). The quantile of a result variable Y can be approximated by using Eqs. (36) and (37) which is explained in Ref. []. In this paper, permilles (1000-quantiles) are used to approximate the cumulative probability plots of result variables.

(36)\begin{array}{l}{Y}^{\prime }(q)=Q(q)+{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}[{\{Q(q)\}}^{\mathrm{2}}-\mathrm{1}]\hspace{0.17em}{C}_{\mathrm{3}}\\ +{\displaystyle \frac{\mathrm{1}}{\mathrm{24}}}[{\{Q(q)\}}^{\mathrm{3}}-\mathrm{3}\{Q(q)\}]\hspace{0.17em}{C}_{\mathrm{4}}\\ -{\displaystyle \frac{\mathrm{1}}{\mathrm{36}}}[\mathrm{2}{\{Q(q)\}}^{\mathrm{3}}-\mathrm{5}\{Q(q)\}]\hspace{0.17em}{C}_{\mathrm{3}}^{\mathrm{2}}\\ +{\displaystyle \frac{\mathrm{1}}{\mathrm{120}}}[{\{Q(q)\}}^{\mathrm{4}}-\mathrm{6}{\{Q(q)\}}^{\mathrm{2}}+\mathrm{3}]\hspace{0.17em}{C}_{\mathrm{5}}\\ -{\displaystyle \frac{\mathrm{1}}{\mathrm{24}}}[{\{Q(q)\}}^{\mathrm{4}}-\mathrm{5}{\{Q(q)\}}^{\mathrm{2}}+\mathrm{2}]\hspace{0.17em}{C}_{\mathrm{3}}{C}_{\mathrm{4}}\\ +{\displaystyle \frac{\mathrm{1}}{\mathrm{324}}}[\mathrm{12}{\{Q(q)\}}^{\mathrm{4}}-\mathrm{53}{\{Q(q)\}}^{\mathrm{2}}+\mathrm{17}]\hspace{0.17em}{C}_{\mathrm{3}}^{\mathrm{3}}\hspace{0.17em},\end{array}

(37)Y=(Y^{\prime }){\rho }_Y+{\mu }_Y.

Computational procedures

In order to accomplish PLF in the sensitivity matrix-based model, the computational procedures of MCS and ECM are discussed underneath.

Monte-Carlo simulation

Transmission system

1)The correlated samples pertaining to input RVs as per the given correlation coefficient matrix are incorporated using polynomial normal transformation technique [].

2)The samples of result variables are obtained by simulating each set of samples of input RVs using Newton-Raphson iterative method.

3)Statistical moments as well as probability distributions of result variables can be established from the samples as obtained in 2).

Distribution system

1)The number of PHEVs and BEVs being charged at an instant is generated randomly using Eq. (23) for the charging station and Eq. (24) for the residential community. Then for each case, follow 2) to 8).

2)Select the parameters of PHEV and BEV class as per their market shares.

3)For PHEV, generate values of P_EV and B_CAP using Eq. (17). For BEV, P_EV=1, hence randomly generate B_CAP as the Gaussian distribution with a mean value of m_CAP and a standard deviation of s_CAP. The values of B_CAP exceeding the limiting range as provided in Table 1 are set to the limiting values.

4)Calculate M_E using Eq. (18) and randomly generate M_D using Eq. (19).

5)Calculate E_D using Eq. (21) for PHEV and E_D= M_DM_E for BEV.

6)Randomly generate charging time T using Eq. (25).

7)For PHEV, I=\min \left[{\displaystyle \frac{E_{\mathrm{D}}}{VT}},\hspace{0.17em}{I}_{\max }\right], where I_max and V depend on the level of charging as described in Section 3. For BEV, I={\displaystyle \frac{E_{\mathrm{D}}}{VT}}.

8)Calculate the overall charging demand using Eq. (26).

9)The correlations between the samples of input RVs as per the given correlation coefficient matrix are incorporated using polynomial normal transformation technique [].

10) Finally, the samples of result variables are obtained by simulating each set of samples of input RVs using the forward-backward sweep iterative method.

11) In the case of the mesh system, the system is first converted into a radial distribution system as discussed in Section 2.2 and then 10) is applied.

12) Statistical moments as well as probability distributions of result variables can be established from the samples as obtained in 10) and 11).

Extended cumulant method

1)Obtain the first five cumulants of (a) EV load demands from the random generations obtained in Section 5.1.2 3) and (b) other input RVs such as PV generations and aggregate bus load power from their distribution parameters.

2)Develop bus voltage and branch power flow sensitivity matrices.

3)The first five cumulants of result variables are obtained using 2) and the formulation developed in Section 4.1.

4)CFM as discussed in Section 4.2 is applied to approximate the cumulative probability plots of result variables.

Case studies and result analysis

The applicability of ECM is verified in both transmission and distribution systems. The performance of the proposed method is compared with that of MCS. Solution accuracy and computational efficiency are the two performance criteria considered.

PLF in transmission systems

ECM is verified on two transmission systems such as a modified Ward-Hale 6-bus system [] and a modified IEEE 14-bus system []. The verification includes the examination of the effect of PV penetration on load flow results in the 6-bus system and the analysis of the effect of various input correlations on probability distributions of result variables in the 14-bus system. The modifications in the systems are done by installing solar parks. The technical details of solar parks are given in Table 1. The PV penetration is based on the local bus load demand. In the case of Gaussian load distributions, the mean values are taken as specified deterministic data and the standard deviation values taken as 7% of the mean value.

Ward-Hale 6-bus system

The Ward-Hale 6-bus system is a simple system consisting of two generating units, seven transmission lines and two transformer branches. The sets (PV₁, PV₂) and (PV₃, PV₄) represent PV generation of units pertaining to the solar parks 1 and 2 respectively. The PV generation is expressed by Eq. (14).

Case 1 Uncorrelated input RVs. In this case, the input RVs are uncorrelated and the effect of PV generation on result variables in the base case and in the increased penetration case is studied. The first five cumulants are obtained using ECM and the results are compared with that of MCS using different combinations of 15000 input samples for each RV. The absolute percentage errore_X is used as the error index. The value of e_X (average of e_µ and e_s) in obtaining mean values and standard deviations of result variables is computed by using Eqs. (38) and (39).

The absolute percentage error in mean value (m) is calculated as

(38)e_{\mu }=\left|{\displaystyle \frac{{\mu }_{X,\hspace{0.17em}\mathrm{ECM}}-{\mu }_{X,\hspace{0.17em}\mathrm{MCS}}}{{\mu }_{X,\hspace{0.17em}\mathrm{MCS}}}}\right|\times \mathrm{100.}

The absolute percentage error in standard deviation (s) is calculated as

(39)e_{\sigma }=\left|{\displaystyle \frac{{\sigma }_{X,\hspace{0.17em}\mathrm{ECM}}-{\sigma }_{X,\hspace{0.17em}\mathrm{MCS}}}{{\sigma }_{X,\hspace{0.17em}\mathrm{MCS}}}}\right|\times \mathrm{100.}

Due to the increase in PV penetration (by 25% of the base case value), the voltage profile of the system buses is found to improve and the system real power loss is reduced by 17.6% from the base case value. The probability density plots of system real power loss both in the base case and after increased penetration is shown in Fig. 4.

Case 2 Correlated input RVs. In this case, the effect of input correlation on PLF result is discussed. The correlation coefficients between input RVs are kept fixed. It is supposed that RVsR₁, R₂, R₃ and R₄ are correlated with RVs T₁, T₂, T₃ and T₄ respectively by a fixed correlation coefficient of 0.3. The remaining RVs are assumed to be correlated as per the correlation coefficient matrix expressed in Table 2.

The real and reactive component of load powers at buses 3, 5 and 6 are completely correlated because of the constant power factor assumption. Since correlation coefficient is independent of change of scale, the correlation coefficient matrix pertaining toQ_D3, Q_D5 and Q_D6 is exactly the same as the one which P_D3, P_D5 and P_D6 have formed in Table 2. The scatter plots of various input RVs in the base case of PV penetration are depicted in Fig. 5.

The effectiveness of ECM using the proposed model under various conditions is compared in Table 3.

IEEE 14-bus system

The IEEE 14-bus system consists of five generating units and twenty lines including three transformer branches.P_D5, P_D13 and P_D14 are assumed as RVs whereas others are presumed to be precisely obtained. It is supposed that RVsR₁, R₂, R₃ and R₄ are independent with that of T₁, T₂, T₃ and T₄. The correlation coefficient matrix for remaining RVs is expressed in Table 4. The first five cumulants of result variables are obtained using the proposed method. The development of the PMCC matrix for reactive load demand follows the same steps as that is considered in the 6-bus system.

In this section, the influence of input correlation on result variables using ECM is discussed. In order to perceive the influence load demand correlation on system real power loss, two extreme values of correlation coefficient such as 0 and 1 are considered. The probability density plots are compared in Fig. 6 where PMCC stands for product moment correlation coefficient. In this case, the correlation between PV generations is considered as per Table 4.

The impact of various PV generation correlations on bus voltage magnitude and angle of bus-13 is depicted in Fig. 7. The load demand correlation is kept fixed as per the base case correlation coefficient matrix of Table 4. The proposed method is capable of handling both uncorrelated and correlated inputs precisely. For both\left|V_{\mathrm{13}}\right| and d₁₃, ECM is closer to MCS.

PLF in distribution systems

The performance of the proposed method is verified on a modified IEEE 69-bus radial and a modified IEEE 33-bus mesh distribution system considering input correlations. The branch and bus data of the systems are adopted from Ref. []. The solution accuracy and computational efficiency are the performance criteria considered. The original loads connected at buses 66 and 68 in the case of radial and at buses 20 and 7 in the case of mesh distribution systems are replaced by PHEV demands at the charging station and the residential community respectively. Similarly, the original loads at buses 50 and 64 in radial and at buses 32 and 24 in mesh distribution systems are replaced by BEV demands at the charging station and the residential community respectively. Aggregate real and reactive bus power demands at buses 11, 16, 17 and 20 in the case of radial and at buses 11, 14, 16, and 17 in the case of mesh distribution systems are assumed to follow Gaussian distribution with mean values as specified deterministic data and the standard deviation values taken as 7% of the mean value.

The details of EV class parameters are provided in Table 5. The EV charging loads at the residential community and the charging station are correlated to each other. It is assumed that the number of PHEVs charged at the residential community is high when the number of PHEVs charged at the charging station is low and vice versa. Hence, a negative value of correlation coefficient is assumed between PHEV demands in residential charging and community charging (see Table 6). Similarly, a negative value of correlation coefficient is assumed between the number of BEVs charged at the residential community and at the charging station. The correlation coefficient between PHEVs as residential community load and BEVs as residential load is assumed to be 0.2 and the correlation between PHEV as residential community load and BEV as charging station load is assumed to be –0.2. The correlation coefficient between PHEV as residential community load and PHEV as charging station load is considered as –0.8. Both ECM and MCS are applied for PLF in both test systems. The cumulative probability plots of\left|V_{\mathrm{68}}\right|, P_L,14–15 and P_L,15–16 with and without considering base case correlation is compared in Figs. 8, 9 and 10 respectively in the case of the radial distribution system. In all the three cases, ECM plots are close to MCS plots. The consideration of the base case correlation has decreased the variance of\left|V_{\mathrm{68}}\right|but increased the variance of P_L,14–15 and P_L,15–16.

In the case of the mesh system, the correlation coefficient matrix is the same as that in Table 6 but is defined among the RVsP_D11, P_D14, P_D16, P_D17, Q_D11, Q_D14, Q_D16, Q_D17, P_D7, P_D20, P_D24 and P_D32 . Cumulative probability plots P_L,16–17 and \left|V_{\mathrm{17}}\right| with and without considering input correlation is compared in Figs. 11 and 12 respectively in the case of the mesh distribution system.

The consideration of the base case correlation has increased the variance of P_L,16–17 and \left|V_{\mathrm{17}}\right|. A close observation of Figs. 8 to 12 indicates that ECM clearly distinguishes uncorrelated and correlated cases. Further, the effect of input correlation has a noticeable effect on the tail region of probability distribution of the result variables of both systems. Depending on the sign of sensitivity coefficients, the effect of input correlation either elongates or shortens the tail region of the distributions of result variables. The computational efficiency and accuracy (error indexe_X) of ECM is compared with that of MCS in Table 7. The consideration of correlation takes slightly more computational time as compared to the uncorrelated case. The computational effort of ECM is nearly 5.9% and 27.5% of the MCS time respectively for the radial and the mesh system without considering input correlation but nearly 6% and 28.5% of the MCS time after considering base case correlation.

Conclusions

An analytical method ECM is effectively implemented on a modified Ward-Hale 6-bus system and an IEEE 14-bus transmission system as well as a modified IEEE 69-bus radial and an IEEE 33-bus mesh distribution system. Sensitivity matrix-based PLF models for both transmission and distribution systems are clearly explained. The uncertainty pertaining to PV generation, PHEV and BEV charging demands for both residential community and charging stations are modeled probabilistically. Multiple input correlation cases are effectively incorporated. The comparison of the results of ECM with those of MCS has clearly established the fact that ECM has accurately approximated the distributions of the result variables in considerably less computational effort.

References

[1] Borkowska B. Probabilistic load flow. IEEE Transactions on Power Apparatus and Systems, 1974, PAS-93(3): 752–759

[2] Prusty B R, Jena  D. A critical review on probabilistic load flow studies in uncertainty constrained power systems with photovoltaic generation and a new approach. Renewable & Sustainable Energy Reviews, 2017, 69: 1286–1302

[3] Zhang P, Lee  S T. Probabilistic load flow computation using the method of combined cumulants and Gram-Charlier expansion. IEEE Transactions on Power Systems, 2004, 19(1): 676–682

[4] Usaola J. Probabilistic load flow in systems with wind generation. IET Generation, Transmission & Distribution, 2009, 3(12): 1031–1041

[5] Fan M, Vittal  V, Heydt G T ,  Ayyanar R . Probabilistic power flow studies for transmission systems with photovoltaic generation using cumulants. IEEE Transactions on Power Systems, 2012, 27(4): 2251–2261

[6] Villanueva D, Feijóo  A E, Pazos  J L. An analytical method to solve the probabilistic load flow considering load demand correlation using the DC load flow. Electric Power Systems Research, 2014, 110: 1–8

[7] Prusty B R, Jena  D. Modeling of correlated photovoltaic generations and load demands in probabilistic load flow. In: 12th IEEE International Conference Electronics, Energy, Environment, Communication, Computer, Control (INDICON 2015), New Delhi, 2015, 1–6

[8] Prusty B R, Jena  D. Combined cumulant and Gaussian mixture approximation for correlated probabilistic load flow studies: a new approach. CSEE Journal of Power and Energy Systems, 2016, 2(2): 71–78

[9] Conti S, Raiti  S. Probabilistic load flow using Monte Carlo techniques for distribution networks with photovoltaic generators. Solar Energy, 2007, 81(12): 1473–1481

[10] Ruiz-Rodriguez F J ,  Hernández J C ,  Jurado F . Probabilistic load flow for photovoltaic distributed generation using the Cornish–Fisher expansion. Electric Power Systems Research, 2012, 89: 129–138

[11] Carpinelli G, Caramia  P, Varilone P . Multi-linear Monte Carlo simulation method for probabilistic load flow of distribution systems with wind and photovoltaic generation systems. Renewable Energy, 2015, 76: 283–295

[12] Vlachogiannis J G . Probabilistic constrained load flow considering integration of wind power generation and electric vehicles. IEEE Transactions on Power Systems, 2009, 24(4): 1808–1817

[13] Pashajavid E, Golkar  M A. Charging of plug-in electric vehicles: stochastic modelling of load demand within domestic grids. In: Proceedings of 20th Iranian Conference on Electrical Engineering (ICEE2012), Tehran, 2012, 535–539

[14] Kong S, Cho  H C, Lee  J U, Joo  S K. Probabilistic modeling of electric vehicle charging load for probabilistic load flow. In: 2012 IEEE Vehicle Power and Propulsion Conference (VPPC 2012), Seoul, Republic of Korea, 2012, 1010–1013

[15] Wu C, Wen  F, Lou Y ,  Xin F. Probabilistic load flow analysis of photovoltaic generation system with plug-in electric vehicles. International Journal of Electrical Power & Energy Systems, 2015, 64: 1221–1228

[16] Papadopoulos P, Skarvelis-Kazakos  S, Grau I ,  Awad B, Cipcigan  L M, Jenkins  N. Impact of residential charging of electric vehicles on distribution networks, a probabilistic approach. In: Proceedings of 2010 45th International Universities’ Power Engineering Conference (UPEC), Cardiff, Wales, 2010, 1–5

[17] Tehrani N H, Wang  P. Probabilistic estimation of plug-in electric vehicles charging load profile. Electric Power Systems Research, 2015, 124: 133–143

[18] Li G, Zhang  X P. Modeling of plug-in hybrid electric vehicle charging demand in probabilistic power flow calculations. IEEE Transactions on Smart Grid, 2012, 3(1): 492–499

[19] Dimitrovski A, Ackovski  R. Probabilistic load flow in radial distribution networks. In: Proceedings of the 1996 14th IEEE Transmission and Distribution Conference, Los Angeles, CA, 1996, 102–107

[20] Golkar M A. A new probabilistic load-flow method for radial distribution networks. European Transactions on Electrical Power Systems, 2003, 13(3): 167–172

[21] Hoese A, Garcés  F. Stochastic correlated simulation: an extension of the cumulant method to include time-dependent energy sources. International Journal of Electrical Power & Energy Systems, 1999, 21(1): 13–22

[22] Cai D, Shi  D, Chen J . Probabilistic load flow computation with polynomial normal transformation and Latin hypercube sampling. IET Generation, Transmission & Distribution, 2013, 7(5): 474–482 doi:10.1049/iet-gtd.2012.0405

[23] Pai M A, Chatterjee  D. Computer Techniques in Power System Analysis. 3rd ed. Noida: McGraw Hill Education (India) Private Limited, 2014

[24] Electrical Engineering, University of Washington. Power systems test case archive, 2016-09-26.

[25] Kumar D, Agrawal  S. Load flow solution for meshed distribution networks. Dissertation for the Bachelor’s Degree. Rourkela: National Institute of Technology Rourkela, 2013

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg