Probabilistic Power Flow Calculation using Non-intrusive Low-rank Approximation Method

TL;DR

This paper introduces a non-intrusive low-rank approximation method for probabilistic power flow analysis, providing accurate results with less computational effort, especially for high-dimensional problems involving renewable energy variability.

Contribution

The paper presents a novel LRA-based approach that efficiently estimates probabilistic power flow solutions, outperforming traditional methods in high-dimensional scenarios.

Findings

Accurately estimates probabilistic characteristics with less computation.

Performs well on large-scale power systems like 1354-bus.

Handles high-dimensional random inputs effectively.

Abstract

In this paper, a novel non-intrusive probabilistic power flow (PPF) analysis method based on the low-rank approximation (LRA) is proposed, which can accurately and efficiently estimate the probabilistic characteristics (e.g., mean, variance, probability density function) of the PPF solutions. This method aims at building up a statistically-equivalent surrogate for the PPF solutions through a small number of power flow evaluations. By exploiting the retained tensor-product form of the univariate polynomial basis, a sequential correction-updating scheme is applied, making the total number of unknowns to be linear rather than exponential to the number of random inputs. Consequently, the LRA method is particularly promising for dealing with high-dimensional problems with a large number of random inputs. Numerical studies on the IEEE 39-bus, 118-bus, and 1354-bus systems show that the proposed method can achieve accurate probabilistic characteristics of the PPF solutions with much less computational effort compared to the Monte Carlo simulations. Even compared to the polynomial chaos expansion method, the LRA method can achieve comparable accuracy, while the LRA method is more capable of handling higher-dimensional problems. Moreover, numerical results reveal that the randomness brought about by the renewable energy resources and loads may inevitably affect the feasibility of dispatch/planning schemes.

Tables

Table 1. TABLE I: Standard Forms of Classical Continuous Distributions and their Corresponding Orthogonal Polynomials [ 38 ]

Distribution	Density Function	Polynomial	Support
Normal	$\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}$	Hermite	(-$\infty$,$\infty$)
Uniform	$\frac{1}{2}$	Legendre	[-1,1]
Beta	$\frac{{\left(1-x\right)}^{\alpha }{\left(1+x\right)}^{\beta }}{2^{\alpha +\beta +1}B(\alpha +1,\beta +1)}$	Jacobi	[-1,1]
Exponential	$e^{-x}$	Laguerre	(0,$\infty$)
Gamma	$\frac{x^{\alpha }{e}^{-x}}{\Gamma \left(\alpha +1\right)}$	Generalized Laguerre	[0,$\infty$)

Table 2. TABLE II: Comparison of the estimated mean of the bus voltage, branch flow and generator reactive power

$V_{im}/S_{ij}$	${\mu }_{mcs}$	${\mu }_{pce}$	${\mu }_{lra}$	$\frac{\Delta {\mu }_{pce}}{{\mu }_{mcs}}\%$	$\frac{\Delta {\mu }_{lra}}{{\mu }_{mcs}}\%$
$V_{8,m}$	0.9804	0.9801	0.9802	-0.0322	-0.0241
$V_{7,m}$	0.9805	0.9802	0.9803	-0.0333	-0.0248
$S_{6-11}$	2.2189	2.2247	2.2436	0.2627	1.1129
$S_{4-5}$	4.2638	4.2692	4.2716	0.1268	0.1824
$S_{10-13}$	4.6271	4.6324	4.6298	0.1145	0.0577
$S_{13-14}$	4.6821	4.6892	4.6853	0.1510	0.0670
$Q_{32}$	2.7018	2.7186	2.7078	0.6214	0.2213
$Q_{36}$	1.1593	1.1584	1.1605	-0.0757	0.1062

Table 3. TABLE III: Comparison of the estimated standard deviation of the bus voltage, branch flow and generator reactive power

$V_{im}/S_{ij}$	${\sigma }_{mcs}$	${\sigma }_{pce}$	${\sigma }_{lra}$	$\frac{\Delta {\sigma }_{pce}}{{\sigma }_{mcs}}\%$	$\frac{\Delta {\sigma }_{lra}}{{\sigma }_{mcs}}\%$
$V_{8,m}$	0.0155	0.0141	0.0154	-8.9862	-0.5056
$V_{7,m}$	0.0158	0.0144	0.0157	-9.0738	-0.5686
$S_{6-11}$	0.8112	0.6962	0.8210	-14.1836	1.1999
$S_{4-5}$	1.4085	1.3331	1.4131	-5.3578	0.3225
$S_{10-13}$	0.7968	0.7652	0.8034	-3.9612	0.8265
$S_{13-14}$	0.8835	0.8446	0.8899	-4.4042	0.7202
$Q_{32}$	0.4899	0.4321	0.4859	-11.8063	-0.8291
$Q_{36}$	0.1057	0.0956	0.1043	-9.5410	-1.3486

Table 4. TABLE IV: A Comparison of the estimated mean of the bus voltage and the branch flow by different methods

$V_{im}/S_{ij}$	${\mu }_{mcs}$	${\mu }_{pce}$	${\mu }_{lra}$	$\frac{\Delta {\mu }_{pce}}{{\mu }_{mcs}}\%$	$\frac{\Delta {\mu }_{lra}}{{\mu }_{mcs}}\%$
$V_{53,m}$	0.9412	0.9412	0.9412	-0.0002	-0.0002
$V_{21,m}$	0.9435	0.9435	0.9435	0.0008	-0.0005
$V_{44,m}$	0.9530	0.9530	0.9530	0.0049	0.0029
$V_{20,m}$	0.9465	0.9465	0.9465	0.0004	-0.0003
$S_{49-69}$	0.6652	0.6631	0.6676	-0.3243	0.3595
$S_{47-69}$	0.7716	0.7684	0.7729	-0.4229	0.1614

Table 5. TABLE V: A Comparison of the estimated standard deviation of the bus voltage and the branch flow by different methods

$V_{im}/S_{ij}$	${\sigma }_{mcs}$	${\sigma }_{pce}$	${\sigma }_{lra}$	$\frac{\Delta {\sigma }_{pce}}{{\sigma }_{mcs}}\%$	$\frac{\Delta {\sigma }_{lra}}{{\sigma }_{mcs}}\%$
$V_{53,m}$	0.0021	0.0021	0.0021	-0.1653	-0.0370
$V_{21,m}$	0.0036	0.0036	0.0036	-1.6728	0.2189
$V_{44,m}$	0.0067	0.0066	0.0068	-2.3151	0.2594
$V_{20,m}$	0.0025	0.0025	0.0025	-1.6444	0.1961
$S_{49-69}$	0.2694	0.2619	0.2697	-2.8034	0.1002
$S_{47-69}$	0.2718	0.2565	0.2719	-5.6293	0.0091

Table 6. TABLE VI: A Comparison of the Computation Time Required by the LHS-based MCS, the PCE and the LRA

Method	$t_{ed}\left(s\right)$	$t_{sc}\left(s\right)$	$t_{es}\left(s\right)$	$t_{total}\left(s\right)$
MCS	–	–	653.94	653.94
PCE	28.84	73.96	0.29	103.09
LRA	28.84	46.53	1.39	76.76

Table 7. TABLE VII: A Comparison of the estimated mean and standard deviation of the branch flow by different methods

