Optimal Scheduling of an Isolated Wind-Diesel-Energy Storage System Considering Fast Frequency Response and Forecast Error

Nowadays, the hybrid wind–diesel system is widely used on small islands. However, the operation of these systems faces a major challenge in frequency control due to their small inertia. Furthermore, it is also difficult to maintain the power balance when both wind power and load are uncertain. To solve these problems, energy storage systems (ESS) are usually installed. This paper demonstrates the effectiveness of using ESS to provide Fast Frequency Response (FFR) to ensure that the frequency criteria are met after the sudden loss of a generator. An optimal day-ahead scheduling problem is implemented to simultaneously minimize the operating cost of the system, take full advantage of the available wind power, and ensure that the ESS has enough energy to provide FFR when the wind power and demand are uncertain. The optimization problem is formulated in terms of two-stage chance-constrained programming, and solved using a Modified Sample Average Approximation (MSAA) algorithm—a combination of the traditional Sample Average Approximation (SAA) algorithm and the k-means approach. The proposed method is tested with a realistic islanded power system, and the effects of the ESS size and its response time is analyzed. Results indicate that the proposed model should perform well under real-world conditions.


Introduction
Frequency is an important criterion in the power system's operation and is related to the instantaneous balance between supply and demand. To ensure stable operation of a power system, the balance between power demand and supply must be kept at all times; the system frequency is only allowed to vary in a tight band around the nominal value. Large frequency disturbances, caused by events such as the sudden loss of a generator, lead to serious active power imbalances and may lead to load shedding or partial or complete blackout. Fortunately, immediately after a frequency disturbance, the kinetic energy stored in the spinning masses of the generators is released into the power system to preserve the power balance, thereby reducing the rate of frequency change. This process is called the Inertial Response (IR) of the generator. The most crucial frequency control activity is the Primary Frequency Control Response (PFR), which is based on the characteristic of the conventional speed governor. PFR automatically starts in the event of a large frequency deviation to adjust the power output of generators. Large synchronous generators can provide both IR and PFR, thereby ensuring the frequency stability of the power system.
In recent years, the penetration of renewable energy sources (RES) such as wind and solar into the power system is increasing rapidly. Unlike conventional power plants, which are based on synchronous machines, RES power plants use inverter-based generators that do not have IR. Therefore, the overall power system inertia is reduced by the increasing penetration of RES. Furthermore, due to the stochastic nature of wind and solar irradiation, the reserved capacity for PFR from RES is also uncertain. Both of these factors contribute to the more complicated frequency control problem for power systems containing large fractions of wind and solar generation [1]. In islanded power systems, frequency control and regulation are even more challenging because the primary resources are diesel generators (DGs) with low inertia and limited operating capability.
With the increasing penetration of renewable generation and energy storage in power systems, the Fast Frequency Response (FFR) method has been introduced as a measure to improve frequency stability. In the Australian electricity market, FFR is defined as "any type of rapid active power increase or decrease by generation or load, in a timeframe of less than two seconds, to correct supply-demand imbalances and assist with managing frequency" [2]. There have been several studies of FFR encompassing a wide range of technologies [3][4][5][6]. In one study [6], it was shown that wind generators (WG) can provide IR for a very short duration (~10 s). Although this method proved to be useful for frequency regulation, the kinetic energy provided by wind turbines is highly dependent on the wind speed; as a result, insufficient support is delivered in the case of low wind speed. Furthermore, in a low inertia system, the frequency can drop below the threshold of the Under-Frequency Load Shedding (UFLS) relay within 1 s, so the response of a WG is not effective. For the case of photovoltaics (PV), another method involves keeping the PV power setpoint below the total available power, at the expense of economic performance [2].
With a very short response time, which can be less than 250 ms depending on the technology, energy storage systems (ESS) are able to instantly increase or decrease their power output to counteract a system power imbalance. There have been many previous studies focusing on the support of an ESS in frequency response such as [7][8][9][10][11][12]. In [9][10][11][12], the authors focused on the size of the ESS, whereas control strategies for an ESS to provide virtual inertia are proposed in [7,8]. The results presented in these articles show the effectiveness of using an ESS for frequency response control.
An interesting research approach in frequency-constrained operation planning is to include frequency constraints in Unit Commitment (UC) models. Ahmadi and Ghasemi [13] proposed a UC formulation including a frequency limit constraint based on the general-order system frequency response model. In contrast, a first-order model for a governor-prime mover system was used in [14]. However, the UC models in [13,14] did not consider the uncertainty in the available wind power, so the results are less reliable for actual operation. Teng et al. [15] develops a sophisticated representation for the frequency dynamics, which includes load damping, but this work did not consider the application of ESS in frequency response. A more recent study [16] developed a UC framework in which the ESS is considered to provide frequency response. However, the wind power uncertainty in this study is described by only three scenarios: the central forecast, the upper bound, and the lower bound. Although this approach helps reduce the computational complexity, it may lead to conservative solutions with higher operating cost.
In this work, we propose a frequency stability-constrained UC model and apply it to a realistic isolated wind-diesel system on Phu Quy Island, Binh Thuan province, Vietnam. The UC model in this work focuses on FFR. The primary goal of ESS in this model is to compensate for the fluctuation of wind and solar generation and to help increase the energy produced from renewable sources (rather than using curtailment). Besides, this ESS provides FFR in large frequency disturbances, such as in the loss of a generator, as an ancillary service. The proposed model has practical implications for isolated power systems with a high penetration level of renewable resources.
To account for the stochastic nature of wind power and demand, the optimal scheduling in this work is formulated as a two-stage chance-constrained optimization problem. The constraints related to uncertain parameters are written as probabilistic constraints with a chosen risk level [17,18]. A common method used to solve chance-constrained problems is the Sample Average Approximation (SAA) algorithm, which involves Monte Carlo simulation to approximate the distribution function of a Energies 2019, 12, 843 3 of 20 random vector using N samples [19][20][21][22][23]. Although SAA is simple and convenient, all N samples are considered to have the same probability regardless of the true distribution of the random vector, so the number of samples must be large to ensure that a feasible solution is found. If there are many uncertain parameters, the size of the optimization problem increases, and a significantly longer computing time is required. For these reasons, it is necessary to improve the algorithm, which is one of the objectives of the present work.
The salient features of the present study include the following.

