*2.1. Preliminaries*

For studying the dynamic response of power systems dealing with frequency variations, traditional system modeling has been based on the following per unit (pu) expression [5]:

$$
\Delta P\_{\overline{G}} - \Delta P\_L = D \Delta f + 2H \frac{d \Delta f}{dt},\tag{1}
$$

where Δ *PG* − Δ *PL* is the power imbalance, Δ *f* is the difference between instantaneous and nominal system frequency and *D*, referred to as the damping factor, is the load dependence on frequency. Assuming that grid frequency can be considered as constant over large interconnected areas, the generating units can be then combined into an equivalent rotating mass *M* (being *M* = 2*H* and *H* the inertia constant expressed in seconds). Likewise, loads are grouped and considered as an equivalent load, being *D* their equivalent damping factor [37].

For multi-area analysis purposes, the previous expression can be extended to different areas. With this aim, Figure 1 shows a general scheme of two-area interconnected power system, including primary and secondary frequency control by conventional generation units. This multi-area general scheme is considered in this paper, where a grid frequency and a system dynamic response characterizes each area. A high-impedance (elastic) transmission line is proposed to connect both areas [5]. Conventional generating units—highlighted with dashed lines—are modeled by considering Generation Rate Constraint (GRC) and speed Governor Dead-Band (GDB). The GRC is modeled by considering an "open loop" method and by adding two limiters bounded by ±0.0017 (pu/s) within the generating units [38]. The GDB is defined as the total magnitude of a sustained speed change within which there is no change in the turbine valve position. The GDB transfer function model can be found in [38], assuming that the GDB of the 0.06% backlash type can be linearized in terms of change and the rate of change in the speed. Under the presence of certain nonlinearities and constraints, such as GRC and GDB, system dynamic responses can present significant overshoots and long settling times for frequency and tie-line power oscillations [39,40].

**Figure 1.** General scheme of two-area interconnected power system.

If one area is under power imbalance, the governor control mechanisms of both areas modify the power system to recover power balance. Grid frequency deviation, Δ *f* , is used as input signal for primary and secondary frequency control. PFC is performed locally at the generator, being the active power increment/decrement proportional to Δ *f* through the speed regulation parameter *R*, defined as the ratio of Δ *P* and Δ *f* . Power flow control at the tie-line and damping of tie-line power oscillations are required for effective active power generation and frequency control. In this respect, SFC involves an integral controller modifying the turbine set-point of each area. A linear combination of both tie-line power errors between neighboring areas, Δ *Ptie*,*ij*, and local frequency variations in each area, Δ *fi*, is used as the input to the corresponding integral controllers, which is called the Area Control Error (ACE) [41]. ACE in each area is then defined as follows [41]:

$$A \mathbb{C} E\_i = \sum\_{j=1}^{n} \Delta P\_{ti \varepsilon, ij} + \beta\_i \Delta f\_{i\prime} \tag{2}$$

where the suffix *i* refers to the CA and *j* to the generator number; *βi* is the bias coefficient; *fi* is the grid frequency of the *i*-CA; and Δ *Ptie*,*ij* is the actual value of the interchange power between *i*-CA and *j*-CA. Further information regarding tie-line bias control applicability to load frequency control for multi-area interconnected power systems can be found in [42]. The dynamic performance of the AGC system thus depends on frequency bias factor *β<sup>i</sup>*, in MW/Hz, and the integral controller gain value *KI*. Optimal values of *KI* and *β* are estimated by means of an Integral Squared Error (ISE) technique,

$$ISE = \int\_{0}^{T} \left(\Delta P\_{tic}^{2} + \left(\beta\_{1}\Delta f\_{1}\right)^{2} + \left(\beta\_{2}\Delta f\_{2}\right)^{2}\right) dt,\tag{3}$$