$V_{im}/S_{ij}$	${\mu }_{mcs}$	${\mu }_{lra}$	$\frac{\Delta {\mu }_{lra}}{{\mu }_{mcs}}\%$	${\sigma }_{mcs}$	${\sigma }_{lra}$	$\frac{\Delta {\sigma }_{lra}}{{\sigma }_{mcs}}\%$
$S_{118-4598}$	8.6613	8.6616	0.0032	0.6782	0.6779	-0.0416
$S_{1754-960}$	7.6526	7.6524	-0.0030	0.5192	0.5197	0.0924
$S_{6901-4874}$	5.6797	5.6797	-0.0004	0.1744	0.1745	0.0891
$S_{7267-6581}$	7.1229	7.1227	-0.0019	0.1766	0.1760	-0.3457
$S_{1001-3580}$	5.7043	5.7041	-0.0027	0.1596	0.1599	0.1500
$S_{1001-516}$	5.7132	5.7131	-0.0011	0.1179	0.1184	0.3920
$S_{1758-1923}$	7.9033	7.9034	0.0004	0.1939	0.1934	-0.2739

Table 8. TABLE VIII: A Comparison of the estimated line flow values at the 10% and 90% confidence levels by different methods

$S_{ij}$	$S_{mcs}^{10}$	$S_{lra}^{10}$	$\frac{\Delta S_{lra}^{10}}{S_{mcs}^{10}}\%$	$S_{mcs}^{90}$	$S_{lra}^{90}$	$\frac{\Delta S_{lra}^{90}}{S_{mcs}^{90}}\%$
$S_{118-4598}$	7.8045	7.8116	0.0911	9.5266	9.5349	0.0873
$S_{1754-960}$	6.9978	7.0051	0.1042	8.3152	8.3198	0.0561
$S_{6901-4874}$	5.4569	5.4584	0.0282	5.9048	5.9057	0.0139
$S_{7267-6581}$	6.8885	6.8941	0.0809	7.3523	7.3501	-0.0306
$S_{1001-3580}$	5.4841	5.4957	0.2099	5.8998	5.9061	0.1067
$S_{1001-516}$	5.5529	5.5564	0.0633	5.8587	5.8627	0.0683
$S_{1758-1923}$	7.6511	7.6560	0.0638	8.1527	8.1526	-0.0009

Table 9. TABLE IX: Comparison of the unknowns to be solved by the PCE and the LRA for different number of inputs

$n$	Sparse PCE				LRA
	$p_{max}^{pce}$	$N_{full}^{pce}$	$N_{trunc}^{pce}$	$N_{solved}^{pce}$	$p_{max}^{lra}$	$r^{lra}$	$N_{solved}^{lra}$
100	3	176851	5251	89	2	1	301
200	3	1373701	20501	161	2	1	601
300	4	348881876	135751	200	2	1	901
400	3	10827401	81001	214	2	1	1201
500	4	2656615626	376251	339	2	1	1501
600	-	-	-	-	2	1	1801
713	-	-	-	-	2	1	2140

Table 10. TABLE X: A Comparison of the estimated mean of the line flow of branch 1758-1923 in different cases

$n$	${\mu }_{mcs}$	${\mu }_{pce}$	${\mu }_{lra}$	$\frac{\Delta {\mu }_{pce}}{{\mu }_{mcs}}\%$	$\frac{\Delta {\mu }_{lra}}{{\mu }_{mcs}}\%$
100	7.9026	7.9030	7.9027	0.0044	0.0008
200	7.9027	7.9022	7.9027	-0.0059	-0.0002
300	7.9029	7.9029	7.9030	0.0006	0.0007
400	7.9030	7.9030	7.9031	0.0001	0.0014
500	7.9032	7.9032	7.9033	-0.0005	0.0011
600	7.9034	-	7.9033	-	-0.0002
713	7.9033	-	7.9034	-	0.0004

Table 11. TABLE XI: A Comparison of the estimated standard deviation of the line flow of branch 1758-1923 in different cases

$n$	${\sigma }_{mcs}$	${\sigma }_{pce}$	${\sigma }_{lra}$	$\frac{\Delta {\sigma }_{pce}}{{\sigma }_{mcs}}\%$	$\frac{\Delta {\sigma }_{lra}}{{\sigma }_{mcs}}\%$
100	0.1546	0.1501	0.1547	-2.9228	0.0594
200	0.1793	0.1731	0.1795	-3.4811	0.1250
300	0.1892	0.1824	0.1890	-3.5640	-0.1099
400	0.1952	0.1900	0.1941	-2.6852	-0.5928
500	0.1948	0.1897	0.1947	-2.6008	-0.0623
600	0.1961	-	0.1962	-	0.0447
713	0.1939	-	0.1934	-	-0.2739

Equations

$$f_V\left(v\right)=\frac{k}{c}{\left(\frac{v}{c}\right)}^{k-1}\exp \left[-{\left(\frac{v}{c}\right)}^k\right]$$

$${{P}_{w}}(v)=\left\{\begin{array}[]{*{35}{l}}0&v\leq{{v}_{in}}\mbox{ or }v>{{v}_{out}}\\ \displaystyle\frac{v-{{v}_{in}}}{{{v}_{rated}}-{{v}_{in}}}{{P}_{r}}&{{v}_{in}}<v\leq{{v}_{rated}}\\ {{P}_{r}}&{{v}_{rated}}<v\leq{{v}_{out}}\\ \end{array}\right.$$

$$f_R\left(r\right)=\frac{\Gamma (\alpha +\beta )}{\Gamma (\alpha )\Gamma (\beta )}{\left(\frac{r}{r_{\max }}\right)}^{\alpha -1}{\left(1-\frac{r}{r_{\max }}\right)}^{\beta -1}$$

$${{P}_{pv}}(r)=\left\{\begin{array}[]{*{35}{l}}\displaystyle\frac{{{r}^{2}}}{{{r}_{c}}{{r}_{std}}}{{P}_{r}}&0\leq r<{{r}_{c}}\\ \displaystyle\frac{r}{{{r}_{std}}}{{P}_{r}}&{{r}_{c}}\leq r\leq{{r}_{std}}\\ {{P}_{r}}&r>{{r}_{std}}\\ \end{array}\right.$$

$$\mathbf{f}(\mathbf{x})=\left[\begin{array}{c}{P}_{Gi}-P_{Li}-P_i(\mathbf{x})\\ Q_{Gi}-Q_{Li}-Q_i(\mathbf{x})\end{array}\right]=0$$

$$\begin{array}{c}{P}_i(\mathbf{x})=V_i\sum _{j=1}^N{V}_j(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij})\\ Q_i(\mathbf{x})=V_i\sum _{j=1}^N{V}_j(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij})\end{array}$$

$$\begin{array}{c}{P}_{Gi}+P_{wi}(v_i)+P_{pvi}(r_i)-P_{Li}(P_{Li})-P_i(\mathbf{x})=0\\ Q_{Gi}+Q_{wi}(v_i)-Q_{Li}(P_{Li})-Q_i(\mathbf{x})=0\end{array}$$

