diff --git a/CHANGELOG.md b/CHANGELOG.md
index 78ae29da..d29d9e71 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 ### Added
 
 - Add flow_rates as explicit argument in PressureDrop curve function.
+- Add function for cascading heat pumps.
+- Added _get_max_power functions to BuildingLoad classes.
 - Inlet and outlet water temperature (issue #271).
 - Added 'type' attribute to temperature profile plotting so also the inlet and outlet temperatures can be shown (issue
   #271).
diff --git a/GHEtool/Borefield.py b/GHEtool/Borefield.py
index 69f3fbbc..14012bc1 100644
--- a/GHEtool/Borefield.py
+++ b/GHEtool/Borefield.py
@@ -2012,19 +2012,17 @@ def get_rb(temperature, limit=None, power=None):
                                    hourly_load) / self.number_of_boreholes / H)
                 if not self.borehole.use_constant_Rb:
                     results._Tf_inlet, results._Tf_outlet = self.calculate_borefield_inlet_outlet_temperature(
-                        self.load.hourly_net_resulting_injection_power, results.peak_injection,
-                        simulation_period=self.load.simulation_period)
+                        hourly_load, results.peak_injection, simulation_period=self.load.simulation_period)
                     if sizing:
                         results._Tf_extraction_inlet, results._Tf_extraction_outlet = self.calculate_borefield_inlet_outlet_temperature(
-                            self.load.hourly_net_resulting_injection_power, results._Tf_extraction,
-                            simulation_period=self.load.simulation_period)
+                            hourly_load, results._Tf_extraction, simulation_period=self.load.simulation_period)
             return results
 
         def calculate_difference(results_old: Union[ResultsMonthly, ResultsHourly],
                                  result_new: Union[ResultsMonthly, ResultsHourly]) -> float:
             return max(
-                np.max(result_new.peak_injection - results_old.peak_injection),
-                np.max(result_new.peak_extraction - results_old.peak_extraction))
+                np.max(np.abs(result_new.peak_injection - results_old.peak_injection)),
+                np.max(np.abs(result_new.peak_extraction - results_old.peak_extraction)))
 
         if isinstance(self.load, _LoadDataBuilding) or \
                 isinstance(self.borehole.fluid_data, TemperatureDependentFluidData):
diff --git a/GHEtool/Examples/sizing_with_building_load_hourly.py b/GHEtool/Examples/sizing_with_building_load_hourly.py
index c3dab5f7..43640803 100644
--- a/GHEtool/Examples/sizing_with_building_load_hourly.py
+++ b/GHEtool/Examples/sizing_with_building_load_hourly.py
@@ -57,6 +57,7 @@ def L3_sizing() -> Tuple[float, float]:
     print(
         f'When sizing with an L3 method, the required borehole length is {length:.2f}m. '
         f'The SCOP is {borefield.load.SCOP_total:.2f}.')
+    borefield.load.SEER
     borefield.print_temperature_profile()
     return length, borefield.load.SCOP_heating
 
diff --git a/GHEtool/Examples/sizing_with_building_load_hourly_limit.py b/GHEtool/Examples/sizing_with_building_load_hourly_limit.py
new file mode 100644
index 00000000..b8188a06
--- /dev/null
+++ b/GHEtool/Examples/sizing_with_building_load_hourly_limit.py
@@ -0,0 +1,174 @@
+"""
+This document contains an example on sizing with an hourly building load and the difference between L3 and L4 sizing.
+It can be seen that sizing with an hourly method not only gives a completely different result, but also that the SCOP
+differs and is more accurate.
+"""
+
+from GHEtool import *
+from GHEtool.VariableClasses.Efficiency._Efficiency import combine_n_heat_pumps
+from typing import Tuple
+
+# initiate ground data
+data = GroundFluxTemperature(2.1, 10, flux=0.07)
+borefield_gt = gt.borefield.Borefield.rectangle_field(10, 12, 6, 6, 110, 1, 0.075)
+
+points_HP500 = np.array([
+    [-4.5, 84.7],
+    [-4.5, 66.3],
+    [-4.5, 38.6],
+    [-1.5, 93.0],
+    [-1.5, 72.9],
+    [-1.5, 42.4],
+    [3.5, 107.8],
+    [3.5, 84.6],
+    [3.5, 49.4],
+    [8.5, 123.5],
+    [8.5, 97.4],
+    [8.5, 56.9],
+    [11.5, 133.2],
+    [11.5, 105.4],
+    [11.5, 61.4],
+])
+eff_HP500 = np.array([
+    3.87, 4.12, 3.78,
+    4.12, 4.39, 4.04,
+    4.57, 4.89, 4.49,
+    5.10, 5.50, 5.08,
+    5.48, 5.92, 5.48,
+])
+points_HP400 = np.array([
+    [-4.5, 67.2],
+    [-4.5, 49.9],
+    [-4.5, 29.1],
+    [-1.5, 73.7],
+    [-1.5, 54.8],
+    [-1.5, 32.0],
+    [3.5, 85.2],
+    [3.5, 63.5],
+    [3.5, 37.2],
+    [8.5, 97.7],
+    [8.5, 73.1],
+    [8.5, 42.7],
+    [11.5, 105.7],
+    [11.5, 79.0],
+    [11.5, 46.0],
+])
+eff_HP400 = np.array([
+    3.91, 4.38, 3.99,
+    4.14, 4.64, 4.27,
+    4.58, 5.16, 4.77,
+    5.12, 5.80, 5.34,
+    5.51, 6.22, 5.75,
+])
+cascaded_system_points, cascaded_system_eff = combine_n_heat_pumps([points_HP500] * 4 + [points_HP400],
+                                                                   [eff_HP500] * 4 + [eff_HP400])
+
+cop_system = COP(cascaded_system_eff, cascaded_system_points, True)
+
+
+def L3_sizing() -> Tuple[float, float]:
+    """
+    Size the borefield with a monthly sizing method.
+
+    Returns
+    -------
+    length, SCOP
+    """
+    # initiate borefield
+    borefield = Borefield()
+
+    # set parameters
+    borefield.set_min_fluid_temperature(0)
+    borefield.set_max_fluid_temperature(25)
+
+    # set ground data in borefield
+    borefield.ground_data = data
+
+    # set Rb
+    borefield.Rb = 0.12
+
+    # set borefield
+    borefield.set_borefield(borefield_gt)
+
+    # load the hourly profile
+    load = HourlyBuildingLoad(efficiency_heating=cop_system)
+    load.load_hourly_profile(FOLDER.joinpath("Examples/hourly_profile.csv"), header=True, separator=";")
+    borefield.load = load
+
+    # size the borefield and plot the resulting temperature evolution
+    length = borefield.size(100, L3_sizing=True)
+    print(
+        f'When sizing with an L3 method (without limit), the required borehole length is {length:.2f}m. '
+        f'The SCOP is {borefield.load.SCOP_total:.2f}.')
+
+    # with limit
+    borefield.load._limit_to_max_heat_pump_power = True
+    length = borefield.size(100, L3_sizing=True)
+    borefield.print_temperature_profile()
+    return length
+
+
+def L4_sizing() -> Tuple[float, float]:
+    """
+    Size the borefield with an hourly sizing method.
+
+    Returns
+    -------
+    length, SCOP
+    """
+    # initiate borefield
+    borefield = Borefield()
+
+    # set parameters
+    borefield.set_min_fluid_temperature(0)
+    borefield.set_max_fluid_temperature(25)
+
+    # set ground data in borefield
+    borefield.ground_data = data
+
+    # set Rb
+    borefield.Rb = 0.12
+
+    # set borefield
+    borefield.set_borefield(borefield_gt)
+
+    # load the hourly profile
+    load = HourlyBuildingLoad(efficiency_heating=cop_system)
+    load.load_hourly_profile(FOLDER.joinpath("Examples/hourly_profile.csv"), header=True, separator=";")
+    borefield.load = load
+
+    # size the borefield and plot the resulting temperature evolution
+    length = borefield.size(100, L4_sizing=True)
+    print(
+        f'When sizing with an L4 method, the required borehole length is {length:.2f}m. '
+        f'The SCOP is {borefield.load.SCOP_total:.2f}.')
+    borefield.load._limit_to_max_heat_pump_power = True
+    length = borefield.size(100, L4_sizing=True)
+    print(
+        f'When sizing with an L4 method (with limit), the required borehole length is {length:.2f}m. '
+        f'The SCOP is {borefield.load.SCOP_total:.2f}.')
+    missing_power = borefield.load.hourly_heating_load_simulation_period - borefield.load.cop._get_max_power(
+        borefield.results.Tf)
+    missing_power = np.maximum(missing_power, 0)
+    print(f'Over the whole simulation period, a maximum of {max(missing_power):.2f} kW cannot be delivered by the '
+          f'cascaded heat pumps at 0°C. This accumulates to {np.sum(missing_power):.2f} kWh over the whole simulation period.')
+
+    borefield.set_min_fluid_temperature(3)
+    length = borefield.size(100, L4_sizing=True)
+    print(
+        f'When sizing with an L4 method (with limit at 3°C), the required borehole length is {length:.2f}m. '
+        f'The SCOP is {borefield.load.SCOP_total:.2f}.')
+    missing_power = borefield.load.hourly_heating_load_simulation_period - borefield.load.cop._get_max_power(
+        borefield.results.Tf)
+    missing_power = np.maximum(missing_power, 0)
+    print(f'Over the whole simulation period, a maximum of {max(missing_power):.2f} kW cannot be delivered by the '
+          f'cascaded heat pumps at 0°C. This accumulates to {np.sum(missing_power):.2f} kWh over the whole simulation period.')
+
+    borefield.print_temperature_profile(plot_hourly=True)
+    borefield.load.SEER
+    return length, borefield.load.SCOP_heating
+
+
+if __name__ == "__main__":
+    L3_sizing()
+    L4_sizing()
diff --git a/GHEtool/Examples/variable_flow_rate.py b/GHEtool/Examples/variable_flow_rate.py
index acb13b1c..f2c115d6 100644
--- a/GHEtool/Examples/variable_flow_rate.py
+++ b/GHEtool/Examples/variable_flow_rate.py
@@ -79,12 +79,30 @@ def sizing():
     plt.legend()
     plt.ylabel('Average fluid temperature [°C]')
     plt.xlabel('Time [years]')
+    plt.title('Borehole fluid temperature')
 
     plt.figure()
     plt.plot(time_array, rb_variable_flow, label='variable flow')
     plt.plot(time_array, rb_constant_flow, label='constant flow')
     plt.ylabel('Effective borehole thermal resistance [mK/W]')
     plt.xlabel('Time [years]')
+    plt.title('Effective borehole thermal resistance')
+    plt.legend()
+
+    plt.figure()
+    plt.plot(time_array[:8760], borefield_variable_flow.results.peak_injection[:8760], label='variable flow')
+    plt.plot(time_array[:8760], borefield_constant_flow.results.peak_injection[:8760], label='constant flow')
+    plt.legend()
+    plt.title('Borehole fluid temperature')
+    plt.ylabel('Average fluid temperature [°C]')
+    plt.xlabel('Time [years]')
+
+    plt.figure()
+    plt.plot(time_array[:8760], rb_variable_flow[:8760], label='variable flow')
+    plt.plot(time_array[:8760], rb_constant_flow[:8760], label='constant flow')
+    plt.ylabel('Effective borehole thermal resistance [mK/W]')
+    plt.title('Effective borehole thermal resistance')
+    plt.xlabel('Time [years]')
     plt.legend()
     plt.show()
 
