When a disturbance occurs, SG will respond according to the rotor motion equation under the action of unbalanced torque (Anderson and Mirheydar, 1990): 2 H s d Δ f t d t = Δ P s t − Δ P e t − D s Δ f t (1)

In the formula, Δ f t , Δ P s t , and Δ P e t are the frequency deviation, mechanical power change, and external power change, respectively; H s and D s are the inertial time constant and damping constant of the system, respectively.

The generator governor will respond to power imbalance, and its characteristics can be described as follows: Δ P s = − 1 R s 1 + s F H T R 1 + s T R Δ f (2)

In the formula, R s is the static adjustment coefficient, F H is the output power ratio of the high-pressure cylinder, and T R is the volume effect time constant of the reheater.

In the frequency response model of multi-SG system, the capacity difference between units makes the disturbance power shared between units different, which will change at different speeds, and the influence weight on the system is expressed by the proportion coefficient k s i (Huang et al., 2020): k s i = S s i / S s y s (3)

In the formula, k s i is the proportion coefficient, which indicates that the rated capacity of the ith SG accounts for the part of the whole system capacity. Where S s i is the rated capacity of the ith SG, and S s y s is the sum of the system capacity.

When the disturbance occurs, the wind turbine will respond according to the rotor motion equation under the action of unbalanced torque (Sudipta et al., 2016; Huang et al., 2022): 2 H w d Δ f t d t = Δ P w t − Δ P e t − D w Δ f t (4)

In the formula, H w and D w are the virtual inertia and virtual damping provided by the WF, and Δ P w t are the output power changes of the WF.

The governor dynamic description is as follows: Δ P w = − 1 R w 1 1 + s T w Δ f (5)

In the formula, R w is the droop coefficient and T w is the response time constant of WF.

Similarly, in the frequency response model of the multi-WF system, the influence weight of wind farms with different capacities on the system is expressed as a proportion coefficient k w j : k w j = S w j / S s y s (6)

In the formula, k w j is the proportion coefficient, which represents the part of the rated capacity of the jth WF relative to the whole system. Where S w j is the rated capacity of the jth WF.

The key to temperature control loads participating in system frequency regulation is to respond to the change in system frequency through its power consumption (Conte et al., 2021).

The temperature control load used in this paper is a typical temperature-controlled household electric water heater, and the electric heating wire is a pure resistive electric load. The purpose of frequency control can be achieved by changing the load power. When a disturbance causes the system’s frequency to fluctuate, the power consumption is adjusted by using the linear relationship between the target temperature change value and the frequency deviation. As a result, the adaptive load regulation of an electric water heater contributes to frequency response. Δ T = K ATL Δ f T max * = T max 0 + Δ T T min * = T min 0 + Δ T (7)

In the formula, K ATL is the frequency regulation coefficient, T max * and T min * are the upper and lower limits of the water temperature after the change. Δ P ATL = 0 , Δ T ⩽ T min * − T max * P atl N atl G on + Δ T T max * − T min * P atl N atl G on , T min * − T max * < Δ T < 0 P atl N atl G on + Δ T T max * − T min * P atl N atl G off , 0 ⩽ Δ T < T max * − T min * P atl N atl , Δ T ⩾ T max * − T min * (8)

In the formula, G on and G off are the proportion of electric water heaters that are on or off during stable operation, Δ P ATL is the power of the electric water heater polymerization load participating in frequency control, N atl is the number of electric water heaters, and P atl is the rated power of a single electric water heater.

The larger the K ATL , the larger the Δ T , the more the electric water heater participates in the frequency regulation and the lower the user’s comfort level. Therefore, the balance between the frequency regulation effect and the user’s comfort level should be considered in the actual control.

The way temperature control load aggregation participates in the primary frequency regulation of the power grid is frequency droop control. The demand response aggregator acts as a virtual generator to respond to the frequency modulation demand, and responds according to the rotor motion equation under unbalanced torque: 2 H a d Δ f t d t = Δ P A t − Δ P e t − D a Δ f t (9)

In the formula, H a and D a are the inertia constant and damping constant of ATL, and Δ P A t is the power change of ATL.

