The role of wind turbines (WT) is to convert wind energy, an environmentally friendly energy source, into a form of electrical energy. It is one of the most promising renewable energy sources and its global installed capacity is increasing year by year (Ren et al., 2021). According to the direction of their rotation axis, WTs can be classified into horizontal axis wind turbines (HAWT) and vertical axis wind turbines (VAWT). HAWT extracts wind energy on the horizontal axis and the blades rotate due to the lift provided by the aerodynamic force from wind. Because of it has higher efficiency than VAWT, the HAWT is more popular in research. However, HAWT needs to always point to the wind direction to work effectively. For unpredictable wind directions, the HAWT needs additional mechanisms to ensure that the blades are always facing the wind direction. The VAWT has its blades perpendicular to the ground and rotating around the vertical axis. VAWT can receive wind power from any direction. Compared with HAWT, VAWT has several advantages, including no need to face to the wind flow, the ability of generating electricity at low wind speed, no need for additional control on yaw and pitch, less maintenance costs, and less noise (Eriksson et al., 2008; Johari et al., 2018). However, VAWT has much lower efficiency at high wind speeds than HAWT, its dynamic stability is relatively poor due to its light weight, and it is vulnerable to backtracking wind. For the above reasons, VAWT is often installed in slow and turbulent wind environments, such as roof top, for low-power generation (Bhutta et al., 2012), while HAWT is often used for high-power generation in isolated and remote areas, such as on-shore or off-shore wind farms (Johari et al., 2018). In recent years, some researchers have combined the two types of WT together to improve its operational capabilities (Govind, 2017).

Although HAWT and VAWT have lots of difference in structure, control strategy and power efficiency, the methods they use to capture power from wind energy are the same, and thus can be expressed by the same equation. The available power output of a WT is proportional to the cube of the wind speed. If the wind speed is too low, the output power of the WT can be ignored. When the wind speed increases, the output power of the WT increases rapidly. Therefore, in order to ensure the normal operation of a WT, three threshold wind speeds usually need to be considered, including the cut-in speed, the rated speed and the cut-out speed (Ren et al., 2016). When the wind speed is lower than the cut-in speed, the WT stops to avoid unnecessary mechanical wear and energy consumption; when the wind speed is higher than the cut-in speed and lower than the rated speed, the WT controls its rotating speed to track the maximum power point referring to the current wind speed; when the wind speed is higher than the rated speed and lower than the cut-out speed, the WT adjusts its pitch angle to limit the captured wind power not to exceed the limit of the generator; when the wind speed exceeds the cut-out speed, the WT shuts down to protect the rotor and the generator from the damage of strong wind.