diff --git a/GHEtool/Validation/enhanced_sizing_method.py b/GHEtool/Validation/enhanced_sizing_method.py
index d51f869e..d7c15220 100644
--- a/GHEtool/Validation/enhanced_sizing_method.py
+++ b/GHEtool/Validation/enhanced_sizing_method.py
@@ -43,18 +43,18 @@ def test_monthly_quadrant_4():
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 93.26359850472193)
-    assert np.isclose(borefield.results.min_temperature, 0.00014450766670659476)
-    assert np.isclose(borefield.results.max_temperature, 15.625961732265106)
+    assert np.isclose(borefield.H, 93.24146806456535)
+    assert np.isclose(borefield.results.min_temperature, 0.0001600069660581127)
+    assert np.isclose(borefield.results.max_temperature, 15.739507592388385)
     print(f'Simulation time with speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = False
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 93.26372520583197)
-    assert np.isclose(borefield.results.min_temperature, 0.00014378529190839373)
-    assert np.isclose(borefield.results.max_temperature, 15.625956914178497)
+    assert np.isclose(borefield.H, 93.243081505459)
+    assert np.isclose(borefield.results.min_temperature, 0.00015381852217721814)
+    assert np.isclose(borefield.results.max_temperature, 15.739453158826853)
     print(f'Simulation time without speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = True
@@ -63,9 +63,9 @@ def test_monthly_quadrant_4():
 
     borefield.size_L3()
     borefield.calculation_setup(approximate_req_depth=True)
-    assert np.isclose(borefield.H, 93.26359850472193)
-    assert np.isclose(borefield.results.min_temperature, 0.00014450766670659476)
-    assert np.isclose(borefield.results.max_temperature, 15.625961732265106)
+    assert np.isclose(borefield.H, 93.24146806456535)
+    assert np.isclose(borefield.results.min_temperature, 0.0001600069660581127)
+    assert np.isclose(borefield.results.max_temperature, 15.739507592388385)
     print(f'Simulation time with speed up and approximate_req_depth {time.time() - start}s')
 
     borefield.calculation_setup(use_explicit_multipole=True)
@@ -73,18 +73,18 @@ def test_monthly_quadrant_4():
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 93.23667934535844)
-    assert np.isclose(borefield.results.min_temperature, 0.00014464998743424218)
-    assert np.isclose(borefield.results.max_temperature, 15.617688373954934)
+    assert np.isclose(borefield.H, 93.21459868840302)
+    assert np.isclose(borefield.results.min_temperature, 0.00016020852361098292)
+    assert np.isclose(borefield.results.max_temperature, 15.728080432710481)
     print(f'Simulation time with speed up and explicit models {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = False
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 93.23667934535844)
-    assert np.isclose(borefield.results.min_temperature, 0.00014464998743424218)
-    assert np.isclose(borefield.results.max_temperature, 15.617688373954934)
+    assert np.isclose(borefield.H, 93.21528026470325)
+    assert np.isclose(borefield.results.min_temperature, 0.00015743055001760098)
+    assert np.isclose(borefield.results.max_temperature, 15.728056277153398)
     print(f'Simulation time without speed up and explicit models {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = True
@@ -93,9 +93,9 @@ def test_monthly_quadrant_4():
 
     borefield.size_L3()
     borefield.calculation_setup(approximate_req_depth=True)
-    assert np.isclose(borefield.H, 93.23667934535844)
-    assert np.isclose(borefield.results.min_temperature, 0.00014464998743424218)
-    assert np.isclose(borefield.results.max_temperature, 15.617688373954934)
+    assert np.isclose(borefield.H, 93.21459868840302)
+    assert np.isclose(borefield.results.min_temperature, 0.00016020852361098292)
+    assert np.isclose(borefield.results.max_temperature, 15.728080432710481)
     print(f'Simulation time with speed up, approximate_req_depth and explicit models {time.time() - start}s')
 
 
@@ -122,18 +122,18 @@ def test_monthly_quadrant_1():
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 83.09881599111853)
-    assert np.isclose(borefield.results.min_temperature, 5.296214968773752)
-    assert np.isclose(borefield.results.max_temperature, 16.00273196708503)
+    assert np.isclose(borefield.H, 85.79127756099108)
+    assert np.isclose(borefield.results.min_temperature, 5.451954106073291)
+    assert np.isclose(borefield.results.max_temperature, 16.009933134246406)
     print(f'Simulation time with speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = False
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 83.17430886882794)
-    assert np.isclose(borefield.results.min_temperature, 5.300760137618011)
-    assert np.isclose(borefield.results.max_temperature, 15.99936830071769)
+    assert np.isclose(borefield.H, 86.04140068480496)
+    assert np.isclose(borefield.results.min_temperature, 5.459660593826807)
+    assert np.isclose(borefield.results.max_temperature, 15.999214960639087)
     print(f'Simulation time without speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = True
@@ -142,9 +142,9 @@ def test_monthly_quadrant_1():
 
     borefield.size_L3()
     borefield.calculation_setup(approximate_req_depth=True)
-    assert np.isclose(borefield.H, 83.09881599111853)
-    assert np.isclose(borefield.results.min_temperature, 5.296214968773752)
-    assert np.isclose(borefield.results.max_temperature, 16.00273196708503)
+    assert np.isclose(borefield.H, 85.79127756099108)
+    assert np.isclose(borefield.results.min_temperature, 5.451954106073291)
+    assert np.isclose(borefield.results.max_temperature, 16.009933134246406)
     print(f'Simulation time with speed up and approximate_req_depth {time.time() - start}s')
 
 
@@ -171,9 +171,9 @@ def test_monthly_quadrant_1_more_data_points():
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 83.09881599111853)
-    assert np.isclose(borefield.results.min_temperature, 5.296214968773752)
-    assert np.isclose(borefield.results.max_temperature, 16.00273196708503)
+    assert np.isclose(borefield.H, 85.79127756099108)
+    assert np.isclose(borefield.results.min_temperature, 5.451954106073291)
+    assert np.isclose(borefield.results.max_temperature, 16.009933134246406)
     print(f'Simulation time with speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = False
@@ -181,9 +181,9 @@ def test_monthly_quadrant_1_more_data_points():
     borefield.borehole._nb_of_data_points = 500
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 83.09936353547558)
-    assert np.isclose(borefield.results.min_temperature, 5.29624824242701)
-    assert np.isclose(borefield.results.max_temperature, 15.99936830071769)
+    assert np.isclose(borefield.H, 85.9223631352438)
+    assert np.isclose(borefield.results.min_temperature, 5.459485430207949)
+    assert np.isclose(borefield.results.max_temperature, 15.999114386583042)
     print(f'Simulation time without speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = True
@@ -192,9 +192,9 @@ def test_monthly_quadrant_1_more_data_points():
 
     borefield.size_L3()
     borefield.calculation_setup(approximate_req_depth=True)
-    assert np.isclose(borefield.H, 83.09936353547558)
-    assert np.isclose(borefield.results.min_temperature, 5.29624824242701)
-    assert np.isclose(borefield.results.max_temperature, 15.99936830071769)
+    assert np.isclose(borefield.H, 85.79127756099108)
+    assert np.isclose(borefield.results.min_temperature, 5.451954394523563)
+    assert np.isclose(borefield.results.max_temperature, 16.004510365189816)
     print(f'Simulation time with speed up and approximate_req_depth {time.time() - start}s')
 
 
@@ -223,18 +223,18 @@ def test_monthly_quadrant_2():
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 116.00524762354382)
-    assert np.isclose(borefield.results.min_temperature, 8.407099988410305)
-    assert np.isclose(borefield.results.max_temperature, 17.001626884710568)
+    assert np.isclose(borefield.H, 117.42360117447696)
+    assert np.isclose(borefield.results.min_temperature, 8.448478389721295)
+    assert np.isclose(borefield.results.max_temperature, 17.000155792690094)
     print(f'Simulation time with speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = False
     start = time.time()
 
     borefield.size_L3()
-    assert np.isclose(borefield.H, 116.0394997074172)
-    assert np.isclose(borefield.results.min_temperature, 8.408106764334239)
-    assert np.isclose(borefield.results.max_temperature, 17.000621479904506)
+    assert np.isclose(borefield.H, 117.40961843257512)
+    assert np.isclose(borefield.results.min_temperature, 8.448073848783803)
+    assert np.isclose(borefield.results.max_temperature, 17.000585066710645)
     print(f'Simulation time without speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = True
@@ -243,9 +243,9 @@ def test_monthly_quadrant_2():
 
     borefield.size_L3()
     borefield.calculation_setup(approximate_req_depth=True)
-    assert np.isclose(borefield.H, 116.00524762354382)
-    assert np.isclose(borefield.results.min_temperature, 8.407099988410305)
-    assert np.isclose(borefield.results.max_temperature, 17.001626884710568)
+    assert np.isclose(borefield.H, 117.42360117447696)
+    assert np.isclose(borefield.results.min_temperature, 8.448478389721295)
+    assert np.isclose(borefield.results.max_temperature, 17.000155792690094)
     print(f'Simulation time with speed up and approximate_req_depth {time.time() - start}s')
 
 
@@ -269,18 +269,18 @@ def test_case_office():
     start = time.time()
 
     borefield.size_L4()
-    assert np.isclose(borefield.H, 141.96534390958553)
-    assert np.isclose(borefield.results.min_temperature, 8.643460326864346)
-    assert np.isclose(borefield.results.max_temperature, 17.00064802237182)
+    assert np.isclose(borefield.H, 141.9659861806117)
+    assert np.isclose(borefield.results.min_temperature, 8.643499473771335)
+    assert np.isclose(borefield.results.max_temperature, 17.000616753506712)
     print(f'Simulation time with speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = False
     start = time.time()
 
     borefield.size_L4()
-    assert np.isclose(borefield.H, 141.968830876313)
-    assert np.isclose(borefield.results.min_temperature, 8.643563349090133)
-    assert np.isclose(borefield.results.max_temperature, 17.00055762922439)
+    assert np.isclose(borefield.H, 141.965075411997)
+    assert np.isclose(borefield.results.min_temperature, 8.643452832019406)
+    assert np.isclose(borefield.results.max_temperature, 17.000654666548748)
     print(f'Simulation time without speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = True
@@ -289,9 +289,9 @@ def test_case_office():
 
     borefield.size_L4()
     borefield.calculation_setup(approximate_req_depth=True)
-    assert np.isclose(borefield.H, 141.9653439076926)
-    assert np.isclose(borefield.results.min_temperature, 8.643460326811507)
-    assert np.isclose(borefield.results.max_temperature, 17.000648022418627)
+    assert np.isclose(borefield.H, 141.9659861806117)
+    assert np.isclose(borefield.results.min_temperature, 8.643499473771335)
+    assert np.isclose(borefield.results.max_temperature, 17.000616753506712)
     print(f'Simulation time with speed up and approximate_req_depth {time.time() - start}s')
 
     borefield.calculation_setup(use_explicit_multipole=True)
@@ -299,18 +299,18 @@ def test_case_office():
     start = time.time()
 
     borefield.size_L4()
-    assert np.isclose(borefield.H, 141.955823246302)
-    assert np.isclose(borefield.results.min_temperature, 8.643335762145801)
-    assert np.isclose(borefield.results.max_temperature, 17.000646440707527)
+    assert np.isclose(borefield.H, 141.95646363193922)
+    assert np.isclose(borefield.results.min_temperature, 8.643374849202972)
+    assert np.isclose(borefield.results.max_temperature, 17.000615221248903)
     print(f'Simulation time with speed up and explicit models {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = False
     start = time.time()
 
     borefield.size_L4()
-    assert np.isclose(borefield.H, 141.96035547742298)
-    assert np.isclose(borefield.results.min_temperature, 8.643467717894024)
-    assert np.isclose(borefield.results.max_temperature, 17.000530338271854)
+    assert np.isclose(borefield.H, 141.9566103028618)
+    assert np.isclose(borefield.results.min_temperature, 8.643378819194513)
+    assert np.isclose(borefield.results.max_temperature, 17.000611682426978)
     print(f'Simulation time without speed up and explicit models {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = True
@@ -319,9 +319,9 @@ def test_case_office():
 
     borefield.size_L4()
     borefield.calculation_setup(approximate_req_depth=True)
-    assert np.isclose(borefield.H, 141.955823246302)
-    assert np.isclose(borefield.results.min_temperature, 8.643335762145801)
-    assert np.isclose(borefield.results.max_temperature, 17.000646440707527)
+    assert np.isclose(borefield.H, 141.95646363193922)
+    assert np.isclose(borefield.results.min_temperature, 8.643374849202972)
+    assert np.isclose(borefield.results.max_temperature, 17.000615221248903)
     print(f'Simulation time with speed up, approximate_req_depth and explicit models {time.time() - start}s')
 
 
@@ -346,18 +346,18 @@ def test_case_auditorium():
     start = time.time()
 
     borefield.size_L4()
-    assert np.isclose(borefield.H, 146.56168945870414)
-    assert np.isclose(borefield.results.min_temperature, 7.627998329436393)
-    assert np.isclose(borefield.results.max_temperature, 17.00060373723548)
+    assert np.isclose(borefield.H, 146.5617348169237)
+    assert np.isclose(borefield.results.min_temperature, 7.627999444626022)
+    assert np.isclose(borefield.results.max_temperature, 17.000602713187273)
     print(f'Simulation time with speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = False
     start = time.time()
 
     borefield.size_L4()
-    assert np.isclose(borefield.H, 146.56251979924573)
-    assert np.isclose(borefield.results.min_temperature, 7.628019277863065)
-    assert np.isclose(borefield.results.max_temperature, 17.000585005188466)
+    assert np.isclose(borefield.H, 146.56434784374707)
+    assert np.isclose(borefield.results.min_temperature, 7.62890725121998)
+    assert np.isclose(borefield.results.max_temperature, 17.00056637799696)
     print(f'Simulation time without speed up {time.time() - start}s')
 
     borefield.USE_SPEED_UP_IN_SIZING = True
@@ -366,9 +366,9 @@ def test_case_auditorium():
 
     borefield.size_L4()
     borefield.calculation_setup(approximate_req_depth=True)
-    assert np.isclose(borefield.H, 146.56168944805887)
-    assert np.isclose(borefield.results.min_temperature, 7.627998329174655)
-    assert np.isclose(borefield.results.max_temperature, 17.000603737475817)
+    assert np.isclose(borefield.H, 146.5617348169237)
+    assert np.isclose(borefield.results.min_temperature, 7.627999444626022)
+    assert np.isclose(borefield.results.max_temperature, 17.000602713187273)
     print(f'Simulation time with speed up and approximate_req_depth {time.time() - start}s')
 
 
diff --git a/GHEtool/VariableClasses/Efficiency/EERCombined.py b/GHEtool/VariableClasses/Efficiency/EERCombined.py
index 4f6cbeb5..86dbe9cf 100644
--- a/GHEtool/VariableClasses/Efficiency/EERCombined.py
+++ b/GHEtool/VariableClasses/Efficiency/EERCombined.py
@@ -186,6 +186,64 @@ def get_SEER(self,
 
         return np.sum(power) / np.sum(w_array)
 
+    def _get_max_power(self,
+                       primary_temperature: Union[float, np.ndarray],
+                       secondary_temperature: Union[float, np.ndarray] = None,
+                       month_indices: Union[float, np.ndarray] = None) -> Union[np.ndarray, float]:
+        """
+        This function calculates the maximum available power for a given primary and secondary temperature.
+
+        Parameters
+        ----------
+        primary_temperature : np.ndarray or float
+            Value(s) for the average primary temperature of the heat pump for the EER calculation.
+        secondary_temperature : np.ndarray or float
+            Value(s) for the average secondary temperature of the heat pump for the EER calculation.
+        month_indices : np.ndarray or float
+            Array with all the monthly indices, after correction for the start month. Should be the same length as the
+            other input parameters
+
+        Raises
+        ------
+        ValueError
+            When secondary_temperature is in the dataset, and it is not provided. Same for power.
+
+        Returns
+        -------
+        EER
+            np.ndarray
+        """
+
+        if isinstance(primary_temperature, (float, int)) and (
+                isinstance(month_indices, (float, int)) or month_indices is None):
+            # check temperature threshold
+            active_cooling_bool = False
+            if self.threshold_temperature is not None and primary_temperature > self.threshold_temperature:
+                active_cooling_bool = True
+
+            if not active_cooling_bool and self.months_active_cooling is not None:
+                if month_indices is None:
+                    raise ValueError('Please provide a month value, for otherwise the system cannot decide if it is '
+                                     'active or passive cooling.')
+                active_cooling_bool = month_indices in self.months_active_cooling
+
+            if active_cooling_bool:
+                return self.efficiency_active_cooling._get_max_power(primary_temperature, secondary_temperature)
+            return self.efficiency_passive_cooling._get_max_power(primary_temperature, secondary_temperature)
+
+        # now for monthly loads
+        active_cooling_power = self.efficiency_active_cooling._get_max_power(primary_temperature, secondary_temperature)
+        passive_cooling_power = self.efficiency_passive_cooling._get_max_power(primary_temperature,
+                                                                               secondary_temperature)
+
+        if month_indices is not None and isinstance(primary_temperature, (float, int)):
+            primary_temperature = np.full(month_indices.shape, primary_temperature)
+
+        active_cooling_bool = self.get_time_series_active_cooling(primary_temperature, month_indices)
+
+        # select correct data
+        return active_cooling_power * active_cooling_bool + passive_cooling_power * np.invert(active_cooling_bool)
+
     def __eq__(self, other) -> bool:
         if not isinstance(other, self.__class__):
             return False
diff --git a/GHEtool/VariableClasses/Efficiency/SEER.py b/GHEtool/VariableClasses/Efficiency/SEER.py
index 3d14b5cc..649a6779 100644
--- a/GHEtool/VariableClasses/Efficiency/SEER.py
+++ b/GHEtool/VariableClasses/Efficiency/SEER.py
@@ -76,5 +76,8 @@ def get_EER(self, *args, **kwargs) -> float:
         """
         return self.SEER
 
+    def _get_max_power(self, *args, **kwargs) -> float:
+        return 1e16
+
     def __export__(self):
         return {'SEER [-]': self.SEER}
diff --git a/GHEtool/VariableClasses/Efficiency/_Efficiency.py b/GHEtool/VariableClasses/Efficiency/_Efficiency.py
index 7a0cb143..70691e9d 100644
--- a/GHEtool/VariableClasses/Efficiency/_Efficiency.py
+++ b/GHEtool/VariableClasses/Efficiency/_Efficiency.py
@@ -1,8 +1,12 @@
+import itertools
+
+import matplotlib.pyplot as plt
 import numpy as np
 
-from scipy.interpolate import interpn, interp1d
-from typing import Union
+from collections import defaultdict
 from GHEtool.VariableClasses.BaseClass import BaseClass
+from scipy.interpolate import interpn
+from typing import Union
 
 
 class _EfficiencyBase(BaseClass):
@@ -161,6 +165,23 @@ def find_value(x, y, z=None):
                         self._data[i, j, k] = find_value(self._range_primary[i],
                                                          self._range_secondary[j],
                                                          self._range_part_load[k])
+            # get max powers per temperature
+            x = coordinates[:, 0]
+            y = coordinates[:, 1]
+            z = coordinates[:, 2]
+
+            ix = np.searchsorted(self._range_primary, x)
+            iy = np.searchsorted(self._range_secondary, y)
+
+            flat_idx = ix * len(self._range_secondary) + iy
+
+            max_z_flat = np.full(len(self._range_primary) * len(self._range_secondary), -np.inf)
+            np.maximum.at(max_z_flat, flat_idx, z)
+
+            max_z = max_z_flat.reshape(len(self._range_primary), len(self._range_secondary))
+            # convert index to part load value
+            self._max_part_load = max_z
+
         elif dimensions == 2:
             self._data = np.empty(
                 (len(self._range_primary), max(len(self._range_secondary), len(self._range_part_load))))
@@ -174,6 +195,17 @@ def find_value(x, y, z=None):
                     for j in range(self._data.shape[1]):
                         self._data[i, j] = find_value(self._range_primary[i],
                                                       self._range_part_load[j])
+
+                # get max powers per temperature
+                x = coordinates[:, 0]
+                y = coordinates[:, 1]
+
+                idx = np.searchsorted(self._range_primary, x)
+
+                max_y = np.full(len(self._range_primary), -np.inf)
+                np.maximum.at(max_y, idx, y)
+
+                self._max_part_load = max_y
         else:
             p = self._range_primary.argsort()
             self._data = data[p]
@@ -256,3 +288,448 @@ def _get_efficiency(self,
         interp = interpn(self._points, self._data, xi, bounds_error=False, fill_value=np.nan)
         if not np.isnan(interp).any():
             return interp
+
+    def _get_max_power(self,
+                       primary_temperature: Union[float, np.ndarray],
+                       secondary_temperature: Union[float, np.ndarray] = None,**kwargs) -> np.ndarray:
+        """
+        This function returns the maximum available power for a certain primary and secondary temperature.
+
+        Parameters
+        ----------
+        primary_temperature : np.ndarray or float
+            Value(s) for the average primary temperature of the heat pump for the efficiency calculation.
+        secondary_temperature : np.ndarray or float
+            Value(s) for the average secondary temperature of the heat pump for the efficiency calculation.
+
+        Raises
+        ------
+        ValueError
+            When secondary_temperature is in the dataset, and it is not provided. Same for power.
+
+        Returns
+        -------
+        Efficiency
+            np.ndarray
+        """
+
+        if not self._has_part_load:
+            return 1e16
+
+            # reuse your existing clipping and array logic
+        _max_length = np.max([
+            len(i) if i is not None and not isinstance(i, (float, int)) else 1
+            for i in (primary_temperature, secondary_temperature)
+        ])
+
+        Tp = np.array(
+            np.full(_max_length, primary_temperature)
+            if isinstance(primary_temperature, (float, int))
+            else primary_temperature
+        )
+
+        Ts = None
+        if self._has_secondary:
+            if secondary_temperature is None:
+                raise ValueError("Secondary temperature is required.")
+            Ts = np.array(
+                np.full(_max_length, secondary_temperature)
+                if isinstance(secondary_temperature, (float, int))
+                else secondary_temperature
+            )
+
+        # clip
+        Tp = np.clip(Tp, np.min(self._range_primary), np.max(self._range_primary))
+        if self._has_secondary:
+            Ts = np.clip(Ts, np.min(self._range_secondary), np.max(self._range_secondary))
+
+        # interpolate directly on precomputed surface
+        if self._has_secondary:
+            xi = list(zip(Tp, Ts))
+        else:
+            xi = Tp
+
+        return interpn(self._points[:1 + self._has_secondary], self._max_part_load, xi, bounds_error=False,
+                       fill_value=np.nan)
+
+
+def plot_heat_pump_envelope(points, eff, ax=None, label_prefix="T"):
+    """
+    Plot heat pump efficiency as a function of power for each primary temperature.
+
+    Parameters
+    ----------
+    points : ndarray of shape (N, 2)
+        Array of (primary_temperature, available_power) pairs.
+    eff : ndarray of shape (N,)
+        Efficiencies corresponding to `points`.
+    ax : matplotlib.axes.Axes, optional
+        Axis to plot on. If None, a new figure and axis are created.
+    label_prefix : str, optional
+        Prefix for legend labels. Default is "T".
+
+    Returns
+    -------
+    ax : matplotlib.axes.Axes
+        Axis containing the plot.
+    """
+
+    if ax is None:
+        fig, ax = plt.subplots()
+
+    # group by temperature
+    grouped = defaultdict(lambda: {"power": [], "eff": []})
+
+    for (T, P), e in zip(points, eff):
+        grouped[T]["power"].append(P)
+        grouped[T]["eff"].append(e)
+
+    # plot each temperature
+    for T in sorted(grouped.keys()):
+        p = np.asarray(grouped[T]["power"])
+        e = np.asarray(grouped[T]["eff"])
+
+        # sort and remove duplicate power values
+        idx = np.argsort(p)
+        p = p[idx]
+        e = e[idx]
+
+        p_unique, idx_unique = np.unique(p, return_index=True)
+        e_unique = e[idx_unique]
+
+        ax.plot(p_unique, e_unique, marker="o", label=f"{label_prefix} = {T}")
+
+    ax.set_xlabel("Power")
+    ax.set_ylabel("Efficiency")
+    ax.grid(True)
+    ax.legend()
+
+    return ax
+
+
+def combine_n_heat_pumps(points_list, eff_list):
+    """
+    Combine the operating envelopes of multiple modulating heat pumps into a
+    single equivalent operating envelope using strict cascade staging.
+
+    At each primary temperature, the heat pumps are ordered by increasing
+    minimum available power and combined according to a deterministic
+    staging strategy:
+
+    Operating logic (per temperature)
+    ---------------------------------
+    1. Single-machine operation
+       Below the sum of the minimum powers of the two smallest heat pumps,
+       only the smallest heat pump may operate.
+
+    2. Single-machine overlap (HP1 vs HP2 only)
+       In the same low-power region, both of the two smallest heat pumps
+       may operate individually. At equal part load, the heat pump with
+       the highest efficiency is selected.
+
+       No overlap regions are allowed beyond this first staging level.
+
+    3. Cascade operation with strict staging
+       Once the combined minimum power of k heat pumps is reached,
+       exactly k heat pumps operate simultaneously.
+
+       For each cascade stage:
+       - All active heat pumps operate at the same part load ratio.
+       - The combined efficiency is computed as a power-weighted average.
+       - Operation with fewer heat pumps is no longer allowed once a
+         higher cascade stage is available.
+
+       This enforces a monotonic staging sequence:
+       1 → 2 → 3 → … → n heat pumps.
+
+    Interpolation is permitted within the operating envelope of each heat
+    pump. Extrapolation outside the envelope is not allowed.
+
+    Parameters
+    ----------
+    points_list : list of ndarray
+        List of arrays, one per heat pump. Each array has shape (Ni, 2) and
+        contains (primary_temperature, available_power) pairs.
+    eff_list : list of ndarray
+        List of efficiency arrays corresponding to `points_list`. Each array
+        has shape (Ni,).
+
+    Returns
+    -------
+    combined_points : ndarray of shape (K, 2)
+        Combined array of (primary_temperature, available_power) pairs
+        representing the equivalent operating envelope.
+    combined_eff : ndarray of shape (K,)
+        Efficiencies corresponding to `combined_points`.
+
+    Notes
+    -----
+    - Multiple power levels per primary temperature are supported.
+    - Each primary temperature is processed independently.
+    - Once a cascade with k heat pumps is possible, operation with fewer
+      heat pumps is strictly disallowed.
+    - The output envelope represents physically allowed operating states
+      under a strict cascade control philosophy.
+    - The output format matches the input format.
+    """
+
+    def group_by_temperature(points, eff):
+        """
+        Group power and efficiency data by primary temperature.
+
+        Parameters
+        ----------
+        points : ndarray of shape (N, 2)
+            Array of (temperature, power) pairs.
+        eff : ndarray of shape (N,)
+            Efficiencies corresponding to `points`.
+
+        Returns
+        -------
+        grouped : dict
+            Dictionary keyed by temperature with values containing sorted
+            power and efficiency arrays.
+        """
+        grouped = defaultdict(lambda: {"power": [], "eff": []})
+
+        for (T, P), e in zip(points, eff):
+            grouped[T]["power"].append(P)
+            grouped[T]["eff"].append(e)
+
+        for T in grouped:
+            p = np.asarray(grouped[T]["power"])
+            e = np.asarray(grouped[T]["eff"])
+            idx = np.argsort(p)
+            grouped[T]["power"] = p[idx]
+            grouped[T]["eff"] = e[idx]
+
+        return grouped
+
+    def interp_eff(P, p_arr, e_arr):
+        """
+        Interpolate efficiency at a given power level.
+        """
+        return np.interp(P, p_arr, e_arr)
+
+    def combine_at_temperature_n(
+            hps,
+            n_pl_single=25,
+            n_pl_cascade=40
+    ):
+        """
+    Combine multiple heat pumps at a fixed primary temperature using strict
+    cascade staging and linspace-based part load discretization.
+
+    Heat pumps are ordered by increasing minimum available power and combined
+    according to the following rules:
+
+    - Below the first cascade threshold, only single-machine operation is
+      allowed.
+    - A single overlap zone exists only between the two smallest heat pumps
+      and only below the first cascade threshold.
+    - Above each cascade threshold, exactly k heat pumps operate
+      simultaneously.
+    - Once a higher cascade stage is available, operation with fewer heat
+      pumps is not permitted.
+
+    Part load behavior is evaluated on uniform linspace grids to produce
+    smooth operating envelopes.
+
+    Parameters
+    ----------
+    hps : list of dict
+        List of heat pump operating envelopes available at this temperature.
+        Each dict contains:
+        - "power" : ndarray
+            Sorted array of available powers.
+        - "eff" : ndarray
+            Efficiencies corresponding to "power".
+    n_pl_single : int, optional
+        Number of part load points used for the single-machine overlap region
+        between the two smallest heat pumps.
+    n_pl_cascade : int, optional
+        Number of part load points used for each cascade stage.
+
+    Returns
+    -------
+    P_comb : ndarray
+        Combined available powers at this temperature.
+    E_comb : ndarray
+        Corresponding combined efficiencies.
+
+    Notes
+    -----
+    - All active heat pumps in a cascade stage operate at the same part load.
+    - Combined efficiencies are computed as power-weighted averages.
+    - Strict staging ensures that exactly one operating mode is valid for
+      each power level.
+    """
+        # sort heat pumps by minimum power
+        hps = sorted(hps, key=lambda hp: hp["power"][0])
+        n = len(hps)
+
+        p_min = [hp["power"][0] for hp in hps]
+        p_max = [hp["power"][-1] for hp in hps]
+
+        P_comb = []
+        E_comb = []
+
+        # -------------------------
+        # zone 1: single HP1 only
+        # -------------------------
+        for P in hps[0]["power"]:
+            if n == 1 or P < p_min[1]:
+                P_comb.append(P)
+                E_comb.append(interp_eff(P, hps[0]["power"], hps[0]["eff"]))
+
+        # -------------------------
+        # overlap zone: HP1 vs HP2 ONLY
+        # -------------------------
+        if n >= 2:
+            hp1 = hps[0]
+            hp2 = hps[1]
+
+            P_overlap_max = p_min[0] + p_min[1]
+
+            pl_grid = np.linspace(0.0, 1.0, n_pl_single)
+
+            for pl in pl_grid:
+                P1 = p_min[0] + pl * (p_max[0] - p_min[0])
+                P2 = p_min[1] + pl * (p_max[1] - p_min[1])
+
+                candidates = []
+
+                if p_min[0] <= P1 <= p_max[0]:
+                    candidates.append((P1, interp_eff(P1, hp1["power"], hp1["eff"])))
+
+                if p_min[1] <= P2 <= p_max[1]:
+                    candidates.append((P2, interp_eff(P2, hp2["power"], hp2["eff"])))
+
+                if not candidates:  # pragma: no cover
+                    continue
+
+                P_best, E_best = max(candidates, key=lambda x: x[1])
+
+                if P_best < P_overlap_max:
+                    P_comb.append(P_best)
+                    E_comb.append(E_best)
+
+        # -------------------------
+        # cascade zones: exactly k+1 machines
+        # -------------------------
+        for k in range(1, n):
+            active = hps[:k + 1]
+
+            P_min_stage = sum(p_min[:k + 1])
+            P_max_stage = (sum(p_min[:k + 2]) if k + 1 < n else np.inf)
+
+            pl_grid = np.linspace(0.0, 1.0, n_pl_cascade)
+
+            for pl in pl_grid:
+                powers = []
+                effs = []
+
+                for i, hp in enumerate(active):
+                    Pi = p_min[i] + pl * (p_max[i] - p_min[i])
+                    if Pi < p_min[i] or Pi > p_max[i]:
+                        break  # pragma: no cover
+                    powers.append(Pi)
+                    effs.append(
+                        interp_eff(Pi, hp["power"], hp["eff"])
+                    )
+                else:
+                    P_tot = sum(powers)
+
+                    # STRICT staging window
+                    if not (P_min_stage <= P_tot < P_max_stage):
+                        continue
+
+                    E_tot = np.dot(powers, effs) / P_tot
+
+                    P_comb.append(P_tot)
+                    E_comb.append(E_tot)
+
+        # -------------------------
+        # cleanup
+        # -------------------------
+        P_comb = np.asarray(P_comb)
+        E_comb = np.asarray(E_comb)
+
+        mask = np.isfinite(E_comb)
+        idx = np.argsort(P_comb[mask])
+
+        return P_comb[mask][idx], E_comb[mask][idx]
+
+    # group each heat pump by temperature
+    hp_groups = [group_by_temperature(points, eff) for points, eff in zip(points_list, eff_list)]
+
+    combined_points = []
+    combined_eff = []
+
+    # all temperatures across all heat pumps
+    all_T = sorted(set().union(*[hp.keys() for hp in hp_groups]))
+
+    for T in all_T:
+        # collect all heat pumps available at this temperature
+        hps_at_T = [hp[T] for hp in hp_groups if T in hp]
+
+        if not hps_at_T:  # pragma: no cover
+            continue
+
+        if len(hps_at_T) == 1:
+            # only one machine available
+            P = hps_at_T[0]["power"]
+            E = hps_at_T[0]["eff"]
+        else:
+            # multiple machines → combine
+            P, E = combine_at_temperature_n(hps_at_T)
+
+        for p, e in zip(P, E):
+            combined_points.append((T, p))
+            combined_eff.append(e)
+
+    return np.asarray(combined_points), np.asarray(combined_eff)
+
+
+def _find_optimal_heat_pump_configuration(heat_pumps: list[_Efficiency], power: float, prim_temp: float,
+                                          sec_temp: float = None) -> list:
+    """
+    This function finds the required combination of different heat pumps to be able to deliver the required power
+    at a certain primary (and secondary) temperature.
+
+    Parameters
+    ----------
+    heat_pumps : list
+        List of efficiency objects
+    power : float
+        Required power [kW]
+    prim_temp : float
+        Primary temperature at which the power is required [°C]
+    sec_temp : float
+        Secondary temperature at which the power is required [°C]
+
+    Returns
+    -------
+    list
+        List with the number of heat pumps required for a certain power.
+    """
+
+    best = None
+    n = len(heat_pumps)
+
+    heat_pump_max_powers = [heat_pump._get_max_power(primary_temperature=prim_temp, secondary_temperature=sec_temp)[0]
+                            for
+                            heat_pump in heat_pumps]
+
+    for mask in itertools.product([0, 1], repeat=n):
+        mask = np.array(mask)
+        total = heat_pump_max_powers @ mask
+        if total >= power:
+            units = mask.sum()
+            overshoot = total - power
+            candidate = (units, overshoot, -total)
+            if best is None or candidate < best:
+                best = candidate
+                best_mask = mask
+
+    return best_mask
diff --git a/GHEtool/VariableClasses/LoadData/Baseclasses/_HourlyDataBuilding.py b/GHEtool/VariableClasses/LoadData/Baseclasses/_HourlyDataBuilding.py
index 2c8f83cb..6775051c 100644
--- a/GHEtool/VariableClasses/LoadData/Baseclasses/_HourlyDataBuilding.py
+++ b/GHEtool/VariableClasses/LoadData/Baseclasses/_HourlyDataBuilding.py
@@ -3,6 +3,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
+from pyparsing.tools.cvt_pyparsing_pep8_names import pre_pep8_method_name
 
 from ._HourlyData import _HourlyData
 from ._LoadDataBuilding import _LoadDataBuilding
@@ -179,9 +180,9 @@ def _get_hourly_cop(self, power: np.ndarray = None) -> Union[float, np.ndarray]:
         if isinstance(self.cop, SCOP) and isinstance(self.eer, SEER) and isinstance(self.cop_dhw, SCOP):
             return self.cop.get_COP(0, power=np.nan_to_num(power))
         if isinstance(self.results, ResultsMonthly):
-            raise TypeError('You cannot get an hourly EER values based on monthly temperature results.')
+            raise TypeError('You cannot get hourly COP values based on monthly temperature results.')
         if isinstance(self.results, tuple):
-            temperature = self.results[1]
+            temperature = self.results[0]
         else:
             temperature = self.results.Tf
 
@@ -209,9 +210,9 @@ def _get_hourly_cop_dhw(self, power: np.ndarray = None) -> Union[float, np.ndarr
         if isinstance(self.cop, SCOP) and isinstance(self.eer, SEER) and isinstance(self.cop_dhw, SCOP):
             return self.cop_dhw.get_COP(0, power=np.nan_to_num(power))
         if isinstance(self.results, ResultsMonthly):
-            raise TypeError('You cannot get an hourly EER values based on monthly temperature results.')
+            raise TypeError('You cannot get an hourly COP values based on monthly temperature results.')
         if isinstance(self.results, tuple):
-            temperature = self.results[1]
+            temperature = self.results[0]
         else:
             temperature = self.results.Tf
 
@@ -247,6 +248,88 @@ def _get_hourly_eer(self, power: np.ndarray = None) -> Union[float, np.ndarray]:
 
         return self.eer.get_EER(temperature, power=np.nan_to_num(power), month_indices=self.month_indices)
 
+    @property
+    def _get_max_power_heating(self) -> Union[float, np.ndarray]:
+        """
+        This function returns the maximum power available in heating, taking into account the maximum power of the
+        heat pump (when available). When the attribute '_limit_to_max_heat_pump_power' is False, the input array
+        is simply returned.
+
+        Returns
+        -------
+        maximum heating power : float | np.ndarray
+            Array with the maximum heating power
+
+        Raises
+        ------
+        TypeError
+            When the results are provided on a monthly basis
+        """
+        if isinstance(self.cop, SCOP) or not self._limit_to_max_heat_pump_power:
+            return self.hourly_heating_load_simulation_period
+        if isinstance(self.results, ResultsMonthly):
+            raise TypeError('You cannot get hourly COP values based on monthly temperature results.')
+        if isinstance(self.results, tuple):
+            temperature = self.results[0]
+        else:
+            temperature = self.results.Tf
+        return np.minimum(self.hourly_heating_load_simulation_period, self.cop._get_max_power(temperature))
+
+    @property
+    def _get_max_power_dhw(self) -> Union[float, np.ndarray]:
+        """
+        This function returns the maximum power available in dhw, taking into account the maximum power of the
+        heat pump (when available). When the attribute '_limit_to_max_heat_pump_power' is False, the input array
+        is simply returned.
+
+        Returns
+        -------
+        maximum dhw power : float | np.ndarray
+            Array with the maximum dhw power
+
+        Raises
+        ------
+        TypeError
+            When the results are provided on a monthly basis
+        """
+        if isinstance(self.cop_dhw, SCOP) or not self._limit_to_max_heat_pump_power:
+            return self.hourly_dhw_load_simulation_period
+        if isinstance(self.results, ResultsMonthly):
+            raise TypeError('You cannot get an hourly COP values based on monthly temperature results.')
+        if isinstance(self.results, tuple):
+            temperature = self.results[0]
+        else:
+            temperature = self.results.Tf
+        return np.minimum(self.hourly_dhw_load_simulation_period, self.cop_dhw._get_max_power(temperature))
+
+    @property
+    def _get_max_power_cooling(self) -> Union[float, np.ndarray]:
+        """
+        This function returns the maximum power available in cooling, taking into account the maximum power of the
+        heat pump (when available). When the attribute '_limit_to_max_heat_pump_power' is False, the input array
+        is simply returned.
+
+        Returns
+        -------
+        maximum cooling power : float | np.ndarray
+            Array with the maximum cooling power
+
+        Raises
+        ------
+        TypeError
+            When the results are provided on a monthly basis
+        """
+        if isinstance(self.eer, SEER) or not self._limit_to_max_heat_pump_power:
+            return self.hourly_cooling_load_simulation_period
+        if isinstance(self.results, ResultsMonthly):
+            raise TypeError('You cannot get an hourly EER values based on monthly temperature results.')
+        if isinstance(self.results, tuple):
+            temperature = self.results[1]
+        else:
+            temperature = self.results.Tf
+        return np.minimum(self.hourly_cooling_load_simulation_period,
+                          self.eer._get_max_power(temperature, month_indices=self.month_indices))
+
     @property
     def hourly_injection_load_simulation_period(self) -> np.ndarray:
         """
@@ -258,6 +341,10 @@ def hourly_injection_load_simulation_period(self) -> np.ndarray:
             Hourly injection values [kWh/h] for the whole simulation period
         """
         part_load = self.hourly_cooling_load_simulation_period
+        if self._limit_to_max_heat_pump_power:
+            return np.multiply(
+                self._get_max_power_cooling,
+                self.conversion_factor_secondary_to_primary_cooling(self._get_hourly_eer(part_load)))
         return np.multiply(
             self.hourly_cooling_load_simulation_period,
             self.conversion_factor_secondary_to_primary_cooling(self._get_hourly_eer(part_load)))
@@ -285,6 +372,10 @@ def _hourly_extraction_load_heating_simulation_period(self) -> np.ndarray:
             Hourly extraction values [kWh/h] for the whole simulation period
         """
         part_load = self.hourly_heating_load_simulation_period
+        if self._limit_to_max_heat_pump_power:
+            return np.multiply(
+                self._get_max_power_heating,
+                self.conversion_factor_secondary_to_primary_heating(self._get_hourly_cop(part_load)))
         return np.multiply(
             self.hourly_heating_load_simulation_period,
             self.conversion_factor_secondary_to_primary_heating(self._get_hourly_cop(part_load)))
@@ -300,6 +391,10 @@ def _hourly_extraction_load_dhw_simulation_period(self) -> np.ndarray:
             Hourly extraction values [kWh/h] for the whole simulation period
         """
         part_load_dhw = self.hourly_dhw_load_simulation_period
+        if self._limit_to_max_heat_pump_power:
+            return np.multiply(self._get_max_power_dhw,
+                               self.conversion_factor_secondary_to_primary_heating(
+                                   self._get_hourly_cop_dhw(part_load_dhw)))
         return np.multiply(
             self.hourly_dhw_load_simulation_period,
             self.conversion_factor_secondary_to_primary_heating(self._get_hourly_cop_dhw(part_load_dhw)))
@@ -328,6 +423,45 @@ def monthly_baseload_cooling_simulation_period(self) -> np.ndarray:
         """
         return self.resample_to_monthly(self.hourly_cooling_load_simulation_period)[1]
 
+    @property
+    def _monthly_baseload_heating_simulation_period_with_limit(self) -> np.ndarray:
+        """
+        This function returns the monthly heating baseload in kWh/month actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        baseload heating : np.ndarray
+            Baseload heating for the whole simulation period
+        """
+        return self.resample_to_monthly(self._get_max_power_heating)[1]
+
+    @property
+    def _monthly_baseload_cooling_simulation_period_with_limit(self) -> np.ndarray:
+        """
+        This function returns the monthly cooling baseload in kWh/month actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        baseload cooling : np.ndarray
+            Baseload cooling for the whole simulation period
+        """
+        return self.resample_to_monthly(self._get_max_power_cooling)[1]
+
+    @property
+    def _monthly_baseload_dhw_simulation_period_with_limit(self) -> np.ndarray:
+        """
+        This function returns the monthly domestic hot water baseload in kWh/month actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        baseload domestic hot water : np.ndarray
+            Baseload domestic hot water for the whole simulation period
+        """
+        return self.resample_to_monthly(self._get_max_power_dhw)[1]
+
     @property
     def monthly_peak_heating_simulation_period(self) -> np.ndarray:
         """
diff --git a/GHEtool/VariableClasses/LoadData/Baseclasses/_LoadDataBuilding.py b/GHEtool/VariableClasses/LoadData/Baseclasses/_LoadDataBuilding.py
index 7de5712e..a3bb1750 100644
--- a/GHEtool/VariableClasses/LoadData/Baseclasses/_LoadDataBuilding.py
+++ b/GHEtool/VariableClasses/LoadData/Baseclasses/_LoadDataBuilding.py
@@ -47,8 +47,12 @@ def __init__(self,
         self._cop_dhw = None
         self._results = None
         self._results_fixed = (0, 17)
-        self.exclude_DHW_from_peak: bool = False  # by default, the DHW increase the peak load. Set to false,
+
+        # False by default, the DHW increase the peak load. Set to false,
         # if you only want the heating load to determine the peak in extraction
+        self.exclude_DHW_from_peak: bool = False
+
+        self._limit_to_max_heat_pump_power: bool = False
 
         # set variables
         self.cop = efficiency_heating
@@ -110,6 +114,45 @@ def monthly_baseload_dhw_simulation_period(self) -> np.ndarray:
             Baseload domestic hot water for the whole simulation period
         """
 
+    @property
+    def _monthly_baseload_heating_simulation_period_with_limit(self) -> np.ndarray:
+        """
+        This function returns the monthly heating baseload in kWh/month actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        baseload heating : np.ndarray
+            Baseload heating for the whole simulation period
+        """
+        return self.monthly_baseload_heating_simulation_period
+
+    @property
+    def _monthly_baseload_cooling_simulation_period_with_limit(self) -> np.ndarray:
+        """
+        This function returns the monthly cooling baseload in kWh/month actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        baseload cooling : np.ndarray
+            Baseload cooling for the whole simulation period
+        """
+        return self.monthly_baseload_cooling_simulation_period
+
+    @property
+    def _monthly_baseload_dhw_simulation_period_with_limit(self) -> np.ndarray:
+        """
+        This function returns the monthly domestic hot water baseload in kWh/month actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        baseload domestic hot water : np.ndarray
+            Baseload domestic hot water for the whole simulation period
+        """
+        return self.monthly_baseload_dhw_simulation_period
+
     @abc.abstractmethod
     def set_results(self, results: Union[ResultsMonthly, ResultsHourly]) -> None:
         """
@@ -398,7 +441,7 @@ def _get_monthly_cop(self, peak: bool, power: np.ndarray = None) -> Union[float,
         elif peak:
             temperature = self.results.peak_extraction
         else:
-            temperature = self.results.monthly_extraction
+            temperature = self.results.baseload_temperature
 
         return self.cop.get_COP(temperature, power=np.nan_to_num(power))
 
@@ -424,7 +467,7 @@ def _get_monthly_cop_dhw(self, peak: bool, power: np.ndarray = None) -> Union[fl
         elif peak:
             temperature = self.results.peak_extraction
         else:
-            temperature = self.results.monthly_extraction
+            temperature = self.results.baseload_temperature
 
         return self.cop_dhw.get_COP(temperature, power=np.nan_to_num(power))
 
@@ -454,6 +497,64 @@ def _get_monthly_eer(self, peak: bool, power: np.ndarray = None) -> Union[float,
 
         return self.eer.get_EER(temperature, power=np.nan_to_num(power), month_indices=self.month_indices)
 
+    def _get_max_power_heating_monthly(self) -> Union[float, np.ndarray]:
+        """
+        This function returns the maximum power available in heating, taking into account the maximum power of the
+        heat pump (when available). When the attribute '_limit_to_max_heat_pump_power' is False, the input array
+        is simply returned.
+
+        Returns
+        -------
+        maximum heating power : float | np.ndarray
+            Array with the maximum heating power
+        """
+        if isinstance(self.cop, SCOP):
+            return self.monthly_peak_heating_simulation_period
+        if isinstance(self.results, tuple):
+            temperature = self.results[0]
+        else:
+            temperature = self.results.peak_extraction
+        return np.minimum(self.monthly_peak_heating_simulation_period, self.cop._get_max_power(temperature))
+
+    def _get_max_power_dhw_monthly(self) -> Union[float, np.ndarray]:
+        """
+        This function returns the maximum power available in dhw, taking into account the maximum power of the
+        heat pump (when available). When the attribute '_limit_to_max_heat_pump_power' is False, the input array
+        is simply returned.
+
+        Returns
+        -------
+        maximum dhw power : float | np.ndarray
+            Array with the maximum dhw power
+        """
+        if isinstance(self.cop_dhw, SCOP):
+            return self.monthly_peak_dhw_simulation_period
+        if isinstance(self.results, tuple):
+            temperature = self.results[0]
+        else:
+            temperature = self.results.peak_extraction
+        return np.minimum(self.monthly_peak_dhw_simulation_period, self.cop_dhw._get_max_power(temperature))
+
+    def _get_max_power_cooling_monthly(self) -> Union[float, np.ndarray]:
+        """
+        This function returns the maximum power available in cooling, taking into account the maximum power of the
+        heat pump (when available). When the attribute '_limit_to_max_heat_pump_power' is False, the input array
+        is simply returned.
+
+        Returns
+        -------
+        maximum cooling power : float | np.ndarray
+            Array with the maximum cooling power
+        """
+        if isinstance(self.eer, SEER):
+            return self.monthly_peak_cooling_simulation_period
+        if isinstance(self.results, tuple):
+            temperature = self.results[1]
+        else:
+            temperature = self.results.peak_injection
+        return np.minimum(self.monthly_peak_cooling_simulation_period,
+                          self.eer._get_max_power(temperature, month_indices=self.month_indices))
+
     @staticmethod
     def conversion_factor_secondary_to_primary_heating(cop_value: Union[int, float, np.ndarray]) -> Union[
         float, np.ndarray]:
@@ -560,6 +661,10 @@ def monthly_peak_injection_simulation_period(self) -> np.ndarray:
             Peak injection for the whole simulation period
         """
         part_load = self.monthly_peak_cooling_simulation_period
+        if self._limit_to_max_heat_pump_power:
+            return np.multiply(
+                self._get_max_power_cooling_monthly(),
+                self.conversion_factor_secondary_to_primary_cooling(self._get_monthly_eer(True, part_load)))
         return np.multiply(
             self.monthly_peak_cooling_simulation_period,
             self.conversion_factor_secondary_to_primary_cooling(self._get_monthly_eer(True, part_load)))
@@ -591,6 +696,10 @@ def _monthly_peak_extraction_heating_simulation_period(self) -> np.ndarray:
             Peak extraction for the whole simulation period
         """
         part_load = self.monthly_peak_heating_simulation_period
+        if self._limit_to_max_heat_pump_power:
+            return np.multiply(
+                self._get_max_power_heating_monthly(),
+                self.conversion_factor_secondary_to_primary_heating(self._get_monthly_cop(True, part_load)))
         return np.multiply(
             self.monthly_peak_heating_simulation_period,
             self.conversion_factor_secondary_to_primary_heating(self._get_monthly_cop(True, part_load)))
@@ -606,6 +715,10 @@ def _monthly_peak_extraction_dhw_simulation_period(self) -> np.ndarray:
             Peak extraction for the whole simulation period
         """
         part_load_dhw = self.monthly_peak_dhw_simulation_period
+        if self._limit_to_max_heat_pump_power:
+            return np.multiply(
+                self._get_max_power_dhw_monthly(),
+                self.conversion_factor_secondary_to_primary_heating(self._get_monthly_cop_dhw(True, part_load_dhw)))
         return np.multiply(
             self.monthly_peak_dhw_simulation_period,
             self.conversion_factor_secondary_to_primary_heating(self._get_monthly_cop_dhw(True, part_load_dhw)))
@@ -693,6 +806,40 @@ def yearly_cooling_load_simulation_period(self) -> np.array:
         return np.sum(np.reshape(self.monthly_baseload_cooling_simulation_period, (self.simulation_period, 12)),
                       axis=1)
 
+    @property
+    def _yearly_heating_load_simulation_period_with_limit(self) -> np.array:
+        """
+        This function returns the yearly heating demand in kWh/year actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        yearly heating : np.ndarray
+            yearly heating for the whole simulation period
+        """
+        if not self._limit_to_max_heat_pump_power:
+            return self.yearly_heating_load_simulation_period
+        return np.sum(
+            np.reshape(self._monthly_baseload_heating_simulation_period_with_limit, (self.simulation_period, 12)),
+            axis=1)
+
+    @property
+    def _yearly_cooling_load_simulation_period_with_limit(self) -> np.array:
+        """
+        This function returns the yearly cooling demand in kWh/year actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        yearly cooling : np.ndarray
+            yearly cooling for the whole simulation period
+        """
+        if not self._limit_to_max_heat_pump_power:
+            return self.yearly_cooling_load_simulation_period
+        return np.sum(
+            np.reshape(self._monthly_baseload_cooling_simulation_period_with_limit, (self.simulation_period, 12)),
+            axis=1)
+
     @property
     def yearly_heating_peak_simulation_period(self) -> np.array:
         """
@@ -872,6 +1019,22 @@ def yearly_dhw_load_simulation_period(self) -> np.array:
         """
         return np.sum(np.reshape(self.monthly_baseload_dhw_simulation_period, (self.simulation_period, 12)), axis=1)
 
+    @property
+    def _yearly_dhw_load_simulation_period_with_limit(self) -> np.array:
+        """
+        This function returns the yearly DHW demand in kWh/year actually provided by the
+        heat pump taking into account the maximum power for the whole simulation period.
+
+        Returns
+        -------
+        yearly DHW : np.ndarray
+            yearly DHW for the whole simulation period
+        """
+        if not self._limit_to_max_heat_pump_power:
+            return self.yearly_dhw_load_simulation_period
+        return np.sum(np.reshape(self._monthly_baseload_dhw_simulation_period_with_limit, (self.simulation_period, 12)),
+                      axis=1)
+
     @property
     def yearly_electricity_consumption(self) -> np.ndarray:
         """
@@ -894,7 +1057,7 @@ def yearly_electricity_consumption_cooling(self) -> np.ndarray:
         -------
         Yearly cooling electricity consumption : np.ndarray
         """
-        return -self.yearly_cooling_load_simulation_period + self.yearly_injection_load_simulation_period
+        return -self._yearly_cooling_load_simulation_period_with_limit + self.yearly_injection_load_simulation_period
 
     @property
     def yearly_electricity_consumption_heating(self) -> np.ndarray:
@@ -905,7 +1068,7 @@ def yearly_electricity_consumption_heating(self) -> np.ndarray:
         -------
         Yearly heating electricity consumption : np.ndarray
         """
-        heating = np.sum(np.reshape(self.monthly_baseload_heating_simulation_period,
+        heating = np.sum(np.reshape(self._monthly_baseload_heating_simulation_period_with_limit,
                                     (self.simulation_period, 12)), axis=1)
         extraction = np.sum(np.reshape(self._monthly_baseload_extraction_heating_simulation_period,
                                        (self.simulation_period, 12)), axis=1)
@@ -920,7 +1083,7 @@ def yearly_electricity_consumption_dhw(self) -> np.ndarray:
         -------
         Yearly DHW electricity consumption : np.ndarray
         """
-        dhw = np.sum(np.reshape(self.monthly_baseload_dhw_simulation_period,
+        dhw = np.sum(np.reshape(self._monthly_baseload_dhw_simulation_period_with_limit,
                                 (self.simulation_period, 12)), axis=1)
         extraction = np.sum(np.reshape(self._monthly_baseload_extraction_dhw_simulation_period,
                                        (self.simulation_period, 12)), axis=1)
@@ -936,7 +1099,8 @@ def SEER(self) -> float:
         SEER : float
         """
         # negative sign to get a positive value
-        return self.yearly_average_cooling_load / np.mean(self.yearly_electricity_consumption_cooling)
+        return np.mean(self._yearly_cooling_load_simulation_period_with_limit) / np.mean(
+            self.yearly_electricity_consumption_cooling)
 
     @property
     def SCOP_total(self) -> float:
@@ -947,7 +1111,8 @@ def SCOP_total(self) -> float:
         -------
         SCOP : float
         """
-        return (self.yearly_average_heating_load + self.yearly_average_dhw_load) / \
+        return (np.mean(self._yearly_heating_load_simulation_period_with_limit) + np.mean(
+            self._yearly_dhw_load_simulation_period_with_limit)) / \
             np.mean(self.yearly_electricity_consumption_heating + self.yearly_electricity_consumption_dhw)
 
     @property
@@ -959,7 +1124,8 @@ def SCOP_heating(self) -> float:
         -------
         SCOP : float
         """
-        return self.yearly_average_heating_load / np.mean(self.yearly_electricity_consumption_heating)
+        return np.mean(self._yearly_heating_load_simulation_period_with_limit) / np.mean(
+            self.yearly_electricity_consumption_heating)
 
     @property
     def SCOP_DHW(self) -> float:
@@ -970,7 +1136,8 @@ def SCOP_DHW(self) -> float:
         -------
         SCOP : float
         """
-        return self.yearly_average_dhw_load / np.mean(self.yearly_electricity_consumption_dhw)
+        return np.mean(self._yearly_dhw_load_simulation_period_with_limit) / np.mean(
+            self.yearly_electricity_consumption_dhw)
 
     @property
     def yearly_SEER(self) -> np.ndarray:
@@ -982,7 +1149,7 @@ def yearly_SEER(self) -> np.ndarray:
         SEER : np.ndarray
         """
         # negative sign to get a positive value
-        return self.yearly_cooling_load_simulation_period / self.yearly_electricity_consumption_cooling
+        return self._yearly_cooling_load_simulation_period_with_limit / self.yearly_electricity_consumption_cooling
 
     @property
     def yearly_SCOP_total(self) -> np.ndarray:
@@ -993,7 +1160,8 @@ def yearly_SCOP_total(self) -> np.ndarray:
         -------
         SCOP : np.ndarray
         """
-        return (self.yearly_heating_load_simulation_period + self.yearly_dhw_load_simulation_period) / (
+        return (
+                self._yearly_heating_load_simulation_period_with_limit + self._yearly_dhw_load_simulation_period_with_limit) / (
                 self.yearly_electricity_consumption_heating + self.yearly_electricity_consumption_dhw)
 
     @property
@@ -1005,7 +1173,7 @@ def yearly_SCOP_heating(self) -> np.ndarray:
         -------
         SCOP : np.ndarray
         """
-        return self.yearly_heating_load_simulation_period / self.yearly_electricity_consumption_heating
+        return self._yearly_heating_load_simulation_period_with_limit / self.yearly_electricity_consumption_heating
 
     @property
     def yearly_SCOP_DHW(self) -> np.ndarray:
@@ -1016,7 +1184,7 @@ def yearly_SCOP_DHW(self) -> np.ndarray:
         -------
         SCOP : np.ndarray
         """
-        return self.yearly_dhw_load_simulation_period / self.yearly_electricity_consumption_dhw
+        return self._yearly_dhw_load_simulation_period_with_limit / self.yearly_electricity_consumption_dhw
 
     @property
     def simulation_period(self) -> int:
diff --git a/GHEtool/test/general_tests/test_GHEtool.py b/GHEtool/test/general_tests/test_GHEtool.py
index 4acc4665..031991e2 100644
--- a/GHEtool/test/general_tests/test_GHEtool.py
+++ b/GHEtool/test/general_tests/test_GHEtool.py
@@ -499,4 +499,4 @@ def test_case_issue_390():
 
     borefield.borehole = borehole
     borefield.size_L3()
-    assert np.isclose(borefield.H, 86.20617039094199)
+    assert np.isclose(borefield.H, 86.3465965114525)
diff --git a/GHEtool/test/general_tests/test_multiple_years.py b/GHEtool/test/general_tests/test_multiple_years.py
index 498c3e2f..f761ad58 100644
--- a/GHEtool/test/general_tests/test_multiple_years.py
+++ b/GHEtool/test/general_tests/test_multiple_years.py
@@ -29,7 +29,7 @@ def test_multiple_years_L4():
                        load["cooling"].clip(0) * load_factor - load["heating"].clip(0) * load_factor)
     h = borefield.size_L4(150)
     assert np.isclose(h, 115.26627395156956, rtol=0.001)
-    load["heating"][8760 * 25:] = 0
+    load.loc[8760 * 25:, "heating"] = 0
     hourly_load = HourlyGeothermalLoadMultiYear(load["heating"].clip(0) * load_factor,
                                                 load["cooling"].clip(0) * load_factor)
     borefield.load = hourly_load
@@ -73,7 +73,7 @@ def test_multiple_years_L3():
                        monthly_cooling_load - monthly_heating_load)
     h = borefield.size_L3(150)
     assert np.isclose(h, 110.18056686276053, rtol=0.001)
-    load["heating"][8760 * 25:] = 0
+    load.loc[8760 * 25:, "heating"] = 0
     hourly_load = HourlyGeothermalLoadMultiYear(load["heating"].clip(0) * load_factor,
                                                 load["cooling"].clip(0) * load_factor)
     borefield.load = hourly_load
@@ -128,7 +128,7 @@ def test_multiple_years_L3_monthly_data():
     assert borefield.load.simulation_period == 50
     h = borefield.size_L3(150)
     assert np.isclose(h, 110.18056686276053, rtol=0.001)
-    load["heating"][8760 * 25:] = 0
+    load.loc[8760 * 25:, "heating"] = 0
     monthly_heating_load = np.array([np.mean((load["heating"].clip(0) * load_factor)[i - 730:i]) for i in
                                      range(730, len(load["heating"]) + 1, 730)])
     monthly_cooling_load = np.array([np.mean((load["cooling"].clip(0) * load_factor)[i - 730:i]) for i in
diff --git a/GHEtool/test/methods/method_data.py b/GHEtool/test/methods/method_data.py
index 00aa152c..58da11e7 100644
--- a/GHEtool/test/methods/method_data.py
+++ b/GHEtool/test/methods/method_data.py
@@ -826,7 +826,7 @@
 borefield.calculation_setup(size_based_on='outlet')
 borefield.flow_data = ConstantDeltaTFlowRate(extraction=4, injection=4)
 list_of_test_objects.add(
-    SizingObject(borefield, L3_output=89.45718988593278, L4_output=73.49659696896121, quadrant=1,
+    SizingObject(borefield, L3_output=96.13895180044284, L4_output=81.15952849537456, quadrant=1,
                  name='BS2023 Auditorium, (var temp and flow, outlet)'))
 borefield.calculation_setup(size_based_on='average')
 borefield.flow_data = ConstantDeltaTFlowRate(extraction=4, injection=4)
@@ -1199,9 +1199,9 @@
 borefield.calculation_setup(size_based_on='average')
 borefield.fluid_data = TemperatureDependentFluidData('MPG', 25)
 borefield.flow_data = ConstantDeltaTFlowRate()
-list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 24.214674441763602, 25.753989837964735,
-                                                   41.116952098990666, 34.79674508776775, 481.21353320134585,
-                                                   643.2765132735545,
+list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 24.212056052962087, 25.750678943889987,
+                                                   41.1119786760319, 34.791392610197995, 481.2201644319575,
+                                                   643.28161087124,
                                                    name='Optimise load profile (power, average, var flow)', power=1,
                                                    hourly=False))
 list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 29.469091502295186, 38.59024834792813,
@@ -1209,15 +1209,15 @@
                                                    642.3540174513246,
                                                    name='Optimise load profile (energy, average, var flow)', power=2,
                                                    hourly=False))
-list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 13.698017382258854, 23.47795628412292,
-                                                   22.175763502822953, 31.189147914587412, 506.46845132956946,
-                                                   646.7123201051548,
+list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 13.696484501705566, 23.474869195054694,
+                                                   22.173081167540246, 31.18435036199417, 506.47202777661306,
+                                                   646.7168892028627,
                                                    name='Optimise load profile (balance, average, var flow)', power=3,
                                                    hourly=False))
 borefield.calculation_setup(size_based_on='inlet')
-list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 17.642440693076743, 15.266828929915876,
-                                                   29.077975024702067, 19.09436975853684, 497.26550263373065,
-                                                   658.2311564442506,
+list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 17.477913367656956, 15.1900163253187,
+                                                   28.790074281883236, 18.988083079232872, 497.6493702908224,
+                                                   658.3323818531115,
                                                    name='Optimise load profile (power, inlet, var flow)', power=1,
                                                    hourly=False))
 list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 28.4206346064365, 29.7813899030333,
@@ -1225,8 +1225,8 @@
                                                    654.9042157277862,
                                                    name='Optimise load profile (energy, inlet, var flow)', power=2,
                                                    hourly=False))
-list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 0.10398122055215273, 0.17763821367920854,
-                                                   0.12212612538444671, 0.17867050911056828, 535.8733011661541,
-                                                   676.2461081103709,
+list_of_test_objects.add(OptimiseLoadProfileObject(borefield, load, 146, 0.10201576191504061, 0.17673320915080462,
+                                                   0.11971975824374738, 0.17771150017503234, 535.876509655675,
+                                                   676.2470214522142,
                                                    name='Optimise load profile (balance, inlet, var flow)', power=3,
                                                    hourly=False))
diff --git a/GHEtool/test/test_examples.py b/GHEtool/test/test_examples.py
index 14c53954..50df96e2 100644
--- a/GHEtool/test/test_examples.py
+++ b/GHEtool/test/test_examples.py
@@ -55,15 +55,22 @@ def test_sizing_with_building_load(monkeypatch):
         size_with_variable_ground_temperature, \
         size_with_part_load_data
     assert np.allclose(size_with_scop(), (95.94350214711557, 4.072466974615784))
-    assert np.allclose(size_with_variable_ground_temperature(), (95.27809127378154, 4.135583487407119))
-    assert np.allclose(size_with_part_load_data(), (97.78277317589445, 4.627161451859563))
+    assert np.allclose(size_with_variable_ground_temperature(), (94.91607600186884, 4.140667366915088))
+    assert np.allclose(size_with_part_load_data(), (97.40238165884419, 4.633232724880249))
 
 
 def test_sizing_with_building_load_hourly(monkeypatch):
     monkeypatch.setattr(plt, 'show', lambda: None)
     from GHEtool.Examples.sizing_with_building_load_hourly import L3_sizing, L4_sizing
-    assert np.allclose(L3_sizing(), (128.14624698264709, 6.291201457558123))
-    assert np.allclose(L4_sizing(), (153.91434358641897, 6.3879487848981364))
+    assert np.allclose(L3_sizing(), (128.0601612984388, 6.292983056938187))
+    assert np.allclose(L4_sizing(), (153.81812162652483, 6.387522978931683))
+
+
+def test_sizing_with_building_load_hourly_limit(monkeypatch):
+    monkeypatch.setattr(plt, 'show', lambda: None)
+    from GHEtool.Examples.sizing_with_building_load_hourly_limit import L3_sizing, L4_sizing
+    assert np.allclose(L3_sizing(), 111.8121680540246)
+    assert np.allclose(L4_sizing(), (144.11894710530277, 5.027207193115156))
 
 
 def test_combined_active_and_passive_cooling(monkeypatch):
diff --git a/GHEtool/test/unit-tests/test_efficiency.py b/GHEtool/test/unit-tests/test_efficiency.py
index e1acfc0d..15b5a749 100644
--- a/GHEtool/test/unit-tests/test_efficiency.py
+++ b/GHEtool/test/unit-tests/test_efficiency.py
@@ -6,6 +6,101 @@
 import numpy as np
 
 from GHEtool.VariableClasses.Efficiency import *
+from GHEtool.VariableClasses.Efficiency._Efficiency import plot_heat_pump_envelope, combine_n_heat_pumps, \
+    _find_optimal_heat_pump_configuration
+
+points_HP300 = np.array([
+    [-4.5, 53.6],
+    [-4.5, 39.8],
+    [-4.5, 23.2],
+    [-1.5, 58.7],
+    [-1.5, 43.7],
+    [-1.5, 25.5],
+    [3.5, 67.9],
+    [3.5, 50.7],
+    [3.5, 29.6],
+    [8.5, 77.6],
+    [8.5, 58.3],
+    [8.5, 33.9],
+    [11.5, 83.9],
+    [11.5, 62.9],
+    [11.5, 36.6],
+])
+eff_HP300 = np.array([
+    3.86, 4.28, 3.93,
+    4.08, 4.60, 4.18,
+    4.53, 5.12, 4.70,
+    5.04, 5.72, 5.30,
+    5.41, 6.17, 5.63,
+])
+points_HP300_new = np.array([
+    [-1.5, 58.7],
+    [-1.5, 43.7],
+    [-1.5, 25.5],
+    [3.5, 67.9],
+    [3.5, 50.7],
+    [3.5, 29.6],
+    [8.5, 77.6],
+    [8.5, 58.3],
+    [8.5, 33.9],
+    [11.5, 83.9],
+    [11.5, 62.9],
+    [11.5, 36.6],
+])
+eff_HP300_new = np.array([
+    4.08, 4.60, 4.18,
+    4.53, 5.12, 4.70,
+    5.04, 5.72, 5.30,
+    5.41, 6.17, 5.63,
+])
+points_HP400 = np.array([
+    [-4.5, 67.2],
+    [-4.5, 49.9],
+    [-4.5, 29.1],
+    [-1.5, 73.7],
+    [-1.5, 54.8],
+    [-1.5, 32.0],
+    [3.5, 85.2],
+    [3.5, 63.5],
+    [3.5, 37.2],
+    [8.5, 97.7],
+    [8.5, 73.1],
+    [8.5, 42.7],
+    [11.5, 105.7],
+    [11.5, 79.0],
+    [11.5, 46.0],
+])
+eff_HP400 = np.array([
+    3.91, 4.38, 3.99,
+    4.14, 4.64, 4.27,
+    4.58, 5.16, 4.77,
+    5.12, 5.80, 5.34,
+    5.51, 6.22, 5.75,
+])
+points_HP500 = np.array([
+    [-4.5, 84.7],
+    [-4.5, 66.3],
+    [-4.5, 38.6],
+    [-1.5, 93.0],
+    [-1.5, 72.9],
+    [-1.5, 42.4],
+    [3.5, 107.8],
+    [3.5, 84.6],
+    [3.5, 49.4],
+    [8.5, 123.5],
+    [8.5, 97.4],
+    [8.5, 56.9],
+    [11.5, 133.2],
+    [11.5, 105.4],
+    [11.5, 61.4],
+])
+eff_HP500 = np.array([
+    3.87, 4.12, 3.78,
+    4.12, 4.39, 4.04,
+    4.57, 4.89, 4.49,
+    5.10, 5.50, 5.08,
+    5.48, 5.92, 5.48,
+])
 
 
 def test_SCOP():
@@ -19,6 +114,7 @@ def test_SCOP():
 def test_SEER():
     seer = SEER(50)
     assert seer.get_SEER(12, test=5) == 50
+    assert seer._get_max_power()==1e16
     assert seer.get_EER(12, test=5) == 50
     with pytest.raises(ValueError):
         SEER(0)
@@ -49,6 +145,7 @@ def test_COP_basic():
     assert np.array_equal(cop_basic.get_COP(np.array([7.5, 10]), power=5), np.array([3.5, 4]))
     assert np.array_equal(cop_basic.get_COP(np.array([7.5, 10, 0]), power=np.array([5, 6])),
                           np.array([3.5, 4, 3]))
+    assert np.allclose(cop_basic._get_max_power(np.array([1, 1.5, 2, 2.5, 3])), 1e16)
 
 
 def test_COP_secondary():
@@ -68,6 +165,7 @@ def test_COP_secondary():
     assert np.array_equal(cop_sec.get_COP(np.array([1.5, 3]), np.array([3, 4]), np.array([1, 2])),
                           np.array([1.25, 3.5]))
     assert np.array_equal(cop_sec.get_COP(1.5, np.array([2.5, 4.5])), np.array([1, 2]))
+    assert np.allclose(cop_sec._get_max_power(np.array([1, 1.5, 2, 2.5, 3])), 1e16)
 
 
 def test_COP_part_load():
@@ -79,6 +177,8 @@ def test_COP_part_load():
         assert cop_part.get_COP(5) == 3
 
     assert cop_part.get_COP(1.5, power=3.5) == 1.5
+    assert cop_part.get_COP(1.5, power=5.5) == 2
+
     assert np.array_equal(cop_part.get_COP(np.array([1.5, 2]), power=np.array([3, 4])), np.array([1.25, 2.625]))
     assert np.array_equal(cop_part.get_COP(np.array([1.5, 3]), power=np.array([3, 4])), np.array([1.25, 3.5]))
 
@@ -90,6 +190,24 @@ def test_COP_part_load():
         np.array([1.25, 3.5]))
     assert np.array_equal(cop_part.get_COP(1.5, power=np.array([2.5, 4.5])), np.array([1, 2]))
 
+    assert cop_part._get_max_power(1.5) == 4.5
+    assert cop_part._get_max_power(1) == 4.5
+    cop_part = COP(np.array([1, 2, 2, 4]), np.array([[1.5, 2.5], [2.5, 2.5], [1.5, 4.5], [2.5, 5.5]]), part_load=True)
+    assert np.allclose(cop_part._get_max_power(np.array([1, 1.5, 2, 2.5, 3])), [4.5, 4.5, 5, 5.5, 5.5])
+
+
+def test_COP_part_load_real():
+    cop = COP(np.array(
+        [4.42, 5.21, 6.04, 7.52, 9.5, 3.99, 4.58, 5.21, 6.02, 6.83, 3.86, 4.39, 4.97,
+         5.62, 6.19, 3.8, 4.3, 4.86, 5.44, 5.9, 3.76, 4.25, 4.79, 5.34, 5.74]),
+        np.array([[-5, 1.06], [0, 1.25], [5, 1.45], [10, 1.66], [15, 1.9], [-5, 2.05], [0, 2.42], [5, 2.81], [10, 3.2],
+                  [15, 3.54], [-5, 3.05], [0, 3.6], [5, 4.17], [10, 4.73], [15, 5.18], [-5, 4.04], [0, 4.77], [5, 5.54],
+                  [10, 6.27], [15, 6.82], [-5, 5.03], [0, 5.95], [5, 6.9], [10, 7.81], [15, 8.46]]),
+        part_load=True)
+    assert np.isclose(cop.get_COP(0, power=1.25), 5.21)
+    assert np.allclose(cop.get_COP(np.array([0, 5]), power=np.array([1.25, 1.45])), [5.21, 6.04])
+    assert np.allclose(cop._get_max_power(np.array([0, 5])), [5.95, 6.9])
+
 
 def test_COP_full():
     cop_full = COP(np.array([1, 2, 2, 4, 2, 4, 4, 8]),
@@ -111,6 +229,16 @@ def test_COP_full():
     assert np.array_equal(cop_full.get_COP(1.5, secondary_temperature=np.array([2.5, 4.5]), power=4.5),
                           np.array([1, 2]))
 
+    with pytest.raises(ValueError):
+        assert cop_full._get_max_power(1.5)
+    assert cop_full._get_max_power(1.5, 4.5) == 8.5
+    cop_full = COP(np.array([1, 2, 2, 4, 2, 4, 4, 8]),
+                   np.array([[1.5, 2.5, 4.5], [2.5, 2.5, 4.5], [1.5, 4.5, 3.5], [2.5, 4.5, 5.5],
+                             [1.5, 2.5, 8.5], [2.5, 2.5, 10.5], [1.5, 4.5, 7.5], [2.5, 4.5, 9.5]]),
+                   secondary=True, part_load=True)
+    assert np.allclose(cop_full._get_max_power(np.array([1.5, 1.5, 2.5, 2.5]), np.array([2.5, 4.5, 2.5, 4.5])),
+                       [8.5, 7.5, 10.5, 9.5])
+
 
 def test_COP_get_SCOP():
     cop_full = COP(np.array([1, 2, 2, 4, 2, 4, 4, 8]),
@@ -351,6 +479,53 @@ def test_EERCombined():
     assert np.isclose(eer.get_SEER(np.array([10, 10, 10]), np.array([1, 15, 20]), month_indices=np.array([5, 6, 7])),
                       30 / 4.5)
 
+def test_EERCombined_max_power():
+    # with threshold
+    eer = EERCombined(20, 5, 10)
+    assert eer._get_max_power(1, 0, 0) == 1e16
+    assert np.allclose(eer._get_max_power(np.array([1, 10, 20])), np.array([1e16, 1e16, 1e16]))
+
+    # with month array
+    eer = EERCombined(20, 5, months_active_cooling=np.array([7, 8, 9]))
+    assert eer._get_max_power(1, 0, 0) == 1e16
+    assert np.allclose(eer._get_max_power(np.array([1, 10, 20]), month_indices=np.array([6, 7,8])), np.array([1e16, 1e16, 1e16]))
+
+    # with threshold and month array
+    eer = EERCombined(20, 5, 10, months_active_cooling=np.array([7, 8, 9]))
+    assert eer._get_max_power(1, 0, 0) == 1e16
+    assert np.allclose(eer._get_max_power(np.array([1, 10, 20]), month_indices=np.array([6, 7,8])), np.array([1e16, 1e16, 1e16]))
+
+    eer_part = EER(np.array([1, 2, 2, 4]), np.array([[1.5, 2.5], [2.5, 2.5], [1.5, 4.5], [2.5, 4.5]]), part_load=True)
+    # with threshold
+    eer = EERCombined(20, eer_part, 10)
+    assert eer._get_max_power(1, 0, 0) == 1e16
+    assert eer._get_max_power(11, 0, 0) == 4.5
+    assert np.allclose(eer._get_max_power(np.array([1, 11, 20])), np.array([1e16,4.5,4.5]))
+
+    eer = EERCombined(eer_part, 20, 10)
+    assert eer._get_max_power(1, 0, 0) == 4.5
+    assert eer._get_max_power(11, 0, 0) == 1e16
+    assert np.allclose(eer._get_max_power(np.array([1, 11, 20])), np.array([4.5,1e16,1e16]))
+
+    # with month array
+    eer = EERCombined(20, eer_part, months_active_cooling=np.array([7, 8, 9]))
+    with pytest.raises(ValueError):
+        eer._get_max_power(1)
+    with pytest.raises(ValueError):
+        eer._get_max_power(np.array([1, 10, 20]))
+    assert eer._get_max_power(1, 0, 0) == 1e16
+    assert eer._get_max_power(11, 0, 0) == 1e16
+    assert np.allclose(eer._get_max_power(np.array([1, 11, 20]), month_indices=np.array([6, 7,8])), np.array([1e16, 4.5,4.5]))
+    eer = EERCombined(eer_part, 10, months_active_cooling=np.array([7, 8, 9]))
+    assert eer._get_max_power(1, 0, 0) == 4.5
+    assert eer._get_max_power(11, 0, 0) == 4.5
+    assert np.allclose(eer._get_max_power(np.array([1, 11, 20]), month_indices=np.array([6, 7,8])), np.array([4.5, 1e16, 1e16]))
+    assert np.allclose(eer._get_max_power(1, month_indices=np.array([6, 7,8])), np.array([4.5, 1e16, 1e16]))
+
+    # with threshold and month array
+    eer = EERCombined(20, eer_part, 10, months_active_cooling=np.array([7, 8, 9]))
+    assert eer._get_max_power(1, 0, 0) == 1e16
+    assert np.allclose(eer._get_max_power(np.array([1, 10, 20]), month_indices=np.array([6, 7,8])), np.array([1e16, 4.5, 4.5]))
 
 def test_eq_eer_combined():
     eer_combined = EERCombined(20, 5, 10)
@@ -399,3 +574,27 @@ def test_repr_():
     eer_combined = EERCombined(20, 5, months_active_cooling=[5, 6])
     assert {'Active cooling': {'SEER [-]': 5}, 'Passive cooling': {'SEER [-]': 20},
             'With active cooling in months [-]': [5, 6]} == eer_combined.__export__()
+
+
+def test_graph_efficiency(monkeypatch):
+    import matplotlib.pyplot as plt
+    monkeypatch.setattr(plt, "show", lambda: None)
+    plot_heat_pump_envelope(points_HP300, eff_HP300)
+
+
+def test_combine_heat_pumps():
+    combine_n_heat_pumps([points_HP300, points_HP300, points_HP300], [eff_HP300, eff_HP300, eff_HP300])
+    combine_n_heat_pumps([points_HP300], [eff_HP300])
+    combine_n_heat_pumps([points_HP300, points_HP400], [eff_HP300, eff_HP400])
+    combine_n_heat_pumps([points_HP300, points_HP300_new], [eff_HP300, eff_HP300_new])
+
+
+def test_find_optimal_heat_pump_configuration():
+    hp_300 = COP(eff_HP300, points_HP300, part_load=True)
+    hp_400 = COP(eff_HP400, points_HP400, part_load=True)
+    hp_500 = COP(eff_HP500, points_HP500, part_load=True)
+
+    assert np.allclose(_find_optimal_heat_pump_configuration([hp_300, hp_400, hp_500], 30, prim_temp=-1), [1, 0, 0])
+    assert np.allclose(_find_optimal_heat_pump_configuration([hp_300, hp_400, hp_500], 70, prim_temp=-1), [0, 1, 0])
+    assert np.allclose(_find_optimal_heat_pump_configuration([hp_300, hp_400, hp_500], 80, prim_temp=-1), [0, 0, 1])
+    assert np.allclose(_find_optimal_heat_pump_configuration([hp_300, hp_500, hp_400], 120, prim_temp=-1), [1, 0, 1])
diff --git a/GHEtool/test/unit-tests/test_loads/test_hourly_load_building.py b/GHEtool/test/unit-tests/test_loads/test_hourly_load_building.py
index 96db0b8a..d9da3daa 100644
--- a/GHEtool/test/unit-tests/test_loads/test_hourly_load_building.py
+++ b/GHEtool/test/unit-tests/test_loads/test_hourly_load_building.py
@@ -118,6 +118,10 @@ def test_baseload_heating():
     assert np.array_equal(load.monthly_baseload_heating, np.linspace(0, 11, 12) * 730)
     assert np.array_equal(load.monthly_baseload_heating / 730, load.monthly_baseload_heating_power)
     assert np.array_equal(load.monthly_baseload_heating_power, load.monthly_peak_heating)
+    load._limit_to_max_heat_pump_power = True
+    assert np.array_equal(load.monthly_baseload_heating, np.linspace(0, 11, 12) * 730)
+    assert np.array_equal(load.monthly_baseload_heating / 730, load.monthly_baseload_heating_power)
+    assert np.array_equal(load.monthly_baseload_heating_power, load.monthly_peak_heating)
 
 
 def test_baseload_cooling():
@@ -127,6 +131,10 @@ def test_baseload_cooling():
     assert np.array_equal(load.monthly_baseload_cooling, np.linspace(0, 11, 12) * 730)
     assert np.array_equal(load.monthly_baseload_cooling / 730, load.monthly_baseload_cooling_power)
     assert np.array_equal(load.monthly_baseload_cooling_power, load.monthly_peak_cooling)
+    load._limit_to_max_heat_pump_power = True
+    assert np.array_equal(load.monthly_baseload_cooling, np.linspace(0, 11, 12) * 730)
+    assert np.array_equal(load.monthly_baseload_cooling / 730, load.monthly_baseload_cooling_power)
+    assert np.array_equal(load.monthly_baseload_cooling_power, load.monthly_peak_cooling)
 
 
 def test_peak_heating():
@@ -138,6 +146,13 @@ def test_peak_heating():
     load.hourly_heating_load = np.repeat(np.linspace(1, 12, 12), 730)
     assert np.array_equal(load.monthly_peak_heating, np.linspace(1, 12, 12))
     assert np.isclose(load.max_peak_heating, 12)
+    load._limit_to_max_heat_pump_power = True
+    assert np.array_equal(load.monthly_peak_heating, np.linspace(1, 12, 12))
+    assert np.isclose(load.max_peak_heating, 12)
+    cop_part = COP(np.array([1, 2, 2, 4]), np.array([[1.5, 2.5], [2.5, 2.5], [1.5, 4.5], [2.5, 4.5]]), part_load=True)
+    load.cop = cop_part
+    assert np.array_equal(load.monthly_peak_heating, np.linspace(1, 12, 12))
+    assert np.isclose(load.max_peak_heating, 12)
 
 
 def test_peak_cooling():
@@ -149,6 +164,13 @@ def test_peak_cooling():
     load.hourly_cooling_load = np.repeat(np.linspace(1, 12, 12), 730)
     assert np.array_equal(load.monthly_peak_cooling, np.linspace(1, 12, 12))
     assert np.isclose(load.max_peak_cooling, 12)
+    load._limit_to_max_heat_pump_power = True
+    assert np.array_equal(load.monthly_peak_cooling, np.linspace(1, 12, 12))
+    assert np.isclose(load.max_peak_cooling, 12)
+    eer_part = EER(np.array([1, 2, 2, 4]), np.array([[1.5, 2.5], [2.5, 2.5], [1.5, 4.5], [2.5, 4.5]]), part_load=True)
+    load.eer = eer_part
+    assert np.array_equal(load.monthly_peak_cooling, np.linspace(1, 12, 12))
+    assert np.isclose(load.max_peak_cooling, 12)
 
 
 def test_load_simulation_period():
@@ -437,6 +459,69 @@ def test_hourly_injection_load_simulation_period_hourly_data():
                            np.tile(np.concatenate((np.full(4380, 6.25), np.full(4380, 10.5))), 10))[0])
 
 
+def test_hourly_injection_load_simulation_period_hourly_data_limit():
+    load = HourlyBuildingLoad(np.zeros(8760), np.zeros(8760), 10, scop, seer)
+    load.hourly_cooling_load = test_load
+    load._limit_to_max_heat_pump_power = True
+    assert np.allclose(load.hourly_injection_load_simulation_period, test_load_sim_per * 6 / 5)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 6 / 5)[1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 6 / 5)[0])
+    load.eer = eer_basic
+    load.reset_results(0, 1)
+    assert np.allclose(load.hourly_injection_load_simulation_period, test_load_sim_per * 3 / 2)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 3 / 2)[1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 3 / 2)[0])
+    load.reset_results(0, 10)
+    assert np.allclose(load.hourly_injection_load_simulation_period, test_load_sim_per * 21 / 20)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 21 / 20)[1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 21 / 20)[0])
+    load.set_results(results_hourly_test)
+    assert np.allclose(load.hourly_injection_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 7.5), np.full(4380, 11))), 10))
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(np.tile(np.concatenate((np.full(4380, 7.5), np.full(4380, 11))), 10))[
+                           1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(np.tile(np.concatenate((np.full(4380, 7.5), np.full(4380, 11))), 10))[
+                           0])
+    load.eer = eer_pl
+    load.reset_results(0, 1)
+    assert np.allclose(load.hourly_injection_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 6.25 / 5), np.full(4380, 12.5 / 10))), 10))
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 6.25 / 5), np.full(4380, 12.5 / 10))), 10))[
+                           1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 6.25 / 5), np.full(4380, 12.5 / 10))), 10))[
+                           0])
+    load.reset_results(0, 10)
+    assert np.allclose(load.hourly_injection_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 5.125 / 5), np.full(4380, 10.25 / 10))), 10))
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 5.125 / 5), np.full(4380, 10.25 / 10))), 10))[1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 5.125 / 5), np.full(4380, 10.25 / 10))), 10))[0])
+    load.set_results(results_hourly_test)
+    assert np.allclose(load.hourly_injection_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 6.25 / 5), np.full(4380, 10.5 / 10))), 10))
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 6.25 / 5), np.full(4380, 10.5 / 10))), 10))[1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 6.25 / 5), np.full(4380, 10.5 / 10))), 10))[0])
+
+
 def test_hourly_extraction_load_simulation_period_hourly_data():
     load = HourlyBuildingLoad(np.zeros(8760), np.zeros(8760), 10, scop, seer)
     load.hourly_heating_load = test_load
