Minimizing Wind Power Curtailment and Carbon Emissions by Power to Heat Sector Coupling—A Stackelberg Game Approach

This paper proposes a Stackelberg game approach to minimize both the wind power imbalances and carbon emissions by harnessing demand response (DR) of residential heating loads fed by electricity or district heating (DH) options. The problem is formulated as a bilevel optimization model. An aggregator (leader) at the upper level aims at minimizing the mismatch between electricity demand and wind power while mitigating the carbon emissions arising from the DH system. The aggregator owns a wind farm and is responsible for controlling the DH generation through a combustion-based source and a deep well heat pump system (DWHP), which converts power into heat by utilizing wind power and replacing combustion-based DH. The aggregator submits bonuses to the households (followers) at the lower level incentivizing them to modify their consumption profile. The households receive the bonus offers and consequently decide how to optimize their net electricity and DH energy payments via an upward or downward DR strategy. The uncertainties associated with wind power and heating loads are considered using a stochastic programming framework. Long term thermal performance of the DWHP is studied separately. Results prove that the proposed bilevel framework enables significant reductions in wind power imbalances, carbon emissions, and energy payments.


Index and sets t, T
Time step index and set. N 1 , N 2 Set for electrically heated houses and district heated houses, respectively. N = N 1 + N 2 Set of all houses. n ∈ N House index.

sw, SW
Wind power scenario index and set.

sc, SC
Outdoor temperature scenario index and set. i, I Breakpoint index and set for electrically heated houses. j, J Breakpoint index and set for district heated houses. k, K slope index and set for piece-wise linearization of (P dwhp sw,sc,t ) 2 The associate editor coordinating the review of this manuscript and approving it for publication was Chenghong Gu . Break points for SOS2 variables ranging from 0 to 1 c 1 , c 2 Bonus coefficient (e/W 2 ) r k Slope of the k th segment of P Wind power curtailed in scenario (sw,sc) and time t P lc sw,sc,t Electrical load curtailed in scenario (sw,sc) and time t P l sw,sc,t Total electricity demand in scenario (sw,sc) and time t B1 sw,sc,t,n Bonus earned by household n ∈ N 1 in scenario (sw,sc) and time t B2 sw,sc,t,n Bonus earned by household n ∈ N 2 in scenario (sw,sc) and time t P dwhp sw,sc,t

