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"""
11
12from GHEtool import Borefield, GroundConstantTemperature, HourlyGeothermalLoadMultiYear
13
14import pandas as pd
15import numpy as np
16import matplotlib.pyplot as plt
17from skopt import gp_minimize
18
19
20def active_passive_cooling(location='Active_passive_example.csv'):
21
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 = [-3.916, 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
61 def update_load_EER(temp_profile: np.ndarray,
62 EER: np.ndarray,
63 threshold_active_cooling: float,
64 load_profile: np.ndarray) -> np.ndarray:
65 """
66 This function updates the geothermal load for cooling based on a threshold for active/passive cooling,
67 and a variable EER.
68
69 Parameters
70 ----------
71 temp_profile : np.ndarray
72 Temperature profile of the fluid
73 EER : np.ndarray
74 Variable EER i.f.o. temperature
75 threshold_active_cooling : float
76 Threshold of the temperature above which active cooling is needed
77 load_profile : np.ndarray
78 Cooling load of the building
79
80 Returns
81 -------
82 Geothermal cooling load : np.ndarray
83 """
84 EER_array = temp_profile * EER[0] + EER[1]
85 passive: np.ndarray = temp_profile < threshold_active_cooling
86 active = np.invert(passive)
87 return active * load_profile * (1 + 1/EER_array) + passive * load_profile
88
89
90 costs = {"C_elec": 0.2159, # electricity cost (EUR/kWh)
91 "C_borefield": 35, # inv cost per m borefield (EUR/m)
92 "DR": 0.0011, # discount rate(-)
93 "sim_period": SIMULATION_PERIOD}
94
95
96 def calculate_costs(borefield: Borefield, heating_building: np.ndarray, heating_geothermal: np.ndarray,
97 cooling_building: np.ndarray, cooling_geothermal: np.ndarray, costs: dict) -> tuple:
98 """
99 This function calculates the relevant costs for the borefield.
100
101 Parameters
102 ----------
103 borefield : Borefield
104 Borefield object
105 heating_building : np.ndarray
106 Heating demand for the building
107 heating_geothermal : np.ndarray
108 Heating demand coming from the ground
109 cooling_building : np.ndarray
110 Cooling demand for the building
111 cooling_geothermal : np.ndarray
112 Cooling demand coming from the ground
113 costs : dict
114 Dictionary with investment cost for borefield/m, electricity cost, annual discount rate
115
116 Returns
117 -------
118 investment cost borefield, operational cost heating, operational cost cooling, total operational cost:
119 float, float, float, float
120 """
121 # calculate investment cost
122 investment_borefield = costs["C_borefield"] * borefield.H * borefield.number_of_boreholes
123
124 # calculate working costs
125 elec_heating = heating_building - heating_geothermal
126 elec_cooling = cooling_geothermal - cooling_building
127
128 discounted_cooling_cost = []
129 discounted_heating_cost = []
130 for i in range(SIMULATION_PERIOD):
131 tempc = costs["C_elec"] * (elec_cooling[730 * 12 * i:730 * 12 * (i + 1)])
132 tempc = tempc * (1 / (1 + costs["DR"])) ** (i + 1)
133
134 temph = costs["C_elec"] * (elec_heating[730 * 12 * i:730 * 12 * (i + 1)])
135 temph = temph * (1 / (1 + costs["DR"])) ** (i + 1)
136 discounted_cooling_cost.append(tempc)
137 discounted_heating_cost.append(temph)
138 cost_cooling = np.sum(discounted_cooling_cost)
139 cost_heating = np.sum(discounted_heating_cost)
140
141 return investment_borefield, cost_heating, cost_cooling, cost_heating+cost_cooling
142
143 borefield = Borefield()
144 borefield.simulation_period = SIMULATION_PERIOD
145 borefield.set_max_avg_fluid_temperature(17)
146
147 borefield.create_rectangular_borefield(12, 12, 6, 6, 100)
148 borefield.set_ground_parameters(GroundConstantTemperature(2.1, 11))
149 borefield.Rb = 0.12
150
151 ### PASSIVE COOLING
152 depths = [0.9, 0]
153
154 # set initial loads
155 cooling_ground = cooling_building.copy()
156 heating_ground = heating_building.copy()
157
158 while abs(depths[0] - depths[1]) > 0.1:
159 # set loads
160 load = HourlyGeothermalLoadMultiYear()
161 load.hourly_heating_load = heating_ground
162 load.hourly_cooling_load = cooling_ground
163 borefield.load = load
164
165 # size borefield
166 depth_passive = borefield.size_L4()
167 depths.insert(0, depth_passive)
168
169 # get temperature profile
170 temp_profile = borefield.results.peak_heating
171
172 # recalculate heating load
173 heating_ground = update_load_COP(temp_profile, COP, heating_building)
174
175
176 ### ACTIVE COOLING
177 depths = [0.9, 0]
178
179 # set initial loads
180 cooling_ground = cooling_building.copy()
181 heating_ground = heating_building.copy()
182
183 borefield.set_max_avg_fluid_temperature(25)
184 while abs(depths[0] - depths[1]) > 0.1:
185 # set loads
186 load = HourlyGeothermalLoadMultiYear()
187 load.hourly_heating_load = heating_ground
188 load.hourly_cooling_load = cooling_ground
189 borefield.load = load
190
191 # size borefield
192 depth_active = borefield.size_L4()
193 depths.insert(0, depth_active)
194
195 # get temperature profile
196 temp_profile = borefield.results.peak_heating
197
198 # recalculate heating load
199 heating_ground = update_load_COP(temp_profile, COP, heating_building)
200 cooling_ground = update_load_EER(temp_profile, EER, 16, cooling_building)
201
202
203 ### RUN OPTIMISATION
204
205 # initialise parameters
206 operational_costs = []
207 operational_costs_cooling = []
208 operational_costs_heating = []
209 investment_costs = []
210 total_costs = []
211 depths = []
212
213
214 def f(depth: list) -> float:
215 """
216 Optimisation function.
217
218 Parameters
219 ----------
220 depth : list
221 List with one element: the depth of the borefield in mm
222
223 Returns
224 -------
225 total_cost : float
226 """
227 # convert to meters
228 depth = depth[0] / 1000
229 borefield._update_borefield_depth(depth)
230 borefield.H = depth
231 depths.append(depth)
232
233 # initialise
234 heating_ground = heating_building.copy()
235 cooling_ground = cooling_building.copy()
236
237 heating_ground_prev = np.zeros(len(heating_ground))
238 cooling_ground_prev = np.zeros(len(cooling_ground))
239
240 # iterate until convergence in the load
241 while np.sum(cooling_ground + heating_ground - heating_ground_prev - cooling_ground_prev) > 100:
242 # set loads
243 load = HourlyGeothermalLoadMultiYear()
244 load.hourly_heating_load = heating_ground
245 load.hourly_cooling_load = cooling_ground
246 borefield.load = load
247
248 # get temperature profile
249 borefield.calculate_temperatures(depth, hourly=True)
250 temp_profile = borefield.results.peak_heating
251
252 # set previous loads
253 heating_ground_prev = heating_ground.copy()
254 cooling_ground_prev = cooling_ground.copy()
255
256 # recalculate heating load
257 heating_ground = update_load_COP(temp_profile, COP, heating_building)
258 cooling_ground = update_load_EER(temp_profile, EER, 16, cooling_building)
259
260 # calculate costs
261 investment, cost_heating, cost_cooling, operational_cost = calculate_costs(borefield,
262 heating_building, heating_ground,
263 cooling_building, cooling_ground,
264 costs)
265 total_costs.append(investment + operational_cost)
266 operational_costs.append(operational_cost)
267 operational_costs_cooling.append(cost_cooling)
268 operational_costs_heating.append(cost_heating)
269 investment_costs.append(investment)
270 return investment + operational_cost
271
272 # add boundaries to figure
273 # multiply with 1000 for numerical stability
274 f([depth_active * 10 ** 3])
275 f([depth_passive * 10 ** 3])
276
277 res = gp_minimize(f, # the function to minimize
278 [(depth_active * 10 ** 3, depth_passive * 10 ** 3)], # the bounds on each dimension of x
279 acq_func="EI", # the acquisition function
280 n_calls=30, # the number of evaluations of f
281 initial_point_generator="lhs",
282 n_random_starts=15, # the number of random initialization points
283 # noise=0, # the noise level (optional)
284 random_state=1234) # the random seed
285
286 # plot figures
287 fig = plt.figure()
288 ax1 = fig.add_subplot(111)
289 ax1.plot(depths, [x/1000 for x in total_costs], marker = 'o', label = "TC")
290 ax1.plot(depths, [x/1000 for x in investment_costs], marker = 'o', label="IC")
291 ax1.plot(depths, [x/1000 for x in operational_costs], marker = 'o', label="OC")
292 ax1.plot(depths, [x/1000 for x in operational_costs_cooling], marker='o', label="OCc")
293 ax1.plot(depths, [x/1000 for x in operational_costs_heating], marker='o', label="OCh")
294 ax1.set_xlabel(r'Depth (m)', fontsize=14)
295 ax1.set_ylabel(r'Costs ($k€$)', fontsize=14)
296 ax1.legend(loc='lower left', ncol=3)
297 ax1.tick_params(labelsize=14)
298 plt.show()
299
300
301if __name__ == "__main__": # pragma: no cover
302 active_passive_cooling()