From 9c0623b31d6d7703c13bee4fab07f3ec0076f81c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Pedro=20L=2E=20Magalh=C3=A3es?=
 <pedro.magalhaes@uni-bremen.de>
Date: Fri, 8 Mar 2024 02:46:21 +0100
Subject: [PATCH] Comment edits.

---
 src/topupopt/data/buildings/dk/bbr.py    |  76 +++++++-------
 src/topupopt/data/buildings/dk/heat.py   |  51 +++-------
 src/topupopt/data/dhn/network.py         |   4 +-
 src/topupopt/data/finance/invest.py      |  52 +++++-----
 src/topupopt/data/gis/identify.py        |   6 +-
 src/topupopt/data/gis/modify.py          |  22 ++--
 src/topupopt/data/gis/osm.py             |   8 +-
 src/topupopt/data/misc/units.py          |  12 +--
 src/topupopt/data/misc/utils.py          |  28 +++---
 src/topupopt/problems/esipp/converter.py |  80 +++++++--------
 src/topupopt/problems/esipp/dynsys.py    | 116 ++++++++++-----------
 src/topupopt/problems/esipp/signal.py    | 122 +++++++++++------------
 src/topupopt/problems/esipp/system.py    |  68 ++++++-------
 src/topupopt/problems/esipp/utils.py     |  98 +++++++++---------
 src/topupopt/solvers/interface.py        |  72 ++++++-------
 15 files changed, 398 insertions(+), 417 deletions(-)

diff --git a/src/topupopt/data/buildings/dk/bbr.py b/src/topupopt/data/buildings/dk/bbr.py
index ca1aa03..d0d2fc1 100644
--- a/src/topupopt/data/buildings/dk/bbr.py
+++ b/src/topupopt/data/buildings/dk/bbr.py
@@ -1,5 +1,5 @@
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # standard
 
@@ -13,8 +13,8 @@ from geopandas import GeoDataFrame, GeoSeries
 from pandas import DataFrame, merge
 from shapely.geometry import Point
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # URLs
     
@@ -34,8 +34,8 @@ url_prefix_buildings = 'https://api.dataforsyningen.dk/bbrlight/bygninger/'
 
 url_prefix_building_point = 'https://api.dataforsyningen.dk/bbrlight/bygningspunkter/'
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # labels
 
@@ -228,8 +228,8 @@ BBR_CONTAINER_LABELS = {
     'koordinater': list,
     }
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # BBR codes
 
