Flexibility Requirement when Tracking Renewable Power Fluctuation with Peer-to-Peer Energy Sharing

$\mathcal{I}$

Set of consumers.

$\mathcal{J}$

Set of prosumers.

$l\in\mathcal{L}$

Line $l$ in set $\mathcal{L}$.

$f_{k}(.)$

Disutility function of user $k\in\mathcal{I}\cup\mathcal{J}$.

$I$

Number of consumers.

$J$

Number of prosumers.

$d_{k}$

Contract demand of user $k\in\mathcal{I}\cup\mathcal{J}$.

$w_{k}$

Forecast renewable output of prosumer $k\in\mathcal{J}$.

$\Delta w_{k}$

Real-time renewable output deviation of prosumer $k\in\mathcal{J}$.

$\tau$

Parameter in user’s objective function.

$\pi_{kl}$

Line flow distribution factors.

$F_{l}$

Power flow limit of line $l\in\mathcal{L}$.

$\alpha_{k},\beta_{k},\zeta_{k}$

Parameters of the disutility function $f_{k}(.)$.

$\underline{D}_{k},\overline{D}_{k}$

Lower/Upper bound of demand adjustable range.

$\Delta d_{k}$

Demand adjustment of user $k\in\mathcal{I}\cup\mathcal{J}$.

$q_{k}^{c}$

Actual sharing amount of user $k\in\mathcal{I}\cup\mathcal{J}$.

$q_{k}$

Expected sharing amount of user $k\in\mathcal{I}\cup\mathcal{J}$.

$\delta_{k}$

The gap between actual and expected sharing amount of user $k\in\mathcal{I}\cup\mathcal{J}$.

$r_{k}^{d}$

Flexibility requirement of user $k\in\mathcal{I}\cup\mathcal{J}$.

The structure of the proposed energy sharing mechanism is shown in Fig. 1. First, each user enters its private information into the smart meter it links to, including its disutility function $f_{k}(\Delta d_{k})$, the contract demand $d_{k}$, the demand adjustable range $[\underline{D}_{k},\overline{D}_{k}]$, and the predicted output of renewable generators $w_{k}$ if it is a prosumer. Then the smart meter will determine the optimal bid $q_{k}$ by solving (2) or (3) and submit the bid to the operator. Upon receiving all the bids, the operator decides on the optimal sharing schedule to meet users’ expectation as much as possible while satisfying the physical constraints by solving (1). Then the obtained optimal schedule $q_{k}^{c}$ and unfulfilled quantity $\delta_{k}:=q_{k}^{c}-q_{k}$ are returned to the smart meter. The smart meter will adjust its bid and submit it to the operator again. This happens until the optimal schedule $q_{k}^{c},\forall k\in\mathcal{I}\cup\mathcal{J}$ changes little between the last two iterations. Finally, the users will execute the sharing schedule $q_{k}^{c},\forall k\in\mathcal{I}\cup\mathcal{J}$. If $q_{k}^{c}>0$, user $k$ is a buyer; if $q_{k}^{c}<0$, user $k$ is a seller; if $q_{k}^{c}=0$, it does not take part in sharing. The key of the energy sharing mechanism design lies in how each user chooses the optimal bid $q_{k}$ and how the operator decides on the optimal sharing schedule $q^{c}_{k},\forall k\in\mathcal{I}\cup\mathcal{J}$.

For the operator, given users’ bids $q_{k},\forall k\in\mathcal{I}\cup\mathcal{J}$, it solves problem (1) to set the optimal sharing schedule $q_{k}^{c},\forall k\in\mathcal{I}\cup\mathcal{J}$:

The objective function (1a) minimizes the gap between the actual and expected sharing amount. Constraint (1c) means the energy sold equals the energy bought. The exchanged energy is transmitted via the power network limited by constraint (1d).

The smart meter decides on its optimal bid $q_{k}$ by solving

For consumer $k\in\mathcal{I}$:

For prosumer $k\in\mathcal{J}$:

