1"""
2This file contains an example on how GHEtool can be used to size a borefield
3using a combination of active and passive cooling.
4This example is based on the work of Coninx and De Nies, 2021.
5Coninx, M., De Nies, J. (2022). Cost-efficient Cooling of Buildings by means of Borefields
6with Active and Passive Cooling. Master thesis, Department of Mechanical Engineering, KU Leuven, Belgium.
7It is also published as: Coninx, M., De Nies, J., Hermans, L., Peere, W., Boydens, W., Helsen, L. (2024).
8Cost-efficient cooling of buildings by means of geothermal borefields with active and passive cooling.
9Applied Energy, 355, Art. No. 122261, https://doi.org/10.1016/j.apenergy.2023.122261.
10"""
11import copy
12
13from GHEtool import Borefield, GroundConstantTemperature, HourlyGeothermalLoadMultiYear
14
15import pandas as pd
16import numpy as np
17import matplotlib.pyplot as plt
18import optuna
19
20
21def active_passive_cooling(location='Active_passive_example.csv'):
22 # load data
23 columnNames = ['HeatingSpace', 'HeatingAHU', 'CoolingSpace', 'CoolingAHU']
24 df = pd.read_csv(location, names=columnNames, header=0)
25 heating_data = df.HeatingSpace + df.HeatingAHU
26 cooling_data = df.CoolingSpace + df.CoolingAHU
27
28 # variable COP and EER data
29 COP = [0.122, 4.365] # ax+b
30 EER = [-0.3916, 17.901] # ax+b
31 threshold_active_cooling = 16
32
33 # set simulation period
34 SIMULATION_PERIOD: int = 50
35 heating_building: np.ndarray = np.tile(np.array(heating_data), SIMULATION_PERIOD)
36 cooling_building: np.ndarray = np.tile(np.array(cooling_data), SIMULATION_PERIOD)
37
38 def update_load_COP(temp_profile: np.ndarray,
39 COP: np.ndarray,
40 load_profile: np.ndarray) -> np.ndarray:
41 """
42 This function updates the geothermal load for heating based on a variable COP
43
44 Parameters
45 ----------
46 temp_profile : np.ndarray
47 Temperature profile of the fluid
48 COP : np.ndarray
49 Variable COP i.f.o. temperature
50 load_profile : np.ndarray
51 Heating load of the building
52
53 Returns
54 -------
55 Geothermal heating load : np.ndarray
56 """
57 COP_array = temp_profile * COP[0] + COP[1]
58 return load_profile * (1 - 1 / COP_array)
59
60 def update_load_EER(temp_profile: np.ndarray,
61 EER: np.ndarray,
62 threshold_active_cooling: float,
63 load_profile: np.ndarray) -> np.ndarray:
64 """
65 This function updates the geothermal load for cooling based on a threshold for active/passive cooling,
66 and a variable EER.
67
68 Parameters
69 ----------
70 temp_profile : np.ndarray
71 Temperature profile of the fluid
72 EER : np.ndarray
73 Variable EER i.f.o. temperature
74 threshold_active_cooling : float
75 Threshold of the temperature above which active cooling is needed
76 load_profile : np.ndarray
77 Cooling load of the building
78
79 Returns
80 -------
81 Geothermal cooling load : np.ndarray
82 """
83 EER_array = temp_profile * EER[0] + EER[1]
84 passive: np.ndarray = temp_profile <= threshold_active_cooling
85 active = np.invert(passive)
86 return active * load_profile * (1 + 1 / EER_array) + passive * load_profile
87
88 costs = {"C_elec": 0.2159, # electricity cost (EUR/kWh)
89 "C_borefield": 35, # inv cost per m borefield (EUR/m)
90 "DR": 0.0011, # discount rate(-)
91 "sim_period": SIMULATION_PERIOD}
92
93 def calculate_costs(borefield: Borefield, heating_building: np.ndarray, heating_geothermal: np.ndarray,
94 cooling_building: np.ndarray, cooling_geothermal: np.ndarray, costs: dict) -> tuple:
95 """
96 This function calculates the relevant costs for the borefield.
97
98 Parameters
99 ----------
100 borefield : Borefield
101 Borefield object
102 heating_building : np.ndarray
103 Heating demand for the building
104 heating_geothermal : np.ndarray
105 Heating demand coming from the ground
106 cooling_building : np.ndarray
107 Cooling demand for the building
108 cooling_geothermal : np.ndarray
109 Cooling demand coming from the ground
110 costs : dict
111 Dictionary with investment cost for borefield/m, electricity cost, annual discount rate
112
113 Returns
114 -------
115 investment cost borefield, operational cost heating, operational cost cooling, total operational cost:
116 float, float, float, float
117 """
118 # calculate investment cost
119 investment_borefield = costs["C_borefield"] * borefield.H * borefield.number_of_boreholes
120
121 # calculate working costs
122 elec_heating = heating_building - heating_geothermal
123 elec_cooling = cooling_geothermal - cooling_building
124
125 discounted_cooling_cost = []
126 discounted_heating_cost = []
127 for i in range(SIMULATION_PERIOD):
128 tempc = costs["C_elec"] * (elec_cooling[730 * 12 * i:730 * 12 * (i + 1)])
129 tempc = tempc * (1 / (1 + costs["DR"])) ** (i + 1)
130
131 temph = costs["C_elec"] * (elec_heating[730 * 12 * i:730 * 12 * (i + 1)])
132 temph = temph * (1 / (1 + costs["DR"])) ** (i + 1)
133 discounted_cooling_cost.append(tempc)
134 discounted_heating_cost.append(temph)
135 cost_cooling = np.sum(discounted_cooling_cost)
136 cost_heating = np.sum(discounted_heating_cost)
137
138 return investment_borefield, cost_heating, cost_cooling, cost_heating + cost_cooling
139
140 borefield = Borefield()
141 borefield.simulation_period = SIMULATION_PERIOD
142 borefield.set_max_avg_fluid_temperature(17)
143
144 borefield.create_rectangular_borefield(12, 12, 6, 6, 100)
145 borefield.set_ground_parameters(GroundConstantTemperature(2.1, 11))
146 borefield.Rb = 0.12
147
148 ### PASSIVE COOLING
149 depths = [0.9, 0]
150
151 # set initial loads
152 cooling_ground = cooling_building.copy()
153 heating_ground = heating_building.copy()
154
155 while abs(depths[0] - depths[1]) > 0.1:
156 # set loads
157 load = HourlyGeothermalLoadMultiYear()
158 load.hourly_extraction_load = heating_ground
159 load.hourly_injection_load = cooling_ground
160 borefield.load = load
161
162 # size borefield
163 depth_passive = borefield.size_L4()
164 depths.insert(0, depth_passive)
165
166 # get temperature profile
167 temp_profile = borefield.results.peak_extraction
168
169 # recalculate heating load
170 heating_ground = update_load_COP(temp_profile, COP, heating_building)
171
172 ### ACTIVE COOLING
173 depths = [0.9, 0]
174
175 # set initial loads
176 cooling_ground = cooling_building.copy()
177 heating_ground = heating_building.copy()
178
179 borefield.set_max_avg_fluid_temperature(25)
180 borefield.gfunction_calculation_object.store_previous_values = False
181 while abs(depths[0] - depths[1]) > 0.1:
182 # set loads
183 load = HourlyGeothermalLoadMultiYear()
184 load.hourly_extraction_load = heating_ground
185 load.hourly_injection_load = cooling_ground
186 borefield.load = load
187
188 # size borefield
189 depth_active = borefield.size_L4()
190 depths.insert(0, depth_active)
191
192 # get temperature profile
193 temp_profile = borefield.results.peak_extraction
194
195 # recalculate heating load
196 heating_ground = update_load_COP(temp_profile, COP, heating_building)
197 cooling_ground = update_load_EER(temp_profile, EER, threshold_active_cooling, cooling_building)
198
199 ### RUN OPTIMISATION
200
201 # initialise parameters
202 operational_costs = []
203 operational_costs_cooling = []
204 operational_costs_heating = []
205 investment_costs = []
206 total_costs = []
207 depths = []
208
209 def f(depth: float) -> float:
210 """
211 Optimisation function.
212
213 Parameters
214 ----------
215 depth : float
216 Depth of the borefield in meters
217
218 Returns
219 -------
220 total_cost : float
221 """
222 borefield._update_borefield_depth(depth)
223 borefield.H = depth
224 depths.append(depth)
225
226 # initialise
227 heating_ground = heating_building.copy()
228 cooling_ground = cooling_building.copy()
229
230 heating_ground_prev = np.zeros(len(heating_ground))
231 cooling_ground_prev = np.zeros(len(cooling_ground))
232
233 # iterate until convergence in the load
234 while np.sum(cooling_ground + heating_ground - heating_ground_prev - cooling_ground_prev) > 100:
235 # set loads
236 load = HourlyGeothermalLoadMultiYear()
237 load.hourly_extraction_load = heating_ground
238 load.hourly_injection_load = cooling_ground
239 borefield.load = load
240
241 # get temperature profile
242 borefield.calculate_temperatures(depth, hourly=True)
243 temp_profile = borefield.results.peak_extraction
244
245 # set previous loads
246 heating_ground_prev = heating_ground.copy()
247 cooling_ground_prev = cooling_ground.copy()
248
249 # recalculate heating load
250 heating_ground = update_load_COP(temp_profile, COP, heating_building)
251 cooling_ground = update_load_EER(temp_profile, EER, 16, cooling_building)
252
253 # calculate costs
254 investment, cost_heating, cost_cooling, operational_cost = calculate_costs(borefield,
255 heating_building, heating_ground,
256 cooling_building, cooling_ground,
257 costs)
258 total_costs.append(investment + operational_cost)
259 operational_costs.append(operational_cost)
260 operational_costs_cooling.append(cost_cooling)
261 operational_costs_heating.append(cost_heating)
262 investment_costs.append(investment)
263 return investment + operational_cost
264
265 # add boundaries to figure
266 MIN_BOUNDARY = depth_active
267 MAX_BOUNDARY = depth_passive
268
269 def objective(trial: optuna.Trial):
270 depth = trial.suggest_float('depth', MIN_BOUNDARY, MAX_BOUNDARY)
271 return f(depth)
272
273 study = optuna.create_study()
274 study.optimize(objective, n_trials=100)
275
276 depths_sorted = copy.deepcopy(depths)
277 depths_sorted.sort()
278 depths_old_new = {}
279 for idx, depth in enumerate(depths_sorted):
280 depths_old_new[idx] = depths.index(depth)
281
282 # plot figures
283 fig = plt.figure()
284 ax1 = fig.add_subplot(111)
285 ax1.plot(depths_sorted, [total_costs[depths_old_new[idx]] / 1000 for idx, _ in enumerate(depths_sorted)],
286 marker='o', label="TC")
287 ax1.plot(depths_sorted, [investment_costs[depths_old_new[idx]] / 1000 for idx, _ in enumerate(depths_sorted)],
288 marker='o', label="IC")
289 ax1.plot(depths_sorted, [operational_costs[depths_old_new[idx]] / 1000 for idx, _ in enumerate(depths_sorted)],
290 marker='o', label="OC")
291 ax1.plot(depths_sorted,
292 [operational_costs_cooling[depths_old_new[idx]] / 1000 for idx, _ in enumerate(depths_sorted)], marker='o',
293 label="OCc")
294 ax1.plot(depths_sorted,
295 [operational_costs_heating[depths_old_new[idx]] / 1000 for idx, _ in enumerate(depths_sorted)], marker='o',
296 label="OCh")
297 ax1.set_xlabel(r'Depth (m)', fontsize=14)
298 ax1.set_ylabel(r'Costs ($k€$)', fontsize=14)
299 ax1.legend(loc='lower left', ncol=3)
300 ax1.tick_params(labelsize=14)
301 plt.show()
302
303
304if __name__ == "__main__": # pragma: no cover
305 active_passive_cooling()