Optimal day-ahead operation of user-level integrated energy system considering dynamic behaviour of heat loads and customers' heat satisfaction

: This study proposes a novel optimal day-ahead operation method for the user-level integrated energy system (IES), which can take into account the dynamic behaviour of heat loads and customers' heat satisfaction. First, based on the energy hub, the dynamic model of heat loads is added into IES. Additionally, the customers' heat satisfaction model is constructed by the variance of water temperature and punishment coefficient. Then, the optimal day-ahead operation model aimed at minimising the total cost is proposed. Constraints of energy purchase, operation of energy converters, loads, temperature and the variation of temperature are all contained. Piecewise linearisation and a standard modelling method are utilised to simplify the initial non-convex and non-linear problem to a quadratic programming problem. Finally, case studies are carried out on a practical calculation example. Results indicate that the total cost can be reduced by 5.6% when considering the dynamic behaviour of heat loads. Also, customers' heat satisfaction exerts a remarkable influence on the total cost.


Nomenclature
Indices D typical day e electricity g gas h heat I set of customer types i customer type M set of energy converters M α set of energy converters whose inputs are α m energy converter max , min maximum and minimum values t time scale α, β, γ energy type αβ energy conversion path from α to β Parameters C heat capacity of residential hot water system C R heat conductance of residential hot water system energy input E α, purchase energy α purchased from the upstream energy network E α, achieve energy α achieved from the upstream energy network L load L α load of energy α T hot water temperature, °C T¯average temperature of hot water, °C T start hot water temperature at the start of a typical day, °C T end hot water temperature at the end of a typical day, °C P heat heating power of residential hot water system P m power of converter m P m, in input power of converter m P m, out output power of converter m r variation rate of hot water temperature λ αm dispatch coefficient of energy α distributed to converter m

Introduction
Nowadays, the conflict between economic development and energy consumption is becoming more and more serious [1,2]. Although the wide use of distributed generations [3], energy storages [4] and electric vehicles [5] has optimised the current energy structure to some degree, scientists are still seeking a novel way to achieve the overall higher energy efficiency and better economy. Benefitted from the development of energy conversion technology, the appearance of integrated energy system (IES) is expected to soft the contradiction more effectively [6]. Compared with the traditional separate energy system, IES is more complex, coupling electricity, gas and heat [7] and bringing the flexibility of energy supply. The processes of energy production, distribution, conversion and consumption are all contained. IES technology has been utilised in many different fields, such as regular energy transmission and distribution system [8][9][10], smart building [11][12][13], sailing [14], flying [15], data centre [16][17][18] and so on. Based on the geographical scope, IES can be further divided into three levels, which are, respectively, cross-region level, region level and user level. Due to the tight integration of different energies, optimising the operation of IES to achieve more economic profit becomes possible. Lots of scholars have made a profound study in this aspect. In Ref. [19], an optimal day-ahead scheduling method for the integrated urban energy system is proposed, which considers the reconfigurable capability of distribution network. In Ref. [20], the economic dispatch model of a combined-heat-and-power (CHP) plant in the wholesale energy markets is constructed to optimise the power production. In Ref. [21], an integrated stochastic optimal operation model is presented. A comprehensive framework for the natural gas transportation network is considered to address the dispatchability of a fleet of fuel-constrained natural gas-fired units. In Ref. [22], a bi-level programming optimisation model is formulated. The upper level represents the coordinated operation, and the lower level simulates the day-ahead market clearing process. Ref. [23] studies the role of hourly demand response in the optimisation of the stochastic day-ahead scheduling of power systems with natural gas transmission constraints. Ref. [24] proposes a coordinated operation strategy for the gaselectricity integrated distribution system, considering alternating current power flow in the power network and the gas hydraulic calculation in gas network. Ref. [25] optimises the conflicting benefits of the electricity network and gas network for daily operation of the IES, while satisfying the operation constraints.
Ref. [26] shows an optimal dispatch strategy for IES with combined cooling heating and power (CCHP) and wind power. Natural gas system is modelled and its security constraints are integrated into the optimal dispatch model.
According to the existing studies, the optimal operation model of IES has been relatively perfect. However, they are mainly based on the steady-state analysis. More practical scenes and details can be further enriched. For instance, loads in the above models are all treated as rigid physical quantities. Any interruption of energy supply will lead to the unavailability. The dynamic behaviour of loads, especially some heat loads such as residential hot water system, is totally ignored. Actually, for these loads, the purpose of heat consumption is to obtain or maintain an acceptable range of temperature for customers. From this side, due to the thermal inertia, some heat loads can be still satisfied even experience a brief interruption of energy supply. This feature will expand the feasible region of the optimal operation model. The initial results of the traditional model may not be the actual global optimal solution. Additionally, the customers' heat satisfaction has not been taken into account. Once the interruption of heat supply is permitted, it will cause the fluctuation of temperature. For customers, different temperature and temperature variation rate mean different qualities of service. This phenomenon will also have a considerable impact on the optimal operation. This paper attempts to address the dynamic behaviour of heat loads and customers' heat satisfaction into the optimal operation model of the user-level IES. First, the user-level IES is modelled as a dual port network based on the energy hub. The dynamic behaviour of heat loads is analysed in the form of differential equations. The model of customers' heat satisfaction is established based on the temperature variance and punishment coefficient. Then, a novel optimal day-ahead operation model is proposed. The objective function is to minimise the total cost, which is composed of energy purchase cost and punishment cost. Constraints of upstream energy network, energy converters, loads, temperature as well as temperature variation rate are all taken into consideration. Due to the existence of differential equations and the appearance of multiplication between optimal variables, discretisation treatment and a standard modelling method are utilised to transform the initial non-convex and non-linear model to a quadratic programming problem. Finally, the comparison between three different scenes demonstrates the correctness and feasibility of the proposed optimal day-ahead operation method.
Rest of the paper is organised as follows. Section 2 introduces the user-level IES model considering the dynamic behaviour of heat loads and customers' heat satisfaction. Section 3 constructs an optimal day-head operation model and discusses its solving method. Case studies are carried out in Section 4 and conclusions are made in Section 5.

