A receding horizon data-driven chance-constrained approach for energy ﬂexibility trading in multi-microgrid distribution network

Unexpected imbalances between supply and demand during real-time operation in the wake of high penetration of renewable energy resources necessitate system operators to improve ﬂexibility of the system. Moreover, increasing growth in ﬂexible distributed resources installation in distribution networks with multiple MicroGrids (MGs) such as fuel-based distributed generators, storage units, and controllable loads makes distribution system operators be able to manage energy imbalances locally through exploiting potentials of aforementioned resources in ﬂexibility enhancement. This paper proposes a receding horizon data-driven chance-constrained energy ﬂexibility trading framework for active distribution networks with multiple MGs. Speciﬁcally, a distance-based conﬁdence set is constructed in a data-driven manner to apply net load uncertainty into this modelling. According to this proposed data-driven approach, the energy ﬂexibility is modelled using ﬂexibility envelope notion to demonstrate its dynamics in different scenarios in an imminent horizon. Furthermore, MGs’ ﬂexibility trading problem is modelled as an exact potential game that possesses a pure strategy Nash Equilibrium (NE). An iterative decentralized algorithm based on alternating direction method of multipliers is developed to achieve the NE strategy of MGs. Numerical results are also represented to show the performance efﬁciency of the proposed energy ﬂexibility trading framework.


