Chance-constrained OPF: A Distributed Method with Confidentiality Preservation

We start by providing a description of the notations used in the CC-OPF formulation. Let $\mathcal{M}$ be the set of regions in the grid. The dispatchable generator set, non-dispatchable wind farm set, load set, and transmission line set of the whole grid are represented by $\mathcal{G}$, $\mathcal{W}$, $\mathcal{D}$, and $\mathcal{L}$, respectively. Correspondingly, the sets of generators, wind farms, loads, and lines in area $n$ are denoted by $\mathcal{G}_{n}$, $\mathcal{W}_{n}$, $\mathcal{D}_{n}$, and $\mathcal{L}_{n}$, respectively. This paper uses $\mathcal{T}$ as the set of time slots, and uses $|\cdot|$ to express the cardinality calculation of any set. Besides, the active power output of dispatchable generators, active wind power outputs, reactive wind power outputs, active loads, reactive loads, and active line flows in area $n$ at time $t$ are respectively denoted by $\boldsymbol{g}_{n,t}\in\mathfrak{R}^{|\mathcal{G}_{n}|}$, $\boldsymbol{w}_{n,t}\in\mathfrak{R}^{|\mathcal{W}_{n}|}$, $\boldsymbol{\widehat{w}}_{n,t}\in\mathfrak{R}^{|\mathcal{W}_{n}|}$, $\boldsymbol{d}_{n,t}\in\mathfrak{R}^{|\mathcal{D}_{n}|}$, $\boldsymbol{\widehat{d}}_{n,t}\in\mathfrak{R}^{|\mathcal{D}_{n}|}$, and $\boldsymbol{l}_{n,t}\in\mathfrak{R}^{|\mathcal{L}_{n}|}$. We further introduce $N_{a}=|\mathcal{M}|$, $T=|\mathcal{T}|$, $H_{n}=|\mathcal{G}_{n}|\times|\mathcal{T}|$, and $H=|\mathcal{G}|\times|\mathcal{T}|$.

The objective of the CC-OPF is to minimize the overall generation cost across the entire grid given the non-dispatchable wind power:

where $\boldsymbol{\mathcal{I}}_{n}\in\mathfrak{R}^{|\mathcal{G}_{n}|\times|\mathcal{G}_{n}|}$ is a diagonal matrix composed of the quadratic cost coefficients of dispatchable generators in region $n$, whose non-zero elements are all positive, rendering (1) strictly convex. Vector $\boldsymbol{\mathcal{K}}_{n}\in\mathfrak{R}^{|\mathcal{G}_{n}|}$ consists of the linear cost coefficients of dispatchable generators in region $n$, whose elements are assumed to be positive.

Due to the non-dispatchable wind power, the supply-demand constraint at time $t$ is formulated as an SCC, as in [29, 30]:

This constraint enforces the supply to be greater than or equal to the demand with a probability of at least $1-\epsilon_{b}$. Vector $\mathds{1}_{|\mathcal{G}_{n}|}$ is an all-ones vector of size $|\mathcal{G}_{n}|$. Similarly, $\mathds{1}_{|\mathcal{W}_{n}|}$ and $\mathds{1}_{|\mathcal{D}_{n}|}$ are also all-ones vectors with appropriate dimensions. Besides, the wind power outputs $\boldsymbol{w}_{n,t}$ are random variables whereas $\boldsymbol{d}_{n,t}$ is constant, that is, this paper ignores the uncertainty of load because it is relatively small compared to wind power uncertainty [31].

Constraints for the dispatchable generators in area $n$ at time $t$ include the capacity and ramp rate constraints:

where $\boldsymbol{g}_{n}^{+}$ and $\boldsymbol{g}_{n}^{-}$ represent the upper and lower generation bounds of generators in area $n$. Similarly, $\boldsymbol{r}_{n}^{+}$ and $\boldsymbol{r}_{n}^{-}$ are the ramp rate bounds of the generators in area $n$. All of these bounds are vectors in $\Re^{|\mathcal{G}_{n}|\times 1}$.

