Stochastic Optimization of Microgrids with Hybrid Energy Storage Systems for Grid Flexibility Services Considering Energy Forecast Uncertainties

0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 1 Stochastic Optimization of Microgrids with Hybrid Energy Storage Systems for Grid Flexibility Services Considering Energy Forecast Uncertainties Felix Garcia-Torres, Carlos Bordons, Senior Member,IEEE, Javier Tobajas, Rafael Real-Calvo, Member,IEEE, Isabel Santiago and Stephane Grieu, Member,IEEE. Abstract—This paper presents a stochastic framework for the optimization of microgrids that has the functionality of providing flexibility services to System Operators (SOs) considering uncertainties in the energy forecast. The methodology is developed with the aim of being applied to complex microgrids composed of different distributed energy resources and hybrid energy storage systems (ESS). The associated optimization problem is operated in two stages: the first one performs a stochastic optimization of the microgrid in order to reserve an up/down regulation capacity with which to deal with the energy forecast uncertainties of the microgrid. The different microgrid devices are optimized by considering their operational costs in order to achieve their optimal operation in the Day-Ahead Market (DM). The second stage is used to re-schedule the initial planning according to the signal request and an economic offer from the SO. The control problem is developed using Stochastic Model Predictive Control (SMPC) techniques and Mixed-Integer Quadratic Programming (MIQP), owing to the presence of logic, integer, mixed and probabilistic variables. The simulation results show that the proposed methodology reduces the risk of undergoing up/downpenalty deviations in the Regulation Service Market (RM), also being able to provide flexibility services to the SOs, despite being subject to uncertainties in the energy forecast carried out for the microgrid. Index Terms—Stochastic Systems, Stochastic Optimal Control, Energy Storage, Optimization Methods, Hydrogen NOMENCLATURE CCapacity (Wh) CC Capital Cost ( N ) Cost Hourly economic Cost ( N /h) Manuscript received September 22, 2020; revised November 29, 2020 and March 13, 2021; accepted March 27, 2021. This work has been carried out with the financial support of the European Regional Development Fund (ERDF) under the Interreg SUDOE SOE3/P3/E0901 (Project IMPROVEMENT) program. F. Garcia-Torres and J. Tobajas are with the Applications Unit at the Centro Nacional del Hidrogeno. 13500 Puertollano, Ciudad Real, Spain. (corresponding author, phone: +34-926-420-682; e-mail: [email protected]; javier[email protected]). C. Bordons is with Laboratory of Engineering for Energy and Environmental Sustainability, Escuela de Tecnica Superior de Ingenieria, Universidad de Sevilla. 41092 Seville, Spain. (e-mail: [email protected]) R. Real-Calvo and I. Santiago are with the Electronic and Computer Department, Escuela Polit´ ecnica Superior, Universidad de C´ ordoba, Campus de Rabanales, Edificio Leonardo da Vinci, 14071 C´ ordoba, Spain (e-mail: [email protected]; [email protected]) S. Grieu is with the Processes, Materials and Solar Energy (PROMESCNRS) laboratory, University of Perpignan Via Domitia, Rambla de la thermodynamique, Tecnosud, 66100 Perpignan, France (e-mail : [email protected]) Cycles Number of life cycles Hours Number of life hours LOH Level of Hydrogen (Nm3) mMinimum possible value in the control hardware MMaximum possible value in the control hardware NsNumber of scenarios PElectric power (W) ˆ PPredicted electric power (W) SOC State of Charge (p.u) w[S]Weighting factor for the probability of the scenario zElectric power formulated as an MLD variable (W) βUnitary Benefit δOn/Off state εMinimum tolerance provided to the controller ηEfficiency (p.u) χLogical degradation state ϑMLD power variation in degradation state (W) σLogical variable start up state ΓCost of energy( N ) Subscripts: bat Battery ch Charge dis Discharge down Down-regulation elz Electrolyzer exch Exchange ext External fc Fuel Cell flex Flexibility global Global grid Main grid load Load local Local pur Purchase of energy pv Photovoltaic system rem Remaining power sale Sale of energy un Uncertainty up Up-regulation wt Wind turbine Superscripts: Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply. 0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 2 [S]Scenario DM/RM Day-Ahead/Regulation Service Market req Required sch Schedule I. INTRODUCTION THE introduction of microgrids into the wholesale electricity market has two principal interactions [1], the first of which takes place for the energy exchanges in the day-ahead and intraday electrical markets and has to be agreed upon with the Market Operator (MO). The second, meanwhile, affects RM, which is managed by SOs. The role of the SO in the wholesale electricity market is to manage the security of the power system in real time and coordinate the supply of and demand for electricity. There are reserves at the SO’s disposal that can be employed for fast load variations and unforeseen problems with production capacity. The foreseen increase in renewable power installed in the electric power system will require further resources with which to carry out the regulation between power generation and consumption. As stated in [2], the hybridization of energy storage technologies can provide additional flexibility and competitiveness to microgrids. Several energy storage options are currently available, and it is important to highlight the great potential of hydrogen as an energy carrier. As occurs with batteries, hydrogen ESS has advantages and drawbacks. Start-up/shut-down cycles, a rapid transient in the operation conditions or the number of working hours affect electrolyzers and fuel cells, while batteries are affected by the number of charge/discharge cycles or a high stress current in the charging/discharging process. A. Literature Review The control problem of microgrids ranges from their participation in electric markets (which concerns their long-term schedule for periods longer than one day) to the management of power quality (which require sampling periods of millisecond or microseconds). The computational requirements of the control platform utilized by the microgrid Energy Management System (EMS) are usually relaxed by implementing a threelevel hierarchical control structure [3]. An introduction to the importance of flexibility in electric power systems can be found in [4]. In [5], a review and classification of existing DER as flexibility providers and a breakdown of trading platforms for DER flexibility in electricity markets are presented. The authors of [6] develop a comprehensive centralized scheduling model in order to exploit demand flexibility from residential devices by analyzing the monetary value of households for demand response. The study carried out in [7] explores two model-free reinforcement learning techniques with which to schedule the operation of a battery and a heat pump in a residential microgrid as a flexibility resource. The authors of [8] propose a flexibilityoriented microgrid optimal scheduling model with which to efficiently schedule microgrid devices in order to support the SO’s flexibility requirements. In [9], the ability of microgrid components to provide flexibility is analyzed. The model developed is capable of evaluating both the impact of variable energy production and consumption and the impact of energy balancing tariffs. The work developed in [10] comprises an optimization model for EMS from the aggregator’s point of view by managing PV panels, batteries and thermal ESS. The different resources are managed in order to participate in the day-ahead energy and local flexibility markets. The resulting model is a Mixed-Integer Linear Programming (MILP) problem in which the objective is to minimize day-ahead operation costs as regards the DM and local flexibility requests. Several papers have dealt with the introduction of uncertainties into the EMS of microgrids. A chain-constrained optimization is used in [11] to optimize islanded microgrids. In [12], an innovative scheduling model with which to island hybrid AC/DC microgrids by involving sourceload uncertainties is developed. A tri-level two-stage mixedinteger robust optimization model with which to meet the diverse operational requirements of microgrids with multiuncertainties is, meanwhile, presented in [13]. In [14], an innovative decomposition Stochastic MILP is carried out. The proposed method demonstrates the effectiveness of finding the optimal battery energy storage system (BESS) power and energy sizes. The degree of complexity of developing an EMS for a microgrid that is able to interact with MOs and SOs while considering uncertainty in the energy forecast requires advanced control methodologies. Model Predictive Control (MPC) has been widely used as a framework in which to solve complex problems in industry and has recently drawn the attention of researchers in microgrids [1]. The methodology, which exists as a family of control methods, makes it possible to optimize a multi-objective cost function into which the future behavior of the microgrid devices, energy forecast and price prediction can be easily integrated. Its stochastic formulation [1] allows the management of uncertainties, while its hybrid formulation makes it possible to deal with logic and continuous variables. A two-stage chance-constrained SMPC is developed in [15] with the purpose of minimizing the cost of energy and greenhouse gas emissions in the microgrid. In [16], three SMPC strategies are compared on a hydrogenbased microgrid, while the authors of [17] develop an adaptive constrained multi-scenario SMPC that achieves a lessconservative coordination between the ESS and the introduced uncertainty of the Renewable Energy System (RES). In [18], a Bender’s decomposition-based Stochastic MILP controller is developed with the objective of minimizing the environmental and social impacts of the microgrid. The study includes a flexible load in order to reduce the requirements of load curtailment and energy generation. In the work presented in [19], the electric and thermal demands in the microgrid are considered for demand response programs using a two-stage chance-constrained SMPC. A Multi-objective Stochastic MILP is developed in [20] in order to minimize operational costs and pollution in the microgrid. In [21], scenario-based stochastic Mixed-Integer Non-Linear Programming (MINLP) is carried out, which reduces the expected costs of energy and the associated cost of emissions. Furthermore, in [22], a Dynamic Response Stochastic MILP is developed to maximize the net social benefit of the microgrid by including a responsive Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply. 0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 3 demand side in order to provide a more flexible power systems. Another approach is carried out in [23], in which a complete framework for the schedule of a hydrogen-based microgrid in order to allow it to be introduced into the electric market is developed without considering the possibility of providing flexibility services. Uncertainties in generation and demand are not considered in the schedule. The deviations in the energy forecast are corrected in a second stage, which is carried out in [24]. B. Main Contributions Despite the existence of a large number of works focusing on microgrids, the functionality of exchanging a deterministic energy profile with an external agent such as the MO in the day-ahead or the SO in real-time markets by considering a stochastic energy balance in the microgrid is an issuethat has not previously been developed and is the principal innovation of this paper. As stated in [2], the hybridization of different ESS can bring benefits to the final costs of the flexibility services. Nevertheless, there are lacks in the current state of the art as regards control strategies applied to the EMS of microgrids in order to support flexibility services by considering uncertainties, in that hybrid ESS are not employed to improve the final cost of the services provided. Very few authors have studied the potentiality of hydrogen as a technology for grid flexibility services. The proposed control structure for the hydrogen-based microgrid proposed in [23] considers a deterministic energy profile for the schedule of the microgrid, which is corrected in the real-time control of the microgrid using the controller described in [24]. This means of controlling microgrids requires the use of the hybrid ESS to absorb the uncertainties in the energy forecast, which may allow them to attain their maximum or minimum capacity when this forecast is not sufficiently accurate. This may also reduce the opportunity of using the hybrid ESS to provide flexibility services. The principal aim of this paper is to improve the development carried out in [23] and [24] in order to integrate uncertainties into the energy forecast, thus enabling the microgrid controller to exchange a deterministic energy profile with the main grid in the DM, despite the probabilistic nature of renewable sources and consumer behavior. This feature minimizes the use of the hybrid ESS to correct deviations in the energy forecast, which allows its use as a resource for the flexibility market. A second stage is, therefore, integrated into the EMS of the microgrid, which includes the functionality of changing the schedule carried out in the DM as a response to a certain power profile variation requests from the SO, including uncertainties in the energy forecast. This variation in the power profile will be linked to an economic benefit offered by SOs or agents external to the microgrid. In order to optimize the associated cost by considering uncertainties in the forecast, different scenarios in the energy forecast are considered and the optimization problem is formulated as a multi-scenario stochastic problem. Some of the features of the proposed algorithm considered as the main highlights of the present work are the following: •A multi-objective optimization problem for the tertiary control level of microgrids with hybrid ESS based on an SMPC formulation that includes the degradation mechanisms and operation costs related to the different ESS of the microgrid. •The optimal schedule of the microgrid in the DM and RM, considering uncertainties in the energy forecast of the microgrid. •The interaction with the SO for the flexibility services provided, depending on the benefit and flexibility power profile requirement offered. The different criteria integrated into the controller developed in this paper are summarized in the first column of Table I and compared to the existing literature in the other columns. As can be seen in the table, there is no reference that accomplishes all the criteria incorporated into the controller. TABLE I SUMMARY OF MAIN MICROGRID SCHEDULE METHODS COMPARED WITH THE MAIN FEATURES INTEGRATED INTO THE PROPOSED CONTROLLER SHOWN IN THIS PAPER Criterion [6]- [11]- [15]- [19]- [23], [10] [14] [18] [22] [24] Uncertainty of RESs X X X Uncertainty of demand X X X Grid Energy Exchange X X X X Flexibility Services X X Hybrid ESS X X ESS Operation Cost X C. Structure and Organization of the Paper The remainder of the paper is organized as follows: The controller is developed in Section II, in which the associated operation cost of each storage technology included in the microgrid is justified. The controller is formulated as an SMPC in order to introduce the uncertainties into the energy forecast. The results provided by the proposed controller are exposed and discussed in Section III. Finally, the main conclusions are highlighted in Section IV. The microgrid that is the object of this study is shown in Fig. 1. As can be seen, it integrates two kinds of renewable generation (photovoltaic (PV) panels and wind turbine) and two different kinds of ESSs (batteries and hydrogen). The different components of the microgrid are interconnected with the EMS of the microgrid by means of a local area network (LAN). The EMS interacts with SOs and the MO through the use of a Wide Area Network (WAN). II. SMPC-CONTROLLER DESIGN The microgrid controller is designed according to the following criteria: 1) Operational Cost of the Microgrid Components: The controller will integrate the economic costs related to lifetime, operation or degradation issues. 2) Uncertainties Integration: The controller will have the capacity to confront energy forecast uncertainties for the periods considered in the schedule horizon without acting in the RM. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply. 0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 4 Fig. 1. Microgrid and Interaction with MO and SOs 3) SO’s Signal Response: The controller will have the capacity to quantify the closest power exchange, which can be provided according to a power profile request and an economic benefit. In order to integrate the uncertainties, the control algorithm will be formulated as an SMPC problem. A stochastic framework that considers not only MPC control problems integrates uncertainties into the decision making process. With regard to MPC, these considerations affect both the dynamic model and the cost function to be minimized. [1]. Two approaches are principally used in MPC in order to integrate uncertainties into the optimization process (further details about SMPC techniques can be found in [1]): •Scenario-based approaches: The uncertainties in the system dynamic are characterized by a finite set of scenarios utilized to solve the MPC control problem. •Analytic approaches: These methods are based on probabilistic functions. Cost functions and chance are reformulated using statistical elements. A simplistic way in which to formulate the MPC controller by considering the uncertainty in the energy forecast process and avoiding complex statistical formulation could be to integrate an uncertainty band into the energy forecast, as shown in Fig. 2. A positive and negative uncertainty band (∆Pun) is added to a deterministic profile of the remaining power forecast (ˆ Prem =ˆ Ppv +ˆ Pwt−ˆ Pload), which is defined as the difference between generation and consumption, resulting in two energy scenarios: optimistic (S= 1) and pessimistic (S= 2).                     Fig. 2. Forecast remaining power considering an uncertainty band Keeping all of the aforementioned features in mind, the EMS control algorithm is structured in two steps. The first optimizes the energy exchange with the main power grid in                                                   Fig. 3. Controller Block Diagram the DM, integrating the internal operation cost of the microgrid components. It also reserves an uncertainty margin in order to prevent the microgrid from participating in the RM by using a stochastic multi-scenario formulation. The second step is based on an MPC formulation of the controller, which responds to the benefit provided by the SO according to a power exchange request with a certain profile. The controller block diagram is detailed in Fig. 3. Throughout this paper, the upper index notation [S]is used to link the different variables of the controller with their corresponding scenario. The different variables and sub-blocks that are integrated into the block diagram of the controller studied are detailed in the following sections. A. Forecast Module The forecast algorithm integrates, according to the database of meteorological measures, the array of forecast disturbances (d= [ ˆ Ppv,ˆ Pwt,ˆ Pload]) or non-manipulated variables of the microgrid (further details on the methodology employed to obtain these forecast variables can be found in reference [23]). It is composed of the energy generated by photovoltaic solar panels and wind turbine generators, along with the load consumption. A number of scenarios Nswith different probabilities can be used in order to integrate the uncertainties into these forecasts. This issue is simplified by considering just two scenarios of disturbances: assigning a positive and negative uncertainty band ∆Pun to the predicted values in the array of forecast disturbances. (d[1] = [ ˆ Ppv + ∆Pun,ˆ Pwt + ∆Pun,ˆ Pload −∆Pun]) and (d[2] = [ ˆ Ppv −∆Pun,ˆ Pwt − ∆Pun,ˆ Pload + ∆Pun]). These are denominated as the optimistic scenario S= 1 and the pessimistic scenario S= 2. Although several statistical methods can be used to calculate ∆Pun, this value is calculated using the standard deviation between the forecast remaining energy and the measured value for each hour along one complete day throughout this paper. The average value for the last year is taken by following the Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply. 0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 5 equation provided in (1). ∆Pun =1 365 1 24 day=365 X day=1 h=24 X h=1 |ˆ Prem(day, h)−Pmeas rem (day, h)| (1) where day identifies the day when the measure is carried out, hstands for the time at which the measure was carried out, ˆ Prem(day, h)is the forecast value of the remaining energy in the microgrid for each hour and day of the year, and Pmeas rem (day, h)is the measured value. The second part of the Forecast Module provides the array of predicted values for the energy prices in the schedule horizon for the processes of purchasing and selling energy with the main power grid (Γ(t) = [Γpur(t),Γsale(t)]). B. Plant Model The dynamic model of the microgrid, considering the two aforementioned energy forecast scenarios and the presence of logic and dynamic variables, is developed by using its representation in the state space shown in the expressions (2)- (4), where the state variables xare those whose values depend on their past values and on the inputs. Two kinds of input variables can be found in the state-space model of the microgrid: non-manipulated inputs or disturbances, represented by d, continuous control variables u, logic control variables δand mixed control variables z. Finally, the values of the output variables are represented by y. Note that, although two different optimization scenarios are used, the system outputs are maintained with similar values, independently of the scenario considered, see Fig. 3. x[S](t+ 1) = Ax(t) + Buu[S](t) + Bδδ[S](t) +Bzz[S](t) + Bdd[S](t)(2) y[S](t) = Cx[S](t) + Duu[S](t) + Dδδ[S](t) +Dzz[S](t) + Ddd[S](t)(3) Eδδ[S](t)+Ezz[S](t)≤Cx[S](t)+Euu[S](t)+Edd[S](t)(4) The array of continuous variables u u u[S]is defined by employing the variables contained in expression (5). u u u[S]=hP[S] dis, P[S] ch , P[S] elz , P[S] fc , P[S] pur, P[S] salei u u umin =0,0, Pmin elz , Pmin fc ,0,0 u u umax =Pmax bat ,−Pmin bat , Pmax elz , Pmax fc , Pmax grid ,−Pmin grid  (5) where P[S] ch ,P[S] dis,P[S] elz and P[S] fc are the reference power values provided by the microgrid controller to the internal controllers for the charging or discharging processes of the batteries, the electrolyzer and the fuel cell. The terms P[S] pur and P[S] sale represent the purchasing and selling power with the main grid. The array of logic variables δ[S]is composed of the components defined in (6). δ[S]= [δ[S] elz, δ[S] fc , σ[S] elz, σ[S] fc , χ[S] elz, χ[S] fc ]T(6) The electrolyzer and fuel cell have logic signals associated with their on/off-state, which are defined by employing δ[S] elz and δ[S] fc .σ[S] elz and σ[S] fc are auxiliary variables linked to the start-up processes of the electrolyzer and the fuel cell. χ[S] elz and χ[S] fc are auxiliary variables that define the instants at which the electrolyzer and the fuel cells are in on-state but the electrolyzer or the fuel cell are being switched on/off. The array of mixed product variables is formed of the variables expressed in (7). z[S]=hz[S] elz, z[S] fc , ϑ[S] elz, ϑ[S] fc i(7) where z[S] elz and z[S] fc are the mixed product of the continuous and logic control variables of the electrolyzer and the fuel cell. ϑ[S] elz and ϑ[S] fc are the power increments in the electrolyzer and the fuel cell in all those instants at which they are active but are switched on or off. The state variables of the microgrid are defined by the level of energy stored in the battery SOC[S]and in the hydrogen tank LOH[S](8). x[S]=hSOC[S], LOH[S]i(8) The outputs of the plant are defined by employing the power exchanged with the main grid P[S] grid (9). y[S]=hP[S] gridi(9) The expression provided in (2) regarding the state variables of the plant can be particularized for this controller with (10) and (11), where Cbat is the capacity of the battery, and the coefficients ηch and ηdis represent the efficiencies as regards the charging and discharging processes of the batteries. Finally, ηfc and ηelz are the efficiencies in the conversion processes carried out by the electrolyzer and the fuel cell. SOC[S](t+ 1) = SOC[S](t) +Ts P[S] ch (t)·ηch Cbat −P[S] dis(t)/ηdis Cbat !(10) LOH[S](t+ 1) = LOH[S](t) + Ts z[S] elz(t)·ηelz −z[S] fc (t) ηfc ! (11) The constraint expressed in (3) is obtained by the difference between the purchased and the sold power from/to the main grid (12), which can be obtained from the energy balance defined by the pessimistic scenario (13) and the optimistic scenario (14). P[S] grid(t) = P[S] pur(t)−P[S] sale(t)(12) P[1] pur(t)−P[1] sale(t) + P[1] dis(t)−P[1] ch (t) + z[1] elz(t) −z[1] fc (t) + ˆ Ppv(t) + ˆ Pwt(t)−ˆ Pload(t) + ∆Pun(t) = 0 (13) P[2] pur(t)−P[2] sale(t) + P[2] dis(t)−P[2] ch (t) + z[2] elz(t) −z[2] fc (t) + ˆ Ppv(t) + ˆ Pwt(t)−ˆ Pload(t)−∆Pun(t) = 0 (14) The expressions provided in (4) can be obtained by transforming the logic relationships between the variables u,δand z, following the procedure defined in [25]. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply. 0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 6 Low power values in the electrolyzer and the fuel cell can lead to degradation aspects [23]. The manufacturers of these devices, therefore, recommend switching in a value lower than Pmin elz and Pmin fc . This aspect is introduced into the controller by means of the mixed product defined in (15), which can be transformed into the linear constraints defined in (16)-(19), but replacing αby {elz, fc}. z[S] α(t) = P[S] α(t)·δ[S] α(t)(15) m≤zα(t)−Pmax αδ[S] α(t)≤0(16) 0≤z[S] α(t)−Pmin αδ[S] α(t)≤M(17) P[S] α(t)−Pmax α(1 −δ[S] α(t)) −z[S] α(t)≤0(18) 0≤P[S] α(t)−Pmin α(1 −δ[S] α(t)) −z[S] α(t)(19) where mand Mrepresent the minimum/maximum possible value in the control hardware and εrepresents the tolerance provided to the controller. The start-up states σ[S] elz and σ[S] fc can be related to δ[S] elz and δ[S] fc with the logic relationships expressed in (20), which can be linearized as constraints with the expressions (21)-(23). σ[S] α(t) = δ[S] α(t)∧ ∼ δ[S] α(t−1)|α=elz,fc (20) m≤ −δ[S] α(t) + σ[S] elz(t)≤0|α=elz,fc (21) m≤ −(1 −δ[S] α(t−1)) + σ[S] α(t)≤0|α=elz,fc (22) m≤δ[S] α(t) + (1 −δ[S] α(t−1) −σ[S] α(t)≤1|α=elz,fc (23) where the symbols ∧and ∼represent the logic operations AND and NOT, respectively. The relation between χelz(t)and δelz(t)is defined in (24), which can be introduced into the controller with the constraints (25)-(27): χ[S] α(t) = δ[S] α(t)∧δ[S]aα(t−1)|α=elz,fc (24) m≤ −δ[S] α(t) + χ[S] α(t)≤0|α=elz,fc (25) m≤ −δ[S] α(t−1) + χ[S] α(t)≤0|α=elz,fc (26) m≤δ[S] α(t) + δ[S] α(t−1) −χ[S] α(t)≤1|α=elz,fc (27) Finally, ϑα(t)is defined as a mixed product defined in (28), which is integrated into the controller by employing the linear constraints provided in (29)-(32), where α={elz, fc}. ϑα(t)=(Pα(t)−Pα(t−1)) ·χα(t)(28) m≤ϑα(t)−∆Pmax αχα(t)≤0(29) 0≤ϑα(t)−∆Pmin αχα(t)≤M(30) m≤ϑα(t)−∆Pα(t) + ∆Pmin α(1 −χα(t)) ≤0(31) 0≤ϑα(t)−∆Pα(t) + ∆Pmax α(1 −χα(t)) ≤M(32) C. Day-Ahead Market MPC Controller The stochastic and hybrid MPC-controller is based on minimizing the cost function expressed in (33) for a number Nsof considered scenarios. Note that y(t), as output of the system, has to be considered as a deterministic variable, regardless of the scenario, and its value must be similar in all of them. The cost function expressed in (33) is integrated into the controller for a schedule horizon of 24 hours. Note that the DM closure is carried out at 10.00 a.m. on the day before it is executed. The schedule horizon is, therefore, between the instants SH0= 14 and SH = 38. The notation tk, which stands for t+k, is used for reasons of simplification. As can be seen in (33), each term is underbraced, indicating its contribution to the cost function. In (33), Nsidentifies the number of scenarios considered, where w[S]is the weighting factor corresponding to the assigned probability for each scenario. The notation CC is used for the capital cost of acquisition for each device. The batteries have a number of life cycles (Cycles) as provided by the manufacturer. This value is included in order to quantify the cycling cost of the batteries. As specified in [23], high current values for the charging and discharging processes can lead to degradation issues in batteries. These issues are incorporated into the term related to the degradation of the batteries, where the cost of these aspects is quantified with Costdegr,α. The electrolyzer and the fuel cell have a limited number of working hours Hours, and an hourly operation cost associated with their maintenance Costo&m,α. Fuel cells and electrolyzers have degradation mechanisms associated with their start-up cycles, but also incur in aspects of degradation when fluctuations are applied to their power reference once they have been started up. These aspects are included in the terms concerning the Hydrogen ESS Degradation. The last terms in (33) are included in order to maintain the energy stored in each ESS in a reference value, which confronts future uncertainties in the next schedule horizon. Jlocal(t) = Ts Ns X S=1 w[S] k=SH X k=SH0−ΓDM sale (tk)z[S] sale(tk)+ΓDM pur (tk)z[S] pur(tk) | {z } Grid Exchange Revenue&Cost(Jgrid) +CCbat 2·Cyclesbat X α=ch,dis (z[S] bat,α(tk)) | {z } Batteries Cycling Cost +X α=ch,dis (Costdegr,α ·(z[S] bat,α(tk))2) | {z } Batteries ESS Degradation +X α=elz,fc  CCα Hoursα +Costo&m,αδ[S] α(tk) | {z } Hydrogen ESS Hourly Cost Use +Coststartup,α ·σ[S] α(tk) + Costdegr,α ·(ϑ[S] α(tk))2 | {z } Hydrogen ESS Degradation +wSOC ·(SOC[S](tSH )−SOCref )− | {z } F uture Uncertainties +wLOH ·(LOH[S](tSH )−LOHref )− | {z } F uture Uncertainties (33) subject to: xmin ≤x[S]≤xmax;umin ≤u[S]≤umax 0≤δ(t)≤1; 0 ≤z[S](t)≤zmax ymin ≤y(t)≤ymax;y[1] =y[2] =... =y[Ns] Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply. 0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 7 D. Flexibility Market MPC Controller One of the design criteria for the EMS is the ability to supply flexibility services. The SO may require the microgrid EMS to provide a variation in the energy exchange, with the main power grid as a flexibility demand, at any sample instant. The nomenclaure local is utilized to represent the microgrid acting as a single system that exchanges energy only with the grid in the DM, while the terminology global represents the operation point at which the microgrid EMS responds to a request from the SO, who offers an economic benefit. The day-ahead optimization by the EMS is carried out at 10 a.m. on the day before the 24 hour-optimization period is executed. The uncertainties in the energy prediction can be considered with a lower uncertainty band once the optimization process is being carried out, owing to the reduction in the schedule horizon of the EMS. It is in this time scale that the flexibility services are demanded by the SO. At the sample instant at which the SO makes its request, the solver calculates its optimal operation, providing the following set of optimized control variables [u[S],opt local (t),δ[S],opt local (t)] as a solution, which are associated with an optimal value of Jlocal. Note that, if the term of the cost function underbraced as Grid Exchange Revenue&Cost is executed after the DM closure, the expression (34) has to be replaced. As can be seen, the RM penalties are included in expression (34) if more energy than that scheduled is exchanged (δup = 1) or less energy than that planned is acquired (δdown = 1). The prices associated with these penalty costs are the constants Γup and Γdown. Jgrid = k=SH X k=SH0X α=up,down ΓαPgrid(tk)−PDM,sch grid (tk)δα | {z } Regulation Market P enalties (34) Jglobal =Jlocal + Ns X S=1 k=SH X k=SH0 wgrid P[S] up (tk) + P[S] down(tk)−Preq SO (tk)2 | {z } System Operator Exchange Revenue (35) Jlocal(u[S],opt global ,δ[S],opt global −Jlocal(u[S],opt local ,δ[S],opt local )) + k=SH0 X k=SH βup(tk)Pup(tk)−βdown(tk)Pdown(tk)≤Bext P[1] α(tk) = P[2] α(tk)|α=up,down;P[S] up ≥0; P[S] down ≤0 (36) The global cost function (35) can be formulated by adding a term to the cost function in order to optimally track the power profile requested by the SO Preq SO and a constraint to guarantee that the cost of satisfying this energy provided to the SO is lower than the economic benefit provided by it. The constraint shown in (36) has to be added to the controller. The terms P[S] α(tk)|α=up,down are related to the energy provided to the SO as a flexibility service. βup(tk)and βdown(tk)are the unitary benefit per Wh of flexibility energy provided, while Bext represents the external benefit offered by the SO. The upper and lower limits for P[S] down and P[S] up are ”0” for those instants at which there is no flexibility requirement from the SO. Note also that they are considered as system outputs and that their values are similar for all the scenarios considered, in order to manage the risk related to uncertainty in the energy forecast. The flexibility services can be activated at any instant of the schedule horizon. In order to compare the variation in the cost function, the same number of sample instants and schedule horizons as those used in the DM optimization are utilized. Nevertheless, the values of the optimization variables have to be imposed as constraints in these sample instants, which take place prior to those at which the SO communicates request for a flexibility service to the microgrids. The final exchange with the main grid in this step of the algorithm is defined by the following expression: Pgrid(t) = PDM,sch grid (t) + Pup(t) + Pdown(t)(37) TABLE II VALUES OF THE CONTROLLER Renewable Energy PV Panels: 30 kWp, Wind Turbine: 10 kW Hydrogen ESS Electrolyzer: 30 kW, Tank: 35 Nm3, Fuel Cell: 20 kW Costdegr,elz = 0.0577 N /W, Hour= 10000 h, LOHmax = 35Nm3LOHmin = 5Nm3 ς= 0.23 Nm3/kW h, CC= 8.22 N /kW, Costo&m,elz = 0.002 N /h Coststartup,elz = 0.123 N Costdegr,fc = 0.0018 N /W, Hours= 10000 h ς= 1.320 kWh/Nm3, CC=30 N /kW, Costo&m,f c = 0.001 N /h Coststartup,fc = 0.01 N , N ,wLOH = 1 Batteries ESS: 15 kW, 55 kWh ηch = 0.90,ηdis = 0.95, CC=125 N /kWh, Cycles=3000, wSOC = 10 Costdegr,dis = 10−9 N /W2h, Costdegr,ch = 10−9 N /W2h SOCmax = 1 SOCmin = 0.25 Data based on reference [23] Fig. 4. Day-ahead Energy Prices Prediction III. RESULTS In this section, the controller is developed and validated by means of numerical simulations carried out using the TOMLAB® as solver tool. The execution time was 10.78 s for an uncertainty band of ±2000 Wand 14.54 s for an uncertainty band of ±5000 W. The computation was carried out by employing a PC with an Intel® Core™ i7-4510U @ 2.60 GHz as a computer processor. The values integrated into the controller are shown in Table II. The sample time employed is Ts= 1 hour and the schedule horizon is 24 hours, Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply. 0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 8 corresponding to a complete day. Although the method can be used by integrating a higher number of considered scenarios and by assigning different probabilities through the use of the weighting factor w[S], here, the scenarios used are Ns= 2, and a similar probability of w[S]= 0.5is applied, for the sake of simplicity. Fig. 5. Validation of the proposed schedule algorithm in different scenarios A. Day-Ahead MPC Controller The energy prices predicted using the forecast algorithm for energy selling in the DM are shown in the graph depicted in Fig. 4. These prices are based on the Iberian Market Operator history data executed for 05/04/20. It has been considered that ΓDM pur = 2ΓDM sale . The benefit required for flexibility services is considered to be βup = 3ΓDM sale and βdown = 1.5ΓDM sale . The different scenarios of the remaining power (Prem) in the microgrid, defined as the difference between generation and consumption, for which the algorithm will be validated are shown in Fig. 5. The base forecast scenario is labeled with legend P[0] rem. The optimistic and the pessimistic scenarios are calculated in the base scenario by adding an uncertainty band of ∆Pun =±5000Wto the base scenario. The algorithm presented in Section II is employed in order to execute these optimistic and pessimistic scenarios (P[1] rem(t) = P[0] rem(t) + ∆Pun, P[2] rem(t) = P[0] rem(t)−∆Pun). Scenarios [3] and [4] are defined by employing the following mathematical expressions P[3] rem(t) = P[0] rem(t) + (−1)t·∆Pun ∀t∈[1,24] and P[4] rem(t) = P[0] rem(t)−(−1)t·∆Pun ∀t∈[1,24]. Scenarios P[5] rem and P[6] rem are generated by adding ±4000Wto the scenario P[0] rem and changing the uncertainty-band sign only at the sample instant corresponding to 12 a.m. The multi-scenario algorithm is executed with the corresponding scenarios P[1] rem and P[2] rem. The remaining scenarios are used only in order to validate the algorithm. A similar algorithm is used for the other scenarios, but considering only the corresponding validation scenarios (S= 0, S = 3 −6), in which the values for the power exchange with the grid at all the sample instants are imposed as constraints with similar values to those obtained after executing the multi-scenario algorithm applied to S=1 and S=2. As can be seen in the aforementioned figure, despite obtaining different energy forecasts using the flexibility provided by the hybrid ESS, it is possible to exchange the power profile obtained for S= 1 and S= 2 with the main power grid for all the sample instants in the remaining scenarios. It will be also observed that, independently of the scenario, the operation of the hybrid ESS is carried out, and their associated cost is optimized, the number of cycles used for the batteries is minimized, and the values of the charging and discharging power applied to the batteries are minimized. The optimization of the fuel cell and the electrolyzer is also carried out, thus minimizing their operation cost. The number of working hours is, therefore, optimized, as are the start-up and shut-down states. The fluctuating conditions associated with these devices are also avoided. The exchange of energy with the main grid is also carried out at those times at which the forecast prices are most beneficial to the microgrid. B. Flexibility Market MPC Controller The results obtained for the second case study, concerning the up-regulation and down-regulation offered to the SO as a flexibility service, are shown in Fig. 6 and Fig. 7. In the first simulation, shown in Fig. 6, the SO requires a certain power profile for the sample instants corresponding to 11-13 hours. Note that until the instant of the requirement, the previous sample instants are considered as past sample instants and both remaining power scenarios are, therefore, similar. After the sample instant at which the down-regulation is requested Ts = 10 a.m., an uncertainty band of ±2000 W is introduced. The four graphs in the upper half of Fig. 6 represent the exchange when a benefit B= 3 N is sent. The operation cost of the microgrid acting as a single system for these scenarios is, respectively, 64.53 N for the optimistic scenario and 63.04 N for the pessimistic scenario (using the cost factors shown in Table II). The costs of carrying out the down-regulation are 70.57 N and 67.64 N , and the EMS will, therefore, agree to carry out the down-regulation in its totality for offers higher Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply. 0885-8950 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TPWRS.2021.3071867, IEEE Transactions on Power Systems IEEE TRANSACTIONS ON POWER SYSTEMS 9 Fig. 6. Down-Regulation with different offers from the SO Fig. 7. Up-Regulation under different offers from the SO than 6.07 N , as can be observed in the four last graphs of Fig. 6. A second simulation for the up-regulation can be observed in Fig. 7, when the SO demands a requirement at 19.00. The four upper graphs show the maximum regulation that can be carried out by the microgrid EMS because the minimum levels of energy stored have been reached. The last four graphs displayed in Fig. 7 show the regulation results for an offer of 4 N . As can be seen in Fig. 6 and Fig. 7, the criteria of Fig. 8. Optimization of the Microgrid in the DM using the proposed methodology in [23] Fig. 9. Optimization of the Microgrid in the RM using the proposed methodology in [23] after a change in the energy scenario minimizing the operation cost of the ESS are also followed by the controller, despite the modification made to the energy exchanged with the main grid. Note that, although the sample time used for the flexibilty services for the sake of simplicity is also Ts= 1 hour, it is the same as that used in the day-ahead market optimization. This sample period could be reduced in the case of flexibility services by simply considering that the values shown in Table II are based on this sample period, and their proportional values should, therefore, be considered for lower sample periods than that selected in all the terms of the cost function with the exception of those related to the start-up states, which would be similar to the independence of the sample period selected. C. Comparison with previous existing methods in literature In order to demonstrate the potentiality of the proposed controller, it is compared with the methodology proposed in [23], in which only an energy scenario was used. Fig. 8 shows the results of this controller for the scenarios considered, S= 1 and S= 2. As can be seen, there is a different schedule for the energy exchanged with the main grid. Fig. 9 shows the results of the optimization of the microgrid if a change in the scenario forecast occurs. The schedule of the energy exchanged with the main grid in S= 1 is, therefore, imposed as a constraint to the optimization for S= 2 and vice-versa. As what is stored in the ESSs is not sufficient to carry out Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 28,2021 at 09:33:09 UTC from IEEE Xplore. Restrictions apply.