Dynamic Modeling and Simulation of Non-Interconnected Systems under High-RES Penetration: The Madeira Island Case

: The defossilization of power generation is a prerequisite goal in order to reduce greenhouse gas emissions and transit for a sustainable economy. Achieving this goal requires increasing the penetration of renewable energy sources (RESs) such as solar and wind power. The gradual shrinking of conventional generation units in an energy map introduces new challenges to the stability of power systems as there is a considerable reduction of stored rotational energy in the synchronous generators (SGs) and the capability to control their power output, which has been taken for granted until today. Inertia and primary reserve reduction have a substantial e ﬀ ect on the ability of the power system to maintain its security and self-resilience during contingency events. Such issues become more evident in the case of non-interconnected islands (NII) as they have unique features associated with their small size and low inertia. The present study examines in depth the NII system of Madeira, which is composed of thermal, hydro, solid-waste, wind and solar generation units, and additional RES integration is planned for the near future. Electromagnetic transient (EMT) simulations are performed for both the current and future states of the system, including the installation of planned variable RES capacities. To alleviate the stability issues that occurred in the high-RES scenario, the introduction of a utility-scale battery energy storage system (BESS), capable of mitigating the active power imbalance due to the power system’s disturbances resultant of RES penetration, is examined. In addition, a comparison between a ﬂywheel energy storage system (FESS) and BESS is shortly investigated. The grid has been modeled and simulated utilizing the open-source, object-oriented modeling language Modelica. The dynamic simulation results proved that battery storage is a promising technology that can be a solution for transitioning to a sustainable power system, maintaining its self-resilience under severe disturbances such as rapid load changes, the tripping of generation units and short-circuits.


Introduction
In 2014, the EU raised the targets of reducing greenhouse gas (GHG) emissions by at least 40% and below the levels of 1990 by 2030, increasing the renewable energy share to 32% and improving energy efficiency by 32.5% [1]. Accomplishing such a transformation is strongly related to the capability of local power systems to incorporate renewable energy technologies while keeping sufficient levels of security

•
The transient modeling of a complex NII power system with a wide range of renewable and conventional power generation units, utilizing exclusively open-source software. This feature enhances the innovative character of the current work, as the value of open-source software for power systems' simulations has already been emphasized by the research community [37].
There are no publications of such complex real systems without using a commercial package, making this paper innovative. • A novel methodological approach for the impact assessment of high-RES penetration on the island's frequency stability. This approach involves examining both the current and future states in terms of stability and considering all the additional RES capacities with the simultaneous reduction of fossil-based plant capacity. Furthermore, power flow analysis is incorporated as an ancillary step in the scenario determination procedure to identify extreme cases in terms of the online inertia and primary frequency reserves. The final part of the methodological approach regards the conduction of transient simulations for disturbances such as rapid load changes and the tripping of generation units and short-circuits. The comparative analysis of the results reveals the impact of increased vRES capacities, decreased inertia and limited primary frequency reserves on the system's stability.

•
The consideration of energy storage technologies (BESS, FESS) to address the revealed instabilities in future low-inertia scenarios.
The following approach revealed the presence of stability issues in the future high-RES scenario. These issues are tackled by considering a BESS that contributes to power system inertia and primary frequency ancillary services, ensuring frequency stability. Simulations demonstrate the crucial role that such systems will play, along with the continuous increase of vRES deployment in order to support the frequency in contingency events. Additionally, a short comparison between BESS and a flywheel energy storage system (FESS) is investigated in terms of response time.
This paper is structured as follows: In Section 2, a description of the methodology used, scenarios considered and proposed BESS is presented. Section 3 presents a detailed description of Madeira's power system as it is currently and in the near future after the scheduled installation of new RES plants. Finally, Section 4 includes the dynamic simulation results accompanied by the relevant discussion.