Energy hub model
To describe the cooperative operation mechanism of different energy flows in the IES, the energy hub model is introduced in Ref. [27]. Based on the energy hub model, IES can be abstracted as a dual port network containing multiple forms of energy input, energy converters and various forms of energy output. The coupling relationship can be described as follows: The coupling matrix C is composed of the energy conversion efficiency and the dispatch coefficient. Actually, it can be expanded to include demand responses [28], energy storages and renewable energy sources. Different from the cross-region-level and region-level IES, the user-level IES is the lowest layer in the energy transmission and distribution system. It connects the customers and the upstream energy network directly through energy converters. Therefore, the user-level IES is suitable for the end-use customers in a small district such as a residential building or a factory.
Its physical structure is relatively simple. The input can be viewed as energy purchased from the upstream energy network, including upstream distribution network (UDN), upstream gas network (UGN) and upstream heat network. Energy converters play a significant role of converting the input energy to other types. The regular converters include transformer, air condition, gas furnace, electrical boiler, heat exchanger and CHP. Loads are composed of electricity, gas and heat. Most loads are measured in watts while the performance of part heat loads is finally determined by temperature.

Dynamic behaviour of heat loads
Loads in the user-level IES can be divided into two types, rigid and flexible loads. The first type has an exact energy demand. It cannot bear any interruption of energy supply, otherwise the load shedding occurs, such as electricity. For the natural gas system, the diameter of gas pipeline is relatively smaller. Gas stored in the pipeline is not enough to support gas loads when gas supply fails. Hence, the dynamic behaviour of gas system can be also ignored.
The second type of loads is characterised by energy inertia. Even if the energy supply fails, loads can be still satisfied for a while. The purpose of this type of loads is not to achieve a very exact power or energy, but to utilise the energy achieved to satisfy customers' actual demand. In the user-level IES, the most typical example is the residential hot water system where the hot water is consumed for washing, cooking, bathing and so on [29].
In the residential hot water system, customers' energy demand is to achieve or maintain a satisfying range of water temperature. As a continuous quantity, temperature cannot change suddenly. So, even though the heating power becomes zero, the temperature of stored hot water will only decrease slowly. It can still satisfy customers' demand for some time before losing desirable performance. The energy inertia of residential hot water system is shown in Fig. 1. Fig. 1 simulates the trend of temperature variation when the state of energy supply changes. T max and T min are, respectively, the maximum and the minimum temperature that customers can accept. When the energy supply is reliable, the heating power is greater than the heat dissipation. The temperature rises from the initial value as time goes by. By the time that water temperature reaches the maximum value, the input energy supply will be adjusted to maintain the temperature. As the energy supply fails, the temperature gradually decreases until the energy supply works again. In conclusion, no matter the state of energy supply is, the customers' demand can be always viewed as satisfied as long as the temperature is in the accepted range. In order to delve the physical essence of the dynamic behaviour, the linearised energy balance concept is presented [30,31]. The schematic diagram is shown in Fig. 2.
In Fig. 2, the residential hot water system can be seen as a close system aggregated with water tank, thermal insulation layer, cold water supplementary device and the heating component. Water tank is covered by thermal insulation layer to prevent the fall of water temperature. The heating component generates heat to satisfy the hot water loads. Cold water supplementary device guarantees that the tank is always full. Even though the thermal insulation layer can reduce the heat loss, the heat exchange between hot water and surrounding air cannot be totally eliminated. Also, the added cold water absorbs heat and drops the inside water temperature in the mixing process. In a word, apart from the hot water loads, the added cold water and the surrounding air also consume the heat generated by the heating component.
If we assume that the added cold water and initial hot water can be well mixed, the mathematical expression of thermal balance can then be formulated as a differential equation as (2).
The left side represents the energy required for temperature variation at the rate of dT /dt in a unit time interval. The right side consists of three parts. The first part denotes the power of heating component. The second part represents the heat dissipation in the surrounding air. It is the product of the heat conductance and the difference between the hot water temperature and the ambient temperature. The heat conductance can indicate the performance of thermal insulation layer. The better the thermal insulation, the less the heat dissipation. The last part is the energy utilised to heat the added cold water. Due to the assumption of uniformly mixing, it is related to the hot water consumption flow and the temperature difference between output hot water and input cold water. With the above differential equation, temperature at any time can be predicted with the initial temperature and other parameters.

