Combination of active and passive cooling

  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()