diff --git a/src/topupopt/data/buildings/dk/bbr.py b/src/topupopt/data/buildings/dk/bbr.py
index f5c1185e4fab8520bfa5665b115c08fb646bffd8..9f41972b8e21e041af48789be89d403a60f7eebd 100644
--- a/src/topupopt/data/buildings/dk/bbr.py
+++ b/src/topupopt/data/buildings/dk/bbr.py
@@ -518,25 +518,19 @@ list_labels_bbr = [
 # *****************************************************************************
 # *****************************************************************************
 
-
 def get_bbr_building_data_geodataframe(
     building_entrance_ids: list,
     selected_bbr_bdg_entrance_labels: list = SELECT_BBR_BDG_ENTR_LABELS,
     selected_bbr_building_labels: list = SELECT_BBR_BDG_LABELS,
     selected_bbr_building_point_labels: list = SELECT_BBR_BDG_POINT_LABELS,
 ) -> Tuple[GeoDataFrame, list]:
-    # *************************************************************************
-    # *************************************************************************
 
     # get data about building entrances
-
     dict_building_entrances, list_failures = fetch_building_entrance_data(
         building_entrance_ids
     )
-
     if selected_bbr_bdg_entrance_labels == None:
         # includes all labels
-
         selected_bbr_bdg_entrance_labels = BBR_BDG_ENTR_LABELS
 
     list_entries = [
@@ -547,7 +541,6 @@ def get_bbr_building_data_geodataframe(
         ]
         for key, value in dict_building_entrances.items()
     ]
-
     df_building_entrances = DataFrame(
         data=list_entries,
         columns=selected_bbr_bdg_entrance_labels,
@@ -560,19 +553,15 @@ def get_bbr_building_data_geodataframe(
     # get data about buildings
 
     dict_buildings = fetch_building_data(df_building_entrances.index)
-
     if selected_bbr_building_labels == None:
         # includes all labels
-
         selected_bbr_building_labels = BBR_BDG_LABELS
 
     # create dataframe with building data
-
     list_entries = [
         [value[bbr_key] for bbr_key in value if bbr_key in selected_bbr_building_labels]
         for key, value in dict_buildings.items()
     ]
-
     df_buildings = DataFrame(
         data=list_entries,
         columns=selected_bbr_building_labels,
@@ -583,12 +572,10 @@ def get_bbr_building_data_geodataframe(
     # *************************************************************************
 
     # get building point data
-
     if selected_bbr_building_point_labels == None:
         # includes all labels
 
         selected_bbr_building_point_labels = BBR_BDG_POINT_LABELS
-
     dict_buildings_points = {
         building_entrance_id: dict_buildings[building_entrance_id][  # (
             label_bbr_bygningpunkt
@@ -599,7 +586,6 @@ def get_bbr_building_data_geodataframe(
     }
 
     # create dataframe with building point data
-
     list_entries = [
         [
             value[bbr_key]
@@ -608,7 +594,6 @@ def get_bbr_building_data_geodataframe(
         ]
         for key, value in dict_buildings_points.items()
     ]
-
     df_building_points = DataFrame(
         data=list_entries,
         columns=selected_bbr_building_point_labels,
@@ -637,15 +622,11 @@ def get_bbr_building_data_geodataframe(
     # *************************************************************************
 
     # create a geodataframe whose geometry is that of building points
-
     # specify the coordinate system
-
     coordinate_system = "EPSG:4326"  # latitude, longitude
-
     key_bbr_coordinate_system = 5  # WGS 84 = EPSG:4326
 
     # raise an error if different coordinates systems are being used
-
     for building_entrance_id in dict_building_entrances:
         if (
             dict_buildings[building_entrance_id][label_bbr_bygningpunkt][
@@ -656,7 +637,6 @@ def get_bbr_building_data_geodataframe(
             raise NotImplementedError("Only WGS 84 coordinates can be used.")
 
     # create a dictionary with the building point geometries (i.e. points)
-
     dict_building_point_geometry = {
         building_entrance_id: Point(
             dict_buildings[building_entrance_id][label_bbr_bygningpunkt][
@@ -671,121 +651,70 @@ def get_bbr_building_data_geodataframe(
     }
 
     # create geodataframe
-
     gdf_buildings = GeoDataFrame(
         data=df_buildings,
         geometry=GeoSeries(data=dict_building_point_geometry),
         crs=coordinate_system,
     )
-
+    
     return gdf_buildings, list_failures
 
-    # *************************************************************************
-    # *************************************************************************
-
-
 # *****************************************************************************
 # *****************************************************************************
 
-
 def fetch_building_entrance_data(building_entrance_ids: list) -> Tuple[dict, list]:
-    # *************************************************************************
-    # *************************************************************************
 
     # retrieve data about each node identified through OSM
-
     dict_building_entrances = {}
-
     list_failures = []
-
     # # determine the number of osm entries
-
     # number_building_entrance_ids = len(building_entrance_ids)
-
     # for each building entrance id
-
     for building_entrance_id in building_entrance_ids:
         # compose the url from which to get bbr data associated with the id
-
         _url = url_prefix_entrance + building_entrance_id
-
         try:
             # retrieve the building entrance data
-
             with urllib.request.urlopen(_url) as response:
                 # parse the data
-
                 bbr_entrance_data_json = json.loads(response.read().decode("utf-8"))
-
                 # store the data
-
                 if len(bbr_entrance_data_json) != 0:
                     for bbr_entry in bbr_entrance_data_json:
                         dict_building_entrances[
                             bbr_entry[label_bbr_building_id]
                         ] = bbr_entry
-
                 else:
                     list_failures.append(building_entrance_id)
-
             response.close()
-
         except Exception:
             response.close()
 
-    # *************************************************************************
-    # *************************************************************************
-
     return dict_building_entrances, list_failures
 
-    # *************************************************************************
-    # *************************************************************************
-
-
 # *****************************************************************************
 # *****************************************************************************
 
-
 def fetch_building_data(building_codes: list):
-    # *************************************************************************
-    # *************************************************************************
 
     # get data about each specific building
-
     dict_buildings = {}
-
     # for each building id
-
     for building_id in building_codes:
         # compose a url with it
-
         _url = url_prefix_buildings + building_id
-
         # try statement
-
         try:
             # retrieve that data
-
             with urllib.request.urlopen(_url) as response:
                 # parse the data
-
                 bbr_data = json.loads(response.read().decode("utf-8"))
-
                 dict_buildings[building_id] = bbr_data
-
             response.close()
-
         except Exception:
             response.close()
 
-    # *************************************************************************
-    # *************************************************************************
-
     return dict_buildings
 
-    # *************************************************************************
-    # *************************************************************************
-
-
 # *****************************************************************************
 # *****************************************************************************
diff --git a/src/topupopt/data/buildings/dk/heat.py b/src/topupopt/data/buildings/dk/heat.py
index 66d60eea9004bd3a371a9c0719fe39af818d1b98..76bb32275c25d73446a0d3c6787f304bd26584a9 100644
--- a/src/topupopt/data/buildings/dk/heat.py
+++ b/src/topupopt/data/buildings/dk/heat.py
@@ -54,7 +54,6 @@ label_osm_entrance_id = "osak:identifier"
 # *****************************************************************************
 # *****************************************************************************
 
-
 def heat_demand_dict_by_building_entrance(
     gdf_osm: GeoDataFrame,
     gdf_buildings: GeoDataFrame,
@@ -68,39 +67,30 @@ def heat_demand_dict_by_building_entrance(
     avg_state: list = None,
     state_correlates_with_output: bool = False,
 ) -> dict:
+    
     # initialise dict for each building entrance
-
     demand_dict = {}
 
     # for each building entrance
-
     for osm_index in gdf_osm.index:
         # initialise dict for each building consumption point
-
         heat_demand_profiles = []
 
         # find the indexes for each building leading to the curr. cons. point
-
         building_indexes = gdf_buildings[
             gdf_buildings[key_bbr_entr_id] == gdf_osm.loc[osm_index][key_osm_entr_id]
         ].index
 
         # 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
-
             if type(avg_state) == type(None):
                 # ignore states
-
                 heat_demand_profiles.append(
                     np.array(
                         discrete_sinusoid_matching_integral(
@@ -109,9 +99,6 @@ def heat_demand_dict_by_building_entrance(
                             min_to_max_ratio=bdg_ratio_min_max[building_index],
                             phase_shift_radians=(
                                 bdg_demand_phase_shift[building_index]
-                                # bdg_demand_phase_shift_amplitude*np.random.random()
-                                # if (type(bdg_demand_phase_shift_amplitude) ==
-                                #     type(None)) else None
                             ),
                         )
                     )
@@ -119,7 +106,6 @@ def heat_demand_dict_by_building_entrance(
 
             else:
                 # states matter
-
                 heat_demand_profiles.append(
                     np.array(
                         create_profile_using_time_weighted_state(
@@ -134,21 +120,15 @@ def heat_demand_dict_by_building_entrance(
                     )
                 )
 
-            # *****************************************************************
-
         # add the profiles, time step by time step
         if len(heat_demand_profiles) == 0:
             final_profile = []
         else:
             final_profile = sum(profile for profile in heat_demand_profiles)
 
-        # *********************************************************************
-
         # store the demand profile
         demand_dict[osm_index] = final_profile
 
-        # *********************************************************************
-
     # return
     return demand_dict
 
@@ -178,4 +158,4 @@ def total_heating_area(
 
 
 # *****************************************************************************
-# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/model.py b/src/topupopt/problems/esipp/model.py
index b19f63b27cccf8fd0a699baac15fd16f05a2f0c8..1b37e1640920faa213eb7b1e178c9e492fc15d1e 100644
--- a/src/topupopt/problems/esipp/model.py
+++ b/src/topupopt/problems/esipp/model.py
@@ -14,7 +14,7 @@ def create_model(
     enable_validation: bool = True,
     enable_initialisation: bool = True,
 ):
-    # TODO: make default values, validation, and initialisation optional
+    # TODO: test initialisation
 
     # *************************************************************************
     # *************************************************************************
@@ -41,8 +41,6 @@ def create_model(
 
     # TODO: try to use blocks to improve readability and modularity
 
-    # TODO: add piecewise constraints for the import and export flows/prices
-
     # *************************************************************************
     # *************************************************************************
     # *************************************************************************