Cooperative Optimization of Networked Microgrids for Supporting Grid Flexibility Services Using Model Predictive Control

The transition towards fully renewable energy-based power systems will require to increase the number of reserves at the System Operators’ (SOs) disposal to provide flexibility on the energy management process. The microgrid’s ability of integrating distributed energy resources, loads and energy storage systems (ESS) appears as a powerful flexibility tool. Nevertheless, the associated control problem of microgrids increases with the number of connected devices. A structuration of the distribution grids in networks of microgrids is proposed, focusing on their ability to provide flexibility services. The complexity of the associated optimization algorithm is faced using Distributed Model Predictive Control (MPC). The algorithm is divided in two steps. The first one is applied to the cooperative participation of microgrids in the day-ahead market. The second step covers the interaction with the SO offering flexibility services in exchange for a financial benefit. The financial benefit is optimally shared between the networked microgrids to satisfy the power profile requested by the SO at the lowest cost. As the proposed control algorithm presents both continuous and binary variables, its associated optimization problem is formulated using the Mixed Logic Dynamic (MLD) framework, which results in a Mixed Integer Quadratic Programming problem (MIQP).


LOH
Level of Hydrogen (Nm 3   T HE EXPECTED penetration of distributed energy resources (DER) in the following years will mainly affect two actors: Market and System Operators. Market Operators (MOs) are entities responsible of managing day-ahead and intraday markets [1], receiving buying and selling bids from the participants and planning the long-term operation of the electric system for the next day. While the complexity of these long-term electrical markets is currently being increased with the presence of DER, the wide schedule horizon available for the optimization of these markets, which is usually done the day before they are executed, makes it not so difficult to deal with the expected renewable energy penetration. SOs are the entities responsible for guarantying the real-time operation of electric power systems. They are in charge of maintaining a continuous balance between electricity supply and demand, ensuring the provision of reserves. This reserve provision is carried out through the short-term and real-time electrical markets management. For the current fossil-fuelsbased power system, a reduced number of flexibility reserves based on hydropower and gas turbines have been enough to balance the power system. Nevertheless, there will be a huge problem concerning flexibility of power systems in the envisaged paradigm of full penetration of renewable energy systems. Microgrids appear as a key solution to increase the flexibility for a higher penetration of DER. The transformation of the smart grids towards more structured systems based on microgrids networks opens the path to increase the number of flexibility reserves. Cooperative control of microgrids will also be an important issue to consider in order to base the energy consumption of large facilities or installations such as university campuses or military bases on renewable energy generation.

A. Literature Review and Paper Contributions
The control of microgrids is a complex issue where several aspects have to be addressed. The control problem starts with their economic optimization through their participation in electric markets (which concerns to their long-term schedule for periods longer than one day) having to solve aspects related to power quality (which require of sampling periods of millisecond or microseconds). In order to relax the computational requirements of the control platform utilized by the microgrid EMS, a three-level hierarchical control structure is usually implemented. Aspects related with the long-term optimization of microgrids, such as their participation in electricity market, are executed in the tertiary control level which sends the power-reference set points to the secondary and primary control level. Both the control horizon, as well as the sample period are greater in the tertiary control than in the rest of control levels [1], [2].
The cooperative control between networked microgrids is introduced by Lasseter in [3]. A review of different flexibility products, as well as their importance to manage the growing presence of renewable generation is detailed in [4]. Microgrids as flexibility resource has been object of study in several papers such as [5], [6]. A recent study concerning the different Energy Management System (EMS) applied to interconnected multi-microgrids is found in [7]. The authors of [8] focus their research on the energy management problem for microgrids, considering that the internal demand may exceed internal power supply and energy stored. A three level hierarchical coordination strategy is proposed. The top level is responsible for energy coordination between the SO and the microgrid EMS. The SO purchases/sells energy from/to a microgrid that has surplus/deficit energy at a slow sampling period. In [9], a multi-agent-based rolling optimization is developed for restoration scheduling of distribution grids with distributed generation. A cooperative market mechanism is proposed in [10] to define the energy transactions and market price in multi-microgrids. In [11], a cooperative multi-objective optimization for the networked microgrids is carried out. A multi-microgrid operator is introduced in order to reduce the total daily costs. Microgrids can participate in demand response programs through the use of flexible loads which can response to price signals. Studies related with the benefits of a network configuration for microgrids can be found in the review carried out in [12]. The associated control problem of microgrids networks as system of subsystems presents even higher complexity than EMS applied to microgrids which act as single systems. While their configuration as networked system can bring several benefits to integrate flexibility in the power system, as well as to reduce the final energy cost, the development of complex EMS for network of microgrid requires an advanced control methodology. In the last years, the MPC framework, as a widely control structure to solve complex problems in industry, has been widely applied to microgrids control problems [1]. The methodology of this family of control methods allows to optimize a multi-objective cost function where aspects as the future behavior of the microgrid devices, the energy forecast and the price prediction can be easily integrated. The most direct system decomposition in a network of microgrids consists of considering each one of the microgrids as a subsystem, with the same state variables, inputs, and disturbances as considered when modeling single microgrids, but adding new input or disturbance variables: the energy exchange between neighbor subsystems [1]. In distributed MPC, each subsystem is controlled by a local MPC controller, but the key point is the communication and what information is exchanged among the different controllers, and also how and when the transmission is performed. This formulation allows the interaction between different agents searching a common objective [13]. DMPC controllers has been recently applied to optimize networked microgrids as can be seen in [14], [15]. An overview of DMPC techniques applied to microgrids can be found in [1]. In the aforementioned references, there is a lack of studies regarding the use of networked microgrids as a tool available to the SO to increase or decrease the energy exchange with the main power grid for real-time system balancing. Sharing the storage resources of the networked microgrids can be an effective way to reduce the associated costs of providing flexibility services to SOs in a more effective way than acting as single microgrids.

B. Main Contributions
The main aim of this article is to develop an optimization algorithm to be applied to the tertiary control level of networked microgrids with the functionality of changing the schedule carried out in the day-ahead market as response to a certain power profile variation request from the SO. This variation of power profile will be linked to an economic benefit offered by the SO. In order to optimize the associated cost, the economic benefit is shared between the different microgrid agents depending on their contribution to the supply of the requested power profile. The best relation between power variation and associated cost is searched by the control algorithm. The following are some of the features of the proposed algorithm for the cooperative optimization of networked microgrids considered as main highlights of the present work: • In order to increase the offered flexibility by the microgrids, the problem is formulated for microgrids with hybrid ESS composed by batteries and hydrogen. • The EMS considers the different degradation mechanisms of the hybrid ESS components, as well as their operation costs. It is formulated as a distributed MPC-controller which includes the energy exchange between microgrids as a flexibility resource. • The EMS interacts with the SO for flexibility services as response to an offered benefit and flexibility power profile requirement, it distributes the power variation associated to the service among the cluster of microgrids in order to find the optimal cost. In order to provide a complete framework which can be applied to microgrids whose optimization problem contains a cost function formulated as mixed integer quadratic or nonlinear problems, the developments given in [16] are used for the scheduling of each single microgrid. The contribution of the present work regards the expansion of the methodology to carry out the scheduling of a network of microgrids including the provision of flexibility services. Once the schedule is carried out, the procedure given in [17] for real-time control of the microgrids can be applied.

C. Structure and Organization of the Paper
The remain of this article is organized as follows: Section II develops the algorithm to optimize the network of microgrids under a flexibility service request from the SO. In order to introduce the interaction between the different microgrid agents, it is formulated as a Distributed MPC. In Section III, the results of the proposed controller are exposed and discussed based on different case studies developed. Finally, the main conclusions are outlined in Section IV.
An overview of a smart grid with centralized reserves at disposal of the SO, as well as distributed reserves through the structuration of the distribution grid in network of microgrids with hybrid ESS which interact with the SOs, as proposed in the present work is shown in Fig. 1. The mentioned figure also depicts the communication among the different electrical market agents. As can be seen, the microgrids integrate two kinds

II. DMPC CONTROLLER DESIGN
The optimization algorithm is designed to incorporate the following features: 1) Operational Cost of the Microgrid Components: The algorithm integrates the different economic cost costs related to lifetime, operation and degradation issues. 2) Cooperative Day Ahead Operation: The algorithm computes the best operation of the whole network of microgrids in the day-ahead market participation through cooperative control actions between the different microgrid agents.

