Robust Scheduling for Networked Microgrids Under Uncertainty

This paper proposes a framework that facilitates the energy exchange between networked microgrids (MGs) in the electricity market. An alternating direction method of multipliers (ADMM)-based robust optimization algorithm is proposed to derive the optimal energy exchange strategy for the networked MGs considering the uncertainties of the electrical load, intermittent generation, and electricity prices in the external market. The proposed method naturally lends itself to a classical market-clearing problem between two hierarchical levels comprising (i) main-grid-to-MG and (ii) MG-to-MG, aiding in the result interpretation and practical realization. Leveraging from the decentralized organization, the operational autonomy and information privacy of each MG is protected. The proposed framework is tested on a modified 144-node network with 3 MGs. The numerical experiments demonstrate the convergence of the proposed market clearing to the market equilibrium in different grid operational scenarios with different conservatism parameter settings for MG operators.


INTRODUCTION
Conventionally, power systems contain a few large-scale generation facilities that supply energy to the passive loads. Lately, this arrangement has been slowly transforming into a multi-layered, cyber-physical intelligent grid with active loads and small-scale generators. These small-scale generators based on intermittent renewable sources are presently not involved in grid operations and the major concerns for rejecting their participation often include the small size and unreliability of Distributed Generatorss (DGs) (Saraiva and Gomes, 2010). By using the concept of a Microgrid (MG), the DGs and aggregated Flexible Loads (FLs) can be clustered to form larger market players and provide ancillary services (Guerrero et al., 2013). The MGs can either be connected to the main grid or be interconnected to form a networked-MG system (Figure 1), where the new decentralized structure contributes to establishing a competitive market and benefits the social welfare of the underlying system. Some practical examples include the networked MGs on the campus of the Illinois Institute of Technology and its adjacent community Bronzeville (Shahidehpour et al., 2017), where the optimal energy management mechanisms for interconnected MGs are exploited. The networked MG system is different from the traditional power distribution grid, because the power flow is bidirectional from one MG to another MG or to the distribution grid. The change in network topology can be frequent in networked MGs. If implemented and managed properly, networked MGs can provide a variety of benefits to electric power utilities and consumers in terms of supply efficiency, security, and reliability (Alam et al., 2019). In the scope of their market participation, the literature proposes the networked MGs participate in the local Distribution System Operator (DSO) market (Hirsch et al., 2018;Alam et al., 2019), wherein the networked MG system can act as virtual power plant with each MG's controller having direct control over its resources. In regard to this, MGs can be owned by utilities, end customers, or consortia of prosumers. For the latter cases that the Distributed Energy Resourcess (DERs) are not the assets of the MG, direct control of these assets might not be easily exercised in the context of a deregulated market. To this end, the indirect control over DERs can be achieved by the introduction of a virtual layers of Mirogrid Operator (MGO)s. Each MG is assumed to be managed by the local MGO that is an independent entity to collect the local prosumers' bids from FLs, and DGs and to coordinate with the market of the external power system to provide energy and ancillary services. With the networked MGs geographically scattered over multiple communities, the network constraints of the networked MGs can be respected at the same time.

