Stochastic Predictive Energy Management of Multi-Microgrid Systems

: Next-generation power systems will require innovative control strategies to exploit existing and potential capabilities of developing renewable-based microgrids. Cooperation of interconnected microgrids has been introduced recently as a promising solution to improve the operational and economic performance of distribution networks. In this paper, a hierarchical control structure is proposed for the integrated operation management of a multi-microgrid system. A central energy management entity at the highest control level is responsible for designing a reference trajectory for exchanging power between the multi-microgrid system and the main grid. At the second level, the local energy management system of individual microgrids adopts a two-stage stochastic model predictive control strategy to manage the local operation by following the scheduled power trajectories. An optimal solution strategy is then applied to the local controllers as operating set-points to be implemented in the system. To distribute the penalty costs resulted from any real-time power deviation systematically and fairly, a novel methodology based on the line ﬂow sensitivity factors is proposed. Simulation and experimental analyses are carried out to evaluate the effectiveness of the proposed approach. According to the simulation results, by adopting the proposed operation management strategy, a reduction of about 47% in the average unplanned daily power exchange of the multi-microgrid system with the main grid can be achieved.


Introduction
In recent years, the multi-microgrid (MMG) concept has emerged to improve the operational and economic performance of distribution networks [1,2]. Coordinating the operation of interconnected microgrids (IMGs) will result in more efficient energy usage and less frequent power exchange with the main grid (upstream network). Moreover, distribution system reliability will be enhanced by microgrids (MGs) involvement in supporting activities including sharing surplus power with the adjacent MGs. However, to exploit these potential benefits, developing an efficient energy management system (EMS) for the entire MMG network is required [3]. The most important issues in establishing appropriate models and control strategies are related to the complexity that is resulted from interaction of multiple MGs as well as the uncertainties that are caused by the intermittent nature of the renewable energy (RE)-based distributed energy resources (DERs) and variability of loads. An overview of the integration challenges of REs as well as smart grid control and reliability issues can be found in [4,5].
The development of efficient EMSs for IMGs has drawn the attention of many researchers. Especially, centralized energy management approaches have been widely studied during the last years. In [6], a two-stage scenario-based programming approach is developed in the centralized form to manage the operation of IMGs in a cooperative manner. In [7], model predictive control (MPC) is adopted to coordinate the operation of a renewable-based MMG system. However, no uncertainty management strategy is used. Accommodating the intermittent nature of renewable energy sources (RESs) production and load fluctuations robust techniques are used in [8,9]. In all of the reviewed papers, a central controller is responsible for devising the operation strategy of the whole system considering the internal dynamics of individual MGs and interacting parts. Although adopting this optimization approach could result in an optimal solution, the problem may not be always computationally tractable especially when the number of system players (e.g., MGs) increases.
Accordingly, distributed and hierarchical approaches have been developed in which decisions are made at different control levels. Although the final solution of distributed approaches will be sub-optimal compared to the centralized strategies, a substantial reduction in computational time will be acquired. Multi-agent systems [10,11], game theory [12,13], decentralized optimization techniques based on Alternating Direction Method of Multipliers (ADMM) [14], and distributed MPC [15][16][17] are among the most common methodologies that have been used to decompose the energy management problem of a network of IMGs into several smaller optimization problems. In [18], power transfer between multiple AC microgrids is scheduled by using distributed cooperative control. Two inter-cluster and intra-cluster schemes are proposed to manage power sharing and regulate DER's voltage and frequency, respectively. In [19], smart MGs are modeled as a team of cooperative agents. Decentralized MPC is used to manage the energy level of storage devices and power exchange of MGs around pre-defined reference levels. In [20], a two-level hierarchical EMS is proposed for an MMG system connected to the main grid. Karush-Kuhn-Tucker optimality conditions are used to decompose the problem into smaller subproblems. In [2], a cooperative distributed MPC is proposed to achieve an economic performance of MMG systems. Supply and demand uncertainties are handled through chance constraints. However, MGs can not exchange power directly and the distribution system operator (DSO) plays an intermediary role. In [21], line agents are responsible for determining optimal amount of power to be transferred through connecting lines of cooperating MGs. The distributed robust technique is utilized to manage the uncertainty of generation and consumption. In [22], a hierarchical MPC methodology is proposed for energy management of MG clusters while taking into account different sources of uncertainties. In [23], an EMS based on chance-constrained MPC is proposed in which the cooperation of different MGs is modeled using a joint probabilistic constraint. In [24], the energy exchange of MMG systems with the main grid is studied from the perspective of the DSO. Reducing the demand-side peak-to-average ratio (PAR) and maximizing the profit from selling energy to the MMGs is the main goal of the DSO. Authors in [25] present a decentralized-distributed adaptive robust optimization method to address the distributed scheduling of AC/DC hybrid MMG systems considering accidental communication failures. An approximate reinforcement learning (RL) methodology is used in [26] for bi-level power management of networked MGs considering the limited knowledge of MGs behavior. The RL-based learning method is used to discover the relationship between the retail price and the exchanged power of MGs. Recent review studies on energy management and control of MMG networks can be found in [10,27,28].
According to the current literature, the majority of the studies in this area concern development of the high-level controllers. The common assumption is the existence of efficient local controllers that are capable of appropriately following the scheduled set points. In other words, no analysis is conducted to indicate the applicability of the developed controllers to be incorporated in the full control hierarchy, which can guarantee power system technical requirements. Moreover, there are only a few experimental results in the context of multiple IMGs [29,30]. Besides, the potential of the two-stage stochastic MPC in operation management of MMG systems has not been considered. Finally, considering the significant role of the MMG system in the next-generation power systems, designing systematic and fair cost allocation methodologies to motivate MGs participation in the cooperative energy management strategies is of vital importance.
In this regard, a hierarchical EMS is developed to tackle complexities in the integrated operation management of interconnected RE-based MGs in this paper. MGs in a neighborhood area desire to keep the power exchange level of the MMG system with the main grid around a reference trajectory. Considering the hierarchical structure of MGs control [31], the proposed methodology includes three control levels. The central energy management system (CEMS) at the system-level, is responsible for coordinating the operation of MGs and determining day-ahead optimal trajectory for power transactions with the main grid considering system technical requirements. During the operation day, in a 60-min scheduling cycle, local EMSs (LEMSs) at the MG-level are responsible for managing the charging/discharging process of storage devices based on the power references received from CEMS. Two-stage stochastic programming with recourse variables is adopted to mitigate the uncertainty of non-controllable generators production including wind turbine (WT) and PVs as well as the variability of loads. Finally, local controllers (LCs) are designed to follow the provided set points through LEMS of each MG. To distribute any penalty cost resulted from real-time power deviation systematically and fairly, a new method based on the line flow sensitivity factors is proposed in which penalty cost is allocated by following the MGs contribution in the power deviation. The effectiveness of the proposed approach is verified by hardware-in-the-loop (HIL) simulation study.
The main contributions of this paper can be summarized as follows: 1. Modeling integrated energy management of an MMG system in a hierarchical framework; 2. Proposing an efficient control strategy using MPC and two-stage stochastic programming with recourse; 3. Developing a novel methodology for cost allocation among MGs using the line flow sensitivity factors.
The rest of this paper is organized as follows. Problem formulation is presented in Section 2. Two-stage stochastic programming with recourse and stochastic MPC are represented in Sections 3 and 4, respectively. The proposed cost allocation procedure is introduced in Section 5. Sections 6 and 7 are related to simulation and experimental results, respectively. Finally, Section 8 concludes the paper.