Frequency Control Ancillary Services
Frequency stability concerns the ability of the system to maintain a steady frequency after a severe system upset, leading to a considerable imbalance between power production and demand. Traditionally, the subject of frequency stability was closely related to the operation of SGs, as they were predominant generation units. SGs present particular features that have greatly affected power system planning and operation since the introduction of electric machinery. From a physics standpoint, the turbine system and rotating components of each machine present mechanical inertia, and therefore, they are able to store kinetic energy in this rotating mass. Because energy can be withdrawn by or directed to these rotating masses during system disturbances, a system of machines is capable of adequately responding to fluctuations in net load and generation [38]. Furthermore, the active power of a SG can be controlled through the shaft torque. To ensure the resilience of a power system, its frequency is tightly adjusted by a combination of fast-acting closed-loop controllers at each machine (primary frequency response through governor controls), which are constantly controlling the machine's mechanical power [39]. Figure 1 presents a typical frequency response of generation loss. Three main areas can be distinguished, namely inertia response, primary frequency control and secondary frequency control. Each area corresponds to different transient phenomena which require the provision of corresponding ancillary services.
Energies 2020, 13, x FOR PEER REVIEW 4 of 26 and therefore, they are able to store kinetic energy in this rotating mass. Because energy can be withdrawn by or directed to these rotating masses during system disturbances, a system of machines is capable of adequately responding to fluctuations in net load and generation [38]. Furthermore, the active power of a SG can be controlled through the shaft torque. To ensure the resilience of a power system, its frequency is tightly adjusted by a combination of fast-acting closed-loop controllers at each machine (primary frequency response through governor controls), which are constantly controlling the machine's mechanical power [39]. Figure 1 presents a typical frequency response of generation loss. Three main areas can be distinguished, namely inertia response, primary frequency control and secondary frequency control. Each area corresponds to different transient phenomena which require the provision of corresponding ancillary services. Based on this categorization, conventional power plants (e.g., hydro, gas, coal) can offer the following frequency control ancillary services: Inertia response [40]: Because of the synchronous coupling of the machines with the grid, their rotational speed (i.e., ωm) is aligned with the angular velocity of the electromagnetic field (i.e., ωe). When a disturbance occurs leading to an imbalance between the two opposing torques, their sum on the rotor is nonzero, resulting in acceleration or deceleration according to the electromechanical swing equation (Equation (1)): where J represents the combined moment of inertia of the generator and the turbine (kg·m 2 ), Tm and Te represent mechanical and electrical torque, respectively (N·m) and Ta represents acceleration/deceleration torque (N·m). From the above equation, it can be deduced that mitigating the effect of power imbalances in terms of the rate of change of frequency (RoCoF) can be realized by enhancing the system's rotational inertia, utilizing fast reserves with RoCoF-based control. RoCoF is a key variable in the inertia response, defined as the time derivative of the power system frequency (df/dt). This variable was considered of low importance for systems with a generation mix dominated by synchronous generators due to the high inertia presence, which inherently compensates load imbalances and checks RoCoF in these cases. Recently, RoCoF has become important during considerable load-generation imbalances (attributed to the disconnection of either large loads or generators or system splits), as higher RoCoF values may be detected. This is due to low system inertia in cases of high instantaneous penetration of inverter-based units that are not synchronously connected to generation facilities. High df/dt values may threaten reliable system operation because of mechanical limitations of individual synchronous machines (inherent capability), protection devices triggered by a particular RoCoF threshold value or timing issues related to load-shedding schemes [41]. Based on this categorization, conventional power plants (e.g., hydro, gas, coal) can offer the following frequency control ancillary services: Inertia response [40]: Because of the synchronous coupling of the machines with the grid, their rotational speed (i.e., ω m ) is aligned with the angular velocity of the electromagnetic field (i.e., ω e ). When a disturbance occurs leading to an imbalance between the two opposing torques, their sum on the rotor is nonzero, resulting in acceleration or deceleration according to the electromechanical swing equation (Equation (1)): where J represents the combined moment of inertia of the generator and the turbine (kg·m 2 ), T m and T e represent mechanical and electrical torque, respectively (N·m) and T a represents acceleration/deceleration torque (N·m). From the above equation, it can be deduced that mitigating the effect of power imbalances in terms of the rate of change of frequency (RoCoF) can be realized by enhancing the system's rotational inertia, utilizing fast reserves with RoCoF-based control. RoCoF is a key variable in the inertia response, defined as the time derivative of the power system frequency (df /dt). This variable was considered of low importance for systems with a generation mix dominated by synchronous generators due to the high inertia presence, which inherently compensates load imbalances and checks RoCoF in these cases. Recently, RoCoF has become important during considerable load-generation imbalances (attributed to the disconnection of either large loads or generators or system splits), as higher RoCoF values may be detected. This is due to low system inertia in cases of high instantaneous penetration of inverter-based units that are not synchronously connected to generation facilities. High df /dt values may threaten reliable system operation because of mechanical limitations of individual synchronous machines (inherent capability), protection devices triggered by a particular RoCoF threshold value or timing issues related to load-shedding schemes [41].
Inverter-based generation lacking any control does not exhibit such inherent characteristics; consequently, high inverter penetration could result in higher df /dt values in a power system. The relationship between inverter penetration and RoCoF is, however, not straightforward and countermeasures, mostly in the form of control algorithms, need to be implemented carefully.
Primary control-droop control [42]: The coupling of active power to frequency forms the basis of frequency control where active power is adjusted according to linear characteristics, based on the following control equation (Equation (2)): where f is the system frequency, f 0 the nominal frequency, r p the droop gain, P is the active power of the unit and P 0 to the initial active power. Secondary control: The secondary control is responsible for restoring the frequency value to nominal (i.e., 50 or 60 Hz). This operation is needed, as the new operating point gained by the primary control has a steady-state frequency error, which must be reset. This method of control is called secondary control and is achieved through an integrator.

BESS Primary Frequency Ancillary Services
As already explained, the replacement of SGs with inverter-based RESs results in low inertia, as inverters provide small or zero inertia to the system due to their nonsynchronous connection to the grid and primary reserves reduction. The replacement of conventional units by inverter-based RES units reduces the number of units contributing to frequency support. To compensate for these effects, battery storage systems can assist primary frequency regulation with the appropriate control implementation in the inverter units. The general aim of the implemented inverter control through which the BESS interfaces with the grid is to modulate battery behavior to act as a synchronous generator in the event of a frequency change.
Since BESSs are capable of absorbing and delivering real power to the grid, they can offer ancillary services by properly controlling the inverter interfaced to the grid. This paper focuses only on inertial response and droop control, which are power-based services. The main challenge is to maintain a steady-state frequency in an acceptable value after disturbances. Figure 2 presents the schematic diagram for the BESS control that can provide the required services. Besides the setpoint for the active power (P ref ), the control also receives a signal for the system frequency (f ), which is compared with the nominal value (f 0 ) to compute the error. The error signal is routed through proportional and derivative sides, implementing droop and synthetic inertia control, respectively. For parameters, it receives two gains: one for the droop control (R) and one for the inertia response (k) [43], as well as the limits for the rated power exchange (P max and P min ). Inverter-based generation lacking any control does not exhibit such inherent characteristics; consequently, high inverter penetration could result in higher df/dt values in a power system. The relationship between inverter penetration and RoCoF is, however, not straightforward and countermeasures, mostly in the form of control algorithms, need to be implemented carefully.
Primary control-droop control [42]: The coupling of active power to frequency forms the basis of frequency control where active power is adjusted according to linear characteristics, based on the following control equation (Equation (2)): where f is the system frequency, f0 the nominal frequency, rp the droop gain, P is the active power of the unit and P0 to the initial active power. Secondary control: The secondary control is responsible for restoring the frequency value to nominal (i.e., 50 or 60 Hz). This operation is needed, as the new operating point gained by the primary control has a steady-state frequency error, which must be reset. This method of control is called secondary control and is achieved through an integrator.

