Simplified Event-Based Load Shedding Scheme for Frequency Stability in an Isolated Power System With High Renewable Penetration. El Hierro: A Case Study

This study presents a novel approach to calculate the load to be shed in El Hierro isolated power system in generation tripping events. The proposed shedding law is based on a linear regression model. The regression model is obtained by means of offline dynamic simulations in a set of representative test cases. The proposed shedding law has been compared to the under-frequency load shedding scheme currently utilized in El Hierro, considering three different implementation schemes. The results of this study demonstrate that the proposed shedding law can contribute to reduce the average amount of load shed in generation tripping events and meet the requirements of the system’s frequency with a moderate investment in smart power management devices.


INTRODUCTION
Frequency stability is a capital issue in isolated power systems with high penetration of renewable generation. Power imbalances usually produced by the outage of a single generation unit may produce substantial frequency deviations (Sigrist et al., 2018). Two are the main causes that explain this: the limited system's rotating inertia in such small power systems and the fact that each generation unit generates an important percentage of the total demand, so the outage of any generation unit provokes an important imbalance. Since most variable renewable generation units are connected to the grid through power converters, substituting conventional thermal generators by, e.g., wind power plants contributes to decreasing the system's rotating inertia. Moreover, in many cases, non-dispatchable energy sources do not contribute to primary regulation reserve (Hafiz and Abdennour, 2016). This kind of problem does not often affect large interconnected power systems thanks to their strong power-frequency relation (ENTSO-E, 2019).
Theoretically, primary frequency control should arrest frequency deviations caused by power imbalances. However, this control action is in some cases not able to maintain the frequency within acceptable limits after large disturbances in isolated power systems (Padrón et al., 2016). In these cases, in order to avoid severe damage in rotating machines associated with low frequency levels, the frequency decay is stopped by load disconnection or load shedding, with subsequent consequences on the economic activity, among others (Delille et al., 2012). For example, from 2005 to 2010, in the seven Canary Islands, there were more than 200 generation tripping events per year and load shedding was activated in 50 of these events (Padrón et al., 2016).
Automatic Under-Frequency Load Shedding (UFLS) schemes are widely adopted in isolated power systems (Sigrist et al., 2018). Static UFLS schemes (as referred to in (Sigrist et al., 2018)) consist of continuously measuring the system's frequency and shedding a predefined amount of load when the frequency falls below a certain threshold. Though widely used, static UFLS schemes have an obvious limitation since, as it is well known, it is not possible to estimate the exact load that must be shed only from the system's frequency. Note that the same frequency deviation can be a consequence of diverse incidents and can occur for different system's operating conditions. In the 80's and 90's of the past century, there were several attempts to enhance static UFLS schemes by the use of the Rate of Change of Frequency (RoCoF) as an input variable to the UFLS scheme logic (Girgis and Ham, 1982), (Thompson and Fox, 1994). The initial RoCoF allows qualitatively estimating the severity of the incident that has caused the frequency decay and quantitatively estimating (in an approximate manner) the frequency NADIR from the system's inertia constant. However, as demonstrated in (Anderson and Mirheydar, 1992), the RoCoF is subject to significant distortions by local dynamics. In addition, the inertia constant is continuously changing, as discussed in (Chuvychin, et al., 1996). To mitigate the RoCoF distortions, a multi-stage RoCoF filtering procedure is proposed in (Sodin, et al., 2020), enabling monitoring RoCoF in real-time for UFLS purposes. In order to accurately estimate the frequency of NADIR, it would be necessary to consider as well the sensitivity of the loads to the system's frequency, the settings of the speed governors of the spinning units (Elkateb and Dias, 1993), which also change pretty often. UFLS schemes using RoCoF relays are referred to in (Sigrist, et al., 2018) as semiadaptive. UFLS schemes using the response of the system (frequency, RoCoF, voltage, etc.) to the contingency as input variable are referred to in (Seyedi and Sanaye-Pasand, 2009) as response-based schemes. Event-based Load Shedding (ELS) schemes differ from UFLS schemes in the variable which activates the load shedding and, consequently, the criteria used to define the amount of load to be shed. In these load shedding schemes, the state of specific elements of the system, such as, e.g., the status of certain breakers, are used to trigger the load shedding and to determine the load to be shed so as to get ahead of very severe disturbances (Dai et al., 2012). In ELS schemes, the corrective action is usually taken almost immediately after the event has occurred, and there is no need to wait until the frequency falls below a certain value. The main benefits of ELS schemes can be summarized as follows: 1) they are much faster than static UFLS schemes and 2) the control efficiency is stronger. The RoCoF in isolated systems with high penetration of renewable generation is very high and hence the quickness of the control action is key to increasing load shedding effectivity and to minimizing the amount of load to be shed.
The most important part of ELS schemes is a decision table (Dai et al., 2012), where the amount of load to be shed is determined from the measured operational variables. This decision table can be defined a priori by the offline simulation of a complete set of representative scenarios or in real-time by the use of dynamic simulations.
The main contribution of this study is to propose a novel approach to determine the amount of load to be shed in El Hierro power system in generation tripping events. El Hierro is a very special isolated power system in the Canary Archipelago. With the aim of becoming free from greenhouse gases emissions, a hybrid wind pumped-storage hydropower plant was committed in June 2014. During last years, the share of electricity supplied by this particular power plant has monotonically increased. So much so that during the 2019 summer, for 596 h, the electricity demand was continuously and exclusively supplied by renewable energy resources. The increase in the penetration of renewable sources in the El Hierro power system supposes, however, an important challenge for the system's frequency control. In (Girgis and Ham, 1982), it is reported that load shedding is needed not only when a generator outage occurs but also to cope with abrupt wind speed down ramps. According to (Taveira, et al., 2018), this kind of incident gives rise to several load shedding events per month. In this study, it is proposed to use a set of offline simulations to determine a linear regression model, which gives the amount of load to be shed from a set of operational variables so that the system's frequency meets the operational requirements imposed by the TSO. The novelty of the proposed approach does not lie only in the particularity of the power system where it has been applied. To the authors' knowledge, there are several examples of applying event-based load shedding in the literature (Dai et al., 2012;Xu et al., 2017;Seyedi and Sanaye-Pasand, 2009;Shi et al., 2019). All of them are applied to very big power systems compared with El Hierro one. In addition, all of these publications propose the use of some kind of optimization method for determining the optimal load shedding. Some approaches utilize, among others, interpretable machine learning methods, such as decision trees (Genc et al., 2010). Although such rule-based approaches are effective, there remain some challenges. For the online ELS control process with limited response time after the contingency occurs, the additional  computational burden will be introduced by the optimization part (Li et al., 2021). These computation processes require hundreds of milliseconds, which are manageable in such complex systems. However, this is not a proper solution in small isolated systems, where RoCoF can overpass 1.5 Hz/s when a sudden disconnection of a generation unit takes place (Sarasúa et al., 2019). As far as we know, the approach we propose has not been previously proposed in any other study. The effectiveness of the regression model is demonstrated by means of a set of out-of-sample simulations of generation tripping events. All simulations are performed assuming the load is shed in a single step, analogously to (Tofis, et al., 2017).
The proposed approach is fully aligned with the findings of previous research articles, as summarized below.
-The simulation model used to perform the offline simulations is based on the well-known swing equation (Terzija, 2006;Mohamad et al., 2013;Mohamad et al., 2014;Alhelou et al., 2019). -The regression model does not consider a single local measurement of the system's RoCoF, but the RoCoF of the so-called center of inertia (RoCoF-COI) (Seyedi and Sanaye-Pasand, 2009;Rudez and Mihalic, 2011a;Njenda et al., 2018). -The regression model does not use only the RoCoF-COI as an independent variable but also other complementary measurements (Rudez and Mihalic, 2011b). -The simulation model considers the deployment of the frequency containment reserve (Concordia et al., 1995;Lin et al., 2008;Karimi et al., 2012;Chang-Chien et al., 2012;Rudez and Mihalic, 2011c;Rudez and Mihalic, 2016;Parizad et al., 2017). -The simulation model considers the sensitivity of the loads to changes in the system's frequency (Terzija, 2006;Lin et al., 2008;Rudez and Mihalic, 2016;Parizad et al., 2017). -The regression model is linear so as to reduce the processing time that can be critical in weak power systems such as the one of El Hierro island (Seyedi and Sanaye-Pasand, 2009;Karimi et al., 2012;Rudez and Mihalic, 2016).
It is important to note that the proposed approach can be easily replicated in similar isolated hybrid power systems currently under development in, e.g., the Aegean islands or the Azores archipelago. The proposed shedding law has been compared with the shedding law implemented in the current UFLS scheme by means of dynamic simulations carried out in MatLab Simulink.
This study is organized as follows. In the Description of El Hierro Power System, the main characteristics of El Hierro power system are presented. The proposed approach is described in the Proposed Approach. In the Discussion of Results, the effectiveness of the proposed approach is demonstrated by means of simulations. Finally, the main conclusions of this study are drawn in Section 5.