$$\begin{array}{c}{P}_{Gi}+P_{wi}(v_i)+P_{pvi}(r_i)-P_{Li}(P_{Li})-P_i(\mathbf{x})=0\\ V_i=V_{i0}\\ Q_{Gi}=-Q_{wi}(v_i)+Q_{Li}(P_{Li})+Q_i(\mathbf{x})\\ Q_{min,i}\le Q_{Gi}\le Q_{max,i}\end{array}$$

$$\begin{array}{c}\mathbf{f}\left(\mathbf{x},\mathbf{U}\right)=0\end{array}$$

$$Y\approx \hat{Y}=\hat{g}\left(\mathbf{\xi }\right)=\sum _{l=1}^r{b}_l{\omega }_l\left(\mathbf{\xi }\right)$$

$${\omega }_l(\mathbf{\xi })=\prod _{i=1}^n{v}_l^{(i)}({\xi }_i)$$

$$\hat{Y}=\hat{g}(\mathbf{\xi })=\sum _{l=1}^r{b}_l\left[\prod _{i=1}^n\left(\sum _{k=0}^{p_i}{z}_{k,l}^{\left(i\right)}{\varphi }_k^{(i)}\left({\xi }_i\right)\right)\right]$$

$$\begin{array}{c}{\omega }_r(\xi )=v_r^{(1)}({\xi }_1)\times v_r^{(2)}({\xi }_2)\times v_r^{(3)}({\xi }_3)\\ v_r^{(1)}({\xi }_1)=z_{0,r}^{(1)}{\varphi }_0^{(1)}({\xi }_1)+z_{1,r}^{(1)}{\varphi }_1^{(1)}({\xi }_1)+z_{2,r}^{(1)}{\varphi }_2^{(1)}({\xi }_1)\\ v_r^{(2)}({\xi }_2)=z_{0,r}^{(2)}{\varphi }_0^{(2)}({\xi }_2)+z_{1,r}^{(2)}{\varphi }_1^{(2)}({\xi }_2)+z_{2,r}^{(2)}{\varphi }_2^{(2)}({\xi }_2)\\ v_r^{(3)}({\xi }_3)=z_{0,r}^{(3)}{\varphi }_0^{(3)}({\xi }_3)+z_{1,r}^{(3)}{\varphi }_1^{(3)}({\xi }_3)+z_{2,r}^{(3)}{\varphi }_2^{(3)}({\xi }_3)\end{array}$$

$$\begin{array}{c}\sqrt{b_{k+1}}{\tilde{\pi }}_{k+1}\left({\xi }_i\right)=\left({\xi }_i-d_k\right){\tilde{\pi }}_k\left({\xi }_i\right)-\sqrt{b_k}{\tilde{\pi }}_{k-1}\left({\xi }_i\right)\\ {{\tilde{\pi }}_{-1}\left({\xi }_i\right)=0,\tilde{\pi }}_0\left({\xi }_i\right)=1\end{array}$$

$$\begin{array}{c}{d}_k=\frac{⟨\xi {\pi }_k,{\pi }_k⟩}{⟨{\pi }_k,{\pi }_k⟩},k\ge 0\\ b_0=⟨{\pi }_0,{\pi }_0⟩,b_k=\frac{⟨{\pi }_k,{\pi }_k⟩}{⟨{\pi }_{k-1},{\pi }_{k-1}⟩},k\ge 1\end{array}$$

$${\omega }_r(\mathbf{\xi })$$

$$=\arg \underset{\omega \in W}{\min }\sum _{m=1}^{M_C}{\left[y^{(m)}-{\hat{g}}_{r-1}({\mathbf{\xi }}^{(m)})-\omega ({\mathbf{\xi }}^{(m)})\right]}^2$$

$${\mathbf{z}}_r^{\left(i\right)}=\arg \underset{\mathbf{\zeta }\in R^{(p_i+1)}}{\min }{\Vert e_{r-1}-C_i\left(\sum _{k=0}^{p_i}{\zeta }_k{\varphi }_k^{\left(i\right)}\right)\Vert }_{{\mathbf{\xi }}_C}^2$$

$$C_i=\prod _{j\ne i}{v}_r^{(j)}({\xi }_j)=\prod _{j\ne i}\left(\sum _{k=0}^{p_j}{z}_{k,r}^{\left(j\right)}{\varphi }_k^{\left(j\right)}\left({\xi }_j\right)\right)$$

$$\mathbf{b}=\arg \underset{\mathbf{\beta }\in R^r}{\min }{\Vert g-\sum _{l=1}^r{\beta }_l{\omega }_l\Vert }_{{\mathbf{\xi }}_C}^2$$

$${\hat{e}}_r=\frac{{\Vert e_{r-1}-{\omega }_r\Vert }_{{\mathbf{\xi }}_C}^2}{\mathbb{V}({\mathbf{y}}_C)}$$

$${\mu }_y=\mathbb{E}\left[\hat{g}\left(\mathbf{\xi }\right)\right]=\sum _{l=1}^r{b}_l\left(\prod _{i=1}^n{z}_{0,l}^{\left(i\right)}\right)$$

$${\sigma }_y^2=\sum _{l=1}^r\sum _{m=1}^r{b}_l{b}_m\prod _{i=1}^n\left[\left(\sum _{k=0}^{p_i}{z}_{k,l}^{\left(i\right)}{z}_{k,m}^{\left(i\right)}\right)-z_{0,l}^{\left(i\right)}{z}_{0,m}^{\left(i\right)}\right]$$

$$\begin{array}{c}\frac{\Delta {\mu }_{\ast }}{{\mu }_{mcs}}\%=\left|\frac{{\mu }_{\ast }-{\mu }_{mcs}}{{\mu }_{mcs}}\right|\times 100\%\\ \frac{\Delta {\sigma }_{\ast }}{{\sigma }_{mcs}}\%=\left|\frac{{\sigma }_{\ast }-{\sigma }_{mcs}}{{\sigma }_{mcs}}\right|\times 100\%\end{array}$$

Full text

I Introduction

Power flow analysis is undoubtedly a fundamental and essential tool in power system operation and planning. It is employed to determine the system operating state, to compare different dispatching/planning schemes, and to provide the initial condition for other advanced applications. The growing integration of wind farms and utility-scale solar photo-voltaic (PV) power plants in transmission systems inevitably results in increasing degree of uncertainties to the operating state of a power system. Therefore, adopting a probabilistic framework is essential to ensure realistic and accurate power flow analysis results.

Indeed, the proposal of probabilistic power flow (PPF) problem can be tracked back to 1970s, the aim of which is to incorporate uncertainties into the power flow problem [1], [2]. Many efforts were made thereafter. These methods can be generally divided into two categories: simulation methods and analytical methods. Monte Carlo simulations (MCS) and Monte-Carlo like simulations [3], [4], [5] are the most widely used in the first category due to its simplicity, yet its computational effort grows exponentially as the random inputs increases.

Attempting to release the computational burden, some analytical methods are developed, which utilize mathematical approximations and assumptions. The convolution method proposed in [2] applies the convolution to obtain the probability density function (PDF) of the power flow solutions. However, the computational burden is still high due to the time-consuming convolution. Besides, the cumulant method was proposed in [6], which can indirectly estimate the moments of responses. Nevertheless, its accuracy will degrade if the random inputs have large variations. Another representative method is the point estimation method [7], [8], which is computationally efficient yet its accuracy will deteriorate as the dimension of input random variables or the order of moments increases. In addition, both the cumulant method and point estimation method can not provide PDF and cumulative distribution function (CDF) of the responses directly, and hence series expansions were proposed to work together with them including Gram-Charlier series [6], Cornish-Fisher series [9], and Edgeworth series [10]. However, it has been shown in [10] that the convergence of the series expansions can not always be guaranteed.

