Privacy-Preserving Distributed Parameter Estimation for Probability Distribution of Wind Power Forecast Error

TL;DR

This paper introduces a privacy-preserving distributed algorithm based on Gaussian mixture models to estimate joint and conditional wind power forecast error distributions without sharing raw data.

Contribution

It develops novel privacy-preserving distributed expectation-maximization algorithms using local computations and neighbor communications for wind power error modeling.

Findings

Algorithms successfully estimate joint distribution parameters.

Preserves data privacy during distributed estimation.

Validated with NREL wind data set.

Abstract

Building the conditional probability distribution of wind power forecast errors benefits both wind farms (WFs) and independent system operators (ISOs). Establishing the joint probability distribution of wind power and the corresponding forecast data of spatially correlated WFs is the foundation for deriving the conditional probability distribution. Traditional parameter estimation methods for probability distributions require the collection of historical data of all WFs. However, in the context of multi-regional interconnected grids, neither regional ISOs nor WFs can collect the raw data of WFs in other regions due to privacy or competition considerations. Therefore, based on the Gaussian mixture model, this paper first proposes a privacy-preserving distributed expectation-maximization algorithm to estimate the parameters of the joint probability distribution. This algorithm consists of two original methods: (1) a privacy-preserving distributed summation algorithm and (2) a privacy-preserving distributed inner product algorithm. Then, we derive each WF's conditional probability distribution of forecast error from the joint one. By the proposed algorithms, WFs only need local calculations and privacy-preserving neighboring communications to achieve the whole parameter estimation. These algorithms are verified using the wind integration data set published by the NREL.

Figures

Tables

Table 1. TABLE I: The RSEs between the distributions built by the centralized EM algorithm and the proposed algorithm

Wind Farm	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$
PDF ($\times {10}^{-3}$)	2.4	2.4	2.4	2.4	2.4	2.4	2.4	2.4	2.4
CDF ($\times {10}^{-5}$)	4.8	4.8	4.8	4.8	4.8	4.8	4.8	4.8	4.8

Table 2. TABLE II: KL divergences between the distributions built by the 1st WF and those built by other WFs

Wind Farm	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$
KLD ($\times {10}^{-15}$)	0	0.02	0.19	2.19	1.09	0.33	1.11	2.06	1.01

Equations

$$P\big{(}{[\mathbf{X},\mathbf{Y}]}|\boldsymbol{\theta}\big{)}=\sum\nolimits_{j=1}^{J}w_{j}\mathcal{N}\big{(}{[\mathbf{X},\mathbf{Y}]}|\boldsymbol{\mu}_{j},\boldsymbol{\Sigma}_{j}\big{)}$$

$${\mathbf{\mu }}_j=\left[\begin{array}{cc}{\mathbf{\mu }}_{j,x}& {\mathbf{\mu }}_{j,y}\end{array}\right],{\mathbf{\mu }}_{j,x}=\left[\begin{array}{ccc}{\mu }_{j,1}& \cdots & {\mu }_{j,M}\end{array}\right]$$

$${\mathbf{\mu }}_{j,y}=\left[\begin{array}{ccc}{\mu }_{j,M+1}& \cdots & {\mu }_{j,2M}\end{array}\right],{\mathbf{\Sigma }}_j=\left[\begin{array}{cc}{\mathbf{A}}_j& {\mathbf{B}}_j\\ {\mathbf{B}}_j^T& {\mathbf{C}}_j\end{array}\right]$$

$${\mathbf{B}}_j=\left[\begin{array}{ccc}{\sigma }_{j,1,M+1}& \cdots & {\sigma }_{j,1,2M}\\ ⋮& ⋮& ⋮\\ {\sigma }_{j,M,M+1}& \cdots & {\sigma }_{j,M,2M}\end{array}\right]=\left[\begin{array}{c}{\mathbf{b}}_{j,1}\\ ⋮\\ {\mathbf{b}}_{j,M}\end{array}\right]$$

$${\mathbf{C}}_j=\left[\begin{array}{ccc}{\sigma }_{j,M+1,M+1}& \cdots & {\sigma }_{j,M+1,2M}\\ ⋮& ⋮& ⋮\\ {\sigma }_{j,2M,M+1}& \cdots & {\sigma }_{j,2M,2M}\end{array}\right]$$

$$P({\mathbf{z}}_m|{\mathbf{y}}^0)={\sum }_{j=1}^J{\alpha }_{j,m}{\mathcal{N}}_j({\mathbf{z}}_m+{\mathbf{y}}_m|{\mathbf{y}}^0;{\lambda }_{j,m},{\sigma }_{j,m})$$

$${\alpha }_{j,m}=\frac{w_j{\mathcal{N}}_j({\mathbf{y}}^0;{\mathbf{\mu }}_{j,y},{\mathbf{C}}_j)}{{\sum }_{j=1}^J{w}_j{\mathcal{N}}_j({\mathbf{y}}^0;{\mathbf{\mu }}_{j,y},{\mathbf{C}}_j)}$$

$${\lambda }_{j,m}={\mu }_{j,m}+{\mathbf{b}}_{j,m}{\mathbf{C}}_j^{-1}({\mathbf{y}}^0-{\mathbf{\mu }}_{j,y})$$

$${\delta }_{j,m}={\sigma }_{j,m,m}-{\mathbf{b}}_{j,m}{\mathbf{C}}_j^{-1}{\mathbf{b}}_{j,m}^T$$

$$Q_{j,n}^k=\frac{w_j^{k-1}\mathcal{N}({\mathbf{\varsigma }}_n;{\mathbf{\mu }}_j^{k-1},{\mathbf{\Sigma }}_j^{k-1})}{{\sum }_{j=1}^J{w}_j^{k-1}\mathcal{N}({\mathbf{\varsigma }}_n;{\mathbf{\mu }}_j^{k-1},{\mathbf{\Sigma }}_j^{k-1})}$$

$$w_j^k=\frac{1}{N}\sum _{n=1}^N{Q}_{j,n}^k$$

$${\mathbf{\mu }}_j^k=\frac{{\sum }_{n=1}^N{Q}_{j,n}^k{\mathbf{\varsigma }}_n}{{\sum }_{n=1}^N{Q}_{j,n}^k}$$

$${\mathbf{\Sigma }}_j^k=\frac{{\sum }_{n=1}^N{Q}_{j,n}^k({\mathbf{\varsigma }}_n-{\mathbf{\mu }}_j^k{)}^T({\mathbf{\varsigma }}_n-{\mathbf{\mu }}_j^k)}{{\sum }_{n=1}^N{Q}_{j,n}^k}$$

$$\mathcal{N}({\mathbf{\varsigma }}_n|\mathbf{\mu },\mathbf{\Sigma })=\frac{exp[-\frac{1}{2}({\mathbf{\varsigma }}_n-\mathbf{\mu }){\mathbf{\Sigma }}^{-1}({\mathbf{\varsigma }}_n-\mathbf{\mu }{)}^T]}{\sqrt{(2\pi {)}^{2M}det(\mathbf{\Sigma })}}$$

$${ϵ}_n=({\mathbf{\varsigma }}_n-\mathbf{\mu }){\mathbf{\Sigma }}^{-1}({\mathbf{\varsigma }}_n-\mathbf{\mu }{)}^T$$

$$={\sum }_{i=1}^{2M}(s_{n,i}-{\mu }_i){\sum }_{m=1}^{2M}{\sigma }_{m,i}(s_{n,m}-{\mu }_m)$$

$${\tau }_{n,i}={\sum }_{m=1}^{2M}{\sigma }_{m,i}(s_{n,m}-{\mu }_m)={\sum }_{m=1}^M{d}_{n,i,m}$$

$${ϵ}_n={\sum }_{i=1}^{2M}{\tau }_{n,i}(s_{n,i}-{\mu }_i)={\sum }_{m=1}^M{e}_{n,m}$$

$$d_{n,i,m}={\sigma }_{m,i}(s_{n,m}-{\mu }_m)+{\sigma }_{M+m,i}(s_{n,M+m}-{\mu }_{M+m})$$

$$e_{n,m}={\tau }_m(s_{n,m}-{\mu }_m)+{\tau }_{M+m}(s_{n,M+m}-{\mu }_{M+m})$$

$$G={\sum }_{m=1}^M{l}_m,l_m=d_{n,i,m}\text{or}{l}_m=e_{n,m}$$

$$\displaystyle\mu_{m}=\sum\nolimits_{n=1}^{N}Q_{n}s_{n,m}\bigg{/}\sum\nolimits_{n=1}^{N}Q_{n}=\sum\nolimits_{n=1}^{N}c_{n}s_{n,m}$$

$$\displaystyle c_{n}=Q_{n}\bigg{/}\sum\nolimits_{n=1}^{N}Q_{n}$$

$${\sigma }_{m,m}$$

$$={\sum }_{n=1}^N{c}_n(s_{n,m}-{\mu }_m{)}^2$$

$${\sigma }_{m,i}$$

$$={\sum }_{n=1}^N{c}_n(s_{n,m}-{\mu }_m)(s_{n,i}-{\mu }_i)$$

$$=⟨{\mathbf{p}}_m,{\mathbf{p}}_i⟩,m,i=1,2,...,2M$$

$${\mathbf{p}}_m$$

$${\mathbf{p}}_i$$

$${\alpha }_{m,i}=\{\begin{array}{ll}\frac{1}{\max (|{\mathbf{N}}_m|,|{\mathbf{N}}_i|)+1}& {\nu }_i\in {\mathbf{N}}_m,\\ 1-\sum _{{\nu }_i\in {\mathbf{N}}_m}\frac{1}{\max (|{\mathbf{N}}_m|,|{\mathbf{N}}_i|)+1}& {\nu }_i={\nu }_i,\\ 0& \text{others}\end{array}$$

$$G_m^{t+1}$$

$$={\alpha }_{m,m}{G}_m^t+{\sum }_{{\nu }_i\in N_m}{\alpha }_{m,i}{G}_i^t$$

$$G_m^0$$

$${\Xi }_m^t={\sum }_{{\nu }_i\in N_m}{\alpha }_{m,i}{G}_i^t$$

$$\mathbf{ct}=E_{pk}(\mathbf{pt},H)$$

$$\mathbf{pt}=D_{sk}(\mathbf{ct},H)$$

$$D_{sk}\left(\prod _{m=1}^M{\mathbf{ct}}_m\text{mod}{H}^2\right)=\sum _{m=1}^M{\mathbf{pt}}_m\text{mod}H$$

$$c{t}_i=E_{pk}({\alpha }_{m,i}{G}_i^0,H),\forall {\nu }_i\in {\mathbf{N}}_m$$

$$c{t}_m=E_{pk}(R_0,H),\text{for}{\nu }_m$$

$$c{t}_d=c{t}_m\times {\prod }_{i=1}^{|{\mathbf{N}}_m|}c{t}_i$$

$$D_{sk}\left(c{t}_d\text{mod}{H}^2\right)={\sum }_{i=1}^{|{\mathbf{N}}_m|}{\alpha }_{m,i}{G}_i^0+R_0$$

Full text

I Introduction

As the penetration of wind power continues to increase in multi-regional interconnected grids [1], a better understanding of wind power forecast error is highly desirable [2]. Building the probability distribution of wind power forecast error can benefit both wind farms (WFs) and regional independent system operators (ISOs). For the former, WFs can perform better in market bidding by quantifying stochastic features of their forecast errors [3]. For the latter, regional ISOs can make optimal decisions regarding stochastic economic dispatch [4] or can schedule enough reserves to meet the demand [5]. Note that, if we consider the correlation between WFs when estimating the probability distribution’s parameters, the reserve cost will be significantly reduced since the distributions of forecast errors will be more precise [6].

To take the wind power correlation into account, one should first establish the joint probability distribution of the wind power and the corresponding forecast data from correlated WFs. Then, one can derive the conditional probability distribution of the forecast error under a given forecast value from the joint one [7, 8]. An accurate distribution model is the prerequisite for the above. Therefore, we choose the Gaussian mixture model (GMM) as the distribution model [7], as the GMM can characterize multivariate random variables subject to an arbitrary distribution with remarkable performance [8].

The expectation-maximization (EM) algorithm is the most commonly used method to estimate the parameters of GMM [9, 10]. The centralized EM algorithm requires collecting all the historical wind power and forecast data to train the GMM-based joint probability distribution. After that, the conditional probability distribution of forecast error under a given forecast value can be derived from it. However, the centralized EM algorithm might not be practical. For example, in a multi-regional interconnected grid, the whole grid is actually controlled by multiple regional ISOs, respectively [11]. A central operator with access to all data of the whole grid may not exist for political and technical reasons [12]. Moreover, regional ISOs are unable to collect the raw data from other areas due to privacy considerations [13]. Furthermore, WFs with different stakeholders will not share their forecast values with others, as this may leak their commercial secrets in market bidding [8]. Except for the problem of the possible data barriers, the centralized architecture also faces the risk of single-point failures [14] and the requirement for high bandwidth [15]. To deal with the aforementioned dilemmas, a privacy-preserving distributed (PPD) EM algorithm is a better choice. Specifically, ‘distributed’ means that each WF only needs local calculations and neighboring communications with surrounding WFs, and “privacy-preserving” means that the raw data of a WF cannot be deduced by others during the whole estimation process.

In the data-mining field, many efforts have been made on PPD EM algorithms [16, 17, 18, 19]. The authors of [16] and [17] have proposed a PPD EM algorithm based on the secure sum technique. The algorithm can accurately calculate the sum of data without revealing the data privacy of any parties. A cyclic communication topology is adopted to perform the algorithm. Using the additive homomorphic encryption technique, Kaleb et al. present their PPD EM algorithm by encoding the raw data into cryptographic messages [18]. To prevent the leakage of data privacy when an adversary controls multiple parties, Kaleb et al. enforce one-way communication by a ring topology to guarantee the corruption-resistant feature of the proposed algorithm. Similar to [18], Yang et al. also utilize the additive homomorphic encryption technique to protect the raw data [19]. The differences lie in that (1) Yang et al. design a local-global secure summation protocol, and (2) the cryptographic messages are sent through a spanning tree communication topology in [19].

However, the above PPD EM algorithms mentioned in [16, 17, 18, 19], even including the privacy-free distributed EM algorithms in [20, 21, 22, 23], cannot be applied to the joint probability distribution estimation of correlated WFs because all these algorithms are designed for “horizontally partitioned data”. In fact, wind power and its forecast data are “vertically partitioned” among correlated WFs. Take three parties and 100 samples of 3-dimensional random variables as an example. The data being horizontally partitioned refer to the situation where party A owns 30 samples of 3-dimensional data, party B owns 40 samples, and party C owns 30 samples, while the data being vertically partitioned means that all parties own 100 samples but each of them only has one dimension of the 3-dimensional data. Since each WF has only its own historical wind power and forecast data, i.e., two dimensions of the full multidimensional data, wind power and forecast data are actually vertically partitioned among correlated WFs.

Moreover, the PPD EM algorithms in [16, 17, 18, 19] are not fully distributed. Both the preselected cyclic communication topology in [16, 17, 18] and the preselected spanning-tree communication topology in [19] need a global perspective for preselection. Besides, the failure of any communication line on the preselected path will make the whole algorithm fail.

In this paper, we aim to solve the above problems and develop a PPD EM algorithm to deal with vertically partitioned wind power and forecast data. This algorithm should enable each WF to obtain the GMM-based joint probability distribution in a fully distributed and privacy-preserving manner, and should be robust to communication line failures. After that, using the PPD derivation algorithm proposed in [24], each WF can eventually achieve its conditional probability distribution of forecast error considering the correlation of all WFs. The original contributions of this paper are as follows:

We formulate a distributed framework for the EM algorithm by reformulating the algorithm into local and global parts. Then, the keys to developing a PPD EM algorithm for vertically partitioned wind power and forecast data are pointed out — to find PPD methods to calculate summations and inner products among the statistics of the correlated WFs. 2. 2.

We propose a distributed summation algorithm by leveraging the average consensus algorithm. Then, this algorithm is modified with an additive homomorphic encryption technique to become a PPD summation algorithm. Moreover, we also propose a PPD inner product algorithm on the basis of the randomized binary hash mapping and the average consensus algorithm. 3. 3.

Combining the proposed PPD summation and inner product algorithms, we finally propose a PPD EM algorithm to estimate the GMM-based joint probability distribution. This algorithm is fully distributed and strictly protects the raw data of the correlated WFs from leakage. Meanwhile, this algorithm is robust to communication line failures.

The remainder of this paper is organized as follows. In Section 2, the GMM-based joint and conditional probability distributions are described. In Section 3, the PPD framework for the EM algorithm is formulated to point out the keys to the realization of the PPD EM algorithm. The PPD summation algorithm is developed in Section 4, and the PPD inner product algorithm is designed in Section 5. In Section 6, the PPD EM algorithm is finally proposed. Case studies are provided in Section 7. Section 8 concludes this paper.

II GMM-based Probability Distributions

This section introduces the GMM-based joint and conditional probability distributions. For M spatially correlated WFs, their wind power and wind power forecast constitute a 2M-dimensional random variable $[\mathbf{X},\mathbf{Y}]\in\mathbb{R}^{2M}$, which is defined as $[\mathbf{X},\mathbf{Y}]=[\boldsymbol{x}_{1},...,\boldsymbol{x}_{M},\boldsymbol{y}_{1},...,\boldsymbol{y}_{M}]$. Elements $\boldsymbol{x}_{m}$ and $\boldsymbol{y}_{m}$ represent the wind power and the forecast of the m-th WF, respectively. The GMM-based joint probability distribution function (PDF) of $[\mathbf{X},\mathbf{Y}]$ is actually a convex combination of J multivariate Gaussian distributions with its weighted coefficient $w_{j}\in\mathbb{R}$, mean vector $\boldsymbol{\mu}_{j}\in\mathbb{R}^{2M}$ and covariance $\boldsymbol{\Sigma}_{j}\in\mathbb{R}^{2M\times 2M}$, as given in (1):

[TABLE]

where $\mathcal{N}(\cdot)$ is a 2M-dimensional Gaussian distribution and $\boldsymbol{\theta}=\{w_{j},\boldsymbol{\mu}_{j},\boldsymbol{\Sigma}_{j}|j=1,2,...,J\}$ is the parameter set of GMM-based joint PDF. The details of $\boldsymbol{\theta}$ are given as follows:

[TABLE]

Define the forecast error of the $m$-th WF as $\boldsymbol{z}_{m}=\boldsymbol{x}_{m}-\boldsymbol{y}_{m}$. Given a forecast value $\boldsymbol{y}^{0}\in\mathbb{R}^{M}$, the conditional PDF of $\boldsymbol{z}_{m}$ can be derived from (1) and detailed in (6), where its weighted coefficient $\alpha_{j,m}$ is given in (7), its mean vector $\lambda_{j,m}$ is given in (8) and its variance $\delta_{j,m}$ is given in (9).

[TABLE]

III PPD Framework for the EM Algorithm

In this section, the centralized EM algorithm for estimating the parameters of GMM is first introduced. Then, we formulate a distributed framework for the EM algorithm by reformulating the algorithm into local and global parts. Specifically, the local parts can be calculated by each WF, and the global parts require the results of all local parts to be computed. To design a PPD EM algorithm, the keys lie in how to achieve the global parts in a PPD manner. Therefore, the distributed framework points out the direction for the PPD EM algorithm design.

III-A Centralized EM Algorithm

The training set $\mathbf{S}\in\mathbb{R}^{N\times 2M}$ consists of N historical observations of $[\mathbf{X},\mathbf{Y}]$. The n-th observation is described as $\boldsymbol{\varsigma}_{n}=[s_{n,1},...,s_{n,2M}]$, where $s_{n,m}$ is the n-th wind power observation of the m-th WF, while $s_{n,M+m}$ is its n-th forecast observation. The closed-form expression of the centralized EM algorithm consists of the expectation step (E-step) and maximization step (M-step). In the k-th iteration, the centralized E-step is given in (10) and the centralized M-step in (11):

[TABLE]

[TABLE]

where T represents the transpose of a vector or matrix. After convergence, the estimation of $\boldsymbol{\theta}$ is achieved, and GMM-based joint PDF is established. For detailed derivation and proof, please refer to [25]. It should be emphasized that, since the calculation processes are the same for every Gaussian component in every iteration, we will omit the subscripts k and j in later derivations when it does not cause ambiguity.

III-B Distributed Framework for the E-step

The E-step aims to calculate the statistics $Q_{n}$ in (10) using the parameter $\boldsymbol{\theta}$ updated by each WF from the last iteration. Therefore, in the E-step, $\boldsymbol{\theta}$ becomes public knowledge for all WFs. The part that actually involves $\boldsymbol{\varsigma}_{n}$ exists only in the exponential term of (12), as given in (13).

[TABLE]

Equation (13) shows that for vertically partitioned data, the E-step can be divided into two summations among WFs: the first one is given in (14), and the second one is given in (15).

[TABLE]

The local part of the first summation for the $m$-th WF is defined as $d_{n,i,m}$ in (16), where the $m$-th WF can compute it with the known $\boldsymbol{\theta}$ and its own data. The global part of the first summation is (14) itself.

[TABLE]

The local part of the second summation for the $m$-th WF is defined as $e_{n,m}$ in (17), while the global part of the second summation is (15) itself as well.

[TABLE]

In fact, the relationships between the local and global parts of (14) and (15) are the same. Therefore, we provide a unified form of the relationships in (18).

[TABLE]

To achieve the global part $G$, one should collect $l_{m}$ of others. However, for the m-th WF, sharing $l_{m}$ with others might leak the information of its raw data because the wind power $s_{n,m}$ is close to its forecast value $s_{n,M+m}$; thus, other WFs can estimate the m-th WF’s data to some extent from $d_{n,i,m}$ or $e_{n,m}$. Therefore, how to calculate (18) in a distributed manner under the premise of data privacy preservation is the key to realizing the PPD E-step.

III-C Distributed Framework for the M-step

The M-step aims to update $\boldsymbol{\theta}$ in (11) using $Q_{n}$ calculated from the E-step. For (11a), since $Q_{n}$ is already obtained by all WFs as public knowledge, every WF can directly compute (11a) to update $\omega$. For (11b), its m-th element is reformulated in (19). The m-th WF can compute $\mu_{m}$ and $\mu_{M+m}$ by itself. Meanwhile, since no WF can deduce N observations from the result of (19), the m-th WF can share its $\mu_{m}$ and $\mu_{M+m}$ with other WFs as no data privacy is sacrificed. Finally, each WF can update $\boldsymbol{\mu}$ using the results of (19) from others.

[TABLE]

For (11c), the diagonal and nondiagonal elements of $\boldsymbol{\Sigma}$ are reformulated in (20) and (21), respectively. The m-th WF can directly compute $\sigma_{m,m}$ and $\sigma_{M+m,M+m}$ in (20). Besides, neither $\sigma_{m,m}$ nor $\sigma_{M+m,M+m}$ contain private information because no WF can deduce raw data from them. Thus, the m-th WF can share $\sigma_{m,m}$ and $\sigma_{M+m,M+m}$ with others.

[TABLE]

However, the situation is completely different when calculating $\sigma_{m,i}$ in (21), which requires computing an inner product between vector $\boldsymbol{p}_{m}\in\mathbb{R}^{N}$ and $\boldsymbol{p}_{i}\in\mathbb{R}^{N}$ in (22). The local part of the inner product is defined in (23) for the $m$-th WF and (24) for the $i$-th WF, where the $m$-th or $i$-th WF can directly compute $\boldsymbol{p}_{m}$ or $\boldsymbol{p}_{i}$ with its own data. The global part of the inner product is (22). As can be observed, to obtain the global part, one should collect the local part, i.e., the vector $\boldsymbol{p}_{m}$ of all WFs. However, since $c_{n}$ and $\mu_{m}$ are all public knowledge after the calculations of (10) and (11b), collecting $\boldsymbol{p}_{m}$ is essentially collecting the raw data $s_{n,m}$ of the $m$-th WF, which reveals privacy. Therefore, how to calculate (22) for any two WFs in a distributed manner under the protection of data privacy is the key for realizing the PPD M-step.

IV PPD Summation Algorithm

This section proposes a PPD summation algorithm to calculate (18) in a fully distributed and privacy-preserving manner.

IV-A Average Consensus Algorithm

To calculate (18) in a fully distributed manner, the average consensus algorithm is an effective approach [26]. Some definitions are presented before the demonstration. The communication topology of M spatially correlated WFs is represented by a graph $(\boldsymbol{\nu},\boldsymbol{\xi})$, where $\boldsymbol{\nu}$ denotes the set of nodes $\boldsymbol{\nu}=\{\nu_{1},\nu_{2},...,\nu_{M}\}$ and $\boldsymbol{\xi}$ denotes the set of edges $\boldsymbol{\xi}\subseteq\nu\times\nu$. Once the distance between two nodes is less than a preset distance threshold $\eta$, the two nodes are connected. The neighbors of node m are denoted by $\boldsymbol{N}_{m}=\{\nu_{i}\in\boldsymbol{\nu}:(\nu_{m},\nu_{i})\in\xi\}$. The weighted adjacency matrix is represented by $\boldsymbol{A}\in\mathbb{R}^{M\times M}$ with adjacent coefficient $\{\alpha_{m,i}|m,i=1,...,M\}$ as given in (25), where $\left|\boldsymbol{N}_{m}\right|$ and $\left|\boldsymbol{N}_{i}\right|$ denote the degree of nodes m and i. $\boldsymbol{A}$ is a symmetric matrix, and $\boldsymbol{A}\boldsymbol{1}=\boldsymbol{1}$, where $\boldsymbol{1}=\left[1,...,1\right]^{T}\in\mathbb{R}^{M}$.

[TABLE]

The discrete form of the average consensus algorithm is presented in (26). After convergence, each WF can obtain the average value of (18) in a distributed manner, as given in (27).

[TABLE]

To compute (26), the $m$-th WF only needs to collect $\alpha_{m,i}G_{i}^{t}$ (for $\nu_{i}\in N_{m}$) from its neighbors to calculate a local summation during each iteration, i.e.,

[TABLE]

However, in the first iteration, $G_{i}^{0}=l_{i}$ (for $\nu_{i}\in N_{m}$) is revealed to the $m$-th WF. Thus, the average consensus algorithm is not privacy-preserving.

IV-B PPD Summation Algorithm

To achieve the local summation in (28) under the premise of protecting privacy, we leverage an additive homomorphic encryption technique named Paillier cryptosystem. The Paillier cryptosystem is favored by many researchers for the analysis of social networks [19] or the Internet [27].

Let $\boldsymbol{pt}\in\mathbb{Z}$ denote a plaintext, $\boldsymbol{ct}\in\mathbb{Z}$ denote a ciphertext and $H$ denote a prespecified large prime integer. The encryption process with a public key $pk$ is given in (29), and the decryption with a secret key $sk$ is given in (30).

[TABLE]

To compute the sum of M plaintexts, the decrypter only needs to multiply the corresponding M ciphertexts and then decrypt the multiplication result, as given in (31). The entire process strictly protects data privacy. See [28] for more details.

[TABLE]

Inspired by the secure summation protocol in [19], we utilize the Paillier cryptosystem to compute (28), which helps us to realize a privacy-preserving average consensus algorithm. Specifically, in the first iteration of the average consensus algorithm, the neighbors of the m-th WF, which are numbered from 1 to $|N_{m}|$, encrypt their initial value using the 1st neighbor’s $pk$ in (32). Meanwhile, the m-th WF encrypts a random and secret number $R_{0}$ by (33). Then, these neighbors send their ciphertexts to the m-th WF. After that, the m-th WF performs the multiplication calculation in (34) and sends the result to the 1st neighbor. Thereafter, the 1st neighbor decrypts $ct$ into the summation in (35) and sends it back to the m-th WF. Finally, the m-th WF subtracts $R_{0}$ to obtain the result of (28). For the subsequent iterations of the average consensus algorithm, no encryption is needed. Details of the PPD summation algorithm are given in Algorithm 1.

[TABLE]

Please note that although the Paillier cryptosystem is introduced in the average consensus algorithm, it does not affect the convergence of the average consensus algorithm. In fact, the PPD summation algorithm and the average consensus algorithm are mathematically equivalent. See [29] for the convergence proof of the original average consensus algorithm.

V PPD Inner Product Algorithm

In this section, a PPD inner product algorithm is proposed to calculate (22) for any two WFs in a fully distributed manner considering privacy protection.

For the inner product calculation in (22), once the angle $\beta_{m,i}$ in (36) is obtained, by sharing the norm $\|\boldsymbol{p}_{m}\|\in\mathbb{R}$ and $\|\boldsymbol{p}_{i}\|\in\mathbb{R}$ among WFs, which will not reveal any raw data, the inner product can be directly calculated by every WF. Therefore, the problem of computing (22) becomes into this one: how to compute the angle between two vectors in (36) without revealing any raw element of the vectors?

[TABLE]

In an N-dimensional space, the probability of finding a random hyperplane separating two vectors $\boldsymbol{p}_{m}$ and $\boldsymbol{p}_{i}$ is proportional to the angle $\beta_{m,i}$ [30]. For calculating the probability, a publicly known random vector set $\boldsymbol{\Gamma}\in\mathbb{R}^{N\times L}$ is first defined, where each column is a random vector $\boldsymbol{\gamma}_{l}\in\mathbb{R}^{N}$. Then, the probability can be computed by (37).

[TABLE]

For a further demonstration, the randomized binary hash mapping function $h\colon\mathbb{R}^{N}\mapsto\mathbb{R}^{L}$ is defined in (38),

[TABLE]

where the $sign$ function can encode an L-dimensional real vector into an L-dimensional binary vector according to the sign of each element in the real vector. Thus, $h\left(\boldsymbol{p}_{m}\right)$ actually represents the sign information of the multiplication results between $\boldsymbol{p}_{m}$ and $\forall\boldsymbol{\gamma}_{l}\in\boldsymbol{\Gamma}$.

Once $h\left(\boldsymbol{p}_{m}\right)$ and $h\left(\boldsymbol{p}_{i}\right)$ are obtained, (37) can be easily computed by counting the number of different sign pairs. Note that the counting process is essentially calculating the Hamming distance between $h\left(\boldsymbol{p}_{m}\right)$ and $h\left(\boldsymbol{p}_{i}\right)$, i.e., $H_{am}[h\left(\boldsymbol{p}_{m}\right),h\left(\boldsymbol{p}_{i}\right)]$ [31]. Therefore, based on the randomized binary hash mapping function and Hamming distance calculation, the angle $\beta_{m,i}$ can be calculated by (39).

[TABLE]

It should be emphasized that, our goal is not only to calculate the inner product of two vectors under the premise of protecting privacy but also to obtain all the inner products between any two WFs through a fully distributed manner. For computing all the inner product values, the set $\{\|\boldsymbol{p}_{m}\|,h\left(\boldsymbol{p}_{m}\right)|m=1,...,2M\}$ is required. Thus, based on (39) and the average consensus algorithm, the PPD inner product algorithm is proposed in Algorithm 2.

VI PPD EM Algorithm

Since the two keys mentioned in Section 3 are solved by the proposed PPD summation and inner product algorithm, the PPD EM algorithm for estimating the GMM-based joint PDF of multiple spatially correlated WFs is eventually developed in Algorithm 3.

The advantages of the proposed algorithm are as follows:

• •

Strict privacy-preserving. For the summation calculation in the PPD E-step, the Paillier cryptosystem is used to protect the raw data; for the inner product calculation in the PPD M-step, the randomized binary hash mapping function is used to prevent data privacy disclosure. The two techniques that we utilized can strictly protect the privacy of each WF.

• •

Fully distributed. As we introduce the average consensus algorithm into the PPD E-step and M-step, each WF only needs to communicate with its neighbors. Thus, we avoid the assumption made in [27, 16, 32] that any two nodes are connected, and we improve the scalability of the proposed algorithm. Meanwhile, the preselected communication paths in [16, 17, 18, 19] are no longer required. Thus, the proposed algorithm is fully distributed.

• •

Robust. As the communication between neighbors may fail, robustness to communication failure is necessary. Since the proposed PPD EM algorithm is developed based on the average consensus algorithm, as long as the communication topology is still connected, the communication failure basically will not affect the final estimation results due to the consensus feature [33].

Please note that once the joint PDF is established, each WF can derive its conditional PDF of the forecast error in (6) via the PPD derivation algorithm presented in [24].

VII Case Study

The historical data of wind power and forecast value are from the “eastern wind integration data set” published by the National Renewable Energy Laboratory (NREL), the U.S., where we choose 9 WFs in Maryland, numbered 4401, 5405, 6211, 6359, 6526, 6812, 6931, 7187, and 7460. Their communication topology is shown in Fig. 1(a). Each WF has 20 days of hourly wind power and forecast data. Thus, $M=9$ and $N=480$. Then we build the joint PDF of the wind power and its forecast of the nine spatially correlated WFs by the proposed EM algorithm. After that, by leveraging the PPD derivation algorithm presented in [24], we also derive the conditional PDF of the forecast error of each WF from the joint PDF. Since the privacy-preserving feature of the proposed algorithm has already been discussed in the previous sections, this section mainly aims to verify the correctness and robustness of the proposed algorithm.

VII-A Correctness of the PPD Summation Algorithm

We use the wind power data of 9 WFs at 2004/1/1-01:00 as input, and use the proposed PPD summation algorithm to estimate the summation of the 9 data points. To show more details, we illustrate the iterative process of each WF’s estimation for the summation in Fig. 1(b). It can be observed that after 30 iterations, all WFs achieve consensus on the exact summation value, showing the correctness of the proposed PPD summation algorithm.

VII-B Correctness of the PPD Inner Product Algorithm

To verify the correctness of the proposed PPD inner product algorithm, we use it to calculate the inner products between every two WFs’ private vectors. The private vector of a WF consists of its 20-day historical wind power data. Then, compared to the real inner product results, we provide the average relative errors of the proposed algorithm under different values of $L$ in Fig. 2. It can be observed that the error decreases significantly as $L$ increases, but it nearly stabilizes when $L$ reaches $2^{11}$. Thus, we finally choose $L$ as $2^{11}$ bit $=0.25$ kb. Furthermore, the average relative error is $3.5\times 10^{-3}$ when $L=2^{11}$, which proves the correctness of the proposed algorithm.

VII-C Verification of the PPD EM Algorithm

We use the Bayesian information criterion (BIC) to set the number of Gaussian components $J$ as 5. After that, we build the joint probability distribution of wind power and forecast of the nine spatially correlated WFs using the proposed PPD EM algorithm. The distribution constructed by the centralized EM algorithm is also given as the benchmark. Since the 18-dimensional distribution cannot be drawn for illustration, we derive several 1-dimensional and 2-dimensional distributions from the 18-dimensional one based on the linear invariance property of GMM [8].

First, the 1-dimensional PDF and the 1-dimensional cumulative distribution function (CDF) are shown in Fig. 3. Only the first dimension is provided.

In Fig. 3, the empirical distributions are obtained from the corresponding original historical data, the benchmarks are built by the centralized EM algorithm, and the distribution of each WF is constructed by the proposed PPD EM algorithm. It can be observed that (1) the benchmark and the distributions built by the WFs match the corresponding empirical distributions; (2) the distributions built by the WFs are coincident with each other, indicating that the consensus among WFs is achieved by the proposed algorithm; and (3) the distribution built by each WF is coincident with the benchmark, indicating the correctness of the proposed algorithm.

Then, further comparisons between the benchmark and the distribution of each WF are made using the relative standard error (RSE) as defined in (40), where $f(\cdot)$ represents the PDF or CDF built by each WF, $f_{0}(\cdot)$ represents the benchmark PDF or CDF, and $\overline{f_{0}(\cdot)}$ represents the mean value. The RSE results are provided in Table. I. First, all the RSE values are less than $2.4\times 10^{-3}$, which means that the difference between the benchmark and the distribution of each WF is tiny. Second, the consensus effect of the proposed algorithm is obvious as the RSE values are identical correspondingly. Third, the RSE values of the CDFs are much smaller than those values of the PDFs. Note that the CDF is what we ultimately want for optimal decisions, e.g., calculating the quantile from the CDF. Thus, the RSE values of the CDFs eventually represent the accuracy of the proposed algorithm.

[TABLE]

Furthermore, we also choose two dimensions from the 18-dimensional joint distribution to form 2-dimensional PDF and CDF. The 2-dimensional benchmarks built by the centralized EM algorithm and the 2-dimensional joint distributions built by the 1st WF are illustrated in Fig. 4.

Thereafter, the Kullback–Leibler divergence (KLD) between the joint distribution built by the 1st WF and by other WFs are given in Table II to illustrate the differences between them. Since all the KLDs are less than $2.19\times 10^{-15}$, the distribution built by the 1st WF and those built by others are basically the same. Thus, using the result of the 1st WF as a representative is reasonable and acceptable. We can observe that, in Fig. 4, the 2-dimensional PDF and CDF built by the 1st WF match the 2-dimensional benchmarks well. Therefore, the correctness of the proposed PPD EM algorithm is eventually verified.

VII-D Robustness of the PPD EM Algorithm

Since the proposed PPD EM algorithm is developed based on the average consensus algorithm, as long as the communication network topology is still connected, the communication failure basically will not affect the final estimation results due to the consensus feature.

To verify the robustness of the proposed PPD EM algorithm, we first cut off communication lines to simulate communication failures. Then, we inspect the variations of the CDFs built by the proposed algorithm after the failures. Since the consensus of the proposed algorithm has already been verified, we still use the estimation result of the 1st WF as the representative. The CDFs built by the 1st WF after the communication failures are shown in Fig. 5. In this figure, the benchmarks are the CDFs built by the centralized EM algorithm. Besides, legend ‘Original’ represents the CDFs built by the 1st WF when no failure occurs, and legend ‘line m–i’ represents the CDFs built by the 1st WF when the communication between the m-th WF and the i-th WF fails. For example, legend ‘line 1-3’ means that the communication between the 1st WF and the 3rd WF fails while other communication lines operate normally. As we can see, in Fig. 5, the CDFs built by the 1st WF under different communication failures still coincide with the corresponding benchmarks and original CDFs. This proves that the proposed PPD EM algorithm can still maintain high accuracy after communication failures. Therefore, the robustness of the proposed algorithm is verified.

VII-E Verification of the Conditional Distribution

Based on the established joint distribution via the proposed PPD EM algorithm, we derive each WF’s conditional distribution of forecast error in (6) by the PPD derivation algorithm proposed in [24].

The PPD derivation algorithm can enable each WF to obtain its conditional distribution from the joint one in the PPD manner. Here we illustrate the conditional CDF obtained by the 1st, 3rd, 5th, and 7th WFs in Fig. 6, where the benchmarks are all built in a centralized manner. The matches between the benchmarks and the CDFs built by WFs show that the joint distribution obtained via the proposed PPD EM algorithm is correct, so the conditional distributions derived from the joint distribution are correct as well.

VIII Conclusion

Under the consideration of wind power correlation, estimating the conditional probability distribution of WF’s forecast error requires the historical wind power and forecast data of all WFs. However, for the multi-correlated WFs spreading over a multi-regional interconnected grid, data barriers among the WFs belonging to different regions may exist. Therefore, we propose a PPD EM algorithm to estimate the joint probability distribution of the vertically partitioned wind power and forecast data. Then, we derive each WF’s conditional probability distribution of forecast error from the joint one. To achieve this, we first formulate a distributed framework for the EM algorithm by reformulating the algorithm into local and global parts. Hereafter, the keys to developing a PPD EM algorithm are pointed out: calculating summations and inner products among the statistics of the correlated WFs in a PPD manner. After that, we design a PPD summation algorithm based on the additive homomorphic encryption and the average consensus algorithm. Thereafter, we develop a PPD inner product algorithm by leveraging the randomized binary hash mapping function and the average consensus algorithm. Combining the PPD summation and inner product algorithms, we finally propose the PPD EM algorithm. This algorithm can enable each WF to estimate the joint probability distribution of the wind power and forecast data of all the WFs in a fully distributed and privacy-preserving manner.

Compared with the centralized EM algorithm, the proposed algorithm not only has high accuracy but also is fully distributed because it only needs local communication between neighboring WFs. Moreover, it strictly protects the data privacy of every WF during communications. Furthermore, its robustness to communication failure is guaranteed by the consensus feature.