Customers' heat satisfaction model
For energy operators, the dynamic behaviour of heat loads can be utilised to achieve a more economic operation schedule. However, the economic operation based on the dynamic behaviour is usually accompanied with the temperature variation, which affects customers' satisfaction. For customers, service satisfaction is one of the most concerning indexes that must be carefully thought out.
To quantify the influence of temperature vibration on customers' satisfaction, the punishment cost is constructed as (3).
where it is proportional to the variance of water temperature. The punishment coefficient is a constant that depended on customers' consumption habits.

Optimal day-ahead operation model
The aim of optimal day-ahead operation is to achieve the most economic energy purchase scheme. Hence, the objective function is established as (4) to minimise the total cost, which contains the energy purchase cost and the punishment cost.
Detailed constraints are listed as follows taking comprehensive consideration of every link in the user-level IES.

(1) Constraints of upstream energy network:
Constraints of upstream energy network describe the process of energy input. As shown in (5), the input efficiency η input connects the energy purchased and the energy achieved together. Equation (6) defines the time-of-use (TOU) price of each kind of energy. The upper and lower limits of energy purchased and achieved are both included in (7) and (8).

(2) Constraints of energy converters:
Constraints of energy converters are the reflection of energy conversion, covering the energy dispatch and conversion efficiency. Equation (9) represents the upper and lower limits of energy converters' power. In (10) and (11), the dispatch and transformation of energy achieved is embodied by the dispatch coefficient and efficiency. Equation (10) describes the equality constraint of energy conversion. Equation (11) indicates the relationship between the input and output power of energy converters.

