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