Here, the governor of the first-order inertia link is considered, so the equivalent transfer function of temperature control load participating in primary frequency regulation can be expressed as: Δ P A = − 1 R a 1 1 + s T a Δ f (10)

In the formula, Δ P A is the power of ATL participating in frequency modulation, R a is the droop coefficient, and T a is the response time constant of ATL.

Similarly, in the frequency response model of multi ATL system, the influence weight of polymerization units with different capacities on the system is represented by the proportion coefficient k a k : k a k = S a k / S s y s (11)

In the formula, k a k is the proportion coefficient, which represents the part of the rated capacity of the kth ATL relative to the whole system. Where S a k is the rated capacity of the kth ATL.

BESS can achieve millisecond response, and the time scale of primary frequency regulation is typically seconds, so the frequency response function of BESS can be expressed by a delay function and a first-order inertial link: P e s s = P ES s e − s T ES 1 1 + s T e s (12)

In the formula, P ES is the power of BESS, T ES is the communication delay of BESS, and T e s is the emergency control response time of BESS.

The control function of the load can be expressed by a delay function with time constant T L : P l s = P L s e − s T L (13)

In the formula, P L is the power reduction of the load, and T L represents the time interval between the fault occurrence time and the load shedding action time.

They are considering the frequency response characteristics of various resources, including synchronous machines, wind farms, battery energy storage systems, temperature control loads, and conventional loads. The multi-resource SFR model is established, as shown in Figure 1. The frequency variation in the power grid can usually be regarded as the unity of the whole system, so the inertia of each unit of SG, WF, and ATL can be superimposed according to the unified standard during synchronous operation, which is equivalent to one unit providing the overall inertia of the system, and the response time constant is very small, which is simplified after ignoring. The parameters of the reduced order model are as follows: κ s i = k s i / R s i , κ w j = k w j / R w j , κ a k = k a k / R a k 1 / R = ∑ i = 1 N SG κ s i + ∑ j = 1 N WF κ w j + ∑ k = 1 N ATL κ a k λ = κ s i / 1 / R , 0 ≤ i ≤ N SG κ w j / 1 / R , 0 ≤ j ≤ N WF κ a k / 1 / R , 0 ≤ k ≤ N ATL H = ∑ i = 1 N SG λ i H s i + ∑ j = 1 N WF λ j H w j + ∑ k = 1 N ATL λ k H a k D = ∑ i = 1 N SG λ i D s i + ∑ j = 1 N WF λ j D w j + ∑ k = 1 N ATL λ k D a k F H = ∑ i = 1 N SG λ i F H i , T R = ∑ i = 1 N SG λ i T R i (14)

In the formula, N SG and N WF are the number of SGs and WFs, H s i and D s i are the inertia and damping constants of the ith SG, and H w j and D w j are the virtual inertia and virtual damping constants of the jth WF. κ s i and κ w j are the branch equivalent gains of SG and WF. λ is the weighted coefficient of the branch. H and D are the system inertia and damping constants after polymerization, F H is the output power ratio of the high-pressure cylinder after polymerization, and T R is the volume effect time constant of the reheater after polymerization.

Since the response time of BESS and load is fast, it is assumed that the time constant can be ignored. The parameters of the reduced model are as follows: P ES = ∑ m = 1 N ES P ES m P L = ∑ n = 1 N L P L n (15)

In the formula, P ES m and P L n are the power of the mth BESS and the nth load, and N ES and N L are the number of BESS and load involved in emergency control.

In summary, the reduced-order model is shown in Figure 2.

The change of the power shortage of the system can be expressed as: Δ P e t = L − 1 − Δ P lost s + P ES s + P L s = − Δ P lost + P ES + P L (16)

In the formula, Δ P lost is the active power shortage in the system. P ES and P L are the power of BESS and load.

According to the transfer function in Figure 2, the expression of Δ f can be written as: Δ f = R ω n 2 D R + 1 1 + T R s P e s 2 + 2 ζ ω n s + ω n 2 (17)

Considering that the step type disturbance in the power system is the most common and has the greatest impact, such as generator cut-off, sudden load increase, et cetera, assume that P e is the unit step function: P e s = P e s (18)

The expression of frequency change Δ f t can be obtained by combining the inverse Laplace transform of formula (17) and formula (18): Δ f t = P e ⋅ U t U t = R D R + 1 1 + α e − ζ ω n t ⁡ sin ω r t + φ (19)