Related Work and Contributions
For the market integration of networked MGs under the transactive energy paradigm (Rahimi et al., 2016), Sabounchi and Wei (2017) propose a block-chain based mechanism to enable the electricity exchange for energy prosumers in an networked MG system environment, where actual grid operation challenges for networked MGs remain undiscussed. Various studies Wu et al., 2019) propose a bi-level-programming model for the MGs' bidding strategy in the distribution grid market based on the assumption that the prosumers have full knowledge of the distribution grid information, such as topology and loading profile. This may be unrealistic for some countries where the grid information is considered to be confidential (see, e.g., Trpovski et al., 2018). To this end, market design and frameworks for networked MGs can be developed by considering the statistical distribution of available data, such as the market-clearing prices and its underlying uncertainty (e.g., Kim and Poor, 2011). Various approaches including stochastic programming, chance-constrained models, robust optimization, and distributionally robust optimization exist in the literature to develop the optimal exchange framework for the networked MG operation in coordination with the main grid.
Chance-constrained models are employed to tackle the uncertainties for the networked MG scheduling problem in works (Bazmohammadi et al., 2018(Bazmohammadi et al., , 2019aDaneshvar et al., 2020a), where Gaussian distribution and beta distribution are assumed to describe the uncertain parameters for electricity prices and Photovoltaics (PV) generation to relax the chance constraints while ensuring that the problem is numerically tractable. From a practical perspective, the exact distribution of uncertain variables is difficult to obtain for given empirical data as we later demonstrate in a case study for electricity prices and load profile, where machine learning-based approaches may deliver a better estimation for these uncertainties. To remove this drawback, the sampling technique can be adopted to transform the chance-constrained model to stochastic programming problems, where it provides an alternative way to handle uncertainties by generating discrete probabilistic scenarios of the uncertain parameter and integrating them into the optimization problems (see, e.g., Weron, 2005;Daneshvar et al., 2020b). Though this approach does not rely on specific types of distribution, both the solution quality and computational efficiency can be highly impacted by the scenarios that are utilized in the problem formulation. Similar to stochastic programming, another approach that does not rely on the exact knowledge of the distribution of uncertain parameters is robust optimization. Instead of optimizing the expected objective value as in stochastic programming, robust optimization essentially provides a conservative solution by considering the worstcase scenario realization of the stochastic parameters that are modeled with continuous uncertainty sets. This, on the one hand, may require less computational power, but, on the other hand, leads to a conservative solution by considering the worstcase scenario. It is worth pointing out that robust optimization based solution methodologies allow the tuning of conservatism parameters to reflect the attitude of energy prosumers. Some examples are shown in  and Correa-Florez et al. (2020), where robust optimization for load aggregation in a DSO market and DG generation for single MG operation are proposed, respectively.
In the center of the solution methodology design for networked MGs operation, an important aspect is the operational autonomy. Networked MGs naturally have a decentralized operation structure, where the ownership of the MGs can vary, being utility-owned, customer-owned, or third-party operated MGs (Alam et al., 2019). Centralized approaches for managing the networked MGs require each MG operator to reveal their complete local system information to a central entity, which may introduce both privacy and scalability concerns in the system. A large body of existing literature has focused on a centralized dispatch model (e.g., Bazmohammadi et al., 2019a;Daneshvar et al., 2020a;among others). In contrast, with a decentralized operation structure, competition can be introduced between the MGs to establish a competitive market at the distributiongrid level. A recent exposition (Liu et al., 2020) uses a bilevel model for the market-clearing purpose between DSO and networked MGs, where bi-level optimization models introduce two hierarchical levels to regulate the market-clearing procedure. In this work, we focus on the interconnected MGs exchange problem, wherein the distribution network is clustered into MGs with a plain operation structure for the MGOs. For this structure, decomposition methods can be applied (Zhang et al., 2020). There are some works considering the decomposition method for the networked MG management problem. For example, Ma et al. (2018) proposes an online Alternating Direction Method of Multipliers (ADMM)-based algorithm to handle real-time dispatch problems for networked MGs. Work by Liu et al. (2017) proposes the decomposition methods for deterministic day-ahead networked MGs scheduling. Although these works concentrate on the decomposable operation structure to preserve information privacy, the uncertainty handling for day-ahead trading within networked MGs has not been addressed yet.
In view of the existing works and research gap, we focus in this work on the day-ahead energy exchange problem for networked MGs considering the grid constraints and the uncertain parameters preserving a decentralized Peer-topeer (P2P) operation structure. In terms of essential solution methodology choices, we use the combination of the ADMMbased decomposition method and robust optimization. The reasons are 3-fold. First, ADMM is well-suited to distributed optimization and large-scale distributed computing systems. By using ADMM, we do not need to have a centralized computation center and the computation can be distributed across the MG operators. Specifically, ADMM consists of several inner loop subroutines that can be assigned to each MG operator. Each subroutine is an optimization problem that can be solved in a parallel and independent fashion. Second, ADMM can be interpreted as a form of tâtonement process, which essentially adjusts the prices to achieve the market equilibrium, following the Walras theory of general equilibrium (Boyd et al., 2011). This plays an essential part in the market design for the networked MGs. Third, the combination of robust optimization and ADMM allows individual parameter settings by each MGO to set their conservatism, avoiding an over-conservative solution based on their preference. Moreover, full operational autonomy is ensured with only information regarding the selling/buying quantity exchanged. Despite the different estimation of uncertainties and conservatism parameter settings, given the solution feasibility, a unique market equilibrium can be attained.
The main contributions of this work are summarized as follows: • To enable fully decentralized operation, we propose an ADMM-structured robust optimization method for the networked MG energy exchange problem, while maintaining the autonomy and robustness of the decentralized grid organization considering local MG contingencies (e.g., congestion between MGs). The combination of robust optimization and ADMM allows individual parameter settings by each MGO, ensuring full operational autonomy with only information regarding the selling/buying quantity exchanged. It also enables individual parameter settings of each MGO to set their conservatism, avoiding an over-conservative solution. • The market equilibrium is derived for the energy exchange prices between MGs under different grid conditions, whereas the heterogeneous preferences for parameter settings for MGOs in robust optimization are exploited to analyze the impact on the equilibrium point.
The paper is structured as follows. In section 2, we provide the modeling choices for DER assets, networked MGs, grid models including voltage and line flows, and the uncertainty models. In section 3, both the deterministic and robust form for the networked MG energy exchange problem is described. The ADMM-structured robust optimization solution methodology is presented in section 4. In particular, we analyze the market equilibrium using the exchange price decomposition in section 4.2. Finally, the simulation results and case studies for different grid conditions and parameter settings are presented in section 5. Section 6 summarizes this work with observations and outlook.

Networked MG System
In this work, we focus on the networked MG system, where each MG is assumed to be operated by a MGO. MGOs are supposed to be independent entities for coordinating, controlling, and monitoring the local MG system. They obtain buying and selling bids from the local prosumers (small-scale generators, FLs) and can coordinate with other MGOs in the same networked MG system. The networked MG system can also exchange electricity through Power Common Coupling (PCC) with the main grid. The collective objective of the networked MG system is to minimize the total MG operational cost. First, local prosumers submit their bids to the respective MGOs. Second, the welfare maximization problem of networked MGs is solved by all MGOs using the proposed distributed algorithm with information exchange regarding proxy variables of power exchange quantity.
The energy exchange schedule and market-clearing prices of the networked MG are determined during the market clearing.

Distributed Energy Resource
We assume a generic energy prosumer model for DERs that manages the optimal set points of its assets, i.e., generators, energy storage, and loads. To this end, the DERs are categorized into (i) DGs with active/reactive power injection capabilities and (ii) FLs with flexible active power consumption capabilities. All market participants are assumed to be economically rational, hence they seek the maximization of their individual economic surplus. The prosumer nodes are labeled as DGs and FLs and described by set G : = {1, 2, . . . , g} for the DGs and F : = {1, 2, . . . , fl} for the FLs. The set T : = {1, 2, . . . , T} is defined to include all the market intervals.

DER Model
For DERs, the electric power injection or consumption is limited by the system constraints, e.g., the curtailment capacity of a commercial building is limited by the comfort of its occupants; the reactive power injection of a DG is limited by the inverter capacity. Such system constraints are expressed with as box constraints: DERs within a MG are assumed to form their bids based on the cost function and the system constraints and submit the bids to the local MGO. MGs are assumed to be equipped with Battery Energy Storage Systems (BESSs). Let set B : = {1, 2, ..., b} contain the BESSs in all MGs and soc t ∈ R b denote the State of Charge (SOC) of the BESSs at time step t. Q ∈ R b is used represent the battery capacity and p d , p ch ∈ R b denote the discharge and charge power. SOC of BESSs is described with a linearized model: where η c , η d ∈ R represent the charge and discharge efficiency and t represents the market period duration. Hence, the nodal injections of BESS at time t from is given as

Cost Function
It is further assumed that the cost for each MGO to procure electric power generation/load curtailment from the DERs can be expressed with a quadratic function where p x ≥ 0 ∈ R x is the electric power injection or consumption and c x , d x are the coefficients for the linear and quadratic terms. The quadratic cost function is strictly convex with zero-crossing. The bid formulation for all the DER types is summarized in Figure 2. For example, the FL bid consists of a baseline load p fl t and a price-sensitive part [p fl t , p fl t ] with the marginal cost c fl t + diag(d fl t ) · p fl t representing their willingness to be curtailed in dependence with the market-clearing price.
For batteries, the rational behind the price sensitivity is to capture the battery aging effect from the charging and discharging process. Interested readers are encouraged to refer to Liu et al. (2017) for the derivation of the parameter settings for batteries, which gives a quadratic cost function:

Grid Model
Consider the networked MG system in Figure 1; we assume the individual MGs can be considered as a portion of a low-voltage (LV) distribution grid with a reference node indexed by 0. In the grid-connected mode, the MG is connected to the distribution grid via PCC, where the reference is modeled as a slack node with fixed voltage. The rest of the nodes are modeled as PQ nodes. All MG nodes are included in set N : = {0, 1, 2, ..., n}.
In addition, we define the set H : = {1, 2, ..., h} containing all the lines connecting the nodes. Let set MG : = {1, 2, ..., mg} denote the set of all MGs, and we obtain for MG i ∈ MG the regional version of its respective nodes set and lines set as

Linearized Nodal Injection Model
The complex nodal injections and voltages are, respectively, represented for all nodes as s : = p +  q and u : = v · e  θ all with dimension C n+1 . The individual node i ∈ N defines s i : = p i +  q i and u i = v i e  θ i . The overall grid is modeled using the admittance matrix Y : = {y ij } 1≤i,j≤n+1 ∈ C (n+1)×(n+1) with its regional version for MG i ∈ MG denoted by Y i ∈ C (n i +1)×(n i +1) . Considering only PQ (constant power) nodes exist in the grid, we have the nodal injection model at each time step t ∈ T (Bolognani and Dörfler, 2015): The nodal injection model is complex and non-linear. We adapt the first-order linearization for eq. (6) as in Bolognani and Dörfler (2015, Proposition 1). It assumes that the nodal angle differences are very small throughout the network, so the flat voltage, i.e., θ ≈ 0, v ≈ 1, is used as the linearization point. We obtain the following linearized power injection model (Linearized DistFlow Model) (Farivar et al., 2013).
Frontiers in Energy Research | www.frontiersin.org Note that the approximation based on Linearized DistFlow has been shown to perform reasonably well in voltage estimation [0.25% error for 5% voltage variation from the nominal voltage (Farivar et al., 2013)] while neglecting the higher order active/reactive loss terms. For the operation of geographically closed networked-MGs, the aim is to minimize the overall operating cost, without focusing on the losses that are small compared to the amount of energy exchange among the MGs.

Linearized Line Current Model
To constrain the line flow so as not to exceed the thermal limit, we consider the following linearized line current model. Defining the line current between node i and node j as l ij ∈ C, thermal limits are usually defined as the line current magnitude |l ij |, which is directly related to the conductor's temperature. For brevity, we obtain the following line current magnitude by assuming no existence of shunt admittance: Based on the same assumption of small angle variance, i.e., θ i − θ j = 0, we obtain FIGURE 3 | Line current model with a from/to direction.
It essentially provides a linear projection from the voltage difference to the line current magnitude. We now rewrite eq. (9) in a compact form using the admittance matrix Y line ∈ C h×h and the from/to mapping matrix C f , C t ∈ R h×(n+1) . The basic idea of representing in from/to quantities is illustrated in Figure 3, where each line segment is assigned a from/to direction from the start/end node to the end/start node. Hence, C f /C t maps each line segment to the corresponding node index (see also Zimmerman, 2020). Y line is defined as a matrix with its diagonal terms equal to the line admittance and the off-diagonal terms equal 0. We obtain the line current magnitude vector |l| ∈ R h at time step t ∈ T represented as where |Y line | ∈ R h×h extracts the magnitude of the matrix Y line . Consider the thermal constraint that |l| ≤ |l|; an equivalent box constraint at time step t ∈ T is obtained as.

Uncertainty Model
Since the networked MG system is exposed to the external electricity market in the main grid, uncertainties in the electricity prices may need to be explicitly considered. For example, the forecast for electricity prices in the Singapore real-time spot market is released 36 h prior to the actual market period. Other uncertainties are from the local MG demand and DG generations, e.g., solar PV generators. These uncertainties need to be explicitly considered and modeled to ensure the robustness of the solution. The motivation for modeling the electricity prices is to have an adequate representation of spot price trends on a weekly and monthly basis, with time resolution consistent with the underlying market periods for the prosumers to decide on their market participation preference. Since the demand and the spot prices are correlated, models for both shall be obtained simultaneously. Note that the renewable generation is neglected in the following. This is because, for the studied data set from Singapore's wholesale market, the renewable injections are still very small compared to the demand. Nevertheless, the renewable generation can be taken into account by following the same procedures as in modeling the demand and their influence on the spot prices.
We adopt a dynamic regression model with Autoregressivemoving-average (ARMA) errors to model demand and price time series. First, the ARMA model for the demand time series is given as The first three terms show the constant components of the load model, representing trends for yearly κ y t , monthly κ m t , and weekly κ w t cycles. The last three terms are the stochastic components in the ARMA model including a Gaussian-distributed term ǫ D t . Based on the demand model, the price model consists of a deterministic portion and a stochastic portion. The deterministic portion is determined by a linear regression model to represent long-term trends, e.g., seasonality, and the demand's influence on the spot prices. On the other hand, the stochastic component is determined by an ARMA model to represent the underlying evolution of the error terms. The log-transformed spot price model is described as Herein, b t is the vector of dummy variables values for the respective time step including a single variable representing the index for the linear trend. Entries of vector β are the corresponding coefficients, estimated by the regression model. Furthermore, a is the offset determined by the regression. Term δ l log(d t−l ) is the influence of the demand time series, in which coefficient δ l has been determined by the regression as well. Accordingly, all mentioned terms so far are the deterministic component of the model. In addition, the remaining terms represent the ARMA model, with φ c l and θ c l being the coefficient of the AR and MA terms, respectively. Finally, ǫ c t denotes the error terms of the ARMA model. One example is given for the generated scenarios for electricity prices in Singapore in Figure 4. Similar modeling procedures can be followed for PV generations (see also Huang et al., 2012).
Based on the demand model in eq. (12) and the spot pricing model in eq. (13), we need to obtain the lower and upper bound for these uncertain quantities. To this end, to achieve a suitable representation of the uncertainty, a high number of scenarios are needed, which are generated using Monte Carlo simulation. The upper and lower bound of the confidence interval can be used to bound the uncertain price sequences by accessing the generated scenarios, where whereγ +,− t is the maximal allowed deviation to the estimated value given a predefined confidence level.

NETWORKED MICROGRID OPERATION UNDER UNCERTAINTY
Since the MGs are interconnected, each MG has the coupled nodes shared with other MGs to exchange power. All the coupled nodes for MG i are included in set M i : = {1 i , 2 i , ...m i } ∈ R m i . Furthermore, the power exchange vectors p exc i,t , q exc i,t ∈ R m i are associated with the coupled nodes. Hence, the nodal injections p i,t , q i,t ∈ R n i +1 in MG i are the total of the electricity withdrawn from the distribution grid p pcc i , DG injections p g i,t , q g i,t ∈ R g i , energy exchange p exc i,t ∈ R m i and flexible/non-elastic loads p fl ∈ R fl i , p D , q D ∈ R n i +1 . It is obtained as where a )×fl i is the incidence matrix that maps the PCC node, DG nodes, power exchange nodes, and FL nodes to the respective local MG nodes. For an MG i, all its physically connected neighboring MGs are included in set E i ⊆ MG. In particular, the power imported by MG i from MG j ∈ E i is defined as M ij p exc j . The incidence matrix M ij ∈ R m i ×m j projects the same coupled nodes in MG j to MG i. Furthermore, the active/reactive power exchange between the MGs should hold the energy conservation rule: Frontiers in Energy Research | www.frontiersin.org