Problem Formulation
Adopting the dynamical system approach, each subsystem i can be represented as a stack variable with associated inflows and outflows. The inflow is related to the net MG production including the aggregated power generation of WTs (P i WT,r (t)) and PV units (P i PV,r (t)) as well as the power consumption (L i (t)) as represented below.
In (1), P i WT,r (t) and P i PV,r (t) represent the output power of the rth renewable unit in the ith MG at hour t. Besides, N i WT and N i PV are related to the number of WTs and PVs in the ith MG, respectively. On the other hand, according to the following equation, the outflow of the ith MG is the total power that is transferred from it to the neighboring MGs and the main grid where N i is the set of the ith MG neighbors including the main grid.
Assume that IF i (t) represents the net power balance of the ith subsystem which could be forecasted based on the historical data during the scheduling horizon T. MGs provide the CEMS with the day-ahead predicted values of IF i (t) with 1-hour time resolutions. Utilizing this abstracted form of information, each MG is modeled through its net power balance from the CEMS point of view, disregarding further details and complexities. Afterward, system-wide power flow analyses are conducted by CEMS to specify the amount of power to be transferred through power lines considering voltage and frequency constraints as well as system losses. Thus, the system is expected to be in balance based on the day-ahead predicted values, which means Here, M is the total number of subsystems including the main grid with the assumption that the highest index is assigned to the upstream network.
After deciding day-ahead scheduling, the entire MMG system will be represented to the main grid as an aggregated load or generation unit through its net energy content (see (3)).
Considering day-ahead power scheduling, IMGs are required to exchange power with the main grid according to the pre-specified trajectory of I M (t). However, the intrinsic uncertainty of RESs production and the variability of loads as well as the forecasting errors will result in uncertain behavior of MGs' real-time power balance quantities. It is assumed that any real-time power deviation will be technically compensated by the main grid to guarantee the stability of the system. However, a large penalty cost will be imposed on the MGs. Consequently, MGs are required to devise operation strategies to compensate for the uncertainty locally at the MGs level. In this paper, to reduce the computational complexity, it is proposed that a local variable named power mismatch can be defined for the individual MGs as follows, for which the desired value is zero.
In other words, it is desired that the realized outflow of each MG (ÔF i (t)) follows the estimated trajectory for the outflow of the MG (OF i (t)). However, accounting for the inherent uncertainty of the MGs, this constraint cannot be fully guaranteed before the realization of the uncertain parameters. Thus, it is assumed that each MG is equipped with an energy storage system (ESS), which can be used as an appropriate tool to mitigate the uncertainty. The cost function of each MG is to minimize the operation cost of ESS unit as represented in the following.
In (5), c i ( $ kW ) represents operational cost of the ESS in the ith MG and P b,i (t) is the amount of charging/discharging power, which is assumed to be positive/negative during charging/discharging intervals.
The LEMS of each MG is responsible to efficiently devise charging/discharging strategies of ESSs in the associated subsystem. Considering the uncertain nature of generation and demand parameters, LEMS tries to keep the possible power mismatch between realized and predicted trajectories of the MG outflow around zero while minimizing the ESS operational cost. Consequently, the control problem of each MG can be considered as a decentralized optimization problem with the following constraints in which (10) is a random constraint.
In (6), SOC i represents state of charge of the ESS in the ith MG and E b,i is its nominal capacity. Besides, ∆ refers to the time length between two consecutive steps, which is assumed to be equal to 1h in this study. In (9),P i WT,r (t),P i PV,r (t), andL i (t) refer to the realized power values of the rth WT and PV units and the total power consumption in the ith MG at time slot t, respectively.
Taking into account the dynamical equation of the storage devices, successive time steps will be related to each other. Accounting for the system's dynamics and accommodating the stochastic nature of the tracking constraint, a two-stage stochastic MPC approach is adopted to achieve the desired solution strategy. The optimal operating set points designed at this level are followed by the local controllers of each MG. Figure 1 illustrates the proposed hierarchical control structure along with the required input data and associated outputs of each control layer.