INTRODUCTION
Traditionally, distribution networks were highly passive with uncontrollable loads and distributed energy resources that could only work based on the fit-and-forget approach [1]. This fact made them unable to work independently in balancing supply and demand during real-time operation from the upstream network or the balancing market. Under such circumstances, Distribution System Operators (DSOs) most likely face high market prices, leading to a substantial increase in their total cost. However, in Active Distribution Networks (ADNs), the DSO can reliably and independently operate the distribution network. Despite passive distribution networks, DSO in an ADN is able to make use of controllable distributed energy resources such as can support local balancing in real-time by utilizing the great potential of the mentioned controllable resources in providing energy flexibility as a reaction to the uncertainty of its net load [2,3]. Generally, energy flexibility in power system is defined as the ability of components of the power system in changing the level of energy that they generate or consume in response to an incentive signal. The incentive signal is sent to agents due to encountering expected as well as unexpected changes in power system condition. These changes can be variations in net load or any equipment contingencies. Accordingly, for a reliable operation, improving the flexibility of the power system regarding the rapid growth of renewable energy resources penetration is unavoidable [4,5]. Furthermore, the great need for reliable, economic, and environmentally friendly energy supply leads to increasing growth in the formation of MicroGrids (MGs).
In such environments in which information privacy is an issue of vital importance, decentralized market mechanisms can act more efficiently than conventional centralized ones [6,7].
Additionally, centralized energy market mechanisms cannot provide energy that procured through flexible behaviour of resources such as FLs [4]. Consequently, the emergence of decentralized market mechanisms that facilitate energy flexibility transaction is becoming unavoidable. From this aspect, DSO can manage energy balance through trading with MGs rather than real-time balancing market when MGs experience energy shortage or surplus [8]. Specifically, this paper focuses on local management of unexpected energy imbalances by utilizing the flexible behaviour of distributed resources in ADNs with multiple MGs. Several previous works in the scope of flexibility management in power systems can be found in the literature. An investigation on flexibility reserve in day-ahead operation was presented in [9] in which a stochastic optimization was introduced for optimizing both energy and flexibility reserve. This paper aims to effectively plan reserve capacity to handle energy imbalances and ramping requirements of uncertain net load in real-time. In [10], authors developed a receding horizon stochastic control framework for ADNs. In this framework, controllable resources were utilized to ensure the ADN's reliable operation in case of facing uncertainties that arise from loads and renewable generations. Authors in [11] developed an optimization framework for managing flexibility services that are procured from distributed energy resources. In the mentioned study, a two-stage stochastic programming model was performed to handle day-ahead uncertainties and intra-day fluctuations, aiming at maintaining balance and solving optimal power flow problem in the ADN. In [12], a framework for distributed energy flexibility management and optimization in distribution networks was performed. In this study, DSO offers aggregated distributed flexibility procured by FLs at the day-ahead energy market. [13] studied day-ahead energy and reserve market clearing under wind generation penetration aiming at balancing energy using ramp-capability reserves modelling. By applying a two-stage optimization approach, energy flexibility was procured by controllable loads in [14]. Here, home energy management systems were utilized to optimize controllable load energy consumption patterns. In [15], authors applied a multi-stage distribution expansion planning in which two different regulatory policies are proposed to incentivize investment in flexible resources such as dispatchable DGs, energy storage systems, and demand response programs. In [4], monopolistic and game-based approaches were implemented to manage flexibility transaction in distribution network. In this study, the flexibility provided by FLs is transacted among the end-users, aggregators and DSO. A fully decentralized energy trading framework based on a non-hierarchical analytical target cascading algorithm was applied in [6] to manage energy trading in ADNs among DSO and MGs. The literature contains several studies investigating uncertainty modelling in power systems as well as those that ignore the probabilistic nature of wind generation and loads which leads to deterministic optimization. Stochastic programming, robust optimization, and chance-constrained programming are highly popular approaches in uncertainty modelling in power system optimization problems [16][17][18]. In [19], a stochastic bilevel model was presented for energy management of ADNs with multiple MGs. In the mentioned paper, the authors introduced a decision-making framework for microgrids and DSO in case of facing uncertainties caused by wind energy and load in real-time. Authors in [16] introduced a stochastic mixedinteger programming model for unit commitment problem in MGs. This predictive-stochastic energy management framework addresses the uncertainty of intermittent energy resources based on a two-stage decision-making mechanism. A study on day-ahead scheduling and control of networks with microgrids was performed in [20]. This study proposed a stochastic programming model to manage uncertainty of different sources. Moreover, the authors considered risk-neutral and risk-averse models in their proposed optimization framework. In stochastic programming, uncertain parameters such as wind generation are assumed to follow a particular probability distribution. In other words, the true probability distribution is supposed to be known which is rarely the case in practice. The probability distribution usually is determined by curve fitting over the historical data due to the limited number of historical data [21,22]. In robust optimization, the objective function often comprises minimization of worst-case cost which results in a relatively conservative solution. Despite stochastic programming, robust optimization can be implemented regardless of the true probability distribution of random parameter [22,23]. In chance-constrained programming, a chance constraint is satisfied within a specified confidence level. Furthermore, the solution is not as conservative as robust optimization solution [24]. In addition to the above common approaches, the envelope concept can be used to facilitate uncertainty management. This way, the realization of random parameters can be restricted by envelope segments. Moreover, in flexibility management applications, potential dynamics of power system resources can also be modelled by envelopes [25,26]. A flexibility envelope model can be found in [27] which represents flexibility management of energy-limited resources as deployable energy reserve using energy-centric perspective.
This paper proposes a receding horizon data-driven chanceconstrained approach for energy flexibility trading in an ADN. We characterize the energy flexibility of resources using flexibility envelopes to demonstrate their dynamics during the prediction horizon. By using flexibility envelope enclosures, output trajectories of flexible resources are restricted in different scenarios. Besides, the energy flexibility trading problem is modelled as an exact potential game. As discussed, several previous works investigated flexibility management in power systems. However, there is a distinct gap in locally energy flexibility provision in distribution networks with multiple MGs equipped with energy-constrained resources such SUs and FLs. The contributions of this paper are as follows: 1) We propose a novel multi-interval receding horizon energy flexibility trading framework in an ADN with multiple MGs in which DSO locally trades energy with MGs to accommodate net load variations. In our proposed framework, MGs optimize their strategy based on updated net load forecast data at any time step which brings out net load uncertainty reduction. Accordingly, flexibility provision cost reduces by limiting the consequences of net load forecast error. With regard to multi-interval modelling, we characterize FLs as shiftable loads that enables them to optimize their consumption during the prediction horizon with constant total energy demand. 2) We propose a novel data-driven chance-constrained programming model to capture and manage the net load's uncertainty. Specifically, we construct a distance-based confidence set to apply net load uncertainty using available historical data. Then, we reformulate chance constraint into its equivalent closed-form. Furthermore, we characterize energy flexibility of resources by flexibility envelope concept to demonstrate their dynamics in various scenarios during the prediction horizon. From this aspect, we schedule energy flexibility for times ahead, aiming at managing the net load uncertainty. Using flexibility envelope enclosures, we limit output trajectories of distributed flexible resources in different scenarios within a predefined confidence level. 3) We prove that MGs' energy trading problem is a noncooperative game named as exact potential game. Afterward, we show that the proposed potential game possesses a pure strategy Nash equilibrium (NE). In order to achieve the NE point, an iterative decentralized algorithm based on the Alternating Direction Method of Multipliers (ADMM) is developed.
The remainder of this paper is organized as follows: Section 2 contains a system model that presents MGs' constraints and objective function as well as distribution network constraints. A data-driven chance-constrained approach is applied in Section 3 to model uncertainty of the net load. In this section, chanced constraint is reformulated to its equivalent deterministic form. Section 4 performs MGs energy flexibility trading problem. This section contains the proposed decentralized game-based algorithm. Performance evaluation of the proposed framework using several case study simulation results are demonstrated in Section 5 and the whole paper is concluded in Section 6.