@@ -453,7 +538,7 @@ def test_hourly_extraction_load_simulation_period_hourly_data():
                        load.resample_to_monthly(test_load_sim_per * 1 / 2)[1])
     assert np.allclose(load.monthly_peak_extraction_simulation_period,
                        load.resample_to_monthly(test_load_sim_per * 1 / 2)[0])
-    load.reset_results(0, 10)
+    load.reset_results(10, 10)
     assert np.allclose(load.hourly_extraction_load_simulation_period, test_load_sim_per * 19 / 20)
     assert np.allclose(load.monthly_baseload_extraction_simulation_period,
                        load.resample_to_monthly(test_load_sim_per * 19 / 20)[1])
@@ -482,7 +567,7 @@ def test_hourly_extraction_load_simulation_period_hourly_data():
                        load.resample_to_monthly(
                            np.tile(np.concatenate((np.full(4380, 3.75), np.full(4380, 7.5))), 10))[
                            0])
-    load.reset_results(0, 10)
+    load.reset_results(10, 10)
     assert np.allclose(load.hourly_extraction_load_simulation_period,
                        np.tile(np.concatenate((np.full(4380, 4.875), np.full(4380, 9.75))), 10))
     assert np.allclose(load.monthly_baseload_extraction_simulation_period,
@@ -528,7 +613,7 @@ def test_hourly_extraction_load_simulation_period_hourly_data():
                        load.resample_to_monthly(test_load_sim_per * 1 / 2)[1])
     assert np.allclose(load.monthly_peak_extraction_simulation_period,
                        load.resample_to_monthly(test_load_sim_per * 1 / 2)[0])
