Optimal Operation Strategy of a Community Integrated Energy System Constrained by the Seasonal Balance of Ground Source Heat Pumps

: Ground source heat pumps (GSHPs) are now widely used in community integrated energy systems (CIES) because of their high e ﬃ ciency in energy conversion. However, the coe ﬃ cient of performance (COP) of GSHPs is unstable if the soil temperature changes with seasonal imbalanced cooling and heating loads, thus downgrading the overall performance of CIES. In this paper, an annual optimization model for CIES that considers the seasonal balance of GSHPs is established. Then, a day-ahead operation strategy based on the pre-allocated load of the GSHP in the yearly balance is proposed while considering the uncertainties in daily conditions. The proposed strategy is validated on a practical CIES in China and assessed on a year-round time scale. The results show that the operation of CIES can be stable, economical and sustainable while ensuring the seasonal balance of the


Introduction
With the wide integration of renewable energy resources, various energy forms have been coupled through interconnection technologies [1]. Many countries around the world have conducted research on the influence of renewable energy on traditional power systems [2,3]. Challenges have been put forward for the actual operation of multiple energy systems, such as the coordination and conversion of diverse energy sources [4,5]. A community integrated energy system (CIES) is regarded as a potential solution to the challenges on the consumer side, which will plan, operate and manage different energy subsystems under a unified framework [6]. As an important part of CIES, a ground source heat pump (GSHP) utilizes the relatively stable temperature of underground soil or groundwater to achieve energy transfer [7]. The mechanism of heat exchange through geothermal energy enables the GSHP to offer higher conversion efficiency, fewer carbon emissions [8], lower operational and maintenance costs and a longer service life compared to ordinary air conditioners and boilers. Entchev et al. [9] compared the performance of a GSHP system and a conventional system that utilizes boiler and chiller to supply the thermal loads of a detached house and a small office building. By consuming a small amount of electricity, the GSHP can replace the electric boiler to obtain the heat stored under the ground and supply the heating load in winter. In the summer, it can also transfer the excess heat from the building to the ground instead of air conditioning. At the same time, the GSHP can also supply domestic hot water [10]. Hassam et al. [11] presented a design for a solar thermal system operating with an Sustainability 2020, 12 integrated GSHP and a seasonal borehole storage to provide domestic hot water (DHW) and space heating (SH) for community-sized demand. In most regions, the GSHP experiences cooling, heating and non-operational periods throughout the year, as shown in Figure 1. The energy flow directions of the GSHP for cooling and heating are opposed, thus, the soil temperature will rise or fall accordingly [12]. Ji et al. [13] integrated ground temperature distribution with the coefficient of performance (COP) of GSHPs and built a model that can accurately and quickly predict the long-term ground temperature distribution. Ideally, the soil will naturally maintain the relative balance of energy reception and divergence, which can be regarded as a huge dynamic energy balance system [14]. During the cooling/heating period, heat pumps supply high-temperature/low-temperature water to buildings on demand through the supply pipes. After the exchange, the water is fed back to the heat pumps in the energy station through the return pipes. Under ideal conditions, the effect of the underground heat exchanger on the soil temperature during the cooling and heating periods can be recovered via natural heat transfer with the soil during non-operational periods. After the annual dispatch, the soil temperature is expected to be maintained at its initial value.
Sustainability 2020, 12, x FOR PEER REVIEW 2 of 24 thermal system operating with an integrated GSHP and a seasonal borehole storage to provide domestic hot water (DHW) and space heating (SH) for community-sized demand. In most regions, the GSHP experiences cooling, heating and non-operational periods throughout the year, as shown in Figure 1. The energy flow directions of the GSHP for cooling and heating are opposed, thus, the soil temperature will rise or fall accordingly [12]. Ji et al. [13] integrated ground temperature distribution with the coefficient of performance (COP) of GSHPs and built a model that can accurately and quickly predict the long-term ground temperature distribution. Ideally, the soil will naturally maintain the relative balance of energy reception and divergence, which can be regarded as a huge dynamic energy balance system [14]. During the cooling/heating period, heat pumps supply high-temperature/low-temperature water to buildings on demand through the supply pipes. After the exchange, the water is fed back to the heat pumps in the energy station through the return pipes. Under ideal conditions, the effect of the underground heat exchanger on the soil temperature during the cooling and heating periods can be recovered via natural heat transfer with the soil during non-operational periods. After the annual dispatch, the soil temperature is expected to be maintained at its initial value.  However, there are often obvious differences between the cooling and heating loads in actual operation due to the climates in different regions. As shown in Figure 2, the regions can be divided into three types according to the length of the cooling/heating period: normal regions, hot regions and cold regions. There may be significant differences in the load demand between these regions, and the soil temperature changes caused by the operation of the GSHP will also be very different. As the three conditions on the left side of Figure 2 show, it is possible to keep the thermal balance of soil and maintain reliable cooling and heating over long periods of operation under ideal conditions. However, in most cases, the difference between the heat extracted from, and released into, the soil can be relatively large, and it is difficult to balance the heat exchange in the soil throughout the year. Ref. [15] established a theoretical ground heat exchange model, and Beijing, Harbin and Zhengzhou were used in numerical simulations as three typical cities (a cold region, a severely cold region and a hot summer and cold winter region, respectively). The simulation results showed that the imbalance efficiencies and annual average soil temperature decreases were not the same in the different regions. Therefore, GSHPs are rarely adopted in areas that have extraordinarily imbalanced heating and cooling loads [16]. However, there are often obvious differences between the cooling and heating loads in actual operation due to the climates in different regions. As shown in Figure 2, the regions can be divided into three types according to the length of the cooling/heating period: normal regions, hot regions and cold regions. There may be significant differences in the load demand between these regions, and the soil temperature changes caused by the operation of the GSHP will also be very different. As the three conditions on the left side of Figure 2 show, it is possible to keep the thermal balance of soil and maintain reliable cooling and heating over long periods of operation under ideal conditions. However, in most cases, the difference between the heat extracted from, and released into, the soil can be relatively large, and it is difficult to balance the heat exchange in the soil throughout the year. Ref. [15] established a theoretical ground heat exchange model, and Beijing, Harbin and Zhengzhou were used in numerical simulations as three typical cities (a cold region, a severely cold region and a hot summer and cold winter region, respectively). The simulation results showed that the imbalance efficiencies and annual average soil temperature decreases were not the same in the different regions. Many previous studies have focused on the analysis of the cooling and heating load balance of a GSHP and its influence on the COP. Choi et al. [21] discussed possible future problems that may be caused by the long-term operation of a GSHP under load imbalance conditions, such as thermal buildup or depletion and performance degradation. Ref. [22] studied the long-term performance of a GSHP system for an office building in Suzhou, which is a hot summer and cold winter climate region of China. The results proved that the GSHP system produces ground heat accumulation, which reduces the system's COP and increases energy consumption. In [23], the influences of the underground thermal imbalance rate, soil type, borehole layout and groundwater advection on the underground soil temperature distribution were analyzed. According to the results, the better the soil thermal conductivity and diffusivity, the faster the soil heat diffusion. Refs. [24,25] analyzed the changes in soil temperature due to the imbalance of GSHP cooling and heating loads. Considering the building environment in the UK, ref. [26] showed that during the application of a GSHP for building energy supply, the COP of the GSHP decreased or even failed due to the imbalance of cooling and heating loads. Cui et al. [27] explored the effects of the soil thermal properties on its temperature, and the COPs of the GSHP are investigated and compared between the first year and tenth year of operations. Miglani et al. [28] concluded that the ground temperature will continue to fall due to unbalanced cooling and heating loads and suggested the following solutions: Retrofitting the building's facade, roofs and floors can reduce CO2 emissions and the cooling demand. Similarly, GSHPs equipped with heat recovery are also an efficient approach compared to conventional systems. This type of system can reject less heat into the soil and extract more heat from it, which helps reduce the soil's thermal imbalance [29]. However, these methods produce additional investment costs.
In order to maintain a constant soil temperature and ensure the stable operation of the GSHP unit, the GSHP usually works in a seasonal intermittent operation mode (continuous operation in winter and summer) [30]. However, in most cases, seasonal balance cannot be simply achieved even with this operation mode. Therefore, it is necessary to limit the working time of the GSHP or adopt an intermittent operation mode with auxiliary equipment [31,32]. In this way, the temperature of the soil can gradually recover its initial temperature during non-operational periods [33,34], and the heat The long-term imbalanced operation of a GSHP will change the soil temperature around the underground heat exchanger [17]. The closer the soil is to the heat exchanger, the more obvious the temperature change is [18]. The three conditions on the right reflect the changes in soil temperature that may occur when the cooling and heating power of the GSHP cannot reach a balance. When the difference between the temperature of the soil and the return pipe is too small, the heat transfer efficiency will significantly decrease, which means a reduction in the COP of the GSHP unit, which may even cause a shutdown of the GSHP units [19,20].
Many previous studies have focused on the analysis of the cooling and heating load balance of a GSHP and its influence on the COP. Choi et al. [21] discussed possible future problems that may be caused by the long-term operation of GSHPs under load imbalance conditions, such as thermal buildup or depletion and performance degradation. Ref. [22] studied the long-term performance of a GSHP system for an office building in Suzhou, which is a hot summer and cold winter climate region of China. The results proved that the GSHP system produces ground heat accumulation, which reduces the system's COP and increases energy consumption. In [23], the influences of the underground thermal imbalance rate, soil type, borehole layout and groundwater advection on the underground soil temperature distribution were analyzed. According to the results, the better the soil thermal conductivity and diffusivity, the faster the soil heat diffusion. Refs. [24,25] analyzed the changes in soil temperature due to the imbalance of GSHP cooling and heating loads. Considering the building environment in the UK, ref. [26] showed that during the application of GSHPs for building energy supply, the COP of the GSHP decreased or even failed due to the imbalance of cooling and heating loads. Cui et al. [27] explored the effects of the soil thermal properties on its temperature, and the COPs of the GSHP are investigated and compared between the first year and tenth year of operations. Miglani et al. [28] concluded that the ground temperature will continue to fall due to unbalanced cooling and heating loads and suggested the following solutions: Retrofitting the building's facade, roofs and floors can reduce CO 2 emissions and the cooling demand. Similarly, GSHPs equipped with heat recovery are also an efficient approach compared to conventional systems. This type of system can reject less heat into the soil and extract more heat from it, which helps reduce the soil's thermal imbalance [29]. However, these methods produce additional investment costs.
In order to maintain a constant soil temperature and ensure the stable operation of the GSHP unit, the GSHP usually works in a seasonal intermittent operation mode (continuous operation in winter and summer) [30]. However, in most cases, seasonal balance cannot be simply achieved even with this operation mode. Therefore, it is necessary to limit the working time of the GSHP or adopt an Sustainability 2020, 12, 4627 4 of 24 intermittent operation mode with auxiliary equipment [31,32]. In this way, the temperature of the soil can gradually recover its initial temperature during non-operational periods [33,34], and the heat transfer capacity of the underground heat exchanger and the operating efficiency of the GSHP can be guaranteed.
A constant soil temperature is easier to achieve under reasonable operation dispatching, so the selection of auxiliary equipment in the intermittent operation mode is important. Many studies have been conducted on a hybrid ground source heat pump (HGSHP) system that adds auxiliary equipment to the GSHP. Hou et al. [35] proved, by long-term simulation results, that the HGSHP has a better performance compared to a conventional GSHP system. Ref. [36] presented a review of the research and applications of the GSHP integrated with a thermal energy storage (TES) system, including: ice storage tank, solar collector, soil, water tank and phase change materials (PCM). Ma et al. [37] proposed a new distributed energy system that integrates cogeneration, photovoltaic energy and GSHPs for a large office building in Beijing, which is an initial form of integrated energy system. A HGSHP with auxiliary heat source was proposed in [38] to solve the problem of the unbalanced loads of most buildings in the cold climate zone. Refs. [39,40] focused on the performance of building-integrated solar energy systems. These studies examined the coupling between GSHPs and various types of energy equipment, but the effect of the GSHPs' operation on soil temperature has not been considered in detail. The obtained results are likely more optimistic than the actual performance.
From the perspective of short-term operation, various types of auxiliary equipment can shorten the continuous operation time of GSHPs, which is beneficial for the recovery of soil temperature. However, the best manner to maintain a stable soil temperature under long-term operation is to ensure a balanced energy supply between the cooling and heating periods. Therefore, in the long-term operation of the CIES, the seasonal balance of a GSHP under the year-round time scale has to be considered. At the same time, on the premise of ensuring the stable operation of the CIES [41], the operation cost needs to be optimized through the coupling and complementation of multiple devices [42]. This paper proposes an optimization model for the optimal operation of the CIES considering the seasonal balance of the GSHP on a year-round time scale. Moreover, in order to reduce the impact of day-ahead uncertainties caused by load demand and environmental factors in practice, a day-ahead scheduling strategy considering the uncertainties in daily conditions is presented. Compared with traditional CIES optimization methods, the application of the scheduling strategy proposed in this paper can realize a stable, economic, and sustainable operation of the CIES in different types of regions in order to ensuring the GSHP's seasonal balance. The main contributions of this paper are summarized as follows: (1) A year-round optimal scheduling model for the CIES that considers the seasonal balance of the GSHP is established. The GSHP works in seasonal intermittent operation mode and cooperates with various types of energy supply and energy storage equipment in the CIES. By optimizing the allocation of the cooling and heating loads among the GSHP and the other equipment, an optimal operation plan can be developed. This plan can guarantee the balance of the total cooling and heating supply of the GSHP throughout the year and can realize a long-term stable operation.
(2) Control variables are introduced to build a CIES day-ahead optimal scheduling strategy considering the daily uncertainties of demand. To lower the effect of the uncertainties from renewable energy generation and loads in daily operation, the scheduling scheme based on a pre-allocated working plan is adjusted according to the day-ahead forecasting of environmental factors and load demands. By introducing uncertainty constraints, the daily operation plan will be adjusted based on the forecasted load and the output of distributed generation (DG), while ensuring the annual seasonal balance.
The remainder of this paper is organized as follows. Section 2 builds the annual optimization model for the CIES considering the seasonal balance of the GSHP. Section 3 presents the corresponding day-ahead optimal operation strategy considering uncertain factors. Case studies are given in Section 4 to verify the effectiveness of the proposed optimization methods, and conclusions are provided in Section 5.