The objective of each user $k\in\mathcal{I}\cup\mathcal{J}$ consists of two parts: its disutility $f_{k}(\Delta d_{k})$ and the gap between its bid $q_{k}$ and operator’s schedule $q_{k}^{c}$. Parameter $\tau$ indicates the trade-off between the above two parts and help evaluate the gap in monetary terms. Constraints include the power balance conditions for consumers (2b) and for prosumers (3b), and the demand adjustable range limits (2c) and (3c).

The proposed energy sharing mechanism well fits the practice in that the users do not need to know the network constraint (1d) available only to the operator, and the operator does not require local constraint (2c) and (3c) private to each user. Compared to centralized dispatch, our mechanism has two advantages: 1) Regarding privacy. Centralized dispatch needs to gather information from all users to make the central decision. However, the users may be unwilling to provide this sort of private information. In the proposed peer-to-peer energy sharing, each user only needs to submit a bid to the operator with its private data input only to its own smart meter, so privacy can be preserved. 2) Less communication burden. Under peer-to-peer energy sharing, the only information exchanged will be the bid, which is a scalar. However, under centralized dispatch, the operator needs to collect information such as cost coefficients, demands, and demand adjustable boundaries, from all users, which incurs a high communication burden overhead. Meanwhile, this mechanism can achieve a social optimal outcome (proved in Proposition 1 later).

Remark 1: Our proposed energy sharing mechanism focuses on the real-time market [33], which will be conducted hourly. In real-time, when the uncertainty is realized, we resort to a mechanism similar to the economic dispatch [34] (happens every several mins or an hour) to balance the system. Specially, instead of balancing the real-time power via centralized dispatch, we develop an energy sharing mechanism that allows the users to adjust their elastic demands and exchange energy with each other. AGC that happens at a smaller time resolution and usually at the transmission level is not considered in this paper as most of the economic dispatch works [34, 4].

Remark 2: In this paper, deciding on the optimal bids and finally the sharing quantities involve solving optimization problems (2) or (3) and two-way communication with the operator. With the current technology, advanced smart meters have certain data processing abilities and enable two-way communication between the meter and the central system [35]. Therefore, it is reasonable to assume that with the assist of a smart meter, the user would be able to decide on the optimal bids. When taking into account various kinds of appliances and facilities, the decision-making may become more complex. In that case, we may need a more sophisticated energy management system to handle all the devices.

The above peer-to-peer energy sharing market can be regarded as a generalized Nash game [36] with its elements summarized in TABLE I. For simplicity, we use $\mathcal{G}=\{\mathcal{K},\mathcal{X},\Pi\}$ to denote the sharing game in an abstract form, where $\mathcal{K}$ is the set of players, $\mathcal{X}$ the action set, and $\Pi$ the collection of cost functions. It’s worth noting that, the action sets of users depend on the strategies of the operator (the $q_{k}^{c},\delta_{k},\forall k\in\mathcal{I}\cup\mathcal{J}$); the action sets of the operator depends on the strategies of the users (the $\Delta d_{k},q_{k},\forall k\in\mathcal{I}\cup\mathcal{J}$). Therefore, it constitutes a generalized Nash game, whose equilibrium is harder to analyze than the standard Nash game due to the complex coupling [37]. The equilibrium of the generalized Nash game, from which every player has no incentive to deviate, is defined below.

In the following proposition, we prove the existence of the equilibrium and offer a centralized optimization problem to compute it, which casts down to multi-parametric program parameterized in the renewable output deviations. Denote $\mathcal{Q}:=\{q^{c}~{}|~{}\mbox{\eqref{eq:operator.3} and \eqref{eq:operator.4} are met}\}$, $\mathcal{D}_{k}:=\{\Delta d_{k}~{}|~{}\underline{D}_{k}\leq d_{k}+\Delta d_{k}\leq\overline{D}_{k}\}$ for all $k\in\mathcal{I}\cup\mathcal{J}$.