-    load.reset_results(0, 10)
+    load.reset_results(10, 10)
     assert np.allclose(load.hourly_extraction_load_simulation_period, test_load_sim_per * 19 / 20)
     assert np.allclose(load.monthly_baseload_extraction_simulation_period,
                        load.resample_to_monthly(test_load_sim_per * 19 / 20)[1])
@@ -557,7 +642,7 @@ def test_hourly_extraction_load_simulation_period_hourly_data():
                        load.resample_to_monthly(
                            np.tile(np.concatenate((np.full(4380, 3.75), np.full(4380, 7.5))), 10))[
                            0])
-    load.reset_results(0, 10)
+    load.reset_results(10, 10)
     assert np.allclose(load.hourly_extraction_load_simulation_period,
                        np.tile(np.concatenate((np.full(4380, 4.875), np.full(4380, 9.75))), 10))
     assert np.allclose(load.monthly_baseload_extraction_simulation_period,
@@ -588,6 +673,157 @@ def test_hourly_extraction_load_simulation_period_hourly_data():
     assert np.allclose(load._monthly_peak_extraction_heating_simulation_period, np.zeros(120))
 
 
+def test_hourly_extraction_load_simulation_period_hourly_data_max_power():
+    load = HourlyBuildingLoad(np.zeros(8760), np.zeros(8760), 10, scop, seer)
+    load.hourly_heating_load = test_load
+    load._limit_to_max_heat_pump_power = True
+    assert np.allclose(load.hourly_extraction_load_simulation_period, test_load_sim_per * 5 / 6)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 5 / 6)[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 5 / 6)[0])
+    load.cop = cop_basic
+    load.reset_results(0, 1)
+    assert np.allclose(load.hourly_extraction_load_simulation_period, test_load_sim_per * 1 / 2)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 1 / 2)[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 1 / 2)[0])
+    load.reset_results(10, 10)
+    assert np.allclose(load.hourly_extraction_load_simulation_period, test_load_sim_per * 19 / 20)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 19 / 20)[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 19 / 20)[0])
+    load.set_results(results_hourly_test)
+    assert np.allclose(load.hourly_extraction_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[
+                           1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[
+                           0])
+    load.cop = cop_pl
+    load.reset_results(0, 1)
+    assert np.allclose(load.hourly_extraction_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 1 * 3.75 / 5), np.full(4380, 1 * 7.5 / 10))), 10))
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 7.5 / 10))), 10))[
+                           1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 7.5 / 10))), 10))[
+                           0])
+    load.reset_results(10, 10)
+    assert np.allclose(load.hourly_extraction_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 4.875 / 5), np.full(4380, 9.75 / 10))), 10))
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 4.875 / 5), np.full(4380, 9.75 / 10))), 10))[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 4.875 / 5), np.full(4380, 9.75 / 10))), 10))[0])
+    load.set_results(results_hourly_test)
+    assert np.allclose(load.hourly_extraction_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))
+    assert np.allclose(load._hourly_extraction_load_heating_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))
+    assert np.allclose(load._hourly_extraction_load_dhw_simulation_period, np.zeros(87600))
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[0])
+    assert np.allclose(load._monthly_baseload_extraction_heating_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[1])
+    assert np.allclose(load._monthly_peak_extraction_heating_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[0])
+    assert np.allclose(load._monthly_baseload_extraction_dhw_simulation_period, np.zeros(120))
+    assert np.allclose(load._monthly_peak_extraction_dhw_simulation_period, np.zeros(120))
+
+    # now for only DHW
+    load.hourly_heating_load = np.zeros(8760)
+    load.dhw = test_load
+    load.cop_dhw = scop
+    assert np.allclose(load.hourly_extraction_load_simulation_period, test_load_sim_per * 5 / 6)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 5 / 6)[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 5 / 6)[0])
+    load.cop_dhw = cop_basic
+    load.reset_results(0, 1)
+    assert np.allclose(load.hourly_extraction_load_simulation_period, test_load_sim_per * 1 / 2)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 1 / 2)[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 1 / 2)[0])
+    load.reset_results(10, 10)
+    assert np.allclose(load.hourly_extraction_load_simulation_period, test_load_sim_per * 19 / 20)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 19 / 20)[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(test_load_sim_per * 19 / 20)[0])
+    load.set_results(results_hourly_test)
+    assert np.allclose(load.hourly_extraction_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[
+                           1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[
+                           0])
+    load.cop_dhw = cop_pl
+    load.reset_results(0, 1)
+    assert np.allclose(load.hourly_extraction_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 7.5 / 10))), 10))
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 7.5 / 10))), 10))[
+                           1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 7.5 / 10))), 10))[
+                           0])
+    load.reset_results(10, 10)
+    assert np.allclose(load.hourly_extraction_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 4.875 / 5), np.full(4380, 9.75 / 10))), 10))
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 4.875 / 5), np.full(4380, 9.75 / 10))), 10))[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 4.875 / 5), np.full(4380, 9.75 / 10))), 10))[0])
+    load.set_results(results_hourly_test)
+    assert np.allclose(load.hourly_extraction_load_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))
+    assert np.allclose(load._hourly_extraction_load_dhw_simulation_period,
+                       np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))
+    assert np.allclose(load._hourly_extraction_load_heating_simulation_period, np.zeros(87600))
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[0])
+    assert np.allclose(load._monthly_baseload_extraction_dhw_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[1])
+    assert np.allclose(load._monthly_peak_extraction_dhw_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[0])
+    assert np.allclose(load._monthly_baseload_extraction_heating_simulation_period, np.zeros(120))
+    assert np.allclose(load._monthly_peak_extraction_heating_simulation_period, np.zeros(120))
+
+
 def test_hourly_injection_load_simulation_period_monthly_data():
     # these results are the same as with hourly resolution, since the baseload and the peak load have the same
     # temperature result and are also the same power, since the hourly load is a constant, so it reduces to the same