Another class of methods for uncertainty quantification is meta-modeling, which aims to build a statistically-equivalent analytical representation for the desired responses (e.g., the PPF solutions in this study) using a small number of model evaluations. Particularly, the polynomial chaos expansion (PCE) is the most popular one because (i) it has strong mathematical basis; (ii) it works with deterministic tools in a non-intrusive fashion; (iii) it can provide accurate and comprehensive statistical properties of responses with low computational effort. The PCE method has been applied in the context of power systems to study the probabilistic power flow [11], [12], the load margin problem [13], and the available delivery capability problem [14].

An emerging method, alternative to PCE, is the canonical low-rank approximations (LRA), which employs the canonical decomposition to express the desired response as a sum of rank-one functions [15]. The original idea of canonical decomposition dates back to 1927 [16] and becomes an attractive method for uncertainty quantification of structural vibration problems recently [17]-[18]. By exploiting the tensor-product structure of the multivariate polynomial basis, LRA can provide polynomial representations in highly compressed formats. The outstanding advantage of the LRA over the PCE is that the number of coefficients to be solved grows linearly rather than exponentially with the number of inputs, making it more powerful on dealing with high dimensional problems.

In this paper, we apply the LRA approach to solve the probabilistic power flow problem in which various uncertainties from the wind power, the solar PV and the loads are incorporated. To the knowledge of authors, the paper seems to be the first attempt to apply the LRA method in the context of power system. The contributions of the paper are as below:

• •

A computationally efficient yet accurate algorithm is proposed to evaluate the PPF problem. Particularly, in the proposed LRA method:

• –

Random variables with diverse marginal distributions can be accommodated using proper polynomial basis. The correlations among them can be incorporated using the Nataf transformation.

• –

A sequential correction-updating scheme is employed to build up the LRA model, making the total number of unknown coefficients to be linear rather than exponential to the number of random inputs.

• –

The mean and variance of the bus voltages and line flows can be obtained analytically without evaluating a large number of samples.

• •

Accurate probabilistic characteristics of the bus voltages and line flows can be achieved by the proposed methodology with much less computational effort compared to the Latin hypercube sampling (LHS)-based MCS.

• •

A higher capability of dealing high-dimensional problems can be reached compared to the PCE method because the number of unknown coefficients grows only linearly with the number of random inputs.

It should be noted that a relevant problem ”optimal power flow (OPF)”, originating from the economic dispatch problem [19], [20] is fundamentally different from the power flow or PPF study considered in this paper. The power flow aims to determine the operating state, i.e., bus voltages and line flows, of the power system under given loads, generations and network conditions. It is formulated as a set of nonlinear equations. The PPF extends the power flow to assess the probability characteristics of the operating state for a range of loads, generations and network conditions. In contrast, the OPF aims to determine the optimal operating state of the power system by simultaneously minimizing a given objective function and satisfying certain physical and operating constraints. It is typically formulated as an optimization problem. The probabilistic OPF [21] and chance-constrained OPF [22] are two relative terms of PPF in the probabilistic framework. They have been tailored for various applications in both transmission and distribution systems, e.g., decentralized energy trading [23]. In the rest of the paper, we focus only on PPF.

The remaining of the paper is organized as follows. Section II introduces the probabilistic power system models. Section III describes the mathematical formulation of PPF problem. Section IV elaborates the low-rank approximation method and its implementation in probabilistic power flow analysis. The detailed algorithm to assess the PPF is summarized in Section V. The simulation results on the modified IEEE 39-bus, 118-bus and 1354-bus systems are shown in Section VI. Conclusions and perspectives are given in Section VII.

II The Probabilistic Power System Model

Power systems are subjected to various uncertainties and randomness from different sources, ranging from renewable generations, load variations, topology changes, to unexpected outages and faults. In this study, we focus on the randomness brought about by the renewable generators and loads. The other sources of uncertainty can also be integrated into the formulation if needed. Correspondingly, the random inputs of interests are wind speed ${v}$, solar radiation ${r}$ and load power ${P_{L}}$. Each random input typically can be expressed as a random variable associated with a probabilistic density function (PDF) ${X\sim f_{X}(x)}$.

II-A Wind Generation

In long-time scale, the wind speed in many locations around the world can be represented by a Weibull distribution [24], [25], [26], the probability density function of which is

[TABLE]

where $v$ is the wind speed, $k$ and $c$ are the shape and scale parameters, respectively. Several methods have been proposed to extract the $k$ and $c$ from historical wind speed data, such as the moment method [27] and the maximum likelihood method [28]. In short-time scale, however, the wind speed can be modeled as a Gaussian distribution [29] in which the mean value is equal to the forecasted wind speed and the variance represents the forecasted error.

As a result, the active power output ${{P}_{w}}$ can be calculated by the piece-wise wind speed-power output relation [30]

[TABLE]

where ${{v}_{in}}$, ${{v}_{out}}$ and ${{v}_{rated}}$ are the cut-in, cut-out, and rated wind speed ($m/s$), ${{P}_{r}}$ is the rated wind power ($kW$). Once the active power is obtained, the reactive power can be determined according to the speed control type of the wind turbine [31] since the wind turbine can be modelled as either a constant P-Q bus or a constant P-V bus with given Q-limits.

II-B Solar Generation

Typically, the solar radiation in long time scale is represented by the Beta distribution [24]

[TABLE]

where $\alpha$ and $\beta$ are the shape parameters of the distribution, which can be fitted from historical solar radiation data by several methods, such as the moment method [32] and the maximum likelihood method [33]. $\Gamma$ denotes the Gamma function, $r$ and ${{r}_{\max}}$ (${W/m}^{2}$) are the respective actual and maximum solar radiations. Similar to the wind speed, the solar radiation in short-time scale can be modeled as a Gaussian distribution provided that accurate mean value is available [29].

The active power ${{P}_{pv}}$ corresponding to the solar radiation $r$ is determined by the following piece-wise function [24]

[TABLE]

where ${{r}_{c}}$ is a certain radiation point set usually as 150 $W/m^{2}$, ${{r}_{std}}$ is the solar radiation in the standard environment, ${{P}_{r}}$ is the rated power of the solar PV. The reactive power ${{Q}_{pv}}$ is assumed to be zero in this study according to the solar generation is required to inject into the power grid at unit power factor [34]. However, similar to the wind farm, the PV power plant can be modeled as P-Q bus or P-V bus with Q-limits in PPF once similar interconnection requirements are formulated and applied in the near future [35].

II-C Load Variation

By nature, load demand is uncertain in power systems. It is a common practice to model the load uncertainty by Gaussian distribution with specified mean value ${{\mu}_{P}}$ and variance ${{\sigma}_{P}}$ [1], [36], which are typically provided by the load forecaster and historical data, respectively. Traditionally, only the active power is predicted, whereas the reactive power is calculated under the assumption of constant power factor [37].

III Mathematical Formulation of Probabilistic Power Flow

The power flow equations can be represented as