I. INTRODUCTION
The prevailing usage of fossil fuels in electricity as well as the district heating (DH) generation sectors resulted in a drastic increase in carbon emissions. To mitigate its harmful effects, the European Union (EU) has set up a roadmap towards carbon neutrality in which zero net emissions need to be achieved by the year 2050 [1]. On this basis, member states are striving hard to achieve individual emission targets defined for them. Similarly, the Paris Climate Accord aims to limit the global temperature rise to 1.5 • C. Such ambitious goals require to balance carbon emissions and carbon sinks by the middle of this century. Moreover, the Finnish government program also aims to achieve carbon neutrality already by the year 2035 [2]. It implies huge emission reductions in the following years which calls for a massive energy transition ahead.
To accelerate the decarbonization objective of energy systems, the penetration of renewable energy sources (RESs) needs to be substantially increased in the electricity generation mix worldwide. Although there are clean renewables available, for example, hydro-power, a downside is that it is geography dependent and requires a long time to build the required infrastructure. As a consequence, other renewable generation technologies such as wind and solar power, are developing rapidly for vast deployment in power systems. Such clean energy sources replace a high percentage of conventional fossil fuel-based power plants and mitigate emissions at almost negligible operating costs. However, they are un-dispatchable since they have highly volatile and intermittent power generation levels. Their integration into the current electrical grid poses a big challenge that requires smart grid solutions. Without such solutions, the power system operation would suffer from enormous renewable generation curtailments that undermine the optimal development of such renewable technologies [3].
In order to address the aforementioned issue of renewablebased generation curtailments, sector coupling of electricity and heat has been extensively proposed in the literature. This solution has simultaneous benefits of renewable energy integration as well as carbon emission reductions. In this area, the authors in [4] demonstrated that power to heat (P2H) conversion by employing renewable energy is the least cost strategy for mitigating the carbon emissions in energy systems. The contribution also pointed out the need for heat storage in addition to the hydro-power. Likewise, in [5], a high wind power generation for the case of Helsinki was considered while studying the impact of P2H coupling in the energy system. It was established that even employing only heat pumps would significantly reduce the emissions of the DH. Performing macro and micro-level energy system optimizations in [6] for the case of Helsinki, revealed that wind power utilization could be enhanced from 20% to 37% by P2H coupling using electric boilers only. An energy system optimization in [7], considering a 50% share of wind power for Helsinki city, demonstrated the importance of combining the P2H and thermal storage in absorbing all the surplus wind generation. The aforementioned works targeted the Nordic climate since the DH load is high due to the extreme cold conditions. The promising results obtained from such studies endorse the P2H conversion and utilization of the DH system as a power sink to increase the deployment of renewable energy and, consequently, reduce emissions.
In the literature, P2H coupling is mostly suggested to perform using electric boilers and ground source heat pumps. There exist only a few studies that investigated the scope of deep borehole heat exchangers (BHE) for P2H applications. Most of the ground-coupled heat pumps employ geothermal resource that is less than 300m deep. Contrarily, deep BHEs have 1-2 km depth where the temperature is between 20 and 40 • C. The heat source still requires a heat pump to raise its temperature further for DH applications. The whole system is then called a deep well heat pump system (DWHP). In this respect, a field test of three 2 km deep BHEs located in Xian was conducted in [8]. An average Coefficient of Performance (CoP) of 4.6 was obtained for the whole system. Recently, a comprehensive study of a 2 km deep BHE was performed in [9]. The operation was simulated with varying parameters for the case of Finland for a period of 25 years. It was proved that a 2 km deep coaxial borehole can produce about 110kW (55W/m) heat continuously in steady state, which was reached after 15-20 years of operation.
Despite the firm benefits of P2H coupling, demand response (DR) is also considered a cost-efficient solution to address renewable energy integration. The DR is an efficient and promising tool that aims to provide ramping or balancing power for variable and unpredictable renewable generations [10]. The DR programs are also designed to motivate consumers to shift demand to off-peak hours by offering some monetary incentives. Hence, the DR can contribute to flattening the load profile as well as the deferral of the infrastructure upgrade. The advent of smart grid infrastructure and technologies like home energy management system (HEMS), automatic meter reading (AMR), etc., have reinforced how end-users participate in energy scheduling, provide requested DR actions, and respond to real-time prices. A variety of domestic appliances in a smart home can be managed using HEMS to provide a fast response during critical hours [11].
Among the DR loads, thermostatically controlled loads (TCLs) such as space heating, ventilation and air conditioning, heater for domestic hot water and electric iron, etc. capture prime positions. The main reason is their inherent flexibility in providing load preponement or postponement features that are least noticeable to the users. In other words, it can provide high volumes of power sink and power resource capability by consuming more available power from RESs and deferring consumption, respectively. In the case of TCLs, the user comfort level is a function of temperature deadband only. Whereas, defining acceptable user comfort limits for other loads is not an easy task. Moreover, only a smart thermostat is required to unleash the flexibility of a TCL [12]. Hence, TCLs can effectively react to variable prices and monetary incentives received from the aggregator or service provider.
A variety of works have acknowledged the benefits offered by the DR options of TCLs. For instance, the potential of domestic electric water heater and space heating loads in both centralized and distributed control management was realized in [13], aiming to achieve a desirable load profile in the system operator's context. The country-wide flexibility potential of domestic space heating loads was considered in [14] in order to quantify the optimal size of centralized storage for renewable energy integration. In [15], [16], the DR capability of space heating and electric water heater loads was studied to minimize the power curtailments in a residential microgrid.
In [17], the flexibility of space heating loads was envisioned for distribution network management during high penetrations of wind power. In [18], the authors probed the potential of space heating load in coordination with building inertia to propose an optimal bidding strategy for day-ahead and real-time markets in a stochastic environment. Authors in [19] demonstrated the benefits of DR of space heating integrated with a real-time thermal rating of the distribution network for wind generation balancing. Similarly, in [20], significant improvements in end user's energy procurement costs and wind power curtailments in a distribution network were obtained by unleashing the DR of space heating loads. The work [21] attempted to maximize the utilization of on-site solar power production with the aid of TCLs. In [22], the aggregator traded the flexibility of residential heating loads in the day-ahead Nordic market to minimize its cost while analyzing both the price-taker and price-maker strategies.
Almost all of the aforementioned works studying DR considered direct load control mechanism where only one entity, such as load aggregator or service provider, was authorized to adjust the demand of end-users and the same was involved in the decision-making process. Such a framework is best suited to network management like frequency containment reserves (FCR) or balance management. It not only compromises the privacy of end-users but also sacrifices their comfort levels since the decision process is more tailored towards the aggregator objectives. Contrarily, the indirect load control method is based on interacting with the end-users to activate the DR in exchange for some monetary incentive. This methodology tends to respect the privacy of end-users. Hence, a proper interaction between the two players needs to be studied as it is in the best interest of both. This approach leads to the game theory concept, where rational decision-makers interact with each other based on their strategies. Among game theory concepts, one of the most applied techniques in power system studies is the Stackelberg game theory, which is suitable for solving hierarchical problems. The two players of this game are called leader and follower, and they compete with each other. A classic example is the interaction of the aggregator (leader) and the end-users (several followers) in a smart grid. The two players have different objectives and strategies, and they iteratively strive to reach an equilibrium that fits best for both players.
Game-theoretic approaches have been applied in the power system context to integrate prosumers in the decision-making process. For instance, in [23], a pricing mechanism based on Stackelberg game model was proposed aiming to maximize the aggregator's profit at the upper level while minimizing the charging cost of EVs at the lower level. The problem was solved iteratively until achieving an equilibrium point. Similarly, a real-time pricing strategy based on the Stackelberg game was proposed in [24] in which equilibrium was obtained by trading energy between prosumers and consumers aiming to maximize the prosumer's profit and minimize the consumer's bill. The interaction between the regional control center and village microgrids was studied in [25] to achieve Stackelberg equilibrium. In that problem, the control center's strategy was the price offering and the microgrid was established to schedule the energy consumption. The load models considered in these works were too simple and could not encompass practical aspects, such as user comfort priority. A game-based model between a retailer and residential users was proposed in [26] to minimize the balancing power costs and risks of a retailer by incentivizing the end-users. Notwithstanding, a direct load control approach, which sacrifices the privacy of the end-users, was used. In [27], a Stackelberg equilibrium was obtained between multiple utilities and end-users via a distributed algorithm that aimed to maximize the benefit of each entity. The proposed methodology required communication links between end-users and utilities that suffered from privacy issues as well as security of supply.
In this work, most of the above-mentioned shortcomings or drawbacks are addressed by proposing a Stackelberg game between an aggregator and households through a bonus-based DR. It is assumed that the aggregator, as a leader, owns a wind farm and operates a DWHP already coupled with the DH network. There are several households with direct electricity heating (DEH) or DH demand, that act as followers in a medium-scale microgrid (MG). The objective of the aggregator is to maximize wind power utilization while minimizing the carbon emissions originating from the DH network. To accomplish these goals, the aggregator offers a bonus to households to adjust their load according to the wind power profile. The excess wind power, which would be curtailed otherwise, can also be utilized by DWHP for P2H conversion. The households at the lower level focus on their energy procurement costs and thermal comfort levels. The households' objective is to optimize their energy payments individually by modifying their heating demand against the offered incentive by interacting with the aggregator in an iterative routine to reach an equilibrium.
The primary reason for choosing space heating loads for this work is their high share in the annual energy consumption of Northern Europe. In 2018, space heating alone formed 26% of all energy used in Finland, and it coincides with the Finnish peak demand as well. In Northern Europe, the DH accounts for the largest part of the heating, and in large cities, it is still heavily based on fossil fuels. For instance, about 90% of all the DH in Helsinki is produced by CHP plants powered by natural gas, coal, and wood pellets. The specific carbon emission factor for heating in Helsinki escalated to 198g/kWh in 2019 compared to 158g/kWh in 2018 [28]. Thus motivated by the increasing heat energy consumption and its consequences in the energy sector, this work investigates the impact of space heating load flexibility in the energy sector.
The unique aspects of this work are listed as under: • A Stackelberg Game is proposed to consider the interaction between an aggregator as the leader and several households as followers. The problem is formulated as a bi-level optimization framework that simultaneously minimizes the power imbalances, mitigates carbon emissions aiming at reducing the energy procurement costs for households.
• In order to avoid the iterative process to reach equilibrium for both participants, the strong duality condition is applied to recast the proposed bilevel problem into an equivalent single-level problem which is then solved to obtain the optimal global solution. VOLUME 8, 2020 • The P2H conversion is performed using a novel DWHP. The thermal behavior of the DWHP is analyzed separately in Finnish conditions by simulating over a 25-year period.
• The heating load is modeled using a two-capacity building thermal model, which not only takes into account the thermal comfort priority of each end-user but also transforms it into offered DR.
• We consider variable electricity prices and flat DH tariff according to current practices. Hence, we also consider the event when high wind generation and high electricity prices co-exist.
• The wind power and heating loads are driven by uncertain parameters such as wind speed and outdoor temperature, respectively. A finite number of discrete realizations for each uncertain parameter is considered in a stochastic framework. The remainder of the paper is organized as follows. Section II discusses the system description and models used. Moreover, a description of the DWHP performance is also provided. Section III presents the mathematical optimization model of the proposed Stackelberg game. Simulation parameters and results are described in Section IV. Finally, Section V concludes the paper.