@@ -615,6 +851,47 @@ def test_hourly_injection_load_simulation_period_monthly_data():
                            np.tile(np.concatenate((np.full(4380, 6.25), np.full(4380, 10.5))), 10))[0])
 
 
+def test_hourly_injection_load_simulation_period_monthly_data_limit():
+    # these results are the same as with hourly resolution, since the baseload and the peak load have the same
+    # temperature result and are also the same power, since the hourly load is a constant, so it reduces to the same
+    # power for both the peak and baseload
+    load = HourlyBuildingLoad(np.zeros(8760), np.zeros(8760), 10, scop, seer)
+    load.hourly_cooling_load = test_load
+    load._limit_to_max_heat_pump_power = True
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(np.tile(np.concatenate((np.full(4380, 6),
+                                                                        np.full(4380, 12))), 10))[0])
+    load.set_results(results_monthly_test)
+    load.eer = eer_basic
+    with pytest.raises(TypeError):
+        load.hourly_injection_load_simulation_period
+
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(np.tile(np.concatenate((np.full(4380, 7.5),
+                                                                        np.full(4380, 11))), 10))[1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(np.tile(np.concatenate((np.full(4380, 7.5),
+                                                                        np.full(4380, 11))), 10))[0])
+    assert np.isclose(load.max_peak_injection, 11)
+    load.eer = eer_pl
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 6.25), np.full(4380, 10.5))), 10))[1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 6.25 / 5), np.full(4380, 10.5 / 10))), 10))[0])
+
+    load.reset_results(0, 20)
+    load.eer = eer_basic
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       load.resample_to_monthly(np.tile(np.concatenate((np.full(4380, 5.25),
+                                                                        np.full(4380, 10.5))), 10))[1])
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       load.resample_to_monthly(np.tile(np.concatenate((np.full(4380, 5.25),
+                                                                        np.full(4380, 10.5))), 10))[0])
+    assert np.isclose(load.max_peak_injection, 10.5)
+
+
 def test_hourly_extraction_load_simulation_period_monthly_data():
     # these results are the same as with hourly resolution, since the baseload and the peak load have the same
     # temperature result and are also the same power, since the hourly load is a constant, so it reduces to the same