(3) Constraints of loads:
Constraints of loads reflect the power and energy balance of energy converters' output and actual loads. For rigid loads, the relationship between converters' output power and loads is shown in (12), which is a strict equality constraint. The flexible heat loads must satisfy the extra constraints from (13) to (16). Equation (13) indicates the thermal balance of the residential hot water system. Equations (14) and (15) are, respectively, the upper and lower limits of water temperature and temperature variation. Equation (16) ensures that the starting state and the final state in a typical day are the same.
T min ≤ T ≤ T max .
T end = T start .

Solution
Due to the thought of dynamic behaviour of heat loads, the temperature constraints in the form of the differential equations increase the difficulty of solving. Besides, the appearance of the multiplication of dispatch coefficient and energy achieved makes the programming harder. The initial optimisation problem is highdimensional, non-linear and non-convex. For differential equations of temperature, piecewise linearisation can be used to simplify the problem. The first step is to transform the heat loads whose unit is kilowatt into residential hot water consumption flow whose unit is litre per second. In the steady state, the temperature variation rate should be zero, and the hot water temperature should be equal to the expected temperature. Make the left side of (13) to be zero and replace T with T d and replace P heating with k f L heat (t). Based on the forecast heat loads L heat (t) of each time interval, the hot water consumption flow at the expected temperature can be calculated by (17)- (19).
Then, take q 0 as a constant value and discretise (13). Constraints of temperature can be simplified as (20) and (21).
After the discretisation, T(t + 1) can be calculated by T(t) and dT(t)/dt even though the hot water temperature is a continuous variable.
For the multiplication of variables in (10), a standard modelling method based on the graph theory can be taken to make the constraint linear [32]. Owing to the single-layer structure of userlevel IES, the energy conversion occurs only one time. All energy converters are connected to customers directly. Hence, the dispatch of energy achieved in the model can be replaced by addition, instead of using the dispatch coefficient.
To summarise, by the simplification above, the initial problem can be transformed into a quadratic programming problem. The final optimal model is as follows: The commercial solver CPLEX is employed to solve the above programming problem.

Case studies
To illustrate the validity of the proposed method, three scenes with different objective functions are compared. Case 1 ignores the dynamic behaviour of heat loads. Case 2 ignores the customers' satisfaction of heat loads. Case 3 establishes the proposed model above, taking comprehensive consideration of dynamic behaviour and customers' heat satisfaction. Case 1 and Case 2 are both linear programming problem, while Case 3 is a quadratic programming problem.
The structure of the user-level IES is shown in Fig. 3. Four different types of energy converters are contained, which are, respectively, electrical boiler, transformer, CHP and gas pipeline. Each converter corresponds to a different energy dispatch coefficient and different energy conversion efficiency.
Parameters of these converters are listed in Table 1, including maximum power, minimum power and energy conversion efficiency.
In Table 1, there are two efficiencies of CHP which, respectively, represents the efficiency of electricity and heat generated by CHP. The heat-electricity ratio is 1.286.
The hourly loads of electricity, heat and natural gas as well as the forecast wind power in a typical winter day in the north of China are shown in Fig. 4 [6].
It can be summed up that the electricity loads and gas loads are higher at noon (12-16 h) and night (19-23 h), and lower in other time. The heat loads are higher at morning and night (1-9, 19-24 h), and lower at noon. The power of wind turbine fluctuates around 45 kW over time.  The TOU prices of electricity are, respectively, 0.3144, 0.6144 and 0.8422 RMB/kWh in different periods. The peak period is 11-12 and 20-21 h. The valley period is 1-8 h. The rest time is normal period. The price of gas is 2.5 RMB/m 3 . The punishment coefficient of all customers' dissatisfaction is 1. The flexible coefficient of heat loads is 10%. The parameters of the residential hot water system are listed in Table 2.
Based on the parameters of residential hot water system above and (13), it can be calculated that when the energy supply fails and the water consumption flow is 0.15 l/s, the water temperature fall will be 5.07°C/h.