II. SYSTEM DESCRIPTION A. TWO-CAPACITY BUILDING MODEL
The heating or cooling load in a detached house is modeled using a two-capacity thermal model [3], [16]. The model has good accuracy in capturing the indoor air dynamics relative to the outdoor temperature changes. As the name suggests, the model uses two capacitances. One of the capacities is distributed to the building structure, C m , and the other is allocated to the indoor air, C a . There are two unknown temperatures, building mass temperature T m and indoor ambient temperature T a . The temperature node points are attached through heat conductance, or heat capacity in case of air flows. The infiltration air is connected between outdoor temperature and indoor temperature and there is no warming of the air in the building mass. The windows have an insignificant mass compared to the rest of the envelope. The heat capacity of building mass C m is much higher than that of indoor air C a , but C a still has a significant role in studying indoor air dynamics. The heating or cooling power of the HVAC unit is of convective type and therefore allocated to indoor air node. The air heating unit is set to operate at constant temperature T x . The model is depicted in Figure 1.
The HVAC unit, depending on the indoor and outdoor temperature, consumes electrical power in order to maintain indoor ambient temperature T a equal to a set point temperature, say 21 • C. The DR of space heating load means the optimal adjustment of the air heating unit's consumption. During higher power consumption, building thermal masses store energy and well-insulated buildings act as a small storage buffer. Alternatively, during lower consumption, building thermal masses release this stored energy. This adjustment  in the power consumption is also reflected in T a . The twocapacity model accurately captures such dynamics of indoor and building mass temperatures. This way, the heating load in coordination with thermal inertia enables to effectively cope with the variable renewable generation. This flexibility comes at the cost of a slight sacrifice of thermal comfort level. Usually, heating loads are used in tandem with thermal storage which needs charging and prevents the loss of thermal comfort.
The building parameters of the two-capacity model are calibrated using the dynamic building energy simulation tool IDA-ICE [11], which is used for studying indoor climate and energy consumption of buildings. The considered house is a single-family, two-floor with the medium weight passive structure according to the Finnish guidelines. More details of the house can be found in [29]. The indoor temperature preference was 21 • C [30]. For calibration, the heating power is discontinued for six hours. The variance between the response obtained from IDA-ICE and the derived two-capacity model is then minimized. The best combination of parameters was derived using conventional search approach. The calibration is performed at three different outdoor temperatures i.e., +10 • C, 0 • C, and -10 • C. The calibration method is depicted in Figure 2.

B. DEEP WELL HEAT PUMP SYSTEM
The deep-heat energy system consists of a combination of a deep-well and a heat pump system. The deep-well comprises a 1-2 km drilled borehole with a coaxial pipe having an outer and inner pipe with two opposite flows as shown in Figure 3.
Here the insulation is placed between the two pipes to reduce thermal short-circuiting, i.e., unwanted heat exchange between forward and return flows which also contributes to heat losses. The heat extraction works as follows: -a fluid with low temperature is injected into the well along the outer pipe, and it flows back upwards through the central pipe. During this cycle, the fluid heats up and provides thermal energy to the heat pump, which raises the temperature of the fluid to the desired level, e.g., suitable for a district heating system. The deep-heat energy system can be used as a heat generation system and/or as a thermal storage. When the well is used as a thermal storage, hot fluid is injected to the well downwards the middle pipe, and it returns back cooled fluid along the outermost pipe. In this way the bedrock surrounding the borehole can be charged.
The advantage of a deep well over a normal 300m borehole is that the deep well provides significantly more heat. In previous work [9], the annual heat energy output of a deep well was analysed in Finnish geological conditions, showing that the deep well produces up to 30 times more heat than a conventional 300m deep borehole.
Because a deep well is significantly deeper than a conventional borehole, the thermal short-circuiting issue needs to be effectively addressed. Compared to a conventional borehole employing a U-tube, a coaxial pipe has a larger cross-sectional area allowing also higher mass flow rates. With a higher mass flow rate, the thermal shunting between the forward and return flows also decreases. However, the mass flow cannot be indefinitely increased, because the pressure drop and hence the pumping power required increases with the flow rate. Figure 4 shows the pressure drop and power required by the pump at different mass flow rates [9].
The thermal short circuiting between the upstream and downstream flows can be reduced by reducing the thermal conductivity between the pipes. For pressure losses and pumping capacities to be reasonable, thermal short-circuiting can be reduced by using a vacuum tube in a deep well. A conventional 300m deep borehole typically employs a plastic U-tube, with a thermal conductivity of 0.42 W/mK. The thermal conductivity of the vacuum tube is 0.02 W/mK so the thermal conductivity is significantly reduced.  The long-term thermal performance of the DWHP system is modeled with the COMSOL Multiphysics programme, which numerically solves the coupled heat and mass transfer problem describing the deep-well physics. The deep-well is modeled in two parts consisting of the coaxial borehole exchanger and the surrounding bedrock. Both models are physically coupled and simulated for a 25-year period with a time step of 7.5 hours. A detailed description of the simulation model and physical equations is given in [9]. The physical parameters of the borehole are given in Table 1.
In this work, the operation sequence is six months of heat extraction, followed by six months of heat injection. This sequence is adopted because the Finnish peak load occurs in the winter season when heating demand is high. Contrarily, there is no space heating load in the summer season, but at the same time, there is a relatively higher excess wind generation. The heat injection during the summer season regenerates the ground, which is realized through power to heat schemes making use of excess wind power. Such thermal regeneration improves the heat output in the extraction mode, i.e., during the winter season. The interested readers are referred to [9] for detailed comparisons of BHE with and without heat injection. Figure 5 demonstrates the dynamics of the fluid temperature in the pipes of a BHE with a 6-month extraction and 6-month injection cycle. Figure 5(a) shows the extraction phase: the fluid is injected at 6 • C through the outer pipe, returning back through the central pipe at a higher temperature, on average at about 14 • C. In Figure 5(b), during the heat injection from P2H, the fluid is injected at 70 • C through the middle pipe and the fluid returns back through the outer pipe at 41-46 • C. This temperature drop demonstrates that heat is VOLUME 8, 2020  absorbed by the surrounding rock. From Figure 5(a), one can see that more heat is extracted at the beginning of the extraction cycle with a declining trend over time. Similarly, more heat is absorbed by the BHE at the beginning of the injection cycle, in Figure 5(b). This kind of saturation phenomenon is explained by the limited thermal conductivity of the rock. The temperature profiles depicted in Figure 5 can easily be transformed into respective thermal powers.
The operational characteristics of the DWHP is demonstrated in Figure 6 using input electrical and output thermal power, which is available from the condenser of the heat pump. The performance curve is obtained from simulations with nine different mass flow rates ranging from 0.5kg/s to 15kg/s. The curve has been optimized in terms of finding an optimal trade-off between the mass flow rate, which affects the pumping power, and the output heat power. In Figure 6, the input electrical power is the sum of the pumping power and the power required by the heat pump. The output heat power showcased in Figure 6 corresponds to the average steady-state values attained in the last six months of a 25-year simulation period against different mass flow rates. Figure 6 shows that the thermal output increases linearly with input electrical power up to about 77kW, which corresponds to a fluid flow rate of about 5kg/s. After 77kW, the system coefficient of performance (CoP) starts to decline.
Assuming a CoP of 3.5 at the beginning, the DWHP CoP decreases from 3.5 to 2.31, when the mass flow rate increases from 0.5kg/s to 15kg/s. This corresponds to a decline of 33.98% in performance relative to the maximum CoP. Up to a mass flow rate of 6kg/s (input power 85.7kW), the decline is just 4.66%.