Constraints for system states (e.g., nodal voltages or line flows) are also formulated as SCCs. Let $s_{n,i,t}$ denote state $i$ in area $n$ at time $t$, which could be either a nodal voltage or a line flow. The constraint for $s_{n,i,t}$ is given as follows:

where $s_{n,i}^{+}\in\mathfrak{R}$ is the upper bound of $s_{n,i,t}$, and set $\mathcal{S}_{n}$ consists of the states in region $n$. Constraint (3) enforces $s_{n,i,t}$ to be within its bound with a probability of at least $1-\alpha_{S}$.

Converting the SCCs in (2) and (3) into their deterministic versions is the prerequisite for solving the CC-OPF problem. To this end, this paper first adopts a linear but sufficiently accurate power flow model, namely the decoupled linearized power flow (DLPF) model proposed in [32], whose performance has been verified in [33, 34]. Accordingly, $s_{n,i,t}$, i.e. line flow or voltage, can be expressed as a linear function of power injections across regions:

where $\boldsymbol{\mathcal{U}}_{n,i,m}\in\mathfrak{R}^{|\mathcal{G}_{m}|}$ consists of the mapping coefficients (derived from the DLPF model) that map $\boldsymbol{g}_{m,t}$ to $s_{n,i,t}$. Similarly, $\boldsymbol{\mathcal{O}}_{n,i,m}\in\mathfrak{R}^{|\mathcal{W}_{m}|}$, $\boldsymbol{\widehat{\mathcal{O}}}_{n,i,m}\in\mathfrak{R}^{|\mathcal{W}_{m}|}$, $\boldsymbol{\mathcal{Y}}_{n,i,m}\in\mathfrak{R}^{|\mathcal{D}_{m}|}$, and $\boldsymbol{\widehat{\mathcal{Y}}}_{n,i,m}\in\mathfrak{R}^{|\mathcal{D}_{m}|}$ also consist of the corresponding mapping coefficients. Note that $\xi_{n,i,t}$ is a constant value computed from all known system states, e.g., the voltage magnitude/angle of the slack node, voltage magnitudes of $PV$ nodes, etc. For more details regarding the DLPF model, please see [32].

Using (II-B), the SCCs in (2) and (3) can be converted into deterministic constraints using the quantile calculation method proposed in [13]. For (2), the resulting deterministic constraint is given by

where

and $\mathbb{P}_{\omega_{t}}^{-1}(\epsilon_{b})$ is the $\epsilon_{b}$-quantile of $\omega_{t}$.

Similarly, the SCC in (3) is transformed into:

where

and $\mathbb{P}_{\Omega_{n,i,t}}^{-1}\left(1-\alpha_{S}\right)$ is the $(1-\alpha_{S})$-quantile of $\Omega_{n,i,t}$.

After the above conversions, the CC-OPF problem becomes a quadratic programming problem with linear constraints. Its compact formulation, defined as ($\rm P0$), is given as follows:

Matrix $\boldsymbol{Q}$ in the objective function is composed of the quadratic cost coefficients of dispatchable generators:

where $\boldsymbol{\Lambda}_{n}$ is detailed by

Furthermore, $\boldsymbol{c}$ in the objective function consists of the linear cost coefficients of dispatchable generators of the whole grid, i.e.

where $\boldsymbol{c}_{n}$ is defined as follows:

The decision variable $\boldsymbol{x}$ includes the active dispatchable generations of the whole grid for all time steps:

where