DESCRIPTION OF EL HIERRO POWER SYSTEM
El Hierro is an island belonging to the Canary Islands archipelago, which was declared as a biosphere reserve by UNESCO in 2000. The power system of El Hierro comprises 11.18/11.2/11.5 MW of diesel/hydro/wind generation. The power system also includes a pump station equipped with both variable and fixed speed pumps (Sarasúa et al., 2019). Figure 1 shows the simplified single line diagram of the El Hierro power system.

PROPOSED APPROACH
This study proposes the use of a linear regression model to determine the load to be shed in the El Hierro power system in generation tripping events. The regression model is obtained from a set of offline representative simulations.
The block diagram of the simulation model used to perform the simulations is shown in Figure 2. The simulation model is a single busbar model. The model calculates the frequency (f) from Eq. 1 in p. u. values.
The frequency variation described by Eq. 1 is the result of the imbalance between the power supplied by the wind (p w ), diesel (p d ), and Pelton units (p h ) and the power consumption of the pumping units (p p ) and loads (p dem ). The sensitivity of the loads consumption to frequency variations is considered by means of the term Dnet. Eq. 1 is computed in the Power System Block in Figure 2. Load demand and wind power are both input variables of the simulation model. The system's inertia is represented in Eq. 1 by H d and H h , which correspond to the inertia constant of the spinning diesel and Pelton units, respectively.
The diesel generation units are modeled by means of two transfer functions, as shown in Figure 3. The reader is referred to (Kundur, 1994) for details about the model of the diesel units. The equations used to model the hydropower plant and the pumping station have been used previously by the authors in . More precisely, equations 2-12 from  have been used in this study. The pump station and the hydropower plant share the upper and lower reservoirs, but each hydraulic circuit uses a different pipe. Because of the length of both conduits, an elastic water column model is required for modeling their dynamic response. Authors have adapted the lumped parameters approach to the penstocks characteristics (Martínez-Lucas et al., 2015). The pump station model considers the eight hydraulic pumps and the hydropower plant model considers the four Pelton units and their corresponding joints and manifolds. In addition, the influence of the frequency variations in pumps power consumption is considered in the model .
Frequency deviations are initially arrested by the online generation units accordingly to their respective speed droop (primary frequency regulation). Then, the frequency is restored to 50 Hz by the secondary regulation action, which is commanded by the system's Automatic Generation Control (AGC). The AGC "continuously" sends power set-point signals to the online generation units (Martínez -Lucas et al., 2020). The AGC system is modeled as in (Pérez-Díaz, et al., 2014), specifically equations (22)-(25). It is supposed that all online generation units analogously participate in secondary regulation and that their participation factors are obtained as a function of each unit's speed droop.
Although modern wind generators and variable-speed pumping units can contribute to frequency regulation and/or provide inertia emulation , in this study, no contribution to frequency regulation or inertia emulation from the wind generators or the variable-speed pumping units is assumed.
Electromagnetic transients are supposed to be very rapid compared to the dynamics of interest for frequency control and are therefore neglected.
As is stated in (Xu, et al., 2011), the effectiveness and robustness of ELS schemes are closely related to the proper selection of representative scenarios. In order to make a proper selection of the scenarios for the offline simulations, the historical data of electricity demand and wind power generation in El Hierro corresponding to 2018 have been divided, respectively, into four and six intervals (D i , i 1,. . .4; W j , j 1,. . .,6). These data have been taken from (Red Electrica de España, 2019) and have a 10 min time resolution.
The 10 min generation data available in (Red Electrica de España, 2019) are grouped by technology. In order to define a set of representative scenarios of hourly generation schedules (including the status and power consumption of the pumping units), we have used the unit commitment model presented in (Fernández-Muñoz, et al., 2019) (referred therein to as DAUC). With the help of this model and using as input data the historical hourly demand and wind power corresponding to 2018, a large database of 2,976 realistic hourly generation schedules have been built. We have then grouped these generation schedules according to the interval D i -W j they "belong to." For each D i -W j interval, we have then selected the five most representative hourly generation schedules, comprising only diesel units (DIE), only Pelton units (HYD), diesel and Pelton units (DIE + HYD), diesel and pumping units (DIE + PUMP), and Pelton and pumping units (HYD + PUMP). This latter operation mode is referred to as hydraulic short-circuit operation in (Fernández-Muñoz and Pérez-Díaz, 2020). Not all above combinations (DIE, HYD, DIE + HYD, DIE + PUMP, and HYD + PUMP) are present in the hourly generation schedules provided by the unit commitment model for every D i -W j interval. For this reason, the total number of scenarios is 60 instead of 120 (24x5). For the selection of these 60 hourly generation scenarios, we have used clustering techniques as suggested in (Sigrist, et al., 2010). More precisely, we have used the well-known K-means algorithm in Matlab.
Once the scenarios were defined, we used the above-described simulation model to simulate for each scenario the sudden loss of the generating unit injecting the highest value of active power to  September 2021 | Volume 9 | Article 698081 the system. By means of an iterative procedure, we have calculated for each scenario the exact amount that should be shed 200 ms after the sudden loss of the generating unit for the system's frequency to meet the requirements specified in (Sigrist, 2015). According to (Sigrist, 2015), the frequency in the Spanish insular power systems must not fall below 47 Hz and not stay below 48 Hz for more than 2 s. The load to be shed calculated in this way is referred to hereinafter as P ex shed . The time considered for the load shedding to be realized has been previously considered in other articles (Dai et al., 2012) and is compatible with currently available industrial power management devices (ABB Group, 2019).
Then, a linear regression model (2) has been obtained to approximate the "upper envelope" of P ex shed (hereinafter referred to as P env shed ) as a function of the lowest number of system variables x that best capture the behaviour of P env shed .
The coefficients of the regression model have been obtained from the solution of a constrained least squares problem. The difference between P ex shed and P env shed has been minimized subject to In order to determine the independent variables x of the regression model, a trial and error approach has been followed. The initial candidate variables were chosen consistently with the literature review summarized in Introduction and the authors' experience. The initial candidate variables considered in this study were H t 0 − and H t 0 + , which represent the system's inertia right before and after the tripping event, respectively, P p , P h , P fsp(t 0 − ) (which represents the power consumption of the fixed speed pumps right before the tripping event) and B p and B d , (binary variables which indicate whether one or more Pelton units and Diesel units, respectively, are operating, 0 if yes, one if not). The reader is referred to Regression Model for more details on the regression model. Finally, the response of the system's frequency using Eq. 2 as load shedding law has been compared to that obtained with the current UFLS scheme in the El Hierro power system, both in the 60 above-mentioned scenarios and the other scenarios provided by the unit commitment model. The current UFLS, which is summarized in Table 1, is a conventional step-wise shedding scheme.
Three different implementation schemes IS1-3 for the proposed load shedding law have been considered, namely: -IS1. The load provided by Eq. 2 is shed 200 ms after the generator outage. This implementation scheme would require the use of a Wide Area Monitoring, Protection and Control (WAMPAC) system (Terzija et al., 2010) and can be considered an optimistic implementation scheme. -IS2. The load is assumed to be evenly distributed across four Medium Voltage (MV) feeders responding to relay tripping signals. P The implementation schemes IS2 and IS3 are more realistic than IS1, even though the assumption that the load consumption is evenly distributed across the system's feeders can be criticized as discussed in (Sigrist and Rouco, 2014); (Potel, et al., 2019). Nonetheless, it should be noted that the simulations of the current UFLS scheme have been performed using an analogous assumption. IS2 is fully compatible with the current structure of El Hierro's 20-kV MV distribution network as described in (Cabrera Quintero, et al., 2016). The entire electricity consumption in the system flows through four 20-kV overhead lines connected to the substation of the diesel power plant.
The methodology followed in this study is summarized in Figure 4.