@@ -682,6 +959,73 @@ def test_hourly_extraction_load_simulation_period_monthly_data():
     assert np.allclose(load._monthly_peak_extraction_heating_simulation_period, np.zeros(120))
 
 
+def test_hourly_extraction_load_simulation_period_monthly_data_limit():
+    # these results are the same as with hourly resolution, since the baseload and the peak load have the same
+    # temperature result and are also the same power, since the hourly load is a constant, so it reduces to the same
+    # power for both the peak and baseload
+    load = HourlyBuildingLoad(np.zeros(8760), np.zeros(8760), 10, scop, seer)
+    load.hourly_heating_load = test_load
+    load._limit_to_max_heat_pump_power = True
+    load.cop = cop_basic
+    load.set_results(results_monthly_test)
+    with pytest.raises(TypeError):
+        load.hourly_extraction_load_simulation_period
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[0])
+    assert np.allclose(load._monthly_baseload_extraction_heating_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[1])
+    assert np.allclose(load._monthly_peak_extraction_heating_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[0])
+    assert np.allclose(load._monthly_baseload_extraction_dhw_simulation_period, np.zeros(120))
+    assert np.allclose(load._monthly_peak_extraction_dhw_simulation_period, np.zeros(120))
+    assert np.isclose(load.max_peak_extraction, 9)
+    assert np.isclose(load.imbalance, -50370.0)
+
+    load.cop = cop_pl
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75), np.full(4380, 9.5))), 10))[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[0])
+
+    # now for only DHW
+    load.hourly_heating_load = np.zeros(8760)
+    load.dhw = test_load
+    load.cop_dhw = cop_basic
+    with pytest.raises(TypeError):
+        load._hourly_extraction_load_dhw_simulation_period
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[
+                           1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 2.5), np.full(4380, 9))), 10))[
+                           0])
+    load.cop_dhw = cop_pl
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75), np.full(4380, 9.5))), 10))[1])
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[0])
+    assert np.allclose(load._monthly_baseload_extraction_dhw_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75), np.full(4380, 9.5))), 10))[1])
+    assert np.allclose(load._monthly_peak_extraction_dhw_simulation_period,
+                       load.resample_to_monthly(
+                           np.tile(np.concatenate((np.full(4380, 3.75 / 5), np.full(4380, 9.5 / 10))), 10))[0])
+    assert np.allclose(load._monthly_baseload_extraction_heating_simulation_period, np.zeros(120))
+    assert np.allclose(load._monthly_peak_extraction_heating_simulation_period, np.zeros(120))
+
+
 def test_time_array():
     load = HourlyBuildingLoad(np.zeros(8760), np.zeros(8760), 10, scop, seer)
     assert np.allclose(load.month_indices, np.tile(np.repeat(np.arange(1, 13), load.UPM), 10))