The proof of Proposition 1 can be found in Appendix -A. It offers a more convenient way to calculate and analyze the energy sharing game equilibrium by solving (6). As problem (6) minimizes the total disutility of all users, $\Delta d^{*}$ is social optimal. Take the real-time renewable output deviation $\Delta w$ as parameters, problem (6) can be regarded as a quadratic multi-parametric optimization problem.

For the sake of analysis, the objective function is linearized via a convex combination approach [38]. Take a univariate convex function $z=g(\xi)$, $\xi_{l}\leq\xi\leq\xi_{u}$ for example. The function is evaluated at some discrete points $\xi_{1},\cdots,\xi_{K}$, and $z_{1}=g(\xi_{1}),\cdots,z_{K}=g(\xi_{K})$. By introducing variables $\sigma_{1},\cdots,\sigma_{K}\geq 0$ and $\sum_{k=1}^{K}\sigma_{k}=1$, $\xi$ and $g(\xi)$ in the optimization problem can be replaced with linear functions $\sum_{k=1}^{K}\sigma_{k}\xi_{k}$ and $\sum_{k=1}^{K}\sigma_{k}z_{k}$ in $\sigma$, respectively. Applying this technique to (6a), problem (6) turns into a multi-parametric linear program (MP-LP) with the compact form as:

where $x\in\mathbb{R}^{n}$ and $\theta\in\mathbb{R}^{p}$ are vectors of decision variables and parameters, respectively; $A\in\mathbb{R}^{m\times n}$, $t\in\mathbb{R}^{m}$, $B\in\mathbb{R}^{m\times p}$ and $c\in\mathbb{R}^{n}$ are input data. The decision variables $x$ is bounded due to the lower and upper bound constraints. Denote $v(\theta)$ and $x^{*}(\theta)$ as the optimal value and optimal solution associated with given parameter $\theta$, respectively.

Based on the above MP-LP (8), in the following, we develop an adaptive vertex generation (AVG) algorithm to approximate the optimal cost and the demand adjustment at equilibrium under different renewable output deviations and quantify the flexibility requirements of individual users.

Let $\bar{\Theta}$ be the maximum set of $\theta$ such that problem (8) has at least one feasible solution $x$, known as the dispatchable region [22]. We assume $\theta\in\Theta$, a bounded polyhedral subset of $\bar{\Theta}$. The flexibility requirement of each user is defined below.

The flexibility in this paper is a kind of reserve capacity [34]. It defines the backup adjustable capacity the demand can provide in the occurrence of renewable generation output deviations due to uncertainty. Since $\Delta d^{*}(\theta)$ involves solving problem (6), getting $[\underline{r}_{k}^{d},\overline{r}_{k}^{d}]$ directly by definition can be difficult. In this paper, an MP-LP based approach is proposed to obtain the flexibility requirement. To be specific, we divide $\Theta$ into non-overlapping critical regions CR₁, $\cdots$, CR_N such that in each region CR_i, the optimal value $v(\theta)$ and optimal solution $x^{*}(\theta)$ are linear functions in the parameter vector $\theta$. Then (9)-(10) is equivalent to finding the maximum/minimum point of a piecewise linear function, which is easier to solve.

The majority of current MP-LP related works identifies the critical regions by leveraging the graph of optimal bases based on which the expression of $x^{\star}(\theta)$ and $v(\theta)$ are obtained [28]. However, degeneracy that indicates the existence of multiple primal or dual optimal solutions may cause difficulties in building the critical regions [39]. In this paper, we develop an adaptive vertex generation algorithm to approximate $v(\theta)$ via dual variables without the information of critical regions first, and then retrieve the critical region and $x^{*}(\theta)$ accordingly. In this way, the impact of degeneracy can be avoided. The details of the algorithm are as follows:

1) Lower bound of $v(\theta)$ from dual problem

Suppose the optimal value function $v(\theta)$ can be represented by a piecewise linear function