Before describing $\boldsymbol{\mathbb{A}}$ and $\boldsymbol{\mathbb{B}}$ in $(\rm{P0})$, several parameter matrices are first defined, including $\boldsymbol{I}_{n}\in\mathfrak{R}^{|\mathcal{G}_{n}|\times|\mathcal{G}_{n}|}$, $\boldsymbol{J}_{n}\in\mathfrak{R}^{H_{n}\times H_{n}}$, $\boldsymbol{E}_{n}\in\mathfrak{R}^{(H_{n}-|\mathcal{G}_{n}|)\times H_{n}}$, $\boldsymbol{K}_{n}\in\mathfrak{R}^{T\times H_{n}}$, $\boldsymbol{\Phi}_{n}\in\mathfrak{R}^{|\mathcal{S}|\times|\mathcal{G}_{n}|}$, and $\boldsymbol{U}_{n}\in\mathfrak{R}^{T|\mathcal{S}|\times H_{n}}$, where $\mathcal{S}=\bigcup_{n\in\mathcal{M}}\mathcal{S}_{n}$. Among these parameter matrices, $\boldsymbol{I}_{n}$ is an identity matrix; other parameter matrices are specified as follows:

Additionally, the following parameter vectors are introduced:

Based on the above parameter matrices and vectors, $\boldsymbol{\mathbb{A}}$ and $\boldsymbol{\mathbb{B}}$ are given by

Ideally, to formulate $(\rm P_{0})$, each regional ISO or a centralized entity should first form $\boldsymbol{Q}$, $\boldsymbol{c}$, $\boldsymbol{\mathbb{A}}$, and $\boldsymbol{\mathbb{B}}$. However, formulating these matrices or vectors requires knowledge of system data of the entire grid. For example, if computations are done by regional ISOs, then each of them has to collect the generation cost information of generators in other areas to form $\boldsymbol{Q}$ and $\boldsymbol{c}$, the capacities and ramp rates of generators in other regions to build $\boldsymbol{\mathbb{B}}$, the system topologies and line parameters of other regions to form $\boldsymbol{U}_{n}$ ($\forall n$) in $\boldsymbol{\mathbb{A}}$. Such information may however be confidential.

Besides, once regional ISOs solve $(\rm P_{0})$, they would infer additional information from the solution, which includes the dispatchable generation of each area. This information, however, might be confidential as well.

Ignoring the quadratic term in ($\rm P0$) leads to a linear programming problem $(\rm{L0})$:

The basic idea of the TE technique is to utilize random matrices to encrypt the decision variable (i.e., dispatchable generation) of each region, which, in the meantime, also masks the confidential information of this area. Specifically, regional ISO $n$ first randomly generates an invertible matrix $\boldsymbol{M}_{n}\in\mathfrak{R}^{H_{n}\times H_{n}}$. This matrix, called the encryption matrix, is only known by ISO $n$. Then, ISO $n$ transforms its own decision variable $\boldsymbol{x}_{n}$ into the encrypted decision variable $\boldsymbol{\bar{x}}_{n}\in\mathfrak{R}^{H_{n}\times 1}$ using $\boldsymbol{M}_{n}$. The relationship between $\boldsymbol{x}_{n}$ and $\boldsymbol{\bar{x}}_{n}$ is

If we define

then the following equation holds:

Substituting (67) into $(\rm L0)$ yields:

Noticeably, after the transformation, $\boldsymbol{x}_{n}$ becomes $\boldsymbol{\bar{x}}_{n}$. This means that, once the transformed problem is solved, only ISO $n$ can recover $\boldsymbol{x}_{n}$ from $\boldsymbol{\bar{x}}_{n}$ using (62). Also, the generation cost $\boldsymbol{c}_{n}$ is encrypted into $\boldsymbol{c}_{n}\boldsymbol{M}_{n}$ while $\boldsymbol{\mathbb{A}}_{n}$ is encrypted into $\boldsymbol{\mathbb{A}}_{n}\boldsymbol{M}_{n}$ after the transformation. Consequently, ISO $n$ can freely share $\boldsymbol{c}_{n}\boldsymbol{M}_{n}$ and $\boldsymbol{\mathbb{A}}_{n}\boldsymbol{M}_{n}$ with others, as no one can deduce $\boldsymbol{c}_{n}$ or $\boldsymbol{\mathbb{A}}_{n}$ from the exchanged information. Overall, $\boldsymbol{x}_{n}$, $\boldsymbol{c}_{n}$, and $\boldsymbol{\mathbb{A}}_{n}$ in $(\rm{L0})$ can be effectively masked by the encryption matrix $\boldsymbol{M}_{n}$.