Two-Stage Stochastic Programming with Recourse
In the two-stage stochastic programming with recourse variables, decision variables include the first and the second-stage variables. The first-stage decision variables are determined before the realization of uncertainty while the second-stage variables are set after uncertain parameters are realized and considered as recourse variables (actions). Taking into account the second-stage variables for compensating any real-time infeasibility, the first-stage decision variables are determined with the allowance of constraint violation as a result of changing random parameters [32,33]. In other words, in this model, the constraints which include random parameters are modeled as soft constraints. To have a more quantitative point of view, consider the optimization problem illustrated in the following.
In (11)- (14), x ∈ R n 1 represents the first-stage decision variables and µ ∈ R r contains the uncertain parameters. In (12), A and b belong to the R m 1 ×n 1 and R m 1 , respectively. While the constraint represented in (12) is a deterministic constraint, the random constraint of (13) includes both deterministic and random parameters where T ∈ R m 2 ×n 1 and H ∈ R m 2 ×r are coefficient matrices.
Since the only information about the uncertain parameters is their past behavior which could be incorporated to extract some statistical characteristics, these kinds of constraints could not be guaranteed to hold before the realization of random variables. In the two-stage stochastic programming, some slack variables γ + and γ − are introduced into the model as the second-stage decision variables as represented below.
In (15) and (16), γ + and γ − are slack variables, which represent corrective actions for compensating real-time infeasibilities. To deal with the violated random constraints, a second-stage or penalty function is included in the objective function. In general, the penalty function is considered in the form of q + γ + + q − γ − in which q + and q − are penalty coefficients.
However, since the penalty function could only be evaluated after the realization of uncertain parameters, the second-stage cost could be replaced with an appropriate measure of risk. Thus, the first-stage decision variables should be selected in a process that encompasses possible realizations of µ represented by a finite set of scenarios. Considering the expected value for evaluating the probable penalty cost of different realizations of the uncertain variables and assuming finite discrete probability distribution for µ, the two-stage stochastic optimization problem can be stated as follows.
In (17)- (21), S represents the number of scenarios and π j shows related discrete probability of the jth scenario. Thus, in a two-stage programming problem, the objective is to decide the first-stage decision variables such that the total cost including the operation cost and expected penalty cost resulted from real-time constraint violation is minimized. Although increasing the number of scenarios will result in improving quality of the solution, a long computational time might be required. Consequently, the solution quality should be traded against tractability of the problem.