Where ω n , ζ , ω r , α and φ are the coefficients of the governor, as follows: ω n = D R + 1 2 H R T R ζ = D R T R + 2 H R + F H T R 2 D R + 1 ω n ω r = ω n 1 − ζ 2 α = 1 − 2 T R ζ ω n + T R 2 ω n 2 1 − ζ 2 φ = arctan ω r T R 1 − ζ ω n T R − arctan 1 − ζ 2 − ζ (20)

According to the superposition principle of linear system, the frequency deviation Δ f can be obtained as follows: Δ f = Δ P e + Δ P A ⋅ U t (21)

Combined with Formulas (16–20), the analytical formulas of the frequency nadir f nadir and the time to reach frequency nadir t nadir can be derived as follows: f nadir = f 0 ⋅ 1 + Δ P e ⋅ U t nadir t nadir = 1 ω r tan − 1 ω r T R ζ ω n T R − 1 (22)

In the formula, f 0 is the frequency of the steady-state operation of the system.

After the frequency of the system passes through the lowest frequency point, the primary frequency regulation makes the system return to the quasi-steady state frequency. The expression of quasi-steady state frequency is: f qss = f 0 ⋅ 1 + Δ P e ⋅ U t ∞ (23)

In the formula, f qss is the quasi-steady state frequency, Δ f ′ is the frequency deviation at the quasi-steady state frequency, and t ∞ is the time when the system enters the quasi-steady state, usually tens of seconds.

The analytical expressions of frequency minima and quasi-steady state frequencies have a linear relationship with the decision variables. Frequency constraints include frequency minima and quasi-steady state frequency constraints. f nadir ≥ f nadir min f qss ≥ f qss min (24)

In the formula, f nadir min and f qss min are the safety thresholds of the frequency nadir and the quasi-steady state frequency.

During regular operation, the system frequency deviation is required to be controlled at ±0.2 Hz to set an appropriate constraint boundary.

To prevent the control amount of each resource exceeding its own limits, it is necessary to constrain the control amount of each resource. 0 ≤ P ES , m ≤ P ES , m max − P ES , m 0 , ∀ m ∈ N ES 0 ≤ P L , n ≤ P L , n 0 − P L , n min , ∀ n ∈ N L (25)

In the formula, P ES , m max and P ES , m 0 are the maximum output power and pre-fault output power of the mth BESS. P L , n 0 and P L , n min are the pre-fault load and the minimum reserve after load shedding of the nth load.

The power flow constraints are set on all transmission lines to ensure the power flow on transmission lines does not exceed their limits, which are as follows: ∑ m = 1 N ES P ES , m λ l ES , m + ∑ n = 1 N L P L , n λ l L , n + P l 0 ≤ P l max , ∀ m ∈ N ES , ∀ n ∈ N L (26)

In the formula, P l 0 is the existing power flow on lth line, and P l max is the maximum power flow on lth line. λ l ES , m and λ l L , n are the power flow transmission ratio of the mth BESS and the nth load on the lth line.

The EFC decision optimization aims at the minimum control cost, which is a mixed integer linear programming (MILP) problem considering frequency constraints. The decision variables are controllable electric energy storage output power P ES and load shedding P L . min ⁡ F = a ES ∑ m = 1 N ES C ES , m P ES , m + a L ∑ n = 1 N L C L , n P L , n over P ES , m , ∀ m ∈ N ES P L , n , ∀ n ∈ N L f nadir , f qss s . t .   f nadir ≥ f nadir min f qss ≥ f qss min 0 ≤ P ES , m ≤ P ES , m max − P ES , m 0 , ∀ m ∈ N ES 0 ≤ P L , n ≤ P L , n 0 − P L , n min , ∀ n ∈ N L ∑ m = 1 N ES P ES , m λ l ES , m + ∑ n = 1 N L P L , n λ l L , n + P l 0 ≤ P l max , ∀ m ∈ N ES , ∀ n ∈ N L (27)

In the formula, P ES , m and P L , n are the power of the mth BESS and the nth load. a ES and a L are the weight coefficients of BESS and load. C ES , m and C L , n are the control cost per kW of the mth BESS and the nth load.