However, $\boldsymbol{M}_{n}$ cannot obscure the information in $\boldsymbol{\mathbb{B}}$. But the way to protect each ISO’s confidential data in $\boldsymbol{\mathbb{B}}$ is straightforward. For instance, to protect $\boldsymbol{G}_{n}^{+}$ in $\boldsymbol{\mathbb{B}}$ from disclosure, ISO $n$ can equivalently modify the corresponding inequality constraint from $\boldsymbol{J}_{n}\boldsymbol{x}_{n}\leq\boldsymbol{G}^{+}_{n}$ to:

where $\boldsymbol{\bar{J}}_{n}=\boldsymbol{\varphi}_{n}\boldsymbol{J}_{n}$ and $\boldsymbol{\bar{G}}^{+}_{n}=\boldsymbol{\varphi}_{n}\boldsymbol{G}^{+}_{n}$. Matrix $\boldsymbol{\varphi}_{n}\in\mathfrak{R}^{H_{n}\times H_{n}}$, randomly generated and privately hold by ISO $n$, is a diagonal matrix whose non-zero elements are positive. After introducing $\boldsymbol{M}_{n}$, (68) further becomes:

As a result, the capacity information of generators in region $n$ is encrypted into $\boldsymbol{\bar{G}}^{+}_{n}$. Similar operations can be done to $\boldsymbol{G}^{-}_{n}$, $\boldsymbol{R}^{+}_{n}$, and $\boldsymbol{R}^{-}_{n}$ to encrypt them into $\boldsymbol{\bar{G}}^{-}_{n}$, $\boldsymbol{\bar{R}}^{+}_{n}$, and $\boldsymbol{\bar{R}}^{-}_{n}$, respectively. Simultaneously, $\boldsymbol{J}_{n}$ and $\boldsymbol{E}_{n}$, i.e., the coefficients in the corresponding inequality constraints, are also modified to $\boldsymbol{\bar{J}}_{n}$ and $\boldsymbol{\bar{E}}_{n}$, respectively. For distinction, this paper uses $\boldsymbol{\bar{\mathbb{A}}}$ and $\boldsymbol{\bar{\mathbb{B}}}$ to separately denote $\boldsymbol{\mathbb{A}}$ and $\boldsymbol{\mathbb{B}}$ after the above modifications.

Finally, the original linear problem $(\rm L0)$ is transformed into its encrypted but equivalent version $(\rm L1)$:

Consequently, ISO $n$ can share its encrypted information, including $\boldsymbol{c}_{n}\boldsymbol{M}_{n}$, $\boldsymbol{\bar{\mathbb{A}}}_{n}\boldsymbol{M}_{n}$, $\boldsymbol{\bar{G}}^{+}_{n}$, $\boldsymbol{\bar{G}}^{-}_{n}$, $\boldsymbol{\bar{R}}^{+}_{n}$, and $\boldsymbol{\bar{R}}^{-}_{n}$, with other ISOs to allow each to build $(\rm L1)$. Once $(\rm L1)$ is solved, only ISO $n$ can derive its desired optimal solution, $\boldsymbol{x}_{n}^{*}$, from the encrypted optimal solution $\boldsymbol{\bar{x}}^{*}$. Clearly, no confidential data is leaked.

Although the TE technique is a confidentiality-preserving optimization method without parameter tunings, this method cannot be directly applied to the CC-OPF problem, i.e., $(\rm P_{0})$. More specifically, there are three unsolved challenges standing in the way:

1. 1.

   As mentioned earlier, the TE technique is only proven to be effective for linear programming so far, but $(\rm P0)$ is a quadratic programming problem. Hence, the first challenge is to prove the method’s effective application to quadratic problems.

2. 2.

   To use the TE technique, ISO $n$ must know $\boldsymbol{\bar{\mathbb{A}}}_{n}$. In $(\rm P0)$, however, ISO $n$ is unable to generate the complete $\boldsymbol{\bar{\mathbb{A}}}_{n}$ independently, because the sub matrix $\boldsymbol{U}_{n}$ in $\boldsymbol{\bar{\mathbb{A}}}_{n}$, known as the coupled coefficient matrix among dispatchable generations across regions, consists of the mapping coefficients that map nodal power injections to all the system states. These mapping coefficients are not locally obtainable unless global information, including the system topology and line parameters, is known. ISO $n$ does not have the global information, making $\boldsymbol{U}_{n}$ unobtainable and thereby the TE technique unusable. Given this, how to enable ISO $n$ to obtain $\boldsymbol{U}_{n}$ under the protection of regions’ confidentiality, is the second unsolved challenge.

3. 3.

   Both $\boldsymbol{\Upsilon}$ and $\boldsymbol{\Delta}$ in $\boldsymbol{\bar{\mathbb{B}}}$ are aggregations of the regions’ confidential information, including load data, line parameters, etc. Hence, the third challenge is to allow ISO $n$ to acquire $\boldsymbol{\Upsilon}$ and $\boldsymbol{\Delta}$ without having access to other ISOs’ confidential data.

The next section will address these challenges and provide solutions to each of them.

After introducing the encryptions mentioned in (67) and (68), the original CC-OPF problem $(\rm P0)$, a typical quadratic programming problem, is transformed into the encrypted version $(\rm P1)$:

where

Proving the TE technique’s adaptability to $(\rm P0)$, is equivalent to proving the following two propositions.

As $\boldsymbol{M}_{n}$ is invertible, we can form

meaning that $\boldsymbol{M}$ is also invertible. Accordingly, $\rm rank(\boldsymbol{M}^{-1})=\textit{H}$ holds, where $\rm rank(\cdot)$ denotes the rank function. Therefore, in the space of $\mathfrak{R}^{H\times 1}$, matrix $\boldsymbol{M}^{-1}$ fully maps $\forall\boldsymbol{x}\neq\boldsymbol{0}$ to $\forall\boldsymbol{\bar{x}}\neq\boldsymbol{0}$ using the relationship $\boldsymbol{\bar{x}}=\boldsymbol{M}^{-1}\boldsymbol{x}$. That is to say, $\forall\boldsymbol{x}\neq\boldsymbol{0}$ yields $\forall\boldsymbol{\bar{x}}\neq\boldsymbol{0}$. As a result, (69) is equivalent to

To summarize, $\boldsymbol{\bar{Q}}$ is a real symmetric matrix and (74) holds. Then, according to the definition of being positive definite, $\boldsymbol{\bar{Q}}$ is a positive definite matrix. Consequently, ($\rm P1$) is strictly convex, thereby ensuring the existence of a unique global optimal solution.$\hfill\blacksquare$

If $\boldsymbol{\bar{x}}^{*}$ is the global optimal solution of ($\rm P1$), then $\boldsymbol{\bar{x}}^{*}$ satisfies the Karush-Kuhn-Tucker (KKT) condition of ($\rm P1$):

where $\boldsymbol{\lambda}_{L}$ is the vector of lagrangian multipliers. Note that $\boldsymbol{M}$ is invertible, hence $\boldsymbol{M}^{\top}\!(\boldsymbol{Q}\boldsymbol{M}\boldsymbol{\bar{x}}^{*}+\boldsymbol{c}+\boldsymbol{\bar{\mathbb{A}}}^{\top}\!\boldsymbol{\boldsymbol{\lambda}_{L}})=0$ in (75) is equivalent to:

Substituting $\boldsymbol{\bar{x}}^{*}=\boldsymbol{M}^{-1}\boldsymbol{x}^{*}$ into (75) and (76) leads to

Noticeably, (77) are exactly the KKT conditions of ($\rm P2$). Most importantly, $\boldsymbol{x}^{*}$ satisfies this KKT conditions, thereby making $\boldsymbol{x}^{*}$ the global optimal solution of ($\rm P2$).

Conversely, if $\boldsymbol{x}^{*}$ is the global optimal solution of ($\rm P2$), then $\boldsymbol{x}^{*}$ meets the KKT conditions of ($\rm P2$), i.e., (77). Substituting $\boldsymbol{x}^{*}=\boldsymbol{M}\boldsymbol{\bar{x}}^{*}$ into (77) generates:

Since $\boldsymbol{M}$ is invertible, multiplying the left and right sides of the first equation in (78) by $\boldsymbol{M}$ yields

Replacing the first equation in (78) with (79) results in the KKT condition of ($\rm P1$), i.e., (75). Clearly, $\boldsymbol{\bar{x}}^{*}$ fulfills these KKT conditions, thus ensuring $\boldsymbol{\bar{x}}^{*}$ to be the global optimal solution of ($\rm P1$).

The above proves that if $\boldsymbol{\bar{x}}^{*}$ is the global optimal solution of ($\rm P1$), then $\boldsymbol{x}^{*}=\boldsymbol{M}\boldsymbol{\bar{x}}^{*}$ is the global optimal solution of ($\rm P2$), and vice versa. Note that ($\rm P0$) and ($\rm P2$) are equivalent. Consequently, Proposition 2 holds. $\hfill\blacksquare$

The above two propositions guarantee the TE technique’s adaptability to quadratic programming.

According to (30), the coupled coefficient matrix $\boldsymbol{U}_{n}$ is actually formed by $\boldsymbol{\Phi}_{n}$. The latter, defined by (25), consists of the mapping coefficients that map the dispatchable generations in region $n$ to the system states of the whole grid. Therefore, calculating $\boldsymbol{U}_{n}$ corresponds to determining the mapping coefficients in $\boldsymbol{\Phi}_{n}$.

To this end, the DLPF model of the whole grid is needed. This model is given by

where $\boldsymbol{V}_{n}$ consists of the voltage magnitudes of the $PQ$ nodes in region $n$, while $\boldsymbol{\widetilde{Q}}_{n}$ is formed by the reactive power injections at the same nodes. Besides, $\boldsymbol{\theta}_{n}$ is composed of the voltage angles of the $PQ$ and $PV$ nodes in region $n$, while $\boldsymbol{\widetilde{P}}_{n}$ includes the active power injections at the same nodes. In addition, $\boldsymbol{G}_{Q,n}$ denotes the line conductances between the $PQ$ nodes in region $n$ and all the other nodes in the grid, while $\boldsymbol{G}_{P,n}$ includes the corresponding line susceptances. Apparently, ISO $n$ only knows $\boldsymbol{G}_{Q,n}$ and $\boldsymbol{G}_{P,n}$ in (95), as the related lines are either the inner lines within region $n$ or the tie lines between region $n$ and other regions.

We further define $\boldsymbol{G}$ as:

Then

where $\boldsymbol{G}^{-1}$ is denoted by

Clearly, once ISO $n$ obtains $\boldsymbol{B}_{Q,n}$ and $\boldsymbol{B}_{P,n}$, it can directly form the desired matrix $\boldsymbol{\Phi}_{n}$. Computing $\boldsymbol{B}_{Q,n}$ and $\boldsymbol{B}_{P,n}$ requires the inverse calculation of $\boldsymbol{G}$. However, ISO $n$ only has partial information of $\boldsymbol{G}$, preventing ISO $n$ from calculating the inverse of $\boldsymbol{G}$. Actually, calculating the inverse of $\boldsymbol{G}$ can be considered as solving the following system of linear equations:

where $\boldsymbol{I}$ is an identity matrix with the appropriate dimension, and $\boldsymbol{B}_{Q,1}\ \cdots\boldsymbol{B}_{P,N_{a}}$ are the solution of this system. Note that the linear equation system in (100) has a special feature: ISO $n$ only knows some equations of the system, i.e., ISO $n$ only has $\boldsymbol{G}_{Q,n}$ and $\boldsymbol{G}_{P,n}$ at hand. This is a typical problem in the field of distributed computing, which has drawn much attention recently [35]. However, the state-of-the-art methods to solve this problem usually require a lot of iterations, leading to low computational efficiency. Therefore, we propose a fast distributed method that can solve (100) in a confidentiality-preserving way.

First, ISO $n$ randomly generates two invertible and private matrices: $\boldsymbol{W}_{Q,n}$ and $\boldsymbol{W}_{P,n}$. These two matrices are used to adapt $\boldsymbol{B}_{Q,n}$ and $\boldsymbol{B}_{P,n}$ into their encrypted versions, i.e., $\boldsymbol{\bar{B}}_{Q,n}$ and $\boldsymbol{\bar{B}}_{P,n}$:

Further, we define $\boldsymbol{W}$ as:

Then the following equation holds:

Substituting (108) into (100) generates

Since only ISO $n$ knows $\boldsymbol{G}_{Q,n}$, $\boldsymbol{G}_{P,n}$, $\boldsymbol{W}_{Q,n}$, and $\boldsymbol{W}_{Q,n}$, ISO $n$ can freely share the masked information $\boldsymbol{G}_{Q,n}^{\top}\boldsymbol{W}_{Q,n}$ and $\boldsymbol{G}_{P,n}^{\top}\boldsymbol{W}_{P,n}$ with other ISOs, as none of the other ISOs can deduce $\boldsymbol{G}_{Q,n}$ and $\boldsymbol{G}_{P,n}$ from the exchanged information. Once all the masked information is received from the other ISOs, each ISO can solve (113) and obtain the encrypted solution. Eventually, ISO $n$ can recover $\boldsymbol{B}_{Q,n}$ and $\boldsymbol{B}_{P,n}$ from the encrypted solution using (101). Consequently, the whole process is confidentiality-preserving.

It is worth mentioning that in the above process, only the sharing of masked information requires communications between ISOs — the actual calculations can be implemented independently by each ISO. To make the whole process a distributed procedure, i.e. not all ISOs having to communicate with all other ISOs, an accelerated average consensus (AAC) algorithm [36] is leveraged to adapt the information sharing process into a fully distributed calculation.

The AAC algorithm is a graph-theory-based distributed method, where in the considered case ISOs are the vertices in the graph while communication lines between ISOs constitute the edges. During the iterations of this algorithm, each ISO only needs to communicate with its one-hop neighbors, yet the calculation result of every ISO will converge to the mean value of all ISOs’ initial values. Hence, to realize the information sharing between ISOs by the AAC algorithm, ISO $n$ should first form its initial value $\boldsymbol{y}_{n}(0)$ using its encrypted information:

where all elements are zero except for the ($2n-1$)-th element $\boldsymbol{G}_{Q,n}^{\top}\boldsymbol{W}_{Q,n}$ and the $2n$-th element $\boldsymbol{G}_{P,n}^{\top}\boldsymbol{W}_{P,n}$. Then, ISO $n$’s initial value is updated according to:

and $k$ is the iteration counter. All the parameters above, including $\beta$, $\varpi_{n,n}$, $\varpi_{n,m}$, $\eta_{1}$, and $\eta_{2}$, can be easily set using the generalized rules mentioned in [36]. With this process, the iterative result $\boldsymbol{y}_{n}(k)$ of ISO $n$, will converge, i.e.,

Multiplying the converged result by $N_{a}$, ISO $n$ finally obtains