Created new methods for problem statistics.

src/topupopt/data/buildings/dk/heat.py

@@ -4,7 +4,9 @@
 import numpy as np
 from geopandas import GeoDataFrame
 from ...misc.utils import discrete_sinusoid_matching_integral
-from ...misc.utils import create_profile_using_time_weighted_state
+# from ...misc.utils import create_profile_using_time_weighted_state
+from ...misc.utils import generate_manual_state_correlated_profile
+from ...misc.utils import generate_state_correlated_profile
 from .bbr import label_bbr_entrance_id, label_bbr_housing_area
 
 # *****************************************************************************
@@ -66,6 +68,8 @@ def heat_demand_dict_by_building_entrance(
     key_bbr_entr_id: str = label_bbr_entrance_id,
     avg_state: list = None,
     state_correlates_with_output: bool = False,
+    deviation_gain: float = 1.0,
+    solver: str = 'glpk'
 ) -> dict:
     
     # initialise dict for each building entrance
@@ -84,11 +88,8 @@ def heat_demand_dict_by_building_entrance(
         # for each building
         for building_index in building_indexes:
             # get relevant data
-            # base_load_avg_ratio = 0.3
-            # specific_demand = 107 # kWh/m2/year
             area = gdf_buildings.loc[building_index][label_bbr_housing_area]
-
-            # estimate its demand
+            # generate the profile
             if type(avg_state) == type(None):
                 # ignore states
                 heat_demand_profiles.append(
@@ -104,18 +105,46 @@ def heat_demand_dict_by_building_entrance(
                     )
                 )
 
+            elif type(deviation_gain) == type(None):
+                # states matter but the gain must be determined
+                # heat_demand_profiles.append(
+                #     np.array(
+                #         create_profile_using_time_weighted_state(
+                #             integration_result=(
+                #                 bdg_specific_demand[building_index] * area
+                #             ),
+                #             avg_state=avg_state,
+                #             time_interval_durations=time_interval_durations,
+                #             min_to_max_ratio=bdg_ratio_min_max[building_index],
+                #             state_correlates_with_output=state_correlates_with_output,
+                #         )
+                #     )
+                # )
+                heat_demand_profiles.append(
+                    np.array(
+                        generate_state_correlated_profile(
+                            integration_result=(
+                                bdg_specific_demand[building_index] * area
+                            ), 
+                            states=avg_state, 
+                            time_interval_durations=time_interval_durations,
+                            states_correlate_profile=state_correlates_with_output,
+                            min_max_ratio=bdg_ratio_min_max[building_index],
+                            solver=solver
+                            )
+                        )
+                    )
             else:
-                # states matter
+                # states matter and the gain is predefined
                 heat_demand_profiles.append(
                     np.array(
-                        create_profile_using_time_weighted_state(
+                        generate_manual_state_correlated_profile(
                             integration_result=(
                                 bdg_specific_demand[building_index] * area
                             ), 
-                            avg_state=avg_state,
+                            states=avg_state, 
                             time_interval_durations=time_interval_durations, 
-                            min_to_max_ratio=bdg_ratio_min_max[building_index],
-                            state_correlates_with_output=state_correlates_with_output,
+                            deviation_gain=deviation_gain
                             )
                         )
                     )

src/topupopt/data/misc/utils.py

@@ -368,7 +368,7 @@ def generate_state_correlated_profile(
     states_correlate_profile: bool,
     min_max_ratio: float,
     solver: str
-) -> list:
+) -> tuple:
     """
     Returns a profile observing a number of conditions.
     
@@ -379,6 +379,23 @@ def generate_state_correlated_profile(
     
     This method relies on linear programming. An LP solver must be used.
     
+    The profile for interval i is defined as follows:
+
+    P[i] = (dt[i]/dt_mean)*( (x[i]-x_mean)*gain + offset)
+
+    where:
+
+    dt[i] is the time interval duration for interval i
+    
+    dt_mean is the mean time interval duration
+    
+    x[i] is the state for interval i
+    
+    x_mean is the mean state for the entire profile
+    
+    The offset is defined as the integration result divided by the number 
+    time intervals, whereas the gain is determined via optimisation.
+
     Parameters
     ----------
     integration_result : float
@@ -402,14 +419,16 @@ def generate_state_correlated_profile(
 
     Returns
     -------
-    list
-        A profile matching the aforementioned characteristics.
+    tuple
+        A tuple containing the profile, the gain and the offset.
 
     """    
         
     # *************************************************************************
     # *************************************************************************
     
+    # TODO: find alternative solution, as this is most likely overkill
+    
     # internal model
     
     def model(states_correlate_profile: bool) -> pyo.AbstractModel:
@@ -477,10 +496,8 @@ def generate_state_correlated_profile(
     # *************************************************************************
     # *************************************************************************
     
-    
     dt_total = sum(time_interval_durations)
     dt_mean = mean(time_interval_durations)
-    # x_mean  = mean(states)
     x_mean = sum(
         deltat_k*x_k 
         for deltat_k, x_k in zip(time_interval_durations, states)
@@ -510,7 +527,11 @@ def generate_state_correlated_profile(
     problem = fit_model.create_instance(data=data_dict) 
     opt.solve(problem, tee=False)
     # return profile
-    return [pyo.value(problem.P_i[i]) for i in problem.I]
+    return (
+        [pyo.value(problem.P_i[i]) for i in problem.I], 
+        pyo.value(problem.alpha),
+        beta
+        )
 
 # *****************************************************************************
 # *****************************************************************************

tests/test_data_utils.py

 # imports
 
 # standard
-
 import random
-
 import math
-
 from statistics import mean
 
 # local, internal
-
 from src.topupopt.data.misc import utils
 
-# ******************************************************************************
-# ******************************************************************************
-
-
 class TestDataUtils:
-    # *************************************************************************
-    # *************************************************************************
 
     def test_profile_synching2(self):
         integration_result = 10446
@@ -390,7 +380,7 @@ class TestDataUtils:
         states_correlate_profile = True
         min_max_ratio = 0.2
         
-        profile = utils.generate_state_correlated_profile(
+        profile, a, b = utils.generate_state_correlated_profile(
             integration_result=integration_result, 
             states=states, 
             time_interval_durations=time_interval_durations, 
@@ -400,6 +390,7 @@ class TestDataUtils:
             )
         
         # test profile 
+        assert a > 0 and b > 0
         assert len(profile) == number_time_intervals
         assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)
         assert math.isclose(min(profile), max(profile)*min_max_ratio, abs_tol=1e-3)
@@ -413,7 +404,7 @@ class TestDataUtils:
         states_correlate_profile = False
         min_max_ratio = 0.2
         
-        profile = utils.generate_state_correlated_profile(
+        profile, a, b = utils.generate_state_correlated_profile(
             integration_result=integration_result, 
             states=states, 
             time_interval_durations=time_interval_durations, 
@@ -423,6 +414,7 @@ class TestDataUtils:
             )
         
         # test profile 
+        assert a < 0 and b > 0
         assert len(profile) == number_time_intervals
         assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)
         assert math.isclose(min(profile), max(profile)*min_max_ratio, abs_tol=1e-3)
@@ -464,7 +456,6 @@ class TestDataUtils:
         # correlation: direct, inverse
         # states: positive, negative
         # time intervals: regular irregular
-        # 
         
         # profile with positive correlation, positive states, regular intervals
         number_time_intervals = 10