The power output of WT captured from wind is presented as P W T ( h ) = { 0 V w s ( h ) < V c u t − i n o r V w s ( h ) ≥ V c u t − o u t 1 2 ρ A r C p V w s 3 ( h ) V c u t − i n ≤ V w s ( h ) < V r a t e P W T , r a t e V r a t e ≤ V w s ( h ) < V c u t − o u t (1) where ρ is the air density, A r is the swept area of WT that depends on the turbine radius or length of blades, V ws ( h ) is the hourly mean wind speed, and C p is the power conversion coefficient that is determined by the aerodynamic structure of the WT. V cut − in , V rate , and V cut − out denote the cut-in wind speed, rated wind speed and cut-out wind speed.

When designing a renewable energy system, it is necessary to determine the suitable capacity of the WT according to the rated power of the WT and the local wind resources. During the operation, it is necessary to consider whether the current wind speed has reached the cut-in speed, rated speed and cut-out speed, and calculate the power that the WT can generate according to the wind speed.

Solar energy is another well-known clean energy. Since most parts of the world have sufficient solar irradiance, solar power generation technology has attracted a large number of researchers. Photovoltaic (PV) is a device that uses the photoelectric effect to directly and conveniently convert solar energy into electricity (Liu et al., 2017). In recent years, with the development of PV technology and the rapid decline of PV costs, large-scale grid-connected solar PV power stations have been built all over the world (Mandal et al., 2012).

However, PV systems are highly dependent on solar radiation and weather characteristics. Due to the variability of solar radiation and ambient temperature, the power output of the PV system is uncertain and random (Mandal et al., 2012). The uncertainty is from the alternation of day and night, effect of seasonal changes, and random cloud movement (Jamal et al., 2017; Jamal et al., 2019). Compared with the distributed PV system, the power output of the relatively concentrated and single format PV power generation system in a small area can change very quickly. An isolated community can be completely or partially obscured by fast-moving clouds in the span of a few seconds to a few minutes (Schmidt et al., 2017). The movement of the cloud will affect the performance of the PV system, so prediction must be made to avoid undesirable technical problems and cost losses (Jamal et al., 2019). Predicting the power output of a PV system is a challenging task.

For a photovoltaic array composed of PV panels, the maximum power output can be expressed as (Tao et al., 2010) P P V ( h ) = η p v A p v S T ( h ) ⋅ [ 1 − 0.005 ( T a i r ( h ) − 25 ) ] (2) where η pv is the photoelectric conversion efficiency of the PV array, A pv is the total area of the PV array, S T is the solar radiation incident on the panel, and T air is the air temperature of the ambient in Celsius degree.

The power output of PV system mainly depends on the amount of solar radiation and the operating temperature of the installation location. In addition, as the equipment ages, the photoelectric conversion efficiency and power conversion efficiency will gradually decrease. But their changes in the system life cycle are very small, so they can be ignored. Therefore, as long as the solar radiation and panel temperature are known, the power output of the PV system can be predicted.

As the main renewable energy sources, solar and wind energy are intermittent and fluctuating weather resources. One of the possible solutions is to use a hybrid renewable energy system to integrate various renewable energy sources in the best combination (Yang et al., 2018a). A good example is the complementarity of solar and wind energy (Prasad et al., 2017). However, this approach still has large uncertainties. If the distributed resources cannot provide sufficient flexibility in power generation, network operators will have to limit the penetration rate of renewable energy to a certain level, which will adversely affect the use of renewable energy systems (Jamal et al., 2019). For example, the power generation of the PV system usually exceeds the residential electricity load during the day but cannot meet the demand in the late afternoon and evening when the residential load tends to increase (O'Shaughnessy et al., 2018). Therefore, in order to highly integrate variable wind and solar energy and reduce the impact of weather changes on renewable power generation, energy storage systems have been widely accepted as one of the potential solutions (e Silva and Hendrick, 2017; Babacan et al., 2017; Zhang et al., 2018; Ren et al., 2020). The main function of energy storage system is to compensate the imbalance between power generation and demand, transfer energy from periods of high power generation to periods of low power generation, so that the system can maintain full functionality under various operating conditions (Dragicevic et al., 2014).

ESS includes a variety of technologies, such as pumped water storage, compressed air energy storage, hydrogen storage, flywheels, supercapacitors, and batteries (da Silva Lima et al., 2021). The batteries have higher energy density than supercapacitors, higher charging rate than pumped water storage and less power loss than flywheels. Therefore, the battery storage system (BSS) is the most suitable storage system to store the electricity energy for such localised small-scale applications due to its advantage of higher charging rate, higher energy density and shorter response time. The response time of BSS is between milliseconds and seconds, while that of PHS and CAES are from a few seconds to minutes (Castillo and Gayme, 2014). Due to their durability, low maintenance and low social environmental impact, it is expected that the development and use of BSS will increase significantly in the next decades (da Silva Lima et al., 2021). Among batteries, the most used are lithium-ion batteries (LIB) and vanadium redox flow batteries (VRB). The advantages of LIB are its high energy density, high efficiency, long lifecycle, and more environment friendly characteristics. However, due to these advantages of LIB, they are also widely used in a variety of other applications, including small electronic products and electric vehicles. VRB is a good substitute for LIB because of its safety, long service life, better scalability and high recyclability (da Silva Lima et al., 2021). However, the disadvantage of VRB is its low energy density (Castillo and Gayme, 2014). Therefore, a VRB requires more floor space than a LIB at the same capacity.

In the charging and discharging behaviour of a battery, the state of charge (SOC) is often used to reflect the ratio of the battery’s remaining capacity to its rated capacity (Yang et al., 2018b). To protect battery lifecycle, it is necessary to limit the SOC within a certain range. During battery charging, the SOC can be given as (Zhang et al., 2018) S O C B S S ( h ) = S O C B S S ( h − 1 ) ⋅ ( 1 − σ ) + P B S S ( h ) ⋅ Δ h E B S S , c a p ⋅ η c h (3) and during battery discharging, the SOC can be given as S O C B S S ( h ) = S O C B S S ( h − 1 ) ⋅ ( 1 − σ ) + P B S S ( h ) ⋅ Δ h E B S S , c a p / η d c h (4) where σ denotes the self-discharge rate, P BSS denotes the power of BSS generated or absorbed, E BSS , cap denotes the energy capacity of BSS, η ch and η dch denote the battery charge and discharge efficiency, respectively. Δ h represents the time step, in this article, the time step is calculated in hours.

Heat pumps are a group of energy systems that can extract heat energy from lower temperature sources such as ambient air, soil, water in lakes and rivers, and then upgrade and deliver it at a higher temperature for heating applications (Watzlaf and Ackman, 2006). If the heat source is outside air, the system is called an air source heat pump (ASHP). If the heat source is the ground, soil or groundwater, it is called a geothermal heat pump (GHP). GHP is one of the fastest growing renewable energy sources in the world (Arat and Arslan, 2017). Compared with traditional heating fuels, such as natural gas, heating oil and propane, the effectiveness of operational cost of GHPs is directly related to the cost of electricity that drives the HP (Self et al., 2013). But the installation costs of GHPs can be an order of magnitude higher compared to a gas boiler.

Heat energy delivered by the HP normally exceeds its electricity consumption. Coefficient of performance (COP) defined as the ratio of heat output to input power is used to describe the energy performance of heat pumps (Casasso and Sethi, 2014). The COP of heat pumps can reach 3–5 or more depending on temperature conditions (Renaldi et al., 2017), which means the heat pump can provide 3–5 kW of heat by consuming 1 kW of electricity. Furthermore, the smaller temperature difference between the heat source and the radiator, the greater the COP of the HP system. In cold weather, the COP of an ASHP is usually about 2 or lower. Therefore, GHP have higher efficiency than ASHP because the underground temperature is closer to the required indoor conditions than the outdoor air temperature in winter and the GHP usually use water which reduces the size of the evaporator heat exchanger due to the higher energy density compared to air. However, as GHP requires higher installation costs and larger spaces, ASHP is more suitable for heating smaller spaces.

In HP applications, its COP can be affected by different variables, such as external temperature, water supply temperature, inlet water temperature and load factor. Simplifications can be made to reduce this complexity. The COP can be calculated as (Ozgener and Hepbasli, 2007): C O P H P = Q ˙ H P P H P , c o m p (5) where Q ˙ HP is the rate of heat transfer in condenser of HP, P HP , comp is the input power to the compressor.

The rate of heat transfer from the condenser is obtained as Q ˙ H P = m ˙ r e f ( T H P − T r e t u r n ) (6)

In this study, the air source heat pump is considered. For simplification, the COP of the HP is modelled as a function of temperature lift, or the temperature lift is known as the difference between the supply water temperature T wat and the external air temperature T temp (Renaldi et al., 2017). C O P H P = a ⋅ ( T w a t − T t e m p ) + b (7) where a and b are the parameters from the linear regression fits of manufacturer’s data (Renaldi et al., 2017). However, the linear regression is just an initial fitting for the easy implementation of modelling. A nonlinear model can be more accurate. Since this paper aims to provide a feasible solution for investment of heat pump rather than operational control, the linear regression model is sufficient for estimating the COP of HP.

The outdoor temperature changes between summer and winter during the year can produce huge changes of heat demand (Xu et al., 2018). In this case, a lot of heat will be wasted in the summer when the heating demand is minimal, and there will be insufficient heat in the winter. In addition, the mismatch between energy demand and supply is often observed. Among several storage technologies, thermal energy storage (TES) seems to be one of the most promising technologies that can compensate for the intermittent heat generation of new energy sources. The advantages of using TES in energy systems include improved overall efficiency, better reliability, bring better economy, reduced investment, and operating costs, and less pollution and CO 2 emission to the environment (Sarbu and Sebarchievici, 2018).

In recent decades, the TES technology has been widely studied. Considering the required storage time, the TES system can be divided into short-term storage and long-term storage. The short-term storage TES has a charging and discharging cycle of hours to several days, which is called diurnal thermal energy storage. The long-term storage TES can store heat for several months and is called seasonal TES (Dahash et al., 2019; Renaldi and Friedrich, 2019; Maximov et al., 2021). The principle of seasonal TES is to store heat energy through heating equipment in summer, and then discharge it for space heating in winter (Kubiński and Szabłowski, 2020). The use of TES can significantly increase the flexibility and self-consumption of renewable energy for end users (Waser et al., 2018). Therefore, this concept makes a significant contribution to the efficient use of renewable energy in district heating systems and the decarbonization of the building sector.

The most traditional TES design is provided in the form of a water tank. When there is excess thermal energy on the supply side, these tanks are charged (heated up), and when the demand for thermal energy exceeds the supply, these tanks are discharged (Abdelsalam et al., 2020). In the thermal system, the charging and discharging process of the storage tank can be carried out by direct heat exchange. The energy stored in the TES at time h is expressed as (Renaldi et al., 2017) Q T E S , s t o ( h ) = Q T E S , s t o ( h − 1 ) + [ Q ˙ T E S , c h ( h ) − Q ˙ T E S , d c h ( h ) − Q ˙ l o s s ] ⋅ Δ h (8) where Q TES , sto is the energy stored in TES, Q ˙ TES , ch and Q ˙ TES , dch denote the charge the discharge power, Q ˙ loss denotes the standing losses, which normally can be found from manufacturer’s datasheet.

Among the performance indicators used to evaluate TES, the most common is TES efficiency, which is the ratio of the total heat energy recovered from the heat accumulator at the discharge temperature to the total heat input at the charging temperature as (Dahash et al., 2019) η T E S = Q T E S , d c h Q T E S , c h (9)

The purpose of this study is to find the most effective investment plan to reduce CO 2 emissions in building heating under a given financial budget. Therefore, the optimisation problem is transferred to minimising the total CO 2 emission in the period between 2020 and 2050. The objective function is expressed as min ∑ y i = 2020 2050 S C O 2 ( y i ) (10) where S CO 2 ( y i ) denotes the annual CO 2 emission in year y i , which is expressed as the product of specific CO 2 emission factor for nature gas, g gas , and heat energy produced by gas boiler, Q gas .

In this paper, the RHS uses the electricity from local renewable energy. Thus, it is assumed that electricity power generation has no CO 2 emissions in the RHS, and the CO 2 emissions are entirely generated by using the natural gas in the traditional heating system. This assumption can intuitively indicate the reduction in the use of traditional heating system and nature gas as well as the total CO 2 emission during operation. The factor of S CO 2 ( y i ) can be used as an indicator to quantify the decarbonization performance, which is defined as follows. S C O 2 ( y i ) = g g a s ⋅ Q g a s ( y i ) (11)

This equation only considers the CO 2 generated from gas and assumed that the electricity generation has no CO 2 emission in the RHS. This assumption is simple and directly react the CO₂ in operation. And use this factor as an index to justify the performance of decarbonisation.

And Q gas can be presented as Q g a s ( y i ) = ∑ h = 1 8760 [ Q ˙ h e a t ( h ) − Q ˙ H P ( h ) + Q ˙ T E S , c h ( h ) − Q ˙ T E S , d c h ( h ) ] (12) where Q ˙ heat ( h ) is the hourly heat demand of buildings and Q ˙ HP ( h ) is the hourly heat produced by heat pump, h is the time factor of hours in a whole year, which in total are 365 days times 24 h. The Q ˙ HP ( h ) is determined by both the electricity generated from PV, WT and BSS and its COP as Q ˙ H P ( h ) = C O P H P ( T a i r ( h ) ) ⋅ ( P P V ( h ) + P W T ( h ) − P B S S ( h ) ) (13)

In practice, the output of power generator and heating equipment cannot exceed its rated value, and the energy stored by storage equipment cannot exceed its specified capacity. In addition, the annual installed capacity cannot be negative, which means that the total capacity of each component cannot be lower than the previous year. The constraints of all types of components are given as 0 ≤ P P V ≤ P P V , r a t e ( y i ) (14) 0 ≤ P W T ≤ P W T , r a t e ( y i ) (15) 0 ≤ Q ˙ H P ≤ Q ˙ H P , r a t e ( y i ) (16) S O C min < S O C B S S ( h ) < S O C max (17) 0 ≤ Q T E S , s t o ≤ Q T E S , c a p ( y i ) (18) P P V , r a t e ( y i ) ≥ P P V , r a t e ( y i − 1 ) (19) P W T , r a t e ( y i ) ≥ P W T , r a t e ( y i − 1 ) (20) Q ˙ H P , r a t e ( y i ) ≥ Q ˙ H P , r a t e ( y i − 1 ) (21) E B S S , c a p ( y i ) ≥ E B S S , c a p ( y i − 1 ) (22) Q T E S , c a p ( y i ) ≥ Q T E S , c a p ( y i − 1 ) (23)

In addition to the above technical constraints, another constraint is the financial budget. In this study, it is assumed that the investment plan has an Annual Investment Cost (AIC) limit. The investment cost in each year can be used to increase the capacity of any component, but the total cost cannot exceed AIC as A I C ≥ Δ P P V , r a t e ( y i ) ⋅ U P V ( y i ) + Δ P W T , r a t e ( y i ) ⋅ U W T ( y i ) + Δ Q ˙ H P , r a t e ( y i ) ⋅ U H P ( y i ) + Δ E B S S , c a p ( y i ) ⋅ U B S S ( y i ) + Δ Q T E S , c a p ( y i ) ⋅ U T E S ( y i ) (24) where Δ P P V , c a p , Δ P W T , c a p , Δ Q ˙ H P , r a t e , Δ E B S S , c a p , and Δ Q T E S , c a p denote the increased rated power of PV, WT and HP and capacity of BSS and TES of each year. U PV ( y i ) , U WT ( y i ) , U HP ( y i ) , U BSS ( y i ) , and U TES ( y i ) denote the unit cost of each component in year y i .

The forecasted component unit price in future years is obtained based on different financial assumptions in recently published papers (Keiner et al., 2019; Yue et al., 2020; Rinaldi et al., 2021). In this study, the financial forecast for the unit cost of components for 2020–2050 is in British pounds, as shown in Figure 2. The prediction shows that the unit prices of PV, WT, HP, BSS, and TES will continue to fall in the next few decades, while the unit price of natural gas and CO₂ will increase rapidly in the coming decades.

From the literature in the past decades, there are four main optimisation approaches that are commonly used, including the direct search, calculus-based optimisation, genetic algorithm (GA), and particle swarm optimisation (PSO). Direct search is a straightforward optimisation method that does not consider time delay and derivatives. Thus, the direct search method can be applied in optimising many nonlinear functions that have less dependence on its derivatives (Kolda et al., 2003). The most commonly used Calculus-based Optimisation method is to set the gradient of the objective function to zero (Kheiri, 2018). Another calculus-based optimisation method is Newton’s method, which is similar to the Steepest Descent Method that uses an iterative process from an initial guessed starting point to finally converge to the optimum point. The GA is a population-based algorithm to search for the global-optimum solution. The iterative process of GA will converge to better solutions based on the breeding of the parents with higher performance (Elbeltagi et al., 2005; Kheiri, 2018). However, the GA has shortages that it requires a large amount of non-optimal data as the requirement of global search (Kheiri, 2018). The PSO is a population-based metaheuristic algorithm that attempts to find the global optimal solution of the optimisation problem by simulating the social gathering behaviour of animals (Eberhart and Kennedy, 1995; Kheiri, 2018). This paper chooses the PSO method to find the optimal sizing of RHS due to its strength in finding the global optimal solution and easy to implement.

PSO method uses a large number of the swarm to search the optimum point globally and share the information among all swarm for the next search step. In the PSO algorithm, each feasible solution is called a particle, which is specified by a vector containing problem variables (Maleki et al., 2016). Particles can remember their previous optimal position and share it with others (Shi, 2001). The motion of each particle is composed of two randomly weighted influences and an initial random velocity (Shi, 2001). The PSO algorithm simulates the sociality of particle, combines the best position of the entire swarm and the trend of its own movement, effectively avoids the particles and swarm falling into the local optimal solutions (Baniasadi et al., 2020).

In the PSO algorithm, the state of each particle is represented by its position x j and velocity v j . The velocity update equation of the PSO algorithm has three key parameters, including the inertia constant w , the acceleration constant c 1 that controls the direction of the particles toward its best position in history, and the acceleration constant c 2 that attracts the particles to the best position of the swarm. The formula for updating the velocity and position of each particle in the space can be expressed as (Maleki et al., 2016; Mekontso et al., 2019) v j k + 1 = w k ⋅ v j k + c 1 ⋅ r 1 k ⋅ ( p b e s t , j k − x j k ) + c 2 ⋅ r 2 k ⋅ ( g b e s t k − x j k ) (25) x j k + 1 = x j k + v j k + 1 (26) where v j k and x j k are the velocity and position of the j th particle in the k th iteration, p best , j k is the best position achieved by the j th particle and g best k is the best position of the swarm, r 1 and r 2 are the random factors between 0 and 1.

In the previous simulation case, it can be found that the use of HP will bring additional power load to the grid. This shows that although the use of HP can greatly reduce CO 2 emissions, there is a potential risk of overloading or even fault occurring in the power grid. Therefore, the most promising alternative is to use local renewable resources to generate electricity to drive HP to heat buildings. This will reduce the impact of the large amount of electrical energy required by HP on the grid. In this case, it is assumed that the heating system is independent of the grid, and only uses renewable energy for power supply to the HP. The heating system combines the RHS based on local renewable resources and the traditional heating system using gas boiler. When renewable resources can generate enough power, the heat demand is completely supplied by HP, and the unconsumed power is stored in the BSS. When renewable resources driven HP cannot provide enough heat for the building, traditional gas boilers will be used to compensate for the remaining heat demand by consuming natural gas.

Figure 4 shows an example of monthly heat production when renewable generation capacity is 10 MW, where the rated power of PV and WT are both 5 MW. Since this case aims to consider the role of the BSS rather than its actual size, the capacity of the BSS is assumed to be large enough to store all unused energy from renewable energy sources. The blue part, the green part, and the yellow part in the bar chart represent the heat generated by the HP driven by power from PV, WT, and BSS, respectively. It can be seen from the results that, due to the large changes in the illumination time and the angle of solar radiation in high latitude areas such as Scotland, PV power generation is higher in summer and lower in winter. WT has the opposite trend as it will produce more power in winter and less power in summer. From the results of each month, the heat demand of the building in July, August, and September is less, so the heat demand can be fully supplied by HP. In addition, due to the low heat demand, the unconsumed power generated by renewable resources is stored in the BSS. Therefore, in October and November, most of the electricity is provided by the BSS that has stored electricity in the past few months. From December to April, the heat demand is high, but the power from renewable energy and BSS are low, so a large amount of heat is provided by gas boilers for building heating.

It can be seen from the above results that the natural gas consumption of gas boilers is related to the installed capacity of renewable energy. Therefore, the relationship between natural gas consumption and renewable energy installed capacity is shown in Figure 5A. Without the use of BSS, the gas consumption decreases slowly as the installed capacity of PV or WT increases, as shown by the dashed lines. Under the same installed capacity, WT-driven HP performs better than PV-driven HP in reduction of gas consumption. The reason is that WT’s power generation is large in winter and small in summer, which is more in line with building heat demand. However, PV power generation has a small amount of power generation in winter, so the power generation of PV in winter is not enough to support HP consumption. However, after using BSS with installed capacity in Figure 5B, the time distribution of renewable resources has less impact, and the total power generation is more important. As the installed capacity of renewable energy increases, the consumption of natural gas decreases faster since unconsumed electricity is stored in the BSS and used to assist the renewable energy generation. With the same installed capacity, PV can generate more electricity than WT, and thus results in less gas consumption. In this case, HP driven by PV can reduce gas consumption faster than HP driven by WT. However, due to the time mismatch between heat demand and PV power generation, the BSS capacity required for PV is much higher than that of WT. The required capacity of BSS corresponding to different renewable energy installation capacity is shown in Figure 5B. When the PV installed capacity reaches 15 MW, the gas consumption can be reduced to zero. The corresponding required BSS is nearly 5 GWh, which is a very high capacity, while the same installation capacity of WT requires the BSS for just 1 GWh, which is still a massive capacity. Therefore, this result only considers the technical requirements and does not consider the financial feasibility.

The results in Figure 5 are for cases with only a single renewable generation technology, i.e., either PV or WT but not at the same time. The result shows that PV driven HP is better in terms of total power generation but requires more BSS for energy storage. For this reason, another method is to install PV and WT at the same time to combine their advantages. The gas consumption and the required BSS corresponding to the PV and WT at different rated powers are shown in the 2D surfaces of Figures 6A,B. When the rated power of PV and WT is greater than 7 MW, only 2.5 GWh of BSS can reduce natural gas consumption to zero, which is only half the installed capacity of PV alone.

The previous case only considered the technical requirements for the selection of components in RHS. This section considers the financial cost in reducing CO 2 emissions and optimises the investment capacity of each component for every year. The financial assumptions for the unit cost of each component in future years use the predictions given in Figure 2. Due to the limited AIC, more investment on one component means less investment on other components. Thus, its sizing needs to be optimised to find the optimal scheme to achieve the minimum CO 2 emission. The variables for optimisation include the capacity of PV, WT, BSS, HP, TES considering their constraints described in Optimisation problem. For this multi-parameter Optimisation problem, the optimisation method uses the PSO method described in Optimisation approach to find the most suitable component size to minimize the total CO 2 emissions. Since the financial assumption for the unit cost of each component in future years predicted in Figure 2 is in every 5 years, the optimisation purpose is the total installed amount of each component in every 5 years. The indicator for the optimisation is the total amount of CO 2 emissions in 5 years. After that, the optimal size of each period is based on the installation amount of each component at the previous period and re-run the optimisation approach according to AIC. Therefore, the optimisation results of every 5 years are based on the optimal solution for the current period. It is ensured that the size of the components in the RHS will increase, and the CO 2 emissions will decrease year by year, and finally achieve the optimal result.

The optimisation result is shown in Figure 7, which shows the annual investment capacity of each component in the RHS based on different investment budgets. For example, if the university plans to spend 100,000 pounds per year, as shown by the blue line, the optimal capacity of PV, WT, HP, BSS, and TES for investment according to the optimisation results are given in Figures 7A–E. Using the optimal scheme, the least CO 2 emissions can be obtained, as shown in Figure 7F. In the blue line, the CO 2 emission will still not be reduced to zero by 2050 but will be reduced from 2.7 tons/year to 0.7 tons/year. If the university plans to use the grid independent RHS to completely eliminate CO 2 emissions by 2050, it will need at least 400,000 pounds per year. The optimisation results give a preliminary idea of how to construct the RHS. From the results, in the first 10 years, HP, as the main component of RHS, has the fastest growth in its investment capacity. At least after 2030, the investment of TES will have an effective impact on reducing CO 2 emissions. The investment capacity of PV, WT, and BSS is increased according to the local renewable energy resources and the unit cost of each component to achieve the optimal application of renewable energy.

Compared with different AIC, the more investment, the less CO 2 emissions are obtained, as shown in Figure 8A. Numerically, if RHS is not invested, the traditional heating system will generate at least 8 × 10 5  tons of CO 2 from 2020 to 2050. The higher the AIC of RHS, the less total CO 2 emissions are obtained. If the AIC reaches £300,000 per year, the total CO 2 emissions from 2020 to 2050 can be reduced to 1.6 × 10 5  tons, which is only 20% of continuous using the traditional heating system.

The estimated CO 2 emissions are based on the output of a district heating network model, which has been calibrated based on data from the past few years. The actual CO 2 emissions of the entire university are much more complicated and include commuting and business travel (Duncan and Haydon, 2019). Therefore, the estimated CO 2 emissions in this article are only for preliminary verification of the feasibility and effectiveness of RHS. In addition, the investment cost shown in this paper only considers the predicted unit price of devices. This makes the price and investment amount of the RHS system look much lower than expected. The actual cost will also include other cost, such as the installation fee, maintenance fee, ground rent, etc. These will greatly increase the complexity of prediction and optimisation, so this paper does not consider these factors.

In addition to the total CO 2 emissions, the operating cost of natural gas and CO 2 is another cost of the heating system, as shown in Figure 8B. If a traditional heating system is used, the total operating cost of natural gas and CO 2 is approximately 260 million pounds. And if the AIC of RHS reaches £300 k per annum for 30 years from 2020 to 2050, the total operating cost will be reduced to 36 million pounds, which is just 14% of the operating cost of the traditional heating system. This shows that the investment of RHS can reduce operating costs as well, which is a further benefit and profit for this investment.