Structure of the Studied CIES
To solve the problem of unbalanced cooling and heating energy demand during the operation of the CIES on a year-round time scale, this paper proposes an optimal operation model that considers the seasonal balance of the GSHP. We also select a practical CIES in Tianjin [43] as the case study. As shown in Figure 3, the CIES consists of multiple energy subsystems: the regenerative electric boiler (EB) subsystem, GSHP subsystem, conventional water-cooled chiller (CWC) subsystem and ice-storage (IS) subsystem. The electric load of the CIES is satisfied by external power grid and photovoltaic (PV) energy. A variety of energy supply and energy storage devices cooperate with the GSHPs to meet the cooling/heating loads of the CIES and keep the temperature of buildings within a suitable and comfortable range [44]. The GSHP subsystem consists of three GSHP units, and GSHPs cooperate with a cold water tank during the cooling period. The water-cooled chiller subsystem is equipped with two CWCs. The IS consists of two double-duty chillers (DCs), two chilled water pumps (CWPs), two ethylene glycol pumps and an ice-storage tank (IT). The electric boiler subsystem is equipped with four electric boilers and three hot water tanks.

Structure of the Studied CIES
To solve the problem of unbalanced cooling and heating energy demand during the operation of the CIES on a year-round time scale, this paper proposes an optimal operation model that considers the seasonal balance of the GSHP. We also select a practical CIES in Tianjin [43] as the case study. As shown in Figure 3, the CIES consists of multiple energy subsystems: the regenerative electric boiler (EB) subsystem, GSHP subsystem, conventional water-cooled chiller (CWC) subsystem and icestorage (IS) subsystem. The electric load of the CIES is satisfied by external power grid and photovoltaic (PV) energy. A variety of energy supply and energy storage devices cooperate with the GSHPs to meet the cooling/heating loads of the CIES and keep the temperature of buildings within a suitable and comfortable range [44]. The GSHP subsystem consists of three GSHP units, and GSHPs cooperate with a cold water tank during the cooling period. The water-cooled chiller subsystem is equipped with two CWCs. The IS consists of two double-duty chillers (DCs), two chilled water pumps (CWPs), two ethylene glycol pumps and an ice-storage tank (IT). The electric boiler subsystem is equipped with four electric boilers and three hot water tanks.