where $m_{1},\cdots,m_{N}$ are constant scalars; $n_{1},\cdots,n_{N}$ are constant vectors. For any fixed $\theta$, the dual problem of (8) is

where $\gamma$ is the dual variable, and its feasible region is denoted by $\Gamma=\left\{\gamma\middle|A^{\top}\gamma=c,\gamma\leq 0\right\}$. According to strong duality, the optimal value of (12) is bounded and equals to $v(\theta)$ for any given $\theta$. Since $v(\theta)$ is finite, the dual optimal solution $\gamma$ must belong to the set of extreme points of $\Gamma$ [vert($\Gamma$) for short], although $\Gamma$ may be unbounded and contain extreme rays. In this regard, by enumerating the vertices of $\Gamma$, $v(\theta)$ can be equivalently expressed as

Comparing (13) and (11), coefficients $m_{i}$ and $n_{i}$ and dual variable $\gamma_{i}$ have the following relation

Nevertheless, vertex enumeration is an exhaustive task, and only a small fraction of vertices correspond to valid pieces in (11) and (13), but we do not know which $\mbox{one in vert}(\Gamma)$ appears in (13) in advance. Let $\underline{\Gamma}$ be a subset of $\mbox{vert}(\Gamma),$ and

Hence, $\underline{v}(\theta)$ is an underestimator of $v(\theta)$. The approximation quality depends on the choice of $\underline{\Gamma}$. In the following, an adaptive vertex generation algorithm is developed to construct $\underline{\Gamma}$ so that the gap between $\underline{v}(\theta)$ and $v(\theta)$ can be small.

2) Recover the critical regions

Suppose (15) is a non-redundant representation. Intuitively, if the value of any piece $\gamma^{\top}_{i}t+\gamma^{\top}_{i}B\theta$ can reach maximum for some $\theta\in\Theta$, (15) is non-redundant. The exact definition and redundancy elimination methods are in Appendix -B. Then the critical region associated with each piece can be retrieved from (15). Recall the original piecewise linear form (11) and the point-wise maximum form (15), the value of $m_{i}+n_{i}^{\top}\theta$ achieves maximum in CR_i, i.e.:

where $[-i]=\{1,\cdots,i-1,i+1,\cdots,n\}$ stands for the set of complementary indexes of $i$, so CR_i is described by $n-1$ inequalities in (16). It can be represented in a matrix form (17), where $H\in\mathbb{R}^{(n-1)\times p}$, $h\in\mathbb{R}^{p}$.

For notation conciseness, index $i$ of critical region for matrix $H$ and vector $h$ are omitted; $H_{j}\theta\leq h_{j}$, $j=1,\cdots,(n-1)$ represent the individual constraints in CR_i.

3) Overall algorithm

Suppose we already have an approximation (15) for the optimal value function, as well as critical regions retrieved from (16), then we have to investigate the quality of this approximation in each CR_i. The task comes down to evaluate the maximal approximation error through

where $v(\theta)$ is the true optimal value function of problem (8), which is unknown. $\varepsilon_{i}$ must be nonnegative because $\underline{v}(\theta)$ is an underestimator. First, due to strong duality, we have

Problem (19) is nonconvex due to the product term $\gamma^{\top}B\theta$ in the objective function. Since the parameter space has a low dimension, we can solve (19) by numerating the vertex of $\mbox{CR}_{i}=\{\theta|H\theta\leq h\}$, which is

Another approach to solve (19) is the big-M method [40], but our approach is more robust since it avoids the difficulty of choosing a proper $M$. The proposed adaptive vertex generation algorithm is summarized in Algorithm 1. The output is a $\varepsilon$-optimal underestimator of the true optimal value function.

Remark: Because the set $\mbox{vert}(\Gamma)$ has finite elements, Algorithm 1 must terminate in a finite number of steps. Nonetheless, as only a small fraction of elements in $\mbox{vert}(\Gamma)$ corresponds to the pieces in the optimal value function, Algorithm 1 is efficient in practice. In contrast to the existing approaches where critical regions are obtained from the analysis of optimal bases invariancy, in Algorithm 1, critical regions are retrieved from the optimal value function. Therefore, model degeneracy does not affected the proposed algorithm.