C. STACKELBERG GAME
The leader in this game is the aggregator who owns a wind farm and controls the DH generation. The DH network is assumed to be fed by two sources, i.e., combustion-based generation and DWHP, both controllable by the aggregator. The objective of the aggregator is two-fold, i.e., maximizing the wind power utilization and minimizing the carbon emissions originating from the DH sector. The aggregator forecasts wind power and aims to optimally balance the generation with electrical demand and DWHP operation. Since the aggregator has no control over wind generation, and the DWHP operation alone is not flexible enough, households need to be offered a monetary benefit to modify their energy consumption according to the wind generation profile.
At the lower level, there are several households as followers with either electricity or both electricity and DH demand. All households have two types of load, i.e., critical and heating loads. The critical load can neither shift nor curtailed against any price or incentive. Contrarily, the heating load is flexible and participates in the bonus framework. Households are further segregated depending on direct electric heating (DEH) and DH loads. Houses with DEH loads are also equipped with a partial thermal storage, which is integrated into the heating unit. This thermal storage space heating system is considered a ubiquitous element in the Nordic household sector [14], [20]. For electricity supply, households are assumed to follow a spot-price based contract, which is a common practice in the Nordic countries. Note that in the Nordic electricity market, hourly day-ahead spot prices are announced at 12 noon each day. Hence, the households know the actual spot prices 12 hours in advance. Contrarily, DH is subject to a flat tariff irrespective of the heat source. The operation of the DWHP not only contributes towards the aggregator objectives, but also enables DH loads to interact with the aggregator and optimize their energy costs despite flat DH tariff.
Before the Stackelberg game interaction begins, households with DEH optimally schedule their heating loads according to the day-ahead electricity spot prices and outdoor temperature forecast whereas DH loads depend on temperature forecast only. The households then further reduce their energy payments during Stackelberg game interaction. The aggregator needs to compete with electricity prices when it comes to DEH households. All households are assumed to be equipped with HEMS to receive bonus information from the aggregator and convey the corresponding DR actions.
In this game, the strategy set of players are defined at first. The aggregator's strategy is the bonus offered to households, and the households' strategy is to offer upward or downward DR by deviating from day-ahead scheduling in relation to the offered bonus. This DR is in the form of electrical or heating power depending on the type of heating. The interaction between players begins with the aggregator's strategy by sending a bonus signal to each household. Each household reacts to the bonus and turns up with a strategy that is optimal for its own energy usage. This strategy is then communicated to the aggregator through HEMS, and the household waits for the aggregator's response. At a subsequent phase, the aggregator updates its bonus offering on the basis of the DR action received from each household and sends this revised bonus offer to each household. This process of interaction continues iteratively until an equilibrium is reached.
In this work, the DEH loads are indexed over subscript n ∈ N 1 while the DH loads are indexed over subscript n ∈ N 2 such that N 1 + N 2 = N and N 1 ∩ N 2 = {}. Let P sw,sc,t,n , ∀n ∈ N 1 and Q sw,sc,t,n , ∀n ∈ N 2 be the strategy sets of DEH and DH loads for each scenario and time step, respectively. The corresponding strategy sets of the aggregator are B1 sw,sc,t,n = c 1 P h sw,sc,t,n , ∀n ∈ N 1 and B2 sw,sc,t,n = c 2 Q h sw,sc,t,n , ∀n ∈ N 2 . These two terms represent the bonus given by the aggregator to the DEH and DH households, respectively and these bonuses depend on the DR volume offered by the corresponding households. It implies that the household earns more bonuses by contributing more to the DR framework.