Regression Model
Our first attempt was to find a linear regression model that approximates P env shed as a function of the system's RoCoF, following the recommendations of (Sigrist, 2015). However, as can be seen in Figure 5, the particular characteristics of the El Hierro power system make it difficult to establish an accurate linear relationship between the system's RoCoF and P ex shed or P env shed . Our second attempt was to use the power imbalance (hereinafter referred to as P imb ) as a unique independent variable. However, as can be seen in Figure 6, it is not possible either to establish an accurate linear relationship between P imb and P ex shed or P env shed . The lack of linear relationships is probably due to the important dynamic differences between the generation mixes present in the El Hierro power system. Graphical simulation results, Figures 8, 9 (see In-Sample Test Cases), confirm this hypothesis.
By means of a trial and error approach, we have obtained a 4-variable linear regression model (6), where H t 0 + , P fsp(t 0 − ) and B p . The use of these regression variables demands proper communication between generation units and pumps and the control center. Furthermore, the correlation coefficient R between these variables is 0.9564 and the p-value is near zero. The average value of the sixty residuals resulting from the linear regression is 0.182 MW and the standard deviation is 0.159 MW. The graphical results of the regression model in the 60 representative scenarios are shown in Figure 7.

Implementation Results
In-Sample Test Cases Table 2 summarizes the average results obtained in the 60 representative scenarios when P ex shed , P env shed , P 4f shed and P 8f shed are shed 200 ms after the trip of the generator with the highest power output, as well as with the current UFLS scheme. As can be seen in Table 2, the proposed shedding law, with the three implementation schemes, helps reduce considerably the average amount of load shed after the trip of the generator with the highest power output in comparison to the current UFLS scheme.
The results obtained for each scenario are included in Table 3. As can be seen in Table 3, with the current UFLS scheme, the frequency meets the requirements specified in (Sigrist, 2015) in 58 out of the 60 scenarios, whereas with the proposed shedding law and the three implementation schemes, the system's frequency meets these requirements in all scenarios. Figures 8, 9 show simulation results corresponding to two different scenarios. The first one, 51 in Table 3, presents an intermediate demand situation (4.5 MW) with high wind power participation (7.1 MW), two diesel units, two fixed speed pumps, and two variable-speed pumps operating. One of the five wind  turbines, which is initially supplying 1.42 MW, trips at t 10 s P ex shed 0.39 MW are in this case necessary for the frequency to stay below 48 Hz for exactly 2 s. In the second example, 49 in Table 3 represents a high demand situation (6.1 MW), with an intermediate wind power participation (6.7 MW). One of the five wind turbines, which are supposed to supply the same power (1.34 MW), trips in this case at t 10 s P ex shed 0.55 MW are necessary for the frequency to stay below 48 Hz for exactly 2 s.
The differences in the dynamic response between the Diesel (thermal) and Pelton (hydroelectric) units can be observed by comparing Figures 8, 9 to each other. The high water inertia in the hydropower conduits makes the provision of frequency containment reserves by the hydropower plant slower than that of the diesel units. This explains the need to use B p as an independent variable of the regression model (6).
The influence of fixed speed pumps in frequency regulation is visible since the power consumed by the pump station is an important percentage of the total demand in many scenarios and the asynchronous machines reduce their energy consumption when the frequency decay.
As can be seen in Table 3, the proposed shedding law, with the three considered implementation schemes, also contributes to increasing the frequency of NADIR (in 80% of the scenarios with IS1 and 100% with IS2 and IS3). Taking into account that with the proposed shedding law and the three implementation schemes IS1, IS2, and IS3, the load shed is lower than or equal to the one corresponding to the current UFLS scheme in 54 out of the 60 analyzed scenarios, it is obvious that the improvement in the frequency NADIR is mainly due to the speed with which the corrective action is executed.
The proposed load shedding is executed in all scenarios 200 ms after the trip of the generator with the highest output. By contrast, with the current UFLS scheme, in scenario 16 (49), the first/second load step is shed 500/700 (400/1,000) ms after the trip of the generator. This is a critical issue in small isolated power systems or parts of larger systems that become isolated due to the trip of a bus coupler or a bus tie breaker, as emphasized by some manufacturers of smart power management devices (ABB Group, 2019). The importance of the speed with which the load shedding is performed is confirmed by the measurement of the system's RoCoF. As can be seen in Table 3, in scenario 16/49, the RoCoF reaches the value of 3.93/1.94 Hz/s. These values, unthinkable in large interconnected power systems, indicate the need for a fastacting UFLS scheme in the system under study.