Results analysis
The differences among models can be reflected in the total cost, energy purchase schedule and the power of energy converters. The total cost of a typical day in Case 1 is 5508.14 RMB. Meanwhile, for Cases 2 and 3, the total costs are, respectively, 5199.65 and 5315.76 RMB. The decrease of total cost in Cases 2 and 3 shows that the dynamic behaviour of heat loads expands the feasible region and can earn more economic profit. Also, it must be noted that Case 2 brings a decrease of 5.6%, which is quite greater than the decrease of 3.5% in Case 3. It indicates the considerable impact of customers' satisfaction on the total cost. Electricity and gas purchased from upstream energy network function as the energy input of the user-level IES. The optimal energy purchase schemes of different cases are shown in Fig. 5.
In Fig. 5, the general trend of energy purchase schemes is similar. The gas purchased from UGN is far more than the electricity purchased from UDN. It is because that CHP based on micro-turbine has a better economy performance. Based on the low heating value, the unit price of gas can be calculated as 0.19 RMB/ kWh. Considering the loss in the energy conversion process, the cost of electricity generated by CHP is 0.5428 RMB/kWh, which is still lower than the electricity price in normal and peak period. The cost of heat generated by CHP is 0.4222 RMB/kWh. Therefore, gas purchase should have a higher priority over electricity purchase.
Taking Case 1 as an example to analyse in detail, the electricity purchase mainly occurs in 1-8, 12-16 and 18-24 h. This phenomenon is generated by the constraint of the heat-electricity ratio. Although gas purchase is more economic, it must match the electricity or heat loads. The ratio of electricity loads and heat loads during these periods fluctuates from 1.318 to 4.161, which is larger than the determined heat-electricity ratio 1.286. It means that CHP cannot satisfy all the loads independently. IES still needs the support of UDN. In other time, nearly all the electricity loads can be satisfied by gas without any electricity purchase.
The difference among the three cases mainly lies in the period of 12-16 and 18-24 h. Cases 2 and 3 correspond to more gas purchase in this time. It is the optimal results after considering dynamic behaviour of heat loads. From observation, heat loads in 18-24 h are heavy and the gas purchase has arrived the upper limit. In 18-24 h, the time-of-use price of electricity is in the normal or peak period. Hence, preheating the residential hot water by utilising the thermal inertia is a potential way to cut the peak of heat loads. More gas purchase in 12-16 h can preheat the water and reduce the electricity purchase cost in 18-24 h. When the load peak arrives, the residential hot water system can even work without any energy supply because of the preheating. The electricity purchase in 18-24 h can be partly saved.The power of each energy converter supporting heat loads in different cases is shown in Fig. 6. Fig. 6 shows the energy allocation schemes between energy converters after the energy purchase. The heat loads are supplied by both electric boiler and CHP. In 11-16 h, the ratio of heat loads and electricity loads waves from 0.571 to 1.261, which is lower than the heat-electricity ratio of CHP. During this period, the heat loads can be fully met by CHP in all three cases. In Cases 2 and 3, because of the consideration of dynamic behaviour of heat loads, the sum of the power of CHP and electric boiler no longer equals to the product of the forecast heat loads and flexible coefficient  Table 2 Parameters of residential hot water system Volume, m 3 exactly. The heating power of residential hot water system becomes one of the optimal variables, affecting the energy purchase cost and total cost by preheating. Water can be preheated before the peak of loads arrives so long as the temperature is in the acceptable range. It is the embodiment of dynamic characteristics of flexible heat loads, which is different from electricity loads and gas loads.The relationship between heating power and residential hot water temperature at different time is shown in Fig. 7.
In Case 1, the heating power equals to the product of the flexible coefficient and the forecast heat loads. The heat loads are treated as rigid loads, and the failure of heating power is not permitted. The hot water temperature is the steadiest, maintained at 65°C without any wave.
In Case 2, when the dynamic behaviour is considered, the temperature fluctuates from 50 to 90°C, instead of maintaining at the expected temperature. Also, it must be pointed out that the heating power in 1-2 and 18-23 h is zero. In 1-4 h, the heating power is lower than the heat dissipation power. So, the temperature decreases as time goes by. In 5-11 h, the temperature is basically around 50°C, which indicates the energy balance between heating and dissipation. In 12-16 h, the temperature rises at the greatest acceptable variation rate. After a brief stability, temperature decreases again. The start temperature and end temperature are 65°C.
In Case 3, the dynamic behaviour of heat loads can also be observed in the curve. What's more, compared with Case 2, the range of temperature variation is quite narrower. It is because that customers' heat satisfaction is related to the variance of temperature. As part of the objective function, the bigger the variance of temperature, the more punishment cost. Hence, the temperature curve of Case 3 is relatively smoother.