BESS Primary Frequency Ancillary Services
As already explained, the replacement of SGs with inverter-based RESs results in low inertia, as inverters provide small or zero inertia to the system due to their nonsynchronous connection to the grid and primary reserves reduction. The replacement of conventional units by inverter-based RES units reduces the number of units contributing to frequency support. To compensate for these effects, battery storage systems can assist primary frequency regulation with the appropriate control implementation in the inverter units. The general aim of the implemented inverter control through which the BESS interfaces with the grid is to modulate battery behavior to act as a synchronous generator in the event of a frequency change.
Since BESSs are capable of absorbing and delivering real power to the grid, they can offer ancillary services by properly controlling the inverter interfaced to the grid. This paper focuses only on inertial response and droop control, which are power-based services. The main challenge is to maintain a steady-state frequency in an acceptable value after disturbances. Figure 2 presents the schematic diagram for the BESS control that can provide the required services. Besides the setpoint for the active power (Pref), the control also receives a signal for the system frequency (f), which is compared with the nominal value (f0) to compute the error. The error signal is routed through proportional and derivative sides, implementing droop and synthetic inertia control, respectively. For parameters, it receives two gains: one for the droop control (R) and one for the inertia response (k) [43], as well as the limits for the rated power exchange (Pmax and Pmin).

Methodology
A custom methodology was utilized to evaluate the impact of increasing RES penetration in the Madeira island system. The approach is schematically depicted in Figure 3 and includes the following steps.

Methodology
A custom methodology was utilized to evaluate the impact of increasing RES penetration in the Madeira island system. The approach is schematically depicted in Figure 3 and includes the following steps. primary reserve. The rationale behind this calculation is to identify the time at which maximum RES penetration (low demand and high RES production) is observed. When the grid's stability is maintained in the worst case, the robustness of the grid in the examined scenario is verified. The next step includes transient grid modeling in Modelica and simulation conduction for each case and the three disturbances examined, namely "load sudden change (2-8 MWe with a step of 2 MWe)," "short-circuit (SC)" of a transmission line 60 kV with the requirement of a 100 ms clearance time and "loss of the island's second-largest production unit online." After the evaluation and comparison of the results, a third scenario is considered, which corresponds to the future scenario with the additional integration of a BESS to mitigate the frequency stability issues that are revealed in the high-RES case. Lastly, the results are thoroughly assessed and taken into account in the decision-making process in the context of the long-term planning of the island's power system. At this point, it is worth noting that the simulated model does not contain any underfrequency shedding relays, which may activate during severe transients. However, a relay is included on every SG, which activates when the generator speed is higher than 1.03 pu or lower than 0.92 pu (overfrequency/underfrequency relay).
A critical security constraint for the Madeira system is the need to avoid automatic underfrequency load-shedding situations [45]. Specifically, it is assumed that load-shedding occurs if the RoCoF is lower than −1.5 Hz/s and the absolute frequency value is lower than 49 Hz [46].

System Description and Scenarios Definition
The present section aims to introduce the investigated power system to the reader by gathering its main operational components for production and consumption, as well as some passive elements, while avoiding a lengthy detailed description of each asset, for which the reader can refer to [47].
Madeira's power system has a yearly load peak of approximately 134 MWe. Its generation fleet is composed of two thermal plants, nine hydro plants, one pump-hydro plant, one solid-waste plant and a number of wind turbines (WTs) and photovoltaic (PV) parks, as well as some distributed PV generation rooftop-mounted systems [48]. Table 1 summarizes the installed power for all generation units per type, along with the number of generating units for conventional power plants. Initially, two scenarios are defined based on the current state of the power system ("reference scenario") and the planned future case of the system ("future scenario"), the latter including the additional variable RES capacities. The two scenarios were modeled in Matpower [44] and simulations were carried out utilizing the optimal power flow (OPF) engine. Through power flow analysis over a one-year period, the resulting dispatches are analyzed. For each one, extreme cases (best and worst) are identified with respect to the online instantaneous rotational inertia and available primary reserve. The rationale behind this calculation is to identify the time at which maximum RES penetration (low demand and high RES production) is observed. When the grid's stability is maintained in the worst case, the robustness of the grid in the examined scenario is verified. The next step includes transient grid modeling in Modelica and simulation conduction for each case and the three disturbances examined, namely "load sudden change (2-8 MW e with a step of 2 MW e )," "short-circuit (SC)" of a transmission line 60 kV with the requirement of a 100 ms clearance time and "loss of the island's second-largest production unit online." After the evaluation and comparison of the results, a third scenario is considered, which corresponds to the future scenario with the additional integration of a BESS to mitigate the frequency stability issues that are revealed in the high-RES case. Lastly, the results are thoroughly assessed and taken into account in the decision-making process in the context of the long-term planning of the island's power system.
At this point, it is worth noting that the simulated model does not contain any underfrequency shedding relays, which may activate during severe transients. However, a relay is included on every SG, which activates when the generator speed is higher than 1.03 pu or lower than 0.92 pu (overfrequency/underfrequency relay).
A critical security constraint for the Madeira system is the need to avoid automatic underfrequency load-shedding situations [45]. Specifically, it is assumed that load-shedding occurs if the RoCoF is lower than −1.5 Hz/s and the absolute frequency value is lower than 49 Hz [46].