@@ -511,8 +511,8 @@ list_labels_bbr = [
     label_bbr_extra_heating
     ]
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def get_bbr_building_data_geodataframe(
         building_entrance_ids: list,
@@ -521,8 +521,8 @@ def get_bbr_building_data_geodataframe(
         selected_bbr_building_point_labels: list = SELECT_BBR_BDG_POINT_LABELS
         ) -> Tuple[GeoDataFrame,list]:
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # get data about building entrances
     
@@ -548,8 +548,8 @@ def get_bbr_building_data_geodataframe(
         index=dict_building_entrances.keys()
         )
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # get data about buildings
     
@@ -577,8 +577,8 @@ def get_bbr_building_data_geodataframe(
         index=dict_buildings.keys()
         )
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # get building point data
       
@@ -626,8 +626,8 @@ def get_bbr_building_data_geodataframe(
                          left_index=True,
                          suffixes=(None,"_y")) # adds "_y" to duplicate columns
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # create a geodataframe whose geometry is that of building points
     
@@ -671,17 +671,17 @@ def get_bbr_building_data_geodataframe(
     
     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
     
@@ -733,21 +733,21 @@ def fetch_building_entrance_data(building_entrance_ids: list) -> Tuple[dict,
             
             response.close()
             
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     return dict_building_entrances, list_failures
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def fetch_building_data(building_codes: list):
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # get data about each specific building
     
@@ -781,13 +781,13 @@ def fetch_building_data(building_codes: list):
             
             response.close()
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     return dict_buildings
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
 
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/data/buildings/dk/heat.py b/src/topupopt/data/buildings/dk/heat.py
index 47941b5..3b35b05 100644
--- a/src/topupopt/data/buildings/dk/heat.py
+++ b/src/topupopt/data/buildings/dk/heat.py
@@ -1,18 +1,14 @@
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 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 .bbr import label_bbr_entrance_id, label_bbr_housing_area
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # labels
 
@@ -55,11 +51,10 @@ selected_bbr_building_labels = [
 
 
 # label under which building entrance ids can be found in OSM
-
 label_osm_entrance_id = 'osak:identifier'
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def heat_demand_dict_by_building_entrance(
         gdf_osm: GeoDataFrame, 
@@ -67,7 +62,7 @@ def heat_demand_dict_by_building_entrance(
         number_intervals: int,
         time_interval_durations: list,
         bdg_specific_demand: dict,
-        bdg_base_load_avg_ratio: dict,
+        bdg_ratio_min_max: dict,
         bdg_demand_phase_shift: dict = None,
         key_osm_entr_id: str = label_osm_entrance_id,
         key_bbr_entr_id: str = label_bbr_entrance_id,
@@ -119,7 +114,7 @@ def heat_demand_dict_by_building_entrance(
                         discrete_sinusoid_matching_integral(
                             bdg_specific_demand[building_index]*area,
                             time_interval_durations=time_interval_durations, 
-                            ratio_min_avg=bdg_base_load_avg_ratio[building_index],
+                            bdg_ratio_min_max=bdg_ratio_min_max[building_index],
                             phase_shift_radians=(
                                 bdg_demand_phase_shift[building_index]
                                 # bdg_demand_phase_shift_amplitude*np.random.random() 
@@ -142,39 +137,33 @@ def heat_demand_dict_by_building_entrance(
                                 ), 
                             avg_state=avg_state, 
                             time_interval_durations=time_interval_durations, 
-                            ratio_min_avg=bdg_base_load_avg_ratio[building_index],
+                            bdg_ratio_min_max=bdg_ratio_min_max[building_index],
                             state_correlates_with_output=state_correlates_with_output
                             )
                         )
                     )
             
-            #******************************************************************
+            # *****************************************************************
             
         # 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
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def total_heating_area(
         gdf_osm: GeoDataFrame, 
@@ -184,27 +173,19 @@ def total_heating_area(
         ) -> float:
     
     area = 0
-    
     for osm_index in gdf_osm.index:
-
         # find the indexes for each building leading to the curr. cons. point
-
         building_indexes = (
             gdf_buildings[
                 gdf_buildings[label_bbr_entrance_id] ==
                 gdf_osm.loc[osm_index][label_osm_entrance_id]
                 ].index
             )
-
         # for each building
-
         for building_index in building_indexes:
-
             # get relevant data
-
             area += gdf_buildings.loc[building_index][label_bbr_housing_area]
-            
     return area
 
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/data/dhn/network.py b/src/topupopt/data/dhn/network.py
index a8839ce..7e0d4d1 100644
--- a/src/topupopt/data/dhn/network.py
+++ b/src/topupopt/data/dhn/network.py
@@ -16,8 +16,8 @@ from ...problems.esipp.network import ArcsWithoutProportionalLosses
 
 from ...data.finance.utils import ArcInvestments
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # scenarios: different datasets/trench objects
 # options: fully specified trenches, including prices
diff --git a/src/topupopt/data/finance/invest.py b/src/topupopt/data/finance/invest.py
index ca682b2..da6f6d0 100644
--- a/src/topupopt/data/finance/invest.py
+++ b/src/topupopt/data/finance/invest.py
@@ -1,12 +1,12 @@
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 from math import prod
 
 from statistics import mean
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # TODO: enable swapping the polarity
 
@@ -129,8 +129,8 @@ class Investment:
             for i in range(self.analysis_period_span+1)
             )
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def add_investment(self, 
                        investment: float,
@@ -169,8 +169,8 @@ class Investment:
             self.net_cash_flows[investment_period] += investment
             self.net_cash_flows[self.analysis_period_span] += -residual_value
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def add_operational_cash_flows(self,
                                    cash_flow: float or int,
@@ -202,19 +202,19 @@ class Investment:
                 
                 self.net_cash_flows[i+start_period] += cash_flow
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def net_present_value(self):
         """Returns the net present value for the investment under analysis."""
         
         return npv(self.discount_rates, self.net_cash_flows)
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def npv(discount_rates: list,
         net_cash_flows: list,
@@ -272,8 +272,8 @@ def npv(discount_rates: list,
             for (ncf_t, df_t) in zip(net_cash_flows, discount_factors)
             )
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
     
 def discount_factor(discount_rates: list) -> float:
     """
@@ -298,8 +298,8 @@ def discount_factor(discount_rates: list) -> float:
     """    
     return prod([1/(1+i) for i in discount_rates])
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def salvage_value_linear_depreciation(
         investment: int or float, 
@@ -355,8 +355,8 @@ def salvage_value_linear_depreciation(
         analysis_period_span
         )*investment/investment_longevity
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def salvage_value_annuity(investment: int or float, 
                           discount_rate: float,
@@ -377,8 +377,8 @@ def salvage_value_annuity(investment: int or float,
         tuple(discount_rate for i in range(analysis_period_span))
         )
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def annuity(investment: int or float, 
             investment_longevity: int,
@@ -392,8 +392,8 @@ def annuity(investment: int or float,
             ))
         )
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def present_salvage_value_annuity(investment: int or float, 
                                   investment_longevity: int,
@@ -497,5 +497,5 @@ def present_salvage_value_annuity(investment: int or float,
             net_cash_flows=net_cash_flows
             )
 
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/data/gis/identify.py b/src/topupopt/data/gis/identify.py
index 33ccd46..97beec8 100644
--- a/src/topupopt/data/gis/identify.py
+++ b/src/topupopt/data/gis/identify.py
@@ -1694,13 +1694,13 @@ def identify_edge_closest_to_node(
 
     """
     
-    #**************************************************************************
+    # *************************************************************************
     
     # 1) ensure that the network crs is correct and convert if not    
     # 2) identify the edges that are nearest to the nodes
     # 3) revert network crs back to the original, if necessary
 
-    #**************************************************************************
+    # *************************************************************************
             
     # if it is a geographic CRS, convert it to a projected CRS
     
@@ -1714,7 +1714,7 @@ def identify_edge_closest_to_node(
         
         projected_network = network
 
-    #**************************************************************************
+    # *************************************************************************
     
     # 2) identify the edges that are nearest to the nodes
     
diff --git a/src/topupopt/data/gis/modify.py b/src/topupopt/data/gis/modify.py
index 52bd2b4..e8cad22 100644
--- a/src/topupopt/data/gis/modify.py
+++ b/src/topupopt/data/gis/modify.py
@@ -483,7 +483,7 @@ def replace_path(
             )
         }
     
-    #**************************************************************************
+    # *************************************************************************
     
     # add edges
     start_node = path[0]  
@@ -827,7 +827,7 @@ def split_linestring(line: LineString,
     # return the geometries for each segment and the relevant points by order
     return line_segments, close_to_start, close_to_end
      
-    #**************************************************************************
+    # *************************************************************************
 
 # *****************************************************************************
 # *****************************************************************************
@@ -1100,7 +1100,7 @@ def connect_nodes_to_edges(
 
     """
         
-    #**************************************************************************
+    # *************************************************************************
     
     # 1) group nodes by the edge they are closest to
     # 2) for each edge, and node that is to be connected to it, find its closest 
@@ -1110,8 +1110,8 @@ def connect_nodes_to_edges(
     # 5) delete the original edges, if they have been split
     # 6) calculate or update the edge attributes
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # 1) group nodes by the edge they are closest to
         
@@ -1124,8 +1124,8 @@ def connect_nodes_to_edges(
         for edge_key in set(edge_keys)
         }
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # 2) for each edge, and node that is to be connected to it, find its closest 
     # point on the edge
@@ -1161,8 +1161,8 @@ def connect_nodes_to_edges(
         # TIP: exclude the points that can be considered close to the start or end nodes
     # TIP: use the shortest line method to obtain the line geometry
             
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # 3) recreate each edge after dividing it at the specified points
     
@@ -1184,8 +1184,8 @@ def connect_nodes_to_edges(
         
     network.remove_edges_from(recreated_edges)
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # 4) connect the nodes to the edges
     
diff --git a/src/topupopt/data/gis/osm.py b/src/topupopt/data/gis/osm.py
index e784105..efc83bf 100644
--- a/src/topupopt/data/gis/osm.py
+++ b/src/topupopt/data/gis/osm.py
@@ -1,5 +1,5 @@
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # OpenStreetMaps
 
@@ -36,7 +36,7 @@ KEY_OSM_BUILDING_ENTRANCE_ID = {
     }
  
 
-#******************************************************************************
+# *****************************************************************************
 
 # osmnx
 
@@ -101,4 +101,4 @@ KEYS_OSMNX_EDGES_ESSENTIAL = {
     KEY_OSMNX_REVERSED
     }
 
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/data/misc/units.py b/src/topupopt/data/misc/units.py
index a60a6bb..a19fb1f 100644
--- a/src/topupopt/data/misc/units.py
+++ b/src/topupopt/data/misc/units.py
@@ -1,7 +1,7 @@
 # constants
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # unit conversion factors
 
@@ -15,12 +15,12 @@ GJ_DIV_MWh = 1000/3600
 
 MWh_DIV_J = 1/(3600*1000*1000)
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # currency conversions
 
 EUR_DIV_DKK = 1/(743.95/100)
     
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/data/misc/utils.py b/src/topupopt/data/misc/utils.py
index 9a70e63..745a410 100644
--- a/src/topupopt/data/misc/utils.py
+++ b/src/topupopt/data/misc/utils.py
@@ -1,5 +1,5 @@
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # standard
 
@@ -11,8 +11,8 @@ from statistics import mean
 
 # local, external
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
                     
 def generate_pseudo_unique_key(key_list: tuple, 
                                max_iterations: int = 10) -> str:
@@ -38,8 +38,8 @@ def generate_pseudo_unique_key(key_list: tuple,
         ' iterations.'
         )
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def discrete_sinusoid_matching_integral(
         integration_result: float,
@@ -124,8 +124,8 @@ def discrete_sinusoid_matching_integral(
         for i in range(number_time_steps)
         ]
          
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def synch_profile(
         profile: list, 
@@ -199,8 +199,8 @@ def synch_profile(
             for i, r in enumerate(reference_profile)
             ]
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def create_profile_using_time_weighted_state(
         integration_result: float,
@@ -307,8 +307,8 @@ def create_profile_using_time_weighted_state(
     n = len(time_interval_durations)
     return [*new_profile[n-_sorted[0][1]:], *new_profile[0:n-_sorted[0][1]]]
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def max_min_sinusoidal_profile(
         integration_result: float or int,
@@ -367,5 +367,5 @@ def max_min_sinusoidal_profile(
         b*time_interval_duration-amplitude
         )
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
diff --git a/src/topupopt/problems/esipp/converter.py b/src/topupopt/problems/esipp/converter.py
index d3f4435..769179a 100644
--- a/src/topupopt/problems/esipp/converter.py
+++ b/src/topupopt/problems/esipp/converter.py
@@ -17,8 +17,8 @@ from .dynsys import DynamicSystem
 
 from .signal import Signal, FixedSignal
           
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # objectives:
 # 1) organise information for optimisation
@@ -60,13 +60,13 @@ class Converter:
                  state_specific_externality_costs: dict = None # one for every output, then another for every time interval
                  ):
                   
-        #**********************************************************************
+        # *********************************************************************
         
         self.key = key
         
         self.sys = sys
           
-        #**********************************************************************
+        # *********************************************************************
         
         # inputs
         
@@ -117,7 +117,7 @@ class Converter:
         
         self.fixed_inputs = self.identify_fixed_signals(self.inputs)
             
-        #**********************************************************************
+        # *********************************************************************
         
         # amplitude costs: one per signal
         
@@ -141,8 +141,8 @@ class Converter:
             output_specific_externality_costs
             )
         
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # identify the signals with specific amplitude costs
         
@@ -182,8 +182,8 @@ class Converter:
                 self.output_specific_amplitude_costs.keys()
                 )
             
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # identify the signals with specific externality costs
         
@@ -229,7 +229,7 @@ class Converter:
                 if value != 0
                 ]
         
-        #**********************************************************************
+        # *********************************************************************
         
         # identify binary signals
         
@@ -249,11 +249,11 @@ class Converter:
         
         # self.identify_amplitude_penalised_inputs()
     
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
-    #**************************************************************************
-    #**************************************************************************  
+    # *************************************************************************
+    # *************************************************************************  
     
     @staticmethod
     def identify_fixed_signals(signals: list):
@@ -262,8 +262,8 @@ class Converter:
             u for u, sig in enumerate(signals) if isinstance(sig, FixedSignal)
             ]
     
-    #**************************************************************************
-    #**************************************************************************  
+    # *************************************************************************
+    # *************************************************************************  
     
     # TODO: number of time intervals in the dynamic system class
     
@@ -352,8 +352,8 @@ class Converter:
         
         return a_innk, b_inmk, c_irnk, d_irmk, e_x_ink, e_y_irk
     
-# #**************************************************************************
-# #**************************************************************************  
+# # *************************************************************************
+# # *************************************************************************  
         
     # def has_dimensionable_inputs(self):
         
@@ -367,8 +367,8 @@ class Converter:
             
     #         return True
             
-    # #**************************************************************************
-    # #**************************************************************************  
+    # # *************************************************************************
+    # # *************************************************************************  
         
     # def has_binary_inputs(self):
         
@@ -382,8 +382,8 @@ class Converter:
             
     #         return True
             
-    # #**************************************************************************
-    # #**************************************************************************  
+    # # *************************************************************************
+    # # *************************************************************************  
         
     # def has_amplitude_penalised_inputs(self):
         
@@ -397,8 +397,8 @@ class Converter:
             
     #         return True
             
-    # #**************************************************************************
-    # #**************************************************************************  
+    # # *************************************************************************
+    # # *************************************************************************  
         
     # def has_externality_inducing_inputs(self):
         
@@ -412,8 +412,8 @@ class Converter:
             
     #         return True
             
-    # #**************************************************************************
-    # #**************************************************************************  
+    # # *************************************************************************
+    # # *************************************************************************  
         
     # def has_externality_inducing_outputs(self):
         
@@ -427,8 +427,8 @@ class Converter:
             
     #         return True
             
-    # #**************************************************************************
-    # #**************************************************************************  
+    # # *************************************************************************
+    # # *************************************************************************  
     
     # def identify_dimensionable_inputs(self):
         
@@ -437,8 +437,8 @@ class Converter:
     #         for i, u in enumerate(self.inputs)
     #         if u.is_dimensionable]
     
-    # #**************************************************************************
-    # #**************************************************************************  
+    # # *************************************************************************
+    # # *************************************************************************  
     
     # def identify_binary_inputs(self):
         
@@ -447,8 +447,8 @@ class Converter:
     #         for i, u in enumerate(self.inputs)
     #         if u.is_binary]
                 
-    # #**************************************************************************
-    # #**************************************************************************
+    # # *************************************************************************
+    # # *************************************************************************
     
     # def identify_externality_inducing_inputs(self):
         
@@ -457,8 +457,8 @@ class Converter:
     #         for i, c in enumerate(self.input_externalities)
     #         if c != 0]
                 
-    # #**************************************************************************
-    # #**************************************************************************
+    # # *************************************************************************
+    # # *************************************************************************
     
     # def identify_externality_inducing_outputs(self):
         
@@ -467,8 +467,8 @@ class Converter:
     #         for i, c in enumerate(self.output_externalities)
     #         if c != 0]
                 
-    # #**************************************************************************
-    # #**************************************************************************
+    # # *************************************************************************
+    # # *************************************************************************
     
     # def identify_amplitude_penalised_inputs(self):
         
@@ -477,8 +477,8 @@ class Converter:
     #         for i, c in enumerate(self.input_amplitude_costs)
     #         if c != 0]
             
-    # #**************************************************************************
-    # #**************************************************************************
+    # # *************************************************************************
+    # # *************************************************************************
     
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/dynsys.py b/src/topupopt/problems/esipp/dynsys.py
index d461f5c..082298b 100644
--- a/src/topupopt/problems/esipp/dynsys.py
+++ b/src/topupopt/problems/esipp/dynsys.py
@@ -15,8 +15,8 @@ from scipy.linalg import expm, inv
 
 # local libraries, internal
                          
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # TODO: integrate only selected outputs on demand, not all of them
 
@@ -52,8 +52,8 @@ class DynamicSystem:
     # TODO: define whether this class holds information about the number of 
     # time intervals or if that information comes from somewhere else
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
              
     def __init__(self, 
                  time_interval_durations: list or float or int,
@@ -68,8 +68,8 @@ class DynamicSystem:
         # 2) compute and/or adjust the data in preparation for discretisation
         # 3) discretise the model
                  
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # cases:
         # 1) time invariant system, regular time steps
@@ -82,8 +82,8 @@ class DynamicSystem:
         # >> multiple continous SS models, each discretised once
         # note: case 3 and 4 are effectively the same
         
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # if time_interval_durations is a number, transform it into a list
                 
@@ -137,7 +137,7 @@ class DynamicSystem:
                  C=C,
                  D=D)
                  
-        #**********************************************************************
+        # *********************************************************************
         
         # determine the number of models
         
@@ -214,8 +214,8 @@ class DynamicSystem:
         
         # determine the number of models
              
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
              
         # get the data into lists
         
@@ -253,18 +253,18 @@ class DynamicSystem:
             
             self.D = [D]
                 
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # discretise the model
             
         self.discretise(integrate_outputs=integrate_outputs)
             
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
             
-    #**************************************************************************
-    #************************************************************************** 
+    # *************************************************************************
+    # ************************************************************************* 
 
     def validate(self,
                  time_interval_durations: list or np.array,
@@ -273,8 +273,8 @@ class DynamicSystem:
                  C: list or np.array,
                  D: list or np.array) -> tuple:
         
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # time interval duration
         
@@ -319,8 +319,8 @@ class DynamicSystem:
             raise TypeError(
                 'Unrecognised time interval duration format.')
         
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
                 
         # if multiple matrices are provided for a given position, then they
         # must all have the same shape (rows and columns)
@@ -542,7 +542,7 @@ class DynamicSystem:
            
             raise TypeError('Unrecognised D matrix format.')
         
-        #**********************************************************************
+        # *********************************************************************
             
         # make sure the matrices have appropriate shapes
         
@@ -591,7 +591,7 @@ class DynamicSystem:
                     'The B and D matrices do not have the same number of columns.'
                     )
             
-        #**********************************************************************
+        # *********************************************************************
         
         # number of inputs, states and outputs
         
@@ -601,7 +601,7 @@ class DynamicSystem:
         
         # number_outputs = D_m # C_m or D_m (better, since C is optional)
                 
-        #**********************************************************************
+        # *********************************************************************
             
         return (
             B_n if D_n == 0 else D_n,
@@ -613,27 +613,27 @@ class DynamicSystem:
             D_is_time_varying
             )
             
-    #**************************************************************************
-    #************************************************************************** 
+    # *************************************************************************
+    # ************************************************************************* 
 
     def discretise(self, integrate_outputs: bool = True):
         
         # cases:
-        #**********************************************************************
+        # *********************************************************************
         # 1) time invariant system, regular time steps
         # >> one continuous SS model discretised once and reused every step
-        #**********************************************************************
+        # *********************************************************************
         # 2) time invariant system, irregular time steps
         # >> one continuous SS model discretised once for each time step
-        #**********************************************************************
+        # *********************************************************************
         # 3) time-varying system, regular time steps
         # >> multiple continous SS models, each discretised once
-        #**********************************************************************
+        # *********************************************************************
         # 4) time-varying system, irregular time steps
         # >> multiple continous SS models, each discretised once
-        #**********************************************************************
+        # *********************************************************************
         # note: case 3 and 4 are effectively the same
-        #**********************************************************************
+        # *********************************************************************
         
         # TODO: do not assume that A, B, C, D lists have the same lengths if time-varying
         
@@ -684,7 +684,7 @@ class DynamicSystem:
                 expm(self.time_interval_durations[0]*A_matrix)
                 for A_matrix in self.A]
             
-        #**********************************************************************
+        # *********************************************************************
     
         # B matrix
         
@@ -746,7 +746,7 @@ class DynamicSystem:
                     )
                 for (A_matrix, B_matrix) in zip(self.A, self.B)]
             
-        #**********************************************************************
+        # *********************************************************************
     
         # C matrix
         
@@ -840,7 +840,7 @@ class DynamicSystem:
                 
                 self.C_line_k = [C_matrix for C_matrix in self.C]
                 
-        #**********************************************************************
+        # *********************************************************************
         
         # D matrix
         
@@ -898,8 +898,8 @@ class DynamicSystem:
                     self.D[i if self.D_is_time_varying else 0]
                     for i, dt in enumerate(self.time_interval_durations)]
                 
-            #******************************************************************
-            #******************************************************************
+            # *****************************************************************
+            # *****************************************************************
             
         elif len(self.D) == 1: 
             
@@ -955,8 +955,8 @@ class DynamicSystem:
                     self.D[i if self.D_is_time_varying else 0]
                     for i, dt in enumerate(self.time_interval_durations)]
                 
-            #******************************************************************
-            #******************************************************************
+            # *****************************************************************
+            # *****************************************************************
             
         else: # the number of time steps and models do not match
         
@@ -1014,11 +1014,11 @@ class DynamicSystem:
                 
                 self.D_line_k = [D_matrix for D_matrix in self.D]
                 
-            #******************************************************************
-            #******************************************************************
+            # *****************************************************************
+            # *****************************************************************
             
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # store information about the matrices
         
@@ -1068,8 +1068,8 @@ class DynamicSystem:
             
             self.outputs_integrated = []
         
-    #**************************************************************************
-    #**************************************************************************  
+    # *************************************************************************
+    # *************************************************************************  
     
     def simulate(self, U: np.ndarray, X0: np.ndarray):
         
@@ -1113,7 +1113,7 @@ class DynamicSystem:
         #m = self.number_inputs
         r = self.number_outputs
         
-        #**********************************************************************
+        # *********************************************************************
                     
         if (self.A_line_is_time_varying is not None and 
             self.B_line_is_time_varying is not None):
@@ -1259,11 +1259,11 @@ class DynamicSystem:
                 
                 return None, Y[:,1:]
     
-    #**************************************************************************
-    #**************************************************************************  
+    # *************************************************************************
+    # *************************************************************************  
         
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
         
 class StatelessSystem(DynamicSystem):
     "A class for dynamic systems without states."
@@ -1284,16 +1284,16 @@ class StatelessSystem(DynamicSystem):
             D=D,
             integrate_outputs=integrate_outputs) # outputs can be integrated
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def simulate(self, U: np.ndarray):
         "Simulate how the system responds to a set of input signals."
         
         return DynamicSystem.simulate(self, U=U, X0=None)
         
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
         
 class OutputlessSystem(DynamicSystem):
     "A class for dynamic systems without outputs."
@@ -1312,13 +1312,13 @@ class OutputlessSystem(DynamicSystem):
             D=None,
             integrate_outputs=False) # no outputs to integrate
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def simulate(self, U: np.ndarray, X0: np.ndarray):
         "Simulate how the system responds to a set of input signals and initial conditions."
         
         return DynamicSystem.simulate(self, U=U, X0=X0)
 
-#******************************************************************************
-#******************************************************************************   
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************   
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/signal.py b/src/topupopt/problems/esipp/signal.py
index 8029595..e134ea4 100644
--- a/src/topupopt/problems/esipp/signal.py
+++ b/src/topupopt/problems/esipp/signal.py
@@ -11,8 +11,8 @@ Created on Wed Nov 17 13:04:53 2021
 
 # local libraries, internal
                  
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # TODO: consider mentioning that the classes are actually for discretised systems and signals
     
@@ -49,8 +49,8 @@ class Signal:
         
         self.is_bounded = self.is_upper_bounded or self.is_lower_bounded
             
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def validate_inputs(self,
                         number_samples: int,
@@ -152,8 +152,8 @@ class Signal:
 
         return is_fixed, has_upper_bounds, has_lower_bounds             
             
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def has_upper_bounds(self) -> bool:
         "Returns True if the signal has predetermined upper bounds."
@@ -178,8 +178,8 @@ class Signal:
             
         return False
             
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def has_lower_bounds(self) -> bool:
         "Returns True if the signal has predetermined lower bounds."
@@ -204,8 +204,8 @@ class Signal:
             
         return False
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def is_signal_bounded(self) -> bool:
         "Returns True if the signal has predetermined lower or upper bounds."
@@ -214,8 +214,8 @@ class Signal:
         
         return self.has_lower_bounds() or self.has_upper_bounds()
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def is_signal_fixed(self) -> bool:
         "Returns True if the signal is predetermined or has been fixed."
@@ -243,8 +243,8 @@ class Signal:
             
             return False
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def violates_bounds(self, tolerance: float = 0.0) -> bool:
         "Returns True if the signal violates its bounds or if it has none."
@@ -285,8 +285,8 @@ class Signal:
             
             return False
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def set_signal(self, samples: list):
         "Defines a signal using externally-provided samples."
@@ -309,8 +309,8 @@ class Signal:
         
         self.is_fixed = True
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def is_signal_integer_only(self, 
                                integrality_tolerance: float = 0.0) -> bool:
@@ -331,8 +331,8 @@ class Signal:
             
             return False
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def is_signal_binary_only(self, 
                               integrality_tolerance: float = 0.0) -> bool:
@@ -353,11 +353,11 @@ class Signal:
             
             return False
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def are_elements_integer(elements: list, 
                          integrality_tolerance: float or int = 0.0) -> bool:
@@ -391,8 +391,8 @@ def are_elements_integer(elements: list,
     
     return True
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def are_elements_binary(elements: list, 
                         integrality_tolerance: float or int = None) -> bool:
@@ -418,7 +418,7 @@ def are_elements_binary(elements: list,
             
         return True
     
-    #**************************************************************************
+    # *************************************************************************
     
     # a specific integrality tolerance was provided
     
@@ -460,8 +460,8 @@ def are_elements_binary(elements: list,
         
     return True
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class FixedSignal(Signal):
     "A class for signals that are predetermined, measured or decided on."
@@ -481,8 +481,8 @@ class FixedSignal(Signal):
                         lower_bounds=lower_bounds,
                         upper_bounds=upper_bounds)
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class NonNegativeRealSignal(Signal):
     "A class for non-negative real signals."
@@ -525,8 +525,8 @@ class NonNegativeRealSignal(Signal):
             
             raise ValueError('The bounds are not non-negative real.')
                 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
         
     def is_nnr(self, tolerance: float = 0.0) -> bool:
         "Returns True if the samples are all non-negative real."
@@ -551,8 +551,8 @@ class NonNegativeRealSignal(Signal):
             
             return False
             
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def are_bounds_nnr(self, tolerance: float = 0.0) -> bool:
         "Returns True if the bounds are non-negative real."
@@ -598,11 +598,11 @@ class NonNegativeRealSignal(Signal):
             
             return False
        
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class FixedNonNegativeRealSignal(NonNegativeRealSignal):
     "A class for free and controllable signals."
@@ -627,8 +627,8 @@ class FixedNonNegativeRealSignal(NonNegativeRealSignal):
             
             raise ValueError('At least one sample is not non-negative real.')
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class FreeSignal(Signal):
     "A class for undetermined signals."
@@ -644,8 +644,8 @@ class FreeSignal(Signal):
                         lower_bounds=lower_bounds,
                         upper_bounds=upper_bounds)
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class FreeUnboundedSignal(FreeSignal):
     "A class for free and unbounded signals."
@@ -658,8 +658,8 @@ class FreeUnboundedSignal(FreeSignal):
             lower_bounds=None,
             upper_bounds=None)
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class BinarySignal(NonNegativeRealSignal):
     "A class for binary signals."
@@ -677,8 +677,8 @@ class BinarySignal(NonNegativeRealSignal):
             lower_bounds=[0 for _ in range(number_samples)],
             upper_bounds=[1 for _ in range(number_samples)])
         
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class AmplitudeConstrainedSignal(Signal):
     "A class for amplitude-constrained signals."
@@ -731,7 +731,7 @@ class AmplitudeConstrainedSignal(Signal):
         
         self.set_negative_amplitude(negative_amplitude)
             
-    #**************************************************************************
+    # *************************************************************************
 
     def violates_amplitude_limits(
             self, 
@@ -821,7 +821,7 @@ class AmplitudeConstrainedSignal(Signal):
             
             return False
         
-    #**************************************************************************
+    # *************************************************************************
     
     def set_positive_amplitude(self, positive_amplitude: float or int = None):
         """
@@ -853,7 +853,7 @@ class AmplitudeConstrainedSignal(Signal):
             
             raise TypeError('Unknown positive amplitude type.')
                 
-    #**************************************************************************
+    # *************************************************************************
     
     def set_negative_amplitude(self, negative_amplitude: float or int = None):
         """
@@ -885,7 +885,7 @@ class AmplitudeConstrainedSignal(Signal):
             
             raise TypeError('Unknown negative amplitude type.')
         
-    #**************************************************************************
+    # *************************************************************************
     
     def set_positive_amplitude_limits(self, 
                                       max_pos_amp_limit: float or int,
@@ -939,7 +939,7 @@ class AmplitudeConstrainedSignal(Signal):
         self.has_pos_amp_limits = (
             self.has_max_pos_amp_limit or self.has_min_pos_amp_limit)
 
-    #**************************************************************************
+    # *************************************************************************
     
     def set_negative_amplitude_limits(self, 
                                       max_neg_amp_limit: float or int,
@@ -994,7 +994,7 @@ class AmplitudeConstrainedSignal(Signal):
         self.has_neg_amp_limits = (
             self.has_max_neg_amp_limit or self.has_min_neg_amp_limit)
 
-    #**************************************************************************
+    # *************************************************************************
     
     @staticmethod
     def validate_limits(max_amp_limit: float or int,
@@ -1079,7 +1079,7 @@ class AmplitudeConstrainedSignal(Signal):
                 
         return has_max_amp, has_min_amp
 
-    #**************************************************************************
+    # *************************************************************************
     
     @staticmethod
     def validate_positive_bounds(has_max_pos_amp_limit: bool,
@@ -1125,7 +1125,7 @@ class AmplitudeConstrainedSignal(Signal):
                     'The positive amplitude limits are not compatible with '+
                     ' the bounds provided.')
 
-    #**************************************************************************
+    # *************************************************************************
     
     @staticmethod
     def validate_negative_bounds(has_max_neg_amp_limit: bool,
@@ -1173,7 +1173,7 @@ class AmplitudeConstrainedSignal(Signal):
                     'The negative amplitude limits are not compatible with '+
                     ' the bounds provided.')
 
-    #**************************************************************************
+    # *************************************************************************
     
     def validate_positive_amplitude(self, 
                                     positive_amplitude: float or int = None,
@@ -1254,7 +1254,7 @@ class AmplitudeConstrainedSignal(Signal):
                             'The positive amplitude is below its tolerated '+
                             'minimum.')
 
-    #**************************************************************************
+    # *************************************************************************
     
     def validate_negative_amplitude(self, 
                                     negative_amplitude: float or int = None,
@@ -1335,8 +1335,8 @@ class AmplitudeConstrainedSignal(Signal):
                             'The negative amplitude is below its tolerated '+
                             'minimum.')
             
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class AmplitudeConstrainedNNRSignal(AmplitudeConstrainedSignal,
                                     NonNegativeRealSignal):
@@ -1401,5 +1401,5 @@ class AmplitudeConstrainedNNRSignal(AmplitudeConstrainedSignal,
             
             raise ValueError('The bounds are not non-negative real.')
             
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/system.py b/src/topupopt/problems/esipp/system.py
index d737ee1..fdcb636 100644
--- a/src/topupopt/problems/esipp/system.py
+++ b/src/topupopt/problems/esipp/system.py
@@ -17,8 +17,8 @@ from .converter import Converter
 
 from .network import Network
 
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 class EnergySystem:
     """A class to model energy systems for simulation or optimisation."""
@@ -40,8 +40,8 @@ class EnergySystem:
     # technologies selected
     
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
             
     def __init__(self, 
                  networks: dict = None, 
@@ -49,8 +49,8 @@ class EnergySystem:
                  optional_converters: list = None,
                  selected_converters: list = None):
                  
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # networks should be a dict of nx.MultiDiGraph objects
         
@@ -72,8 +72,8 @@ class EnergySystem:
             
             self.converters = {}
             
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # optional converters
         
@@ -89,8 +89,8 @@ class EnergySystem:
         
             self.optional_converters = []
             
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # selected converters
         
@@ -106,8 +106,8 @@ class EnergySystem:
         
             self.selected_converters = []
             
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         # TODO: explain why some of attributes should stay here
         
@@ -133,19 +133,19 @@ class EnergySystem:
         
         self.network_effects = {}
         
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def add_network(self, network_key, network: Network):
         """Add a new network to the energy system."""
         
         self.networks[network_key] = network
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def add_converter(self, 
                       converter_key, 
@@ -176,8 +176,8 @@ class EnergySystem:
                 
         return True
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def set_converter_as_optional_or_not(self,
                                          converter_key,
@@ -198,8 +198,8 @@ class EnergySystem:
                 
                 self.optional_converters.remove(converter_key)
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def set_converter_selection_status(self, 
                                        converter_key, 
@@ -220,8 +220,8 @@ class EnergySystem:
                 
                 self.selected_converters.remove(converter_key)
                 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def set_converter_network_effects(self,
                                       converter_key,
@@ -269,8 +269,8 @@ class EnergySystem:
                 
         raise NotImplementedError
             
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def make_arc_mandatory(self, network_key, arc_key: tuple):
         """
@@ -323,8 +323,8 @@ class EnergySystem:
                 'Either the network key provided is incorrect or the arc key '+
                 'lacks the proper size.')
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def unmake_arc_mandatory(self, network_key, arc_key: tuple):
         """
@@ -366,8 +366,8 @@ class EnergySystem:
                 'Either the network key provided is incorrect or the arc key '+
                 'lacks the proper size.')
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def set_maximum_number_parallel_arcs(self, 
                                          network_key, 
@@ -384,10 +384,10 @@ class EnergySystem:
                 self.max_number_parallel_arcs[
                     (network_key,node_a, node_b)] = limit
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # TODO: create a method to identify mandatory nodes
             
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/utils.py b/src/topupopt/problems/esipp/utils.py
index d0dc4ff..f191307 100644
--- a/src/topupopt/problems/esipp/utils.py
+++ b/src/topupopt/problems/esipp/utils.py
@@ -1,5 +1,5 @@
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # standard
 
@@ -74,13 +74,13 @@ def run_mvesipp_analysis(problem: InfrastructurePlanningProblem = None,
                          analyse_results: bool = False,
                          analyse_problem: bool = False):
     
-    #**************************************************************************
+    # *************************************************************************
     
     if model_instance != None and analyse_problem:
         
         describe_mves(model_instance)
 
-    #**************************************************************************
+    # *************************************************************************
 
     if model_instance != None and analyse_results:
                
@@ -96,7 +96,7 @@ def run_mvesipp_analysis(problem: InfrastructurePlanningProblem = None,
             flow_in_cost,
             flow_out_revenue)
 
-    #**************************************************************************
+    # *************************************************************************
     
     if problem != None and analyse_results:
         
@@ -104,7 +104,7 @@ def run_mvesipp_analysis(problem: InfrastructurePlanningProblem = None,
         
         describe_solution(problem)
 
-    #**************************************************************************
+    # *************************************************************************
     
 # *****************************************************************************
 # *****************************************************************************
@@ -286,16 +286,16 @@ def compute_cost_volume_metrics(instance: pyo.ConcreteModel,
             for q, p in instance.set_QP
             }
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     return flow_in, flow_in_k, flow_out, flow_in_cost, flow_out_revenue
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def compute_network_performance(solved_problem: InfrastructurePlanningProblem):
     
@@ -342,8 +342,8 @@ def compute_network_performance(solved_problem: InfrastructurePlanningProblem):
         flow_out_revenue
         )
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # provide a summary of the results
 
@@ -354,8 +354,8 @@ def present_summary_results(flow_in: dict,
                             flow_unit: str = 'MWh',
                             currency: str = 'EUR'):
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
         
     if len(flow_in) != 0:
 
@@ -385,8 +385,8 @@ def present_summary_results(flow_in: dict,
             print(
                 'Average price: N/A (no flow imports are set to take place).')
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
         
     if len(flow_out) != 0:
     
@@ -416,11 +416,11 @@ def present_summary_results(flow_in: dict,
             print(
                 'Average price: N/A (no flow exports are set to take place).')
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
         
 def unused_node_key(network: nx.MultiDiGraph):
     """Returns an unused node key."""
@@ -432,8 +432,8 @@ def unused_node_key(network: nx.MultiDiGraph):
             # it doesn't, return it
             return i
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # TODO: document
 
@@ -520,8 +520,8 @@ def compute_gross_network_flows(problem: InfrastructurePlanningProblem) -> dict:
         'gross_demand_g': gross_demand_g,
         }
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def describe_mves(obj: pyo.ConcreteModel):
     
@@ -540,7 +540,7 @@ def describe_mves(obj: pyo.ConcreteModel):
     
     for g in obj.set_G:
     
-        #**********************************************************************
+        # *********************************************************************
         
         # the number of nodes for each network
         
@@ -595,7 +595,7 @@ def describe_mves(obj: pyo.ConcreteModel):
         print('Min. and max. flow needs for other nodes: '+
               str(static_flow_limits))
         
-        #**********************************************************************
+        # *********************************************************************
         
         # arcs
         
@@ -615,9 +615,9 @@ def describe_mves(obj: pyo.ConcreteModel):
             
         # unreacheable nodes
         
-        #**********************************************************************
+        # *********************************************************************
     
-    #**************************************************************************
+    # *************************************************************************
     
     # systems
         
@@ -733,14 +733,14 @@ def plot_networks(ipp: InfrastructurePlanningProblem,
     # plt.axis('off')
     # plt.show()
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     plt.rcParams["figure.figsize"] = [7.50, 3.50]
     plt.rcParams["figure.autolayout"] = True
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     node_shapes = ['o','s','d']
     
@@ -837,7 +837,7 @@ def plot_networks(ipp: InfrastructurePlanningProblem,
         edge_alphas = [(5 + i) / (M + 4) for i in range(M)]
         cmap = plt.cm.plasma
         
-        #**********************************************************************
+        # *********************************************************************
         
         # plot import nodes
         
@@ -958,11 +958,11 @@ def plot_networks(ipp: InfrastructurePlanningProblem,
     # nx.draw(G, pos)
     # plt.show()
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 def is_integer(variable: float, integrality_tolerance: float) -> bool:
     """Returns True if a given number qualifies as an integer."""
@@ -973,12 +973,12 @@ def is_integer(variable: float, integrality_tolerance: float) -> bool:
             )
     return not (abs(round(variable)-variable) > integrality_tolerance)
     
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
     
 def describe_solution(ipp: InfrastructurePlanningProblem):
         
-    #**************************************************************************
+    # *************************************************************************
     
     print('******************************************************************')
         
@@ -1137,9 +1137,9 @@ def describe_solution(ipp: InfrastructurePlanningProblem):
                                 #print(tech_options_sec)
                                 #print('******************')
                     
-            #******************************************************************
+            # *****************************************************************
                     
-        #**********************************************************************
+        # *********************************************************************
         
         # for each node
         
@@ -1271,13 +1271,13 @@ def describe_solution(ipp: InfrastructurePlanningProblem):
                                 #print(tech_options_sec)
                                 #print('******************')
                     
-            #******************************************************************
+            # *****************************************************************
         
-        #**********************************************************************
+        # *********************************************************************
                     
     print('******************************************************************')
             
-    #**************************************************************************
+    # *************************************************************************
 
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/solvers/interface.py b/src/topupopt/solvers/interface.py
index d8178e1..279c858 100644
--- a/src/topupopt/solvers/interface.py
+++ b/src/topupopt/solvers/interface.py
@@ -1,5 +1,5 @@
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
 # key features:
 # 1) same name and units for the same attributes across different solvers
@@ -34,8 +34,8 @@ from pyomo.opt.results.solver import check_optimal_termination
 class SolverInterface(object):
     "A class for interfacing with solvers through Pyomo."
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
 
     SOLVER_GLPK = 'glpk'
     
@@ -56,8 +56,8 @@ class SolverInterface(object):
                SOLVER_CPLEX,
                SOLVER_CPLEX_DIRECT]
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # topupopt options
     
@@ -103,8 +103,8 @@ class SolverInterface(object):
     
     PYOMO_OPTIONS = ['executable','warmstart','solver_io','tee']
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     PROBLEM_LP = 'lp'
     PROBLEM_MILP = 'milp'
@@ -168,8 +168,8 @@ class SolverInterface(object):
         SOLVER_CPLEX_DIRECT: PROBLEMS_COMPATIBLE_CPLEX,
         }
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     # USER -->> SOLVER_INTERFACE -->> PYOMO --> SOLVER
     
@@ -214,8 +214,8 @@ class SolverInterface(object):
         TerminationCondition.licensingProblems
         ]
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
          
     def __init__(self, 
                  solver_name: str, 
@@ -244,8 +244,8 @@ class SolverInterface(object):
         
         self.obj = SolverFactory(solver_name, **kwargs)
         
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def get_solver_handler(self, **kwargs):
             
@@ -276,8 +276,8 @@ class SolverInterface(object):
                 options_dict_format=options_dict_format,
                 options_param_format=options_param_format)
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     @staticmethod
     def get_pyomo_solver_object(solver_name: str, 
@@ -302,8 +302,8 @@ class SolverInterface(object):
         
         return opt
     
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     @classmethod
     def problem_and_solver_are_compatible(cls,
@@ -326,8 +326,8 @@ class SolverInterface(object):
             
             return False
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     def was_optimisation_sucessful(self, results, problem_type) -> bool:
                 
@@ -355,8 +355,8 @@ class SolverInterface(object):
             
             raise self.UnknownSolverStatusError(solver_status)
             
-        #**********************************************************************
-        #**********************************************************************
+        # *********************************************************************
+        # *********************************************************************
         
         expected_solver_status = TerminationCondition.to_solver_status(
             termination_condition)
@@ -381,8 +381,8 @@ class SolverInterface(object):
             
             return False
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     class InconsistentProblemTypeAndSolverError(Exception):
         
@@ -393,8 +393,8 @@ class SolverInterface(object):
                 ') cannot be handled by solver ('+solver_name+').'
                 )
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
 
     class InconsistentSolverStatusError(Exception):
         
@@ -405,8 +405,8 @@ class SolverInterface(object):
                 ') does not match the one obtained ('+solver_status+').'
                 )
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     class UnknownSolverError(Exception):
         
@@ -414,8 +414,8 @@ class SolverInterface(object):
             
             super().__init__('Unknown solver: '+str(solver_name))
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     class UnknownProblemTypeError(Exception):
         
@@ -423,8 +423,8 @@ class SolverInterface(object):
             
             super().__init__('Unknown problem type: '+str(problem_type))
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     class UnknownTerminationConditionError(Exception):
         
@@ -433,8 +433,8 @@ class SolverInterface(object):
             super().__init__(
                 'Unknown termination condition: '+termination_condition)
 
-    #**************************************************************************
-    #**************************************************************************
+    # *************************************************************************
+    # *************************************************************************
     
     class UnknownSolverStatusError(Exception):
         
@@ -442,5 +442,5 @@ class SolverInterface(object):
             
             super().__init__('Unknown solver status: '+solver_status)
     
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
-- 
GitLab