diff --git a/GHEtool/test/unit-tests/test_loads/test_monthly_load_building.py b/GHEtool/test/unit-tests/test_loads/test_monthly_load_building.py
index eb0595a6..b2740534 100644
--- a/GHEtool/test/unit-tests/test_loads/test_monthly_load_building.py
+++ b/GHEtool/test/unit-tests/test_loads/test_monthly_load_building.py
@@ -11,7 +11,10 @@
 cop_basic = COP(np.array([2, 20]), np.array([1, 10]))
 eer_basic = EER(np.array([2, 20]), np.array([1, 10]))
 cop_pl = COP(np.array([2, 20, 4, 40]), np.array([[1, 5], [10, 5], [1, 10], [10, 10]]), part_load=True)
+cop_pl_2 = COP(np.array([2, 20, 4, 40]), np.array([[1, 5], [10, 5], [1, 5], [10, 5]]), part_load=True)
 eer_pl = EER(np.array([2, 20, 4, 40]), np.array([[1, 5], [10, 5], [1, 10], [10, 10]]), part_load=True)
+eer_pl_2 = EER(np.array([2, 20, 4, 40]), np.array([[1, 5], [10, 5], [1, 5], [10, 5]]), part_load=True)
+
 results_monthly = ResultsMonthly(np.linspace(0, 120 - 1, 120),
                                  np.linspace(0, 120 - 1, 120) * 2,
                                  np.linspace(0, 120 - 1, 120) * 3,
@@ -71,6 +74,7 @@ def test_different_start_month():
     assert np.allclose(load.baseload_dhw, result)
     assert np.allclose(load.monthly_baseload_heating_simulation_period, np.tile(result, 20))
     assert np.allclose(load.monthly_baseload_cooling_simulation_period, np.tile(result, 20))
+    assert np.allclose(load._monthly_baseload_cooling_simulation_period_with_limit, np.tile(result, 20))
     assert np.allclose(load.monthly_peak_heating_simulation_period, np.tile(result, 20))
     assert np.allclose(load.monthly_peak_cooling_simulation_period, np.tile(result, 20))
 
@@ -352,6 +356,31 @@ def test_monthly_baseload_injection_simulation_period():
                        np.tile(np.array(
                            [7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 11, 11, 11, 11, 11, 11]), 10))
 