System Description and Scenarios Definition
The present section aims to introduce the investigated power system to the reader by gathering its main operational components for production and consumption, as well as some passive elements, while avoiding a lengthy detailed description of each asset, for which the reader can refer to [47].
Madeira's power system has a yearly load peak of approximately 134 MW e . Its generation fleet is composed of two thermal plants, nine hydro plants, one pump-hydro plant, one solid-waste plant and a number of wind turbines (WTs) and photovoltaic (PV) parks, as well as some distributed PV generation rooftop-mounted systems [48]. Table 1 summarizes the installed power for all generation units per type, along with the number of generating units for conventional power plants. Regarding the passive grid components that are incorporated in the developed model, the following aspects are considered: Forty-six transformers: 6: 60 kV/30 kV, 6: 60 kV/6.6 kV, 21: 30 kV/6.6 kV, 4: 6.6 kV/60 kV and 9: 6.6 kV/30 kV.
The future scenario has been elaborated in communication with the Madeira's system operator, Eletricidade da Madeira, S.A. (EEM), and incorporates specific actions towards the energy transition of the island. An existing power plant will be entirely renovated in this framework, and an additional pumped-hydro plant will be constructed. Furthermore, the installed vRES power will be increased. The existing wind farm capacities will be extended, new solar plants will be created and the energy production from solar PV microproducers will be increased. Additionally, a BESS is planned to be installed in the central substation (60 kV), mainly providing ancillary services (primary, secondary frequency control and synthetic inertia response). Dedicated studies regarding the sizing of the BESS in Madeira have been performed [34][35][36] and their results were adopted in the context of the present study, considering a BESS of 10 MWh e and 15 MW e . In this study, the future scenario is investigated considering the inclusion and nonuse of the BESS, aiming to examine the battery's impact on the grid. Although the battery's location has been indicated by EMM, this option was further examined with the conduction of several simulations with the BESS installed in different busses. The simulations proved that the location of the BESS was irrelevant to its performance during disturbances. Table 2 gathers the new total installed capacities per technology.

Case Definition
Following the definition of reference and future scenarios, this section elaborates on the identification of the examined cases. As mentioned in Section 3, yearly OPF simulations were conducted for both scenarios to reveal the best and worst cases in terms of online instantaneous rotational inertia and available primary reserve. The only difference regarding the optimal dispatch formulation between the two scenarios is that, in the reference case, there is a requirement of 30 MW e minimum thermal production to ensure the system's stability, which translates into a considerable RES curtailment mainly during the night hours. This restriction is not present in the future scenario, as the BESS is responsible for securing the power system operation. Figures 4 and 5 depict the resulting operation for the two scenarios and highlight the extreme points that form the basis for the transient simulations.   Tables 3 and 4 focus on the identified best and worst cases and present units of operation. The purpose of including this information is to clarify which units are active at that time and which units can offer ancillary services. All conventional generation units can offer inertia response to the power system, since they are all composed of SGs, which contain rotating masses. Accordingly, all conventional units except steam plants contribute to primary frequency control according to their droop curve (2% for thermal, 3% for hydro). An important point is that many hydro units (mainly those with smaller installed power) operate close to their rated power; therefore, although they can assist in primary frequency control, this capability is limited by their rated power.  Tables 3 and 4 focus on the identified best and worst cases and present units of operation. The purpose of including this information is to clarify which units are active at that time and which units can offer ancillary services. All conventional generation units can offer inertia response to the power system, since they are all composed of SGs, which contain rotating masses. Accordingly, all conventional units except steam plants contribute to primary frequency control according to their droop curve (2% for thermal, 3% for hydro). An important point is that many hydro units (mainly those with smaller installed power) operate close to their rated power; therefore, although they can assist in primary frequency control, this capability is limited by their rated power.   Table 5 presents the calculated available primary reserves and online inertia for each case under consideration. As shown, the reduction of primary frequency reserves and online inertia among worst cases (26.87 MW e to 11.33 MW e and 1.36 pu to 0.65 pu, respectively) is substantial and renders the system particularly sensitive under disturbances. It is mentioned that due to an increase in the installed hydropower plants' capacity, the power system base is 250.17 MW e and 289.95 MW e for the reference and future cases, respectively.

Dynamic Simulations Results and Discussion
This section presents the results of the simulations for the defined cases. Figure 6 depicts an overview of the developed system's representation in the Dymola environment. A detailed description of the utilized dynamic models for all generation units is given in Appendix A. Their electrical and mechanical parameters have been set with typical values from the literature [49] based on their type and capacity.

Disturbances
Three types of disturbances have been examined, namely (i) load step-change, (ii) generation loss and (iii) short-circuit at the transmission level. The purpose of this section is to illustrate the differences between the future and reference scenarios and demonstrate that the increase of vRES penetration affects network stability and requires solutions to avoid curtailment while maintaining power system security and reliability.

Load Step-Change
To examine the transient behavior of the system under stiff load changes, four simulations are presented, parameterized by the value of the load step increase, i.e., 2 MWe, 4 MWe, 6 MWe and 8 MWe. Figures 7 and 8 show the frequency response and RoCoF for the best and worst cases of the reference scenario, respectively. It can be observed that the system responded successfully to this disturbance and no load-shedding relay activation was required as the load-shedding criteria were not met. Of course, the higher the load increase, the greater the steady-state error of frequency and RoCoF fluctuations. Even if load-shedding actions are avoided, frequency and RoCoF could reach extremely low values of 48.93 Hz and 1.14 Hz/s, respectively. This new low steady-state frequency will negatively affect the power quality provided to all customers, especially industries, as the operation of grid-coupled induction motors present in industrial environments is highly affected by the power system's frequency.

Disturbances
Three types of disturbances have been examined, namely (i) load step-change, (ii) generation loss and (iii) short-circuit at the transmission level. The purpose of this section is to illustrate the differences between the future and reference scenarios and demonstrate that the increase of vRES penetration affects network stability and requires solutions to avoid curtailment while maintaining power system security and reliability.