[TABLE]

where

[TABLE]

and ${\bm{x}}={{[{\theta},{V}]}^{T}}$, e.g., voltage angles and magnitudes for all buses; $G_{ij}$ and $B_{ij}$ are the real and imaginary part of the entry $Y_{ij}$ in the bus admittance matrix.

Let $v$, $r$ and ${{P}_{L}}$ be the random vectors that represent wind speeds, solar radiations and load variations, respectively, the probabilistic power flow (PPF) equations of a $N$-bus system can be described as below. Specifically, for P-Q type buses, the PPF equations are:

[TABLE]

For P-V type buses, the corresponding PPF equations are:

[TABLE]

where $P_{wi}({{v}_{i}})$, $P_{pvi}({{r}_{i}})$, $P_{Li}$ and $P_{Gi}$ are the real power injection from the wind turbine, the solar PV, the load, and the conventional generator at bus $i$; $Q_{wi}({{v}_{i}})$, $Q_{Li}$ and $Q_{Gi}$ are the reactive power injection from the wind power, the load, and the conventional generator at bus $i$. If $Q_{Gi}$ exceeds its limits, e,g., ${{Q}_{min,i}}$ or ${{Q}_{max,i}}$, then the terminal bus switches from P-V to P-Q with $Q_{Gi}$ fixed at the violated limit.

In fact, the set of PPF equations (7)-(8) can be described in the following compact form:

[TABLE]

where ${\bm{x}}={{[{\theta},{V}]}^{T}}$ is the state vector, ${\bm{U}=\left[\bm{v},\bm{r},\bm{{P}_{L}}\right]}$ is the random vector describing the wind speed, the solar radiation, and the load power.

IV Canonical Low-rank Approximation using Polynomial Basis

This section presents the general framework of low-rank approximation of a multivariate stochastic response function. For simplicity, we first consider a scalar response function of independent inputs, the case of dependent inputs will be addressed in Section IV–E.

IV-A Low-rank Approximation with Polynomial Basis

Consider a random vector ${\bm{\xi}}$ = (${{\xi}_{1},{\xi}_{2},...,{\xi}_{n},}$) with joint probability density function (PDF) ${f_{\bm{\xi}}}$ and marginal distribution functions ${f_{\xi_{i}},i=1,...,n}$ (${\xi_{i}}$ is related with the random variables ${U_{i}}$ in (9), see Section III), then the canonical rank-$r$ approximation [15] of the target stochastic response (e.g., probabilistic bus voltages or branch flows in this study) ${Y=g(\bm{\xi})}$ can be represented by:

[TABLE]

in which $b_{l},l=1,...,r$ are normalizing weighting factors, and ${{\omega}_{l}}$ is a rank-one function of $\bm{\xi}$ in the form of

[TABLE]

where ${v_{l}^{(i)}}$ denotes the ${i}$-th dimensional univariate function in the ${l}$-th rank-one function. For most applications, the number ${r}$ of rank-one terms is usually small (under 5), hence (10) and (11) represent a canonical low-rank approximation.

In order to obtain the rank-$r$ approximation, a natural choice is expanding ${v_{l}^{(i)}}$ on a polynomial basis ${\{\phi_{k}^{(i)},k\in\mathbb{N}\}}$ that is orthogonal to ${f_{\xi_{i}}}$, the resulting rank-$r$ approximation takes the form:

[TABLE]

where ${\phi}_{k}^{(i)}$ denotes the $k$-th degree univariate polynomial in the $i$-th random input, $p_{i}$ is the maximum degree of ${\phi}^{(i)}$ and ${z}_{k,l}^{(i)}$ is the coefficient of ${\phi}_{k}^{(i)}$ in the $l$-th rank-one function.

Building the low-rank approximation for desired responses in (12) requires: (i) choose an appropriate univariate polynomial for each random input; (ii) solve the polynomial coefficients ${z_{k,l}^{(i)}}$ as well as the weighing factors ${b_{l}}$. This process relies on a set of samples of $\bm{\xi}$, which are usually termed as experimental design (ED) and their corresponding accurate responses $\bm{Y}$.

To illustrate the construction of $r$-th rank-one function, consider a 3-dimensional input random vector ${\bm{\xi}=(\xi_{1},\xi_{2},\xi_{3})}$ with mutually independent components ${{{\xi}_{1}}\sim\mbox{Weibull}\left(\lambda,k\right)}$, ${{{\xi}_{2}}\sim\mbox{Beta}\left(\alpha,\beta\right)}$ and ${{{\xi}_{3}}\sim N\left(0,1\right)}$ respectively, then the $r$-th rank-one function $\omega_{r}$ that corresponds to ${p_{i}=2,i=1,2,3}$ is:

[TABLE]

where ${z_{0,r}^{(i)}}$, ${z_{1,r}^{(i)}}$ and ${z_{2,r}^{(i)}}$ are the coefficients corresponding to the zero-, first- and second-order univariate polynomial ${\phi_{0}^{(i)}({{\xi}_{i}})}$, ${\phi_{1}^{(i)}({{\xi}_{i}})}$ and ${\phi_{2}^{(i)}({{\xi}_{i}})}$ respectively. Now we have the $r$-th rank-one function ${\omega_{r}(\xi)}$ with known polynomial basis, the next step is to determine the coefficients $\{{z_{0,r}^{(i)}},{z_{1,r}^{(i)}},{z_{2,r}^{(i)}}:i=1,2,3\}$.

IV-B Selection of the Univariate Polynomial Basis

To reduce the computational effort for building an accurate low-rank approximation (12) for each response of interests, it is crucial to choose a proper polynomial ${\phi}_{i}$ for the $i$-th random input ${\xi}_{i},i=1,2,...,n$; otherwise, lower convergence rate and higher degree of polynomial basis may be needed [38], which will hamper the capability of LRA in dealing with high-dimensional problems.

Table I shows a set of typical continuous distributions and the respective optimal Wiener-Askey polynomial basis [38], which can ensure the exponential convergence rate. In case that $f_{\xi_{i}}$ is out of the list in Table I, two options are available: (i) employ the isoprobabilistic transformation [39] to project $\xi_{i}$ to a typical distribution in Table I, and then choose its corresponding polynomial basis; (ii) adopt the discretized Stieltjes procedure to numerically construct a set of univariate orthogonal polynomial basis in the form: ${\tilde{\pi}=\frac{{{\pi}_{k}}}{\left\langle{{\pi}_{k}},{{\pi}_{k}}\right\rangle}}$ for $\xi_{i}$. In this paper, we carry out the second option to recurrently compute the orthogonal polynomial basis [40]:

[TABLE]

[TABLE]

where $k$ denotes the polynomial degree.

It is worth mentioning that the fundamental premise of applying polynomial basis in Table I or equations (14)-(15) is that the PDF of random inputs are known in advance. This assumption is reasonable if abundant historical data (several years) of wind speed and solar radiation are available, from which the wind speed and solar radiation can be well estimated (for wind [27], [28]; for solar [41], [32]). In case the historical data cannot be well fitted by any typical distribution, the data-driven polynomial chaos method offers an alternative choice to construct the polynomial basis from dataset [42].

Now we have chosen an optimal polynomial basis for each random input. The next step is to determine the coefficients ${z_{k,l}^{(i)}}$ as well as the weighing factors ${b_{l}}$ in (12).