Annual Optimization Model for the CIES Considering Seasonal Balance
During the cooling and heating periods, the GSHPs operate in two modes: cooling and heating. As shown in the following constraints, the two modes cannot operate simultaneously: where , represents the ON/OFF state of the n-th GSHP unit at time t, , , and , , represent the cooling/heating mode of the n-th GSHP at time t, is the number of GSHPs. The operation models of the GSHP subsystem and related equipment are described as follows, and the models of the other subsystems in the CIES are summarized as Equations (A1)-(A25) in Appendix A.
(a) Modeling of the cooling period When the GSHP subsystem operates in the cooling season, there are two modes: refrigeration and cooling-storage. Equation (4) indicates that GSHPs can only run in one mode at a time. Equations

Annual Optimization Model for the CIES Considering Seasonal Balance
During the cooling and heating periods, the GSHPs operate in two modes: cooling and heating. As shown in the following constraints, the two modes cannot operate simultaneously: where U HP t,n represents the ON/OFF state of the n-th GSHP unit at time t, U HP,C t,n and U HP,H t,n represent the cooling/heating mode of the n-th GSHP at time t, N HP is the number of GSHPs. The operation models of the GSHP subsystem and related equipment are described as follows, and the models of the other subsystems in the CIES are summarized as Equations (A1)-(A25) in Appendix A.
(a) Modeling of the cooling period When the GSHP subsystem operates in the cooling season, there are two modes: refrigeration and cooling-storage. Equation (4) indicates that GSHPs can only run in one mode at a time. Equations (5)- (7) Sustainability 2020, 12, 4627 6 of 24 describe the relationship between the refrigeration and cooling-storage modes and ON/OFF states of the GSHPs. Equations (8) and (9) restrict the startup/shutdown sequence within GSHP units.
where U HP,R t,n and U HP,S t,n represent the refrigeration mode and cooling-storage mode of the n-th GSHP unit at time t.
Equation (10) describes the cooling power of the n-th GSHP unit at time t, which is the sum of refrigeration and cooling-storage power. Equations (11) and (12) describe the lower/upper limits of the refrigeration and cooling-storage power of the GSHPs under two modes.
where Q HP,R t,n and Q HP,S t,n are the refrigeration and cooling-storage power of the n-th GSHP unit at time t, Q HP,R n,min , Q HP,R n,max , Q HP,S n,min and Q HP,S n,max are lower/upper refrigeration and cooling-storage limits of the n-th GSHP unit.
The GSHPs can operate under the cooling-storage mode during the valley price period to store cooling energy, and then release the stored energy for space cooling during the peak price period. This load shift can not only reduce the cost of purchasing electricity, but also relieve the pressure of peak load [45].
As shown in Equations (13) and (14), when the GSHPs are in the cooling-storage or refrigeration mode, the cold water tank cannot be discharged. Equations (15) and (16) describe the relationship between the cooling mode of the cold water tank and the cooling state of the CWPs.
represents the ON/OFF state of the n-th CWP at time t, N WT,CWP is the number of CWPs.
As shown in Equation (17), the cooling power of the cold water tank at time t can be calculated by the rated flow of the CWPs collateral to the cold water tank and the total flow of the primary CWPs. Equation (18) indicates that the number of CWPs collateral to the cold water tank in the cooling state should meet the demands of the cold water tank.
where F WT,CWP is the rated flow of the CWPs collateral to the cold water tank, L C t is the cooling load of the CIES at time t, F C t represents the total flow of the GSHP's CWP at time t, Q WT,CWP max is the maximum cooling power.
The cold storage W WT,C t of the cold water tank at time t can be obtained from the storage of the last scheduling period t − 1, the cooling-storage power of the GSHPs and the cooling power of the cold water tank. At the same time, the cold storage of the cold water tank should meet the constraints of Equation (20): where ε WT,C is the heat loss rate of the cold water tank, ∆t represents the scheduling interval, W WT,C max represents the upper limit of the energy stored in the cold water tank. The consumed power of the cold water tank P WT,C t at time t includes the electricity demand of the interlocked circulating pumps on both sides of the plate heat exchanger during refrigeration.
The total power consumption of the GSHP subsystem in the cooling season can be obtained by the following constraint: t,n P HP,CP + P HP,CWP +U HP,S t,n P WT,CP1 + P WT,CP2 + P HP,CP where COP HP,C represents the cooling COP of the GSHP unit, P HP,CWP and P HP,CP are the rated power of the CWPs and the cooling water pump interlocked with the GSHPs.
(b) Modeling of the heating period.
The GSHP units provide heating power directly to the CIES during the heating period. The ON/OFF state and power constraints are constrained as follows: where Q HP,H t,n represents the heating power of the n-th GSHP unit at time t, F HP t is the total flow of the CWPs at time t, L H t is the heating load of the CIES at time t, F H t is the rated flow of the CWP of the GSHP at time t, Q HP,H n,min and Q HP,H n,max represent the lower/upper heating limits of the n-th GSHP unit. In the heating season, the total power consumption of the GSHP subsystem can be obtained by: where COP HP,H represents the heating COP of the GSHP unit. The operation of the GSHP couples with various energy supply and energy storage devices within the CIES. When these devices are operated in cooperation, the dispatching scheme of the GSHP will be more complex due to the energy coupling. By optimizing the power allocation scheme of the CIES, the most economical heating and cooling power of the GSHP in a year can be calculated, and the seasonal balance of the GSHP can be realized. If the yearly optimization of the CIES is not considered, and the heating power supply is determined directly by the cooling power, or vice versa, it will cause higher operation costs, and might even fail to be realized in actual operation. Therefore, it is necessary to consider the device operation and power balance constraints in the CIES to establish an optimization model as well as satisfying the seasonal balance of the GSHP.
In the annual optimization model, the objective function is to minimize the annual operation cost of the CIES.
where N day represents the total days of cooling and heating period in a year, N T is the number of scheduling intervals in a day, λ t is the purchase price of electricity at time t, P EG t represents the power purchasing from the external power grid at time t, ∆t represents the scheduling interval.
During the CIES operation, the tie-line power P EG t between the CIES and the external power grid needs to be limited by maximum power constraints: where P EG max is the maximum allowed power of the tie-line. Equations (29)-(31) are the load balance constraints of cooling, heating, electricity supply and demand that ensure the stable operation of the CIES.
where U CWC t,n and U EB t,n represent the ON/OFF state of the n-th CWC and EB at time t, N CWC is the number of CWCs, Q CWC t,n and P CWC t are the cooling power and consumed power of CWCs at time t, Q IS t and P IS t are the cooling power and consumed power of IS at time t, Q WT,H t is the heating power of the hot water tanks at time t, N EB is the number of EBs, Q EB t,n and P EB t are the heating power and consumed power of EBs at time t, P PV t is the power output of the PV system, L E t is the electric load of the CIES at time t, P WT,H t is the consumed power of the hot water tanks at time t, P AWP t is the consumed power of air conditioning hot water pump at time t.
Equation (32) guarantees that the total cooling and heating power of the GSHP throughout the year will be balanced: The annual optimization model of the CIES operation considering the GSHP's seasonal balance can be written as where the constraints include equipment operational constraints, Equations (1)- (26) and (A1)-(A25); the power balance constraints, Equations (28)- (31); and the seasonal balance constraint of GSHPs, Equation (32).
This model is a mixed-integer nonlinear programming (MINLP) model, which can be transformed into a mixed-integer linear programming (MILP) model [46]. Through linearization, this model can be solved by the YALMIP optimization toolbox. According to the forecasting climate and load data of the CIES, the following results can be obtained by solving the model on a year-round time scale: load pre-allocated plan of various equipment in the CIES throughout the year, ON/OFF state of each device at each moment, total operation cost of the system, total optimal cooling and heating energy of GSHPs' Q ann , and pre-allocated cooling/heating energy of GSHPs on the i-th scheduling day Q P,C i , Q P,H i , etc.