3) SO's Signal Response:
The algorithm has the capacity to quantify the best power profile shared by each microgrid in order to guarantee the minimum associated cost to the flexibility services providing. The optimization algorithm is structured in two steps. The first one computes the energy exchange of the whole network with the main power grid in the day-ahead market. The second step searches the optimal power shared by each microgrid as response to a common request of flexibility services from the SO. The controllers are developed using the block diagram shown in Fig. 2.

A. Distributed, Hybrid and Agent-Based Negotiation Approaches Using MPC
MPC main idea is to minimize a cost function J expressed as a function of the control variables which are subject to a set of constraints concerning their physical limits or the relation between them. This relationship is usually modeled with a state space representation. The optimization problem is usually expressed as a quadratic cost function. The aim of the optimization is to minimize the difference of the system variables values respect a desired reference. This minimization of the cost function is calculated over a control horizon of N c predicted sample instants in order to decide the optimal values for the control variables considering a prediction horizon N p for the predicted outputs based on the state representation model of the system [1].
Complex control systems often integrate both bynary and continuous variables. In these cases, the so-called hybrid formulation of MPC controllers is utilized. This methodology is based on the application of the foundations of the MLD framework introduced in [18] by Bemporad and Morari. The main idea is to use the propositions to transform relations between logic and continuous variables in linear constraints which are introduced in the controller. The control variables are separated in three main groups: u as continuous control variables, δ the logic variables and finally, z defined as mixed products between logic and dynamic variables. The hybrid formulation of MPC controllers also allows to replace the non-linear formulation of the cost function by a linear one, replacing the non-linear terms by piecewise expressions [19].
The main concept of DMPC-based controllers [13] is to determine how different local controllers can interact in order to find a global optimum solution which gets a lower global cost than acting independently. Cooperative DMPC controllers search a common objective through the interconnection variables while in the case of non cooperative ones, each controller will find a better optimization point of its local control function through the interaction between the different agents. A methodology for DMPC is the agent-based negotiation which is detailed in [13].
This methodology implies a negotiation phase after which agents take a cooperation or no-cooperation decision. At each sampling instants the different agents make proposals to improve their local cost function, state and model. These proposals are accepted if they improve the global cost corresponding to the current solution.
The global cost function J global is defined in (1). It integrates the sum of the local cost functions J local of each controller (2) adding a term related with the interaction between the different subsystems. The superscript (i) refers to the number used as label for each subsystem. The interaction between the different subsystems through the interconnection variables v (i,j) allows to find a better operation point working as global system than acting as independent systems. The nomenclature v (i,j) is defined as an output of the subsystem (i) which acts as input or action in the subsystem (j). The control problem is completed introducing the constraints given by the state space model of each subsystem. The theoretical formulation of a hybrid and distributed representation of the state space model which is defined by means of the expressions (3)-(5), where x (i) stands for the state variables, which are defined as those ones whose values depend on their past values, and d (i) represent the non-manipulable inputs or disturbances. Finally, y (i) are the outputs of the systems. The subindex ref are the reference values for each variable of the subsystems. The system constraints is given by expressions (6)-(10). The controller finds an optimal operation point for the cooperative DMPC problem defined by a set of optimal control vari- subject to: After having found the optimal operation point as cooperative system, an external agent (the SO in our case study) can request a modification of the operation point of the system of subsystems offering a benefit B ext in exchange. The control problem in this case, considering it has similar constraints as defined in expressions (3)- (10), can be expressed as defined in (11). As it occurs in the previous case, the controller will find the optimal point closest to the request sent by the external . The system of subsystems will agree to this variation in its optimal operation point defined as OPT only if the criterion, expressed as constraint integrated in the controller, is satisfied: where the superscript ext and int are used in order to differentiate between the optimization problem with an external agent, equations (11) and (12), and the optimization problem where no interaction of the microgrids network with an external agent is considered, equations (1) and (2). Notice that the interaction with the external agent will suppose a variation in the operation point which will result as a variation of the internal cost of the microgrids network optimization. This cost variation has to be constrained with the external benefit as indicated in (13). Notice that the current state of software destined to solve optimization problems require a high computation cost for solving MLD problems with quadratic or non-linear constraints using normal computers. For this reason, the proposed methodology requires to linearize quadratic or non-linear terms in the cost function using piecewise approaches based on mixed programming to integrate the expression (13) as linear constraint in the controller.