IV-C Calculation of the Coefficients and Weighing Factors

Based on a experiment design of size $M_{C}$, i.e., a set of samples of $\bm{\xi}_{C}=\{\bm{\xi}^{(1)},\bm{\xi}^{(2)},...,\bm{\xi}^{(M_{C})}\}$ and corresponding responses $\bm{y}_{C}=\{y^{(1)},y^{(2)},...,y^{(M_{C})}\}$ which evaluated by deterministic tools, different algorithms have been proposed in the literature for solving the LRA coefficients and weighing factors in a non-intrusive manner [17], [18], [43]. The sequential correction-updating scheme (Algorithm 1) presented in [18] is employed in this study due to its efficiency and capability of constructing low-rank approximation using less sample evaluations. In the ${r}$-th correction step, a new rank-one function ${\omega_{r}}$ is built, while in the ${r}$-th updating step, the set of weighing factors ${\{b_{1},...,b_{r}\}}$ is determined. This process continues until the applied error index stop decreasing [18].

Correction step: the ${r}$-th correction step aims to find a new rank-one tensor $\omega_{r}$, which can be obtained by solving the following minimization problem:

[TABLE]

where ${W}$ represents the space of rank-one tensors, ${{e}_{r-1}}=(g-{{\hat{g}}_{r-1}})$ is the approximation error of the response $Y$ at the ${(r-1)}$-th step, $\|.\|^{2}$ represents the norm 2 of the residual after the new rank-one tensor $w$ is applied, and the subscript ${\bm{\xi}_{C}}$ indicates that the minimization is carried over the whole set of samples in the experiment design $(\bm{\xi}_{C},\bm{y}_{C})$.

By exploiting the retained tensor-product form of the univariate polynomial basis, as shown in (12), typical scheme for solving equation (16) is the alternated least-square (ALS) minimization, which involves sequential minimization along each dimension ${i=1,...,n}$ to solve the corresponding polynomial coefficients ${\bm{z}_{r}^{(i)}=(z_{0,r}^{(i)},...,z_{p_{i},r}^{(i)})}$. The total number of coefficients to be solved in each correction step is $\sum\nolimits_{i=1}^{n}{({p}_{i}+1)}$, which grows linearly as the number of random inputs $n$ increases. Since ${{\omega}_{r}}$ is the product of $v_{r}^{(i)}({\xi}_{i})$ as shown in (11), $v_{r}^{(i)}({\xi}_{i})$ can be initialized as $1.0$.

In the minimization along the ${i}$-th dimension, the polynomial coefficients corresponding to all other dimensions are ”frozen” at their current values and the coefficients ${\bm{z}_{r}^{(i)}=(z_{0,r}^{(i)},...,z_{p_{i},r}^{(i)})}$ can be determined by:

[TABLE]

where $C_{i}$ is a scalar

[TABLE]

Updating step: After the ${r}$-th correction step is completed, the algorithm proceed to the ${r}$-th updating step to determine the weighing factor $b_{r}$ of the newly solved rank-one function ${{\omega}_{r}(\bm{\xi})}$, meanwhile, the set of existing weighing factors ${\bm{b}=(b_{1},...,b_{r-1})}$ are updated too. The updating step can be achieved by solving the following minimization problem

[TABLE]

Stop criteria: The correction-updating scheme successively adds new rank-one function to improve the accuracy of the approximation (12), hence, error reduction in two successive iterations becomes an natural stop criteria for this process. In this study, the relative empirical error is employed which is given by:

[TABLE]

where $\mathbb{V}(\bm{y}_{C})$ denotes the empirical variance of the desired response over the experimental design.

Remark: (i) It is worth pointing out that the minimization problems in (17) and (19) can be efficiently solved with the ordinary least-squares (OLS) technique because the dimension of unknowns are small. (ii) When the LRA (12) for the desired response is built up, the response of any new samples can be evaluated efficiently by directly substituting to (12) instead of solving the original complex problem (e.g., the power flow problem (9)).

IV-D Selection of optimal rank and polynomial degree

Currently, there is no systematic way to identify the optimal rank $r$ and the polynomial degree $p_{i}$ for individual random input beforehand. In this study, we specify a candidate set, say ${\{1,2,3,4,5\}}$, for the rank and another one ${\{2,3,4,5\}}$ for the polynomial degree, and assume the same degree for all the univariate polynomials. The rank selection is performed by progressively increasing the rank and applying the corrected error to select the best one. The selection of optimal degree can be implemented in a similar way. Further study on how to determine an optimal rank and the degree in a systematic and efficient way is desired.

IV-E Integration of Dependent Random Inputs

So far, the random inputs are assumed to be mutually independent as requested by the LRA method. To accommodate dependent random inputs with correlation matrix $\bm{\rho}$, the Nataf transformation [39], [44] and the isoprobabilistic transformation can be employed to build up a mapping between $\bm{U}$ and the independent standard random variables $\bm{\xi}$: $\bm{u}=T(\bm{\xi})$, where $T$ is invertible. Therefore, the set of samples of ${\bm{\xi}}$ can be transformed back into samples of ${\bm{U}}$ to evaluate the corresponding responses ${\bm{Y}}$, after which the desired response $\bm{y}=g(T^{-1}(\bm{u}))$ can be expanded onto the polynomial basis with $\bm{\xi}$ using the aforementioned method [15].

IV-F Moments of a Low-Rank Approximation

Due to the orthogonality of the univariate polynomials that form the LRA basis (see (12)), the mean and variance of the meta-model can be obtained analytically in terms of the polynomial coefficients and the weighing factors. In particular, the mean and variance of the LRA response are given respectively by [15]:

[TABLE]

and

[TABLE]

Hence, if only mean and variance are of interests, (21) and (22) can be applied directly without evaluating a large size of samples which is required by most of simulation-based methods. If other statistical properties (e.g., higher-order moments, PDF, CDF, etc.) are needed, one can sample the random inputs extensively and apply the solved functional approximation (12) to evaluate the corresponding response $\bm{y}$, from which the the statistics of interest can be obtained.

V Computation of Probabilistic Power Flow

In this section, a step-by-step description of the proposed probabilistic power flow calculation is summarized below, which is also illustrated in Fig. 1:

Step 1: Input the network data, the probability distribution and the parameters of the random inputs $\bm{U}$, i.e., the wind speed, the solar radiation, the load active power, and their correlation matrix ${\bm{\rho}}$.

Step 2: Choose the independent standard variable ${\xi_{i}}$ and the corresponding univariate polynomial ${\phi_{i}}$ for each random input $U_{i}$.

Step 3: Generate an experimental design of size ${{M}_{C}}$:

• •

i) Generate ${M_{C}}$ samples ${\bm{\xi}_{C}=(\xi^{(1)},\xi^{(2)},...,\xi^{({M_{C}})})}$ in the standard space by the LHS.

• •

ii) Transform ${\bm{\xi}_{C}}$ into the physical space by the inverse Nataf transformation ${\bm{u}_{C}}=T^{-1}(\bm{\xi}_{C})$.

• •