To obtain the flexibility requirement of individual users, we need to recover an approximate optimizer $x(\theta)$ as a function of $\theta$. In CR_i, the optimal value is $v(\theta)=\gamma^{\top}_{i}t+\gamma^{\top}_{i}B\theta$, the complementarity and slackness condition of (8) implies

Let $A^{\prime}_{i}/B^{\prime}_{i}/t^{\prime}_{i}$ denote the rows of $A/B/t$ corresponding to the indexes of nonzero elements in $\gamma_{i}^{*}$, we obtain a set of linear equations:

If $A_{i}^{\prime}$ is invertible, the optimizer is obtained by

Otherwise, we can calculate the flexibility requirement $[\underline{r}_{ik}^{d},~{}\overline{r}_{ik}^{d}]$ by the following optimization problems. Without loss of generality, suppose the first $I+J$ terms in $x$ correspond to the decision variables $\Delta d_{k},\forall k=1,...,(I+J)$.

and

Therefore, the flexibility requirement of user $k$ over the entire parameter set $\theta$ is $[\underline{r}_{k}^{d},\overline{r}_{k}^{d}]$, where

In general, the obtained critical regions CR${}_{i},\forall i$ can provide geometry information of how uncertainty be handled under energy sharing, while the flexibility requirement $[\underline{r}_{k}^{d},\overline{r}_{k}^{d}],\forall k$ offers a quantitative metric.

The topology of the 5-bus system is given in Fig. 3 with other parameters in TABLE II-2. There are three users with elastic demands and two renewable generators at nodes C and E, respectively, so $\theta=[\Delta w_{1},\Delta w_{2}]^{\top}$. The maximum set of parameters $\bar{\Theta}$ is the dispatchable region [22], the grey region in Fig. 4. The $\theta$ can vary within the rectangle subregion $\Theta$. Applying Algorithm 1, we can obtain six critical regions as in Fig. 5. Note that the quadratic objective function (6) is linearized by 5 segments. The approximate optimal cost function $\underline{v}(\theta)$ in each critical region is as follows.

Following the procedure in Section III-C, we can get the demand adjustment at equilibrium in each critical region as well. Take $\mbox{CR}_{4}$ as an example, we have $\forall[\Delta w_{1},\Delta w_{2}]^{\top}\in\mbox{CR}_{4}$:

If we choose a scenario inside $\mbox{CR}_{4}$, say $[\Delta w_{1},\Delta w_{2}]=[-10,-20]^{\top}$kW, then according to (27), the approximate optimal cost is $770.96 ($=837.26-2.01\times 10-2.31\times 20$); according to (IV-A), the approximate demand adjustments at equilibrium are $\Delta d_{1}^{*}=11.10$ kW, $\Delta d_{2}^{*}=-20$ kW, $\Delta d_{3}^{*}=-26.10$ kW. Both the optimal cost and demand adjustments are the same as the one we get by solving problem (8) directly with given parameters $\theta=[\Delta w_{1},\Delta w_{2}]=[-10,-20]^{\top}$kW. This validates the effectiveness of our algorithm. The optimal cost of the original quadratic problem (6) is $767.24, so the relative error is only 0.48%, showing the linearization is accurate enough. Moreover, the obtained $\Delta d_{k}^{*},\forall k=1,2,3$ are the same as well. We calculate the flexibility requirement of each user, shown in TABLE IV. To accommodate the volatile renewable generation, all users need to have the ability to both increase and decrease their demand. User 3 requires the most flexibility since the value of $\overline{r}_{3}^{d}-\underline{r}_{3}^{d}$ is the largest.