The weight coefficient is related to the priority of the control resources, and the resources with smaller weight coefficients are scheduled preferentially. Dispatchers can set the weight coefficients of different resources according to the control requirements of the system. In the actual power system, the control cost of BESS is usually low, and the cost of load shedding is the highest.

This section initiates the verification of the accuracy of the frequency response reduction model and conducts a frequency analysis of the multi-resource system through simulation. To facilitate this verification, a three-machine nine-node system, as depicted in Figure 4, is employed as the simulation framework. The system’s initial state is configured to be in normal operation, and a disturbance, in the form of a 0.1 per unit (p.u.) active power deficiency, is introduced at 5 s. The center of inertia (COI) frequency is employed as the representative metric for the system’s frequency behavior during time-domain simulations. The simulation results, presented in Figure 5, provide a visual representation of the system’s response to the introduced disturbance.

Analysis of the obtained results, as illustrated in Figure 5 and summarized in Table 1, reveals a notable congruence between the COI frequency curve derived from time-domain simulations and the frequency curve predicted by the SFR transfer function model developed in this study. Additionally, the results computed using the frequency analytical formula align closely with the outcomes of the time-domain simulations. This consistency underscores the accuracy and precision of both the SFR transfer function model introduced in this paper and the associated frequency analytical formula. These modeling approaches, in conjunction with the time-domain simulation, exhibit minimal discrepancies, thus affirming their reliability and robustness.

The modified IEEE 13 bus distribution system is adopted for the test, and the structure is shown in Figure 6. SG, WF, BESS, and ATL are set, respectively. The system’s initial state is set to operate normally, and then the power distribution system is disconnected from the superior power grid due to the accident, resulting in the disturbance of active power shortage. In order to illustrate the feasibility and effectiveness of the EFC strategy, the influence of the EFC strategy decision on the primary frequency regulation of the system will be compared.

Two test scenarios will be set: (1) Set the total capacity of BESS participating in frequency modulation to 0.03p.u. and the total capacity of removable load to 0.05p.u. The system initially operates normally. At 10 s, the active power shortage of 0.12p.u. occurs suddenly in the stable system. (2) Set the total capacity of BESS participating in frequency modulation as 0.05p.u. and the total capacity of removable load to 0.05p.u. The system initially operates normally. At 10 s, the active power shortage of 0.1p.u. occurs suddenly in the stable system.

Figure 7 and Table 2 provide a clear representation of the system’s frequency transient change process. Without the application of the EFC strategy, the system’s lowest frequency exhibits a significant decline to 49.56 Hz, with the quasi-steady state frequency stabilizing at 49.72 Hz.

In Scenario (1), the EFC strategy is implemented, resulting in an initial energy storage system output of 0.03 p.u. This prioritization of the energy storage system leads to the utilization of its entire output capacity. However, due to the limited energy storage system capacity, this output alone is insufficient to restore the system frequency to a safe operating level. Consequently, a load shedding of 0.035 p.u. becomes necessary. Subsequently, the lowest system frequency rises to 49.80 Hz, and the quasi-steady state frequency increases to 49.87 Hz. In contrast, in Scenario (2) with the EFC strategy, the energy storage system exhibits an output of 0.043 p.u. Since this output remains within the energy storage system’s operational limits, no load shedding is required. Consequently, the lowest system frequency is restored to 49.80 Hz, with the quasi-steady state frequency reaching 49.87 Hz.

The test results confirm the EFC strategy’s effectiveness in raising both the lowest and quasi-steady state frequencies to levels within the defined safety thresholds. Furthermore, it successfully allocates the output of controllable resources in an efficient manner. This empirical validation underscores the feasibility and efficacy of the proposed EFC strategy.

Our approach to emergency frequency control aims to address the limitations of traditional methods that predominantly rely on the frequency modulation capability of synchronous machines or involve various combinations of synchronous machines, wind farms, and energy storage systems. These conventional strategies, while effective to a certain extent, may result in load-shedding if the system’s frequency modulation ability falls short of requirements. In the context of distribution systems, load-shedding can lead to an unsatisfactory user experience and have significant societal implications, especially when essential loads are disconnected. In our proposed method, we have introduced an innovative approach by incorporating temperature control loads and optimizing the coordination of multiple power sources within the distribution system. This innovation allows us to aggregate these resources effectively, thereby enhancing the system’s frequency regulation capabilities. The ultimate goal is to minimize or eliminate load-shedding, which is advantageous for both the distribution system’s overall performance and the satisfaction of its users.