Two-Stage Stochastic Model Predictive Control Formulation
In this section, the operation management problem of IMGs is formulated in the framework of two-stage stochastic MPC with recourse. In MPC, based on a dynamical model of the system under consideration, an optimal control sequence is derived over the desired control horizon. Intrinsic robust characteristics of MPC due to the rolling horizon strategy as well as its capability to account for operational and time-coupling constraints, make it a consistent solution strategy for power system applications.
Based on the two-stage stochastic approach that has been outlined in the previous section, the operation management problem of MGs can be represented as follows: min P b,i (t,...,t+H u −1), In (22)-(28), H P and H u refer to the prediction and control horizons, respectively. In addition, π i,j represents the occurrence probability of the jth scenario in the ith MG and γ − i,j (t + k) and γ + i,j (t + k) are the recourse actions to take if the jth scenario occurs in the ith MG at time instant t + k while staying in t. Furthermore, δ mis,i,j (t + k) is the possible scenario of the uncertain parameter δ mis at the time-step t + k in the jth scenario and the ith MG.
At 60-min time intervals, the LEMS of each MG measures the level of stored energy in the battery and solves the optimization problem (22)- (28) to obtain the optimal control sequence P b,i (t, ..., t + H u − 1) and recourse variables γ i,j (t, ..., t + H p − 1). The first sample of the optimal control sequence is applied to the local controllers and the control and prediction horizons are shifted. Consequently, at each time step there are H u + 2 × S × H p decision variables and 2 × (H u + (2S + 1)H p ) constraints for each MG.

Efficient Cost Allocation Strategy
In this paper, it is assumed that MGs in a neighborhood area desire to represent the whole MMG system as a predictable entity to the main grid. In other words, the final goal is to minimize the net value of unexpected power transactions with the main grid. According to (3), the net exchanged energy of the MMG network with the main grid is equal to the summation of the flow of lines that connect MGs to the main grid. It is assumed that the network of MGs will be penalized according to the real-time power deviation from the scheduled plan as represented in (29) where PCost stands for the penalty cost and λ represents the cost factor ( $ kWh ) imposed to the MMG system for one kWh energy deviation. PCost = (Î M (t) − I M (t))λ. (29) In this paper, it is proposed that the final penalty cost could be distributed among MGs (λ = ∑ M−1 i=1 λ i ) by following their participation level. In other words, one should calculate the effect of power deviation in different MGs in the change of I M (t) from the reference value. Considering an exemplary system diagram for M = 4 as represented in Figure 2, changes in I M (t) is equal to the summation of the line flow changes in different lines that connect the MGs to the main grid as shown below.

