Analytic regularity and stochastic collocation of high dimensional Newton iterates

TL;DR

This paper develops a complex analytic regularity framework for high-dimensional stochastic Newton iterates, enabling efficient computation of stochastic moments using sparse grids, with applications to power flow problems.

Contribution

It introduces a new regularity theory for stochastic Newton methods and demonstrates their efficient computation via sparse grids, with proven convergence rates.

Findings

Convergence rates are subexponential or algebraic.

Sparse grids effectively compute low probability events.

Numerical experiments confirm theoretical predictions.

Abstract

In this paper we introduce concepts from uncertainty quantification (UQ) and numerical analysis for the efficient evaluation of stochastic high dimensional Newton iterates. In particular, we develop complex analytic regularity theory of the solution with respect to the random variables. This justifies the application of sparse grids for the computation of stochastic moments. Convergence rates are derived and are shown to be subexponential or algebraic with respect to the number of realizations of random perturbations. Due the accuracy of the method, sparse grids are well suited for computing low probability events with high confidence. We apply our method to the power flow problem. Numerical experiments on the 39 bus New England power system model with large stochastic loads are consistent with the theoretical convergence rates.

Tables

Table 1. Table 1. Sparse grid approximations formulas for TD and HC.

Approx. space	sparse grid:    $m$, $g$	polynomial space:    $\Lambda \left(w\right)$
Total	$m\left(i_n\right)=i_n$	$\{𝐩\in {\mathbb{N}}_0^N:{\sum }_n{p}_n\le w\}$
Degree (TD)	$g\left(𝐢\right)={\sum }_n\left(i_n-1\right)\le w$
Hyperbolic	$m\left(i\right)=i$	$\{𝐩\in {\mathbb{N}}_0^N:$
Cross (HC)	$g\left(𝐢\right)={\prod }_n\left(i_n\right)\le w+1$	${\prod }_n(p_n+1)\le w+1\}$

Equations

$$\mathbf{f}(\mathbf{x})=0.$$

$${\mathbf{x}}^{v+1}={\mathbf{x}}^v-\mathbf{J}({\mathbf{x}}^v{)}^{-1}\mathbf{f}({\mathbf{x}}^v).$$

$$\Vert \mathbf{J}(\mathbf{x})-\mathbf{J}(\mathbf{y})\Vert \le \lambda \Vert \mathbf{x}-\mathbf{y}\Vert$$

$$\Vert \mathbf{J}({\mathbf{x}}^0{)}^{-1}\Vert \le \varkappa ,\Vert \mathbf{J}({\mathbf{x}}^0{)}^{-1}\mathbf{f}({\mathbf{x}}^0)\Vert \le \delta ,h:=2\varkappa \lambda \delta \le 1,$$

$$L_{\mathbb{P}}^q(\Omega ):=\left\{u:\Omega \to \mathbb{R}|{\int }_{\Omega }|u(\omega ){|}^q\text{mbox}d\mathbb{P}(\omega )<\infty \right\}\text{mbox}and{L}_{\mathbb{P}}^{\infty }(\Omega ):=\left\{u:\Omega \to \mathbb{R}|\mathbb{P}\text{mbox}-\underset{\omega \in \Omega }{\mathrm{esssup}}|u(\omega )|<\infty \right\}.$$

$$\mathbb{P}(\mathbf{W}\in A):=\mathbb{P}({\mathbf{W}}^{-1}(A))={\int }_A\rho (\mathbf{q})\text{mbox}d\mathbf{q}.$$

$$\mathbb{E}[\mathbf{W}]={\int }_{\Gamma }\mathbf{q}\rho (\mathbf{q})\text{mbox}d\mathbf{q}.$$

$$\begin{array}{rl}{L}_{\rho }^q(\Gamma )& :=\left\{u:\Gamma \to \mathbb{R}|{\int }_{\Omega }|u(\mathbf{q}){|}^q\rho (\mathbf{q})\text{mbox}d\mathbf{q}<\infty \right\}\text{mbox}and{L}_{\rho }^{\infty }(\Gamma ):=\left\{u:\Gamma \to \mathbb{R}|\rho \text{mbox}-\underset{\mathbf{q}\in \Gamma }{\mathrm{esssup}}|u(\mathbf{q})|<\infty \right\}.\end{array}$$

$$\hat{\rho }(\mathbf{q})=\prod _{n=1}^N\hat{\rho }(q_n),\mathbf{q}=(q_1,\dots ,q_N)\in \Gamma ,$$

$${\mathcal{I}}_n^{m(i_n)}(u(q_n)):=\sum _{j_n=1}^{m(i_n)}u(q_{j_n}^i)l_{n,j_n}(q_n),$$

$${\Delta }_n^{m(i_n)}:={\mathcal{I}}_n^{m(i_n)}-{\mathcal{I}}_n^{m(i_n-1)}.$$

$${\mathcal{S}}_w^{m,g}[u(\mathbf{q})]=\sum _{\mathbf{i}\in {\mathbb{N}}_0^N:g(\mathbf{i})\le w}\underset{n=1}{\overset{N}{⨂}}({\Delta }_n^{m(i_n)})(u(\mathbf{q})).$$

$${\Lambda }^{m,g}(w)=\{\mathbf{p}\in {\mathbb{N}}^N,g({\mathbf{m}}^{-1}(\mathbf{p}+1))\le w\},$$

$${\mathbb{P}}_{{\Lambda }^{m,g}(w)}(\Gamma )=\text{span}\left\{\prod _{n=1}^N{q}_n^{p_n},\text{with }\mathbf{p}\in {\Lambda }^{m,g}(w)\right\}.$$

$$C^{0}(\Gamma;V):=\{u:\Gamma\rightarrow V\,\,\mbox{is continuous on $\Gamma$ and }\max_{y\in\Gamma}\|u(y)\|_{V}<\infty\}.$$