Load Step-Change
To examine the transient behavior of the system under stiff load changes, four simulations are presented, parameterized by the value of the load step increase, i.e., 2 MW e , 4 MW e , 6 MW e and 8 MW e . Figures 7 and 8 show the frequency response and RoCoF for the best and worst cases of the reference scenario, respectively. It can be observed that the system responded successfully to this disturbance and no load-shedding relay activation was required as the load-shedding criteria were not met. Of course, the higher the load increase, the greater the steady-state error of frequency and RoCoF fluctuations. Even if load-shedding actions are avoided, frequency and RoCoF could reach extremely low values of 48.93 Hz and 1.14 Hz/s, respectively. This new low steady-state frequency will negatively affect the power quality provided to all customers, especially industries, as the operation of grid-coupled induction motors present in industrial environments is highly affected by the power system's frequency. Figures 9 and 10 present the corresponding results of frequency response and RoCoF for the future scenario. As shown in Figure 10, severe frequency fluctuations occur in both 6 and 8 MW e load step-changes. The frequency rapidly reaches 46 Hz where generators have to be disconnected and the system collapses. In addition, in the RoCoF diagrams, it can be observed that the frequency drops sharply due to the presence of reduced inertia. disturbance and no load-shedding relay activation was required as the load-shedding criteria were not met. Of course, the higher the load increase, the greater the steady-state error of frequency and RoCoF fluctuations. Even if load-shedding actions are avoided, frequency and RoCoF could reach extremely low values of 48.93 Hz and 1.14 Hz/s, respectively. This new low steady-state frequency will negatively affect the power quality provided to all customers, especially industries, as the operation of grid-coupled induction motors present in industrial environments is highly affected by the power system's frequency.  Figures 9 and 10 present the corresponding results of frequency response and RoCoF for the future scenario. As shown in Figure 10, severe frequency fluctuations occur in both 6 and 8 MWe load step-changes. The frequency rapidly reaches 46 Hz where generators have to be disconnected and the system collapses. In addition, in the RoCoF diagrams, it can be observed that the frequency drops sharply due to the presence of reduced inertia.   Figures 9 and 10 present the corresponding results of frequency response and RoCoF for the future scenario. As shown in Figure 10, severe frequency fluctuations occur in both 6 and 8 MWe load step-changes. The frequency rapidly reaches 46 Hz where generators have to be disconnected and the system collapses. In addition, in the RoCoF diagrams, it can be observed that the frequency drops sharply due to the presence of reduced inertia.  Lastly, Figures 11 and 12 depict the results of the frequency response for the same disturbances after the BESS integration. As shown in these diagrams, no shedding relay needs to be activated after battery installation, since the load-shedding criteria are not met. In all load step increments, the system maintains stability and reduces both RoCoF deviations and steady-state frequency error as the inverter's time response is fast enough to increase the battery's active power output and emulate inertia.    Lastly, Figures 11 and 12 depict the results of the frequency response for the same disturbances after the BESS integration. As shown in these diagrams, no shedding relay needs to be activated after battery installation, since the load-shedding criteria are not met. In all load step increments, the system maintains stability and reduces both RoCoF deviations and steady-state frequency error as the inverter's time response is fast enough to increase the battery's active power output and emulate inertia.   Figure 13 presents the active power output of BESS during the disturbances. In all cases, BESS operates within the permitted limits (15 MWe). At the moment when the disturbance event occurs, the power of the BESS increases rapidly due to virtual inertia control and, due to the droop control, follows a pattern inversely proportional to the frequency.  Lastly, Figures 11 and 12 depict the results of the frequency response for the same disturbances after the BESS integration. As shown in these diagrams, no shedding relay needs to be activated after battery installation, since the load-shedding criteria are not met. In all load step increments, the system maintains stability and reduces both RoCoF deviations and steady-state frequency error as the inverter's time response is fast enough to increase the battery's active power output and emulate inertia.   Figure 13 presents the active power output of BESS during the disturbances. In all cases, BESS operates within the permitted limits (15 MWe). At the moment when the disturbance event occurs, the power of the BESS increases rapidly due to virtual inertia control and, due to the droop control, follows a pattern inversely proportional to the frequency.  Figure 13 presents the active power output of BESS during the disturbances. In all cases, BESS operates within the permitted limits (15 MW e ). At the moment when the disturbance event occurs, the power of the BESS increases rapidly due to virtual inertia control and, due to the droop control, follows a pattern inversely proportional to the frequency. Tables 6 and 7 summarize the most characteristic parameters of the examined scenarios described above. Inertia reduction is reflected not only to the RoCoF nadirs of the reference and future scenarios but also by the fact that, since the frequency drops sharply, the generator's prime movers fail to increase their power rates and prevent system collapse. In addition, it is observed that the battery installation for primary frequency control and inertia support can lead to a smaller RoCoF and steady-state frequency error than the reference case. For the most severe load step-change of 8 MWe, in the future with BESS scenario, the resulting RoCoF oscillates in the interval [−0.18, −0.015] Hz/s. The corresponding interval for the reference scenario is constrained between [−0.3, −0.03] Hz/s, demonstrating a strong advancement over the current operational mode, also considering the reduction of the physical system inertia from 1.36 pu to 0.65 pu.  Figure 13. Contribution of the BESS to the system response of a load step-change for the best and worst cases.
Tables 6 and 7 summarize the most characteristic parameters of the examined scenarios described above. Inertia reduction is reflected not only to the RoCoF nadirs of the reference and future scenarios but also by the fact that, since the frequency drops sharply, the generator's prime movers fail to increase their power rates and prevent system collapse. In addition, it is observed that the battery installation for primary frequency control and inertia support can lead to a smaller RoCoF and steady-state frequency error than the reference case. For the most severe load step-change of 8 MW e , in the future with BESS scenario, the resulting RoCoF oscillates in the interval [−0.18, −0.015] Hz/s. The corresponding interval for the reference scenario is constrained between [−0.3, −0.03] Hz/s, demonstrating a strong advancement over the current operational mode, also considering the reduction of the physical system inertia from 1.36 pu to 0.65 pu.

. Generation Loss
This scenario investigation is of great importance in the future case because facing great power fluctuations from renewable energy sources is very common. For this purpose, it is important to ensure that the power system will remain balanced under high-RES penetration. In order to test power system Energies 2020, 13, 5786 14 of 25 stability under extreme conditions, the tripping of the second-largest unit is investigated. According to the reference case, this unit is a wind power plant with an installed capacity of 8.32 MW e . According to the future case, the second-largest unit is a solar power plant with a 14.22 MW e maximum power output. Although the system maintained its stability in the best case of the future scenario as depicted in Figure 14, it collapsed in the corresponding worst case of Figure 15. Furthermore, load-shedding relays activated before system collapse as the load-shedding criteria were met. Additionally, BESS installation leads to secure system operation without load-shedding relay activation, even after the complete loss of a significant production unit.

. Generation Loss
This scenario investigation is of great importance in the future case because facing great power fluctuations from renewable energy sources is very common. For this purpose, it is important to ensure that the power system will remain balanced under high-RES penetration. In order to test power system stability under extreme conditions, the tripping of the second-largest unit is investigated. According to the reference case, this unit is a wind power plant with an installed capacity of 8.32 MWe. According to the future case, the second-largest unit is a solar power plant with a 14.22 MWe maximum power output. Although the system maintained its stability in the best case of the future scenario as depicted in Figure 14, it collapsed in the corresponding worst case of Figure  15. Furthermore, load-shedding relays activated before system collapse as the load-shedding criteria were met. Additionally, BESS installation leads to secure system operation without load-shedding relay activation, even after the complete loss of a significant production unit.   Table 8 presents the characteristic values derived from the diagrams described above. In the worst case of the future scenario, the RoCoF's nadir is equal to −3.69 Hz/s before the system collapse. On the contrary, in the worst case of the future with the integration of BESS scenario, the RoCoF and frequency remained below −0.98 Hz/s and 48.87 Hz, respectively, far from the load-shedding activation criteria (−1.5 Hz/s, 49 Hz), revealing a significant response improvement. The proposed control method for the BESS inverter leads not only to reliable system operation but also a severe reduction of the minimum RoCoF and prevention of load-shedding relay activation.   Table 8 presents the characteristic values derived from the diagrams described above. In the worst case of the future scenario, the RoCoF's nadir is equal to −3.69 Hz/s before the system collapse. On the contrary, in the worst case of the future with the integration of BESS scenario, the RoCoF and frequency remained below −0.98 Hz/s and 48.87 Hz, respectively, far from the load-shedding activation criteria (−1.5 Hz/s, 49 Hz), revealing a significant response improvement. The proposed control method for the BESS inverter leads not only to reliable system operation but also a severe reduction of the minimum RoCoF and prevention of load-shedding relay activation. This section examines the impact of a three-phase short-circuit fault on the grid, which is considered one of the major disturbances that can occur in a power system and explores the system's ability to clear the fault and maintain stability. On mainly overhead line systems, such as Madeira's power system, the majority of short-circuit faults (approximately 80-90%) tend to appear on overhead lines and the rest on substation equipment and busbars combined [50]. This is due to the exposure of the transmission and distribution assets to the weather phenomena, making them vulnerable to wind gusts and thunders and, at the same time, causing earlier degradation to the insulating materials. Additionally, the increase of vRES penetration decreases the short-circuit currents due to the inverters' current limiters, introducing a further challenge to the grid as it becomes harder to identify and clear the fault quickly. It is evident that the thorough investigation of this kind of fault is crucial for the resilience of the system under consideration.
In this framework, short-circuit scenarios have been examined for the most crucial point of the system, i.e., the 60 kV overhead transition lines. According to the scenario under consideration, the fault occurred on a 60 kV transition line, which supplies the west side of the island, and was cleared after 100 ms (five cycles) according to the Madeira's power system fault-clearing standards.
As can be observed in Figures 16 and 17, the system responded successfully to this disturbance and maintained its stability in all cases, although in the best case of the reference scenario, the fault led to the desynchronization of a hydro unit. As Figure 18 depicts, the short-circuit caused very high oscillations to the hydro unit's speed and, as a result, to power system frequency, which are depreciated after the successful disconnection of the unit from the power system. Generally, the system has no problem maintaining its stability and is driven to a new frequency almost identical to the rated frequency (50 Hz) without the presence of secondary frequency control. In the future case, frequency deviation is more severe, which is mainly due to the significant reduction of system inertia.
Energies 2020, 13, x FOR PEER REVIEW 16 of 26 In this framework, short-circuit scenarios have been examined for the most crucial point of the system, i.e., the 60 kV overhead transition lines. According to the scenario under consideration, the fault occurred on a 60 kV transition line, which supplies the west side of the island, and was cleared after 100 ms (five cycles) according to the Madeira's power system fault-clearing standards.
As can be observed in Figures 16 and 17, the system responded successfully to this disturbance and maintained its stability in all cases, although in the best case of the reference scenario, the fault led to the desynchronization of a hydro unit. As Figure 18 depicts, the short-circuit caused very high oscillations to the hydro unit's speed and, as a result, to power system frequency, which are depreciated after the successful disconnection of the unit from the power system. Generally, the system has no problem maintaining its stability and is driven to a new frequency almost identical to the rated frequency (50 Hz) without the presence of secondary frequency control. In the future case, frequency deviation is more severe, which is mainly due to the significant reduction of system inertia.     Based on these figures, along with Table 9, it is evident that BESS results in both frequency deviations and the reduction of RoCoF. This contribution is primarily due to the derivative control presented in Section 2.2, which helped to increase the system inertia virtually along with BESS droop control since, after fault clearance, it contributed more than the derivative control. Furthermore, it is observed that in the worst case of the future scenario, the RoCoF oscillates between [−6.69, −4.32] Hz/s and [−4.27, −3.22] Hz/s in the same case of the future with BESS scenario, which is a huge enhancement, also taking into account the reduction of physical system inertia from 1.36 pu to 0.65 pu.