SYSTEM MODEL
In this study, we consider an ADN with multiple MGs. MGs are equipped with distributed energy resources such as fuel-based DGs, FLs, SUs, wind turbines, and also uncontrollable loads in our modelling. Each MG can act as a load or power supply that is connected to a node in the ADN. Hence, an ADN that comprises several MGs which can trade energy with each other is seen as a multi-microgrid distribution network [6]. There is a DSO in an ADN that manages energy trading to achieve an energy balance between supply and consumption in our proposed flexibility trading mechanism. Suppose there are M MGs as players in the potential game. The trading framework scheme is illustrated in Figure 1. According to this figure, MGs use historical and forecasted net load data to make their best decision in the energy flexibility provision mechanism.
With regard to multi-interval receding horizon modelling, net load forecast data are updated successively. Moreover, MGs optimize their strategy based on the performed data-driven chance-constrained approach in Section 3 and send it to DSO according to the developed iterative algorithm in Section 4. It is worth noting that t demonstrate the current time slot and  = {1, … , T − t } is set of ahead time slots in the prediction horizon. At any time step, MGs choose their best strategy with respect to updated wind generation forecasts. The goal is to limit the consequences of net load forecast error. In the next subsections, we present the energy flexibility envelope notion and MG's trading problem formulation.

Flexibility envelope notion
Conventionally, flexibility is modelled as constant capacity in the operational horizon in flexibility provision problems. Considering significant growth in wind generation penetration in power systems, the aforementioned state of the art studies in the scope of flexibility management are not able to satisfy the flexibility requirement of the system [26]. This paper characterizes energy flexibility by the flexibility envelope concept. Hence, the energy flexibility can be scheduled for times ahead, aiming at managing the net load uncertainty. Moreover, we assume that resources' energy flexibility follows their respective interval in hours ahead in the simulation horizon to facilitate the net load uncertainty accommodation. From this aspect, flexibility cannot be modelled as a contracted constant capacity. A flexibility envelope can be constructed by bounding flexible resources feasible trajectories or realizations of a random parameter with respect to their allowed ramp rate. In this way, maneuvers of energy output trajectories of distributed flexible resources in different scenarios are limited by using flexibility envelope enclosures in our proposed framework. Accordingly, we designed power and energy-based envelopes for DGs, FLs, and SUs. Figure 2 represents their corresponding flexibility envelope enclosures. Further details can be found in the following subsections.

Distributed generations
Distributed fuel-based resources such as diesel engines can play an important role in improving distribution networks' flexibility by their ability to change their scheduled output energy. The cost of providing flexibility by DG g in MG i at time t ∈  is presented in (1). f u i,g and f d i,g are upward and downward segments of flexibility envelope of DGs. Here, we have c u i,g and c d i,g as their corresponding cost, respectively.

FIGURE 2 Illustration of envelope enclosures
The flexibility envelope for DGs is illustrated in Figure 2(a). As can be seen, flexibility segments of DG's output are bounded by upward e f u i,g and downward e f d i,g flexibility envelope segments in prediction horizon. From this aspect, constraints (2) and (3) are imposed to DGs. Here, P sch i,g denotes scheduled output power of DG g ∈  in MG i ∈ . Also maximum and minimum output powers of resource g are, respectively, denoted by P max i,g and P min i,g . For instance, a violation in constraint (3) can lead to an infeasible trajectory as demonstrated in Figure 2 Based on Figure 2(a), we have several trajectories ofP i,g due to the net load uncertainty during prediction horizon. These several trajectories are related to the possible realization of the net load. Obviously, feasible trajectories ofP i,g in different scenarios are bounded. This way, (6) is given to constrain DG's output power trajectory in several scenarios. Here,P i,g represents output power trajectory of DG g ∈  in MG i ∈ .
Moreover, changes of the output power of DG g ∈  considering its provided flexibility in two consecutive time slots are restricted by its maximum ramp rate, ΔP i,g , as indicated in (7).

Flexible loads
Implementation of demand response programs in distribution networks has profound effects on power system flexibility enhancement [28] due to the ability of flexible loads in handling both positive and negative imbalances [29]. This paper introduces a mechanism in which shiftable loads as energy flexibility providers can improve system flexibility with regard to their capability in shifting a portion of their demand to another time within the daily horizon. For shiftable loads, it is assumed that their whole energy demand is constant during the day. The cost of providing flexibility for FLs for i ∈  and t ∈  is demonstrated in (8). Upward and downward flexibility segments of FLs are defined by f u i,l and f d i,l and also, c u i,l and c d i,l are their corresponding costs, respectively.
. (8) Dynamics of FLs flexibility in hours ahead of time period are exhibited in Figure 2 and downward e f d i,l flexibility envelope segments. Consequently, constraints (9) and (10) must be considered in FLs formulation. Here, D sch i,l is the scheduled demand of flexible end-users. w i,l and y i,l are percentile of shiftable demands out of the whole scheduled demand which are determined by flexible end-users and can take a value between 0 and 1.
Changes in upward and downward segments of FLs energy flexibility are confined in (11) and (12).
are maximum flexibility ramp up and ramp down of FL l for i ∈  and l ∈ . As can be seen in Figure 2(b), maneuvers of feasible flexibility trajectoriesL i,l in different scenarios are restricted. Consequently, constraint (13) is imposed to FLs.
The pivotal point inferred from Figure 2(b) is that the surface area betweenL i,l and D sch i,l is equal to zero. This is to say that FLs are modelled as shiftable loads, which means that the whole demand of these loads remains constant during the day. In order to take into account the mentioned point in problem modelling, Equations (14) and (15) are given. L f i,l is shifted demand in time slot t ∈  and k ∈ . L f i,l takes positive value when controllable end-user reduces its consumption and takes negative value once controllable end-user increases its consumption. Regarding the implementation of the multi-step receding horizon framework, the first term in (14) represents shifted demand in previous time slots. The second term represents a portion of demand with shift capability to hours ahead in the prediction horizon. It should be noted that only the second term is defined as a variable. For infeasible trajectories shown in Figure 2