Out-Of-Sample Test Cases
The simulation of the trip of the generator with the highest power output has been simulated considering the other 2,916 hourly generation schedules provided by the above-mentioned unit commitment model. shed are shed 200 ms after the trip of the generator with the highest power output, as well as with the current UFLS scheme.
As can be seen in Table 4, the proposed shedding law, with the three implementation schemes, helps reduce considerably the average amount of load shed after the trip of the generator with the highest power output in comparison to the current UFLS scheme. Table 4 also shows the number of scenarios where the requirements specified in (Sigrist, 2015) are met with each shedding scheme.
As can be seen in Table 4, the proposed shedding law, with the implementation schemes IS2 and IS3, guarantees the fulfillment of the requirements specified in (Sigrist, 2015) in all out-ofsample scenarios.
In order to provide the reader with a more detailed view of the load shed in the 2,916 out-of-sample scenarios, Figure 10 shows the probability of exceedance (in the ordinate axis) of different values of load shed (in the abscissa axis) when P ex shed , P env shed , P 4f shed and P 8f shed are shed 200 ms after the trip of the generator with the highest power output, as well as with the currently used UFLS scheme. Finally, Figure 11 shows the cumulative distribution of  the difference between the load shed with the current UFLS scheme and the load shed with the proposed shedding law and the implementation scheme IS2. As can be seen in Figure 11, with the proposed shedding law and the implementation scheme IS2, the load shed is higher than the one shed with the current UFLS scheme in only 166 out of the 2,916 out-of-sample scenarios, by a value ranging from 402 to 867 kW. As mentioned above, the proposed shedding law, with the implementation scheme IS2, guarantees the fulfillment of the requirements specified in (Sigrist, 2015) in all out-of-sample scenarios and is fully compatible with the current MV distribution network of El Hierro. It would only require a modest investment in power management devices similar to the one commercialized by ABB (ABB Group, 2019). With the proposed shedding law and the implementation scheme IS3, the load shed would be lower than the one shed with the current UFLS scheme, and, as mentioned above, the frequency would meet the requirements specified in (Sigrist, 2015) in all out-of-sample scenarios. However, the implementation scheme IS3 would require a higher investment than IS2.