III. OPTIMIZATION FORMULATION
This section presents the bi-level optimization framework for the proposed interaction between the aggregator and house-holds. The aggregator objective is formulated as in equation (1) and it is subject to constraints (2)- (22), (the equations (1)- (19) are shown at the bottom of the next page and equations (20)- (22) are shown at the bottom of the page ten). Whereas, the household level problem is formulated in (7)- (22). Please note that both upper and lower level problems are formulated as minimization problems. The household in lower level aims to simultaneously minimize its energy payment and maximize the bonus.
In the above problem, the aggregator's objective (1) has two parts. The first part accounts for the electrical power imbalances, and the second term deals with costs associated with the bonus offerings and carbon taxes. The terms P wc sw,sc,t and P lc sw,sc,t represent total wind and load curtailment in each scenario and time step t respectively. The terms B1 sw,sc,t,n P h sw,sc,t,n and B2 sw,sc,t,n Q h sw,sc,t,n represent the bonus earned by each DEH and DH household in each scenario and time t respectively. The last term µ.v.Q th sw,sc,t accounts for the cost of the carbon emissions in the DH system. The power balance for electrical demand is maintained by (2). Constraints (3) calculates the total electricity demand in each time slot, and (4) is responsible for satisfying the DH demand with two sources. The input and output relationship of DWHP is devised in (5), while its input operation limits are defined in (6) according to Figure 6. The binary variable ensures that DWHP operates within specified input limits. The objective function of the households is presented in (7) which has two parts. The first part in (7) represents the net energy procurement costs of DEH households while the second part accounts for the net energy payments of DH households. Equations (8) and (9) stand for the discrete version of the two-capacity building model that is used to calculate the space heating load and capture the dynamics of indoor temperature in each house. Constraint (10) specifies the thermal comfort preferences for each household in the form of upper and lower bounds for indoor ambient temperature while heating power is capped in (11). To ensure that the two-capacity model operates in heating mode only, (12) is applied. Constraint (13) establishes that the total heating demand for each household remains the same with and without DR activation over the simulation horizon. In other words, the net flexibility offered by the heating load for each household is zero over the simulation horizon. The magnitude of upward and downward DR of heating load with respect to its day-ahead scheduling is restricted by (14), while (15) applies the same to the absolute value of DR. The evolution of the thermal storage for space heating is captured in (16). The thermal storage is charged using electrical power, and it stores energy in the form of heat [17], [19]. The discharging takes place in accordance with the two-capacity building model. For the sake of simplicity, thermal losses of the storage are assumed to be negligible. The allowable limits of thermal storage capacity are defined in (17). To maintain the longer-term thermal balance, constraint (18) guarantees that the net flexibility offered by each thermal storage over the horizon is zero. Constraint (19) maintains that the charging power of thermal storage is always non-negative, while (20) sets an upper limit to it. Constraint (21) states that the DR of thermal storage is within specified bounds, whereas the corresponding absolute value is limited in (22).
The presented bi-level model (1)- (22) is non-linear due to the quadratic and absolute terms, for instance, in (1), (5) and (7). In order to find the equivalent single level of this bilevel problem, first, we need to get rid of the non-convexities, and a straightforward way is to use some appropriate linearization techniques. The absolute term can be easily linearized using two positive auxiliary variables as in (23)-(24) and (32)-(33) [16]. To linearize the quadratic terms resulting from the multiplication of two absolute value functions such as P h sw,sc,t,n in (1) and (7), special ordered sets of type two (SOS2) variables are employed as in (25)- (31). SOS2 is a set of variables in which not more than two variables take non-zero values and these variables are adjacent to each other. The necessary condition to apply SOS2 requires the non-negativity of the base, which is already fulfilled by reformulating the absolute function using two positive auxiliary variables in (24) [31]. In SOS2 formulation, the continuous variable is first expressed as a sum of finite number of breakpoints weighted with the SOS2 variables as in (26). The breakpoints are chosen in accordance with (22). The square of the continuous variable is then computed as in (27). The constraint (28) requires that the sum of SOS2 variables is equal to unity. Similar formulation is performed for Q h sw,sc,t,n which is given by (32)-(40). Please note that the last terms in both the (25) and (34) the two auxiliary variables always take the value zero which render their product as zero. The quadratic term appearing in (5) is linearized using piece wise approximation as shown in (41)-(44) [32]. In this method, the quadratic curve is first represented using a finite number of linear segments as in (41) ; ∀sw ∈ SW , i z sw,sc,t,n,i = 1; ∀sw ∈ SW , sc ∈ SC, t ∈ T , ; ∀sw ∈ SW , sc ∈ SC, t ∈ T , n ∈ N , j ∈ J (36) j w sw,sc,t,n,j = 1; ∀sw ∈ SW , sc ∈ SC, Now, the above-presented bi-level model (1)-(44) has a linear lower level problem, which enables reformulating the bi-level problem into its equivalent single-level problem. This reformulation is useful to avoid the iterative method to solve the bi-level problem [31], [33]. In order to reformulate (1)-(44) into its equivalent single-level problem, the dual of the lower level problem is determined using the dual variable associated with each of the constraints (8)- (22). P h o sc,t,n + P h sw,sc,t,n ≤ P max n ; (ω sw,sc,t,n ), ∀sw ∈ SW , sc ∈ SC, t ∈ T , ∀n ∈ N 1 (20) − P max n ≤ P h sw,sc,t,n ≤ P max n ; (δ lo sw,sc,t,n , δ up sw,sc,t,n ), ∀sw ∈ SW , sc ∈ SC, t ∈ T , ∀n ∈ N 1 (21) P h sw,sc,t,n ≤ P max n ; (π sw,sc,t,n ), ∀sw ∈ SW , sc ∈ SC, t ∈ T , ∀n ∈ N 1 (22) VOLUME 8, 2020 Then the lower level problem is substituted with primal constraints, dual constraints and the strong duality condition. The equivalent single-level problem is presented in (45)-(58), as shown at the bottom of the page and (59)-(63), as shown at the bottom of the next page.
In the reformulation (45)-(63), the absolute value functions will be substituted with their corresponding auxiliary variables presented in (23)- (24) and (32)  − σ up sw,sc,t,n + ψ sw,sc,t,n − ς sw,sc,t,n /n d = 0; ∀sw ∈ SW , sc ∈ SC, t ∈ T , ∀n ∈ N 1 (54) ς sw,sc,t,n − ς sw,sc,t+1,n + τ lo sw,sc,t,n + τ up sw,sc,t,n = 0; ∀sw ∈ SW , sc ∈ SC, t = 1 . . . T − 1, ∀n ∈ N 1 (55) ς sw,sc,T ,n + τ lo sw,sc,T ,n + τ up sw,sc,T ,n = 0; ∀sw ∈ SW , sc ∈ SC, ∀n ∈ N 1 (56) −n c ς sw,sc,t,n + ξ sw,sc,n + κ sw,sc,t,n + ω sw,sc,t,n + δ lo sw,sc,t,n + δ up sw,sc,t,n + π sw,sc,t,n = λ e t − B1 sw,sc,t,n ; ∀sw ∈ SW , sc ∈ SC, t ∈ T , ∀n ∈ N 1 (57) n c ς sw,sc,t,n − ξ sw,sc,n − κ sw,sc,t,n − ω sw,sc,t,n − δ lo sw,sc,t,n − δ up sw,sc,t,n + π sw,sc,t,n = −λ e t − B1 sw,sc,t,n ; ∀sw ∈ SW , sc ∈ SC, t ∈ T , ∀n ∈ N 1 (58) obtained from [34]. Wind speed scenarios are adopted from the study [35] which simulates 100 yearly scenarios of wind speed at varying heights using a statistical approach without considering any site-specific data for the case of Helsinki, Finland. The corresponding wind power can be generated using cut-in speed and rated speed of the turbine. The rated wind power is considered as 400 kW. The hourly outdoor temperature profile is obtained from [36]. The scenarios of outdoor temperature are generated in such a way that the deviation from the actual profile increases linearly through the horizon and is set to ±20% at the end [11], [16]. We consider four scenarios for wind and two for outdoor temperature, which makes a total of eight scenarios that are all equiprobable. A fixed DH tariff defined for the winter season is used [37]. The input data are shown in Figure 7. It is assumed that there are 100 single-family twofloor detached houses served by the aggregator. The first 50 households use DEH and are also equipped with a small thermal storage. The remaining households use DH without any storage. Apart from that, all houses have critical load that cannot shift in time. The critical load for a typical house is demonstrated in Figure 8 [16]. The building occupancy is considered continuous for the studied horizon. The load and house parameters are uniformly distributed among intervals as listed in Table 2.
The COMSOL simulation results of the DWHP operation is mapped accurately using a quadratic regression curve as shown in Figure 6. The fitting error is negligible. The minimum input required for the DWHP operation is 18.5kW for DH application. For the SOS2 variables, 11 breakpoints are used to calculate the square of continuous variables. The cost of carbon emissions in (45) is calculated using a specific emission factor for wood pellets which is a common fuel in the Nordic region [38]. The carbon tax value applicable to the Nordic region is used to convert resulting emissions into a monetary value [39]. The calibrated two-capacity sw sc α up sw,sc,t,n , χ sw,sc,t,n , σ up sw,sc,t,n , ψ sw,sc,t,n ≤ 0; ∀sw ∈ SW , sc ∈ SC, t ∈ T , n ∈ N (60) α lo sw,sc,t,n , γ sw,sc,t,n , σ lo sw,sc,t,n ≥ 0; ∀sw ∈ SW , sc ∈ SC, t ∈ T , n ∈ N (61) τ up sw,sc,t,n , ω sw,sc,t,n , δ up sw,sc,t,n , π sw,sc,t,n ≤ 0; ∀sw ∈ SW , sc ∈ SC, t ∈ T , n ∈ N 1 (62) τ lo sw,sc,t,n , κ sw,sc,t,n , δ lo sw,sc,t,n ≥ 0; ∀sw ∈ SW , sc ∈ SC, t ∈ T , n ∈ N 1 (63) VOLUME 8, 2020  building parameters are given in Table 3. The time resolution is hourly due to the hourly day-ahead electricity prices. The mixed-integer linear programming model (45)-(63) is simulated in the Matlab-GAMS environment while the optimization problem is solved via CPLEX solver. The model is implemented on a Windows desktop computer with a 3.4GHz Intel Xeon processor and 16GB RAM.