Day-Ahead Optimal Operation Strategy
For an overall optimization throughout the year, the input data are the predicted load value and climate data obtained based on the historical data. Since the forecasting data on a year-round time scale may have a large deviation from the day-ahead forecasting data, using day-ahead scheduling strictly according to the pre-allocated scheme obtained by the annual optimization may not be appropriate. Therefore, it is necessary to allow a deviation between the day-ahead operation scheme of the GSHPs and the pre-allocated cooling/heating energy during annual optimization. Thus, we propose a day-ahead optimal operation strategy for the CIES considering the uncertainties of the loads. By introducing control variables, the energy supply of the GSHPs can be adjusted flexibly around the pre-allocated value on a daily time scale according to the day-ahead load and climate information.
In the process of daily scheduling around the year, the range of uncertainty will be flexibly tightened or relaxed. In this way, the reasonable allocation of daily load in the CIES can be achieved while ensuring the seasonal balance of the GSHPs. This strategy can be expressed as follows.
(1) Establish the day-ahead optimization model: The objective function is to minimize the daily operation cost of the CIES.
The equipment operation constraints and supply and demand balance constraints are the same as the annual optimization model. However, the optimal scheduling in actual operation may deviate from the pre-allocated values due to inaccurate forecasting in a year-round time scale. In order to ensure the seasonal balance and obtain the optimal day-ahead scheduling plan, it is necessary to incorporate the output constraints of the GSHP that take into account the uncertainties in daily conditions: (2) Solve the optimization model (Equation (38)) for the current scheduling day: Read in the day-ahead cooling/heating/electric load demands and climate data of the i-th day. Then, the optimal scheduling model can be solved after linearizing the model into an MILP formulation.
(3) Update the control variables: In day-ahead scheduling, the loads will likely be significantly different from the annual forecasting data due to uncertainties. If this happens for many scheduling days, the seasonal balance will be difficult to achieve as planned. In a year-round operation, the deviation of energy supply on a certain day is not a problem. However, it is necessary to ensure that the overall energy supply is generally consistent with the optimal value of pre-allocation. Therefore, the cooling/heating energy in day-ahead scheduling is limited by the control variable based on the pre-allocated scheme. If the overall energy supply deviates significantly from the pre-allocated value, the control variable will be adjusted constantly to limit the output of the GSHPs, until the annual cooling and heating loads reach equilibrium.
Taking the heating period as an example, the control variables and the day-ahead limited energy need to be continuously updated considering the uncertainties, as shown in the following equations: are the actual and pre-allocated heating energy on the i-th scheduling day, ε is the maximum allowable deviation rate between the pre-allocated and actual values. The control variable is adjusted according to the deviation to tighten or relax the constraints. The closer the control variable is to zero, the tighter the constraints on the energy supply of the GSHPs.
(4) Determine whether the optimization of all scheduling days is completed. If it is, skip to Step 6; otherwise, i = i + 1.
(5) Determine whether the total energy supply has reached the planned value: Equation (45) calculates the difference between the day-ahead total cooling/heating energy up to the (i − 1)th scheduling day and the annual total pre-allocated energy of the GSHPs.
During operation, it is necessary to confirm whether the total current energy supply has reached the annual planned value. If it has, the GSHPs will no longer be in operation. Determine whether the difference is less than the pre-allocated energy output of the i-th scheduling day:  (6) Determine whether the annual cooling/heating energy is less than the total value of pre-allocation: Q TA i < Q ann . If it is, the day-ahead scheduling optimization ends; otherwise, the GSHPs continue to operate until the actual output is equal to the planned value.
If the planned total cooling/heating value has not been reached by the end of the annual scheduling, the GSHPs will continue to provide cooling/heating energy under no-load condition on non-scheduled days until the scheme is completed.
The optimal scheduling strategy for the CIES considering the seasonal balance of the GSHPs and the uncertainties in daily conditions can be organized into the flow chart shown in Figure 4. By adopting the operation strategy mentioned above, the daily energy output of the GSHPs can be adjusted flexibly according to the day-ahead conditions while guaranteeing the total energy supply throughout the year. The closer the values of the parameters (cooling and heating load, weather factors, etc.) in daily operation to the historical-based forecasts, the smaller the deviation between the day-ahead optimization results and the pre-allocated values of the scheduling day throughout the year.