Given $[\Delta w_{1},\Delta w_{2}]=[-10,-20]^{\top}$kW, the equilibrium of the energy sharing game (1)-(1) can be reached via a best-response algorithm. The change of $\Delta d$ and $\delta$ during iterations are shown in Fig.6 and Fig. 7, respectively. We can find that the demand adjustments converge to the above approximate values. The value of dual variables of problem (6) at optimum can be retrieved by $\gamma_{i}$ in $\mbox{CR}_{i}$, i.e. $\eta_{1}^{*}=-1.92$, $\eta_{2}^{*}=-2.08$, $\eta_{3}^{*}=-2.31$, which are very close to the value of $-\delta/\tau$ plotted in Fig. 7. The above analysis validates Proposition 1.

A larger case with a modified 69-bus microgrid [42] is tested with topology in Fig. 8. There are three renewable generators connected to node 9, 30, and 60, respectively. Suppose their real-time output may deviate within $[-30,30]$ kW. Applying Algorithm 1, the change of the maximum approximation error (defined in (18)) during iterations is shown in Fig. 9. The error decreases fast and the algorithm can output 10 critical regions (as in Fig. 10) after three iterations. The approximate optimal cost function $\underline{v}(\theta)$ in each critical region is

To examine the accuracy of the approximation, let’s take $\mbox{CR}_{4}$ as an example. Following the procedure in Section III-C, the demand adjustments at equilibrium are:

For scenario $\theta=[-30,-30,30]$ kW $\in\mbox{CR}_{4}$, according to (29) and (IV-B), the optimal cost is $373.29 and $\Delta d_{1}^{*}=-17$ kW, $\Delta d_{2}^{*}=-9$ kW, $\Delta d_{3}^{*}=-4$ kW, $\Delta d_{4}^{*}=8$ kW, $\Delta d_{5}^{*}=5$ kW, $\Delta d_{6}^{*}=-24$ kW, which are the same as the one obtained by solving problem (6). Similarly, we can calculate the flexibility requirement for each elastic demand as in TABLE V. We notice that while elastic demands 2, 3, 4 still have some redundancy, the adjustable ability of elastic demands 1, 5, 6 has been fully exploited as $d_{k}+\underline{r}_{k}^{d}=\underline{D}_{k}$ and $d_{k}+\overline{r}_{k}^{d}=\overline{D}_{k}$, for $k=1,5,6$. Therefore, in long-term planning, we may put more emphasis on increasing the flexibility of demands 1, 5, 6.

To show the applicability and scalability of the proposed algorithm, we test its performance under different numbers of renewable generators and users in TABLE VI; and under different settings in TABLE VII. The impact of the number of users on the computational time is marginal. The time needed increases with the number of RGs, but in this paper, we only consider the distributed renewable generators with a relatively large scale (say about 500kW) while the uncertainties from small-scale units are neglected. Usually, in a residential community, there are few such relatively large-scale units so our proposed method can still be applied.

To test the practicability of the proposed mechanism, we compare the computational time under different numbers of users in TABLE VIII. The computational time does not change much when the number of users increases since the decision-making for each user can be done in parallel. Our algorithm can handle a large number of users.

Problem (1) can be equivalently written as

Suppose $(\Delta d^{*},q^{*},q^{c*},\delta^{*})$ is an GNE of the sharing game, then $q^{c*}$ is the optimal solution of (-A). Therefore, we have

Also, $(\Delta d^{*},q^{*})$ is the optimal solution of problem (2) or (3), which is equivalent to

where $\mathcal{D}_{k}:=\{\Delta d_{k}~{}|~{}\underline{D}_{k}\leq d_{k}+\Delta d_{k}\leq\overline{D}_{k}\}$. Therefore,

Similarly, let $(\Delta\hat{d},\hat{q}^{c})$ be the optimal solution of problem (6) and $\hat{\eta}$ the value of the dual variable. Then, we have $\forall(\Delta d,q^{c},\eta)\in\prod_{k\in\mathcal{I}\cup\mathcal{J}}\mathcal{D}_{k}\times\mathcal{Q}\times\mathbb{R}^{(I+J)}$: