Virtual Inertia Control of Variable Speed Heat Pumps for the Provision of Frequency Support

: The growth in the integration of converter interfaced renewable energy has reduced the system inertia, which threatens system stability due to high rate of change of frequency (RoCoF) and frequency nadir issues unless steps are taken to mitigate it. There is a need to provide sufﬁcient fast frequency response to maintain adequate inertia in the system. This paper investigates the capabilities of a variable speed heat pump to provide an emulated inertial response. This paper presents a virtual synchronous machine control for a variable speed heat pump that provides support for grid frequency regulation over the inertial response time frame. A small-signal model with the transfer function of the variable speed heat pump is developed to analyse the effectiveness and feasibility of providing virtual inertia at the device and grid level, respectively. Furthermore, the small-signal model is validated using hardware in the loop simulation. Finally, the aggregated frequency response and virtual inertia contribution by a population of the heat pumps are evaluated and quantiﬁed in an urban distribution system.


Introduction
Europe has set an ambitious target to combat global warming by reducing the greenhouse gas emissions by at least 40% by 2030 and 80% by 2050 compared to the 1990 levels [1]. The emissions from the heating sector are about 600 Mt or 32% of Europe's total emissions [2]. Electrification of heating combined with the increased use of renewable generation plays a vital role in reducing these emissions. Electric heat pumps are seen as a leading technology to provide a path towards a low carbon heating sector [3]. At the same time, increased use of variable renewable generation from wind and PV reduces power system inertia, as these generators are connected through a power electronics interface and do not provide inertia. The reduction in system inertia due to the increased penetration of renewable generation has the potential to cause cascading failures and blackouts due to the high rate of change of frequency (RoCoF) unless actions are taken to mitigate it. This paper investigates the role that heat pump technology can play in providing an emulated inertial and fast frequency response in the context of low inertia power systems. In particular, it investigates whether a controlled response from heat pumps can provide a virtual inertia-type response to the system.
The concept of using the demand-side response as a cost-effective solution to provide greater system flexibility in the context of growing variable renewable penetration is well known. Domestic appliances such as electric water and space heaters, air conditioners, and heat pumps make use of inherent thermal inertia and can be built to aid in the system frequency regulation [4] while minimising the impact on end-users. For example, in the U.K., domestic heat pumps have become common, and it has been estimated that they could provide a frequency regulation capacity of 2 GW by 2030 [5].
In California, thermostatically controlled loads (TCLs) such as air conditioners, heat pumps, water heaters, and refrigerators have a frequency regulation capacity of 0.6 GW [6]. In many European countries, the electrification of heating combined with increased renewable generation is seen as an important pathway to energy system decarbonization, so it is likely that heat pumps as a potential resource for frequency regulation will grow in the future.
The participation of heats pumps in frequency regulation has been investigated in previous works. However, most of these works have focused on the response to centrally dispatched direct load control signals (DLC). For example, the work in [7] investigated such a response, but only including the dynamics of a heat pump induction machine in the model. The work in [8] developed a dynamic model of the VSHP, which included the dynamics of a variable speed drive (VSD) controlled induction machine, a heat pump cooling cycle, and an experimental commercial building. This was used to estimate the effect of DLC application to the VSHP on the building indoor air temperature. The work in [5] presented a thermodynamic model of the entire population of localised load controlled heat pumps connected to the U.K. power system, but it only reflected the frequency response based on an on/off control of the heat pump. The works in [9,10] presented a hardware demonstration of VSHP participating in frequency response through a centralised control signal.
Although the concept of providing frequency response from direct centralised control of heat pumps has been studied, the provision of emulated inertia from heat pumps has received little attention. In contrast to previous studies, in this paper, the provision of a controlled inertial response from VSHP is investigated. In the power system, inertia is an inherent property of the synchronous machine, present due to energy stored in its large rotating mass and its electromechanical coupling to the grid. In low inertia systems, the main concern for systems operators is to ensure adequate frequency response that is effective in the first hundreds of milliseconds after an event, which is most likely to impact the system dynamics [11]. Consequently, an emulated inertial response should act on fast localised measurements and not rely on centralised communications. The key to an inertial response is that it is proportional to the rate of change of frequency (RoCoF) and not the frequency deviation.
The provision of emulated inertia from power electronics interfaced storage and generation has gained attention in recent years, through for example virtual synchronous machine-type controls [12,13] of the power electronics converters. Providing an inertial response typically requires short-term energy storage, and it has been investigated for energy storage systems (ESS) such as batteries, supercapacitors, and flywheel storage in [14][15][16][17][18][19] and for photo-voltaic systems operating off their maximum power point [20]. For wind turbines, the stored kinetic energy in the turbine rotor can be used in combination with VSMcontrol, as shown in [21][22][23].
The contribution of this paper is to develop a virtual inertia control for the variable speed heat pump (VSHP) and, using this, estimate the level of inertial response that can be contributed by increasing use of VSHP in a distribution system. To do this, a detailed dynamic model of the VSHP is first developed including the dynamics of the VSD controlled induction machine, heat pump cooling cycle, thermal model of a residential building, phase-locked loop (PLL), and virtual inertia control. The resulting small-signal model based on transfer functions is validated with measurements from hardware. Using this model, the virtual inertia control is designed, and the resulting inertial response is characterised. These results show that the inertial response is limited by the uni-directional nature of the typical VSHP drive topology. Consequently, to achieve increased response, a modification of the drive approach is investigated, which can extract the machine stored kinetic energy to provide maximum support for grid frequency regulation. Finally, simulation studies are performed to evaluate and quantify the aggregated frequency response and virtual inertia capability of a population of such heat pumps in an urban distribution system.

Overview
A simplified schematic of the air source heat pump, which is widely used in the residential sector to exchange heat from the outside air to the building, is shown in Figure 1. The figure shows the VSHP, its electrical drive, and its controls, including virtual inertia control and indoor air temperature control. The variable speed electrical drive consists of a single-phase full-bridge diode rectifier and active power factor correction (PFC) feeding a voltage source inverter driving an induction machine. The single-phase full-bridge rectifier with active boost converter based PFC is used to control the DC voltage to the inverter. The voltage source inverter acts as the variable speed drive for the induction machine and is used to control the power consumed by the machine. The heat pump itself contains an evaporator, a compressor, which is driven by the induction machine, a condenser, and an expansion valve. The evaporator exchanges the heat between the refrigerant and the outside air, then the refrigerant temperature and pressure are further raised using the compressor. The heat collected from the evaporator and the compressor is transferred to the ambient air using the condenser. Finally, the expansion valve reduces refrigerant pressure.
The heat pump heating cycle model is used to determine the operational characteristics in both steady-state and transient operation [8]. An equivalent thermal parameter (ETP) model of the building being heated is included to determine the power setpoints and indoor temperature variations. The power setpoint for the VSHP is usually provided by a temperature controller, which attempts to control the room temperature to the user input setpoint.
The frequency response control is achieved by an outer closed-loop power control of the induction machine, which drives the heat pump compressor. For this work, it is assumed that the induction machine has a 2 kW rated power and 230 V rated voltage. The induction machine has an inner vector controller such that the rotor flux is always maintained constant in the steady-state [24]. An enhanced phase-locked loop (EPLL) [25,26] is used to determine the frequency of the grid voltage. The temperature controller provides the power setpoint P set . The virtual inertia controller provides the modifying power ∆P re f to provide frequency response. The combination of temperature setpoint P set and modifying power ∆P re f from the virtual inertia control provides a reference to the outer power controller of the induction machine drive, which then determines a torque reference for the vector control of the induction machine. The vector control outputs the modulation index, and the inverter is controlled using the sinusoidal pulse width modulation (SPWM) technique.
The thermal model of a residential building, EPLL, and virtual inertia control are implemented as explained in the section below. The goal of the modelling is to determine the dynamics of the variation of the heat pump power consumption, which can be achieved in response to a system frequency and RoCoF change. In order to achieve this, the various parts of the system and their small-signal models are described in detail in the section below.
The electrical drive topology of the VSHP shown in Figure 1, although commonly used, does not allow for energy to be returned to the grid supply due to the use of the diode bridge rectifier. However, the stored kinetic energy from induction machine loads can be extracted and returned to the source by employing regenerative braking with the VSD [27]. Figure 2 presents a variation in design of the VSHP drive system, which replaces the diode bridge and power factor correction circuit with a bi-directional full-bridge converter and passive LCLfilter. This allows bidirectional power exchange with the grid so that kinetic energy stored in the induction machine can if required be returned to the grid. In this approach, the full bridge converter is controlled so as to maintain the DC link voltage constant. The indicated DC bus voltage controller regulates the DC bus voltage V DC under variable load power P DC LPF , and the controller design is implemented as in [24]. The LCL filter is used to minimise the input current harmonics.

Thermal Model of a Residential Building
An equivalent thermal parameter (ETP) model of a residential building is used to provide the temperature feedback signal for the temperature controller. The ETP model of a residential building can be used to evaluate the thermal load capacities and the indoor air temperature variation due to the frequency responsive heating loads [28]. The evaluated thermal load capacities from the ETP model provides knowledge about how long the duration of the space heater power can be reduced to support frequency, without significantly changing the room temperature.
The ETP model is described by the following differential equations: whereQ hp is the heat flow into the building provided by the heat pump. C a , C room are the building indoor air and mass heat capacity, respectively. R 1 , R 2 are the thermal resistance between air to the ambient environment and air to the mass respectively. T in , T room , and T o are the indoor air, building mass, and ambient temperature, respectively. Typically, the C room is very large, and the temperature variation dT room /dt is negligible, Hence, it is assumed that T in =T room and the equivalent heat capacity C = C a + C room , and the equivalent model is reduced as shown in Equation (3). Table 1 presents the parameters for a residential building used in the equivalent electric circuit of an equivalent thermal parameter model, as shown in Figure 1. The time constants R 1 C, associated with the temperature controller, are much greater than the frequency response time constants, and hence, they do not interact, essentially meaning that reducing the space heater's power for a few seconds has a negligible impact on the room temperature. Table 1. Building model parameters [5]. 0.121

Dynamic Model of the Heat Pump
The dynamic model of the heat pump used follows the approach presented in [8], which developed a simplified dynamic model of a VSHP for real-time simulation studies. The characteristics of a VSHP (Equations (4) and (6)) in both steady-state and transient operation are extracted from the experimentally verified coupled nonlinear differential equation presented in [29,30], which relates the mass, momentum, and energy balances of all components in the refrigeration cycle [8]. From the data presented in [10], it is shown that the steady-state relationship between the compressor steady-state mechanical power P m0 and shaft speed ω r , ambient temperatures T o , and indoor air temperatures T in can be very reasonably approximated as a linear one as described in Equation (4). Since in this work, the main purpose of this relationship is to determine the steady state operating point for different indoor and outdoor temperatures, it is not expected that minor deviations from the assumed linear relationship would have any significant impact on the dynamic response of interest. The four parameters k ω , k c , k e , and k o f f set were determined using a multiple polynomial regression algorithm applied to 94 different combinations of input variables and output data obtained from the work conducted in [31].
Perturbing (4) and eliminating the small-signal cross products, DC terms yield the small-signal variation in steady-state mechanical power, P m0 , for a small variation in ω r , K in , and K o .
In terms of dynamics, the time constants associated with variations in temperature T in and T o are much larger than those associated with ω r [8]. Hence, for the purposes of fast inertial response, these can be neglected.
The transient change in power,P m , resulting from a change in mechanical speed is of interested here. In [8], it was shown that this dynamic relationship could be accurately approximated by a second-order system as shown in Equation (6). The work in [8] showed that the dynamic response approximated by the second-order system was virtually indistinguishable from that obtained from the nonlinear differential equations. The four coefficients n w1 , n w0 , d w1 , and d w0 are estimated using the polynomial regression algorithm [8].P m = n w1 s + n w0 The mechanical power variation of the compressor with respect to speed variation and mechanical torque variation can be expressed as: where Γ m0 and ω r0 are the initial operating torque and speed of the compressor, respectively. Substituting Equation (7) in Equation (6) yields the small-signal transfer function of the compressor.
Similarly, the heat flow rate, which is the input to the building model,Q hp =Q hp0 +Q hp , can be expressed in the same form as Equations (4) and (6) with different parameter values estimated using the polynomial regression algorithm.

Enhanced Phase-Locked Loop
A vital component of this control architecture is the PLL, which is required to calculate the angle and frequency of the grid voltage waveform with reasonable accuracy. The method for the calculation of the rate of change of frequency (RoCoF) is based on the derivative of the frequency measurement from the PLL. When measuring the RoCoF, noise is a critical issue. This noise comes from the grid voltage and also from the harmonics of the grid voltage, which can give rise to large fluctuations in measured RoCoF unless appropriate filtering is used. There is thus a trade-off between the accuracy and time response in the measurement of RoCoF. In general, increasing accuracy requires filtering, which also increases the time response.
Here, the enhanced phase-locked loop (EPLL) introduced in [25] is used, which has been shown to be more suitable for single-phase measurement in frequency varying and harmonics conditions. An adaptive notch filter is used in the EPLL to eliminate the double frequency errors and to estimate the voltage amplitude [26]. Low pass filters with time constants τ f LPF and τ RoCoF LPF are used in the EPLL measuring frequency and RoCOF, respectively. The schematic representation of the EPLL is shown in Figure 3.  Figure 3. Single-phase enhanced phase-locked loop (EPLL) system with a direct estimation of amplitude, phase, frequency, and the rate of change of frequency (RoCoF).

Small-Signal Model of EPLL
The EPLL differential equations and output of the EPLL phase detector were presented in [26], which are used to derive the small-signal model of EPLL.
Assuming sin(θ − φ) = θ − φ and representing Equation (9) in the Laplace domain representation yield:ω Substituting (11) in (10) yields: The actual phase angle to actual frequency in the Laplace domain representation is: Substituting Equation (13) in Equation (12) yields the closed loop transfer function of the EPLL without the added filter: The outputs from the EPLL are further low pass filtered, before being used to calculate the power reference to remove the disturbances caused by the harmonics. The low pass filter associated with the EPLL measurement for frequency and RoCoF are taken as first-order filters with transfer functions given as: where τ f LPF and τ RoCoF LPF are the time constant of the low pass filters used in the EPLL measuring frequency and RoCoF, respectively. The overall small-signal transfer functions of the EPLL including the low pass filter are given as: The EPLL control parameter k 1 represents the time constant of the adaptive notch filter, whereas k 2 , k 3 represent the natural frequency and damping ratio of the EPLL, respectively. Higher values for the EPLL parameter (k 2 , k 3 ) provide a faster response, but higher noise while measuring RoCoF. Lower gain values measure RoCoF without noise, but give a slower response. Hence, there is a trade-off between response time and noise attenuation while selecting the parameters. An appropriate set of EPLL parameters is selected using a trial and error approach in order to obtain a critically damped response, and the parameters obtained are as follows: k 1 = 40, k 2 = 1000, and k 3 = 50. The low pass filter time constants τ f LPF and τ RoCoF LPF are selected as 0.05 s and 0.1 s, respectively.

Virtual Inertia Control
The virtual inertia control is a strategy that emulates the dynamics of the synchronous machine, including inertia emulation, and damping power, when used to control converter interfacing distributed generation or energy storage [14]. Here, the control is only used to provide the frequency response by emulating the inertial response and the damping power of a synchronous machine. The power control block is shown in Figure 4. The controller modulates the power reference with a component that is proportional to frequency deviation (droop) and a component that is proportional to RoCoF (inertial). The output power reference of the controller ∆P re f is given by Equation (19) [32].
Here, M V I and D V I represent the virtual inertia and droop factor of the virtual inertia control, respectively, ω PLL LPF is a measurement of the angular frequency of the grid, and ω re f is the reference grid frequency. Note that the parameter D V I provides damping to the response during transients and also provides a steady-state response proportional to a frequency deviation. The saturation block limits the reference power to the rated power and maximum export power to the grid. Both components incorporate a dead-band, and if the frequency deviation and RoCoF are within this band, then no component is generated. The control receives its frequency and RoCoF measurement from the phase-locked loop (PLL) and then outputs the ∆P re f , which modifies the power setpoint from the temperature controller. Perturbing Equation (19) and eliminating small-signal cross-products and DC terms gives the virtual inertia control small-signal transfer function: whereω PLL LPF ω g and dω PLL LPF dt ω g are obtained from Equations (17) and (18), respectively.

Induction Machine Power Control
The induction machine is operated with an outer power controller and an inner torque vector control. The outer power controller receives its reference from the temperature controller and the frequency response controller. It utilises a PI control to output a reference torque for the inner torque controller. Here, we assume that the power reference change results only from the frequency controller as described by Equation (19). The inner torque control uses a standard vector control technique, the details of which can be found in [24]. Since the main goal of this section is to determine the small-signal response for power consumption, the torque control is only very briefly reviewed. The PI power controller outputs the machine electrical torque reference Γ e,re f according to the following relation: where K ppc , K ipc , P re f , P DC LPF are the power controller proportional and integrator gain, power reference, and measured power in the DC bus, respectively. Making use of a rotating reference frame synchronised with the rotor field, the q-axis component of stator current reference i sq,re f can be taken to be linearly related to the machine electrical torque reference: where σ r = ( L r L m − 1), L r , L m , and i mr are the rotor leakage factor, rotor inductance, magnetising inductance, and the magnetising current, and assuming that i mr , the machine magnetising current is held constant. Thus, torque can be controlled using the q-axis component of stator current. With a properly designed current controller, the q-axis stator current can be made to follow its reference with a first-order lag, so that: where τ i is the time constant of the current controller. The electrical torque Γ e is linearly related to the q-axis component of stator current i sq as: Perturbing Equations (21)- (24) and eliminating the small-signal cross products and DC terms yield:Γ e,re f = Combining Equations (25)-(28) yields: The speed variationω r of the induction machine in response to a change in electrical torqueΓ e and mechanical torqueΓ m is given as: where J and B are the inertia and friction coefficient of the induction machine, respectively. Rearranging Equation (30) yields:ω Substituting Equation (29) in Equation (31): Similarly, combining Equations (21)-(23), perturbing, and eliminating the small-signal cross products and DC terms yield:î Neglecting the losses in the AC/DC converter (single phase rectifier with the boost PFC), the power consumed by the induction machine is: where v sd and v sq are the dq-axis stator voltages and i sd and i sq are the dq-axis stator currents. Assuming the power loss is negligible or constant and a dq system aligned with the rotor field, i.e., with v sd set to zero, then the power consumed by the induction machine P DC can be expressed as: The q-axis stator voltage variation is given as [8]: v sq = (L a s + R a )î sq + i sd0 L sωr (36) where L a =L s − L 2 m L r and R a =R s + R r L S L r . In order to measure the DC power for feedback, a low pass filter is used to remove the disturbances while measuring the DC bus power. The measured power with the filter dynamics is given as:P L 2 m i mr , and A 3 = 3L s i sq0 i sd0 .Γ m ω r is the compressor transfer function presented in Equation (8). Figure 5 now shows the complete small-signal frequency domain model of the induction machine power control with the main equations used summarised in Table 2

Number Equation Description
(25)Γ e,re f = K ppc s+K ipc s (P re f −P DC LPF ) PI controller for power to torque reference

Overall Small-Signal Model
For this analysis, it is assumed that the power consumed by the PFC or DC bus voltage controller is negligible, and it is assumed thatP AC ≈P DC . Therefore, the overall closed loop transfer function is given as:P whereP DĈ P re f andP re f ω g are obtained from Equations (39) and (20), respectively. Table 3 summarises the initial operating point and parameter values used for controller design studies. From Equation (37), the input power variation of the VSD controlled VSHP (P DC LPF ) has a dependence on the dynamics of the stator currentî sq and the motor speedω r . With the correct design of the current controller,î sq can be made to change rapidly. The change inω r is limited by the compressor characteristics (inertia, friction, etc.). Theî sq term dominates the initial change so that the initial input power variation is similar for different compressor characteristics. Figure 6a presents the Bode plot for the open loop transfer function from the change inΓ e,re f to the resulting change in measured DC link powerP DC LPF for the operating conditions in Table 3. It is inferred from Figure 6a that the input power of the VSD controlled VSHP P DC LPF can be actively controlled to follow the power reference P re f with large PI controller gains. The cross-over frequency of the controlled system is chosen to be one decade less than the switching frequency (i.e., 1000 Hz) and a phase margin of 60 • , which yields K ppc and K ipc equal to 0.103 and 455.67, respectively. Figure 6b shows the resulting closed loop transfer function from a change inP re f to the resulting change in DC powerP DC LPF . It can be inferred that the closed loop system is stable. Figure 7 shows the step response of VSHP power consumption P AC to 1 W variation in P re f , which responds in 0.5 ms.

Figure 7.
Step response of VSHP power consumption P AC to 1 W variation in P re f .

Virtual Inertia Control Parameter Selection
In order to choose the virtual inertia controller parameters, M V I and D V I , some estimate of the quantity of response available is required. The quantity of demand response available from the VSHP obviously depends on its operating conditions. The speed setpoint for different ambient temperatures is determined from the thermal model of a residential building described in Section 2.2. Manufacturers recommend that the heat pump should not be turned off/on frequently and should operate at a minimum speed of (i.e., 1/3 of the rated speed). Hence, a constraint is introduced such that VSHP speed is always greater than or equal to its minimum speed. For this analysis, the VSHP minimum speed is set as 1/3 of its rated speed. The available demand response is constrained by speed variation ∆ω max < ω T o − ω min , where ω T o is an operating speed for a particular ambient temperature, ω min is the minimum speed of the VSHP, and ∆ω max is the maximum available speed variation. The maximum available demand response or the maximum power variation to maximum speed variation is determined using Equation (6).
The maximum released kinetic energy due to heat pump speed variation can be estimated as: The estimated inertia for heat pump with rated power 0.9 kW and nominal speed 1500 rpm is J = 0.0127 kgm 2 using the formula presented in [33]. Table 4 shows the calculated maximum available speed variation, demand response, and released kinetic energy for different ambient temperatures. These parameters can be used as a basis for setting the values of inertia and droop based on the worst case frequency variation scenario as described in the next section.

Selection Based on the Worst Case Scenario
A constraint for the selection of the control parameter based on the worst case scenario is introduced such that the controller provides its entire available power ∆P re f max plus the stored kinetic energy ∆KE max as estimated in the previous section, for the worst case scenario. The worst case scenario considered is a system with maximum RoCoF of 1 Hz/s and maximum frequency deviation of 1 Hz with all of the available kinetic energy released over 1 s. Substituting these constraints in Equations (19) yields: A fixed droop of 4% is used, i.e., for a 4% change in frequency; the power change is 100%, which yields: Substituting Equation (43) in Equation (42) yields the maximum available virtual inertia.
Based on these assumptions, Table 5 presents the selected control parameters for different outdoor temperatures.

Hardware Validation and Single Device Characterisation
Hardware validation and device characterisation were performed using a VSHP emulator consisting of a 2 kW wound rotor induction machine mechanically coupled to a dynamometer with variable torque generator, which supplied the heat pump compressor load torque. The power conversion system consisted of three LabVolt IGBT Chopper/Inverter 8857-1 modules with the control implemented using the hardware in the loop system OP5600 from OPAL-RT.

Single Device Characteristics
To determine the characteristics of the virtual inertia controlled VSHP, power hardware in the loop (PHIL) simulation was performed with the parameters summarised in Table 3 and for three different sets of virtual inertia control parameters (Case 1: M V I = 0 Ws 2 , D V I = 500 Ws; Case 2 M V I = 150 Ws 2 , D V I = 500 Ws; and Case 3 M V I = 350 Ws 2 , D V I = 500 Ws). Initially, a supply voltage of 230 V and frequency of 50 Hz were supplied to the system. At time t = 42 s, a ramp change in frequency −1 Hz/s was introduced until the frequency reached 49.5 Hz. These initial tests were done with the uni-directional electrical drive containing the diode rectifier as shown in Figure 1. Figure 8 shows the variation in the power demand and compressor speed of the VSHP's system for a change in grid frequency. Figure 8b shows that the variation in compressor speed was smooth. The effect of the increased virtual inertia setting can be clearly observed from Figure 8c,d. In Case 1, the virtual inertia controller setting was set as a simple droop; hence, the change in power ∆P DC LPF contained only a droop response that was proportional to the frequency deviation. In Case 2 and Case 3, the virtual inertia controller responded to the RoCoF of −1 Hz/s from simulation time t = 42 s to t = 42.5 s; hence, the change in power ∆P DC LPF was larger compared to Case 1. After simulation time t = 42.5 s, the RoCoF was 0 Hz/s; hence, the virtual inertial response was zero, and ∆P DC LPF contained only a droop response.
In Figure 8c, for a high virtual inertia setting, i.e., for Case 3, it can be observed that the DC bus power was negative, indicating that the induction machine P DC LPF regenerated back to the supply in this situation. However, the power transferred from the grid to the DC bus P AC was zero, as shown in Figure 8d, since the diode rectifier blocked the negative power being supplied to the grid. A simple replacement of the diode rectifier by the proposed bi-directional full-bridge converter as shown in Figure 2 would allow this extra inertial response to be transferred to the grid. Hence, for more effective virtual inertia support, we proposed this modification to replace the diode rectifier and PFC with a bidirectional full-bridge converter.

Modified VSHP System
Again, this system was set up in hardware as described previously, and the frequency of the VSHP's supply voltage was varied at 42 s and for six different sets of virtual inertia control parameters:  Figure 9 shows the variation in power demand of the VSHP's system for different cases. Figure 9b,c shows that the modified VSHP system had the capability to reduce demand quickly, and in this case, for large virtual inertia settings, the system could transfer the negative power to the grid. Simulation Time (s) (c) P AC Figure 9. Modified VSHP measured test result.

Validation of VSHP Small-Signal Model
The small-signal model of the VSHP system was later used for determining the aggregated response of a large population of VSHP systems in Section 4. In order to validate the small signal model of the VSHP system as expressed by Equation (40), the predicted output from the model in terms of power variation in response to a change in frequency was compared to measured results from the hardware system. The system was set up in hardware as described previously, and the frequency of the VSHP's supply voltage was varied at 42 s for four different sets of virtual inertia control parameters  Figure 10 shows the comparison between the response of the small signal model and hardware setup to a ramp change in frequency for the different cases. It can be seen from these results that the small signal model was a good representation of the hardware result, thus giving confidence that it could be used further in modelling studies.

Case Study: Aggregated Response from a Population of Heat Pumps
In order to understand the aggregate response and quantify the maximum attainable virtual inertia from a large population of VSHP embedded in the distribution grid, electromagnetic transient (EMT) simulations were performed on the real-time simulation platform [34]. For the simulations, each VSHP was represented by the validated transfer function in Equation (40) using the virtual inertia control settings derived using the worst case settings described earlier.

Simulation Results in the Test Distribution System
The distribution grid model used for real-time simulation was based on an urban LV distribution grid in the U.K. [35], a diagram of which is shown in Figure 11. The network contained three feeders with 330 houses, 6350 nodes, three phase, and four wire configuration distribution cables. The system frequency of the network was 50 Hz, and the voltage was fed to the network through an MV (13.8 kV) to an LV (400 V) delta-star transformer, as shown in Figure 11. The challenges and methodology used in modelling this network in the real-time platform were further described in [34]. The network was divided into five groups for parallel processing using the real-time simulator, and the colours illustrate how the network was split into the five groups. The loads in the distribution system were modelled as per the recommendations by the IEEE Task Force on Load Representation for Dynamic Performance [36]. The load model consisted of ZIPterms plus two voltage-/frequency-dependent terms as shown in the following equations: where P 0 , Q 0 ,, and V 0 are the nominal active, reactive load powers, and load voltage, respectively. Table 6 provides the load data description and the values used for the load modelling. Per unit of active or reactive load that is voltage and frequency sensitive (Term 1) 0.100, 0.384 K p2 , K q2 Per unit of active or reactive load that is voltage and frequency sensitive (Term 2) 0.400, 0.572 n pv1 , n qv1 Voltage sensitivity exponent (Term 1) 1, 3 n pv2 , n qv2 Voltage sensitivity exponent (Term 2) 0.100, 0.500 n p f 1 , n q f 1 Frequency sensitivity (Term 1) 1, −2.800 n p f 2 , n q f 2 Frequency sensitivity (Term 2) 1.900, 1.200 In order to observe the effect of the virtual inertia controlled VSHPs, the distribution system was supplied by a synchronous generator with a steam turbine governor model, and a transient increase in load was implemented to perturb the frequency. The parameters of the synchronous generator model and load damping are presented in Table 7. Simulations were performed for a different number of houses with the droop and virtual inertia controller installations and different outdoor temperatures. P nom gen is the nominal power of the synchronous machine. P base is the power base for per unit calculations. P load is the power consumed by the distribution network. D is the damping provided by the distribution network. M is the inertia of the synchronous machine. τ t is the time constant of the synchronous machine turbine. τ g is the time constant of the synchronous machine governor.

Droop vs. Virtual Inertia Controller
The time-domain simulations were performed on the distribution test grid with 300 installations of VSHP with the droop and virtual inertia controllers providing frequency response for an operating point at T o = −10 • C. At 7 s, a bulk three-phase load of 112.500 kW was connected to the network to produce a frequency perturbation. Figure 12a shows the frequency variation of the network for the 300 houses with the droop and virtual inertia controller installations. Note here that for "droop control", the entire available power was provided as a droop response (M V I = 0 Ws 2 D V I = 1022 Ws), whereas for "virtual inertia control", it was split between inertial and droop response. Clearly, the frequency control of the VSHP loads provided an improved frequency response compared to the case where no frequency control was implemented. Figure 12b shows the response from an individual heat pump with the droop and virtual inertia controller installation. It can be seen that the droop controller provided a better reduction in frequency nadir, whereas the virtual inertia controller provided a better reduction in RoCoF. There is potential scope to optimise the relative contribution of droop and inertial responses, depending on system requirements at any time. Dynamic tuning of the virtual inertia control parameter could be implemented to update settings so as to obtain a better trade-off between reduction in RoCoF and frequency nadir.

Increased Virtual Inertia Controller Installations
The time domain simulations were performed on the distribution test grid with different numbers of houses with virtual inertia controlled VSHP providing frequency response for an operating point at T o = −10 • C. At 7 s, a bulk three-phase load of 112.500 kW was connected to the network to produce a frequency perturbation. Figure 13 shows the frequency variation of the network for the different number of houses with the virtual inertia controlled VSHP. The RoCoF was calculated 0.300 s after the disturbance occurs, over a 1 ms time window, and the equivalent aggregated distribution system inertia M sys in pu was calculated from Equation (49) [37].
where ∆P is the variation in power and RoCoF (t=0.300 s) is the RoCoF measured 0.300 s after the disturbance. The calculated total system inertia for the system without virtual inertia controller installations and with 100, 200, and 300 virtual inertia controller installations was 5.200 s, 8.320 s, 11.910 s, and 14.920 s, respectively. Hence, the added virtual inertia from 100, 200, and 300 virtual inertia controller installations was 3.120 s, 6.710 s, and 9.720 s, respectively. Figure 13b shows the response from an individual heat pump for different virtual inertia controller installation. It can be seen that increasing the number of installations reduced the response from individual heat pumps, as they saw less frequency deviation and RoCoF.

Different Outdoor Temperatures
As indicted previously, the available response was dependent on outdoor temperature, and to investigate this, the time domain simulations were performed on the distribution test grid with 300 houses with virtual inertia controlled VSHP providing the frequency response for different operating points at T o = 5 • C, 0 • C, −5 • C, and −10 • C. At 7 s, a bulk three-phase load of 1112.500 kW was connected to the network to produce a frequency perturbation. Figure 14 shows the frequency variation of the network for different outdoor temperatures.
The calculated system inertia for the system without virtual inertia controller installations and with 300 virtual inertia controller installations at operating point T o = 5 • C, 0 • C, − 5 • C, and −10 • C was 5.2 s, 7.280 s, 9.710 s, 12.260 s, and 14.920 s, respectively. Hence, the added virtual inertia from 300 virtual inertia controller installations at operating point T o = 5 • C, 0 • C, −5 • C, and −10 • C was 2.080 s, 4.510 s, 7.060 s, and 9.720 s, respectively.  Figure 14b shows the response from an individual heat pump for different outdoor temperatures. It can be seen that as for lower temperatures, the response from individual heat pumps was higher since the virtual inertia controller setting was set at a higher value and vice versa.
Clearly, the ability of VSHPs to deliver frequency support is not constant, but dependant on consumer behaviour in how they use their heat pumps and is also dependant on temperature. This makes the frequency support resource variable, but given a large population of such VSHPs, the available resource at any given time and for any conditions will be reasonably predictable. The availability and use of the resource would therefore need to be planned and coordinated with other sources of flexibility as part of a wider system operation strategy.

Conclusions
The paper presented the dynamic model of a variable speed heat pump providing frequency response with virtual inertia control. The hardware in the loop simulation studies showed that the heat pump could quickly change its power consumption and also provide power during regenerative braking. The paper presented a modification in variable speed heat pump design to provide power from the heat pump to the grid during regenerative braking; hence providing maximum frequency response during the inertial response time frame. A small signal model was developed and validated with hardware in the loop simulation results. The simulation studies performed on a test distribution grid verified that the virtual inertia controlled variable speed heat pump reduced the frequency nadir and RoCoF and provided significant virtual inertia to the system. The level of virtual inertia obtained from the heat pumps was quite significant, but obviously depended on ambient temperatures and the amount of penetration of these demand response controllers. To select the virtual inertia controller settings, there was a trade-off between virtual inertia gain, which helped in the reduction of RoCoF, and droop gain, which helped in the reduction of frequency nadir.

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