Storage units
SUs can participate in flexibility procurement mechanism due to their capability in changing their state from charging to discharge, and vice versa [30,31]. So, batteries act as load in the charging state and as supplier in the discharging state. Equation (16) demonstrates SUs cost owing to providing energy flexibility for i ∈ , t ∈  and k ∈ . In this equation, SUs upward and downward flexibility segments are represented by f u i,s and f d i,s , respectively. Moreover, c u i,s and c d i,s are their corresponding costs. It should be noted that storage is in the charging state in upward trajectory. Likewise, discharging energy leads to a downward trajectory.
According to (17) SUs energy flexibility envelope is demonstrated in Figure 2(c). In this figure, upward and downward flexibility envelope segments are denoted by e f u i,s and e f d i,s . Similar to DGs and FLs, flexibility envelope is constructed by bounding SUs energy flexibility trajectories in several scenarios. As can be seen, flexibility segments of SUs are restricted by envelope segments. Hence, we apply this restriction to SUs by (18) and (19). The right hand sides of (18) and (19) represent the remaining capacity of storage to be charged or discharged in order to manage uncertainty of the net load. E max i,s shows the maximum energy capacity of storage s. Furthermore, maximum and minimum states of charge of storage s in MG i are respectively denoted by SoC max i,s and SoC min i,s . Upward and downward segments of flexibility in two consecutive time slots are also limited through (20) and (21), where, Δ f u i,s and Δ f d i,s are defined as maximum flexibility ramp up and ramp down of storage s ∈ , respectively.
The need for net load variation management necessitates us to characterize SUs energy flexibility using several trajectories. In Figure 2(c), feasible energy flexibility trajectories for SUs are illustrated. In what follows, energy flexibility trajectory,Ẽ i,s , in different scenarios is restricted by upward and downward flexibility segments.
In (23), E f i,s denotes final stored energy level of SU s considering upward and downward segments of flexibility envelope. Constraint (24) limits final state of charge of storage s in each time slot. Eventually, (25) confines changes in final stored energy level of storage by its maximum ramp rate ΔE i,s . It is worth noting that, as we have considered SU's upward and downward flexibility costs in the objective function, there's no need to keep charging and discharging variables mutually exclusive by using binary variable [32].

Total balancing cost of microgrids
Total balancing cost of MG i ∈  due to procuring energy flexibility through trading with DSO in addition to its flexible resources is indicated in (26). C i,tr is MG's cost due to trading energy with DSO which is defined in (27). Traded energy flexibility between MG i and DSO is denoted by P i2DSO . Here, i,tr is the trade price with DSO. It must be noted that P i2DSO takes positive value when MG i is selling its surplus energy flexibility to DSO. Likewise, P i2DSO takes negative value once MG i faces energy flexibility shortage.

Deterministic power balance constraint
This study demonstrates how MGs balance their generation and consumption through an energy balance constraint given in (28). In the right-hand side of this equation, the forecasted net load of MG i ∈  is denoted by NL i . With regard to the implemented receding horizon framework, net load forecasts for hours ahead in prediction horizon are updated in time step

Chance-constrained power balance
MGs are likely to confront unexpected imbalances between generation and demand as a result of net load uncertainty [33]. Consequently, they need to be capable of adjusting either their generation or consumption level according to changes in the net load with respect to its forecasted value. As mentioned before, MGs manage uncertain net load, D i, , through the exploitation of flexible capacity of their flexible resources in addition to trading energy as much as their surplus or shortage with DSO. It is worth noting that, scheduled demand of flexible end-users, D sch i,l , is supposed to have no uncertainty. Hence, its changes relative to scheduled value are aiming to enhance the flexibility of the distribution network. Accordingly, the uncertainty of MG's forecasted net load NL i is captured by constructing a confidence set using historical net load data in Section 3. Chanced power balance constraint for MG i ∈  is demonstrated in (29). Here, is defined as the uncertain parameter which has a finite support. In other words, there are a finite number of possible realizations (scenarios 1 , 2 , … , N ) for MG's net load [22]. Chance constraint (29) implies that the probability of constraint violation is less than the risk level̄for all possible realizations of the random net load. The confidence level is denoted by 1 −̄, and is allowed power imbalance tolerance.