1.
The research focuses on frequency stability-constrained UC models. The ESS, which is employed to keep power balance and take advantage of wind power, is considered to provide fast frequency response (FFR) in large frequency disturbances, such as loss of a generator. The frequency dynamics is approximated using a first-order representation.

2.
The proposed UC model is based on a two-stage stochastic programming framework which is suitable for the short-term planning of power systems with uncertain sources. The model is formulated as a chance constraint problem to allow a certain risk level in the day-ahead scheduling.

3.
The impact of ESS sizing and its response time on frequency nadir is analyzed.
As the computing time required for solving chance-constrained optimization is usually significant, in this paper, a Modified Sample Average Approximation (MSAA) method is presented and applied to solve the proposed optimization problem. The combination of SAA and k-means clustering approach is proven to be more effective than the original SAA approach.
The rest of the paper is organized as follows. Section 2 demonstrates the application of the ESS for FFR. Section 3 presents the mathematical formulation of the proposed chance-constrained optimization problem to determine the optimal scheduling of the power system. Section 4 presents the MSAA algorithm. The computation results are collected and analyzed in Section 5. Section 6 concludes the paper.

Fast Frequency Response and the Role of the ESS
As discussed in Section 1, IR plays an important role in maintaining frequency after a generator is lost. This process slows down the change of frequency before the governors fully react to provide PFR. We can evaluate this process using two important criteria including the rate-of-change-of-frequency (RoCoF) and the lowest frequency known as the frequency nadir f nadir .
Consider, for example, a power system with I generators. If at time t generator j with power output P t j (kW) is lost, the RoCoF immediately after the contingency event is defined as where M H is the system inertia (kW·s/Hz) after the loss of generation j and M H is a function of the inertia of the online generators: where H i , P maxi , and w i are the inertia constant, maximum capacity, and ON/OFF state of the remaining generators, respectively; w i is a binary variable that is equal to 1 if generator i is online and 0 if it is offline. If the inertia of the system is sufficient, the frequency will stop before reaching the threshold of the UFLS relays. However, in an islanded power system, the primary resources are DGs with low inertia constants. Additionally, modern wind turbines are connected to the grid through an electronic power converter, which does not contribute to the system inertia. It is easy to see that the lower the system inertia, the faster the frequency drops. Therefore, even though DGs can adjust their power output quickly, it is still difficult to arrest the frequency decline before reaching the minimum allowable frequency. The concept of FFR is applied to solve this problem. Unlike both IR and PFR, which slowly adjust the power based on the frequency deviation thus arresting the deviation and restoring frequency, the objective of FFR is to immediately inject active power into the grid to correct the power imbalance. This is implemented based on wind turbines or ESS that can change their power output almost instantaneously. FFR can be considered as a measure to compensate for the "interval" between IR and PFR, during which the frequency is too low and PFR is still not fully active.
Note that FFR cannot completely replace PFR-it is only a support measure while waiting for the DGs to provide PFR. Thus, a sustained FFR duration is not necessary. Riesz et al. [24] show that in ERCOT (Texas) this duration is 10 min, while the EirGrid/SONI (Ireland) only requires an 8-second FFR.
An important requirement for FFR is fast response time, and systems with a higher RoCoF will require a faster response time. Assuming the system has a nominal frequency of 50 Hz and a minimum frequency of 49 Hz, FFR must fully react within 250 ms in the case of a 4 Hz/s RoCoF. The response time depends on not only the detection method but also the type of device used to provide FFR. The times required to detect contingency and to send the control signal, as well as the reaction time and rise time of the ESS, are summarized in Table 1 [3]. It can be seen that, with total response times ranging from 100 ms to 200 ms, the ESS is suitable for providing FFR. In this section, we outline the constraints on the power output of the DGs for each hour and the response of the ESS needed to satisfy the frequency criterion, f min .
Using the first-order model for a governor-prime mover previously presented in [14,25], we arrive at an approximate model of the system's response after a sudden generation loss of amount P t j and the application of ESS for FFR, as described in Figure 1. After the governor's dead time t d , which corresponds to the frequency dead band f db , the power output of the DGs will change due to the governor's response with the system ramp rate K = ∑ i =j k i P maxi / ∑ i =j P maxi . On the other hand, a control signal is also sent to the ESS to increase its output from P ch,t ss or P disch,t ss to P t dispost . The adjustment provided by the ESS is P t adjust = P t dispost − P disch,t SS − P ch,t SS ( Figure 1). To simplify the model, we make two assumptions:

•
The rise time of ESS is negligible.

•
The ESS can fully compensate for the power shortage. This means that the frequency decay stops immediately after the ESS fully responds, at which point the frequency nadir occurs.
The time evolution of the system frequency deviation can be described by where ∆P DGs (t) and ∆P ESS (t) describe the additional power provided by DGs (due to the response of the governors) and ESS, respectively. ∆P DGs (t) and ∆P ESS (t) are formulated as follows Using the model presented in Figure 1, we can find the relationship between the adjustment power provided by the ESS and the time when the ESS can fully respond: Considering Equations (4)-(6), Equation (3) can be integrated between t = 0 and t = t nadir : Assuming that before the contingency event, the system frequency is at the nominal value f norm , we have Using the model presented in Figure 1, we can find the relationship between the adjustment power provided by the ESS and the time when the ESS can fully respond: Considering Equations (4)-(6), Equation (3) can be integrated between = 0 and = : Assuming that before the contingency event, the system frequency is at the nominal value , we have Figure 1. Application of an energy storage system for fast frequency response.
Chavez et al. [14] show that the duration of is related to the governor dead-band (Hz) by equation = 1 . Besides, noting that the frequency nadir should not be below the predefined threshold , we obtain the frequency nadir' s requirement as follows. Chavez et al. [14] show that the duration of t d is related to the governor dead-band f db (Hz) by equation f db = 1 M H P t j t d . Besides, noting that the frequency nadir should not be below the predefined threshold f min , we obtain the frequency nadir' s requirement as follows.
Substituting Equations (6) into (7) and noting that P t adjust = P t dispost − P disch,t SS − P ch,t SS , we obtain the following constraint.
Equation (8) shows that the number of DGs in operation and their power output per hour is limited by the time taken for the ESS to fully react, and this will be used in the optimal scheduling formulation presented in Section 3.

Wind and Demand Models
In the present study, we focus on the optimal schedule of an islanded power system considering FFR. A challenge in this problem is the uncertainty in the expected wind power and demand. The wind power and demand forecast errors will affect the operation of the system, so the optimal scheduling problem must be formulated as a predictive optimization with results expressed as ranges of values that assure reliable operation of the system.
Both wind power and demand can be defined as the sum of the forecasted value and the forecasting error: The errors P t W−error and P t D−error are assumed to follow a normal distribution with zero mean, and σ t W is the standard deviation for wind power and σ t D for demand. This means the maximum error would be 3σ t W for wind power and 3σ t D for demand in correspondence to the confidence level of 99.7%.

The Optimal Scheduling Problem
This study implements day-ahead scheduling of the islanded power system including DGs, WG, and ESS with frequency criteria. The operating cost of DGs is minimized when wind power is used as much as possible. With the given forecasted wind power and demand, the operating mode of the DGs including ON/OFF state and power output is determined every hour for a 24-hour time horizon. The total energy excess or deficiency due to forecasting error is fully compensated by the ESS. Furthermore, the ESS must provide FFR if a generator is tripped. For this reason, the ESS should be scheduled such that it not only has enough energy to discharge but also can recharge to utilize the most energy from wind power.
The process of optimal scheduling comprises two stages, as illustrated in Figure 2. In the first stage, the deterministic day-ahead schedule in the commitment and the dispatch of DGs is decided and sent to the grid operator; the solid blue arrows in Figure 2 illustrate this process. The first-stage problem is a day-ahead UC problem that is implemented at least one day before the actual operation date. At that time, the wind and demand values are long-term forecast results with errors. Thus, only the results of the unit commitment (on/off states) and the power output of DGs are fixed. The ESS charge/discharge power and WG's power output are adjusted after wind power and demand are known with higher accuracy using a very short-term forecast. This two-stage optimal scheduling is formulated as a two-stage chance-constrained optimization model ( Figure 3). The constraints in the first stage refer to the deterministic planning, and the second stage ensures that the power balance and frequency criteria after a contingency event are met with a chosen probability.

Objective Function
The objective function to be minimized comprises the first-stage operating cost and the expected value of the second-stage cost: where represents a random vector including wind and demand; ℎ, ( ) and ℎ, ( ) are the ESS charge/discharge powers, respectively, decided corresponding to the actual values or the very short-term forecast value of .

First-Stage Constraints
The first stage is characterized by the following constraints.

Objective Function
The objective function to be minimized comprises the first-stage operating cost and the expected value of the second-stage cost: where represents a random vector including wind and demand; ℎ, ( ) and ℎ, ( ) are the ESS charge/discharge powers, respectively, decided corresponding to the actual values or the very short-term forecast value of .

First-Stage Constraints
The first stage is characterized by the following constraints.

Objective Function
The objective function to be minimized comprises the first-stage operating cost and the expected value of the second-stage cost: where ξ represents a random vector including wind and demand; P ch,t SS (ξ) and P disch,t SS (ξ) are the ESS charge/discharge powers, respectively, decided corresponding to the actual values or the very short-term forecast value of ξ.

First-Stage Constraints
The first stage is characterized by the following constraints. Active power balance constraint. The total active power output from the DGs P t i , the wind plant P t w , and the storage system (P disch,t SS and P ch,t SS ) must equal the given forecasted load P t D f at any time t: • DG operating constraints. The power output of each DG must be in the operating range between P mini and P maxi , which are specified by the manufacturer. The binary variable w t i in Equation (12) is used to keep the DG power output equal to zero if it is shut down. Equations (13) and (14) describe the minimum uptime τ iON and downtime τ iOFF limitations of each DG.
Equation (14) ensures that if at hour t the generator is offline (w t i = 0), the generator cannot start-up within the duration from In contrast, if at hour t the generator is online (w t i = 1), then this constraint ensures that the generator cannot shutdown within the duration from • Primary reserve constraints for DGs. This constraint shows that each DG can take part in the operating reserve.
• WG operating constraint. For a given forecasted wind power P t W f , the power output of the WG must satisfy Equation (16), where P t Wmin is the minimum wind power output: • ESS constraints. Equation (17) states that the charging and discharging power of the ESS must be smaller than the actual power rating of the storage device P SSmax . The process of charging/discharging the ESS is described by Equation (18). This constraint also imposes that the energy stored in the ESS should be smaller than its rated capacity at all times: • Frequency nadir limit. As presented in Section 2, to ensure that the frequency does not drop below the minimum allowable level, Equation (19) must be satisfied: Energies 2019, 12, 843 9 of 20 Note that the system inertia M H and ramp rate K depend on the number of the online DGs. Thus, this constraint determines not only ESS charge/discharge power, but also the number of DGs in operation and the power output of each DG before the contingency event to ensure the frequency nadir criteria.

•
Postcontingency energy storage capacity constraint. This constraint ensures that the ESS has enough energy to sustain FFR for ∆t FFR .

Second-Stage Constraints
In the second stage of the optimization model, uncertainties in wind power generation and load consumption are considered. The power output of the ESS and the WGs are redispatched as P disch,t SS (ξ), P ch,t SS (ξ), and P t w (ξ), where ξ represents a random vector. When a generator is lost, the ESS will discharge P t dispost (ξ) to decrease the disturbance to P t The second-stage constraints can be summarized as follows.
• Active power balance constraint: • WG operating constraint: • ESS constraints: • Frequency nadir limit: • Postcontingency energy storage capacity constraint: In the above constraints, Equation (21) guarantees that ESS and WG will be adjusted so that the probability of power imbalance is less than a risk level ε. Similarly, Equation (26) is the probability formulation of the Equation (19) and ensures that the frequency criterion will be met after a contingency event with high probability, even if the ESS is redispatched due to the difference between the long-term and very short-term forecast results of ξ.

The Modified Sample Average Approximation
Consider a simple two-stage chance-constrained optimization model where x is the first-stage variable, y is the second-stage variable, and ξ is random input data.
In this method, Monte Carlo simulation is used to approximate the distribution function of the random vector ξ by N samples. The optimization formulation in Equations (28) and (29) then becomes where 1 (0,∞) (G(x, y n , ξ n )) is an indicator function that is equal to one if G(x, y n , ξ n ) ≤ 0 and zero otherwise. It is assumed that the N samples have the same probability (1/N). This assumption helps to simplify the formulation of the optimization; however, a large number of samples are required to guarantee accuracy, which means the CPU time required to solve it increases accordingly.
In the present study, a modified approach to the SAA is proposed, by using a k-means clustering approach to reform the samples. Instead of using all N samples, the k-means clustering divides the samples into M clusters. The probability of each cluster is the sum of the probabilities of the constituent samples. Next, M centroids of the clusters are used as the SAA algorithm input samples, with the probability of each centroid being equal to the probability of the cluster that it represents. Figure 4 illustrates a small example: 1000 samples generated from the standard normal distribution N(0,1) are replaced by 10 centroids. In this method, Monte Carlo simulation is used to approximate the distribution function of the random vector ξ by N samples. The optimization formulation in Equations (28) and (29) then becomes subject to where (0,∞) ( ( , , )) is an indicator function that is equal to one if ( , , ) ≤ 0 and zero otherwise. It is assumed that the N samples have the same probability (1/ ). This assumption helps to simplify the formulation of the optimization; however, a large number of samples are required to guarantee accuracy, which means the CPU time required to solve it increases accordingly.
In the present study, a modified approach to the SAA is proposed, by using a k-means clustering approach to reform the samples. Instead of using all N samples, the k-means clustering divides the samples into M clusters. The probability of each cluster is the sum of the probabilities of the constituent samples. Next, M centroids of the clusters are used as the SAA algorithm input samples, with the probability of each centroid being equal to the probability of the cluster that it represents. Figure 4 illustrates a small example: 1000 samples generated from the standard normal distribution N(0,1) are replaced by 10 centroids.
subject to where is the probability of each centroid ( = 1,2, … ). For M centroids and their corresponding probabilities, Equations (28) and (29) are reformulated as where p m is the probability of each centroid (m = 1, 2, . . . M). Now let x and V, respectively, be the optimal solution and value of the optimal problem in Equations (32)-(33) and check whether this solution is feasible or not. Using Monte Carlo simulation to generate a new set of N samples where N is much larger than N, we find the value of the probability constraint in Equation (31) with solution x is The (1-ε)-confidence lower bound on g(x) is then computed using where Φ −1 is the inverse normal distribution function.
x is a feasible solution of the original problem only if L(x) ≥ 1 − ε. We repeated this process K times according to the flow chart illustrated in Figure 5, and found the maximum value U and minimum value L of the optimal value V. If the optimality gap given by U − L /U × 100% is smaller than a predetermined threshold, the algorithm terminates, and we obtain the optimal solution of the original problem. Now let ̅ and ̅ , respectively, be the optimal solution and value of the optimal problem in Equations (32)-(33) and check whether this solution is feasible or not. Using Monte Carlo simulation to generate a new set of ' samples where ' is much larger than , we find the value of the probability constraint in Equation (31) with solution ̅ is The (1-ɛ)-confidence lower bound on ( ̅ ) is then computed using where −1 is the inverse normal distribution function. ̅ is a feasible solution of the original problem only if ( ̅ ) ≥ 1 − . We repeated this process K times according to the flow chart illustrated in Figure 5, and found the maximum value ̅ and minimum value ̅ of the optimal value ̅ . If the optimality gap given by ( ̅ − ̅ ) ̅ ⁄ × 100% is smaller than a predetermined threshold, the algorithm terminates, and we obtain the optimal solution of the original problem.

Study System
The data used in this study is based on an actual power system on Phu Quy Island, Binh Thuan province, Vietnam. This system includes six 500 kW DGs and two WGs, each with a rated power of 1.8 MW. The purpose of the ESS installation is to utilize as much wind power as possible and to mitigate the uncertainty in demand and wind speed. The optimal size of an ESS taking into account the relevant wind scenarios and the annual load growth factor was analyzed in our previous work [26]. The ESS helps balance the fluctuation in wind power, increase the wind penetration level, and also provide FFR as a supplementary service. We assume that the wind power variable costs are zero and the ESS efficiency factor for the charging and discharging process is 90%. The risk level of the probability constraints in the second stage is 5%, which means these constraints should be met with a probability of greater than 95%. The other parameters are presented in Table 2. The scenario of wind and demand considered in this problem is described in Figure 6. The maximum possible instantaneous penetration of wind power is approximately 45% during the first hour and highest in the fifth hour (49%). However, this ratio is only 15% when the load is highest at the 19th hour. Forecast errors are modeled using a normal distribution with zero mean. The standard deviations are assumed to be 0.05 for both wind power and demand, which means that the maximum forecasting error in the values each hour can be considered to be 15%.

Study System
The data used in this study is based on an actual power system on Phu Quy Island, Binh Thuan province, Vietnam. This system includes six 500 kW DGs and two WGs, each with a rated power of 1.8 MW. The purpose of the ESS installation is to utilize as much wind power as possible and to mitigate the uncertainty in demand and wind speed. The optimal size of an ESS taking into account the relevant wind scenarios and the annual load growth factor was analyzed in our previous work [26]. The ESS helps balance the fluctuation in wind power, increase the wind penetration level, and also provide FFR as a supplementary service. We assume that the wind power variable costs are zero and the ESS efficiency factor for the charging and discharging process is 90%. The risk level of the probability constraints in the second stage is 5%, which means these constraints should be met with a probability of greater than 95%. The other parameters are presented in Table 2. The scenario of wind and demand considered in this problem is described in Figure 6. The maximum possible instantaneous penetration of wind power is approximately 45% during the first hour and highest in the fifth hour (49%). However, this ratio is only 15% when the load is highest at the 19th hour. Forecast errors are modeled using a normal distribution with zero mean. The standard deviations are assumed to be 0.05 for both wind power and demand, which means that the maximum forecasting error in the values each hour can be considered to be 15%. In the present study, the following cases are considered to evaluate the effectiveness of using ESS to provide FFR.

•
Case 1a. This case simulates the original UC problem without frequency constraints. The ESS is used to take advantage of wind power and compensate for the uncertainty in wind speed and demand. • Case 1b. The original UC problem with frequency constraints. The ESS is not used to provide FFR. • Case 2. UC problem with frequency constraints. The ESS is used to provide FFR. Besides comparing the frequency criteria after a contingency event in these cases, the effects of other parameters such as ESS size are also considered. The optimization problem is solved In the present study, the following cases are considered to evaluate the effectiveness of using ESS to provide FFR.

•
Case 1a. This case simulates the original UC problem without frequency constraints. The ESS is used to take advantage of wind power and compensate for the uncertainty in wind speed and demand. • Case 1b. The original UC problem with frequency constraints. The ESS is not used to provide FFR. • Case 2. UC problem with frequency constraints. The ESS is used to provide FFR.
Besides comparing the frequency criteria f nadir after a contingency event in these cases, the effects of other parameters such as ESS size are also considered. The optimization problem is solved using the MSAA approach presented in Section 4 with CPLEX version 12.6 and the YALMIP toolbox [27].

Case 1a: Original Optimal Scheduling Model Without Frequency Criteria
With an ESS rating of 400 kW/800 kWh and the other input data given in Section 5.1, the optimal daily schedule for DGs in Case 1 is presented in Figure 7a. Figure 7b shows the scheduled operation of DGs, wind generators, and the ESS. Note that the DGs operating schedules are first stage variables. The ESS charging/discharging power and wind power are second stage variables, considering the uncertain nature of wind speed and load forecast error. Therefore, they are represented using box plots. It can be seen that although the possible wind capacity and demand are uncertain, the available wind power is still fully utilized most of the time; this is undoubtedly due to the involvement of the ESS in the grid.
The RoCoF immediately after a contingency event in which the DG having the largest power output is lost is presented in Figure 8. During the period from the 17th to the 20th hour, demand is at its highest, whereas available wind power is quite low, so four DGs are online. This means that the stored kinetic energy in this period is higher than that in the rest of the day. However, the RoCoF is still approximately 10 Hz/s. Although DGs can increase their power output very quickly, even from a cold start condition (10-15 s), the frequency declines rapidly to below the minimum threshold. This can be explained by the fact that the inertia constant of the DGs being small (H = 0.8).
Energies 2018, 11, x FOR PEER REVIEW 13 of 20 using the MSAA approach presented in Section 4 with CPLEX version 12.6 and the YALMIP toolbox [27].

Case 1a: Original Optimal Scheduling Model Without Frequency Criteria
With an ESS rating of 400 kW/800 kWh and the other input data given in Section 5.1, the optimal daily schedule for DGs in Case 1 is presented in Figure 7a. Figure 7b shows the scheduled operation of DGs, wind generators, and the ESS. Note that the DGs operating schedules are first stage variables. The ESS charging/discharging power and wind power are second stage variables, considering the uncertain nature of wind speed and load forecast error. Therefore, they are represented using box plots. It can be seen that although the possible wind capacity and demand are uncertain, the available wind power is still fully utilized most of the time; this is undoubtedly due to the involvement of the ESS in the grid.
The RoCoF immediately after a contingency event in which the DG having the largest power output is lost is presented in Figure 8. During the period from the 17th to the 20th hour, demand is at its highest, whereas available wind power is quite low, so four DGs are online. This means that the stored kinetic energy in this period is higher than that in the rest of the day. However, the RoCoF is still approximately 10 Hz/s. Although DGs can increase their power output very quickly, even from a cold start condition (10-15 s), the frequency declines rapidly to below the minimum threshold. This can be explained by the fact that the inertia constant of the DGs being small (H = 0.8).    using the MSAA approach presented in Section 4 with CPLEX version 12.6 and the YALMIP toolbox [27].

Case 1a: Original Optimal Scheduling Model Without Frequency Criteria
With an ESS rating of 400 kW/800 kWh and the other input data given in Section 5.1, the optimal daily schedule for DGs in Case 1 is presented in Figure 7a. Figure 7b shows the scheduled operation of DGs, wind generators, and the ESS. Note that the DGs operating schedules are first stage variables. The ESS charging/discharging power and wind power are second stage variables, considering the uncertain nature of wind speed and load forecast error. Therefore, they are represented using box plots. It can be seen that although the possible wind capacity and demand are uncertain, the available wind power is still fully utilized most of the time; this is undoubtedly due to the involvement of the ESS in the grid.
The RoCoF immediately after a contingency event in which the DG having the largest power output is lost is presented in Figure 8. During the period from the 17th to the 20th hour, demand is at its highest, whereas available wind power is quite low, so four DGs are online. This means that the stored kinetic energy in this period is higher than that in the rest of the day. However, the RoCoF is still approximately 10 Hz/s. Although DGs can increase their power output very quickly, even from a cold start condition (10-15 s), the frequency declines rapidly to below the minimum threshold. This can be explained by the fact that the inertia constant of the DGs being small (H = 0.8).

Case 1b: Optimal Scheduling Problem Considering Frequency Criteria
Accounting for the frequency criteria in the optimal scheduling problem, as many DGs as possible are kept online while keeping their power output at a low level (Figure 9). For example, during peak load hours, there are six DGs in operation even though only four DGs are needed in Case 1a. This operating strategy helps increase the inertia of the system and reduce the RoCoF. The reduction of RoCoF can be seen by comparing the results in Figures 8 and 10. However, the frequency nadir f nadir at almost all hours is still much smaller than the minimum threshold f min (49.2 Hz). Moreover, increasing the power output of the DGs leads to higher operating costs (Table 3). Accounting for the frequency criteria in the optimal scheduling problem, as many DGs as possible are kept online while keeping their power output at a low level (Figure 9). For example, during peak load hours, there are six DGs in operation even though only four DGs are needed in Case 1a. This operating strategy helps increase the inertia of the system and reduce the RoCoF. The reduction of RoCoF can be seen by comparing the results in Figures 8 and 10. However, the frequency nadir at almost all hours is still much smaller than the minimum threshold (49.2 Hz). Moreover, increasing the power output of the DGs leads to higher operating costs (Table 3).   Table 3. Comparison of optimal costs between Cases 1 and 2.

Case 2: Optimal Scheduling Problem Considering Frequency Criteria with FFR Provided by ESS
As in Case 1a,b, an ESS rated at 400 kW/800 kWh is used in this case. However, the ESS is not only used to maintain the power balance but also to provide FFR. The optimal schedule of the DGs and box plots of the ESS and wind power for this case are presented in Figure 11.

Case 1b: Optimal Scheduling Problem Considering Frequency Criteria
Accounting for the frequency criteria in the optimal scheduling problem, as many DGs as possible are kept online while keeping their power output at a low level (Figure 9). For example, during peak load hours, there are six DGs in operation even though only four DGs are needed in Case 1a. This operating strategy helps increase the inertia of the system and reduce the RoCoF. The reduction of RoCoF can be seen by comparing the results in Figures 8 and 10. However, the frequency nadir at almost all hours is still much smaller than the minimum threshold (49.2 Hz). Moreover, increasing the power output of the DGs leads to higher operating costs (Table 3).   Table 3. Comparison of optimal costs between Cases 1 and 2.

Case 2: Optimal Scheduling Problem Considering Frequency Criteria with FFR Provided by ESS
As in Case 1a,b, an ESS rated at 400 kW/800 kWh is used in this case. However, the ESS is not only used to maintain the power balance but also to provide FFR. The optimal schedule of the DGs and box plots of the ESS and wind power for this case are presented in Figure 11.  As in Case 1a,b, an ESS rated at 400 kW/800 kWh is used in this case. However, the ESS is not only used to maintain the power balance but also to provide FFR. The optimal schedule of the DGs and box plots of the ESS and wind power for this case are presented in Figure 11.
As discussed in Section 2, if the DG with the largest power output is lost, the frequency deviation will activate the ESS response, and here we assume that the total response time for FFR is 100 ms. Based on the state of the ESS before the contingency event, the charge/discharge power of the ESS postcontingency at each hour can be in a range, as shown by the box plot in Figure 12. The frequency nadir values for each hour are obtained using Equation (7), as shown in Figure 13.
Energies 2018, 11, x FOR PEER REVIEW 15 of 20 As discussed in Section 2, if the DG with the largest power output is lost, the frequency deviation will activate the ESS response, and here we assume that the total response time for FFR is 100 ms. Based on the state of the ESS before the contingency event, the charge/discharge power of the ESS postcontingency at each hour can be in a range, as shown by the box plot in Figure 12. The frequency nadir values for each hour are obtained using Equation (7), as shown in Figure 13.    As discussed in Section 2, if the DG with the largest power output is lost, the frequency deviation will activate the ESS response, and here we assume that the total response time for FFR is 100 ms. Based on the state of the ESS before the contingency event, the charge/discharge power of the ESS postcontingency at each hour can be in a range, as shown by the box plot in Figure 12. The frequency nadir values for each hour are obtained using Equation (7), as shown in Figure 13.    As discussed in Section 2, if the DG with the largest power output is lost, the frequency deviation will activate the ESS response, and here we assume that the total response time for FFR is 100 ms. Based on the state of the ESS before the contingency event, the charge/discharge power of the ESS postcontingency at each hour can be in a range, as shown by the box plot in Figure 12. The frequency nadir values for each hour are obtained using Equation (7), as shown in Figure 13.    To evaluate the effect of the ESS on FFR, we compare the RoCoF immediately after the contingency event for Cases 2 and 1b. Figure 14 shows that the RoCoF in Case 2 is higher than in Case 1b in a few hours. However, the frequency nadir in Case 2 is ensured, while it is violated in several hours in Case 1b.
Note that Equations (19) and (26) in the optimization formulation limit the power output of each DG. This explains why the number of hours with six DGs in operation in Case 2 is more than that in Case 1b. In contrast, the ESS support in Case 2 helps to ensure the frequency criteria after contingency events, even when the number of online DGs is less than in Case 1, thus ensuring maximum utilization of the available wind power. This can be seen by comparing the box plots of wind power in Figures  9b and 11b. It is interesting to note that when the ESS is able to provide FFR, the UC solution will reduce DGs uptime and increase the wind power/ESS output, which in turn reduces the operating cost. The optimal cost of Case 2 is smaller than that of Case 1b and is not significantly higher than the nonconstrained optimal cost (Table 3). To evaluate the effect of the ESS on FFR, we compare the RoCoF immediately after the contingency event for Cases 2 and 1b. Figure 14 shows that the RoCoF in Case 2 is higher than in Case 1b in a few hours. However, the frequency nadir in Case 2 is ensured, while it is violated in several hours in Case 1b.
Note that Equations (19) and (26) in the optimization formulation limit the power output of each DG. This explains why the number of hours with six DGs in operation in Case 2 is more than that in Case 1b. In contrast, the ESS support in Case 2 helps to ensure the frequency criteria after contingency events, even when the number of online DGs is less than in Case 1, thus ensuring maximum utilization of the available wind power. This can be seen by comparing the box plots of wind power in Figures 9b and 11b. It is interesting to note that when the ESS is able to provide FFR, the UC solution will reduce DGs uptime and increase the wind power/ESS output, which in turn reduces the operating cost. The optimal cost of Case 2 is smaller than that of Case 1b and is not significantly higher than the nonconstrained optimal cost (Table 3).