$$m(i_n)=\{\begin{array}{ll}1& \text{for }{i}_n=1\\ 2^{i_n-1}+1& \text{for }{i}_n>1\end{array}\text{ and }g(\mathbf{i})=\sum _{n=1}^N(i_n-1).$$

$$f(p_n)=\{\begin{array}{l}0,p_n=0\\ 1,p_n=1\\ \lceil {\log }_2(p_n)\rceil ,p_n\ge 2\end{array}.$$

$$q_{j_n}^{i_n}=-\cos \left(\frac{\pi (j_n-1)}{m(i_n)-1}\right),j_n=1,\dots ,m(i_n).$$

$${\mathbb{E}}_{\hat{\rho }}^{\mathbf{p}}[u]=\sum _{k=1}^{N_{\mathbf{p}}}{\omega }_{k}u({\mathbf{q}}^{(k)}),{\omega }_k=\prod _{n=1}^N{\omega }_{k_n}{\omega }_{k_n}={\int }_{{\Gamma }_n}{l}_{n,k_n}^2(q_n){\hat{\rho }}_n(q_n)\text{mbox}d{q}_n,$$

$$\mathbb{E}[{\mathcal{S}}_w^{m,g}[u(\mathbf{q})]]\approx {\mathbb{E}}_{\hat{\rho }}^{\mathbf{p}}[{\mathcal{S}}_w^{m,g}[u(\mathbf{q})]\frac{\rho }{\hat{\rho }}],$$

$$\begin{array}{rl}var[u(\mathbf{q})]& \approx \mathbb{E}[({\mathcal{S}}_w^{m,g}[u(\mathbf{q})]{)}^2]-\mathbb{E}[{\mathcal{S}}_w^{m,g}[u(\mathbf{q})]{]}^2\\ & \approx {\mathbb{E}}_{\hat{\rho }}^{\mathbf{p}}[({\mathcal{S}}_w^{m,g}[u(\mathbf{q})]{)}^2\frac{\rho }{\hat{\rho }}]-{\mathbb{E}}_{\hat{\rho }}^{\mathbf{p}}[{\mathcal{S}}_w^{m,g}[u(\mathbf{q})]\frac{\rho }{\hat{\rho }}{]}^2.\end{array}$$

$$C_{\mathrm{mix}}^k(\Gamma ;\mathbb{R})=\left\{u:\Gamma \to \mathbb{R}:\frac{{\partial }^{{\alpha }_1,\dots ,{\alpha }_N}u}{{\partial }^{{\alpha }_1}{q}_1\dots {\partial }^{{\alpha }_N}{q}_N}\in C^0(\Gamma ;\mathbb{R}),n=1,\dots ,N,{\alpha }_n\le k\right\}$$

$$\Vert u{\Vert }_{C_{\mathrm{mix}}^k(\Gamma ;\mathbb{R})}=\left\{u:\Gamma \to \mathbb{R}:\underset{\mathbf{q}\in \Gamma }{\max }\left|\frac{{\partial }^{{\alpha }_1,\dots ,{\alpha }_N}u(\mathbf{q})}{{\partial }^{{\alpha }_1}{q}_1\dots {\partial }^{{\alpha }_N}{q}_N}\right|<\infty \right\}.$$

$$\Vert u-{\mathcal{S}}_w^{m,g}[u]{\Vert }_{L^{\infty }(\Gamma )}\le C(k,N)\Vert u{\Vert }_{C_{\mathrm{mix}}^k(\Gamma )}{\eta }^{-k}(\log \eta {)}^{(k+2)(N-1)+1},$$

$$\begin{array}{rl}{\mathcal{E}}_{n,{\hat{\sigma }}_n}& =\left\{z\in \mathbb{C}\mathrm{with}\mathrm{Re}z=\frac{e^{{\delta }_n}+e^{-{\delta }_n}}{2}\cos (\theta ),\mathrm{Im}z=\frac{e^{{\delta }_n}-e^{-{\delta }_n}}{2}\sin (\theta ):\theta \in [0,2\pi ),{\hat{\sigma }}_n\ge {\delta }_n\ge 0\right\}\end{array}$$

$$\begin{array}{r}{T}_0(y)=1,T_1(y)=y,\dots ,T_{k+1}(y)=2y{T}_k(y)-T_{k-1}(y),\dots \end{array}.$$

$$u(y)={\alpha }_0+2\sum _{k=1}^{\infty }{\alpha }_k{T}_k(y),$$

$${\alpha }_k=\frac{1}{\pi }{\int }_{-1}^1\frac{u(y)T_k(y)}{1-y^2}dy$$

$$|u(y)-{\alpha }_0-2\sum _{k=1}^m{\alpha }_k{T}_k(y)|\le \frac{2M}{\zeta -1}{\zeta }^{-m}.$$

$$\Vert (I-{\mathcal{I}}_n^{m(i_n)})u{\Vert }_{L_{\rho }^{\infty }({\Gamma }_n)}\le (1+{\Lambda }_{m(i)})\underset{v\in {\mathcal{P}}_{m(i_n)-1}}{\min }\Vert u-v{\Vert }_{C^0(\Gamma ,\mathbb{R})}\le (1+{\Lambda }_{m(i)})\frac{2M(u)}{e^{{\sigma }_n}-1}{e}^{-{\sigma }_{n}m(i_n)},$$

Full text

This material is based upon work supported by the National Science Foundation under Grant No. 1736392. Research reported in this technical report was supported in part by the National Institute of General Medical Sciences (NIGMS) of the National Institutes of Health under award number 1R01GM131409-01.

1. Introduction

Newton iteration is a powerful method for solving many scientific and engineering problems, naturally arising in the context of power flow problems [8], non-linear Partial Differential Equations [15], among others.

Computational predictions often form bases for critical decisions in many areas, including engineering, financial markets, weather forecasting and disaster management. The rapid development of computer hardware has allowed simulation of more and more complex phenomena. But how reliable are these predictions? Can they be trusted? Rigorous uncertainty quantification (UQ) methods have become a basic tool for assessing the validity of such computational predictions.

Uncertainty quantification is a process in which uncertainties in a system are characterized and propagated into the calculation of a given Quantity of Interest (QoI). The characterization of uncertainties arises as a problem often in the solution of so-called inverse problems, while uncertainty propagation is often associated with the solution of forward problems. An additional goal of UQ is to allow for reduction of uncertainties in a QoI by improving the collection of new data, or by directly improving the physical understanding itself. Performing UQ for complex applications such as electric power grids is a daunting task. This is mainly due to the fact that it can be very computationally challenging. Effective UQ processes should be systematic, computationally efficient, and reliable.

One of the most widely used UQ techniques is the Monte Carlo method [11], which is robust and easy to implement. Indeed, a deep analysis or understanding of the underlying stochastic model is not required, making this an attractive approach for the practicing engineer and scientist. However, convergence rates for iterative approximation methods can be very slow. For the application of UQ to the power flow problem, due to the large numbers of generators, loads and transmission lines one potentially faces, the problem can be high dimensional, non-linear, non-Gaussian , and not feasible with current computational resources. An alternative approach is the use of tensor product methods. However, these methods suffer significantly from the curse of dimensionality, thus making them unattractive even for moderate dimensionalities.

If the regularity of a QoI is relatively high with respect to the fundamental random variables, then application of stochastic collocation with Smolyak sparse grids [27, 23, 22] is a good choice. Indeed, this method has become popular in the field of computational applied mathematics and engineering as a surrogate model of stochastic Partial Differential Equations (sPDEs) [9] where the QoI is composed of moderately large numbers of random variables. The method is easy to implement and non-intrusive, i.e., each collocation point corresponds to uncoupled deterministic problems. The stochastic collocation method can be used with non-linear dependence of the QoI on the random variables. However, such grids still suffer from the curse of dimensionality. Alternative adaptive techniques have been developed, including anisotropic sparse grids [22], dimension adaptive quadrature [12] and quasi-optimal sparse grids [7, 21]. Yet these methods are still not feasible for very high dimensional problems and/or low regularity of the QoI. In addition, although quasi-optimal sparse grids lead to exponential convergence rates with respect to numbers of realizations, there is to our knowledge still no systematic way to construct them.

We have particular interest in the application of the Newton iteration to the solution of the power flow equations of electric grids [8]. In practice many of the generators (wind, solar, etc) and loads are stochastic in nature, and thus a traditional deterministic power flow analysis is insufficient. UQ applied to mathematical/statistical modeling of electrical power grids is still in its infancy. A major 2016 National Academy of Sciences report underscores the importance of this new area of research in its potential to contribute to this and the next generation of electric power grids [20]. The incorporation of uncertainty in the grid has gathered interest in the power system community.

In [14] Hockenberry et al. proposed the probabilistic collocation method (PCM) with applications to power systems. Although the results of this work are good, the uncertainty is computed only with respect to a single parameter. In [25] the authors test a polynomial chaos collocation and Galerkin approach to study the uncertainty of power flow on a small 2-bus power system. The results are good for low stochastic dimensions, but it is not clear how this approach will scale for large electric power grids, with higher associated dimensions. Moreover, there is little mathematical theory on the effectiveness of this approach.

More recently, in [16], Tang et al. proposed a dimension-adaptive sparse grid method [12] by using the off-the-shelf Matlab Sparse Grid Toolbox [18, 17]. The numerical results show the feasibility of this approach. However, there are also some weaknesses. The authors did not analyze the regularity of the power flow with respect to the uncertainty, so that the rate of convergence of the sparse grid is not known. Reduced regularity of the stochastic power flow can choke the accuracy. In contrast to Monte Carlo methods, which are robust, sparse grid methods are sensitive to the regularity of the function at hand. Lack of regularity will lead to erroneous results.

This motivates the application of numerical analysis and UQ theory to the Newton iteration for such problems as the electric power grid. Many of the ideas of this paper originate from the numerical solution of stochastic PDEs [9, 23, 22], where UQ is having a large impact. The goals of our paper is to integrate these methods with the theory and practice of power systems, where we believe that UQ will eventually have a strong impact. Furthermore, from the numerical analysis perspective regularity can be determined by complex analytic extensions of the functions of interest [9, 3, 29]. These lead to sharper convergence rates than regularity in terms of derivatives.

In Section 2 the mathematical background for this paper is introduced. In particular, the Newton-Kantorovich Theorem, sparse grids and convergence rates are discussed. Furthermore, complex analysis from the numerical analysis perspective for polynomial approximation is also treated. In Section 3 the complex analytic regularity of the solution of the Newton iteration is developed with respect to the random perturbations. In Section 4 the theory developed in Section 3 justifies the application of sparse grids to the power flow equations. The sparse grids are applied to the power flow equations of the 39 bus, 10 Generator, New England model. Subexponential or algebraic convergence rates are obtained than coincide with the sparse grid convergence rates. These convergence rates make the sparse grid method suitable for computing stochastic moments to high accuracy. Furthermore, they will also be well suited for computing small event tail probabilities with high accuracy.

2. Mathematical background

In this section we introduce the general notation and mathematical background that will be used in this paper. We provide a summary of the three important topics: i) Newton-Kantorovich Theorem, ii) Stochastic spaces and iii) Sparse grids (approximation theory).

2.1. Newton-Kantorovich Theorem

Consider the Fréchet differentiable operator $\mathbf{f}:X\rightarrow Y$ that maps a convex open set $D$ of a Banach space $X$ into a Banach space $Y$. Suppose that we are interested in finding an element $\mathbf{x}\in D$ such that

[TABLE]

Assuming that a solution of equation (1) exists, a series of successive approximations $\mathbf{x}^{v}\in D$, where $v\in\mathbb{N}_{0}:=\mathbb{N}\cup{0}$, can be built. Consider the space of bounded linear operators $L(Y,X)$ from $Y$ into $X$. Let $\mathbf{J}(\mathbf{x}^{v})$ be the Fréchet derivative of $\mathbf{f}(\mathbf{x}^{v})$. Suppose that $\mathbf{J}(\mathbf{x}^{v})^{-1}\in L(Y,X)$ and consider the sequence

[TABLE]

Assumption 1**.**

We assume that that $\mathbf{J}$ satisfies the Lipschitz condition

[TABLE]

for some constant $\lambda\geq 0$, and all $\mathbf{x},\mathbf{y}\in D$. Furthermore assume that $\mathbf{x}^{0}\in D$ and there exist positive constants $\varkappa$ and $\delta$ such that

[TABLE]

$U(\mathbf{x}^{0},t^{*})\subset D$, with $U(\mathbf{x},r)$ the open ball $\{\mathbf{y}:\|\mathbf{y}-\mathbf{x}\|\leq r\}$ and $t^{*}=\frac{2}{h}(1-\sqrt{1-h})\delta$.

If Assumption 1 above is satisfied, then by the Newton-Kantorovich theorem [2], for all $v\in\mathbb{N}_{0}$,

1. (1)

The Newton iterates $\mathbf{x}^{v+1}=\mathbf{x}^{v}-\mathbf{J}(\mathbf{x}^{v})^{-1}\mathbf{f}(\mathbf{x}^{v})$ and $\mathbf{J}(\mathbf{x}^{v})^{-1}$ exist. 2. (2)

$\mathbf{x}^{v}\in U(\mathbf{x}_{0},t^{*})\subset D$. 3. (3)

$\mathbf{x}^{*}=\lim_{v\rightarrow\infty}\mathbf{x}^{v}$ exists, $\mathbf{x}^{*}\in\overline{U(\mathbf{x}^{0},t^{*})}$, and $\mathbf{f}(\mathbf{x}^{*})=\mathbf{0}$ uniquely.

2.2. Stochastic spaces

Let $\Omega$ be the set of outcomes from the complete probability space $(\Omega,\mathcal{F},\mathbb{P})$, where $\mathcal{F}$ is a sigma algebra of events and $\mathbb{P}$ is a probability measure. Define $L^{q}_{\mathbb{P}}(\Omega)$, $q\in[1,\infty]$, as the Banach spaces

[TABLE]

Let $\mathbf{W}:=[W_{1},\dots,W_{N}]$ be an $N$-component random vector measurable in $(\Omega,\mathcal{F},\mathbb{P})$ that takes values on $\Gamma:=\Gamma_{1}\times\dots\times\Gamma_{N}\subset\mathbb{R}^{N}$, with $\Gamma_{n}:=[-1,1]$. Let ${\mathcal{B}}(\Gamma)$ be the Borel $\sigma$-algebra. Define the induced measure $\mu_{\mathbf{W}}$ on $(\Gamma,{\mathcal{B}}(\Gamma))$ as $\mu_{\mathbf{W}}(A):=\mathbb{P}(\mathbf{W}^{-1}(A))$ for all $A\in{\mathcal{B}}(\Gamma)$. Assuming that the induced measure is absolutely continuous with respect to Lebesgue measure on $\Gamma$, there exists a density function $\rho(\mathbf{q}):\Gamma\rightarrow[0,+\infty)$ such that for any event $A\in{\mathcal{B}}(\Gamma)$

[TABLE]

Now, for any measurable function $\mathbf{W}\in[L^{1}_{\mathbb{P}}(\Gamma)]^{N}$ define the expected value as

[TABLE]

Define also the Banach spaces (for $q\geq 1$)

[TABLE]

Note that the above $\rho$-essential supremium is with respect to the measure induced by the density function $\rho(\cdot)$, rather than Lebesgue measure itself.

In general the density $\rho(\cdot)$ will not factorize into independent probability density functions, making higher dimensional manipulations difficult in some cases. In [3] the authors recommend use of an auxiliary probability density function $\hat{\rho}:\Gamma\rightarrow\mathbb{R}^{+}$ that factorizes into $N$ independent ones, i.e.,

[TABLE]

where we will assume that $\left\|\frac{\rho(\mathbf{q})}{\hat{\rho}(\mathbf{q})}\right\|_{L^{\infty}(\Gamma)}<\infty$. Note that in contrast $L^{\infty}(\Gamma)$ (without a subscript) is with respect to the supremium in the Lebesgue measure.

2.3. Sparse grids

Our objective is to efficiently approximate a function $u:\Gamma\rightarrow\mathbb{R}$ defined on high dimensional domains using global polynomials. The accuracy of the approach will directly depend on the regularity of the function. Let ${{\mathcal{P}}_{\mathbf{p}}(\Gamma)}\subset L^{2}(\Gamma)$ be the span of tensor product polynomials of degree at most ${\mathbf{p}=(p_{1},\ldots,p_{N})}$; i.e., ${\mathcal{P}}_{\mathbf{p}}(\Gamma)=\bigotimes_{n=1}^{N}\;{\mathcal{P}}_{p_{n}}(\Gamma_{n})$ with ${\mathcal{P}}_{p_{n}}(\Gamma_{n}):=\text{\rm span}(q_{n}^{m},\,m=0,\dots,p_{n}),$ $\quad n=1,\dots,N$. For univariate polynomial approximation, we define a sequence of levels $i=0,1,2,3,\ldots$ corresponding to increasing degrees $m(i)\in\mathbb{N}_{0}$ polynomial approximation for given coordinates. Here for a given approximation scheme, $m(\cdot)$ is a fixed function. In general our multivariate approximation scheme will assume different levels of approximation $i_{n}$ for different coordinates $n=1,\ldots,N$.

We consider separate univariate Lagrange interpolants in $\Gamma$ along each dimension $n$, given as ${\mathcal{I}}^{m(i_{n})}_{n}:C^{0}(\Gamma_{n})\rightarrow{\mathcal{P}}_{m(i_{n})-1}(\Gamma_{n})$. Specifically, let

[TABLE]

where $\{l_{n,j}\}_{j=1}^{m(i_{n})}$ is a Lagrange basis for the space ${\mathcal{P}}_{p_{n}}(\Gamma_{n})$, the set $\{q^{i}_{j_{n}}\}_{j_{n}=1}^{m(i_{n})}$, represents $m(i_{n})$ discrete locations (the interpolation knots) in $\Gamma_{n}$, the index $i_{n}\geq 0$ is the level of approximation, and $m(i_{n})\in\mathbb{N}_{+}$ is the number of collocation nodes at level $i_{n}\in\mathbb{N}_{+}$ where $m(0)=0$, $m(1)=1$ and $m(i_{n})\leq m(i_{n}+1)$ if $i_{n}\geq 1$.

Remark 1*.*

Since $m(i_{n})$ represents the number of interpolation knots for the Lagrange basis, we have $p_{n}=m(i_{n})+1$.

One of the most common approaches to constructing Lagrange interpolants in high dimensions is the formation of tensor products of ${\mathcal{I}}^{m(i_{n})}_{n}$ along each dimension $n$. However, as $N$ increases the dimension of ${\mathcal{P}}_{p}$ increases as $\prod_{n=1}^{N}(p_{n}+1)$. Thus even for moderate dimensions $N$ the computational cost of a Lagrange approximation becomes intractable. However, in the case of sufficient complex analytic regularity of the function $u$ with respect to the random variables defined on $\Gamma$, the application of Smolyak sparse grids is better suited. In the rest of this section the construction of the classical Smolyak sparse grid (see e.g. [27, 6]) is summarized. More details can be found in [4].

Consider the difference operator along the $n^{th}$ dimension given by

[TABLE]

Given an integer $w\geq 0$, called the approximation level, and a multi-index $\mathbf{i}=(i_{1},\ldots,i_{N})$ $\in{\mathbb{N}}^{N}_{0}$, let $g:{\mathbb{N}}^{N}_{0}\rightarrow{\mathbb{N}}$ be a strictly increasing function in each argument and define a sparse grid approximation of function $u(\mathbf{q})\in C^{0}(\Gamma)$, restricted in order by $g$:

[TABLE]

From the previous expression, the sparse grid approximation is obtained as a linear combination of full tensor product interpolations. However, the constraint $g(\mathbf{i})\leq w$ in (7) restricts the use of tensor grids of high degree.

Let $\mathbf{m}(\mathbf{i})=(m(i_{1}),\ldots,m(i_{N}))$ and consider the set of polynomial multi-degrees

[TABLE]

where $\boldsymbol{1}$ is an $N$ dimensional vector of ones. Denote by $\mathbb{P}_{\Lambda^{m,g}(w)}(\Gamma)$ the corresponding multivariate polynomial space spanned by the monomials with multi-degree in $\Lambda^{m,g}(w)$, i.e.

[TABLE]

For a Banach space $V$ let

[TABLE]

It can be shown that the approximation formula given by $\mathcal{S}^{m,g}_{w}$ is exact in $\mathbb{P}_{\Lambda^{m,g}(w)}(\Gamma)$. We state the following proposition that is proved in [4].

Proposition 1**.**

• a)

For any $u\in C^{0}(\Gamma;V)$, we have $\mathcal{S}_{w}^{m,g}[u]\in\mathbb{P}_{\Lambda^{m,g}(w)}\otimes V$.

• b)

Moreover, $\mathcal{S}_{w}^{m,g}[u]=u\;\;\forall u\in\mathbb{P}_{\Lambda^{m,g}(w)}\otimes V$.

Remark 2*.*

The tensor product space $\mathbb{P}_{\Lambda^{m,g}(w)}\otimes V$ is more easily understood as the space of polynomials with Banach-valued coefficients. Furthermore, $\mathcal{S}_{w}^{m,g}[u]\in\mathbb{P}_{\Lambda^{m,g}(w)}\otimes V$ is interpreted as a sparse grid approximation of a $V$-valued continuous function.

A good choice of $m$ and $g$ is given by the Smolyak sparse grid definitions (see [27, 6])

[TABLE]

Furthermore $\Lambda^{m,g}(w):=\{\mathbf{p}\in{\mathbb{N}}^{N}:\;\;\sum_{n}f(p_{n})\leq w\}$ where

[TABLE]

Other common choices are shown in Table 1.

This choice of $m$ and $g$ combined with the choice of Clenshaw-Curtis abscissas (which are locations of interpolation points given as extrema of Chebyshev polynomials) leads to nested sequences of one dimensional interpolation formulas and a sparse grid with a highly reduced number of nodes compared to the corresponding tensor grid. Another good choice includes Gaussian abscissas [24]. For any choice of $m(i_{n})>1$ the Clenshaw-Curtis abscissas are given by

[TABLE]

In Figure 1 an example of Clenshaw Curtis and Gaussian abscissas are shown for $w=5$.

As previously pointed out, the probability density function $\rho$ does not necessarily factorize in higher dimensions. As an alternative we use the auxiliary distribution $\hat{\rho}$, which factorizes as $\hat{\rho}(\mathbf{q})=\prod_{n=1}^{N}\hat{\rho}_{n}(q_{n})$ and is close to the original distribution $\rho(\mathbf{q})$. Suppose $k$ is a given global index determined by the set of indices $k_{1}\dots,k_{N}$ as $k=k_{1}+p_{1}(k_{2}-1)+p_{1}p_{2}(k_{3}-1)+p_{1}p_{2}p_{3}(k_{4}-1)+\dots$. Given a function $u:\Gamma\rightarrow V$, the quadrature scheme $\mathbb{E}_{\hat{\rho}}^{\mathbf{p}}[u]$ that approximates the integral $\mathbb{E}[u(\mathbf{q})]:=\int_{\Gamma}u(\mathbf{q})\hat{\rho}(\mathbf{q})\,d\mathbf{q}$ can now be computed based on the distribution $\hat{\rho}(\mathbf{q})$ as

[TABLE]

with $\mathbf{q}^{(k)}\in\mathbb{R}^{n}$ the locations of the quadrature knots, and $N_{\mathbf{p}}\in\mathbb{N}$ the number of Gauss quadrature points controlling the accuracy of the quadrature scheme. Recall that $l_{n,k_{n}}^{2}(q_{n})$ define Lagrange polynomials defined in equation (5). The term $\mathbb{E}[u(\mathbf{q})]$ can be approximated as

[TABLE]

and similarly the variance $\hbox{var}[u(\mathbf{q})]$ is approximated as

[TABLE]

Remark 3*.*

The weights $\omega_{k_{n}}$ and node locations $\mathbf{q}^{(k)}$ are computed from the auxiliary density $\hat{\rho}$. For standard distributions of $\hat{\rho}$ such as uniform and Gaussian, these are already tabulated to full accuracy. Otherwise they must be computed by solving for roots of orthogonal polynomials and using a quadrature scheme. However, the integrals involved are only one dimensional. See [3] (Section 2) for details.

We now develop some rigorous numerical bounds for the accuracy of the sparse grid approximation. Let $C^{k}_{\rm mix}(\Gamma;\mathbb{R})$ denote the space of functions with continuous mixed derivatives up to degree $k$:

[TABLE]

and equipped with the following norm:

[TABLE]

Assume that $u\in C^{k}_{\rm mix}(\Gamma;\mathbb{R})$. In [6] the authors show that it must follow that

[TABLE]

where $\eta$ is the number of knots of the sparse grid $\mathcal{S}^{m,g}_{w}$. However, the coefficient $C(k,N)$ is in general not known [6].

If the function $u$ admits a complex analytic extension, a better approach for deriving error bounds for the polynomial approximation arises from exploitation of analysis in the complex plane. In [23] the authors derive $L^{\infty}_{\rho}(\Gamma)$ bounds based on analytic extensions of $u$ on a well defined region $\Psi\subset\mathbb{C}^{N}$ with respect to the variables $\mathbf{q}$. These bounds are explicit, and the coefficients can be estimated and depend on the size of the region $\Psi$.

In [22, 23] the authors derive error estimates for isotropic and anisotropic Smolyak sparse grids with Clenshaw-Curtis and Gaussian abscissas, with error $\|u-\mathcal{S}^{m,g}_{w}[u]\|_{L^{\infty}_{\rho}(\Gamma)}$ exhibiting algebraic or sub-exponential convergence with respect to the number of collocation knots $\eta$ (see Theorems 3.10, 3.11, 3.18 and 3.19 in [23] for more details). However, for these estimates to be valid, $u\in C^{0}(\Gamma,\mathbb{R})$ has to admit an analytic extension to the region defined by the following polyellipse in $\mathbb{C}^{N}$: ${\mathcal{E}}_{\hat{\sigma}_{1},\dots,\hat{\sigma}_{N}}:=\Pi_{n=1}^{N}{\mathcal{E}}_{n,\hat{\sigma}_{n}}\subset\Psi$, where

[TABLE]

and $\hat{\sigma}_{n}>0$ (see the Bernstein ellipse in Figure 2).

We now recall the definition of Chebyshev polynomials, useful in deriving error estimates on sparse grids. Let $T_{k}:\Gamma_{1}\rightarrow\mathbb{R}$, $k=0,1,\dots$, be a $k^{th}$ order Chebyshev polynomial over $[-1,1]$. These polynomials are defined recursively as:

[TABLE]

The following theorem characterizes approximation of analytic functions using Chebyshev polynomials.

Theorem 1**.**

Let $u$ be analytic and absolutely bounded by $M$ on ${\mathcal{E}}_{\log\zeta}$, $\zeta>1$. Then the expansion

[TABLE]

holds for all $y\in{\mathcal{E}}_{\log\zeta}$ where

[TABLE]

and, additionally, $|\alpha_{k}|\leq M/\zeta^{k}$. Furthermore if $y\in[-1,1]$ then

[TABLE]

Proof.

See Theorem 8.2 in [29] ∎

We follow the arguments in [24, 3] and using the fact that the interpolation operator ${\mathcal{I}}^{m(i_{n})}_{n}$ is exact on the space ${\mathcal{P}}_{p_{n}-1}$, i.e. for any $v\in{\mathcal{P}}_{p_{n}-1}$ we have that ${\mathcal{I}}^{m(i_{n})}_{n}(v)=v$, it can be shown that if $u$ is continuous on $[-1,1]$ and has an analytic extension on ${\mathcal{E}}_{\sigma_{n}}$ we have, from Theorem 1,

[TABLE]

where $\Lambda_{m(i_{n})}$ is the Lebesgue constant and is bounded by $2\pi^{-1}(\log{(m-1)}+1)$ (see [3]) and $M=M(u)$ is the maximal value of $u$ on ${\mathcal{E}}_{\sigma_{n}}$ . Thus, for $n=1,\dots,N$ we have

[TABLE]

where $\sigma_{n}=\frac{\hat{\sigma}_{n}}{2}>0$ and $C(\sigma_{n}):=\frac{2}{(e^{\sigma_{n}}-1)}$. Recalling the definition of ${\Delta_{n}^{m(i_{n})}}$ from equation (6), we have that for all $n=1,\dots,N$

[TABLE]

By applying equation (9) to Lemma 3.5 in [23], we are now in a position to slightly modify Theorems 3.10 and 3.11 in [23] and restate them into a single theorem given below. However, the following assumptions and definitions are first needed:

• •

We set $\hat{\sigma}\equiv\min_{n=1,\dots,N}\hat{\sigma}_{n},$ i.e. for an isotropic sparse grid the overall asymptotic sub-exponential decay rate $\hat{\sigma}$ will be dominated by the smallest $\hat{\sigma}_{n}$.

• •

Let

[TABLE]

$\sigma=\hat{\sigma}/2$, $\mu_{1}=\frac{\sigma}{1+\log(2N)}$, and $\mu_{2}(N)=\frac{\log(2)}{N(1+\log(2N))}$,

[TABLE]

[TABLE]

$\mu_{3}=\frac{\sigma\delta^{*}\tilde{C}_{2}(\sigma)}{1+2\log(2N)}$, and

[TABLE]

Theorem 2**.**

Suppose that $u\in C^{0}(\Gamma;\mathbb{R})$ has an analytic extension on ${\mathcal{E}}_{\hat{\sigma}_{1},\dots,\hat{\sigma}_{N}}$ and is absolutely bounded by $\tilde{M}(u)$. If $w>N/\log{2}$ and a sparse grid with Clenshaw-Curtis abscissas is used, then the following bound is valid:

[TABLE]

Furthermore, if $w\leq N/\log{2}$ then the following algebraic convergence bound holds:

[TABLE]

Proof.

This is proved by applying the inequality of equation (9) to the proof of Theorems 3.10 and 3.11 in [23]. ∎

In many practical cases not all dimensions of $\Gamma$ are equally important. In these cases the dimensionality of the sparse grid can be significantly reduced by means of anisotropic sparse grids. It is straightforward to build related anisotropic sparse approximation formulas by having the function $g$ act differently on different input random variables $q_{n}$. Anisotropic sparse stochastic collocation [22] combines the advantages of isotropic sparse collocation with those of anisotropic full tensor product collocation.

3. Analyticity of the Newton iteration

It is profitable here for purposes of clarification to consider to consider the Newton iteration in a general function space context. Specifically, let $X$ and $Y$ (see section 2) be Banach spaces. Consider the following problem: Find $\mathbf{x}\in D\subset X$ such that

[TABLE]

where $\mathbf{q}\in\Gamma$ and $\mathbf{f}:D\times\Gamma\rightarrow Y$. Equation (12) is then solved using the Newton iteration under the conditions of the Newton-Kantorovich Theorem (see section 2).

The convergence rate of the Newton iterates based on a sparse grid approximation as a function of grid size is directly affected by the regularity properties with respect to parameters $\mathbf{q}\in\Gamma$. Regularity is characterized in terms of an analytic extension in $\mathbb{C}^{N}$ of the iterates.

In the sequel we will treat $\mathbf{q}\in\Gamma$ as a random parameter, which will be suppressed occasionally. Thus for example, below we will write $\mathbf{f}\equiv\mathbf{f}_{\mathbf{q}}$, $\mathbf{J}\equiv\mathbf{J}_{\mathbf{q}}$, $M_{v}\equiv M_{v,\mathbf{q}}$, etc., and we can write $\mathbf{f}(\mathbf{x},\mathbf{q})\equiv\mathbf{f}_{\mathbf{q}}(\mathbf{x})$. For any $\mathbf{q}\in\Gamma$ (which we fix for now) consider the Newton sequence

[TABLE]

where $\mathbf{J}:X\rightarrow Y$ is the Fréchet derivative of $\mathbf{f}:D\rightarrow Y$. Assume that $\mathbf{x}_{0}\in D_{0}\subset D$ and for all $v\in\mathbb{N}$ let $D_{v}\subset D$ be the successive images under the map $M_{v}$, so that $M_{v}(D_{v-1})=D_{v}$. Note at this point the random parameter $\mathbf{q}\in\Gamma$ is unchanging throughout the iteration; the iterated domains $D_{v}$ however depend on $\mathbf{q}$.