Networked MG Scheduling: Deterministic Form
The objective of each MGO is to minimize operating costs while taking into account the energy balancing, exchange, voltage stability, and congestion management. This objective is extended to the collective networked MGs based on the assumption of collaborative, non-strategic behaviors of MGOs. Let p g i , q g i,t ∈ R g i , p fl i,t ∈ R fl i denote the injected power from DGs/FLs in MG i ∈ MG, t ∈ T , giving the total cost of the MG: where γ i,t ∈ i,t which gives the total energy cost from main grid and α ⊺ i p exc i +β ⊺ i q exc i is the energy exchange payoff between MGs. Term x∈{g,fl,b} C x t (p x t ) provides the total cost for the energy and curtailment procurement from DGs/FLs/BESSs, respectively.
The centralized optimization problem for networked MG scheduling is to minimize the total operational cost subject to all MG constraints and consensus constraints eqs. (17) and (18). The deterministic version of the problem is given as: eqs. (7) and (15) to (19) where constraints eqs. (20b) and (20c) are the power balance constraints indicating that the power supply is satisfying the demand. Note that the power balance equations contain the linearized power flow equation eq. (7) implicitly. Constraints eqs. (20e) and (20g) to (20i) are the box constraints for power dispatches, voltage magnitude, and line current magnitude. Note that the box constraints eqs. (20e) and (20g) can be derived in a straightforward way by considering the capacities of the DG units and FLs. For a typical inverter-interfaced DG units, the upper and lower bounds are derived from a delimited area by a rectangle in the actual active/reactive power operational plan (Cuffe et al., 2014). For FLs, these limits can be given by estimating the maximal shiftable power in the energy consumption (Hanif et al., 2017). Constraints eq. (17) and eq. (18) ensure that the power supply from neighboring MGs meets the local import demand. For constraints eqs. (15), (20c), (20e) and (20g) to (20i), the corresponding dual variables are listed on the right side of the equation.