Impact of ESS Size and Response Time
In this section, we consider the effect of the ESS size and the total response time from the moment the contingency event occurs until the ESS fully responds. Table 4 shows the smallest possible value of the frequency nadir for two total response times-100 ms and 200 ms-and several ESS sizes, which are defined by the power rating and the full charge/discharge duration ⁄ (from 0.5 h to 4 h). It can be seen that the ESS size must be larger than 200 kW/400 kWh to ensure the problem has a feasible solution. It should also be noted that the forecast errors are assumed to be ±15%, so too small an ESS will not be able to compensate for the mismatch between the predicted and actual values of the load and wind power. However, even if the optimization problem has a feasible solution, depending on the size of the ESS, there will still be a nonzero probability that the frequency nadir is lower than the minimum threshold ( Table 4). The reason for this is that the frequency nadir constraint in Equation (26) is formulated as a probabilistic constraint with a risk level of 5%. On the other hand, Equation (26) shows that a longer response time requires a lower power output from each DG or more DGs in operation to provide enough kinetic energy; consequently, increasing the ESS power rating is necessary. It can be seen from Table 4 that when the response time is 200 ms, the power rating of the ESS must be larger than 600 kW to maintain the frequency nadir above 49.2 Hz, whereas an ESS with rated power 400 kW is acceptable if the response time is 100 ms.

Impact of ESS Size and Response Time
In this section, we consider the effect of the ESS size and the total response time from the moment the contingency event occurs until the ESS fully responds. Table 4 shows the smallest possible value of the frequency nadir for two total response times-100 ms and 200 ms-and several ESS sizes, which are defined by the power rating P SSmax and the full charge/discharge duration E SSmax /P SSmax (from 0.5 h to 4 h). It can be seen that the ESS size must be larger than 200 kW/400 kWh to ensure the problem has a feasible solution. It should also be noted that the forecast errors are assumed to be ±15%, so too small an ESS will not be able to compensate for the mismatch between the predicted and actual values of the load and wind power. However, even if the optimization problem has a feasible solution, depending on the size of the ESS, there will still be a nonzero probability that the frequency nadir is lower than the minimum threshold ( Table 4). The reason for this is that the frequency nadir constraint in Equation (26) is formulated as a probabilistic constraint with a risk level of 5%. On the other hand, Equation (26) shows that a longer response time requires a lower power output from each DG or more DGs in operation to provide enough kinetic energy; consequently, increasing the ESS power rating is necessary. It can be seen from Table 4 that when the response time is 200 ms, the power rating of the ESS must be larger than 600 kW to maintain the frequency nadir above 49.2 Hz, whereas an ESS with rated power 400 kW is acceptable if the response time is 100 ms. Case 2 was solved by the traditional SAA algorithm and the MSAA algorithm to compare their computational efficiency. To solve Case 2, the SAA algorithm must be repeated at least 50 times per loop (K = 50) and needs at least 100 samples per loop (N = 100). On the other hand, the original set of 1000 samples can be replaced with five centroids, and the MSAA algorithm needs 50 iterations to obtain the results. Interestingly, five centroids in the MSAA algorithm are equivalent to five samples in the SAA algorithm; thus, it is easy to see that the computing time required for the MSAA algorithm is much smaller than that of the SAA algorithm ( Table 5).
The performance of the proposed MSAA is also compared with that of the Robust chance-constrained formulation, which is also a popular approach.
In this comparison, a two-stage robust chance-constrained model, solved by column-and-constraint generation algorithm (CCG) [28][29][30], is implemented. The constraints related to power balance and frequency criteria are also formulated as probability constraints with the same risk level. Besides, the results obtained with a two-stage robust model is also shown. The results in Table 5 clearly demonstrate the compromise between CPU time and economic performance: although the required CPU time for MSAA is longer than the robust method, the optimal cost obtained by MSAA is significantly lower than both robust models-with or without chance constraints.

Conclusions
In this paper, an optimal day-ahead scheduling problem concerning the application of ESS for FFR is considered and analyzed in detail. The optimization problem is formulated within a two-stage chance-constrained framework, in which the load and the maximum possible wind power are uncertain. In this model, power balance and frequency criteria constraints are formulated as probability constraints with a certain risk level. Based on the first-order model of frequency dynamic, the relationships between the power output of each DG, the ESS charge/discharge power, and the response time are studied. The impact of the size and response time of the ESS on the frequency nadir after the sudden loss of a DG is also analyzed. It is also noteworthy that an MSAA approach was proposed in the present study to solve a chance-constrained problem, and the effectiveness of this method was demonstrated.
The results obtained in two cases-with and without FFR provided by ESS-demonstrate the effectiveness of FFR in arresting frequency deviations after a contingency event. The proposed method ensures that the minimum frequency threshold is not violated, even when the actual values of wind power and demand are different from the predicted values incorporating the predetermined maximum errors. The results also show that a slower FFR will lead to a larger ESS to ensure frequency criteria.
The proposed approach can be extended to consider multiple contingencies such as line outages or load interruptions as well as equipment failures. The model can also be readily adapted to include other uncertain factors, such as solar power generation or electricity prices. These topics are left for future work.

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