Case Study and Results
The case study selected was a typical CIES in Tianjin, China. The energy demands in the CIES include cooling, heating and electricity. Renewable energy, such as solar energy and geothermal energy, is utilized according to the actual structure and equipment. The electricity demand in the CIES is satisfied by the external power grid and the PV system together. The cooling demand of the CIES is met by GSHPs and other types of cooling and cooling-storage equipment during the cooling period. The heating load during the heating period is supplied by GSHPs, EBs and hot water tanks. The GSHPs work in intermittent operation mode, including seasonal intermittent operation (continuous operation in winter and summer) and short-term intermittent operation (configuring auxiliary equipment).
Taking the heating period as an example, the heating energy is provided by the GSHPs during the valley price period. During this period, the regenerative electric boiler subsystem stores hot water. After that, during the peak price period, the hot water tanks release the stored energy to meet the heating demand, and the vacancy can be replenished by the EBs or GSHPs.
The cooling/heating energy generated is sent to the demand-side buildings in the form of airconditioning cold/hot water. After the energy exchange, the water returns to the water collector in the energy station through the return pipe. During the delivery of air conditioning water, the transmission time delay characteristics of the pipelines [47] are considered. The cooling/heating

Case Study and Results
The case study selected was a typical CIES in Tianjin, China. The energy demands in the CIES include cooling, heating and electricity. Renewable energy, such as solar energy and geothermal energy, is utilized according to the actual structure and equipment. The electricity demand in the CIES is satisfied by the external power grid and the PV system together. The cooling demand of the CIES is met by GSHPs and other types of cooling and cooling-storage equipment during the cooling period. The heating load during the heating period is supplied by GSHPs, EBs and hot water tanks. The GSHPs work in intermittent operation mode, including seasonal intermittent operation (continuous operation in winter and summer) and short-term intermittent operation (configuring auxiliary equipment).
Taking the heating period as an example, the heating energy is provided by the GSHPs during the valley price period. During this period, the regenerative electric boiler subsystem stores hot water. After that, during the peak price period, the hot water tanks release the stored energy to meet the heating demand, and the vacancy can be replenished by the EBs or GSHPs.
The cooling/heating energy generated is sent to the demand-side buildings in the form of air-conditioning cold/hot water. After the energy exchange, the water returns to the water collector in the energy station through the return pipe. During the delivery of air conditioning water, the transmission time delay characteristics of the pipelines [47] are considered. The cooling/heating power of the m-th energy source Q m can be calculated as follows: where c and ρ are the specific heat capacity and density of water, respectively, F m denotes the rated flow rate of the m-th primary pump, T R and T S are the temperature of the water in the return pipe and the supply pipe, respectively. According to the actual operation, the total scheduling period is set to 365 days a year, and the scheduling interval is 1 h. The cooling period extends from 21 May to 30 September (133 days in total); the heating period is from 1 January to 16 April and from 24 October to 31 December (175 days in total). The annual and day-ahead forecasting cooling/heating loads of the CIES are shown in Figure 5. The rated peak power of the PV is set to 823 kW, and the maximum allowable power of the tie-line is 10 MW. The specific operating parameters of the various pieces of equipment are shown in Table 1.
The electricity consumption on the demand side will be affected by the operation mechanisms of the electricity market [48], and the electricity price purchased from the external power grid in this paper is based on the time-of-day tariff scheme [49]. The peak periods are 8:00-11:00 and 18:00-23:00, and the electricity price is 1.35 CNY/kW·h. The flat periods are 7:00-8:00 and 11:00-18:00, and electricity price is 0.89 CNY/kW·h. The valley periods are 00:00-7:00 and 23:00-00:00, and the electricity price is 0.47 CNY/kW·h.  Table  1.     According to the historical data of the CIES, the annual heating energy of the GSHP is much larger than the cooling power. If the optimization of each scheduling day during the whole year is carried out without considering the seasonal balance constraints, the annual operation cost should be minimized, but it is too ideal and impractical for the operation of the GSHP. Therefore, five scenarios were compared, as follows, to study the operation of the CIES under various feasible optimization methods: Scenario 1: The annual operation is optimized regardless of the seasonal balance of the GSHP. Scenario 2: To realize the seasonal balance, the annual total energy supply of GSHPs is limited by supply constraints. The constraint value is determined by the annual total cooling/heating energy under annual operation without considering the seasonal balance of the GSHP. Scenario 3: Considering the constraint of seasonal balance, the annual optimization model built in Section 2 is solved to pre-allocate the annual dispatching scheme. The optimization is based on the annual load demand and environmental forecasting data, without taking into account the day-ahead uncertainties.
Scenario 4: Considering the seasonal balance of the GSHP, the day-ahead scheduling is carried out on each scheduling day during the year. The optimization is constrained by the pre-allocated daily cooling/heating energy obtained by the annual optimization in Scenario 3, but the day-ahead uncertainties are not considered.
Scenario 5: Considering the seasonal balance and uncertainty constraints, the day-ahead optimal scheduling is performed on each scheduling day over the year according to the operation strategy proposed in Section 3. Figure 6 shows the final cooling and heating power allocation among the various pieces of equipment in Scenarios 1 to 5. It can be seen from Figure 6a that, regardless of the seasonal balance of the GSHP, the heat energy of the system during the heating period is mainly provided by the GSHPs. However, most of the cooling loads are satisfied by CWCs, and the GSHPs only provide a portion of the cooling energy, which causes a large difference between the annual total cooling and heating power.  In order to balance the cooling and heating load throughout the year, as in Scenario 2, the total heating energy supply of the GSHPs is set equal to the value of cooling energy, according to the initial optimization results in Scenario 1. The optimization results are shown in Figure 6b. It can be seen that the heating power of the GSHPs is significantly reduced in Scenario 2.
Next, Scenario 3 is implemented according to the annual optimization model proposed in Section 2, to obtain a better scheduling scheme under the premise of seasonal balance. After considering the constraints of annual seasonal balance, it can be seen from Figure 6c that the total energy output of the GSHPs during the cooling period increases significantly compared to Scenario 1. It should be noted that the total cooling energy output of the GSHP includes the cooling energy of the GSHP units and the energy released by the cold water tank. This condition is the same for Scenarios 4 and 5, as shown in Figure 6d,e. Figure 7 shows the annual energy supply of the GSHPs in Scenarios 1, 2 and 3. The total heating supply of the GSHPs in Scenario 1 is the largest among all scenarios. The total heating supply in Scenario 2 is the smallest, which is equal to the cooling supply in Scenario 1 and Scenario 2. The cooling energy supply of the GSHPs in Scenario 3 is obviously larger than that in Scenarios 1 and 2.  Figure 8 shows the day-ahead energy supply of the GSHPs in Scenarios 4 and 5. In Scenario 4, the day-ahead cooling/heating energy supply of the GSHPs for each dispatch day is limited according to the pre-allocated value obtained through annual optimization. Thus, the output of the GSHPs under Scenario 4 (shown in Figure 8) is the same as that under Scenario 3 (shown in Figure 7). At the same time, the output in Scenario 5 revolves around the pre-allocated value of Scenario 4 within a certain interval and maintains a dynamic balance after considering uncertainty. Table 2 compares the output of the GSHPs and annual operation costs in different scenarios. Compared with Scenario 1, the operation cost is increased after considering the seasonal balance of the GSHPs in Scenarios 2 to 5, where Scenario 2 has the highest operation cost.   Figure 8 shows the day-ahead energy supply of the GSHPs in Scenarios 4 and 5. In Scenario 4, the day-ahead cooling/heating energy supply of the GSHPs for each dispatch day is limited according to the pre-allocated value obtained through annual optimization. Thus, the output of the GSHPs under Scenario 4 (shown in Figure 8) is the same as that under Scenario 3 (shown in Figure 7). At the same time, the output in Scenario 5 revolves around the pre-allocated value of Scenario 4 within a certain interval and maintains a dynamic balance after considering uncertainty.  Figure 8 shows the day-ahead energy supply of the GSHPs in Scenarios 4 and 5. In Scenario 4, the day-ahead cooling/heating energy supply of the GSHPs for each dispatch day is limited according to the pre-allocated value obtained through annual optimization. Thus, the output of the GSHPs under Scenario 4 (shown in Figure 8) is the same as that under Scenario 3 (shown in Figure 7). At the same time, the output in Scenario 5 revolves around the pre-allocated value of Scenario 4 within a certain interval and maintains a dynamic balance after considering uncertainty. Table 2 compares the output of the GSHPs and annual operation costs in different scenarios. Compared with Scenario 1, the operation cost is increased after considering the seasonal balance of the GSHPs in Scenarios 2 to 5, where Scenario 2 has the highest operation cost.   Table 2 compares the output of the GSHPs and annual operation costs in different scenarios. Compared with Scenario 1, the operation cost is increased after considering the seasonal balance of the GSHPs in Scenarios 2 to 5, where Scenario 2 has the highest operation cost.