B. Day-Ahead Cooperative MPC Controller for a Network of Microgrids
The problem formulated in previous section can be particularized for the day-ahead operation of a network of microgrids as shown in Fig. 1 following the block diagram shown in Fig. 2. The different arrays of control variables for each microgrid are defined with expressions (14)-(20).
where P (i) bat , P (i) elz and P (i) fc are the reference power values given by the microgrid controller to the internal controllers of the batteries, electrolyzer and fuel cell.
where P (i)→(j) exch (t) are the power exchange among the different microgrids of the network N .
where δ where z The outputs of the plant are defined with the power exchanged with the main grid P (i) grid .
Finally the disturbance array d is obtained by the prediction carried out in the Forecast Module of the controller using the methodology defined in [16], the next disturbance array is incorporated in the controller: where the forecasted energy generation for photovoltaic and wind turbine generators are defined by (P pv (t),P wt (t)), and on the other hand, the load prediction for each sample instant in the microgridP load (t). This module also predicts the energy price for the processes of purchasing and selling energy with the main grid (ˆ pur (t),ˆ sale (t)).
Hydrogen ESS Degradation The local cost function can be defined with the expression (21) based on the case study presented in [16]. In this reference, the cost function has two quadratic terms. As commented in Section II-A, the proposed methodology requires the use of linear cost functions. For this reason, the quadratic terms of the cost function, as initially formulated in [16], are linearized by piecewise expressions formulated within the MLD framework.
In order to simplify, the nomenclature t k = t + k is going to be used in (21) and the rest of expressions along the paper, being SH 0 the first instant of the schedule horizon, SH the last one and T s the sampling period of the controller.
A ch,n · z ch,n (t) + B ch,n · δ ch,n (t) + A dis,n · z dis,n + B ch,n · δ dis,n As indicated in (21), each term is underbraced indicating its contribution to the cost function. The first term of (21) is related with the associated cost or benefit of selling (z sale ) or purchasing energy (z pur ) to the main power grid in the dayahead market. These terms are related with the total energy exchanged with the main power grid through the expressions given in (22)-(25) which are introduced as linear constraints in the controller. DM pur and DM sale are the corresponding energy prices in the day-ahead market. The batteries have two main degradation processes: the cycle aging and the application of high-power values to the charging/discharging processes. The cycling aging cost is introduced in the cost function quantifying the cost of each used cycle, CC is the nomenclature used for the capital cost of each device acquisition. The number of life cycles given by the manufacturer is expressed with (Cycles). In order to disperse the charge and discharge of the batteries as much as possible in time, the current applied/demanded values of charging/discharging power are quadratically penalized in [16]. In order to avoid the quadratic values of P bat , it is linearized with the piecewise approaches given in the expressions (26)-(32). These degradation costs are quantified with Cost degr,α .
In order to avoid anomalous working conditions in electrolyzers and fuel cells, manufacturers recommend maintaining these devices within a power range between 10%-100% of their maximum power [1]. For this reason, electrolyzers and fuel cells setpoints are formulated with logics and continuous variables (δ i , P i ), being z i = P i · δ i . Notice that with this control formulation δ i can be set to '0' providing three working limits for the electrolyzer and the fuel cell [0, P min i , P max i ]. The lifetime of the electrolyzers and the fuel cells have a limited number of working hours which are expressed with Hours. These devices also have hourly operation costs associated to their operation and maintenance Cost o&m,α . Start up cycles incur in degradation aspects to the fuel cells and electrolyzers. These aspects are included in the terms concerning the Hydrogen ESS Degradation. The variable σ stands for the startup cycles. Power fluctuations of the electrolyzer and fuel cell also produce degradation issues in these devices [16]. To minimize the power fluctuations of these devices the auxiliary variables defined in (34) and (38) are integrated in the controller. χ is a logic variable which determines the instants in which the electrolyer and fuel cell are active but are not switching on/off. The mixed variables ϑ + and ϑ − quantify the positive and negative values of the power fluctuations in these instants when χ is active. The last terms in (21) are included in order to maintain the level of stored energy in each ESS close to a given reference value, to face future uncertainties in the next scheduling horizon. Only these values whose difference with the reference is negative are penalized. Further description about each term of (21) can be found in [16].  (39) and (40).
In order to find the output array of each microgrid of the network, composed by the energy exchanged with the main power grid y = [P grid (t)], expressed as a controller constraint in (4), the next energy balance has to be integrated as system constraint in the controller: where the remaining power of each microgrid P (i) rem is defined as follows: Finally, when more than one microgrid are considered in the optimization problem, the next constraints for the relationship among the components of the interconnection variables array v (i,j) (expression (43)) and the global balance in the network (expression (44)), as well as the constraint related with the maximum global power exchange with the main grid of the network (45) have to be also included in the optimization problem: The model constraints shown in (5) are composed by the linear constraints obtained by the transformation of the relationships between logic and continuous variables expressed along the paper. While the optimal way to solve the problem would be formulating the global cost function encompassing all the networked microgrids, using normal computing devices the required optimization time is too high when the problem is applied to more than two microgrids. Hence, the procedure displayed in Algorithm 1 where at each iteration is applied to a couple of microgrids is followed.