iii) Evaluate the accurate responses (bus voltages and line flows in this study) ${\bm{y}_{C}=(y^{(1)},y^{(2)},...,y^{({M_{C}})})}$ of ${\bm{u}_{C}}$ by the deterministic power flow solver. Pass the sample-response pairs $({\bm{\xi}_{C}},{\bm{y}_{C}})$ to Step 4.

Step 4: Apply the algorithm in Section IV-C (Algorithm 1) to solve the coefficients and the weighting factors of LRA to build the low-rank approximation (12) for all desired responses $\bm{Y}$. If LRA for all responses have reached the prescribed accuracy, go to Step 6; otherwise, go to Step 5.

Step 5: Generate additional ${\Delta M_{C}}$ new samples and evaluate them, then go back to Step 4 using the enriched experiment design $({\bm{\xi}_{C}}+{\Delta\bm{\xi}_{C}},{\bm{y}_{C}}+{\Delta\bm{y}_{C}})$.

Step 6: Calculate the mean and variance of all responses through (21) and (22), respectively. If other statistical properties (e.g., higher-order moments, PDF, CDF, etc.) are needed, go to Step 7; otherwise, go to Step 8.

Step 7: Sample $\bm{\xi}$ extensively, e.g., ${{M}_{S}}$ samples, and apply the solved functional approximation (12) to evaluate the corresponding responses $\bm{y}$ for all these samples. Compute the statistics of interests for each response.

Step 8. Generate the result report.

Remark: the number of samples ${{M}_{C}}$ in Step 3 is usually much smaller than ${{M}_{S}}$ in Step 6. Unlike MCS, LRA does not solve power flow equations for all ${{M}_{S}}$ samples in Step 6, and hence it is much more efficient. The main computational effort of LRA lies in Step 3.

VI Numerical Studies

In this section, we apply the proposed LRA method to investigate the probabilistic power flow of the modified 39-bus, 118-bus, and 1354-bus systems [45]. The LHS-based Monte Carlos simulation serves as a benchmark for validating the accuracy and the performance of the proposed method. In addition, a comparison between the LRA and the sparse PCE using the UQLab [46] is also presented. In order to compare the accuracy of them, the error indices introduced in [8] are adopted to indicate the distribution accuracy of responses as follows:

[TABLE]

where $\mu$ and $\sigma$ denote the mean and variance of the response, respectively. The star * represents the applied method, e.g., ’lra’ for the LRA method and ’pce’ for the sparse PCE method. It should be noted that the same set of samples was used for building up the LRA and the sparse PCE in all simulation cases.

In this study, we assume that the probability distributions and the associated parameters of all random inputs are available from up-front modeling. Particularly, the wind speed follows Weibull distribution, the parameters of which are $k=9.0$ and $c=2.15$; the solar radiation follows Beta distribution, the parameters of which are $\alpha=0.9$, $\beta=0.9$, $r=0$ and $s=1000$; the load power follows Normal distribution. For each individual load, the mean value is set to be its base case value and the variance is equal to 5% of its mean value. The univariate polynomial basis used in (12) for Weibull, Beta and Normal distributions are chosen to be the numerical polynomial, the Jacobi polynomial and the Hermite polynomial, respectively, according to Table I. It should be noted that the proposed method is not limited to any specific probability distribution. For simplicity, all wind generators share the same set of parameters $v_{rate}=15.0\mbox{ m/s}$, $v_{in}=4.0\mbox{ m/s}$, and $v_{out}=25.0\mbox{ m/s}$. Similarly, all solar PVs have the same values r_{c}=150.0\mbox{ W/m^{2}}, and r_{std}=1000.0\mbox{ W/m^{2}}, respectively. The linear correlation coefficient ${\rho_{ij}}$ between component $i$ and $j$ of wind speed $v$, solar radiations $r$ and load power $P_{L}$ are 0.5053, 0.8040, 0.4000, respectively.

VI-A The Modified 39-Bus System

The 39-bus system is a simplified New-England network which contains 10 generators, 21 loads and 47 branches. The total load of this network is 6254.23 MW and 1387.10 Mvar. In order to assess the impact of uncertainties on the feasibility of power flow solution of this network, 4 solar PV power plants of 120 MW are connected to bus $\{36,37,38,39\}$ and 4 wind farms of 180 MW are connected to bus $\{32,33,34,35\}$. There are totally 29 random inputs including 21 random loads. The renewable capacity penetration level is 17.44%.

We first apply the power flow tool to the deterministic system (i.e., without uncertainty), and the power flow solution is shown to be feasible, i.e., there is no voltage violation or thermal violation in the network. However, the system is prone to thermal violation, since we have scaled up the load power at bus $\{1,3,4,15,16,18,20,21,23,24,25,26,27,28,29\}$ by 10%. In other words, we consider a heavy loading condition, under which the randomness may affect the static security level of the system as shown subsequently.

Next, we exploit the proposed LRA method to assess the probabilistic characteristics (e.g., mean, variance, PDF and CDF) of these branch flows and compare the results with those of the LHS-based MCS and with those of the sparse PCE method. Applying the proposed algorithm, 146 simulations are required in Step 4-5 to build up the LRAs (12) of the desired responses (i.e., the bus voltages, branch flows and generator reactive power), which consist of 1 rank-one functions with the highest polynomial degrees $p_{i}=2$. Once the coefficients and the weighting factors are computed, the mean and the standard deviation of the responses are computed in Step 6 and are compared with those of the LHS-based MCS and with those of the sparse PCE method as shown in Table II and III. Furthermore, 5000 samples are generated in Step 7 to assess the PDF and the CDF of the responses. Fig. 2 shows the PDF and CDF of the branch flow $S_{13-14}$ computed by the the LHS-based MCS, the sparse PCE and the solved LRA, respectively. These results clearly demonstrate that the LRA method can provide accurate estimation for the probabilistic characteristics of the PPF solutions. Particularly, the accuracy of the LRA is comparable to that of the sparse PCE method.

Once we have the statistics of the quantities of interest, say the bus voltages and branch flows, one natural step afterwards is to see how the uncertainty actually affect the operating point of a power system. As shown by the CDF curve in Fig. 2, the branch flow $S_{13-14}$ in the deterministic system (454MW, dot line) is far away from the thermal limit (600 MW, pink line). However, the probability corresponding to the thermal limit 600 MW (the pink line) is 92%, indicating that there is still 8% probability that the corresponding power flow solution will exceed the thermal limit when the uncertainties are considered.

In addition, the convergence rate of the LRA method is traced by increasing the sample size from $0.25n$ to $10n$ ($n$ is the number of random inputs) to verify the robustness of the proposed method. For each sample size, 100 replications are generated and are evaluated by the proposed LRA method. As shown in Fig. 3, the statistics of the branch flow settle down to the values computed by the LHS-based MCS using 5000 samples when the sample size reaches $5n$ (i.e., 146 in this case), after which little improvement can be achieved by increasing the sample size. Similar results have been observed in other PPF solutions, indicating a nice property of the LRA that its simulation time (the required sample evaluations) grows linearly as the number of the random inputs increases.

VI-B The Modified IEEE 118-Bus System

The IEEE-118 bus system is a simplified representation of the Midwest U.S. transmission system in 1962, which contains 19 generators, 35 synchronous condensers, 177 transmission lines, 9 transformers and 91 loads [47]. Six wind farms, each with $100$ MW, are connected to bus $\{10,25,26,49,65,66\}$, and six solar PV parks, each with an installed capacity of $60$ MW, are connected to bus $\{12,59,61,80,89,100\}$. Besides, there are 99 stochastic loads. Therefore, the total number of random inputs is 111. The renewable capacity penetration level is 21.96%.