CONCLUSION
This study presents a novel approach to determine the amount of load to be shed in El Hierro power system in generation tripping events. The proposed approach consists in determining a linear regression model from a set of offline dynamic simulations of representative generation tripping events. The regression model obtained in this study provides the amount of load to be shed as a function of four variables, namely: the power imbalance, the system's inertia right after the tripping event, the power consumption of the fixed speed pumping units, and a binary variable which indicates whether one or more Pelton units are operating.
Three different implementation schemes of the proposed shedding law have been analyzed in this study. The proposed shedding law, with the three implementation schemes, has been compared with the UFLS scheme currently utilized in the El Hierro power system both in a set of 60 representative in-sample test cases and in a large set of 2,196 out-of-sample test cases. The results presented in this study show that the proposed shedding law, with the three implementation schemes, helps reduce the average amount of load shed across the 2,196 out-of-sample test cases. The proposed shedding law, with two of the studied implementation schemes, allows meeting the requirements specified in (Sigrist, 2015) for the system's frequency in all out-of-sample test cases. One of these implementation schemes is fully compatible with the current MV distribution network of El Hierro and would only require a modest investment in smart power management devices currently available in the industry.
The proposed approach (summarized in Figure 4) can be easily replicated in other power systems, both isolated and interconnected. The use of a unit commitment model is not necessary if historical series of hourly generation schedules are available. In interconnected power systems, it would be necessary to consider the power flows through the interconnection in the simulation model.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
JS has contributed to the conceptualization, investigation, methodology, and writing the original draft. JP has contributed to conceptualization, supervision, funding acquisition, and review and editing of the original draft. GM-L has contributed to the methodology and writing the original draft. Finally, DF-M has contributed to the methodology and review of the original draft.

FUNDING
The work presented in this study has been partially funded by the Spanish Ministerio de Economía y Competitividad under the project "Value of pumped-hydro energy storage in isolated power systems with high wind power penetration" of the National Plan for Scientific and Technical Research and Innovation 2013(Ref. ENE 2016.