Distribution network
Evidently, one of the most important roles of DSO is to manage the supply-demand balance in the distribution network. In our proposed framework, DSO manage supply-demand balance in time slot t ∈  and k ∈  in the whole distribution network by imposing (30). Note that, in our modelling, system loss has been ignored. Distribution system constraints modelling can be found in [6,8,34].

DATA-DRIVEN CHANCE-CONSTRAINED SOLUTION APPROACH
Contrary to stochastic programming approaches, in our proposed data-driven chance-constrained solution framework, the true probability distribution of random parameters is assumed to be unknown. In other words, we assume that the net load as the only random parameter in our proposed framework does not follow any particular distribution. Hence, we design a confidence set for the ambiguous true probability distribution of the net load using available historical data samples [22]. Afterward, we reformulate chance constraint (29) into its equivalent deterministic form to be solvable by standard solvers.

Confidence set construction
According to the aforementioned discussion, true realizations of random net load are likely to be different from its forecasted value. Consequently, uncertainty modelling of MGs' net load is unavoidable. From this aspect, we have utilized a distance-based confidence set which is constructed based on historical net load data. Here, we have a reference distribution which is aimed to converge to true probability distribution of random parameter as the number of observations increases [23]. Here, we construct a histogram to estimate reference distribution. Other approaches such as utilizing Kernel Density functions can also be used [24]. Assume there are A data samples for MG's net load. The number of scenarios represents the number of histogram bins which is denoted by N . Furthermore, the number of data samples in each bin are A 1 , A 2 , … , A N with A = ∑ N n=1 A n . Accordingly, the estimated reference distribution for uncertain net load is represented by P re f = (p 1 re f , p 2 re f , … , p N re f ) in which, p 1 re f = A 1 ∕A, p 2 re f = A 2 ∕A, … , and p N re f = A N ∕A are reference probabilities. The designed confidence set under L 1 norm is given below [22] Here, P is ambiguous true distribution, P re f is estimated reference distribution for random net load, and  represents the set of all possible probability distributions. In addition, the tolerance level is denoted by , which takes value based on the number of historical data and the given confidence level. As the number of net load data increases, the estimated reference distribution gets closer to the true distribution function which brings out the divergence tolerance level reduction. Obviously, the divergence tolerance level converges to zero when data numbers go to infinity. This fact can also be inferred from (34). For the distance between two distributions, √ P re f ) 2 ) 1 2 = 2Θ H (P, P re f ), in which Θ H (P, P re f ) denotes Hellinger distance between two distributions [35]. It has been proved in [36] that 4AΘ 2 H (P, P re f ) converges to a 2 N −1 random variable in distribution, as the number of historical data A increases. In (34) The degrees of freedom is represented by N − 1. Based on (34), we can obtain the tolerance level. It is worth noting that the distance between reference and true distributions can be measured by other methods such as KL divergence [37] and Φ-divergence

Chance constraint reformulation
Here, we reformulate chanced power balance constraint (29) using the approach given in [23] to get its identical closed-form to be solvable by common solvers such as CPLEX. In order to consider the worst-case distribution for net load, chance constraint (29) is rewritten in (35) for i ∈ , t ∈  and k ∈ . In other words, we have considered the largest distance between true and reference distributions in designed ambiguous set  1 to ensure that the chanced power balance constraint is satisfied.
Reformulation of constraint (35) with regard to big-M approach is demonstrated in (36) and (37). Here, z n i is binary variable. According to constructed confidence set  1 definition in (32) we have (38)- (40).
We have reformulated (38) as (41)-(43) by introducing an auxiliary variable y n i for i ∈  to remove absolute operation, where, To eliminate minimization operation in constraint (37), the duality theorem is implemented. Hence, equations (44)-(47) represent the dual problem of objective function (37) and its associated constraints (39)- (43). Note that z n i is considered as constant parameter in (37). Here, corresponding dual variables for i ∈  are denoted by i , n i , n i , and i .
As discussed before, MGs are likely to confront unexpected imbalances between generation and demand due to net load uncertainty. In these circumstances, they can sell or buy energy as much as their surplus or shortage. Therefore, the final formulation of MGs' energy flexibility trading problem with DSO can be summarized in P1. P1: MG energy flexibility trading problem

DECENTRALIZED GAME-BASED SOLUTION METHODOLOGY
In this section, we develop a decentralized algorithm based on ADMM to solve interdependent MGs' optimization problem with respect to their information privacy. Additionally, we model MGs' optimization problem as a non-cooperative game named as potential game. Potential game modelling ensures both overall and individual optimality [38].

Potential game approach with ADMM
As we discussed in Section 3, MGs' energy trading problems are interdependent. From this aspect, the optimization problem can be modelled as a game. In game theory, games are divided into two non-cooperative and cooperative categories. There are different types of non-cooperative games, such as potential game. This paper models the MGs energy trading problem as an exact potential game. Specifically, potential game approach implementation facilitates investigation of NE strategy existence [39,40] and leads to the global optimality. What follows is the potential game definition.