Likewise, we first apply the power flow tool to the deterministic system. It is shown that the power flow solution is feasible. However, this system is prone to voltage violation since 160% of load powers at bus $\{19,20,21,43,44,45,50,51,52\}$ are assumed, i.e., a heavy loading condition is considered to investigate the impacts of randomness on the feasibility of power flow solutions.

Next, we exploit the proposed LRA method to assess the probabilistic characteristics of the bus voltages and the branch flows and compare the results with those of the LHS-based MCS and with those of the sparse PCE method. Applying the proposed algorithm, 441 simulations are needed in Step 4-5 to build up the LRAs (12) of the PPF solutions. Moreover, the mean and standard deviation of the responses computed in Step 6 are compared with those of the LHS-MCS and with those of the sparse PCE method as shown in Table IV and V. Particularly, the PDF and the CDF of the voltage magnitude at Bus 53 assessed in Step 7 are compared with those of the MCS and the sparse PCE method using 10000 samples as presented in Fig. 4. All these results and comparisons clearly demonstrate that the LRA method can provide accurate estimation for the probabilistic characteristics of the PPF solutions. Particularly, the accuracy of the LRA is comparable to that of the sparse PCE method.

In terms of efficiency, the LHS-based MCS needs to run 10000 simulations (i.e. solving (9)), while the LRA only requires 441 simulations for solving the coefficients and the weighting factors of (12). Table VI lists a break-down time consumption evaluating the six responses listed in Table IV: the time for the experimental design ${t_{ed}}$, for solving the coefficients and the weighting factors ${t_{sc}}$, for evaluating the statistic samples ${t_{es}}$, and the total time ${t_{total}}$. The LRA method is approximately 8 times faster than the LHS-based MCS, and 1.5 times faster than the sparse PCE method. It is worth noting that the ${t_{sc}}$ and ${t_{es}}$ for the PCE and the LRA will increase linearly with the number of desired responses.

Similar to the previous example, it can be observed that the uncertainty may affect the dispatching/planning scheme of the system. For instance, $|V_{53}|$ without considering the uncertainty is 0.9412 p.u. which is within the feasible voltage range $[0.94,1.06]$. However, the probability corresponding to 0.94 p.u. at the CDF curve (Fig. 4) is 0.28, indicating that there is 28% probability that the voltage at Bus 53 will be below the lower limit. Such high probability of voltage violation implies that a feasible dispatch/planning scheme for the deterministic system may become unfeasible when the randomness of RES and loads are incorporated.

VI-C The Modified European 1354-Bus System

The 1354-bus system is a part of the European high voltage transmission network which contains 1354 buses, 260 generators, and 1991 branches. The total load of this network is 73,060 MW and 13,401 Mvar [48]. To assess the impact of uncertainties on the feasible power flow solution of this network, 20 wind generators, each with 100MW, and 20 solar PVs, each with 80MW, are connected to the system. All parameters of the test system are available online: https://github.com/shenghao/LRA-PPF. There are 713 random inputs in total including the random loads.

This system is prone to thermal violations, hence we study the probabilistic characteristics of the line flows of the set of branches that are most overloaded. The estimated mean and standard deviation of the responses computed by the proposed LRA are compared with those of the LHS-based MCS as shown in Table VII. It can be seen that the LRA method can provide accurate estimations for the statistics of the PPF solutions.

Similar results can be found for the estimation of the PDF and the CDF of the PPF solutions as presented in Fig. 5. The proposed LRA method can provide accurate PDF and CDF approximation for the line flow at branch 1758-1923. Likewise, the deterministic solution renders 40% probability of violating the thermal limit in the stochastic system due to the uncertainties, which evidently demonstrate the necessity of applying a probabilistic framework for power flow study.

In power system planning, the 10% and 90% confidence levels that the line flow will not exceed are essential to planning engineers, since the corresponding line flow values roughly reflect the desired capacity of the transmission path [6]. As an illustration, the MVA values at 10% (e.g., $S_{mcs}^{10}$) and 90% (e.g., $S_{mcs}^{90}$) confidence levels of the lines in Table VII are computed using the MCS and the LRA. These results have shown that the proposed LRA method can provide planning engineers with an accurate estimate of the two indices, as shown in Table VIII.

In terms of simulation effort, the LHS-base MCS requires 50000 simulations to obtain a converged result, while the proposed LRA method requires only 3566 simulations to get comparable accuracy. It should be noted that the sparse PCE is unable to solve this problem (with 713 inputs) due to the large size of the expanded multivariate polynomials. The issue of ‘out of memory’ is met on the computer (DELL OptiPlex 7050 with Intel Core i7-7700 (3.6GHz), 16GB RAM).

To make a more detailed comparison between the sparse PCE and the LRA regarding the capability of handling high-dimensional problems, six reduced systems with an increasing number of random inputs (i.e., 100, 200, 300, 400, 500, 600) are created based on the 1354-bus system by fixing small loads to their mean values. Again, we apply the MCS, the sparse PCE and the LRA to evaluate their PPF solutions.

Table IX shows the number of unknowns to be solved in the sparse PCE and the LRA. Even with the sophisticated hyperbolic truncation scheme [49], the number of unknown coefficients to be sent into the least angle regression algorithm [50] (column 4, $N_{trunc}^{pce}$) grows dramatically as the number of inputs increases. Nevertheless, the unknown coefficients plus weighting factors in LRA (column 8, $N_{solved}^{lra}$) grows only linearly. When the number of inputs reach 600, the spare PCE failed due to ‘out of memory’. These results demonstrate the advantage of the LRA method over the PCE method on handling high-dimensional problems when the number of random inputs is enormous.

Table X and XI show the mean and standard deviation of the line flow at the branch 1758-1923 in the six reduced systems as well as the original system. Clearly, the LRA is able to provided accurate estimation for all of the seven systems.

Remark: Since independent LRA has to be built for each response of interests, the computational effort will increase linearly with the number of desired responses.

VII Conclusion and Perspectives

In this paper, we have developed a novel method to accurately and efficiently assess PPF. Particularly, the proposed method can build up a statistically-equivalent surrogate model (i.e., the low-rank approximations) for the PPF solutions through a small number of power flow evaluations. Numerical studies show that

• •

The proposed LRA method can accurately estimate the probabilistic characteristics of the bus voltages and line flows in the PPF problem.

• •

The proposed LRA method is more computationally efficient compared to the LHS-based MCS. Moreover, once the LRA (12) for the desired response is built up, the response of any new samples can be evaluated efficiently by directly substituting to (12) instead of solving the original PPF equations (9).

• •

The proposed LRA method has higher capability of dealing high-dimensional problems compared to the PCE method. This merit stems from exploiting the retained tensor-product form and employing a sequential correction-updating scheme.

It has been revealed in our study that a feasible power flow may become unfeasible with a relatively high probability when the uncertainties are incorporated. The probability may be even larger if the penetration of RES increases thus need to be carefully investigated. In the future, we plan to develop a grouping scheme to accelerate the LRA calculation by reducing the iterations required in the correction step.