+    load = MonthlyBuildingLoadAbsolute(*load_case(2), 10, scop, seer)
+    load._limit_to_max_heat_pump_power = True
+    load.baseload_cooling = test_load
+    load.peak_cooling = np.full(12, 10)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period, test_load_sim_per * 6 / 5)
+    load.eer = eer_basic
+    load.reset_results(0, 1)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period, test_load_sim_per * 3 / 2)
+    load.reset_results(0, 10)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period, test_load_sim_per * 21 / 20)
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       np.tile(np.array([7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 11, 11, 11, 11, 11, 11]), 10))
+    load.eer = eer_pl
+    load.reset_results(0, 1)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       np.tile(np.array([7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 15, 15, 15, 15, 15, 15]), 10))
+    load.reset_results(0, 10)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       np.tile(np.array([5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5]), 10))
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_baseload_injection_simulation_period,
+                       np.tile(np.array(
+                           [7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 11, 11, 11, 11, 11, 11]), 10))
+
 
 def test_monthly_baseload_extraction_simulation_period():
     load = MonthlyBuildingLoadAbsolute(*load_case(2), 10, scop, seer)
@@ -409,6 +438,62 @@ def test_monthly_baseload_extraction_simulation_period():
                        np.tile(np.array(
                            [2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 9, 9, 9, 9, 9, 9]), 10))
 
+    load = MonthlyBuildingLoadAbsolute(*load_case(2), 10, scop, seer)
+    load._limit_to_max_heat_pump_power = True
+    load.baseload_heating = test_load
+    load.peak_heating = np.full(12, 10)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period, test_load_sim_per * 5 / 6)
+    assert np.allclose(load._monthly_baseload_extraction_dhw_simulation_period, np.zeros(120))
+    assert np.allclose(load._monthly_baseload_extraction_heating_simulation_period, test_load_sim_per * 5 / 6)
+    load.cop = cop_basic
+    load.reset_results(1, 11)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period, test_load_sim_per * 1 / 2)
+    load.reset_results(10, 110)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period, test_load_sim_per * 19 / 20)
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 9, 9, 9, 9, 9, 9]), 10))
+    load.cop = cop_pl
+    load.reset_results(1, 11)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 5., 5., 5., 5., 5., 5.]), 10))
+    load.reset_results(10, 110)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       np.tile(np.array([4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5]), 10))
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       np.tile(np.array(
+                           [2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 9, 9, 9, 9, 9, 9]), 10))
+
+    # now with only DHW
+    load.cop_dhw = scop
+    load.baseload_heating = np.zeros(12)
+    load.dhw = test_load
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period, test_load_sim_per * 5 / 6)
+    assert np.allclose(load._monthly_baseload_extraction_dhw_simulation_period, test_load_sim_per * 5 / 6)
+    assert np.allclose(load._monthly_baseload_extraction_heating_simulation_period, np.zeros(120))
+    load.cop_dhw = cop_basic
+    load.reset_results(1, 11)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period, test_load_sim_per * 1 / 2)
+    load.reset_results(10, 110)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period, test_load_sim_per * 19 / 20)
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 9, 9, 9, 9, 9, 9]), 10))
+    # this will follow the data for the baseload, since the part load data is the same, because peak_heating is not
+    # relevant for the dhw power
+    load.cop_dhw = cop_pl
+    load.reset_results(1, 11)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 5., 5., 5., 5., 5., 5.]), 10))
+    load.reset_results(10, 110)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       np.tile(np.array([4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5]), 10))
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_baseload_extraction_simulation_period,
+                       np.tile(np.array(
+                           [2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 9, 9, 9, 9, 9, 9]), 10))
+
 
 def test_monthly_peak_injection_simulation_period():
     load = MonthlyBuildingLoadAbsolute(*load_case(2), 10, scop, seer)
@@ -437,6 +522,33 @@ def test_monthly_peak_injection_simulation_period():
                        np.tile(np.array(
                            [7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5]), 10))
 
+    load = MonthlyBuildingLoadAbsolute(*load_case(2), 10, scop, seer)
+    load._limit_to_max_heat_pump_power = True
+    load.peak_cooling = test_load
+    load.baseload_cooling = np.full(12, 10)
+    assert np.allclose(load.monthly_peak_injection_simulation_period, test_load_sim_per * 6 / 5)
+    load.eer = eer_basic
+    load.reset_results(0, 1)
+    assert np.allclose(load.monthly_peak_injection_simulation_period, test_load_sim_per * 3 / 2)
+    load.reset_results(0, 10)
+    assert np.allclose(load.monthly_peak_injection_simulation_period, test_load_sim_per * 21 / 20)
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       np.tile(np.array([7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 11, 11, 11, 11, 11, 11]), 10))
+    load.eer = eer_pl_2
+    load.reset_results(0, 1)
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       np.tile(np.array([7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, ]), 10))
+    load.reset_results(0, 10)
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       np.tile(
+                           np.array([5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25]),
+                           10))
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_peak_injection_simulation_period,
+                       np.tile(np.array(
+                           [7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5]), 10))
+
 
 def test_monthly_peak_extraction_simulation_period():
     load = MonthlyBuildingLoadAbsolute(*load_case(2), 10, scop, seer)
@@ -491,6 +603,59 @@ def test_monthly_peak_extraction_simulation_period():
     assert np.allclose(load.monthly_peak_extraction_simulation_period,
                        np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5]), 10))
 
+    load = MonthlyBuildingLoadAbsolute(*load_case(2), 10, scop, seer)
+    load._limit_to_max_heat_pump_power = True
+    load.peak_heating = test_load
+    load.baseload_heating = np.full(12, 10)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period, test_load_sim_per * 5 / 6)
+    assert np.allclose(load._monthly_peak_extraction_dhw_simulation_period, np.zeros(120))
+    assert np.allclose(load._monthly_peak_extraction_heating_simulation_period, test_load_sim_per * 5 / 6)
+    load.cop = cop_basic
+    load.reset_results(1, 11)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period, test_load_sim_per * 1 / 2)
+    load.reset_results(10, 110)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period, test_load_sim_per * 19 / 20)
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 9, 9, 9, 9, 9, 9]), 10))
+    load.cop = cop_pl_2
+    load.reset_results(1, 11)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5]), 10))
+    load.reset_results(10, 110)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       np.tile(np.array([4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75]), 10))
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5]), 10))
+
+    # now with only DHW
+    load.cop_dhw = scop
+    load.peak_heating = np.zeros(12)
+    load.baseload_heating = np.zeros(12)
+    load.dhw = test_load * 730
+    assert np.allclose(load.monthly_peak_extraction_simulation_period, test_load_sim_per * 5 / 6)
+    assert np.allclose(load._monthly_peak_extraction_dhw_simulation_period, test_load_sim_per * 5 / 6)
+    assert np.allclose(load._monthly_peak_extraction_heating_simulation_period, np.zeros(120))
+    load.cop_dhw = cop_basic
+    load.reset_results(1, 11)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period, test_load_sim_per * 1 / 2)
+    load.reset_results(10, 110)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period, test_load_sim_per * 19 / 20)
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 9, 9, 9, 9, 9, 9]), 10))
+    load.cop_dhw = cop_pl_2
+    load.reset_results(1, 11)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5]), 10))
+    load.reset_results(10, 110)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       np.tile(np.array([4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75]), 10))
+    load.set_results(results_monthly_test)
+    assert np.allclose(load.monthly_peak_extraction_simulation_period,
+                       np.tile(np.array([2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5]), 10))
+
 
 def test_max_loads():
     load = MonthlyBuildingLoadAbsolute(*load_case(2), 10, scop, seer)
diff --git a/GHEtool/utils/calculate_load.py b/GHEtool/utils/calculate_load.py
index 4da7107a..ef4b84d2 100644
--- a/GHEtool/utils/calculate_load.py
+++ b/GHEtool/utils/calculate_load.py
@@ -98,7 +98,7 @@ def xor(a, b):
         TMY: pd.DataFrame = pd.read_csv(weather_file, sep=",", header=None, skiprows=8)
 
         # drop first 6 columns
-        TMY.drop(columns=TMY.columns[:5], axis=1, inplace=True)
+        TMY.drop(columns=TMY.columns[:5], inplace=True)
 
         # set dry bulb temperature
         temperature: np.ndarray = np.array(TMY.iloc[:, 1])
@@ -110,7 +110,7 @@ def xor(a, b):
             TMY = pd.read_csv(weather_file, header=None, skiprows=18).iloc[:-10, :]
 
             # drop first 1 column
-            TMY.drop(columns=TMY.columns[:1], axis=1, inplace=True)
+            TMY.drop(columns=TMY.columns[:1], inplace=True)
 
             # set dry bulb temperature
             temperature: np.ndarray = TMY.iloc[:, 0].astype(float).to_numpy()
diff --git a/setup.cfg b/setup.cfg
index 2d100343..2e810864 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = GHEtool
-version = 2.4.1.dev3
+version = 2.4.1.dev4
 author = Wouter Peere
 author_email = wouter@ghetool.eu
 description = Python package for borefield sizing