Power Grid Microgrid 3
Microgrid 2 Microgrid 1 In (30), ∆ f n−m represents the change in power flow of the line n − m, which connects the nth bus to the mth bus. According to [34], it is assumed that any power changes ∆P i at the ith bus of the system will be compensated through an exact opposite change in the reference bus, which is assumed to be the main grid in this study while the power generation in other buses remains fixed. Since power changes in the generation of different MGs do not have the same effect on the power flow of the intended line, the line flow sensitivity factors are used to measure different MGs effect and to share the penalty cost fairly.
Difference between the expected line power flow ( f n−m (t)) and the realized line power flow (f n−m (t)) can be represented as follows: Considering the effects of power changes in different MGs on the line n − m and defining α (n−m),i as the sensitivity of line n − m power flow to a change of ∆P i in the ith MG, ∆ f n−m (t) could be written as below.
Utilizing (30), the change in power flow of the reference line concerning different MGs can be calculated as follows: According to (35), the change in the flow of a line can be written using the phase angle change in the two buses n and m where the line is connected [34]. In (35), x n−m represents the line reactance.
Moreover, using the linear power flow equations, the phase angles and power change of buses are related according to (36) where [X] is derived from the admittance matrix [34].
Hence, the line flow sensitivity factors are derived as follows: It should be noted that if the lth bus is the reference bus, then X l−i = 0, X i−l = 0. Utilizing (34) and (37), the penalty cost coefficient of each MG can be considered as the normalized value of the line flow sensitivity factors.

Simulation Results
In this section, to evaluate the effectiveness of the proposed approach, simulation analyses are carried out in a test system. The model is established in MATLAB/Simulink. The nominal voltage and frequency of the system are set to 230 √ 2V and 50Hz, respectively. Matlab SimPower Systems is used to model the electrical system. The line impedance information is represented in Table 1. Table 1. Electrical system parameters.

Parameter
Symbol Value Frequency f 50 Hz Nominal Voltage V nom 230 √ 2 V Line Impedance 1,2 Z 1,2 R 1,2 = 0.4Ω, L 1,2 = 1mH Line Impedance 1,3 Z 1,3 R 1,3 = 0.4Ω, L 1,3 = 1mH Line Impedance 1,4 Z 1,4 R 1,4 = 0.4Ω, L 1,4 = 1mH Line Impedance 2,3 Z 2,3 R 2,3 = 1Ω, L 2,3 = 1.7mH Line Impedance 2,4 Z 2,4 R 2,4 = 1Ω, L 2,4 = 1.7mH Line Impedance 3,4 Z 3,4 R 3,4 = 0.7Ω, L 3,4 = 1.5mH The simulated system consists of three RE-based MGs that are equipped with battery-based ESSs in size of 350, 300, and 400 kWh, respectively. The initial level of SOC is set to 50% of the nominal capacity, while the min and max SOC values are set to 20% and 80% of the battery nominal capacity, respectively. It is assumed that all MGs and the main grid are electrically connected and information can be exchanged among them to support bidirectional updates. Table 2 represents the predicted aggregated load profile of each MG during the optimization horizon and the forecast data for aggregated RESs production. To generate discrete random scenarios representing the uncertain nature of RESs production and consumer demand, Monte Carlo algorithm is used. Random scenarios are generated from normal distributions where the mean values are set to predicted data and standard deviations is assumed to 5% of the mean values for all MGs. Besides, H P = H u = 3 h and T = 24 h. The cost function is evaluated for a representative case where c i = 0.2 $ kW and λ = 2 $ kWh . Power flow analysis is conducted by CEMS using MATLAB SimPower System to specify the amount of power to be transferred through power lines. All electrical system requirements including voltage and frequency constraints, line capacity, and system losses are considered for deriving optimal power flows. A local PI controller as represented in Figure 3 is designed that controls each inverter output in current control mode. The current proportional and integral terms are set to K pI = 20 and K iI = 50[s −1 ], respectively. For simulation purposes, the main grid is modeled using an inverter, which is controlled in voltage control mode as illustrated in Figure 4. For the experimental tests this part will be replaced by grid simulator. The current and voltage proportional and integral terms are set to K pI = 15, K pV = 20, K iI = 50[s −1 ], and K iV = 50[s −1 ], respectively. Using the proposed strategy for penalty cost allocation, the penalty coefficients have been calculated: s 1 = 0.3184, s 2 = 0.3518, s 3 = 0.3299.  Figure 5 shows the daily ESSs scheduling for the MGs including the normalized SOC and charging/discharging power of the batteries, respectively. Scheduled and realized exchanged power between the MMG system and the utility is also plotted in Figure 6. As it can be noted, the proposed EMS has a good tracking performance in following the power reference trajectories.   To evaluate the effects of the proposed approach on the unscheduled power exchange level of the MGs with the main grid, a single operation mode is assumed. In this mode, each MG is only allowed to exchange power with the main grid. As can be seen in Figure 7, the sum of power deviations of individual MGs is considerably higher than one related to the MMG system and a reduction of about 47% in the average unplanned daily power exchange of the IMGs with the main grid can be achieved. This can be considered as an advantage of the proposed strategy in reducing power system unscheduled power transactions. The penalty cost of MGs in both operating modes is tabulated in Table 3. According to the table, through coordinating the operation of MGs the penalty cost of each MG and the whole system will be reduced. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Time