Suppose that the parameter $\mathbf{q}$ is now extended to a complex parameter $\mathbf{g}$ with $\mathbf{g}\in\Psi\supset\Gamma$, where $\Psi\subset\mathbb{C}^{N}$. We can now form a complex extension of the sequence (13) as follows. We will complexify the pair $(\mathbf{x},\mathbf{q})$ into a pair of complex variables $(\mathbf{z},\mathbf{g})$, with $\mathbf{z}$ the complexification of $\mathbf{x}$. Assume that for $\mathbf{f}(\mathbf{x}^{0})\equiv\mathbf{f}_{\mathbf{q}}(\mathbf{x}^{0}):D_{0}\rightarrow E_{0}$ there exists an analytic extension $\mathbf{f}(\mathbf{x}^{0})\equiv\mathbf{f}_{\mathbf{g}}(\mathbf{z}^{0}):\Theta_{0}\rightarrow\Phi_{0}$, where $\Theta_{0}$ and $\Phi_{0}$ are contained in a suitable complex Banach spaces, which are the respective complex extensions of $X$ and $Y$. Given a function $\mathbf{f}$ on a real linear domain $D$, and a function $\mathbf{f}^{*}$ on a complex linear domain $\Theta\supset D$, we say that $\mathbf{f}^{*}$ is an analytic extension of $\mathbf{f}$ if $\mathbf{f}^{*}$ is analytic on its domain, and the restriction $\mathbf{f}^{*}|_{D}=\mathbf{f}$. When there exists an analytic extension $\mathbf{f}^{*}$ we say that $\mathbf{f}$ can be analytically extended.

Remark 4*.*

Note that as before we write $\mathbf{f}\equiv\mathbf{f}_{\mathbf{g}}$, $\mathbf{J}\equiv\mathbf{J}_{\mathbf{g}}$, $M_{v}\equiv M_{v,\mathbf{g}}$, etc. It is understood from context that the notational equivalence is over the extension of the variables $(\mathbf{x},\mathbf{q})$ into the complex pair $(\mathbf{z},\mathbf{g})$.

Similarly, assume that $\mathbf{J}(\mathbf{x}^{0})\in L(D_{0},E_{0})$ can be extended analytically as $\mathbf{J}(\mathbf{z}^{0})\in L(\Theta_{0},\Phi_{0})$. Here $L(\cdot,\cdot)$ is the space of bounded operators between two spaces. Through the above complexifications, equation (12) defines a complexification of the mapping $M_{v}$. We now repeat the above iteration using the complexified maps defined here. Thus there exists a series of sets $\Theta_{0},\dots,\Theta_{v}$ and $\Phi_{0},\dots,\Phi_{v}$, such that for $\mathbf{f}(\mathbf{x}^{v}):D_{v}\rightarrow E_{v}$ and $\mathbf{J}(\mathbf{x}^{v})\in L(D_{v},E_{v})$ the analytic extensions $\mathbf{f}(\mathbf{z}^{v}):\Theta_{v}\rightarrow\Phi_{v}$ and $\mathbf{J}(\mathbf{z}^{v})\in L(\Theta_{v},\Phi_{v})$ are onto. Thus the sequence (13) is extended in $\mathbb{C}^{N}$ as follows: Let $\mathbf{z}^{0}=\mathbf{x}^{0}$ and for all $v\in\mathbb{N}$ and $\mathbf{g}\in\Psi$ (which we also fix) form the sequence

[TABLE]

where $M_{v}:\Theta_{v-1}\rightarrow\Theta_{v}$.

Remark 5*.*

The domain $\Theta_{0}$ contains the initial condition $\mathbf{z}^{0}$. Under certain assumptions and with a judicious choice of $\Theta_{0}$ it can be shown that the sequence in (14) converges in a pointwise sense inside $\Theta_{0}$. This will be explored in detail in section 3.1.

Suppose that $\mathbf{z}^{v}$ is an analytic extension of $\mathbf{x}_{v}$ on $\Psi\subset\mathbb{C}^{N}$ (see Figure 3). Then the convergence rates of the sparse grid applied to any entry of interest of $\mathbf{z}^{v}$ can be characterized. The size of the set $\Psi$ determines the regularity properties of the solution. From the sparse grid discussion in section 2.3 we embed a polyellipse ${\mathcal{E}}_{\hat{\sigma}_{1},\dots,\hat{\sigma}_{N}}:=\Pi_{n=1}^{N}{\mathcal{E}}_{n,\hat{\sigma}_{n}}$ in $\Psi$. From Theorem 10 the $L^{\infty}_{\rho}(\Gamma)$ convergence rate of the sparse grid is sub-exponential (or algebraic) with respect to the number of sparse grid knots $\eta$. The decay of the sparse grid is dominated by $\sigma=\min_{n=1,\dots,N}\sigma_{n}$. Thus, the larger $\sigma$ is the faster the convergence rate.

Remark 6*.*

For finite dimensional spaces $X=Y=\mathbb{R}^{m}$, $m\in\mathbb{N}_{0}$, the Fréchet derivative $\mathbf{J}$ corresponds to the Jacobian of $\mathbf{f}$. In the rest of the paper it is assumed that $D\subset\mathbb{R}^{m}$ and $\|\cdot\|$ corresponds to the standard Euclidean norm or the standard matrix norm, depending on context. For the case of the power flow equations $m$ will be simply related to the number of nodes of the power system [8]. We will be using the notion of analytic extensions, which can be defined as follows.

We can now prove an important theorem for our purposes. First, denote $\mathbf{f}(\mathbf{z}^{v}):\Theta_{v}\rightarrow\Phi_{v}$ as $\mathbf{f}^{v}$, and $\mathbf{J}(\mathbf{z}^{v})\in L(\Theta_{v},\Phi_{v})$ as $\mathbf{J}^{v}$ (note that this depends on the generic initial $\mathbf{z}_{0}\in\Theta_{0}$ at which the Jacobian is computed).

Theorem 3**.**

Assume that for all $v\in\mathbb{N}_{0}$: $D_{v}\subset X$ and

1. (i)

$\mathbf{f}^{v}:D_{v}\rightarrow E_{v}$* can be analytically extended to $\mathbf{f}^{v}:\Theta_{v}\rightarrow\Phi_{v}$.* 2. (ii)

There exists a coefficient $c_{v}>0$ such that

[TABLE]

where $\sigma_{min}(\cdot)$ refers to the minimum singular value, $\mathbf{J}^{v}_{R}:=\mathop{\text{\rm Re}}{\mathbf{J}^{v}}$ and $\mathbf{J}^{v}_{I}:=\mathop{\text{\rm Im}}{\mathbf{J}^{v}}$.

Then for all $v\in\mathbb{N}_{0}$ there exists an analytic extension of $\mathbf{x}^{v}$ on $\Psi$.

Proof.

The main strategy for this proof is to use the Cauchy-Riemann equations. This avoids having to explicitly show that the inverse of the complex Jacobian matrix $\mathbf{J}^{v}$ is analytic. The existence of an analytic extension for $\mathbf{z}^{v}$ for each separate complex dimension is shown. The Hartog’s Theorem is then used to show analyticity with respect to all the complex dimensions.

For fixed $v$ consider the extension $\mathbf{x}^{v}\rightarrow\mathbf{z}^{v}=\mathbf{x}^{v}+\mathbf{w}^{v}$ in $\Theta_{v}$, where $\mathbf{w}^{v}\in\mathbb{C}^{m}$ and $\mathbf{z}^{v}\in\Theta_{v}$. In complex form $\mathbf{z}^{v}=\mathbf{z}^{v}_{R}+i\mathbf{z}^{v}_{I}$, where $\mathbf{z}^{v}_{R}=\mathop{\text{\rm Re}}\mathbf{z}^{v}$, and $\mathbf{z}^{v}_{I}=\mathop{\text{\rm Im}}\mathbf{z}^{v}$. Furthermore, consider the extension of $\mathbf{q}\rightarrow\mathbf{g}=\mathbf{q}+\mathbf{v}$ in $\Psi$, where $\mathbf{v}\in\mathbb{C}^{N}$ and $\mathbf{g}\in\Psi$. The extension of the iteration (13) on $\Theta_{v}\times\Psi$ leads to the following block form iteration

[TABLE]

where $\mathbf{f}^{v}_{R}:=\mathop{\text{\rm Re}}{\mathbf{f}^{v}}$ and $\mathbf{f}^{v}_{I}:=\mathop{\text{\rm Im}}{\mathbf{f}^{v}}$. From $(ii)$ it follows that equation (15) is well posed and is a valid extension of equation (13) on $\Theta_{v}\times\Psi$. We now show that $\mathbf{z}^{v+1}$ is an analytic extension on $\Theta_{v}\times\Psi$.

We focus our attention on the $k^{th}$ variable of $\mathbf{z}^{v}$ as $z^{v}_{k}$ and write it in complex form as $z^{v}_{k}=s+iw$. By differentiating equation (15) with respect to $s$ and $w$ we obtain

[TABLE]

From assumption $(ii)$ we conclude that $\partial_{s}\mathbf{z}^{v+1}_{R}$, $\partial_{s}\mathbf{z}^{v+1}_{I}$, $\partial_{w}\mathbf{z}^{v+1}_{R}$ and $\partial_{w}\mathbf{z}^{v+1}_{I}$ exist on $\Theta_{v}\times\Psi$. The following step is to show that the Cauchy-Riemann equations for $\mathbf{z}^{v+1}$ are satisfied on $\Theta_{v}\times\Psi$.

Let $P(\mathbf{z}^{v}):=\partial_{s}\mathbf{z}^{v}_{R}-\partial_{w}\mathbf{z}^{v}_{I}$ and $Q(\mathbf{z}^{v}):=\partial_{w}\mathbf{z}^{v}_{R}+\partial_{s}\mathbf{z}^{v}_{I}$, then from equation (16)

[TABLE]

Now, $(i)$ implies that $\mathbf{J}^{v}\in L(D_{v},E_{v})$ can be analytically extended to $\mathbf{J}^{v}\in L(\Theta_{v},\Phi_{v})$. Since $\mathbf{J}(\mathbf{z}^{v},\mathbf{g})$ and $\mathbf{f}^{v}(\mathbf{z}^{v},\mathbf{g})$ are analytic on $\Theta_{v}\times\Psi$ then from the Cauchy-Riemann equations