B. SIMULATION RESULTS
To show the effectiveness of the proposed model, two cases are simulated which are listed below: I. No DR is activated and hence there is no interaction between the players. The aggregator operates the DWHP system to fulfill its objectives. This case serves as the comparison benchmark. II. The whole problem (45)-(63) is simulated.  For Case I, we first generate the heating load profiles using the two-capacity model given in (8)-(9) and the thermal storage in (16). It is assumed that after DA prices are revealed, households forecast outdoor temperature to schedule their heating loads. The heating demands of the DH households follow the outdoor temperature only, whereas the DEH households optimize their storage operation according to the DA electricity prices as well as the outdoor temperature. For this purpose, the hourly indoor temperature for each household is kept at 21 • C. The remaining parameters are the same as reported in Table 2. The resulting heating power and storage charging profiles for the two temperature scenarios are depicted in Figure 9. It is visible in the Figure that the charging is done when the price is low. This charging power is stored as thermal energy in storage, usually a hot water tank. The storage is discharged in the form of heating power required to maintain desired comfort levels, i.e., 21 • C in this case. Unlike electric storage, thermal storage can charge and discharge in the same time slot, since thermal comfort level must be maintained in each time slot. Moreover, total charging power over the horizon is set equal to the total required heating power in proportion to storage efficiency. For the DEH households, the discharging or heating power is the same as that depicted in Figure 9(b).
The simulated total expected energy cost for all the DEH and DH households in Case I is 284.76e and 252.08erespectively. The scenario wise energy cost is  presented in Table 4. In the absence of the proposed framework, the aggregated load profiles portrayed in Figure 9 are considered firm and final, and the aggregator cannot request the households to adjust their demand. The aggregator gets the load information of each household through HEMS, aggregates the demand and schedule DWHP system operation to minimize the power imbalances and carbon emissions. The simulated expected electrical power imbalances and emissions over all scenarios are 2369.82kWh and 609.98kg, respectively. Two out of eight scenarios have been selected to validate the simulation results. These scenarios are (sw = 4, sc = 1) and (sw = 1, sc = 2), herein after called scenario 4 and scenario 5, respectively. Figure 10 illustrates the hourly electrical power balance for the two scenarios. It is worth mentioning that hourly electrical demand consists of the critical electricity demand of both the DEH and DH households as well as the charging power of thermal storages. Similarly, Figure 11 demonstrates the power balance for DH demand which can be satisfied with the DWHP system together with fossil fuel-based thermal generation. In order to account for all the emissions in the DH, it is assumed that there is always ample thermal generation available to satisfy the residual DH demand. By comparing Figures 10 and 11, it can be inferred how the excess wind generation has been utilized by DWHP. Please note that both the DWHP and thermal generation can operate simultaneously within an hour, for instance, in Figure 11(a), hour # 5, 20 and 22. During these hours, the DWHP utilized all the available wind generation and there was no electrical power imbalance, see Figure 10(a), but since the DH demand was higher, so fossil fuel-based generation also needs to be employed.
For Case II, the aggregator can offer a bonus, and households can deviate from their DA schedule portrayed in Figure 9. The individual parameters are followed, as listed in Table 2. It is noteworthy to mention that households are free to choose and tune comfort preferences including the indoor set point temperature. The up or down-regulation power of the DEH households is limited to ± [1.5, 2.1] kW. In order to effectively utilize the thermal comfort band of individual households, the heating power is allowed to mutate in the whole operating range. This feature enables turning the heating loads on or off while building thermal masses store or release heat energy in accordance with the two-capacity model; otherwise, a part of the specified thermal comfort band remains un-utilized. The deviation in heating power is reflected in the indoor ambient temperature. Nevertheless, the indoor temperature is always maintained within the comfort band specified in Table 2.
Simulating the proposed model (45)-(63) brings significant benefits for both the power and heat sectors, compared to Case I. The simulated electrical power imbalances and carbon emissions obtained for each scenario are compared for both cases in Table 5. The proposed interaction mechanism produces remarkable improvement in each aspect and for each scenario, compared to Case I. The expected power imbalances and emissions for all scenarios are 640.56kWh and 325.18kg, respectively. Compared to Case I, the improvement is 72.97% and 46.7%, respectively. This mechanism also allows households to reduce their energy procurement costs by earning a bonus from the aggregator. The numerical VOLUME 8, 2020 results concerning such monetary benefits are summarized in Table 6. Unlike the DH households that have the same total energy cost in scenario sc as in Case I, the DEH households' energy costs increase relative to Case I, which is explained as follows: the sum of regulation powers over the horizon is zero for each DEH household but since such regulation powers are activated in time slots that have different electricity prices and the regulation powers are also subject to these prices. Therefore, these regulation powers result in slight increase in the energy cost. Such additional costs have been compensated in the respective bonus for the DEH households. The expected net energy cost (i.e., energy cost -bonus) for DEH and DH households over all scenarios is reduced to 234.13e and 161.45e, respectively. This leads to a relative expected net cost reduction of 17.78% and 35.95%, respectively, compared to Case I. The DH households obtain relatively more bonus that is explained by the limited participation of thermal storage, flat DH tariff and P2H conversion by the DWHP system. The electrical power balance for Case II is depicted in Figure 12 for the two scenarios (i.e., scenario 4 and 5). It is notable in Figure 12 that the demand closely follows the wind generation, thanks to the DR coordinated with building thermal inertia which also enables more P2H conversion. In scenario 4, the demand exactly matches with wind generation in hour # 1, 5-6, 12, 16 and 18-22. Similarly, scenario 5 secures the exact power balance in hour # 1, 3-6, 9-11, 16-17 and 19-22 as illustrated in Figure 12(b). Figure 13 presents the power balance for the DH obtained in Case II for the two scenarios. By comparing Figures 10 and 13, it can be seen that how excess wind generation has been utilized by offering up-and down-regulating power of the DH demand, which results in reducing the   emissions as well. The up-regulation is activated during hours with excess wind generation, and the DH is satisfied with the DWHP. Contrarily, down-regulation is activated when wind generation was not enough to operate the DWHP or the DWHP alone was unable to satisfy the DH demand. Nonetheless, a widely diversified range of load and comfort parameters have been considered in the simulation. For instance, for hours 7-10 and 17 in Figure 12(a), the excess wind generation could also be utilized by the DWHP if the households had offered even a wider thermal comfort band and rated power of heating units had been higher. The optimal solution requires that the sum of up-and down-regulation powers for each household over the horizon be zero, and it affects the topology of operating the DWHP compared to Case I. Figure 14 demonstrates the total hourly bonus earned by all households in scenario 5. The bonus amount directly expresses their participation in the DR framework at different hours. In scenario 5, the DEH and DH households receive total bonuses of 74.73e and 87.47e, respectively. Moreover, the total daily bonus earned by each household for scenario 5 is depicted in Figure 15. This figure confirms that the DR participation is affected by individual load ratings and the offered thermal comfort band. It is to be noted that house # 1 has the smallest value of physical parameters, which gradually increase until the last house for both the DEH and DH houses, as listed in Table 2. Moreover, for the DEH households, there is relatively less loss of thermal comfort due to considering thermal storage. For these reasons, the daily bonus of the DEH households in Figure 15(a) increases with the size of thermal storage. Contrarily, the DH households' bonus is highly related to their set thermal comfort preferences. The first ten DH houses have a narrow thermal comfort band, i.e., [20], [22] • C, that restricts the DR of heating power which ultimately results in less bonus. The DH houses # 61-80 have the widest comfort band but corresponding load  parameters and house areas are not the highest, as evident in Table 2. Due to this, the bonus of houses # 61-100 has less divergence among themselves.
In order to analyze the simulation results on the household level, we choose one DEH house (house # 22) and one DH house (house # 72) in scenario 4. Their profiles for indoor ambient temperature, thermal charging power, storage level evolution and heating power are illustrated for both cases  in Figures 16 and 17, respectively. The profiles are entirely different for both cases. It is worth to notice that for the DEH household in Figure 16, the indoor ambient temperature follows the corresponding heating power slowly. Contrary to Case I, the heating unit in Case II is turned off in five different time steps, and the ambient temperature drops slowly in those time steps due to insulated house structures. The thermal storage is not allowed to drop below 10% level while the final storage level is always equal to the level before the beginning of the horizon, i.e., 50%. At the first-time step, discharging or heating power is much higher than the charging power in Case II, which drops the storage level from 50% to 22% as depicted in Figure 16(b). Furthermore, the storage capacity considered here is just enough to cover about 25% of the daily heating demand of a house. In Case II, the maximum DR of the DEH houses is limited to ±30% of their rated charging powers, specified in Table 2. For house # 22, it comes out to be ±1.764kW deviation from Case I. In Figure 16(c), the apparent difference between charging powers shows that the household offers full DR capacity in most of the time slots with the intent to reduce its energy payment. In scenario 4, house # 22 receives e1.28 bonus against the daily energy cost of e6.05e.
The heating profile obtained for the DH house # 72 in Figure 17 is more tailored to the DWHP schedule, as illustrated in Figure 13(a). Due to less wind generation, the heating for house # 72 is turned off in the first-time step, which forces the indoor temperature to drop from 21 • C to about 19.5 • C, as shown in Figure 17(a). The heating power is shaped to enable more DWHP operation, which in turn allows to reduce more emissions and earn more bonus. Like Figure 16, Figure 17 also displays that the indoor temperature follows the corresponding heating power slowly. This can be seen in hour # 6-7 and 15-16 when the heating power is transited to its rated power and stays at the same level for a while to reach the upper comfort limit. In scenario # 4, house # 72 earns a total bonus of e2.75 against the daily energy payment of e5.06, i.e., e4.31 for the DH and e0.75 for the critical load.
In the simulation results presented for Case II, the DR participation of the DEH households was limited to ±30% of rated charging powers. Hence, it is crucial to study the impact of varying the allowable DR volume, i.e., P max for the DEH households. Figure 18 characterizes the outcome when P max is altered within a band ranging from zero to full rated capacity, i.e., ±P max in steps of 10%. Figure 18(a) highlights the corresponding effect on the expected carbon emissions and power imbalances computed for all scenarios. Note that 0% capacity of P max does not mean that the DEH households are not equipped with thermal storage, rather it implies that storage charging power is not allowed to change relative to Case I, i.e., the DEH households do not participate in DR framework. However, the corresponding heating power follows the charging power profile according to the framework but without participating in the bonus mechanism. Figure 18(a) also shows that increasing P max for the DEH households decreases the expected power imbalances while increasing the emissions, which is explained as follows: Thermal storage is an additional flexibility besides maintaining thermal comfort in a predefined band. Moreover, unleashing thermal storage flexibility does not always result in comfort loss. This feature makes the DEH households a good candidate for the DR as compared to the DH households when P max is high. Increasing P max enables thermal storages to offer more DR with smaller sacrifice of thermal comfort as compared to the DH loads, which enables thermal storages to contribute more in balancing the wind power. Due to this reason, the power imbalances decrease with increasing levels of P max . Contrarily, at low capacities of P max , DR potential of the DH demand is higher than that of the DEH demand which enables the aggregator to satisfy DH demand mainly with the DWHP system when P max is low. The proportion of the DH demand that is satisfied with the DWHP decreases as P max increases. Such lower operation of the DWHP with increasing P max eventually increases the emissions. It must be emphasized that the maximum expected power imbalances attained for 0% of P max capacity and the maximum emissions obtained for the full P max capacity are still better than that computed in Case I. The relative gain is 11.89% in power imbalances reduction and 21.78% in emission reduction, which are quite noticeable. Increasing P max brings significant reduction in total power imbalances. At 80% capacity, the power imbalances drop to zero with no further rise in emissions beyond this level.
Increasing P max also shifts monetary benefits towards the DEH households as illustrated in Figure 18(b). At 0% capacity in Figure 18(b), no bonus is earned by the DEH households, and the net expected cost is the same as what was obtained in Case I. On the contrary, the DH households capture the maximum monetary benefit at 0% capacity of P max , since only DH demand has flexibility, which results in the maximum corresponding DR potential and leads to a maximum emission reduction as well. At 100% capacity of  P max , the DEH households earn a total bonus of e176.99, which reduces their expected net cost to e137.44. That is a 51.73% reduction, compared to Case I. Additionally, the DH households still manage to earn a bonus of e8.4 at this level of P max , which reduces their expected net cost to e243.68. Due to the share of the DH load, the minimum emission reduction of 21.78% is ensured at any level of P max . The expected total bonus, given by the aggregator to all households, increases from e145.6 to e185.39 as P max rises from zero to the rated capacity. To this end, the combination of the DEH and DH loads in tandem qualifies for both upper and lower level objectives.
So far, we have considered that the thermal storage capacity is just enough to cover 25% of the daily heating demand of a DEH household. Let us examine the whole model with a 50% storage capacity, while all other parameters remain the same as detailed in Table 2. The expected solution for all scenarios obtained in Case I is tabulated in Table 7. Evidently, the expected costs of the DEH households have been decreased as the higher storage capacity results in a higher charging during the lower price period. This point is illustrated in Figure 19 for scenario 5 in order to provide a comparison with Figure 10. Such higher charging in fewer hours reduces the need to charge during most of the horizon which follows higher power imbalances compared to the case when storage capacity was just 25% of the heating demand. Figure 20 displays the outcome with a 50% storage capacity for Case II. In Case II, the higher storage capacity benefits only the DEH households in terms of bonus while the benefits obtained by the aggregator and the DH households remain almost unaffected. As can be seen in Figure 20(b), the DEH households reduce their net energy payment to e94.39 when the full capacity of P max is utilized. It signifies about 65.51% reduction compared to Case I, whereas this relative reduction was limited to 51.73% when storage capacity was 25%.

V. CONCLUSION
This work has proposed a Stackelberg game model to achieve aggregator and households' level objectives. It has been demonstrated that residential heating load has a significant potential in maximizing the wind power utilization and minimizing carbon emissions amid electricity and heat sector coupling. Compared to Case I, when only sector coupling has been performed, the proposed interaction mechanism resulted in about 72.97% improvement in the expected power imbalances and 46.7% in the expected emission reduction. This improvement is accomplished at a cost reduction of 17.78% and 35.95% for the DEH and DH households, respectively. Moreover, it has been shown that the DWHP operation also enables the DH households to interact with the aggregator and offer the DR despite a flat DH tariff. However, increased participation of thermal storage in the DR framework results in more bonuses for the DEH households.
In the future, we aim to quantify the optimal size of thermal storage with the DWHP that is required to account for the days with almost no wind generation or very high wind generation.
ARSLAN AHMAD BASHIR received the bachelor's degree in electrical engineering from the University of Engineering and Technology Lahore, Pakistan, in 2011, and the master's degree in electrical engineering, electrical systems from Aalto University, Finland, in 2016, where he is currently pursuing the Ph.D. degree in electrical engineering. He has worked as the Design Engineer with the Electricity Transmission System Operator of Pakistan from 2012 to 2014 and from 2016 to 2017. His research interests include demand response modeling of thermostatically controlled loads, renewable energy integration, power and heat sector coupling, as well as power system operation, planning, and economics.
ANDREAS LUND received the B.Sc. degree in energy technology from Aalto-University, Espoo, Finland, in 2017, and the M.Sc. degree in electrical engineering from Aalto University, in 2019, where he is currently pursuing the Ph.D. degree in electrical engineering on energy communities with sector coupling to power and heat. Since 2019, he has been an Energy Consultant/Expert at the Granlund Consulting Ltd., on energy efficiency and deep-heat energy systems. His research interest includes renewable energy and how to make heating sector carbon neutral. MATTI LEHTONEN received the master's and Licentiate degrees in electrical engineering from the Helsinki University of Technology, in 1984 and 1989, respectively, and the Doctor of Technology degree from the Tampere University of Technology, in 1992. He was with VTT Energy, Espoo, Finland, from 1987 to 2003. Since 1999, he has been a Professor with the Helsinki University of Technology (now Aalto University), where he is the Head of the Power Systems and High Voltage Engineering. His main activities include power system planning and asset management, power system protection including earth fault problems, harmonic related issues, and applications of information technology in distribution systems. VOLUME 8, 2020