Optimized design and integration of energy storage in Solar-Assisted Ground-Source Heat Pump systems

The integrated use of multiple renewable energy sources to increase the efficiency of heat pump systems, such as in Solar Assisted Geothermal Heat Pumps (SAGHP), may lead to significant benefits in terms of increased efficiency and overall system performance especially in extreme climate contexts, but requires careful integrated optimization of the different system components. In particular, thermal storages take a fundamental role in optimizing the integration of renewable energy sources and the system operation. This work investigates the potential design optimization of a SAGHP system in a mountain site by exploring many different alternatives to optimize the mutual relationship between the solar field, the geothermal field and the water thermal storages. This is done through an original simulation-based multi-objective optimization framework considering energy efficiency and economic feasibility, which allows appraising the impact of the different design alternatives on the overall system performance and on the dynamics of the different system components. Results identify a set of optimized system configurations that optimize the integrated exploitation of the different thermal sources showing a potential increase of the overall system performance leading to 34% lower global cost compared to the initial design. High robustness of the optimal design solutions is reported with respect to the current context of high economic uncertainty.


Introduction
According to latest IEA projections (International Energy Agency 2020), heat pump is a key clean energy technology that will drive decarbonization of the building sector in the next few year s, especially relying on the integration o f multiple renewable energy sources for in creasing overall efficiency and applicability in different building typologies and boundary conditions (Alberti et al. 2018).
In general, it has been demonstrated that the combination of the solar and geothermal energy sources in Solar Assisted Ground-source Heat Pump s (SAGHP) (Ozgener a nd Hepbasli 2007) may po tentially lead to multiple benefits in terms of increased global COP of the system (Wo łoszyn and Gołaś 2017), and geothermal source regeneration for lower energy consumptions and higher ren ewable energy exploitation (Girard et al. 2015).
Thermal energy storages at both source and use sides is important for optimal integration of renewable energy sources considering their flexibil ity (Maturo et al. 2022), climate boundaries and eco nomic feasibility (Osterman and Stritih 20 21).In par ticular, economic feasibility is a key-point constituting a barrier for this technology for the large investment costs (Liu et al. 2021).Therefore, the optimization of the solar i ntegration and of thermal storages is also crucial to reduce the size of the boreholes (Nouri et al. 2019), which are responsible for the higher costs related to this type of systems and may lead to high payback times unless tailored incentives and taxation policies are implemented (Rivoire et al. 2018).
Since the performance of such systems is complex and time variable in the long term (the thermal depletion of the ground can occur in a period of several years), it may not be sufficient to analyze short-term data available from monitoring but there is the need to use detailed transient simulation models for the design and optimization of th e List of symbols C a,i (j) annual cost for component j at the year i C G (τ, P)  various system layouts, uses and operation parameters of the storages (Schellenberg et al. 2018; B uonomano and Guarino 2020; Koşan and Aktaş 2021), which can be optimized according to different objectives (Shah et al. 2020).The TRNSYS ® software has been recognized as very useful for the layout organization an d simulation accuracy of th e different system components (ground heat exchanger, heat pump, thermal stor ages, controllers, building loads, ...), as demonstrated by several studies (Desideri et al. 2011;Buonomano et al. 2015;Calise et al. 2016 ;Marini et al. 2019).Further, the coupling between transient simulation and artificial intelligence-based meth ods, such a s metaheuristics optimization algorithms, has emerged for its ability to drive multiple simulation runs towards the search of optimal design configurations while analysing their performance in details in a complex set of constraints (Schellenberg et al. 2020).However, even though recognized design guidelines exist and many different advanced models for detailed performance assessment of the different s ystem components have been reported (Hein et al. 2016;Luo et al. 2019;Zhang et al. 2022), the o peration and the long-term performance of SAGHP systems are very sensitive to the local condition s (climate, soil conditions, hydraulic parameters, energy storage ratio, etc.) (Dalla San ta et al. 2020) and should be considered for the specific set of constraints characterizing each case (D' Agostino et al. 2020).That is why po pulating the literature with detailed case-by-case analysis for the ir application in multiple different contexts is required to overcome the current barriers to the wide adoption of the technology in different climate scenarios (Ferrara et al. 2015), supported by a standardized methodological framework (Ferrara et al. 2019) for modelling and optimizing th e different components and the complex set of constraints of such systems.

Aim of the work
Starting from the initial design of a system providing thermal energy for heating and DHW in the context of an Alpine ski park, this work aims at investigating the potential design optimization of a SAGHP system in a mountain site from both technical and economic points of view.To do so, the following innovation objectives will be met: Setup of a dedicated transient simulation-based optimization framework and definition of constraints occurring within the different variables related to thermal storages; Definition of key-performance indicators that can be used to set up the system multi-objective optimization to determine technical and financial feasibility; Dynamic exploration and multi-objective optimization of different design alternativ es that are abl e to optimize the mutual r elationship between the solar field, the geothermal field and the thermal storages; Quantification of potential r eduction/increase of energy and cost performance within the defined design space; Sensitivity analysis to identify the stability of the resulting optimal design with the current uncertain financial parameters.

The case study
The case stud y is a SAGHP system for a newly buil t restaurant in a mountain site.It is located at 2500 m above the sea level in the high Alps region in rigid winter climate conditions (the annual degree days are 4524 °Cd, the mean outdoor air temperature is −3.8 °C and the winter design temperature is −20 °C).The energy demand of the building (it is expected to operate from 1 st November to 30 th April) was estimated through a dynamic simulation model based on real historical weather data collected in the location and typical daily occupancy schedules of similar buildings at the site.The total seasonal thermal energy demand resulting from simulation of the building model is 119027 kWh (257 kWh th /m 2 ), of which 38.5% are n eeded for space heating, 50.2% for fresh air heating and 11.3% for domestic hot water production.
The system is designed to meet the building energy needs by providing water at three different thermal levels, based on user requests: (1) high temperature technical water (60 °C) to be used for air heaters and air handling units (AHU); (2) low temperature technical water for s pace heating radiant panels and AHU pre-heating (40 °C ); (3) domestic hot water (52 °C).
As shown in Figure 1 the source-side loop is composed of a geothermal field and a hot thermal storage where a grey DHW recovery system and a solar thermal collector field are integrated with the purpose of increasing t he temperature of the heat transfer fluid (water-ethylene glycol, 40%) at the evaporator side to increase the heat pump efficiency.In the initial design configuration (INI), the geothermal field is composed on 6 vertical boreholes of 200 m connected in parallel, the solar field is composed of 4 thermal collectors (2 m 2 each), the volume of the integration thermal storage is 500 l and the volume of the grey water storage is 4150 l.
The demand-side loop is composed of three water thermal storages, one for each requir ed thermal level at the demand side.In the initial design configuration (INI), 2000 l, 1000 l and 5 00 l are the volumes of the high-temperature thermal storage, the low-temperature thermal storage the DHW storage, respectively.More details on the system initial design configuration can be found in Fabrizio et al. (2015).

Settings of simulation model
The system dynamic simulation model was created in TRNSYS as reported in Figure 2.
To ensure model accuracy, the different model components were selected based on recognized model accuracy validation studies.The largely val idated Type 1b was used for modelling flat plate solar th ermal collectors (Banister et al . 2014).T ype 668 was used to model the water/water heat pump through interpolation of user-defined data, based on performance curves provided by manufacturers (Dott et al. 2013).Type 557 was used to model the geothermal field (Hellstrom 1989;Ruiz-Calvo et al. 2017;Li et al. 2023) with the boreh oles of diameter 40 mm and water flow rate of 0.483 l/s.The soil thermal conductiv ity and the soil th ermal capacity were set to 1.9 W/(m•K) and 2.0 MJ/(m 3 •K), respectively, and the vertical thermal gradient  et al. 2005), was used to mo del the five stratified thermal storages included in the s ystem.More details about the modelling and validation of other system components and control logics can be found in Urone (2015).
Once ensured the appropriate Types for th e different system components were selected, the connections between different types were made based on actual system d esign and according to modeling approaches used in previo us literature reporting the TRNSYS simulation and validation of similar solar assisted grou nd-source heat pump systems (Wang et al. 2 012;Nam et al . 2015;Yu et al. 2016;Li et al. 2018;Sakellariou et al. 2019;Hou et al. 2020).Based on the described TRNSYS code, the model template for the optimization was created for their use in the Genopt ® software.Such template contains information on code locations where the optimization program should input the set of optimization variables, so that a sea sonal TRNSYS simulation with a different system design configuration can be run at each iteration.A tailored "control center" equation type was cr eated to properly link the user-defined set of independent optimization variables to the set of dependen t variables that need to be modified and in putted to th e relevant model types at each iteration (e.g. using the size of the storage volume as indepe ndent user-defined variables means modifying the ent ire set of geo metrical parameters of the storage Type 60 -orig inal equations and constraints were defined for the purpos e, according t o Table 1 and  Table 2).
The optimization process was set to be driven by the binary version of pa rticle swarm optimization (PSO) algorithm, for which th e set of algorith m parameters (particles number = 5; cognitive acceleration c 1 = 3; social acceleration c 1 = 1; max swarm velocity = 4) were demonstrated to be optimal for building system optimization problems (fast convergence with high neighborhood exploration -refer to Ferrara et al. (2017) for details).The maximum number of iteration was set to 2000 for each optimization process.