C. Cooperative MPC Controller for a Network of Microgrids to Provide Flexibility Services
Once that the schedule is carried out, the external variables of the microgrid P grid (t) are committed. Therefore, they are introduced as constraints in the controller, but the established schedule for the energy stored in the batteries and the hydrogen tank, as well as the exchange power among the microgrids, can be modified according to a new objective. In a certain instant t ext of the day-ahead market execution period, the SO (or an external agent to the network) can request a modification of the power exchange planned with the main power grid, demanding for up/down regulation services P ext req (t) and proposing an associate benefit B ext for this action. One way to minimize the global cost of this action is to find the best way to share this benefit and the up/down regulation power (P (α) up (t), P (α) down ) that each microgrid supplies to the global profile requested by the SO or external agent. To carry out this action, the global benefit is divided in a number of parts m, being B u = B ext /m, the unitary benefit used at each iteration. The global cost function (46) is minimized for each pair of microgrids (i, j) belonging to the network, where w ext is a weighting factor that is adjusted in order to have a magnitude order higher than the rest of terms of the cost function, in order to assign the major importance in the cost function.
The constraint concerning the energy balance in the microgrid (41) in this case has to be replaced considering the up/down regulation power provided by each microgrid and considering the schedule carried out for the energy exchange EXIT 5: end if 6: for i = 1, 2, . . . , N MG do 7: for j = i + 1, i + 2, . . . , N MG do 8: Solve J ext global (i, j) s.t. 9: global )) ≤ R(i) + R(j) + 2B u 10: end for 11: end for 12: Select the couple of microgrids (i * , j * ) which minimizes the value (P flex (t) + α=i,j P (α) 14: Update the values P  with the main power grid in the day-ahead market P sch grid (t).
The procedure detailed in Algorithm 2 is followed, where N MG stands for the number of microgrids in the network. Notice that the flexibility services can be activated at whichever instant of the schedule horizon. In order to compare the variation in the cost function the same number of sample instants and schedule horizon as used in the day-ahead market optimization is utilized. Nevertheless, the values of the optimization variables have to be imposed as constraints in these sample instants which are past to the instant when the SO communicates a flexibility service request to the network of microgrids. A second consideration to be taken into account is the fact that the final exchange with the main grid at this step of the algorithm is defined by the following expression: Both P Finally, the unitary benefit B u which is used at each iteration of Algorithm 2 has to be higher than the cost associated to the variation of the logic variables (σ elz , σ fc , δ elz , δ fc )  within the cost function, in order to guarantee an effect over these variables.

III. RESULTS AND DISCUSSION
In this section, the controller is developed and validated by numerical simulations using TOMLAB as solver tool. The different components of the microgrids can be found in Table I. The values integrated in the controller are included in Table II. The average execution time per iteration is 2 s using the aforementioned computer. In order to simplify, the used sample time for both algorithms is T s = 1 hour and the schedule horizon is 24 hours. Both values correspond to the optimization of a complete day and to the sample time which is usually taken in the day-ahead market. Nevertheless, this sample period could be reduced for the case of flexibility services just considering that the values given in Table II are based on this sample period, therefore their proportional values should be considered for lower sample periods than the selected one. The energy prices predicted using the forecast algorithm for energy selling in the DM are exposed in the graph of Fig. 3, which is based on the Iberian Market Operator history data executed for 05/05/20. It has been considered that DM pur = 2 DM sale .   As commented in Section II-B, the first iteration of the algorithm calculates the participation cost of the different microgrids in the day-ahead market acting as single systems. The associated global cost to the operation of each microgrid, acting as single system in the day-ahead market is displayed in the first row of Table III. The rest of iterations correspond to the couples formation among microgrids. The cost-function decrements obtained by the formation of each microgrid couple in the cost function are shown in Table IV, where the couple of microgrids which leads greatest decrement in the cost function is highlighted with bold text. As can be seen, in the iteration X no decrements in the cost functions are obtained with any of the microgrids couples, so the algorithm for the cooperative operation in the day-ahead market is finished. The final results for the energy exchange among all the microgrids in the day-ahead cooperative participation is displayed in Fig. 4.
In Fig. 5, the response of the whole network is displayed, showing a flexibility requirement and offering an external benefit of B ext = 130 e and B ext = 170 e. The distribution of power profile as well as the re-schedule carried out to supply this service by each microgrid can be observed in Fig. 6 and 7. The changes in the local cost of each microgrid for both cases of exchange of energy with the SO are displayed in Table V. The algorithm 'Day-Ahead Cooperative MPC Controller' compares the increment of the benefit that occurs when establishing the different possible associations of microgrids to   find the optimal set of microgrids' couples at each iteration. The required time to solve the problem for each couple of microgrids was lower than 1.2629 seconds using a PC i7-4510U @ 2.60 GHz. For 4 microgrids the problem is solved within 10 iterations of 6 combinatorial possibilities. The total required time is then 75.774 s. It can be therefore stated that the proposed algorithm is practically acceptable for a dayahead optimization. Concerning the algorithm 'Cooperative MPC Controller for a Network of Microgrids to Provide Flexibility Services' is based on an iterative process that increments an unitary benefit, in the case study of 5 e. The average time to solve the problem for each couple of microgrids is 1.5203 seconds. This unitary benefit can be increased or decreased in order to achieve a shorter computational time or a more accurate result of the economic optimization of the network of microgrids. When the unitary benefit is increased simultaneously for two microgrids at each iteration, the required power profile for flexibility services is obtained with a global benefit of 170 e, what means 17 iterations. As for the day-ahead market operation, each iteration has 6 combinatorial possibilities, thus requiring 155.0706 s. This can be considered as an acceptable time for providing flexibility services.

IV. CONCLUSION
The main aim of the present study has been the implementation and validation under numerical simulations of DMPCbased algorithms to enhance the flexibility and economic competitiveness of microgrids by using their cooperative networked operation and the use of hybrid ESS. Two main case studies are analyzed. The first one improves the participation in the day-ahead market through a cooperative behavior. As can be seen in the exposed simulations, the revenue or cost of the participation is improved with our approach with respect to acting as single microgrids. The second one opens the path to serve as flexibility reserves available for SOs reducing their associated cost since the up/down regulation as flexibility service is optimized finding the best distribution of the requested power profile by the SO among the different microgrids belonging to the network. In 2009, he joined the Centro Nacional del Hidrogeno, Spain, where he is currently responsible for the Microgrid Laboratory, Application Unit. He was working in several research centers or companies, such as the Instituto de Automatica Industrial-Consejo Superior de Investigaciones Cientificas, Spain, the University of Seville spin-off GreenPower Technologies, Spain, and the Université Catholique de l'Ouest, Angers, France. He has coauthored the book Model Predictive Control of Microgrids (London, U.K.: Springer-Verlag, 2020). His current research interests are advanced power electronics and control to introduce energy storage technologies in the smart grid.