where <sup>Δ</sup>*f*1 and <sup>Δ</sup>*f*2 are frequency deviations in Areas 1 and 2, respectively [43,44]. According to the specific literature, two different power systems were simulated: *(i)* the two-area interconnected power system proposed in [10]; and *(ii)* a two-area wind farm integrated power system based on [43–45]. Figure 2 depicts the initial two-area interconnected power system considered for the simulation by including the corresponding block diagrams.

**Figure 2.** Block diagram of two-area interconnected power system.

#### *2.2. Wind Power Modeling*

Accounting for the relevance of wind power integration, two WFs containing VSWTs are explicitly considered in the power system model. Averaged time-varying wind speed profiles are applied to both wind power plants to evaluate wind active power fluctuations on grid frequency stability.

Primary frequency response of WFs is improved by including an additional inertia controller. This controller allows VSWTs to contribute to the system inertia through an active power regulation under frequency deviations [25,35],

$$P\_{cl}^{correction} = k\_{\rm WFs} \cdot \frac{I\_{\rm WT}}{2} \frac{df}{dt} \,\tag{4}$$

where *k*WFs is the virtual inertia factor and *I*WT is the WT inertia (pu). Larger *k*WFs-factors would improve WT contribution, maintaining WT rotational speed values within acceptable ranges. Figure 3 shows this additional virtual controller depending on the RoCoF (*d f* /*dt*). An additional control loop that synthesizes virtual inertia is implemented for the VSWTs to contribute to the system's total inertia and thus to the system frequency control. The output power of the WF, *Pref el*(*pu*), is modified by an additional power *Pcorrection el*(*pu*) depending on the RoCoF, *d f* /*dt*, [30,46]. The global power of the WFs (Δ*PWFs*(*pu*)) depends on the number of WT (*NWTs*) and the VSWTs power set-point (*P*0,*WT*) [47,48]. The proposed approach is thus based on a modified inertial control scheme, involving quick response through power-electronic circuits and subsequently providing frequency support from the rotational mass kinetic energy.

**Figure 3.** Virtual Inertia Block (VIB) diagram of the virtual frequency controller.

## *2.3. Demand-Side Modeling*

Even thoguh demand-side was initially considered as non-controllable, nowadays, most authors agree with the idea that the devices can be categorized into controllable or uncontrollable [49]. Controllable loads are then defined as the loads suitable for deferral or shifting their power demand to off-peak periods without harming household's convenience, such as heating, ventilation, and air conditioning [50,51]. Because of their thermal properties, these devices can adapt their operation through on/off signals without adversely affecting their temperature requirements. The rest of the loads are considered uncontrollable, i.e., their use is completely determined by the end user. Assuming uncontrollable loads to be modeled by an equivalent load, their global power depends on frequency excursions through parameter *D* (equivalent damping factor), as discussed in Section 2.1.

In the proposed solution, controllable loads are modeled individually taking into account that their individual power demand can be modified by the following controllers: *(i)* a proportional frequency-dependent controller able to modify the active power demanded by the device and thus emulating PFC from supply-side; and *(ii)* an integral-action controller able to modified the thermostat individual load settings based on the ACE signal. The proportional frequency-dependent controller was previously described by the authors and discussed in [21,25]. This approach modifies the controllable load demand through forced switching-off/on actions. To operate under forced disconnections, frequency excursions exceed a predefined threshold, Δ *f* , for a certain time (see Figure 4). Larger frequency deviations imply faster frequency control responses, thus entering the *Control Region*. Small frequency changes are maintained during longer time intervals, delaying the controller actions. Figure 4 also gives an example of Δ *f-time* feature for a given load and its corresponding demand profiles in the event of a linear under-frequency excursion. As can be seen, when Δ *f-time* boundaries are exceeded, load is forcibly switched-off for a predefined time period (*OFF-forced*). After a short recovery interval of forced connection (*ON*) to maintain the temperature within an acceptable range, and given that the frequency does not recover its rated value thus continuing within the *Control Region*, the load controller restarts a new cycle of forced disconnection. These *ON-* and *OFF-forced* time periods are designed keeping both appliance requirements in view and considering temperature limits that may be almost imperceptible by the customers. Actually, most controls for conditioning spaces and typical thermostat-controlled devices operate on a time scale of 10–15 min [52], which is in line with the SFC timescales. The proposed integral-action controller is thus based on thermostat temperature changes according to comfort level constraints and the expected demand-side contribution to SFC. Further discussion about thermal inertia, heat–cool flow rates and SFC can be found in Section 3.

**Figure 4.** Frequency-responsive load controller: example of demand-side contribution to PFC and Δ*f -time* characteristic.

Figure 5 shows an extended block diagram by including wind power generation frequency control (Δ*PWFs*) and demand-side contribution to PFC and SFC considered in Area 1 (Δ*PSFC DS* ). As can be seen, supplementary power due to wind power generation frequency control, Δ*P*WFs, is added to the conventional generation response (prime movers), Δ*PG*1, to provide an additional generation to the supply-side. The demand-side contributions to PFC and SFC are referred to as Δ*P PFC DS* and Δ*P SFC DS* , respectively. This multi-area power system model is an extension of the block diagram depicted in Figure 2, by including frequency control to both supply-side—wind power plants—and demand-side.

**Figure 5.** Two-area interconnected power system. WF and demand-side contribution to frequency control.

#### **3. Demand-Side Contribution to SFC: Proposed Solution**

## *3.1. General Description*

The implementation of demand-side contribution to SFC is based on adding an integral-action function on the frequency controller proposed by the authors in [25]. By considering two or more areas, this integral-action controller gives a proportional value in line with the accumulated ACE of the corresponding area, i.e., the linear combination of both tie-line power errors (Δ*Ptie*) and local frequency excursions (Δ*f*),

$$
\Delta T\_{set-point} = k \int (\Delta P\_{tie} + \beta \Delta f) dt,\tag{5}
$$

where *k*-parameter (constant of proportionality) is estimated depending on the expected frequency deviations (Δ*f*), the time interval remaining, the frequency excursion and the comfort level ranges allowed by the customers. Figure 6 shows the controller block diagram implemented in *Matlab-Simulink* environment.

**Figure 6.** Demand-side contribution to SFC: controller block diagram.

Assuming that controllable loads operate according to their thermostat settings, TCRLs' normal operation patterns are then modified by changing their temperature set-points, <sup>Δ</sup>*Tset*−*point*. Recovery to grid frequency deviations can be achieved by removing the residual steady-state error according to both AGC from the supply-side and Equation (5) from the demand-side. Communication is thus needed between the dispatch center and appliances when two or more areas are considered in the power system model. Recent contributions evaluate the impact on power system dynamic of errors in AGC systems [53]. Further discussion can be found in Section 3.2.

#### *3.2. Implementation of Integral-Action Controller*

For selecting the optimal value of *k*-parameter, three variables are considered. (i) *frequency deviation*: According to the European Standard EN–50160 [54], it is stated that the nominal value of grid frequency in most European countries, Asia and Africa is 50 Hz. It also provides the acceptable ranges for frequency variations: for interconnected supply systems under normal operating conditions, the averaged value of grid frequency measured over 10 s must be within a range of ±1% (49.5–50.5 Hz) for 99.5% of a week [55]. The proposed frequency control responses proportionally to the severity of the frequency excursion, and thus a proportional number of loads are called to switch-on/off depending on the frequency excursion to avoid undesirable over-frequency values. Subsequently, the number of loads is proportionally distributed taking into account that 100% of loads are required to support frequency control when the deviation is higher than 500 mHz and then grid frequency is kept within a range of ±1% (49.5–50.5 Hz), as previously commented. (ii) *duration of the frequency deviation*: The frequency deviation is considered to remain for at least 5 s, to be in line with the frequency excursion events emulated in the aforementioned contributions. (iii) *maximum allowable temperature variation for the thermostat set points*: Firstly, controllable loads are categorized into different groups according to their usage patterns and load profiles: *(Load-Group I)* including fridges/freezers; *(Load-Group II)* with air-conditioners/heat pumps; and *(Load-Group III)* corresponding to electric water heaters that have a grea<sup>t</sup> potential to store energy in advance [56]. In [57], load profiles of selected major household appliances in the U.S. are discussed in detail. When considering *Load-Group I*, and on the basis of technical specifications from different appliances manufacturers, the maximum allowable

temperature variation range for the fridge/freezer thermostat set point is set to +2 ◦C. For *Load-Group II*, air-conditioners/heat pumps, a variation range of ±3 ◦C is assumed acceptable for the customers, i.e., +3 ◦C for cooling mode and −3 ◦C for heating mode. Examples of cooling devices in demand response can be found in [58]. Finally, −6 ◦C temperature variation range is considered for *Load-Group III*. Hence, *k*-parameters are determined separately for each individual group of loads as follows:

$$\begin{pmatrix} k\_{\text{Load-Group II}} & = -0.8 \, ^\circ \text{C} (\text{Hz} \cdot \text{s})\_\prime \\ k\_{\text{Load-Group II}} & = -1.2 \, ^\circ \text{C} (\text{Hz} \cdot \text{s})\_\prime \\ k\_{\text{Load-Group III}} & = +5.0 \, ^\circ \text{C} (\text{Hz} \cdot \text{s})\_\prime \end{pmatrix} \tag{6}$$

The value of *k* is thus estimated taking into account tolerances of the various frequency excursion phenomena that may occur on the mains (almost 95% of cases included). However, given the potential variability of the frequency deviations, in both magnitude and duration, the output of the integral-action controller needs to be limited. For example, in the case of large generation outages, the frequency deviation occurs too rapidly for the frequency containment reserves to be effectively activated and prevent the frequency from reaching values below and above these tolerances. This situation would cause the controllers to order unacceptable changes for temperature set points. For this reason, upper and lower saturation limits and rising and falling slew rates are included in the implemented model. Additionally, to avoid a synchronized and massive response of the loads, and consequently undesired frequency oscillations, a linear density probability function is implemented on each load controller to decide in a distributed manner when an individual load participates (or not) in the demand-side response. Therefore, the percentage of controllable loads varies gradually depending on the severity of the frequency excursion, and thus emulating proportional response of classical turbine-generator governors subjected to under-frequency excursions. The load decision to be involved in frequency control (or not) is then determined individually by each controllable load attending to: (*i*) the frequency excursion along the time; (*ii*) the severity of such frequency deviations (the more severe is the frequency excursion, the higher is the number of controlled loads to participate in the demand response); and (*iii*) the own load thermal characteristics that are in line with the off-time period ranges allowed by the customers. The decentralized solution can thus be used for practical applications, in which it may need a huge number of appliance switches.

Figure 7 shows both thermal and power demand response of controllable loads subjected to an under-frequency event. Their power demand profiles are also included by considering the load controller proposed in Section 3.1. The thermal set-point value is modified according to Equation (5). In this case, we consider Δ *Ptie* = 0, and, thus, Δ *Tset*−*point* is estimated proportional to the accumulated deviation between the reference grid frequency and the measured value, according to Equation (5). Δ *Tset*−*point* represents the command sent by the individual load controller to the appliance after applying upper/lower and ramp limiters. As can be seen, within the first 25 s after the disturbance, the ramp slope is limited to a maximum rate of 0.08 ◦C per second. At the end of the simulation, the maximum temperature variation is limited to 2 ◦C. Moreover, this command has also been discretized in steps of 0.5 ◦C, i.e., Δ *Tset*−*point discretized*. The larger is the variation of the temperature set point, the lower is the duty-cycle of the load, thus decreasing power consumption and relieving power system reserves.

**Figure 7.** Frequency-response load controller: example of demand-side contribution to SFC.

#### **4. Simulation and Results**