[TABLE]

Since $z^{v}_{k}$ is a linear polynomial of $s+iw$ then $P(\mathbf{z}^{v})=Q(\mathbf{z}^{v})=\mathbf{0}$ on $\mathbb{C}^{N}$ and thus $P(\mathbf{z}^{v+1})=Q(\mathbf{z}^{v+1})=\mathbf{0}$ on $\Theta_{v}\times\Psi$. We conclude that $\mathbf{z}^{v+1}$ is analytic for the $k^{th}$ variable for all $\mathbf{z}^{v}\in\Theta_{v}$ and $\mathbf{g}\in\Psi$. Following a similar argument we can show that for $l=1,\dots,N$ the $l^{th}$ variable extension of $\mathbf{q}$ has leads to an analytic extension of $\mathbf{z}^{v+1}$ whenever $\mathbf{z}^{v}\in\Theta_{v}$ and $\mathbf{g}\in\Psi$. We now extend the analyticity of $\mathbf{z}^{v+1}$ on all of $\Theta_{v}\times\Psi$.

Since $z^{v+1}_{k}$ is analytic for all $k=1,\dots,m$ and the $l^{th}$ variable of $\mathbf{q}$ has an analytic extension for all $l=1,\dots,N$ whenever $\mathbf{z}^{v}\in\Theta_{v}$ and $\mathbf{g}\in\Psi$, then from Hartog’s theorem we conclude that $\mathbf{z}^{v+1}$ is continuous on $\Theta_{v}\times\Psi$. From Osgood’s lemma it follows that $\mathbf{z}^{v+1}$ is analytic on $\Theta_{v}\times\Psi$. From an induction argument and using that fact that the composition of analytic functions is analytic then it follows that $\mathbf{z}^{v+1}$ is analytic in $\Psi$, for all $v\in\mathbb{N}$. ∎

If the assumptions of Theorem 3 are satisfied then $\mathbf{z}^{v}$ is complex analytic in $\Psi$ and it is reasonable to construct a series of sparse grid surrogate models of the entries of the vector $\mathbf{x}^{v}$. Note that in practice we restrict out attention to a subset of the variables of interest of $\mathbf{x}^{v}$. With a slight abuse of notation denote ${\mathcal{S}}^{m,g}_{w}[\mathbf{x}^{v}(\mathbf{q})]$ as the sparse grid approximation of the entries of interest of the vector $\mathbf{x}^{v}$.

From Theorem 11 we observe that the accuracy of the sparse grid approximation is a function of i) the size of the polyellipse ${\mathcal{E}}_{\hat{\sigma}_{1},\dots,\hat{\sigma}_{N}}\subset\Psi$ and ii)

[TABLE]

If the complex sequence (14) does not converge, then the size of the sets $\Theta_{v}$ can become unbounded. In particular, it is possible that $\tilde{M}(\mathbf{z}^{v})\rightarrow\infty$ as $v\rightarrow\infty$ even if $\mathbf{J}(\mathbf{z}^{v})^{-1}\in L(\Phi_{v},\Theta_{v})$ exists for all $v\in\mathbb{N}_{0}$. Our objective now is to analyze under what conditions the complex sequence remains bounded. In particular, for all $v\in\mathbb{N}_{0}$, we ask if it is possible to construct bounded regions $U\subset\mathbb{C}^{m}$ and $\Psi\subset\mathbb{C}^{N}$ such that $\mathbf{z}^{v}$ is contained in $U$ and thus

[TABLE]

To help answer this question we first show that the complex sequence (14) is itself a Newton sequence.

Remark 7*.*

We have to clarify what we mean by the Fréchet derivative of the complex function $\mathbf{f}:\Theta_{v}\rightarrow\Psi_{v}$. The algebraic problem of equation (12) can be complexified as follows: Find $\mathbf{z}\in\Theta_{0}$ such that $\mathbf{f}(\mathbf{z},\mathbf{g})=\mathbf{0}$ for all $\mathbf{g}\in\Psi$. This can be re-written in vector form as: Find $\mathbf{z}=\mathbf{z}_{R}+i\mathbf{z}_{I}\in\Theta_{0}$ such that $\mathop{\text{\rm Re}}\mathbf{f}(\mathbf{z}_{R},\mathbf{z}_{I},\mathbf{g}_{R},\mathbf{g}_{I})=\mathbf{0}$ and $\mathop{\text{\rm Im}}\mathbf{f}(\mathbf{z}_{R},\mathbf{z}_{I},\mathbf{g}_{R},\mathbf{g}_{I})=\mathbf{0}$ for all $\mathbf{g}=\mathbf{g}_{R}+i\mathbf{g}_{I}\in\Psi$. The corresponding Newton iteration is based on

[TABLE]

where $\partial_{\mathbf{z}_{R}}\mathbf{f}^{v}_{R}$ is the Fréchet derivative of $\mathbf{f}^{v}_{R}$ with respect to the variables $\mathbf{z}_{R}$ and similarly for the rest. We refer to the matrix

[TABLE]

as the Fréchet derivative of $\mathbf{f}:\Theta_{v}\rightarrow\Psi_{v}$.

Theorem 4**.**

Suppose assumptions i) and ii) of Theorem 3 are satisfied. Then the complex analytic extension $\mathbf{J}(\mathbf{z}^{v})\in L(\Theta_{v},\Phi_{v})$ of $\mathbf{J}(\mathbf{x}^{v})\in L(D_{v},E_{v})$ is equivalent to the Fréchet derivative of $\mathbf{f}^{v}:\Theta_{v}\rightarrow\Phi_{v}$, i.e.

[TABLE]

Proof.

We first prove this result for $m=1$ dimension. Suppose that $f:D\rightarrow\mathbb{R}$, is a Fréchet differentiable function and let $f:\Xi\rightarrow\mathbb{C}$ be the analytic continuation on the non-empty open set $\Xi\subset\mathbb{C}$. The analytic function $f:\Xi\rightarrow\mathbb{C}$ can be rewritten as $f(x,y)=f_{R}(x,y)+if_{I}(x,y)$ for all $x+iy\in\Xi$. Since $f$ is analytic on $\Xi$, from the the identity theorem [1] (uniqueness of complex analytic extensions) we have that $f_{R}$ and $f_{I}$ are unique in $\Xi$. Furthermore, since $f$ is analytic the Cauchy-Riemann equations are satisfied. Thus

[TABLE]

in $\Xi$ and $f:\Xi\rightarrow\mathbb{C}$ is Fréchet differentiable. Now, $\partial_{x}f(x,y)=\partial_{x}f_{R}(x,y)+i\partial_{x}f_{I}(x,y)$ in $\Xi$, and from the uniqueness property of the Fréchet derivative all the terms are unique. Recall that $\partial_{x}f(x)$ defined in $D$ is the Fréchet derivative of $f:D\rightarrow\mathbb{R}$. Write the analytic extension of $\partial_{x}f(x)$ (defined in $D$) on $\Xi$ as $g(x,y)+ih(x,y)$, with $x+iy\in\Xi$. Since $\partial_{x}f(x)=\partial_{x}f(x,y)=\partial_{x}f_{R}(x,y)+i\partial_{x}f_{I}(x,y)$ for $y=0$ and $x\in D$, from the uniqueness of the analytic extension we conclude $g=\partial_{x}f_{R}$ and $h=\partial_{x}f_{I}$ for all $x+iy\in\Xi$. From equation (18) the conclusion follows.

We can now prove our statement for the general case using a simple extension of the above argument. Since $\mathbf{f}^{v}:\Theta_{v}\rightarrow\Phi_{v}$ is complex analytic, from the identity theorem [1] it is the unique extension of $\mathbf{f}^{v}:D_{v}\rightarrow E_{v}$. (Note that the unique extension of the identity theorem applies in multi-variate case, which includes the variables $\mathbf{x}^{v}$ and $\mathbf{q}$ in the domains $\Theta_{v}$ and $\Psi$ respectively.) From the Cauchy-Riemann equations the functions $\mathbf{f}^{v}:\Theta_{v}\rightarrow\Phi_{v}$ are Fréchet differentiable and unique. Now, with a slight abuse of notation, denote $\mathbf{J}_{\mathbf{z}^{v}_{R}}$ as the Jacobian of $\mathbf{f}^{v}:\Theta_{v}\rightarrow\Phi_{v}$, with respect to the real variables $\mathbf{z}^{v}_{R}$ only. By using the above one dimensional argument we can show that the analytic extension of each entry of $\mathbf{J}(\mathbf{x}^{v})$ matches $\mathbf{J}_{\mathbf{z}^{v}}$ on the real part of $\Theta_{v}$. From the Cauchy-Riemann equations we conclude that $\mathbf{J}(\mathbf{z}^{v})\in L(\Theta_{v},\Phi_{v})$ is equivalent to the Fréchet derivative of $\mathbf{f}^{v}:\Theta_{v}\rightarrow\Phi_{v}$. ∎

From Theorem 4 it follows that the complex sequence (14) is a Newton sequence. We can now apply the Newton-Kantorovich Theorem to study the sequence convergence as $v\rightarrow\infty$.

3.1. Regions of Analyticity

The size of a polyellipse embedded in the domain $\Psi$ and the magnitude of $\mathbf{z}^{v}\in\Theta_{v}$ (for any $v\in\mathbb{N}_{0}$) directly impacts the accuracy of the sparse grid (c.f. Theorem 11). For each $v\in\mathbb{N}_{0}$ the size of the domains $\Theta_{v}$ and $\Psi$ will be characterized by the magnitude of the minimum singular value

[TABLE]