Definition 1. A non-cooperative game in which a single global function
can express all game players' motivation to deviate from their strategy is a potential game. The aforementioned global function represents potential function of the game which is of great importance in existence investigation of the NE strategy [39,41].
In the formulation, strategy profile of players are represented by s = (s i , s −i ) andŝ = (ŝ i , s −i ), where s i andŝ i are two different strategies for player i. Moreover, strategy profile of all game players except for MG i is represented by s −i . We have In (49), Φ(s) is exact potential function. This equation implies that each player's gain from a deviation in its strategy set is equal to the gain in potential function due to an equivalent deviation [41]. U i (s) is MG's payoff function defined as in P2. The vector of the strategy profile of MG i is defined as follows: Interdependency between MGs problem can be inferred from Equation (30). In order to solve the energy flexibility trading problem in a decentralized fashion, we need to reformulate P1 using one of several common decentralized methods such as ADMM. This method has been utilized for optimal scheduling of multi-microgrid ADNs [8,42]. [8,43].

P2: MG energy flexibility trading problem with ADMM
sub ject to (2) Here, the penalty factor is denoted by . represents the Lagrange multiplier. The MG's energy flexibility trading problem which is demonstrated in P2, fulfill exact potential game condition (49) with our proposed potential function given in (51). The MG's payoff function U i (s i , s −i ) is demonstrated in P2's objective function. Also, the proof is performed in Proof 1.
Proof 1:Regarding the exact potential game definition given in Definition 1, proof of MGs trading problem as an exact potential game is given by (52)-(54).
According to (52) and (53) we have, Proof 2: In the potential games category, the exact potential game is considered as a special form of ordinal potential game [38]. In this way, we apply important properties of ordinal potential games in order to investigate the existence of the pure NE strategy. In what follows, some important properties of ordinal potential games are given [41].
• Suppose Φ as an ordinal potential function for strategic game Γ(U 1 , … , U M ). In this strategic game, U i represents the payoff function for each game player. Then, the equilibrium set of strategic game Γ(U 1 , … , U M ) coincides with the equilibrium set of Γ(Φ, … , Φ). Accordingly, if Φ admits a minimal value in the game players' strategy profile, then Γ as a finite ordinal potential game possesses a pure strategy equilibrium. • Every finite ordinal potential game has the finite improvement property.
Based on the above-mentioned first feature, the NE point of exact potential game P2 corresponds to the potential function's minimum point (51). NE point can be achieved using our proposed iterative algorithm. Since every maximal improvement path in a finite ordinal potential game with finite improvement property terminates in an equilibrium point, the convergence and optimality are guaranteed.

Decentralized solution
In this section, we develop an iterative decentralized algorithm based on ADMM. Using our proposed iterative algorithm, the NE point of the given exact potential game in P2 can be achieved. Algorithm 1 illustrates our proposed receding horizon iterative approach. In the first step, both initial value for t , ALGORITHM 1 1: Set: t = 1, T = 24, = 10 −3 .

5:
Each MG receives its net load forecasts for time t + k to T .

7:
DSO and MG i initialize the Lagrange multiplier 0 and traded energy flexibility P 0 i2DSO with DSO, respectively. 8: Repeat

9:
Each MG sends its strategy at t ∈  and k ∈  to DSO based on (55).
which is defined as the current time step, and the stopping criterion are determined. Note that T − t is assumed as the prediction horizon. In line 3, the prediction horizon time counter, k, is set. Then, MG i receives its net load forecast data for hours ahead. In this algorithm, x denotes the iteration number. In line 6, the initial value for x is determined. Afterward, DSO and MG i initialize the Lagrange multiplier and traded energy flexibility with DSO, respectively. Then, MG i sends its best strategy in time slot t ∈  and k ∈  to DSO by solving (55) with respect to P2's constraints. Based on (55), MGs optimization problems are solved sequentially. In iteration (x + 1) th , MG i optimize its strategy s i considering MG r best strategy s r for r ∈  and r > i at iteration x th and MG m best strategy s m for m ∈  and m < i at iteration (x + 1) th .
Next, DSO updates according to (56). The decentralized algorithm's execution stops when the difference between two consecutive Lagrange multipliers becomes less than the stopping criterion. By doing so, the energy balance in the distribution network is satisfied. The procedure mentioned above will continue until the end of the simulation horizon T . Since the objective function in P1 is convex, closed, and proper, convergence conditions hold for the proposed iterative decentralized algorithm. The mathematical proof can be found in [43].
For better illustration, our proposed multi-interval receding horizon decentralized algorithm is summarized in Figure 3.