Discussion
This paper used five scenarios to compare different optimization strategies. According to the results shown in Section 3, we can clearly see the differences between the different scenarios.
The annual operation under Scenario 1 is optimized regardless of the seasonal balance of the GSHP, and the heating energy is mainly provided by the GSHPs according to the principle of economic optimization, since the COP of the GSHP is higher than that of the EB. Therefore, the operation cost of Scenario 1 is the lowest. However, contrary to the heating period, the GSHPs only provide a small part of energy during the cooling period. The annual total cooling and heating power of the GSHPs is quite different, which means that the heat extracted from the soil during the heating period is larger than the heat released during the cooling period. At the end of the year, there will be a certain drop in soil temperature. Long-term operation with this operation scheme will reduce the soil temperature year by year, so the heating COP of the GSHP will continue to decrease. Therefore, from the perspective of long-term operation, the operation cost of Scenario 1 will keep increasing until the GSHP cannot operate normally.
In Scenario 2, in order to make the total cooling energy supply equal to the heating supply, we set the total energy supply for the whole year. The annual heating energy of the GSHP is much larger than the cooling power in this case. If the cooling supply is forced to be equal to the heat supply, a large amount of heat will be wasted. Thus, the total heating energy supply of the GSHPs is set equal to the value of cooling energy obtained through the initial optimization in Scenario 1. The heating power of the GSHPs is significantly reduced in Scenario 2. However, the operation cost under this scheme is the highest of all the scenarios. Obviously, this is not an economic way to achieve seasonal balance.
Next, Scenario 3 is implemented according to the annual optimization model proposed in this paper. Because of the constraints of seasonal balance, the output of GSHPs in the heating period reduces significantly. This portion of the heating load is transferred to the EBs, thereby reducing the heat exchange of the GSHPs in the soil during the heating period. In this way, the GSHPs can operate in intermittent operation mode and naturally achieve a dynamic balance of soil temperatures. These results prove that a reasonable operation plan can be developed by solving the annual optimization model built in this paper. However, the forecasting data on the year-round time scale may have a large deviation from the day-ahead forecasting. The pre-allocated scheme obtained by annual optimization may not be appropriate for day-ahead scheduling.
After considering the uncertainty of the day-ahead cooling/heating load demand and weather conditions, Scenarios 4 and 5 perform an optimal scheduling based on the day-ahead forecasting data. In Scenario 4, the day-ahead cooling/heating energy supply of the GSHPs for each dispatch day is limited according to the pre-allocated value obtained through annual optimization, while Scenario 5 is implemented according to the day-ahead optimal operation strategy proposed in this paper.
Comparing Scenario 4 and Scenario 5, it can be seen that the output of the GSHPs in day-ahead dispatching deviates from the pre-allocated energy to a certain extent due to the difference between the day-ahead and annual forecasting loads. For the seasonal balance, it is necessary to adjust the control variables to constrain the total cooling/heating energy of the GSHPs. Figure 8 shows that the output in Scenario 5 revolves around the pre-allocated value of Scenario 4 within a certain interval and maintains a dynamic balance. When the day-ahead forecasting load differs greatly from the annual forecasting value, considering the current state of seasonal balance, the result of day-ahead optimal dispatching is allowed to make large or small adjustments according to the pre-allocation value. By adopting this optimization strategy, a more flexible day-ahead scheduling can be achieved.
On the other hand, compared to Scenario 4, the optimized power allocation scheme in Scenario 5 also reduces the operation costs slightly. This is because after considering the uncertainties, the cooling/heating energy supply of GSHPs can be more flexibly adjusted to a certain extent based on the day-ahead forecasting data, which is more accurate, so as to obtain a better optimization strategy.
The heating energy supply of GSHPs is more obviously affected by uncertain factors because there are fewer types of energy supply devices in the heating period than in the cooling period. Compared to Scenario 1, by changing the load power allocation, the annual cooling and heating power of the GSHPs in Scenarios 3, 4 and 5 can be maintained at the same level at the expense of operation costs. However, the operation cost is increased after considering the seasonal balance of the GSHPs in Scenarios 2 to 5, where Scenario 2 has the highest running cost. Due to this high cost, the dispatching scheme of Scenario 2 is not feasible.