BESS/FESS Comparison
This paragraph examines the consideration of a flywheel storage system as an alternative solution to the stability problem in Madeira. The FESS model is presented in Appendix A and includes the rotor, a permanent magnet generator and two inverters in a back-to-back topology scheme. The battery in the overall system model has been replaced by the corresponding FESS, and simulations have been conducted for the case of the future scenario. The main parameters include Based on these figures, along with Table 9, it is evident that BESS results in both frequency deviations and the reduction of RoCoF. This contribution is primarily due to the derivative control presented in Section 2.2, which helped to increase the system inertia virtually along with BESS droop control since, after fault clearance, it contributed more than the derivative control. Furthermore, it is observed that in the worst case of the future scenario, the RoCoF oscillates between [−6.69, −4.32] Hz/s and [−4.27, −3.22] Hz/s in the same case of the future with BESS scenario, which is a huge enhancement, also taking into account the reduction of physical system inertia from 1.36 pu to 0.65 pu. Table 9. Results of three-phase short-circuit (3ph) short-circuit for all scenarios and cases.

BESS/FESS Comparison
This paragraph examines the consideration of a flywheel storage system as an alternative solution to the stability problem in Madeira. The FESS model is presented in Appendix A and includes the rotor, a permanent magnet generator and two inverters in a back-to-back topology scheme. The battery in the overall system model has been replaced by the corresponding FESS, and simulations have been conducted for the case of the future scenario. The main parameters include the generator's rated power of 18 MVA and the FESS inertia of 3·10 7 kg·m 2 , which were selected to be analogous to the BESS characteristics. Figure 19 presents the comparison between the battery and flywheel performance in an 8 MWe load step increase. It is observed that the differences between the two storage solutions are negligible. This is expected as the response time of the storage systems for the most part is determined by the inverter selection for the application and overall system design [51]. Therefore, the appropriate storage technology must be selected based on economic criteria regarding the investment and maintenance costs as their technical performance regarding the frequency of ancillary services is identical.

Conclusions
The transition from a fossil-fuel-based power system to a sustainable system with very high instantaneous penetrations of vRES requires the development of efficient methods for matching supply and demand over multiple timescales. Since traditional power system operation and security are based on synchronous generator characteristics, it is of utmost importance to ensure that vRES power systems operate in such a way that does not endanger grid stability. This includes designing an inverter-based system to provide system stability and additional grid services necessary for proper system operations. This paper focused on the primary frequency control of the Madeira island power system under high vRES penetration with the integration of a BESS (15 MWe, 10 MWhe) to provide frequency support under severe disturbances. In island systems, the location of a BESS is of minor importance, and in the context of the present study, the battery is considered to be installed in one of the major nodes of the system.
Through scenario-based fault simulation, the system's response has been examined. It is demonstrated that with the proposed BESS, fault scenarios that would otherwise lead to loadshedding procedures can be mitigated with success. More precisely, three contingency events have been examined, i.e., load step-change, generation loss and short-circuit in a transmission line. It is also important to remark on some of the difficulties that occurred during the conduction of the present study. The most challenging part regards the high level of expertise required to develop

Conclusions
The transition from a fossil-fuel-based power system to a sustainable system with very high instantaneous penetrations of vRES requires the development of efficient methods for matching supply and demand over multiple timescales. Since traditional power system operation and security are based on synchronous generator characteristics, it is of utmost importance to ensure that vRES power systems operate in such a way that does not endanger grid stability. This includes designing an inverter-based system to provide system stability and additional grid services necessary for proper system operations. This paper focused on the primary frequency control of the Madeira island power system under high vRES penetration with the integration of a BESS (15 MW e , 10 MWh e ) to provide frequency support under severe disturbances. In island systems, the location of a BESS is of minor importance, and in the context of the present study, the battery is considered to be installed in one of the major nodes of the system.
Through scenario-based fault simulation, the system's response has been examined. It is demonstrated that with the proposed BESS, fault scenarios that would otherwise lead to load-shedding procedures can be mitigated with success. More precisely, three contingency events have been examined, i.e., load step-change, generation loss and short-circuit in a transmission line. It is also important to remark on some of the difficulties that occurred during the conduction of the present study. The most challenging part regards the high level of expertise required to develop the necessary modules and successfully perform the simulations using open-source software tools. There are many procedures that must be performed manually in contrast to the automated methods of corresponding proprietary tools. Modelica, in particular, requires deep user knowledge to develop and debug models in the power system domain. Additionally, there were conflicts about the compatibility of RMS simulations with high-RES scenarios and if the use of the phasor domain simulation, which is the most common for these studies, can predict all instabilities. Another important issue that is always present in these studies is the gathering and quality of the input data, which is often a sensitive subject, especially in the context of power consumption.
The simulation results proved that using a BESS for this purpose in an islanded non-interconnected island can prevent possible curtailment or load-shedding actions and is a promising solution for even higher RES penetration, ensuring grid security and stability. Similar conclusions were extracted from the analysis of using FESS for the frequency control of the island's system, revealing that both batteries and flywheels are two energy storage technologies that guarantee the grid's stability in a high-RES NII. Acknowledgments: All parameters of the considered power system were kindly provided by EEM. The examined scenarios of the reference and future case were discussed and agreed with EEM. The authors thank Diogo Vasconcelos, Aires Henriques and Henrique Pinto Correia for their support and review.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Dynamic Models Description
The implementation of the power system's dynamic model was carried out in Modelica/Dymola with the use of the open-source PowerSystems library [32] as the basis. Basic models such as grid lines, transformers, loads and synchronous generators were based on the PowerSystems library and others, such as synchronous generator controls, BESS, BESS controls, PV plant, WT plant and inverters, were developed from scratch. In this Appendix, the developed models are presented along with a short description regarding the most important implementation details.