Robust Counterpart for Electricity Prices
To ensure the maximal benefit of the MGs under the worstcase scenario realization for the uncertain prices γ , we utilize the method of robust optimization to solve the day-ahead energy exchange problem. Robust optimization does not rely on the exact distribution of the uncertain parameters, which is generally difficult to obtain based on historical data. Instead, a continuous set in Equation (14) can be used. To further accommodate different levels of conservatism in the solution, the following set is defined to bound the uncertain prices given the conservatism parameter where 0 ≤ Ŵ p ≤ T must hold for the conservatism parameter. Intuitively, the parameter Ŵ p allows the price to reach the upper bound of the interval in no more than T time periods. The objective in eq. (19) can be reformulated for the robust optimization problem as The objective essentially captures the maximization of the cost under uncertain prices (worst-case scenario) in the inner loop and minimization of the cost with respect to the dispatch variables in the outer loop. Note that by replacing the original objective function of eq. (20) with eq. (22), the robustified form has a min-max problem structure. Now consider the term max γ i,t ∈ i,t T t=1 γ i,t p pcc i,t . To derive the robust counterpart of this term, it can be reformulated as the following optimization problem with proxy variables u i,t as max T t=1 (γ mean where u i,t ≤ 1 is used to bound the electricity price within the upper bound for each time step t. Note that only the upper bound deviationγ + i,t is considered in this formulation as the worst-case scenario can only be realized between the upper bound and the mean value (González et al., 2014). The resultant dual problem of eq. (23) is further provided as By substitution of eq. (24) into eq. (22), the min-max problem structure is replaced by a single minimization problem. To this end, the robust counterpart of the objective function in problem eq. (22) for all MGs can be obtained as eqs. (24b) and (24c), where it can be integrated directly into the problem formulation eq. (20).

Robust Counterpart for DGs and Loads
Another source of uncertainty is given by PV and electrical demand in Equation (15). Both quantities can be merged into one to obtain the net load for each market period t as (p g − p D ). Therefore, we consider the uncertain parameter only for one variable p D for the sake of simplicity. Following a similar procedure, for electricity prices the robust counterpart for the electrical load balance constraint eq. (15) can be given as with the upper bound of the load/PV generation given as p D,mean i,t +p D,+ i,t and a predefined the level of conservatism Ŵ D ∈ [0 1] for the electrical load. Note that in contrast to electricity prices, Ŵ D = 1 is the maximal value for the parameter setting, where it corresponds to the worst-case scenario realization for the combined load and DG generations.
Note that the robust optimization problem can be decomposed for each MG with individual parameter settings. Based on the market structure where each MG is operated autonomously by the local MG operator with only P2P communication, we want to design a market-clearing process and solve the above problem in a distributed manner. In the following, we introduce an ADMM-based optimization method to determine the energy exchange price/quantity α i , β i between MGOs and analyze the market equilibrium with the decomposition of the energy exchange prices.

ADMM Solution Methodologies
One of the appealing economic interpretation of ADMM is that it can be interpreted as a process working toward the market equilibrium (Boyd et al., 2011). The ingredients of ADMM are presented as follows: consider the scenario where each local MGO optimizes the local problem and obtains the results for p exc i,t , q exc i,t . Contained by their own set of local constraints, the connected MGs may not satisfy the desired import demand. Hence, a set of proxy variables are needed to temporarily store the values based on consensus between the local MGO and its neighbors. The proxy variables denoted by p exc,+ i,t , q exc,+ i,t ∈ R m i for MG i ∈ MG. Define the following augmented Lagrangian for MG i to include the consensus constraints as: where σ i,t , κ i,t ∈ R m i are the dual variables for the consensus constraints. The second-order regularization terms are included using a penalty factor ρ i ∈ R to ensure the robustness of ADMM convergence. The ADMM steps to solve eq. (28) is presented in Algorithm 1. The algorithm comprises the local maximization step, global variable update and Lagrangian multiplier update. As shown in the algorithm, the local maximization step solves the local version of the revenue maximization problem while the consensus constraints eqs. (17) and (18) are taken into account in the augmented Lagrangian. The global variable update requires the communication between interconnected neighboring MGs. Note that the consensus constraints are only defined for active/reactive power in this work. The voltage magnitude of the coupled nodes is implicitly updated from one MG to another, i.e., after obtaining the voltage magnitude as per the local maximization step in one MG, the rest of the MGs are updated with the result and use it as the reference node voltage. Since the constraints are linear in a problem eq. (20) with a quadratic objective, FIGURE 5 | Simplified schematic of the 144-node test case (Khodr et al., 2008) partitioned into three autonomous MGs.
of interest to pre-condition and change ρ in an effective manner (see, e.g., Xu et al., 2017).

Optimal Exchange Price Decomposition
Consider the robust optimization problem of MGs i in eq. (28). The associated Lagrangian of the local optimization problem is: The Karush-Kuhn-Tucker (KKT) conditions for the local optimization problem comprises the following first-order optimality condition on variable p exc i : together with primal/dual feasibility and complementary slackness (Boyd and Vandenberghe, 2004). We obtain The results make intuitive sense: The price is decomposed into three parts: energy, voltage, and congestion as in the locational marginal prices formulation (see, e.g., Zhang et al., 2020).

SIMULATION SETUP AND RESULTS
We tested the proposed ancillary service market framework on a 144-node AC distribution grid (Khodr et al., 2008), which was partitioned into three autonomous MGs. The topology of the test network is illustrated in Figure 5. Each MGO operated a number of DGs locally. One coupled node (node 7) was the energy exchange node, which connected MGs 1, 2, and 3. The test network had a total fixed nominal load of 11.9 MW and 7.38 Mvar. The cost coefficient for the DGs and BESSs are presented FIGURE 6 | ARMA model for price and load profile based on data from EMC Singapore.
in Table 1. The uncertain price and load profile is depicted in Figure 6, where an ARMA model is used to generate 1,000 scenarios for prices and loads based on historical data taken from the Nation Singapore Electricity Market EMCS from year 2014 to 2015. The upper and lower bounds for the uncertain quantities were determined as the 95% confidence interval for the reduced scenarios. For the sake of simplicity, same load shape was applied to all MGs. Voltage constraints were kept as [0.95, 1.05]; the ADMM parameters were set as: η = 0.1, τ = 0.1, κ = 10, ρ i (0) = 1.0 × 10 3 , λ i (0) = 0, i ∈ MG. The simulations were performed on a personal computer with Intel i7-8665U 3.2 Ghz and 24 GB RAM. The validation of the proposed framework was separated into the robust optimization part and the distributed optimization part.

Case Study of Robust Optimization
In the case study, each DG has a limited generation capacity such that the collective supply from all DGs cannot meet the total demand from all nodes. Therefore, the PCC that connects to MG 1 provides energy to meet the energy demand for all MGs. The first part of the case study is to investigate the robust optimization for the networked MGs, where the performance under different parameter settings for Ŵ D ∈ [0 1] and Ŵ p ∈ [0 48] are tested. To show the impact of conservatism Ŵ p , the worst-case realization of the wholesale market electricity price are obtained for different values of Ŵ p , which is depicted in Figure 7. With the increasing value of conservatism, the worst-case scenarios for the prices tend to locate at the upper bound of the uncertain price. The most conservative decision (Ŵ p = 48) to schedule the energy consumption, however, does not corresponds to the upper bound of the energy prices. This is due to the charge and discharge costs that occurs when using the battery to shift the energy consumption from the peak-price area to the low-price area. More specifically, to charge/discharge the battery at a higher power rate may result in a higher cost for the battery aging that makes this options less cost effective for obtaining the optimal  which does not necessarily increase the cost. Furthermore, a direct correlation between the relative error of the robust prices and the actual prices can be obtained for different conservatism parameter settings. Therefore, a well-chosen parameter for the conservatism setting should well represent the probability that such a worst-case scenario occurs while maintaining the robustness of the solution.

Validation of Distributed Solution
To validate the distributed solution, we solve the robust networked MG exchange problem both in a centralized fashion and distributed fashion based on proposed method. The SOCs results of the BESSs are first presented in Figure 10A for ADMM and centralized solution, respectively. Since the battery aging cost for BESS 3 is the lowest, it can be observed that BESS 3 is utilized at its full capacity to provide energy balance. All three BESSs are discharged when the spot price at the PCC and the total demand is high. One can also conclude from Figure 10A that the centralized solution fully coincides with the distributed solution.
The convergence of the distributed solution is shown in Figure 10B, where the convergence indicator of the ADMM is given as the primal residual r(k), i.e., r(k) = i∈M t∈T p exc i,t (k) − p exc,+ i,t (k) + q exc i,t (k) − q exc,+ i,t (k) . (34) In this part of the case study, the energy exchanged between all three MGs is also supplied by the PCC, where the baseline energy exchange prices are determined by PCC spot prices. The marketclearing price for the energy exchange between the MGs upon convergence of the ADMM is shown in Figure 10C, where the electricity price at the PCC provides the baseline price for the electricity exchange. Since the BESS is utilized to behave as a price arbitrager, the exchange price is flattened. It can also be observed that consensus is reached for the energy import/export on the coupled bus. Due to the existence of BESS as a price arbitrageur in balancing the energy supply and demand, the energy exchange price between multiple MGs are flattened compared to the spot price at the PCC. The validation is further extended under binding thermal constraints for line flows, where additional constraints are added to the line current magnitude in MG 3. Consequently, MG 3 is fully supplied by its local DG generation with a higher marginal cost. We plot the electricity exchange schedule in Figure 11B, where zero electricity import for MG 3 can be observed. Compared to the non-congested scenario, the marketclearing price for energy exchange is increased as shown in Figure 11A. This is because the MGs have to pay the congestion part in the exchange price in addition to the baseline price as discussed in eq. (32).

CONCLUSION AND OUTLOOK
Moving toward the grid decentralization, a market framework for clearing the energy trading in the networked MG system is proposed in this work. The proposed ADMM-structured robust optimization method allows the robust scheduling for energy dispatch for DERs and MGs considering uncertain electricity prices, load, and DG generation variations while each local MG operation preserves its autonomy in the energy trading to achieve the optimality of the overall system. The numerical analysis provides insights for price settling when applying such a framework to different scenarios (cases). Note that from the market design perspective, we assume that each MGO is a price-taker without modeling their strategic behavior. In cases with a limited number of market participants (e.g., only two geographically connected microgrids), market power can be easily exercised, whereby an ADMM-based market mechanism cannot guarantee social-welfare optimality. Therefore, gametheoretical approaches are more appropriate for solving the energy exchange problem and providing the market design for these cases. Another interesting direction is to investigate the co-optimization of reserve provision of the networked MG system, which is considered as one of the future extensions for this proposal.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.