Experimental Results
The test system has been implemented at the AC Microgrid Research Laboratory at Aalborg University as represented in Figure 8 [35]. A real-time platform (dSPACE 1006) is used for real-time simulation. As it can be noted, the setup includes three inverters (2.2 kW, FC302) fed through a stiff bidirectional DC source (0-800 VDC, 0-50 A, 64 kW), a regenerative grid simulator (45 kVA|330 V L-N) as well as line models, and load boxes. Due to the technical limitations, each MG is modeled by a single inverter. However, as the paper aims to show the capability of the MGs to follow reference trajectories, it will not raise any issue. It should be noted that generation/consumption profiles and the time slot have been scaled down such that 1000 W:1 W and 3600 s:20 s. Experimental tests are exemplary carried out in the time interval 18:00-20:00. Figure 9 represents the tracking performance of the local controllers. As can be seen, a good tracking performance has been achieved during the test interval. Moreover, taking into account the time scaling factor, the transient time of MGs is about 6s which approves the efficiency of the proposed approach in calculating and implementing optimal operation strategy and makes it applicable for real EMS. System frequency and voltage profile at PCC are also denoted in Figure 9. As can be seen in this figure, the system frequency has been kept within an acceptable range and a three-phase balanced output voltage regulation has been achieved.

Conclusions
In this paper, a novel control structure for integrated operation management of IMGs was proposed. The main goal was to represent the MMG system as a predictable entity to the main grid. To that end, a decentralized tracking control problem was designed for each MG and a two-stage stochastic MPC approach was adopted to derive the solution strategy. In the proposed scheme, any power deviation of the MMG system from the scheduled trajectory will be penalized. Considering the significant role of the MMG system in the next-generation power systems, designing systematic and fair cost allocation methodologies to motivate MGs participation in the cooperative energy management strategies is of vital importance. In this sense, a novel methodology for allocating the penalty cost among MGs using the line flow sensitivity factors was proposed. Analyzing other systematic cost allocation strategies is under study by the authors. Experimental results demonstrated the capabilities of the proposed approach to be incorporated as an efficient high-level controller in practical applications. According to the obtained results, adopting the proposed operation management strategy, a reduction of about 47% in the average unplanned daily power exchange of the IMGs with the main grid can be achieved. The reduction of unplanned power exchange with the main grid will result in a more predictable power exchange of the MMG system with the main grid as well as less penalty cost. From the viewpoint of the main grid, since the unscheduled power exchanges can be seen as disturbances, reducing it will result in lowering the complexity of power management of the main grid and more efficient reserve management strategies. Future studies of the authors are related to integrating learning-based methodologies in the EMS of MMG systems. Learning the behavior of MGs under different operating conditions will help develop more efficient EMSs in a less uncertain environment. Besides, recent advances in real-time monitoring and situational awareness of MGs will be deployed for early detection of system behavior changes from the predefined trajectories to activate the required responsive actions.