Synchronous Generator and Controls
The SG of the PowerSystems library is a complete electromechanical model. As shown in Figure A1, it receives a mechanical torque and an excitation voltage as inputs and outputs voltage and current.

Automatic Voltage Regulator (AVR)
The purpose of an AVR [52] is to control the voltage at the generator's terminal and keep it at a constant, desired value. AVR IEEE T1 [53], whose diagram is shown in Figure A2, was selected. Power System Stabilizer (PSS) In order to increase the damping of rotor oscillation caused by a disturbance [54], an IEEE type PSS1A power system stabilizer [55] was developed, based on Figure A3. Turbine Governors (TG) To properly simulate the mechanical reactions of synchronous generators, it was necessary to use turbine governor models for all power plant types [56]. Figures A4-A6 show the schematic The purpose of an AVR [52] is to control the voltage at the generator's terminal and keep it at a constant, desired value. AVR IEEE T1 [53], whose diagram is shown in Figure A2, was selected.

Automatic Voltage Regulator (AVR)
The purpose of an AVR [52] is to control the voltage at the generator's terminal and keep it at a constant, desired value. AVR IEEE T1 [53], whose diagram is shown in Figure A2, was selected. Power System Stabilizer (PSS) In order to increase the damping of rotor oscillation caused by a disturbance [54], an IEEE type PSS1A power system stabilizer [55] was developed, based on Figure A3. Turbine Governors (TG) To properly simulate the mechanical reactions of synchronous generators, it was necessary to use turbine governor models for all power plant types [56]. Figures A4-A6 show the schematic diagrams of TG for all types of power plants. Power System Stabilizer (PSS) In order to increase the damping of rotor oscillation caused by a disturbance [54], an IEEE type PSS1A power system stabilizer [55] was developed, based on Figure A3.

Automatic Voltage Regulator (AVR)
The purpose of an AVR [52] is to control the voltage at the generator's terminal and keep it at a constant, desired value. AVR IEEE T1 [53], whose diagram is shown in Figure A2, was selected. Power System Stabilizer (PSS) In order to increase the damping of rotor oscillation caused by a disturbance [54], an IEEE type PSS1A power system stabilizer [55] was developed, based on Figure A3. Turbine Governors (TG) To properly simulate the mechanical reactions of synchronous generators, it was necessary to use turbine governor models for all power plant types [56]. Figures A4-A6 show the schematic diagrams of TG for all types of power plants. Turbine Governors (TG) To properly simulate the mechanical reactions of synchronous generators, it was necessary to use turbine governor models for all power plant types [56]. Figures A4-A6 show the schematic diagrams of TG for all types of power plants.  Figure A4. Hydro turbine governor type GovHydroIEEE0 (adopted by [57]).
Average Inverter Model RES and BESS are interfaced with the three-phase AC grid through an average inverter model presented in Figure A7 [59,60].  . Hydro turbine governor type GovHydroIEEE0 (adopted by [57]).
Average Inverter Model RES and BESS are interfaced with the three-phase AC grid through an average inverter model presented in Figure A7 [59,60].  . Steam turbine governor (adopted by [58]).
Average Inverter Model RES and BESS are interfaced with the three-phase AC grid through an average inverter model presented in Figure A7 [59,60].  . Diesel turbine governor (adopted by [57]).
Average Inverter Model RES and BESS are interfaced with the three-phase AC grid through an average inverter model presented in Figure A7 [59,60].  Figure A7. Modelica's average inverter model.

Photovoltaic Power Plant
As a result of the investigations and discussions, a key simplification that could be applied in a central station PV system model is that the dynamics related to the DC side of the inverter (PV array dynamics, inverter DC link and voltage regulator) shall be ignored [61]. Therefore, according to the above, the developed solar plant consists of a DC voltage source along with an average inverter to connect it to the three-phase grid, as presented in Figure A8. central station PV system model is that the dynamics related to the DC side of the inverter (PV array dynamics, inverter DC link and voltage regulator) shall be ignored [61]. Therefore, according to the above, the developed solar plant consists of a DC voltage source along with an average inverter to connect it to the three-phase grid, as presented in Figure A8.

Wind Power Plant
The double-fed induction generator (DFIG) developed wind turbine model is shown in Figure  A9

Battery Energy Storage System (BESS)
A lithium-ion battery module was developed in order to simulate the BESS operation and calculate important BESS parameters, i.e., state of charge (SOC), terminal voltage, etc. The battery was modeled based on its equivalent electrical circuit in [63], consisting of a voltage source in series with its internal impedance. The latter is represented as a two-series resistor and two RC networks. The developed module is presented in Figure A10.  Figure A9 [62].
connect it to the three-phase grid, as presented in Figure A8.

Wind Power Plant
The double-fed induction generator (DFIG) developed wind turbine model is shown in Figure  A9 [62].

Battery Energy Storage System (BESS)
A lithium-ion battery module was developed in order to simulate the BESS operation and calculate important BESS parameters, i.e., state of charge (SOC), terminal voltage, etc. The battery was modeled based on its equivalent electrical circuit in [63], consisting of a voltage source in series with its internal impedance. The latter is represented as a two-series resistor and two RC networks. The developed module is presented in Figure A10.

Battery Energy Storage System (BESS)
A lithium-ion battery module was developed in order to simulate the BESS operation and calculate important BESS parameters, i.e., state of charge (SOC), terminal voltage, etc. The battery was modeled based on its equivalent electrical circuit in [63], consisting of a voltage source in series with its internal impedance. The latter is represented as a two-series resistor and two RC networks. The developed module is presented in Figure A10.

Flywheel Energy Storage System (FESS)
The developed flywheel storage model is presented in Figure A11. The developed flywheel storage model is presented in Figure A11. The developed flywheel storage model is presented in Figure A11.