Integrated Optimal Control and Borefield Sizing for (Small) Hybrid Heating and Cooling Systems

Proceedings of the 19th IBPSA Conference 0001 Brisbane, Australia, Aug 24-27, 2025 https://doi.org/10.26868/25222708.2025.1442 Integrated Optimal Control and Borefield Sizing for (Small) Hybrid Heating and Cooling Systems Louis Hermans1,2,*, Lone Meertens1,2, Lieve Helsen1,2 1University of Leuven (KU Leuven), Leuven, Belgium 2EnergyVille, Genk, Belgium *Corresponding Author: [email protected] Abstract This study introduces two novel methods for integrated (non-linear) optimal control and borefield sizing for hybrid heating and cooling systems. These methods allow to take into account the optimal control behavior of the operational phase already during the borefield sizing process, avoiding borefield oversizing. Both methods use an iterative approach that combines solving non-linear optimal control problems with borefield sizing algorithms. The key difference between the methods lies in how the borefield is modeled within the optimal control framework: the semi-integrated method uses a borefield fluid temperature profile, while the fully-integrated method incorporates a complete dynamic model. To assess the accuracy and computational efficiency of both methods, they are applied to two use cases: a residential building and an office building. The results reveal that the fully-integrated method generates the smallest borefield sizes, but at the cost of significantly higher computational time. Key Innovations • Integration of non-linear optimal control and borefield sizing in hybrid heating and cooling systems. • Use and comparison of two different dynamic heat transfer models inside boreholes to improve sizing accuracy and avoid oversizing. • Development and comparison of two integrated optimal control and borefield sizing approaches. Practical Implications In hybrid heating and cooling systems with ground-source heat pumps, smart HVAC controllers introduce thermal load uncertainty during the desing phase. This study develops borefield sizing methods that integrate dynamic heat transfer models and optimization, reducing oversizing, improving efficiency, and minimizing total costs (investment + operation). Introduction Hybrid heating and cooling systems, which combine and integrate multiple (renewable) heat and cold sources into one system, have emerged as a promising solution to reduce CO2-emissions in the building sector, as they can achieve a high year-round efficiency by activating the most efficient heat/cold source at every moment. This paper focuses on a specific hybrid heating system that combines an Air-Source Heat Pump (ASHP) with a Ground-Source Heat Pump (GSHP) connected to a borefield. The geothermal part can efficiently supply heat during cold winter days with a relatively high COP and can also provide free cooling in summer (if hightemperature cooling devices are used). Meanwhile, the ASHP is more efficient during intermediate periods when outdoor temperatures are milder. By combining both technologies, the system can reduce overall operational costs and CO₂ emissions, and become future-proof at the same time thanks to a higher intrinsic system flexibility. To keep the investment cost of such a hybrid system acceptable, it is crucial to avoid oversizing its components and in particular the borefield, as a significant portion of the cost arises from drilling the boreholes and installing the associated piping. Accurate borefield sizing is thus essential. Traditionally, borefield sizing is done using software packages such as Earth Energy Designer (EED) or GHEtool, which require an hourly or monthly borefield load profile as an input to subsequently calculate the necessary borefield size (number of boreholes) and depth (Peere & Blanke, 2022). However, during the design stage, this borefield load is uncertain and highly dependent on the eventual control strategy of the system. This is especially the case in hybrid systems, since the control strategy determines which part of the load is covered by the GSHP and which part is covered by the ASHP. Ultimately, a good (optimal) control strategy ensures that the most efficient heat source is activated at all times and takes full advantage of the system's inherent flexibility by anticipating future conditions. Model Predictive Control (MPC) is a promising approach for this, as it uses a dynamic mathematical model of the system, along with forecasts of disturbances such as weather and occupancy behavior, to optimize control inputs over a finite prediction horizon. At each time step, these optimal control inputs are determined by solving a constrained optimal control problem (Drgoňa et al., 2020). Three main groups of MPC formulations exist: (i) purely data-driven methods (black-box), (ii) physics-informed data-driven methods (grey-box), and (iii) completely physics-based methods (white-box) (Drgoňa et al., 2020). Since only white-box MPC formulations do not rely on operational data, only these can be considered during the design phase. Preferably, this white-box MPC formulation also takes into account the most important system nonlinearities (e.g. the temperature-dependent Coefficient-ofPerformance (COP) of heat pumps, cubic relation Proceedings of the 19th IBPSA Conference 0002 Brisbane, Australia, Aug 24-27, 2025 https://doi.org/10.26868/25222708.2025.1442 between mass flow rate and pumping power, …), as the MPC performance is largely dependent on the controller model used. It can typically be enhanced by incorporating additional system dynamics or by refining the existing system dynamics already included (Jansen et al., 2023). In the literature, most studies on borefield sizing, both for purely ground-source and hybrid systems, employ predefined, inflexible building demand profiles, typically generated using building simulation software(Ahmadfard & Bernier, 2019; Blanke et al., 2024; Sharifi et al., 2022). As a result, the interdependence between borefield sizing and (optimal) control is often overlooked. After all, the ground temperature evolution, and thus also the maximum achievable ground-source heat pump’s COP (which is explicitly incorporated in a good optimal controller such as non-linear MPC), depends both on the thermal loads applied to the borefield and on the size of the borefield. The main goal of this paper is to integrate such a nonlinear program-based (NLP), white-box MPC formulation in an optimal borefield sizing strategy such that the overall operational cost of the system and the borefield investment cost are minimized simultaneously. To this end, this paper develops and compares, in terms of accuracy and computational speed, two methods by applying them to two use cases: a residential building and an office building. These cases were selected to represent both heating-dominated (residential) and coolingdominated (office) buildings. High-fidelity, physicsbased, non-linear controller models of these use cases are developed using Modelica, an object-oriented, equationbased, acausal, multi-domain modeling language. Subsequently, TACO (Toolchain for Automated Control and Optimization), an in-house developed Modelicabased toolchain for non-linear white-box MPC, is used to translate these models into NLP-based optimal control problems and efficiently solve them (Jorissen et al., 2019). Both methods follow an iterative procedure where an initial borefield length is set after which the length is updated each iteration until convergence. However, the methods differ in the way the borefield is modeled within the controller model. In the first method (the semiintegrated method), the borefield is represented by a temperature profile, which is an input to the model. In the second method (the fully-integrated method), on the other hand, an optimizable borefield model is added to the controller model such that the ground thermal response is directly taken into account in the control optimization. Use cases and models description In this section, the uses cases and the corresponding controller models are explained in more detail. The hydraulic configuration of both use cases, illustrated schematically in Figure 1, includes: (i) underfloor heating and cooling systems (FH), (ii) a buffer tank, (iii) a modulating air-source heat pump, (iv) a modulating ground-source heat pump, (v) a free cooling heat exchanger (HEX), (vi) a borefield, and (vii) circulation pumps. The corresponding controller model is developed using mostly component and building models of the IDEAS library (Jorissen et al., 2018). However, to satisfy the TACO requirement that all model equations need to be continuous and twice differentiable, the models of the circulation pumps and heat pumps are taken from an inhouse developed component optimization library instead of the IDEAS library. Figure 1: Schematic overview of the hydraulic configuration of the use cases. Building models Linear two-zone models are used to represent the building envelopes. In the office, these zones respectively correspond to the ground and first floor of the building, while in the residential building these zones correspond to a day and a night zone. Each of these zones is equipped with a single underfloor heating model through which the flow is regulated by circulation pumps modeled as pumps with prescribed mass flow rate. Table 1 gives an overview of the most important building parameters. In the office building, the allowable temperature range when people are present is set between 21°C and 24°C for both zones. In the residential building, the temperature limits when people are present are set between 21°C and 26°C in the day zone and between 18°C and 26°C in the night zone. Heating and Cooling System In the heating and cooling system, the following modelling approaches were used: all circulation pumps are modeled as pumps with a prescribed mass flow rate, the free cooling heat exchanger is modeled as a fixed effectiveness heat exchanger with an effectiveness of 0.8 and the buffer tank is modeled as a perfect mixing volume with heat losses to the surrounding air. These heat losses are modeled by a thermal resistance. The 𝑄󰇗𝑑𝑒𝑠𝑖𝑔𝑛 values of the zones are used to calculate the nominal mass flow rates of the underfloor heating and buffer tank, assuming a temperature difference of 7.5°C in the underfloor heating system. The COP of the heat pumps is calculated using Eq. 1. COPHP = COPDef + a (Tcon,out - Tcon,nom) + b (Teva,out – Teva,nom) (1) Where, Tcon,out [°C] and Teva,out [°C] represent the condenser and evaporator outlet temperatures, respectively, while COPDef denotes the default COP corresponding to the nominal outlet temperatures Tcon,nom and Teva. Proceedings of the 19th IBPSA Conference 0003 Brisbane, Australia, Aug 24-27, 2025 https://doi.org/10.26868/25222708.2025.1442 Table 1: Overview of the most important building parameters. Residential Office Day zone Night Zone Total Ground floor 1st Floor Total Floor area 𝑚2 167.0 167.0 334.0 900.0 900.0 1800.0 Air volume 𝑚3 - - 1027.0 - - 5400.0 UA-value building 𝑊𝐾 ⁄ - - 391.7 - - 496.4 𝑸󰇗design 𝑘𝑊 10.3 10.4 20.7 17.7 24.7 42.4 n50 𝑙/ℎ 8.0 8.0 - 3.0 3.0 - Window area 𝑚2 20.1 20.1 40.2 20.1 20.1 40.2 Glazing type - Double Double - Triple Triple - For the air-source heat pump, the default COP is determined at a nominal condenser temperature of 35°C and a nominal evaporator temperature of 12.5°C, resulting in a COPDef of 5.9. Meanwhile, for the ground-source heat pump, the default COP is based on a nominal condenser temperature of 45°C and a nominal evaporator temperature of 4°C, yielding a COPDef of 4. The part-load ratio dependency of the COP is not explicitly taken into account. The nominal mass flow rates of the heat pumps are determined based on their capacities, assuming a temperature difference of 5 °C in the condenser and 3 °C in the evaporator. Finally, Table 2 gives an overview of all relevant borefield parameters. Occupancy and weather conditions The weather conditions in the model are assumed to be the same each year and taken as a Typical Meteorological Year (TMY) in the city of Uccle in Belgium. The corresponding TMY-file can be found in the IDEAS library. In the office building, each zone is occupied by 60 people during office hours. The metabolic heat gains from occupants are modeled as 73 W of sensible heat and 45 W of latent heat per person (ASHRAE, 2009). Office hours are considered during all weekdays from 8 am until 6 pm. In the residential building, no occupants are present during office hours, with occupancy only occurring outside these hours. Moreover, no metabolic heat gains are considered in the residential building. Methodology This section provides a detailed explanation of the two integrated borefield sizing methods developed in this study: the semi-integrated borefield sizing method and the fully-integrated borefield sizing method. Both methods use an iterative approach that couples a non-linear optimal control problem with a borefield sizing strategy. They differ in the way the borefield is modeled within the optimal control problem. First, the two borefield sizing strategies and the considered optimal control problem are explained in detail, followed by a schematic overview and a discussion of the integrated methods. Borefield sizing Accurate borefield sizing requires reliable predictions of both ground and fluid temperatures during operation. Heat transfer processes in borefields occur over multiple temporal and spatial scales. Table 2: Overview of the borefield parameters. Description Unit Value Configuration parameters Borehole radius m 0.066 Borehole buried depth m 0.8 Borehole spacing m 6 Ground parameters Ground conductivity W/(mK) 1.9 Ground specific heat capacity J/(kgK) 1400 Ground density kg/m3 900 Undist. ground temperature °C 10 Ground temperature gradient °C/m 1 Grout parameters Grout conductivity W/(mK) 2 Grout specific heat capacity J/(kgK) 840 Grout density kg/m3 1818 Pipe parameters Heat exchanger type Single U-tube Inner pipe radius m 0.0131 Outer pipe radius m 0.016 Pipe conductivity W/(mK) 0.38 Pipe spacing m 0.043 Temperatures Max. allowed temperature °C 16 Min allowed temperature °C 3 Residential Borefield configuration 2 x 2 Borefield mass flow rate kg/s 0.8 Office Borefield configuration 4 x 3 Borefield mass flow rate kg/s 2.4 On short time scales (i.e. from minutes to hours), internal effects within the borehole, such as those associated with the fluid, pipes, and grout, dominate the heat transfer dynamics. In contrast, over longer periods, thermal conduction in the surrounding ground becomes dominant. Proper borefield sizing ensures that heat transfer does not result in fluid temperatures exceeding specified limits during the operational lifetime of the borefield, which is assumed to be 20 years in this study. Proceedings of the 19th IBPSA Conference 0004 Brisbane, Australia, Aug 24-27, 2025 https://doi.org/10.26868/25222708.2025.1442 The sizing strategies used in this study, detailed further below, build upon the Finite Line Source (FLS) model of Cimmino (Cimmino, 2019) and implemented in the opensource Python package pygfunction (Cimmino & Cook, 2022). In this model, the borehole is represented as a finite line, and the concept of g-functions, as introduced by Eskilson (Eskilson, 1987), is used to characterize heat transfer external to the borehole. The g-functions serve as step-response functions relating the heat extraction rate in the borefield to the effective temperature variation at the borehole walls. Accordingly, the borehole wall temperature, 𝑇𝑏, at time 𝑡𝑛 is determined by superimposing discrete loads 𝑄󰇗𝑖, as shown by Eq. 2 (Cimmino, 2019). 𝑇𝑏(𝑡𝑛)= 𝑇𝑔− 1 2𝜋𝑘𝑠𝐿∑(𝑄󰇗𝑖−𝑄󰇗𝑖−1)𝑔(𝑡𝑛 − 𝑡𝑖−1 𝑡𝑠) 𝑛 𝑖=1 (2) With 𝑇𝑔 the undisturbed ground temperature, L the total borefield length, 𝑘𝑠 the ground thermal conductivity, and 𝑄󰇗𝑖 and 𝑡𝑖 the thermal load and time at instance i. 𝑡𝑠(= 𝐻2/9𝛼𝑠), with H denoting the borefield depth and 𝛼𝑠 the ground thermal diffusivity, is the borehole characteristic time. By introducing the borehole equivalent thermal resistance, 𝑅𝑏∗, which accounts for all steady-state thermal interactions within the borehole, Eq. 3 converts the wall temperature into the average fluid temperature, 𝑇𝑓(𝑡). 𝑇𝑓(𝑡)= 𝑇𝑏(𝑡)−𝑄󰇗𝑝,𝑡𝑅𝑏∗ (3) Where 𝑄󰇗𝑝,𝑡 denotes the thermal peak load at time t. Although this conventional formulation is widely used, its steady-state nature implies an immediate transfer of heat from the fluid to the ground. This assumption neglects the transient behaviour of the fluid, pipe, and grout, leading to an overestimation of the heat exchange and, consequently, an oversized borefield. In light of these limitations, the current study uses and compares two different approaches that account for shortterm dynamics to size the borefield. The first approach, referred to as the GHEtool sizing approach, uses the borefield model described by Meertens et al. (Meertens et al., 2024) and implemented in the open-source GHEtool (Peere and Blanke, 2022), which incorporates an explicit one-dimensional numerical model. This dynamic model accounts for the thermal inertia of the materials between the fluid and the ground and accurately captures the transient heat transfer processes. The second approach, referred to as the Modelica sizing approach, uses a slightly adapted version of the IDEAS borefield model. The primary reason for employing both approaches is to accurately capture short-term dynamic effects within the boreholes. However, each sizing approach incorporates these dynamics differently: the GHEtool sizing approach utilizes the one-dimensional numerical model of Xu and Spitler, while the Modelica sizing approach employs the IDEAS model, which simulates borehole dynamics using axial discretization and a resistance-capacitance network to represent internal thermal resistances between individual pipes and between each pipe and the borehole wall. By incorporating this dynamic behaviour, it becomes necessary to account for the cylindrical geometry of the borehole when calculating the heat transfer outside the borehole on short time scales. To achieve this, the cylindrical correction proposed by Li et al. (2013) is applied to the FLS-generated g-functions according to Eq. 4. 𝑔(𝑡)= 𝑔𝐹𝐿𝑆(𝑡)+(𝑔𝐶𝐻𝑆(𝑡) − 𝑔𝐼𝐿𝑆(𝑡)) (4) Where the additional term represents the difference between the cylindrical heat source (CHS) generated gfunction and infinite line source (ILS) generated gfunction. This correction makes the g-function method applicable for time values below the borehole characteristic time. For longer time scales, the correction converges to Eskilson’s correction factor (Eskilson, 1987), and the adjusted g-function approaches the original FLS formulation. The dynamic average fluid temperature, 𝑇𝑓, computed from the numerical model is subsequently employed in an iterative sizing procedure. The required borehole length is determined iteratively using hourly thermal pulses and the formulation provided by Peere et al. (2023) is given by Eq. 5 𝐻=𝑚𝑎𝑥{𝑚𝑎𝑥(𝑇𝑓)−𝑇𝑔 𝑇𝑚𝑎𝑥−𝑇𝑔 𝑚𝑖𝑛(𝑇𝑓)−𝑇𝑔 𝑇𝑚𝑖𝑛−𝑇𝑔}∙𝐻𝑝𝑟𝑒𝑣 (5) In this formulation, 𝑇𝑚𝑎𝑥 and 𝑇𝑚𝑖𝑛 denote the maximum and minimum allowable average fluid temperature in the borefield, respectively, and 𝐿𝑝𝑟𝑒𝑣 is the borehole length from the previous iteration. This iterative process ensures that the borefield is sized accurately by fully accounting for the transient heat transfer behaviour of the system. Optimal Control Problem The considered optimal control problem is formulated by Eq. 6. 𝑚𝑖𝑛 𝒐(𝑡)∫ (𝐽𝑒𝑙(𝑡)+∑𝑤 2 𝑛=1 (ℎ𝑛2+𝑐𝑛2)) d 𝑡𝑖+𝛥𝑡𝑝𝑟 𝑡𝑖𝑡 (6𝑎) s.t. 𝒅𝒙 ( 𝑡 ) 𝑑𝑡 =𝑭 ( 𝒙 ( 𝑡 ) ,𝒛 ( 𝑡 ) ,𝒐 ( 𝑡 ) ,𝑡 ) (6𝑏) 0 = 𝑯(𝒙(𝑡),𝒛(𝑡),𝒐(𝑡),𝑡) (6𝑐) 𝒙(𝒕𝟎) = 𝒙𝟎 (6𝑑) 𝑇𝑚𝑖𝑛,𝑛(𝑡)−𝑇𝑧,𝑛(𝑡)≤ℎ𝑛 (6𝑒) 𝑇𝑧,𝑛(𝑡)−𝑇𝑚𝑎𝑥,𝑛(𝑡)≤𝑐𝑛 (6𝑓) ℎ𝑛,𝑐𝑛≥0 𝑛=1,2 (6𝑔) Here, x(t) represents the vector of state variables, while o(t) represents the vector of optimization variables, and z(t) corresponds to the vector of remaining algebraic variables. The temperature in the nth-zone is given by Tz,n , with Tmin,n and Tmax,n representing the minimum and maximum permissible temperatures in that zone, respectively. The governing system equations are described by the functions F and H . The considered optimal control variables o(t) are the modulation degree of the heat pumps, and the prescribed mass flow rates of Proceedings of the 19th IBPSA Conference 0005 Brisbane, Australia, Aug 24-27, 2025 https://doi.org/10.26868/25222708.2025.1442 the circulation pumps connected to the free cooling heat exchanger, and underfloor heating models. TACO directly infers all these equations and variables from the Modelica controller model. Subsequently, these equations are manipulated using computer algebra, discretized using a non-equidistant time grid, and solved by a derivative-based NLP solver (Jorissen et al., 2019). The primary objective of the optimal control problem is the minimization of the heat pumps’ and free cooling circulation pump’s electrical energy use, Jel(t) , while also ensuring thermal comfort. This last requirement is enforced through weighted soft constraints introduced by the slack variables hn and cn . To balance thermal comfort and computational efficiency (reasonable convergence time) the weighting factor w of these soft constraints is set to 10000 W/K2. The optimal control problem has a timestep of 1 hour and a prediction horizon of 1 year or 20 years. Semi-Integrated Method Within the semi-integrated method, the borefield is not explicitly modeled in the optimal control problem, but represented by a precalculated fluid temperature profile. The method consists of 3 main blocks: an initialization block, an optimal control block, and a sizing and temperature profile generation block. Figure 2 shows a schematic overview of the method. To initialize the method, a one-year optimal control problem, which assumes a constant ground temperature of 10°C, is solved by TACO and the resulting borefield load profile is used by the GHEtool sizing approach to calculate an initial borefield depth and fluid temperature profile. After this initialization, the initial borefield fluid temperature profile is sent to the method’s optimal control block, where four approaches for calculating an optimized 20-year borefield load profile are considered and compared. The first approach, termed ‘Single’, uses the complete borefield fluid temperature profile to perform a single 20-year optimal control optimization, with the resulting load profile directly sent to the sizing block. The second approach, ‘Multi’, splits the 20-year temperature profile into 20 individual yearly profiles, after which the optimal control problem of each year is solved. The resulting yearly load profiles are then recombined to form a single 20-year borefield load profile. The third approach, ‘Dual’, optimizes only the first and last year of the borefield operation 1 , with linear interpolation applied to the intermediate years to construct a 20-year borefield load profile. Finally, the ‘Simple’ approach calculates an average one-year borefield fluid temperature profile from the 20-year temperature profile and optimizes only this year, repeating the resulting load profile 20 times to create the full 20-year borefield load profile. It is important to stress that full-year or 20-year optimal control problems, with full knowledge of the boundary conditions within that period, are solved in the 1 It is known that the temperature limitation is hit either in the first or last year of operation (Peere et al., 2021). optimal control block. Hence, the resulting optimal controls represent an idealized situation that serves as a proxy for the controls of a real non-linear model predictive controller that typically has a prediction horizon of only a few days. Figure 2: Schematic overview of the semi-integrated sizing method. After generating the 20-year borefield load profile in the optimal control block, the semi-integrated method uses that load profile to size the borefield and generate a borefield fluid temperature profile using either the GHEtool sizing approach or the Modelica sizing approach. The resulting borefield depth H is compared to the previous depth, and if the difference is less than 0.1 m, the convergence criterion is met, and the calculated depth is the final one. However, if the criterion is not met, the temperature profile is sent back to the optimal control block to initiate a new iteration. Fully-Integrated Method The fully-integrated borefield sizing method, schematically shown in Figure 3, follows a similar iterative procedure as the semi-integrated method. However, an initialization block is not needed, as only an initial borefield depth is required, which is set to 100 m. The key difference from the semi-integrated method lies in the optimal control problem: in the fully integrated method, an explicit borefield model is used, allowing the optimal controller to directly account for the borefield’s thermal response within the optimization. Since TACO requires all model equations to be continuous, the default Proceedings of the 19th IBPSA Conference 0006 Brisbane, Australia, Aug 24-27, 2025 https://doi.org/10.26868/25222708.2025.1442 IDEAS borefield model cannot be used, as it relies on an event-based load-shifting algorithm for load aggregation. To overcome this limitation, this work follows the approach of Cupeiro Figueroa et al. (2020), which modifies the IDEAS model by applying the QUICK (Quadratic Upstream Interpolation of Convective Kinetics) method (Cupeiro Figueroa et al., 2020). This adaptation converts the load-shifting algorithm from an event-based approach to a continuous one, thereby ensuring compatibility with TACO. Within the optimal control block, an optimal control problem of 20 years is solved, and the resulting borefield load profile is subsequently sent to the sizing block. Again, the full 20-year optimal control problem is an idealization of the actual model predictive control behavior. Figure 3: Schematic overview of the fully-integrated sizing method. The subsequent steps of the iterative procedure in the fully integrated method follow those of the semi-integrated method, with one extra difference: if the convergence criterion is not met, the calculated borefield depth is sent to the optimal control block instead of a calculated temperature profile. This allows to update the borefield depth within the controller model and start a new iteration. Note that in the fully-integrated method only the GHEtool sizing approach is used, as the results will show that the difference with Modelica sizing approach is small. Results Sizing Accuracy An overview of the resulting borefield sizes for both methods and use cases can be found in Figure 4 and Figure 5, where the percentages next to each bar represent the relative oversizing of the semi-integrated methods compared to the fully-integrated method. The results are presented for two scenarios: (i) a reference scenario, where only a GSHP is used to verify the different sizing approaches, and (ii) a hybrid scenario, where a GSHP is combined with an ASHP. In both scenarios, the total installed heat capacity matches the design heat loss 𝑄󰇗𝑑𝑒𝑠𝑖𝑔𝑛of the respective building. (with each heat pump assumed to contribute to 50% of the capacity in the hybrid scenario). Figure 4: Resulting borefield depths of all sizing approaches for the residential building (incl. the relative oversizing of the semi-integrated methods). Figure 5: Resulting borefield depths of all sizing approaches for the office building (incl. the relative oversizing of the semi-integrated methods). A first interesting finding is that the GHEtool sizing approach mostly results in a slightly smaller borefield depth compared to the Modelica sizing approach, with an average absolute difference of 1.9 m and an average relative difference of 2%. Since this difference is minor, Proceedings of the 19th IBPSA Conference 0007 Brisbane, Australia, Aug 24-27, 2025 https://doi.org/10.26868/25222708.2025.1442 both approaches can be considered reliable, with the variation likely stemming from the slight differences in the modelling approach. Furthermore, Figure 4 clearly shows that when the residential building is equipped with a GSHP only, all considered sizing approaches yield nearly the same borefield size, with a maximum relative oversizing of 2.4%. However, this trend does not hold for the office building, where the relative oversizing in the ‘only GSHP’ scenario ranges from +6.8% (Single) to +17.1% (Simple). The primary reason for this discrepancy between the fully-integrated method and semi-integrated method is the reduced cooling provided in the fully-integrated method for the office. This effect is visible in Figure 6, which illustrates the temperature evolution of one of the office zones in both methods. The significantly higher zone temperature in the fully-integrated case highlights the limited cooling provision. Since the office is a coolingdominated building, this reduced cooling has a notable impact on borefield sizing, as cooling demand, particularly peak loads, determines the required borefield depth. Figure 6: Example of temperature evolution in one of the office zones. An oversizing trend is also observed in the ‘hybrid’ scenario of both use cases, where the semi-integrated method consistently results in an oversized borefield by at least 9.7%. This suggests that to minimize oversizing, the fully-integrated method is the preferred approach. The reason for this difference is that the optimizer in the fully integrated method directly accounts for the ground’s thermal response. As a result, it increases the ASHP to GSHP usage ratio in the residential building while decreasing it in the office building compared to the semiintegrated methods. In the residential building, this adjustment slows the decline of borefield temperature (because less heat is extracted from the ground), improving the GSHP’s COP. Conversely, in the office building, the optimizer prioritizes GSHP operation, reducing long-term borefield temperature rise (because less heat is injected) and enhancing cooling potential. Finally, when comparing the resulting borefield depths of the semi-integrated methods, it becomes clear that the ‘Single’ approach is the most accurate as it has the lowest oversizing (4.7% on average in the ‘onlyGSHP’ scenario and 14.2% on average in the ‘hybrid’ scenario). Computational Time Figure 7 provides an overview of the total computational time in the ‘hybrid’ scenario for each use case. It is immediately clear that the fully integrated method requires significantly more computational time (3 orders of magnitude) than the semi-integrated methods. Figure 7: Overview of the total computational times per iteration for all sizing methods. This substantial difference is due to the fully integrated method requiring the optimal control problem to be recompiled at every iteration. This is necessary because TACO does not allow parameter changes after compilation while each iteration the borefield depth parameter has to be updated. In contrast, the semiintegrated methods only require compilation once during initialization, because the borefield temperature profile is treated as an input to the controller model, which can be changed without recompilation. Moreover, when comparing the computational times of the semi-integrated methods, it becomes clear that computation time increases with the number of optimal control problems solved per iteration. More specifically, the ‘Simple’ approach is the fastest, followed by the ‘Dual’ approach, with the ‘Multi’ approach being the slowest. Additionally, splitting a 20-year optimal control problem into 20 one-year problems significantly reduces computation time, as seen in the lower computational time of the ‘Multi’ approach compared to the ‘Single’ approach (one order of magnitude difference). Finally, the average initialization times (across both use cases and scenarios) for the ‘Simple,’ ‘Dual,’ ‘Multi,’ and ‘Single’ semi-integrated methods are 1514, 1495, 2136, Proceedings of the 19th IBPSA Conference 0008 Brisbane, Australia, Aug 24-27, 2025 https://doi.org/10.26868/25222708.2025.1442 and 34321 seconds, respectively. The initialization time for the ‘Single’ approach is thus an order of magnitude longer than the others, which also mainly explains the difference in overall computational time. This difference is primarily caused by the extensive time required to compile and solve the initial the 20-year optimal control problem. This again highlights the advantage of splitting the single 20-year optimal control problem in 20 1-year optimal control problems. Conclusion This paper introduces two novel methods for integrating non-linear optimal control in borefield sizing in hybrid heating and cooling systems: a fully-integrated method and a semi-integrated method. A comparative analysis demonstrates that the fully-integrated method offers the highest accuracy in minimizing borefield size while ensuring optimal system performance. However, the computational demands of this approach are significantly higher compared to semi-integrated methods. For practical applications, the semi-integrated methods, which yield more conservative sizes, may offer a wellbalanced trade-off between computational efficiency and borefield size accuracy, despite being slightly less precise. Among the four semi-integrated approaches considered, Single, Multi, Dual, and Simple, the Single approach demonstrated the highest accuracy in borefield sizing, but had a significantly higher computational time. In contrast, the Simple approach was the most computationally efficient but resulted in the largest oversizing, making it less suitable. The Multi and Dual approaches provided intermediate results, offering a better balance between accuracy and computational efficiency. The methods have been demonstrated for individual buildings, though they can be easily generalised to collective systems for multiple buildings. For the latter, the borefield investment cost will probably dominate, justifying the use of the fully-integrated method that avoids oversizing. Nevertheless, the semi-integrated methods can be very useful to get first estimates. Acknowledgement The authors are grateful to Filip Jorissen for all TACO developments. The authors acknowledge the funding by the European Union through the SEEDs project under the Horizon Europe Programme. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or CINEA. Neither the European Union nor the granting authority can be held responsible for them. References Ahmadfard, M., & Bernier, M. (2019). A review of vertical ground heat exchanger sizing tools including an intermodel comparison. In Renewable and Sustainable Energy Reviews (Vol. 110, pp. 247–265). Elsevier Ltd. https://doi.org/10.1016/j.rser.2019.04.045 ASHRAE. (2009). ASHRAE Handbook - Fundamentals (Vol. 18). Blanke, T., Born, H., Döring, B., Göttsche, J., Herrmann, U., Frisch, J., & van Treeck, C. (2024). Model for dimensioning borehole heat exchanger applied to mixedinteger-linear-problem (MILP) energy system optimization. Geothermal Energy, 12(1). https://doi.org/10.1186/s40517-024-00301-w Cimmino, M. (2019). Semi-Analytical Method for g-Function Calculation of bore fields with seriesand parallelconnected boreholes. Science and Technology for the Built Environment, 25(8), 1007–1022. https://doi.org/10.1080/23744731.2019.1622937 Cimmino, M., & Cook, J. C. (2022). pygfunction 2.2: New features and improvements in accuracy and computational efficiency. IGSHPA2022, 45–52. Cupeiro Figueroa, I., Cimmino, M., & Helsen, L. (2020). A methodology for long-term model predictive control of hybrid geothermal systems: The shadow-cost formulation. Energies, 13(23). https://doi.org/10.3390/en13236203 Drgoňa, J., Arroyo, J., Cupeiro Figueroa, I., Blum, D., Arendt, K., Kim, D., Ollé, E. P., Oravec, J., Wetter, M., Vrabie, D. L., & Helsen, L. (2020). All you need to know about model predictive control for buildings. In Annual Reviews in Control (Vol. 50, pp. 190–232). Elsevier Ltd. https://doi.org/10.1016/j.arcontrol.2020.09.001 Eskilson, P. (1987). Thermal analysis of heat extraction boreholes. University of Lund, Sweden. Jansen, J., Jorissen, F., & Helsen, L. (2023). Optimal control of a fourth generation district heating network using an integrated non-linear model predictive controller. Applied Thermal Engineering, 223. https://doi.org/10.1016/j.applthermaleng.2023.120030 Jorissen, F., Boydens, W., & Helsen, L. (2019). TACO, an automated toolchain for model predictive control of building systems: implementation and verification. Journal of Building Performance Simulation, 12(2), 180– 192. https://doi.org/10.1080/19401493.2018.1498537 Jorissen, F., Reynders, G., Baetens, R., Picard, D., Saelens, D., & Helsen, L. (2018). Implementation and Verification of the IDEAS Building Energy Simulation Library. Journal of Building Performance Simulation, 11(6), 669–688. https://doi.org/10.1080/19401493.2018.1428361 Meertens, L., Peere, W., & Helsen, L. (2024, May 28). Influence of short-term dynamic effects on geothermal borefield size. IGSHPA Research Conference 2024. Peere, W., & Blanke, T. (2022). GHEtool: An open-source tool for borefield sizing in Python. Journal of Open Source Software, 7(76), 4406. https://doi.org/10.21105/joss.04406 Peere, W., Picard, D., Cupeiro Figueroa, I., Boydens, W., & Helsen, L. (2021, September 1). Validated combined first and last year borefield sizing methodology. https://doi.org/10.26868/25222708.2021.30180 Sharifi, M., Figueroa, I. C., Mahmoud, R., Himpe, E., Helsen, L., & Laverge, J. (2022). Early-stage optimal design of hybrid GEOTABS buildings in terms of costs and CO2 emissions. Energy Conversion and Management, 257. https://doi.org/10.1016/j.enconman.2022.115392