for some $c_{v}>0$. However, constructing the domain $\Psi$ would require imposing inequality conditions for each Newton iteration. This leads to a highly complex coupled problem that is hard to solve. Moreover, if the complex extension $\mathbf{z}^{v}$ grows rapidly with respect to $v$ then the size of the domain $\Psi$ will be most likely severely constrained. In contrast, by applying the Newton-Kantorovich Theorem it is sufficient to impose conditions on the initial Jacobian $(v=0)$ to construct a region of analyticity for $\Psi\subset\mathbb{C}^{N}$. Furthermore, the size of the iteration $\mathbf{z}^{v}$ will be controlled.

Consider the iteration

[TABLE]

where $\mathbf{g}\in\Psi$,

[TABLE]

for all $v\in\mathbb{N}_{0}$.

Remark 8*.*

From Theorem 4 or, alternatively, the Cauchy-Riemann equations, the matrix $\mathbf{J}(\boldsymbol{\alpha}^{v},\mathbf{g})$ corresponds to the Fréchet derivative of $\mathbf{f}(\boldsymbol{\alpha}^{v},\mathbf{g})$. Thus the sequence (20) is an Newton iteration and the Newton-Kantorovich Theorem can be used to analyze its convergence properties.

Assumption 2**.**

For all $\mathbf{q}\in\Gamma$ Assumption 1 is satisfied.

Assumption 3**.**

Assume that $\tilde{D}$, where $D\subset\tilde{D}$, is an open convex set in $\mathbb{R}^{2m}$ and the following Lipschitz condition is satisfied:

[TABLE]

for all $\mathbf{x},\mathbf{y}\in\tilde{D}$, $\mathbf{g}\in\Psi$, and $\lambda_{e}\geq 0$. Furthermore assume that for all $\mathbf{g}\in\Psi$

[TABLE]

and $U(\boldsymbol{\alpha}^{0},t^{*}_{e})\subset\tilde{D}$, where $t^{*}_{e}=\frac{2}{h_{e}}(1-\sqrt{1-h_{e}})\delta_{e}$.

Theorem 5**.**

If Assumptions 2 and 3 are satisfied then for all $\mathbf{g}\in\Psi$

1. (1)

The Newton iterates $\boldsymbol{\alpha}^{v+1}=\boldsymbol{\alpha}^{v}+\mathbf{J}(\boldsymbol{\alpha}^{v},\mathbf{g})^{-1}\mathbf{f}(\boldsymbol{\alpha}^{v},\mathbf{g})$ exist and $\boldsymbol{\alpha}^{v}\in U(\boldsymbol{\alpha}_{0},t^{*}_{e})\subset\tilde{D}$. 2. (2)

$\boldsymbol{\alpha}^{*}:=\lim_{v\rightarrow\infty}\boldsymbol{\alpha}^{v}$* exists, $\boldsymbol{\alpha}^{*}\in\overline{U(\boldsymbol{\alpha}^{0},t^{*}_{e})}$, and $\mathbf{f}(\boldsymbol{\alpha}^{*},\mathbf{g})=\mathbf{0}$.* 3. (3)

$\mathbf{J}(\boldsymbol{\alpha}^{v},\mathbf{g})^{-1}$* exists for all $v\in\mathbb{N}$ and equation (19) is satisfied.*

Proof.

Immediate application of the Newton-Kantorovich Theorem. ∎

Remark 9*.*

From Assumptions 2 and 3 and from the fact that extended Newton iteration is a valid extension of the sequence (13) then we have that $\lambda_{e}\geq\lambda$, $\varkappa_{e}\geq\varkappa$, $\delta_{e}\geq\delta$, $h_{e}\geq h$. This implies that $t^{*}_{e}\geq t^{*}$ and therefore $U(\mathbf{x}_{0},t^{*})\subseteq U(\boldsymbol{\alpha}_{0},t^{*}_{e})$ for all $\mathbf{g}\in\Psi$ (See Figure 4). From the Newton-Kantorovich Theorem it follows that that $\boldsymbol{\alpha}^{v}\in U(\boldsymbol{\alpha}_{0},t^{*}_{e})$ for all $v\in\mathbb{N}_{0}$.

We can construct a region $\Psi$ such that for all $\mathbf{g}\in\Psi$ the extended Newton iteration converges. Let $\mathbf{y}=\begin{bmatrix}\mathbf{y}_{R}\\ \mathbf{y}_{I}\end{bmatrix}=\begin{bmatrix}\mathop{\text{\rm Re}}{\mathbf{g}}\\ \mathop{\text{\rm Im}}{\mathbf{g}}\end{bmatrix}$, and apply the multivariate Taylor theorem for each $k=1,\dots,n$, $l=1,\dots,n$ entry of the Jacobian matrix $\mathbf{J}_{R}(\boldsymbol{\alpha}_{0},\mathbf{g})$. Evaluating $\mathbf{g}$ at $\mathbf{q}+\mathbf{v}$, we have that

[TABLE]

$\mathbf{v}_{R}:=\mathop{\text{\rm Re}}{\mathbf{v}}$, $\mathbf{v}_{I}:=\mathop{\text{\rm Im}}{\mathbf{v}}$, and

[TABLE]

The entries of the remainder term $\mathbf{R}_{k,l}(\mathbf{x}_{0})$ are bounded by

[TABLE]

where $\partial_{y_{\beta}}$ refers to the derivative of the $\beta^{th}$ variable of the vector $\mathbf{y}$. Form the matrix

[TABLE]

and let

[TABLE]

be the $k=1,\dots,n$, $l=1,\dots,n$ entry of the matrix $\mathbf{Q}$. Then

[TABLE]

where $\mathbb{J}:=\begin{bmatrix}[c]\mathbf{J}_{R}(\mathbf{x}_{0},\mathbf{q})&\mathbf{0}\\ \mathbf{0}&\mathbf{J}_{R}(\mathbf{x}_{0},\mathbf{q})\\ \end{bmatrix}$.

Theorem 6**.**

Suppose that $\varkappa_{e}\geq\varkappa$ and

[TABLE]

whenever $\mathbf{g}\in\Psi$ then

[TABLE]

Proof.

First note that

[TABLE]

From Lemma 2.2.3 in [13], if $\|\mathbb{J}^{-1}\mathbf{E}\|_{2}<1$ then $(\mathbf{I}+\mathbb{J}^{-1}\mathbf{E})$ is invertible and

[TABLE]

Given that $\|\mathbb{J}(\mathbf{x}_{0},\mathbf{q})^{-1}\|\leq\varkappa$ (From Assumption 1 ) whenever $\mathbf{q}\in\Gamma$, it follows

[TABLE]

We conclude that if

[TABLE]

whenever $\mathbf{g}\in\Psi$, then from Equations (21) and (22)

[TABLE]

∎

Applying the multivariate Taylor’s theorem for each $k=1,\dots,n$, entry of the vector $\mathbf{f}_{R}(\boldsymbol{\alpha}_{0},\mathbf{g})$ where $\mathbf{g}=\mathbf{q}+\mathbf{v}$, we have

[TABLE]

where $\mathbf{v}_{R}:=\mathop{\text{\rm Re}}{\mathbf{v}}$, $\mathbf{v}_{I}:=\mathop{\text{\rm Im}}{\mathbf{v}}$, and

[TABLE]

The remainder term $\mathbf{S}_{k}(\mathbf{x}_{0},\mathbf{q})$ is bounded by

[TABLE]

where $\partial_{y_{\beta}}$ refers to the derivative of the $\beta^{th}$ variable of $\mathbf{y}$. We can now rewrite the vector $\mathbf{f}(\boldsymbol{\alpha}_{0},\mathbf{g})$ as

[TABLE]

where $\mathbb{F}:=\begin{bmatrix}[c]\mathbf{f}_{R}(\mathbf{x}_{0},\mathbf{q})\\ \mathbf{0}\\ \end{bmatrix}$, $\mathbf{G}:=\begin{bmatrix}\mathbf{P}(\mathbf{x}_{0},\mathbf{q},\mathbf{v}_{R},\mathbf{v}_{I})\\ \mathbf{f}_{I}(\mathbf{x}_{0},\mathbf{q},\mathbf{v}_{R},\mathbf{v}_{I})\\ \end{bmatrix}$ and

[TABLE]

Theorem 7**.**

Suppose that $\varkappa_{e}\geq\varkappa$ and $\delta_{e}\geq\delta$. Then if

[TABLE]

whenever $\mathbf{g}\in\Psi$, it follows

[TABLE]

Proof.

For each of the entries $k=1,\dots,n$ of the vector $\mathbf{P}$.

[TABLE]

Since $\|\mathbf{E}\|_{2}<\frac{1-\frac{\varkappa}{\varkappa_{e}}}{\varkappa}<\varkappa$ (from Theorem 6) from Lemma 2.2.3 in [13] it follows that $(\mathbf{I}+\mathbb{J}^{-1}\mathbf{E})$ is invertible and

[TABLE]

Combining e quations (23) and (24) we have

[TABLE]

The result follows. ∎

From the values of $\varkappa,\delta,\lambda$ and $\varkappa_{e},\delta_{e},\lambda_{e}$ and Theorems 6 and 7, the region of analyticity $\Psi$ can be constructed. From this region, convergence rates from Theorem 11 for the sparse grid interpolation can be estimated. If we are interested in forming a sparse grid for each of the entries of the vector $\mathbf{x}^{v}\in\mathbb{R}^{n}$, then there are potentially $n$ sparse grids. To estimate the convergence rate of the sequence of sparse grids it is sufficient to embed a polydisk ${\mathcal{E}}_{\sigma_{1},\dots,\sigma_{N}}\subset\Psi$ for a suitable set of coefficients $\{\sigma_{1}\dots,\sigma_{N}\}$. Furthermore since $\boldsymbol{\alpha}^{v}\in U(\alpha_{0},t^{*}_{e})$, the maximal coefficient $\tilde{M}(\mathbf{z}^{v})$ can be bounded as