Conclusions
To solve the problem of seasonal balance of the GSHP under long-term operation, this paper proposed an annual optimization model and corresponding day-ahead operation strategy for the CIES. A CIES model was built considering the operation constraints and seasonal balance constraints of the GSHP. By optimizing the power allocation among various devices, the most economical operation scheme of the CIES can be obtained, and the seasonal balance of the GSHP can be realized. Then, a scheduling strategy was given which takes into account the uncertainties under daily conditions to obtain a more flexible and reasonable dispatching plan in the day-ahead scheduling. The case studies were based on the structure and parameters of an actual CIES. Then, the annual scheduling optimization was performed separately under the five scenarios. A comparison of different scenarios showed that the optimization model and scheduling strategy proposed in this paper can realize the most economical operation on the premise of satisfying the seasonal balance. Although the yearly operation cost has a slight increase compared to the unbalanced solution, from the perspective of long-term operation, the proposed strategy not only offers lower costs but also meets environmental and sustainable requirements. At the same time, the issue of the unreasonable pre-allocation scheduling scheme due to day-ahead load demands and weather conditions can be solved after considering the uncertain factors. The final results prove that the optimization strategy proposed in this paper can ensure a stable, economic and sustainable operation of the CIES on the premise of the GSHP seasonal balance.
The optimization methods proposed in this paper can achieve seasonal balance throughout the year, but the uncertainties in real-time operation on an hourly time scale have not been taken into account. The soil temperature may keep changing due to the continuous operation of the GSHP, which could cause real-time dynamic changes of the COP of the GSHP in actual operation. At the same time, real-time changes in load demand and environmental factors will also affect the operation of the CIES. Therefore, in order to achieve a more flexible and precise scheduling of the CIES, future studies can consider the effect of soil temperature recovery characteristics and uncertainties in real-time dispatching through methods such like MPC.

Appendix A
The following are models of all subsystems in the CIES and the operation constraints of related equipment.
(1) Regenerative electric boiler subsystem The regenerative electric boiler subsystem consists of EBs and hot water tanks. The heating power of the subsystem Q REB t is the sum of the heating power of EBs Q EB,H t and heating power of hot water tanks Q WT,H t , as shown in Equations (A1) and (A2), which denotes that the EB and hot water tanks cannot work simultaneously.
where U WT,H t is the ON/OFF state of the hot water tank at time t. The heating power of EBs Q EB t,n can be used for space heating and heating storage, as indicated in Equation (A3). As shown in Equations (A4) and (A5), the electric boilers should operate in sequence, and the heating power of each unit is restricted by the upper limit. U EB t,n ≥ U EB t,n+1 , ∀n = 1, 2, · · · , N EB − 1 (A4) 0 ≤ Q EB t,n ≤ U EB t,n Q EB n,max , ∀n ∈ N EB (A5) where Q EB,H t and Q EB,S t , respectively, are the heating and heating storage power of the EB at time t, Q EB n,max is the maximum power of the n-th EB. The heating power of the hot water tank at time t can be calculated by Equation (A6). Equations (A7) and (A8) restrict the heating energy stored in hot water tanks and the heating power. where W WT,H t represents the total heat storage of the hot water tank at time t, ε WT,H is the heat loss rate of the hot water tank, Q WT,H t denotes the heating power of the hot water tank at time t, ∆t represents the scheduling interval, W WT,H max is the upper limit of the energy stored in the hot water tank, Q WT,H max represents the upper limit of the power supply of the hot water tank. Equation (A9) presents the electric power consumption of the EB, and Equation (A10) presents the electric power consumption of the interlocked water pumps (i.e., circulation water pumps and supply water pumps) on both sides of the plate heat exchanger.
Q EB t,n /η EB n + where η EB n is the efficiency of the n-th EB; P EB,WP is the rated power of the auxiliary circulating water pump, N AWP is the number of interlocked water pumps, U AWP t,n represents the ON/OFF state of the n-th (2) Water-cooled chiller subsystem The cooling power of the n-th CWC at time t can be calculated by Equation (A11): where U CWC t,n represents the ON/OFF state of the n-th CWC at time t, N CWC is the total number of CWC, F CWC t denotes the rated flow of chilled water pumps. Equation (A12) restrains the maximum and minimum cooling powers of CWCs, and the operating sequence of the units is limited by Equation (A13).
U CWC t,n Q CWC n,min ≤ Q CWC t,n ≤ U CWC t,n Q CWC n,max , ∀n ∈ N CWC (A12) where Q CWC n,max and Q CWC n,min represent the lower/upper cooling limits of the n-th CWC. The electric power consumption of the water-cooled chiller subsystem can be obtained as follows: t,n /COP CWC + U CWC t,n · P CWC,CWP + P CWC,CP + P CWC,CT where COP CWC is the COP of CWC and P CWC,CWP , P CWC,CP and P CWC,CT are the rated powers of the interlocked auxiliary equipment.
(3) Ice-storage subsystem The ice storage subsystem is composed of DCs and an IT. Equation (A15) indicates that the cooling and ice-making modes of DCs cannot work in parallel. As shown in Equations (A16)-(A18), the operation modes between DCs are limited, and the operating sequence of the units is constrained by Equations (A19) and (A20).
U DC,C t,n + U DC,I t,n ≤ 1, ∀n ∈ N DC (A15) U DC t,n ≥ U DC,C t,n , ∀n ∈ N DC (A16) U DC t,n ≥ U DC,I t,n , ∀n ∈ N DC (A17) U DC t,n ≤ U DC,C t,n + U DC,I t,n , ∀n ∈ N DC (A18) U DC,C t,n ≥ U DC,C t,n+1 , ∀n = 1, 2, · · · , N DC − 1 (A19) U DC,I t,n ≥ U DC,I t,n+1 , ∀n = 1, 2, · · · , N DC − 1 (A20) where U DC t,n is the ON/OFF state of the n-th DC at time t, U DC,C t,n and U DC,I t,n are the cooling/ice-making mode of the n-th DC at time t, N DC is the number of DCs.
During operation, the power of the DCs in the cooling and ice-making modes should meet the lower and upper power constraints: U DC,C t,n Q DC,C n,min ≤ Q DC,C t,n ≤ U DC,C t,n Q DC,C n,max , ∀n ∈ N DC (A21) U DC,I t,n Q DC,I n,min ≤ Q DC,I t,n ≤ U DC,I t,n Q DC,I n,max , ∀n ∈ N DC (A22) where Q DC,C t,n and Q DC,I t,n represent the cooling and ice-making power of the n-th DC at time t and Q DC,C n,max , Q DC,C n,min , Q DC,I n,max and Q DC,I n,min denote the upper/lower power limits of the n-th DC under the cooling and ice-making modes.
The cooling power can be provided by double-duty chillers and an ice-storage tank, which is depicted in Equation (A23). The cooling storage capacity of the IT W IT t in the ice-storage subsystem at time t can be calculated as Equation (A24).
where Q IT t represents the cooling power of the IT at time t, ε IT is the heat loss rate of ice-storage tank, Q DC,I t,n represents the cooling power of the n-th DC at time t.
The electric power consumption of the ice-storage subsystem can be calculated by Equation (A25): where COP DC,C and COP DC,I denote the COP of the DC in cooling mode and ice-making mode, respectively; P EP , P DC,CP , P DC,CT and P IS,CWP are the rated powers of the ethylene glycol pump, chilled water pump, cooling tower and CWP in the ice-storage subsystem, respectively.