Optimization variables and design space
The set of variables defining the design space for this multienergy system optimization problem is reported in Table 1.
The range an d the ste p for their variation is define d according to feasibility criteria (market availability, available space on th e site for installation), as follows.In fact, th e maximum value of the var iable NColl, indicating the number of solar collectors, is set to the maximum area of the building roof that can be covered by solar collectors (considering a 2 m 2 area fo r each panel, the variation of NColl ca n increase the solar collector area up to 80 m 2 ).Similarly, the maximum number of geothermal boreholes is related to the available space of the geothermal site.The characteristics of thermal storages (volume and serpentine exchange surfaces) were set acco rding to market availability.It has to be remarked that the volume of the integration hot storage is dependent on two optimization variables, one is the ratio between storage volume and solar collecto r area (variabl e STVint) and the other is the number of solar collectors (NColl), from which the solar collector area can be derived considering a 2 m 2 area for each collector.Depending on the values of t he two independent variables, the volume of the integration hot storage may range between 500 l (set as minimum volume constraint) and 6000 l.The res ulting 8-dimensional design space (each variable is a dimension of the design s pace) is composed of 5.4 × 10 7 points, each representing a possible design solution for the multi-energy system.

Cost functions
Cost functions representing the investment and maintenance costs for the variable system components with reference to the defined set of o ptimization variables (Table 1), are reported in Table 2.
An example of the process for determining cost functions is reported in Figure 2, where the total investment cost function for the integration thermal storage was derived by interpolation of real cost val ues for d ifferent storage siz es reported in the price lists (official price list of Valle d'Aosta Region, 2018).The so-called size-dependent costs are related to the size of the storage tank, fixed costs r efer to the set of valves and pipes that are n ecessary for h ydronic system connections, while installation costs are se t to 25% of the sum of fixed a nd size-dependent costs, as suggested by th e price lists.The linear cost func tion is derived as the best fit to the analysed costs (R 2 = 0.97).Regarding thermal storages, the cost analysis based on price lists, and confirmed by manufacturers, has highlighted that only th e volume and the number of se rpentines (two heat exchanger serpentines are included in the integration hot storage, one is included in the grey water storage) have direct impact on the investment cost, if the heat exchange surface area of serpentines is maintained within feasible ranges.Therefore, cost functions of ther mal storages only depend on optimization variables related to storage volume (variables STVint, STVw) and, in case of integration hot storage, on the number of solar collectors, because of th e interdependency between STVint and NColl.That is why the cost function in Figure 3 is defined considering the storage volume as the x independent variable.In the model, as reported in Table 2, such volume is determined as the product of ST Vint and the total area of solar collectors, obtained by multiplying NColl by the single collector area (2 m 2 ), with a minimum value of 500 l set as a constraint.

Objective functions
According to Standard EN 15459 for economic evaluation of energy systems in building s (CEN 2017), the global cost has been considered as the primary objective of this stud y, as follows: where C G (τ, P) represents the global cost r elatively to the starting year τ 0 , considering a number τ of years as the calculation period and the defined set of design variables P, CI is the initial investment cost, C a,i (j) is the annual cost for component j at the year i (including running costs an d periodic or replacement costs), R d (i) is the discount rate for year i, V f,τ (j) is the final value of component j at the end of the calculation period (relatively to the starting year τ 0 ).The calculation period of 30 years and the financial parameters were initially set according to EU guidelines for performing cost-optimal analysis (European Commission 2012) following Directive 2010/31/EU (Market interest rate 4.5%, inflation rate 2%, real interest rate 2.5%, rate o f evolution of energy price 2%).Lower real interest rates (R R = 1.5% and R R = 0.5%) deriving from different combinations of market and inflation rates reflecting the EU macroeconomic context of the last decade were also considered to perform a sensitivity analysis to financial parameters.
With respect to the o ther macroeconomic parameters, energy prices may be subject to much higher variabi lity in the very sh ort and med ium term, due to the unsta ble availability of energy sources and their strict relation with social and political conte xts.Therefore, following official statistics about energy prices in the European Union in view of the current international political context, which is likely to have a very important impact on the current and future energy prices (https://ec.europa.eu/eurostat/web/energy/database), a dedicated sensitivity analysis was  conducted on the rate of development of e nergy prices.In particular, the rate of development of electricity price R e was considered and the sensitivit y of results to the increase of such rate up to R e = 10% was investigated (in the very last years many countries have been experiencing such average increase of energy prices, including Italy).
It is importan t to note that this global co st objective function is directly influenced by the energy performance itself, as lower energy use is related to lower running costs, included in the global cost computation.H owever, energy costs are associated to purchased energy from the grid but may not be representative of the e fficiency of the system.Therefore, the seasonal performance factor of the system was also evaluated to better understand the energy performance in terms of system efficiency.It was defined, as used for multipurpose systems, the ratio of the total useful energy output to the total energy expense of a system (Zottl and Nordman 2012;Biglia et al. 2021) as follows: where Q h , Q v and Q DHW are the useful energy outputs for space heating, ventilation and DHW, wh ile E el,HP and E el,aux are the electrical energy inputs for the heat pump operation and the auxiliary systems (circulation pumps).Given these two objectives, a multi-objective optimization function can be written as follows: where MOF SAGHP Î [0, 1].Weights w 1 and w 2 can be assigned to the two o bjectives according to user-de fined priority, provided that w 1 + w 2 = 1.In this study w 1 = w 2 = 0.5.

System modelling results
The adopted model setup resulted to be able to correctly simulate the system operation according to the initial design specifications.Figure 4 reports several model outputs for four example days of operation (January 14 th to 17 th ).In details, describing Figure 4 from bottom to top, it is possible to see the modelled dynamic operation of the different parts of the system by means of the dynamic trend of several simulated variables.
Regarding the external envir onment, it is shown that the external dry bulb air tem perature (T_ext [°C]) is quite rigid and often below 0 °C, demonstrating the need for heat pump sources with higher temperatures than air, s uch as the geothermal source, of which the average temperature is shown to be around 6 °C [T_ave_geo [°C]).It is also shown that the solar radiation (Solar_rad [W/m 2 ]) is available during opening hours when higher thermal loads occur (Q_tot [kW]), showing the opportunity to use the solar source for further increasing the source temperature of the heat pump.The dynamics of thermal loads for the different uses (Q_AHU -heating load for air handling units operation [kW th ], Q_radiant floor -heating load fo r radiant floor operation [kW th ]; Q_DHW -heating load for domestic hot water [kW th ]) show the accuracy of the model according to the expected operating conditions (smaller heating loads at night are required for the radiant floor operating in antifreeze conditions).
Thermal loads are also reflected in the dynamics of temperatures in the three load-side storages where the variable T_storage_highT represents the average fluid tem perature in the high-temperature storage dedicated to the air handling units loop, T_storage_lowT is the average fluid temperature in the low-temperature storage dedicated to the radiant floor loop and T_storage_DHW is the average temperature in the storage dedicated to domestic hot water.
Such oscillating dynamics of storages show the correct simulation of the degree-minute-based control logic for the operation of t he heat pump.Its intermittent operation is driven by the entity and speed of temperature decrease in the load-side storages and is shown through the dynamics of the electric power demand of the double-stage heat pump (P_ele_HP -[kW el ]) and of thermal capacities absorbed at the evaporator and delivered at the condensator (Q_source_HP and Q_load_HP, respectively [kW th ] -peaks represent the call for second stage of operation).
Regarding the operating con ditions of the storages at the source side, the graph at the top of Fig ure 4 shows in details the different temperature levels at the different points of the loop and the contribution of the solar and the greywater sources to increase the t emperature in the integration hot storage (see system scheme in Figure 1).In details: -T_in_geo -t emperature of the fluid entering the geothermal field [°C]; -T_in_STint -temperature of the fluid coming from the geothermal field and entering the integration hot storage, waiting for gaining additional heat from th e solar and grey water sources [°C]; -T_out_STint -temperature of the flu id at the outlet of the integration hot storage, ready for transferring heat to the heat pump evaporator [°C]; -T_ave_STgw -average wat er temperature in the grey water storage [°C]; -Q_fromHXsol -heat transferred to the integration hot storage from the solar loop [kW]; Fig. 4 Four sample days of dynamic operation of the modelled system in its initial configuration -Q_fromHXgw -heat transferred to the integration hot storage from the grey water loop [kW].These results have been co mpared to previous studies where similar systems were modelled in TRNSYS, as described in Section 2. Despite the d ifferences due to the different boundary conditions and component sizing referring to each particular case study, the obtained sim ulation outputs are co nsistent with expectations ba sed on system design features and with si mulation outputs resulting in other studies reported in the literature.In particular, it is shown that the oscillating tempe rature levels of the fluids entering and exiting the borehole heat exchangers and entering in the heat pump follows the on-o ff operation of the heat pump and related call for storage discharge, similarly to a previous study reporting outputs of a si milar TRNSYS model (Li et a l. 2018).The same study reports a similar behavior of temperatures in water storag es.Coherently with the outputs reported in a previous study regarding the operation of a SAGHP in a cold climate (Rad et al. 2013), the average te mperature increase is aro und 3 °C in the initial configuration, where a recent study demonstrating that the output of such modelling approach may lead to a small underestimation of te mperature of t he fluid exi ting the ground heat exchanger (0.05 °C).However, as mentioned, the reults reported in Figure 4 refer to sample days for the system operation in its initi al design configuration (blue point in Fig ure 5) and should be interpreted as t ypical trends for the all involved and interrelated output variables, rather than absolute results.These are shown as an example of the correct operation and good accuracy of the mod el, but these are subject to variation for the optimiz ed system configurations, which will be show n in the next result subsection.

Optimization results
The entire set of evaluated design alternatives, report ed in Figure 5, wer e obtained within the four single-objective optimization runs to m inimize and max imize the two objectives (global cost and seasonal performance factor, Eqs. ( 1) and ( 2)) and the final multi-ob jective optimization run (Eq.( 3)).Each point of Figure 5 re presents a different design alternative, reported in a gra ph where the global cost is reported on the vertical axis and the SPF is reported on the horizontal axis.Further than the point minimizing the multi-objective function (MOBopt, orange dot), o ther notable points were highlighted, such as the minimum and the maximum achievable S PF (green and red dots for SPFmin and SPFmax, respectively) and the minimum and the maximum achievable glo bal cost (violet and light blue dots for GCmin and GCmax, respectively).The initial design configuration (INI) is r epresented in blue.The d ifferent sets of design variable values for the highl ighted notable points are reported in Table 3, together with the related values of global cost, SPF and non-renewable energy expense (E el,tot , that is the denominator of Eq. ( 2)).
Results show that, within th e design space delimited by the defined set of paramete rs, the global cost may vary between 179023 € and 52342 3 € (−13.7% a nd +152% with respect to the initial design solution), while the S PF may vary between 3.66 and 3.88 (−0.5% and + 5% with respect to the initial design solution) .The minimum global cost is achieved with SPF that is slightly higher than the minimum, while the max imum SPF is only achievable with higher increase of global cost.However, the min imization of th e multi-objective optimization function shows that it is possible to identify a solution leading to nearly optimal SPF C C Fig. 5 Global cost vs SPF of all evaluated design alternatives and notable points but significant reduction of global cost (−34%).
The analysis o f the res ulting optimal values of design variables in the different notable points shows that the global cost is highly influenced b y the investm ent cost, w hich greatly increases when the number and the deep of the borehole increase, without causing a significant increase in energy performance and a consequent significant reduction of operational costs that is able to impact the overall global cost calculation.From the e nergy performance point of view, it is show n that the initial design co nfiguration INI has almost the same low performance level of the resulting cost-optimal solution, but even with a high number of solar field and significant in crease of the geothermal field the achievable increase of SPF is around 6%.
However, such performance increase can be appraised by studying the dynamics of the different system variables.Figure 6 shows the d ynamics of the variables related to the source-side loop in the SPFopt and the MOBopt point for the same sam ple days alrea dy reported in Figure 4.The temperatures of the fluid at the different loop points is reported (T_in_geo, T_in_STint, T_out_STint, T_ave_STgw), showing the ability of the system setup to effectively increase the thermal levels of the fluid thanks to the different heat sources and therefore increase the overall system performance.
If compared to the trends in Figure 4, it is shown that, in both optimized system configurations SPFopt and MOBopt, the temperature of the fluid exiting the integration storage (T_out_STint) is increased of about 3.5°C on average with respect to th e INI configuration (Figure 3).However, the temperature dynamics is different due to greater exploitation of the geothermal field in th e SPFopt configuration, while the MOBopt relies more on the exploitation of a greater solar field, as also shown by the variable Q_fromHXsol in Figure 6.This is also coherent with the resulting optimal value of optimization variabl es reported in Table 3.Both SPFopt and MOBopt are also able to make the contribution of the heat recovery from waste greywater increase with respect to the INI configuration.

Analysis of the optimum neighborhood
Beyond the a bsolute results in terms of en ergy and cost performance (the dependent variables), the relative relationship between one and the other design alternative should be considered, with t he resulting optimal values of the selected optimization variables to be implemented in the real design (the independent variables).In this context, the analysis of the points located in the neighborhood of notable and/or optimal points may give important information about the robustness of the resulting o ptimal design solutions against unce rtainty of modelling inp uts and outputs, as also emerging from other applications of the EDeSSOpt methodology (Ferrara et al. 2018(Ferrara et al. , 2020)).
Therefore, the set of optimization variable s related to different design alternatives falling in the the neighborhoods of the GCmin and the MOBopt points were studied.
The neighborhoods where defined consider ing a 10% maximum distance from opt imum values of the relevant objective functions.The neig hborhood of GCmin resulted to be composed of 216 points, with GC lower than 196925 € and SPF ranging between 3 .66and 3.70, while 440 points were included in the neighborhood of MOBopt, where the 10% deviation from minimum MOB value led to identify a part of the Pareto front where SPF ranges between 3.81 and 3.87 for global cost values ra nging between 298283 € and 361080 € (Figure 7).
The graphs in Figure 8 and Figure 9 show the distribution of the frequency of values of optimization variables occurring in the points composing the neighborhoods.The fact that, for an optimization variable, most points in the neighborhood have the same values as the optimal point means that the variable is significant -if that variable is not optim ized in a design alternative, it is like ly to fall outside the optimu m neighborhood.
Regarding the neighborhood of G Cmin (details in Figure 8), it is shown that the main design variables have high level of robustness as most design alternatives in the neighborhood have a low number of solar collectors and a low number of boreholes (depth appears as not significant).Regarding storages, it is shown that the volume of the integration hot storage (driven by variable STVint) should stay low, as a consequence of the low number of solar collectors.On the other hand, it is shown that the volume of the greywater storage (driven by variable STVgw) is not a driver for the optimality.The proper sizing of the integration hot storage should have priority over the size of the greywater storage in ord er to ensure the cost-optimality of the re al design.
Different results are shown in Figure 9, where the analysis of values of o ptimization variables in th e neighborhood of MOBopt is reported.In fact, only the NBor variable presents a high level of robustness.It is shown that the other variables may assume almost all values in the defined variability range while leading to the neighborhood of the minimum value of the multi-objective function.In other words, it is shown that many different design configurations may lead to similar results and therefore the final decision may be made following other criteria.configurations of optimal and notable points do not change due to variation of real interest rate (the same optimizations that were initially run with R R of 2.5% were then run with R R of 1.5% and 0.5%).These results demonstrate a high stability of the obtained optimal design solutions with respect to variation of financial para meters.However, lower real interest rates may lead to increase the global cost values, for which the present study quantifies up to 16% (for MOBopt solution) or 25% (for GCmin solution).

Sensitivity analysis
Figure 11 reports the results obtained from the analysis of sensitivity to the rate of development of e nergy prices R e when the other financial parameters are set to their initial value (R R set to 2.5% following European guidelines).Also in this case, th e relative position of one point with respect to the others is not subject to substantial variations, meaning that design co nfigurations of optimal and notable points and their related neighborhoods are stable with respect to the rate of development of energy prices (the same optimizations that were initially run with R e of 2% were then run with R e of 6% and 10%).
As in previous case, these results confirm a high level of robustness of the obtained o ptimal design solutions with respect to variation of financial parameters.This means that it is possible to design and implement a system configuration that is op timized in the current economic context and it is likely to stay optimal even in case of quite high variation of financial parameters.However, the higher development shows very high increase of global cost values, for which the present study quantifies up to 228% (for GCmin s olution) or 1 24% (for MOBopt solution).Therefore, in case of high development of energy price rates a much lower relative variation of glo bal cost values with respect to initial design solution can be appraised.In fact, in case of R e = 10%, the global cost may vary between 588534 € and 912244 € in absolute terms, that is −4% an d +48% in relative terms with respect to the initial design solution, with the MOBopt solution leading to increase the global cost by only 14% with respect to the initial design configuration (it was +51% with the basic assum ption of the rate of development of energy prices).
In an economic context when the rate of development of energy price is expected to be subject to very high increase, this approach may lead to change the final design choice in favour of a more energ y efficient design configuration, despite it is clear that, within the defined set of constraints, investment costs will account for the higher weight on the global cost even if the weight of energy cost increases.

Conclusions
The design of a real SAGHP system for a building in the Alps region was modelled in detail and optimized from both energy performance and cost points of view.The simulation model resulted to be effective to support the understanding of the dynamics of the differ ent parts of th e multi-energy systems and a ppraising the impact of desig n variation on the behavior of each component and of the overall system, as well as the mut ual relationship between the solar field , the geothermal field and the thermal storages.In particular, a new dedicated simulation-bas ed optimization framework was created w ith a ta ilored set of optim ization variables, providing accurate definition of modeling constraints occurring within the differen t variables rela ted to ther mal storages.With respect to the initial des ign solution, the study demonstrated a potential reduction of global cost of 13.7% for t he same energ y performance, or a potent ial increase of seasonal performance factor of about 6% at the expense of very high increase of global cost.
This means that such a system configuration can reach valuable performance (similar to that achieved from similar systems in other con texts), if properly optimized and integrated with available renewable sources, also in extreme context of environmental conditions, therefore demonstrating the wide applicability of the heat pump technology towards worldwide decarbonization.Moreover, the performed sensitivity analysis show ed high level of robustness of th e resulting optimal solution with respect to financial parameters and even with the curr ent high development of energy prices, demonstrating the reliability of the performed technical studies even in context of financial and economic uncertainty, while highlighting the need for incentive policies that are able to increase the weight of energy costs in the calculated global cost and the refore increase financial feasibility of this type of system.
In general, beyond the numerical results tha t can drive design of such systems in similar contexts, the paper highlights the potential of dynamic simu lation and optimization to support the complex design of multi-source multi-energy systems and shows the possibility to identify trade-off design solutions driven by KPIs t hat are representative of technical and financial objectives.The presented approach and original setup of simulation-based optimization framework for th e thermal storage design para meters can be replicated for modelling and optimizing the other s imilar multi-energy multi-source systems with different set of external constraints.
Further work should be perf ormed to incl ude control strategies in the set of optimi zation variables, especially regarding control of thermal storages, in o rder to e xplore further possibilities for improving energy performance without causing investment cost increases.

Declaration of competing interest
The authors have no competing interests to declare that are relevant to the content of this article.

Fig. 1
Fig. 1 System scheme in the initial design configuration

Fig. 2
Fig. 2 System modelling network in the TRNSYS Simulation Studio interface

Fig. 3
Fig. 3 Example of cost function for a thermal storage tank

Fig. 6
Fig. 6 Dynamics of source-side variables in four sample days for the SPFopt and MOBopt optimized system configurations (T_in_geotemperature of the fluid entering the geothermal field [°C]; T_in_STint -temperature of the entering the integration hot storage[°C]; T_out_STint -temperature of the fluid at the outlet of the integration hot storage [°C]; T_ave_STgw -average water temperature in the grey water storage [°C]; Q_fromHXsol -heat transferred to the integration hot storage from the solar loop [kW]; Q_fromHXgw -heat transferred to the integration hot storage from the grey water loop [kW], Q_tot -total heating loads)

Figure 10
Figure 10 reports the res ults obtained fro m the analysis of sensitivity to financial p arameters.As shown, design

Fig. 9
Fig. 9 Analysis of the values of optimization variables in t he neighborhood of MOBopt, where red dots indicate the optimal variable values in the MOBopt point

Fig. 11
Fig. 11 All evaluated design alternatives -sensitivity to the rate of development of energy prices

Table 1
Design optimization variables

Table 2
Investment cost, maintenance cost and lifespan of the system components as a function of optimization variables

Table 3
Design configuration and related objective function values in optimal and notable points