[TABLE]

4. Application to power flow

The theory developed in Section 3 can be applied to the computation of the statistics of stochastic power flow. In particular we concentrate on the random perturbations of the generators, loads and admittance uncertainty of the transmission lines. Much of the power system network model presented in this section is based on [8].

Consider a network with $m+1$ mechanical constant power generators. The electrical power injected into the network at each generator is given by

[TABLE]

where the operands of the summation are the power from bus $k$ transmitted to bus $l$ through a line with admittance $Y_{k,l}=G_{k,l}+iB_{k,l}$, phase shift $\varphi_{k,l}$, and voltage $V_{k}$ at the buses. These form the algebraic constraints of the power system. The dynamic constraints at generator $k$ are given by

[TABLE]

where $M_{k}$ is the moment of inertia of generator $i$, $D_{k}$ is the damping factor, $P_{M_{k}}$ denotes the mechanical power, $P_{L_{k}}(\omega)$ is the stochastic load and $P_{I_{k}}(\omega)$ is the intermittent stochastic power applied to bus $i$. Equations (25) and (26) constitute the swing equation model. Since the intermittent power generators and loads are stochastic, the rotor angle $\theta_{k}(\omega)$ and power generation $P_{G_{k}}(\omega)$ will be stochastic as well. A simple example of a 3 bus power system is shown in Figure 5. From the steady state response the power flow equations are given by

[TABLE]

for $k=0,\dots,m$.

It is assumed that at each node the active and reactive power injections (or loads) are given by $P_{1},\dots,P_{m}$ and $Q_{1},\dots,Q_{m}$. The first bus is assumed to be slack bus with known angle $\theta_{0}=0$ and fixed voltage $V_{0}$.

Remark 10*.*

According to power system convention, the numbering of the buses (nodes) starts with 1 instead of 0. To simplify the notation in this section we start from 0. However, for the examples and numerical results we revert to the power system standard.

In this paper we limit our discussion of power flow to the case where the power injections $P_{1},\dots,P_{m}$ and $Q_{1},\dots,Q_{m}$ are assumed to be known, but could be stochastic. The unknowns are formed by the angles $\theta_{1},\dots,\theta_{m}$ and voltages $V_{1},\dots,V_{m}$. The power flow equations are solved with a Newton iteration and posed as

[TABLE]

where

[TABLE]

The Jacobian matrix is given in block form as

[TABLE]

where $\mathbf{J}_{11},\mathbf{J}_{12},\mathbf{J}_{21},\mathbf{J}_{22}\in\mathbb{R}^{n\times n}$. For $k,l=1,\dots,m$ let $\theta_{k,l}:=\theta_{k}-\theta_{l}$ and if $k\neq l$

[TABLE]

otherwise

[TABLE]

It is clear from the structure of $\mathbf{f}$ and the Jacobian $\mathbf{J}$ that they are analytic everywhere except for $V_{k}=0$, for $k=1,\dots,m$. However, in practice the domain $\Theta_{0}$ is chosen such the origin is avoided. Otherwise the analyticity assumptions of Theorem 3 are not satisfied.

There are many forms of uncertainty that can be present in the solution of the power flow equations. We concentrate on the following cases:

• •

Random generators and loads: The power injections $P_{k}$ and $Q_{k}$, $k=1,\dots,m$, will be a function of the random vector $\mathbf{q}\in\Gamma$:

[TABLE]

where $P^{0}_{k}$ and $Q^{0}_{k}$ are the nominal power loads (or generators), $q_{k}\in[-1,1]$ and $c_{k},c_{k+1}\in\mathbb{R}$.

• •

Random admittances: The transmission line admittances $Y_{k,l}$ will be functions of the random vector $\mathbf{q}\in\Gamma$. Let ${\mathcal{A}}$ be the set of network index tuples $(k,l)$ such that the admittance is stochastic. Thus for all $k,l\in{\mathcal{A}}$ let

[TABLE]

where $G^{0}_{k,l}$ and $B^{0}_{k,l}$ are the nominal conductance and susceptance, $q_{k,l,1},q_{k,l,2}\in[-1,1]$ and $c_{k,l,1},c_{k,l,2}\in\mathbb{R}$. Note that with a slight of abuse of notation the vector $\mathbf{q}$ consists of all the stochastic random variables $\{q_{k,l,1},q_{k,l,2}\}_{(k,l)\in{\mathcal{A}}}$.

For sufficiently small coefficients $c_{k},c_{k+1}$ with $k=1,\dots,m$ and $c_{k,l,1},c_{k,l,2}$ for all tuples $(k,l)\in{\mathcal{A}}$ the assumptions of Theorems 5, 6 and 7 are satisfied for some initial condition $\mathbf{x}_{0}$ and thus we can justify the use of the sparse grids. A more detailed analysis of the size of these coefficients is left for a future analysis emphasizing the details of power systems.

We test the sparse grid approximation on the New England 39 Bus, 10 Generator, power system model provided from the Matpower 6.0 steady state simulator [30, 19]. In this model buses 1 - 29 are PQ buses, buses 30, 32-39 are generators and bus 31 is the reference (slack).

Two tests are performed. We randomly perturb either the loads or the admittances of the transmission lines. The mean and variance of the voltage $V_{22}$ at bus 22 are computed. The mean $\mathbb{E}[V_{22}]$ and variance $\hbox{var}[V_{22}]$ are computed with the Clenshaw-Curtis isotropic Sparse Grid Matlab Kit [4, 28] for $N=2,4,12$ dimensions and up to the $w=7$ level. This last level, $w=7$ is taken as the “true” solution. The errors are computed up to level $w=4$ with respect to this solution. Two tests are performed:

• •

Random loads: The loads are considered stochastic and are perturbed by up to $\pm$ 50% of their nominal value. For each $k=1,\dots,N$, the so-called $k^{th}$ PQ bus is stochastically perturbed as

[TABLE]

where $P^{0}_{k}$ and $Q^{0}_{k}$ are the nominal power loads, $q_{k}\in[-1,1]$, $\rho(q_{k})$ has a uniform distribution and the random variables $q_{1},\dots,q_{N}$ are independent. Note that although the load random perturbations are independent, the power flows will be dependent on all the random variables $q_{1},\dots,q_{N}$.

In Figure 6 (a) & (b) the mean and variance convergence error for the stochastic voltage $V_{22}$ of bus 22 are shown. A surrogate model based on the sparse grid operator is formed as $\mathcal{S}^{m,g}_{w}[V_{22}]$ with Clenshaw-Curtis abscissas. Each of the circles corresponds to a sparse grid $\mathcal{S}^{m,g}_{w}$ starting with level $w=1$ up to level $w=4$. The y-axis corresponds to the error of the mean or variance. The x-axis is the number of sparse grid knots needed to form the grid $\mathcal{S}^{m,g}_{w}$. The dimension of the sparse grid is given by $N=2,4,12$.

From Figure 6 (a) & (b) we observe that the error decreases faster than polynomially with respect to the number of knots $\eta$. As we increase the number $w$ of levels, sub-exponential convergence is achieved. This is much faster than the $\eta^{-\frac{1}{2}}$ convergence rate of the Monte Carlo method. However, as the number of dimensions $N$ increases the convergence rate of the sparse grid decreases, as predicted by Theorem 11. Moreover, if the level $w$ is not large enough then the error bound gives algebraic convergence.

• •

Random transmission line admittances: The admittances of the network are assumed to be random with

[TABLE]

where $G^{0}_{k,l}$ and $B^{0}_{k,l}$ are the nominal conductance and susceptance. The coefficients $q_{k,l,1},q_{k,l,2}\in[-1,1]$ have a uniform distribution and are all independent. Figure 6 (c) & (d) indicate sub-exponential convergence of the mean and variance of the voltage $V_{22}$ at bus 22 for a sufficiently large number of knots. However, as the number of stochastic dimensions $N$ increases to 12, the convergence rate decreases and almost approaches polynomial convergence. From Theorem 1 the sufficient condition $w>N/\log{2}$ leads to subexponential convergence. For $N=12$ we have that $w$ has to be larger than $18$ to guarantee sub-exponential convergence. In Figures (c) & (d) the largest level for $w$ is 4.

5. Conclusions

In this paper we have introduced ideas from UQ and numerical analysis for the solution of stochastic PDEs using the Newton iteration. In particular we have developed a regularity analysis of the solution with respect to the random perturbations. Under sufficient conditions based on the Newton-Kantorovich Theorem there exists analytic extensions of the solution of the Newton iteration. These indicate that the application of sparse grids for the computation of the stochastic moments leads to sub-exponential or algebraic convergence. For a moderate number of dimensions the convergence rates are much faster than traditional Monte Carlo approaches ($\eta^{-\frac{1}{2}}$). In addition, numerical experiments applied to the power flow problem confirm these subexponential and algebraic convergence rates.

A weakness in the application of the Newton-Kantorovich Theorem is that is constricts the size of the region of analyticity $\Psi$, thus leading to a conservative convergence rate of the sparse grid. This motivates the application of less restrictive methods such as damped Newton iterates [5]. In addition, if we incorporate the assumption that all the Newton iterations converge for each of the knots of the space grid, then by developing an a posteriori method convergence rates can be further improved.

Future work includes the important application of this method to the security constrained problem [26] from the probabilistic perspective. In other words, given stochastic perturbations of the loads and sources what are the optimal power injections into the grid such that the probability of failure is below a tolerance level. Current approaches rely on simplifications of the stochastic perturbations to deal with the high dimensions. However, this can lead to suboptimal results. The high dimensional stochastic quadrature approach developed in this paper will allow more optimal results.