The modified IEEE 33 bus distribution system is adopted for the test, and the structure is shown in Figure 8. SG, WF, BESS, and ATL are set at two places, respectively. The system initially operates normally. Then, the distribution system is disconnected from the superior power grid due to the accident, resulting in the disturbance of active power shortage. In order to illustrate the feasibility and effectiveness of multi-resource regulation, the effect of EFC strategy with or without multi-resource regulation on the primary frequency regulation of the system will be compared.

The system initially operates normally. However, at 10 s, the active power shortage occurs suddenly in the stable system. To assess the impact of frequency modulation by different resource combinations, we establish the following test scenarios:

Scenario 1

A sudden active power shortage of 0.1 p.u. occurs at 10 s in the stable system. We evaluate the frequency modulation effects using the following combinations:

Case 1

SG only

Case 2

SG and WF

Case 3

SG, WF, and ATL

Case 4

SG, WF, ATL, BESS, and load

Scenario 2

A sudden active power shortage of 0.5 p.u. occurs at 10 s in the stable system. We assess the frequency modulation effects with the same combinations as in Scenario 1.

These scenarios enable a comprehensive examination of the frequency modulation capabilities of the specified resource combinations under varying degrees of active power shortage in the power system.

Analysis of the system’s primary frequency response, as depicted in Figure 9 and summarized in Table 3, reveals the impact of various resource combinations on system frequency during disturbance scenarios. In scenario (1), where only the SG is engaged in frequency modulation, the system’s lowest frequency is observed at 49.59 Hz, with a quasi-steady state frequency of 49.76 Hz. This initial assessment underscores the limitation of relying solely on a single resource for maintaining frequency stability within the distribution system during disturbances. Upon introducing WF in conjunction with SG for frequency modulation, we observe an improvement in the system’s primary frequency response. The lowest frequency in this case is 49.68 Hz, with a quasi-steady state frequency of 49.77 Hz, demonstrating the beneficial impact of combining resources. Expanding the resource pool further to include ATL alongside SG and WF for frequency modulation results in a more robust primary frequency modulation capability. The system’s lowest frequency in this configuration is measured at 49.73 Hz, with a quasi-steady state frequency of 49.78 Hz. However, it is still not enough to ensure that the system frequency is within the safety threshold until SG, WF, ATL, BESS, and load-shedding participation in frequency modulation. The lowest frequency observed in this comprehensive scenario is 49.80 Hz, with a quasi-steady state frequency of 49.84 Hz.

Analysis of the system’s primary frequency response, as depicted in Figure 10 and summarized in Table 4. In scenario (2), where only the SG is engaged in frequency modulation, the system’s lowest frequency is observed at 48.18 Hz, with a quasi-steady state frequency of 48.86 Hz. Upon introducing WF alongside SG for frequency modulation, the lowest system frequency increases to 48.48 Hz, with a quasi-steady state frequency of 48.89 Hz. Furthermore, the addition of ATL to frequency modulation raises the lowest frequency to 48.78 Hz, with the quasi-steady state frequency reaching 48.96 Hz. Subsequently, the involvement of BESS and load-shedding participation in frequency modulation further enhances the system’s frequency stability. In this configuration, where SG, WF, ATL, BESS, and load-shedding collectively engage in frequency modulation, the lowest system frequency significantly improves, measuring 49.80 Hz, while the quasi-steady state frequency closely follows at 49.83 Hz.

The test results have indeed conducted a comparative analysis of various resource combinations, including those found in existing methods, show that the regulation of multi-resource in the distribution system can improve the primary frequency regulation ability of the distribution system in response to disturbances. Moreover, the EFC strategy can improve the system’s lowest frequency and quasi-steady state frequency within the safety threshold, ensuring the safe and stable operation of the distribution system in the event of an accident. These results unequivocally demonstrate that our proposed strategy outperforms the alternatives, delivering superior outcomes in terms of system frequency regulation during emergency scenarios.