Sensitivity analysis
This subsection is set to delve the influence factor of the total cost.
To study the impact of heat conductance and flexible coefficient on the total cost, the sensitivity analysis based on Case 2 is employed. Fig. 8 illustrates that compared with the flexible coefficient of heat loads, the change of heat conductance has little affection on the total cost, which even can be ignored. When the heat conductance is determined, the total cost grows inverse proportion to the flexible coefficient. The bigger coefficient means that more hot water can be preheated. As a result, more energy consumption can be reduced in the peak period. Nevertheless, when the flexible coefficient is beyond 10%, its saving effect is not as remarkable as before. Delving (19), we can find that the flexible coefficient determines the actual hot water consumption flow. High coefficient leads to a huge water consumption which is a big challenge to the initial hot water stored. Thus, most power of the heating component is dispatched to heat the added cold water. The residential hot water system's ability of cutting peak loads is weaken, and the total cost stabilises gradually.
To compare the impact of punishment coefficient and flexible coefficient in Case 3, the sensitivity analysis is presented in Fig. 9.
In Fig. 9, the minimum value of total cost is achieved when the flexible coefficient is biggest and the punishment coefficient is smallest. If we fix the punishment coefficient and continuously increase the flexible coefficient to 80%, the total cost will at first decrease and then will stabilise at around 5316 RMB. In another side, when the flexible coefficient is a determined value, higher punishment coefficient always contributes to a greater total cost. Furthermore, with the increase of flexible coefficient, the impact of punishment coefficient on the total cost is more remarkable.
According to the above-mentioned analysis, both flexible coefficient and punishment coefficient have significant impact on the total cost, but the volume of water tank in the residential hot water system is a fixed value, 8 m 3 . To prove its influence, the relationship between the total cost and water tank volume is shown in Fig. 10.
From Fig. 10, we can conclude that with the water tank volume getting larger, the value of total cost is dropping and finally stabilises. Benefitted from the energy inertia, the residential hot water system can function as an energy storage. It absorbs energy to preheat the stored water in the valley period and reduces the peak loads. The volume is related to the capacity of energy storage. The increase of water tank volume can strength the ability of cutting peak loads. The corresponding total cost is lower. Howbeit, the ability cannot be enhanced unboundedly. When the volume is beyond 40 m 3 , the cost is levelling at around 5100 RMB.

Conclusion
This paper proposes an optimal day-ahead operation method which is suitable for the user-level IES with comprehensive consideration of dynamic behaviour and customers' heat satisfaction. The impacts of heat conductance, flexible coefficient, punishment coefficient and water tank volume on the total cost are, respectively, analysed. Case studies are carried out on three different scenes. The simulation results show that: i. The energy inertia of flexible heat loads provides the abilities of energy storage and peak load cutting. ii. The optimal day-ahead operation considering the dynamic behaviour and customers' heat satisfaction can reduce the energy purchase cost by preheating the flexible heat loads. iii. The flexible coefficient, punishment coefficient and water tank volume are the most remarkable factors impacting the total cost.

Acknowledgments
This work was supported by the National Key R&D Program of China (2018YFB0905000), Shanghai Sailing Program (18YF1411600) and Science and Technology Project of State Grid Corporation of China (5102-201956308A-0-0-00).