PERFORMANCE EVALUATION
In this section, the proposed methodology's effectiveness is evaluated through inclusive case study simulations and result discussions. The problem simulation has been implemented in MATLAB and GAMS environments and solved by the CPLEX solver on a personal computer with processor 2.30-GHz intel Core i5 CPU and 6 GB of RAM memory. There are three MGs as game players in ADN. Each MG is equipped with DGs, FLs, SUs, wind turbines, and uncontrollable loads. The day is divided into 24 equal time slots. Hence, the simulation horizon is 24 h of a day. Energy flexibility provision cost for MGs as well as its corresponding trading price in addition to simulation parameters can be found in Table 1. The uncertain net load D i, (t + k) data samples are generated randomly and the histogram as the reference distribution is constructed using MATLAB. It is worth noting that data samples can be achieved either by using Monte Carlo simulation [23] or gathering real net load datasets. The number of net load data samples are assumed to be A = 1000. Furthermore, the number of scenarios that denotes the number of histogram bins is N = 5. As we discussed before, the estimated distribution contains partial information about the true probability distribution. The confidence level 1 − ′ is supposed to be equal to 95% for each of three MGs which means constructed confidence set  1 contains true probability distribution of MG's net load within 95% confidence level.
The divergence distance between true and reference distributions for 1000 data samples and 95% confidence level is equal to 0.097. We have assumed that the controllable loads' demand in each MG is equal to 25% of each MG's whole net load. That is to say, we have ∑ l ∈ D sch i,l (t + k) = 0.25NL i (t + k). Net load forecast data, NL i (t + k), for MGs during the entire period is illustrated in Figure 4. This paper investigates four different cases. In Case 1, we evaluate the performance of the receding horizon data-driven energy flexibility trading framework under the aforementioned data and parameters assumption. In Case 2, we examine the effect of different confidence levels on total cost and robustness of the model. In Case 3, we vary the number of net load data samples to examine our proposed framework performance in cost reduction. In Case 4, the data-driven energy flexibility trading framework is implemented excluding the receding horizon approach, and in Case 5, the results for energy flexibility trading problem regardless of data-driven chance-constrained approach is represented.

Case 1: Implementation of receding horizon energy flexibility trading framework
In this case, MGs trade energy flexibility with DSO when facing energy shortage or surplus due to variations in net load by reason of wind generation uncertainty. The trading framework is based on the flowchart illustrated in Figure 3. Figure 5 illustrates the convergence of our proposed algorithm for time slot 14. According to this figure, our proposed framework converges after 14 iterations. Based on this figure, MG 3 sells its surplus energy. MG 1 and MG 2 buy energy flexibility as much as their energy shortage. The sum of traded energy flexibility after 14 th iteration is equal to zero, which means the energy balance in the distribution network is fulfilled. Solution times for MGs 1 to 3 for single interval optimization at the end of the day are equal to 0.35, 0.32, and 0.31 s for each iteration, respectively. At the beginning of the day, for multi-interval optimization, these times increase to 2.48, 2.16, and 2.08 s for MGs 1 to 3, respectively. Regarding the fact that our proposed trading mechanism is suggested for short-term operation horizon, this mechanism is quite applicable for real scale systems. To investigate our proposed framework effectiveness in terms of energy flexibility requirement satisfaction for MGs, we have illustrated flexible resources' output, traded energy flexibility, and flexibility requirement for MGs in Figure 6. Note that the energy flexibility requirement for MGs in time slot t is procured when flexible resources output and traded energy flexibility cover reference net load probability distribution within at least the predefined 95% confidence level. In this figure, provided upward and downward energy flexibilities of DGs, FLs, and SUs and traded energy flexibilities of MGs are represented. Besides, the flexibility requirements are plotted using dashed lines. For traded energy flexibility, we have dashed areas and white undashed areas. The dashed area represents surplus upward or downward energy generation of MG. MG faces surplus upward or downward energy generation when unexpected net load reduction or unexpected net load increase occurs, respectively. In these circumstances, MG sells its surplus upward or downward energy generation as flexibility to DSO. Similarly, white areas imply upward or downward energy shortage of MG which occurs when MG experiences unexpected net load increase or unexpected net load reduction, respectively. Hence, MG has to buy energy flexibility as much as needed. This figure implies energy balance in MG.

Case 2: Investigating the effect of different confidence levels
In the previous case, we considered the confidence level 1 − ′ to be equal to 95% for each of three MGs. In this case,  we investigate the effect of different confidence levels on the total cost of energy flexibility procurement in the distribution network. Considering fixed data sample's number A = 1000, we vary ′ from 0.01 to 0.4. In this way, the percentile of the true probability distribution of MG's net load coverage by constructed ambiguity set varies between 60% to 99%. According to (34), changes in confidence level will result in changes in distance between two distributions. Numerical results are depicted in Table 2. According to the results, system cost increases as the confidence level increases. The reason is that a higher confidence level leads to larger confidence set construction which can be inferred from changes in value in Table 2. In this experiment, we also present an out-of-sample assessment in order to study the effect of different confidence levels on problem robustness. The robustness is assessed through determining out-of-sample reliability indexes Expected Energy Not Supplied sub ject to (6), (13), (22) and In this problem, MGs' active power balance equations are relaxed by introducing non-negative slack variables S u i,d (for load shed) and S d i,d (for wind generation curtailment). Here, we generated 5000 data samples for the uncertain net load. The upward and downward flexibility segments procured by flexible resources as well as energy flexibility traded between MGs and DSO are set to their associated values which are obtained by solving the proposed energy flexibility trading framework given in Algorithm 1. The numerical results are represented in Table 3. According to this table, higher confidence levels result in power balance infeasibility reduction. From this aspect, the confidence level can be adjusted according to a trade-off between cost and robustness.

Case 3: Investigating the effect of different number of data samples
The number of data samples A is an important factor in ambiguity set construction. Regarding the fact that the distance between true and reference distributions depends on the number of random parameter samples, different data numbers lead to different reference distributions. Undoubtedly, the more data samples we have, the more information we can get from the net load's ambiguous true probability distribution. In other words, with increasing the number of observations, the reference distribution gets closer to the true probability distribution. In this case, the data number changes between 700 to 5000 and the important parameters are the same as Case 1. The simulation results are demonstrated in Table 4 for 95% confidence level. The total solution time for 24-interval optimization (at the beginning of the day) for convergence is provided in this table. It is evident, the distance between two distributions ( ) decreases as the number of net load data increases. Hence, the problem becomes less conservative. Although this fact leads to a total cost reduction, total solution time as well as the computational complexity increases.

Case 4: Implementation of the energy flexibility trading framework regardless of receding horizon
Regarding the implementation of receding horizon fashion, MGs are able to choose their best strategy with respect to updated wind generation forecasts at any time step. The main purpose is to limit the consequences of net load forecast error, leading to a more reliable energy supply with lower cost. In this case, we have implemented our proposed framework without considering the multi-interval receding horizon concept. Regardless of multi-interval simulation, FLs cannot optimize their consumption during a day by shifting a part of it in time to participate in the flexibility provision mechanism. Accordingly, each FL decides to curtail a part of its demand or increase its consumption at any time slot. Changes in FLs' scheduled demand in 24 h of a day as the simulation horizon is represented in Figure 7. As can be seen in Figure 7(a,c), in multi-interval receding horizon simulation, MG 1 and MG 3 have participated in the energy flexibility provision mechanism by shifting their  demand in time with respect to constant total energy demand in the time period. In other words, the sum of their shifted demand during the day is equal to zero. This point can be easily inferred from Figure 7(a,c). Moreover, in Figure 7(b,d), both MGs 1 and 3 FLs are modelled as curtailable loads. In the sequel, their total consumption during the day may not remain constant. Note that negative changes mean demand reduction and positive changes denote an increase in demand. Another important point of receding horizon fashion implementation is its impact on cost reduction. The numerical results are represented in Table 5. As can be seen, both MGs and overall flexibility provision costs in the distribution network have been increased compared to Case 1. Case 5: Implementation of the energy flexibility trading framework regardless of data-driven chance-constrained approach In order to show the effectiveness of the data-driven chance-constrained trading approach, we compared the result of the Case 1 with a case in which we ignored implementation of our proposed data-driven chance-constrained solution methodology described in Section 3. In this way, we assumed that the net load follows a Gaussian distribution and the flexibility requirement is achieved based on this distribution. The result is demonstrated in Tables 6 and 7 for 95% confidence level. According to the demonstrated results in Table 6, the total cost has been increased about 5% compares to the Case 1 in which the trading is done using the proposed data-driven chance-constrained approach. This is because the results are achieved in a more conservative way. This can be seen in Table 7. According to this table, EENS and ERC, on average, have been decreased 37% and 50% for MGs, respectively.

CONCLUSION
Here, we applied a receding horizon data-driven chanceconstrained energy flexibility trading framework. In the proposed mechanism, DSO trades energy flexibility with MGs in an ADN to accommodate net load variations by exploiting distributed resources' flexible behaviour. We also modelled MGs' energy flexibility trading problem as an exact potential game that possesses a pure strategy NE. According to the results, the proposed decentralized energy trading algorithm converged after a finite number of iterations. Besides, the effectiveness of the envelope-based data-driven model in managing and capturing the net load uncertainty was demonstrated in terms of flexibility requirement satisfaction. Specifically, we investigated the effect of different confidence levels and number of data samples on the proposed framework. It was shown that there is a trade-off between cost and robustness of the problem. In fact, these two cases demonstrated how the conservativeness of the trading problem can be controlled by quantifying the number of given net load data samples and confidence level. Additionally, the results of energy flexibility trading framework implementation regardless of receding horizon concept and the proposed data-driven chance-constrained approach showed that the proposed framework led to about 9% and 5% system cost reduction, respectively.