diff --git a/src/mesido/base_component_type_mixin.py b/src/mesido/base_component_type_mixin.py index c5d2f7c30..9b77c483e 100644 --- a/src/mesido/base_component_type_mixin.py +++ b/src/mesido/base_component_type_mixin.py @@ -20,6 +20,14 @@ def energy_system_components(self) -> Dict[str, str]: """ raise NotImplementedError + @property + @abstractmethod + def energy_system_owners(self) -> Dict[str, str]: + """ + This method return a dict with the components structured by owner. + """ + raise NotImplementedError + def energy_system_components_get(self, list_types: list) -> list: components = [] for component_type in list_types: diff --git a/src/mesido/component_type_mixin.py b/src/mesido/component_type_mixin.py index 8c0d2a3ec..783be8df7 100644 --- a/src/mesido/component_type_mixin.py +++ b/src/mesido/component_type_mixin.py @@ -23,6 +23,7 @@ class ModelicaComponentTypeMixin(BaseComponentTypeMixin): def __init__(self, *args, **kwargs): super().__init__(**kwargs) self.__hn_component_types = None + self.__owners = None def pre(self): """ @@ -347,7 +348,13 @@ def energy_system_components(self) -> Dict[str, Set[str]]: # Find the components in model, detection by string # (name.component_type: type) - component_types = sorted({v for k, v in string_parameters.items()}) + component_types = sorted( + { + v + for k, v in string_parameters.items() + if ("component_type" in k or "component_subtype" in k) + } + ) components = {} for c in component_types: @@ -359,6 +366,28 @@ def energy_system_components(self) -> Dict[str, Set[str]]: return self.__hn_component_types + @property + def energy_system_owners(self): + if self.__owners is None: + # Create the dictionary once after that it will be available + string_parameters = self.string_parameters(0) + + # Find the components in model, detection by string + # (name.component_type: type) + owners = sorted({v for k, v in string_parameters.items() if "owner" in k}) + + ownership_dict = {} + for c in owners: + ownership_dict[c] = sorted( + {k.split(".")[0] for k, v in string_parameters.items() if v == c} + ) + + ownership_dict.pop("__", None) + + self.__owners = ownership_dict + + return self.__owners + @property def energy_system_topology(self) -> Topology: """ diff --git a/src/mesido/esdl/esdl_heat_model.py b/src/mesido/esdl/esdl_heat_model.py index 6306fe60d..b26218327 100644 --- a/src/mesido/esdl/esdl_heat_model.py +++ b/src/mesido/esdl/esdl_heat_model.py @@ -101,6 +101,100 @@ def _rho_cp_modifiers(self) -> Dict: """ return dict(rho=self.rho, cp=self.cp) + def merge_modifiers(self, a: dict, b: dict): + """ + Recursive (not in place) merge of dictionaries. + + :param a: Base dictionary to merge. + :param b: Dictionary to merge on top of base dictionary. + :return: Merged dictionary + """ + b = b.copy() + + for k, v in a.items(): + if isinstance(v, dict): + b_node = b.setdefault(k, {}) + b[k] = self.merge_modifiers(v, b_node) + else: + if k not in b: + b[k] = v + + return b + + def get_owner(self, asset): + return dict( + owner=( + asset.attributes["isOwnedBy"].name + if asset.attributes["isOwnedBy"] is not None + else "NoOwner" + ) + ) + + def get_carrier_id(self, asset, n=0): + if ( + asset.in_ports is not None + and asset.out_ports is not None + and len(asset.in_ports) >= 2 + and len(asset.out_ports) >= 2 + ): # heat pump and heat exchanger + carrier = asset.global_properties["carriers"][asset.in_ports[0].carrier.id] + if "Prim" in asset.in_ports[0].name: + prim_in_id = carrier["id_number_mapping"] + else: + sec_in_id = carrier["id_number_mapping"] + carrier = asset.global_properties["carriers"][asset.in_ports[1].carrier.id] + if "Prim" in asset.in_ports[1].name: + prim_in_id = carrier["id_number_mapping"] + else: + sec_in_id = carrier["id_number_mapping"] + if len(asset.in_ports) == 3: + if "Prim" in asset.in_ports[1].name: + prim_in_id = carrier["id_number_mapping"] + else: + sec_in_id = carrier["id_number_mapping"] + carrier = asset.global_properties["carriers"][asset.out_ports[0].carrier.id] + if "Prim" in asset.out_ports[0].name: + prim_out_id = carrier["id_number_mapping"] + else: + sec_out_id = carrier["id_number_mapping"] + carrier = asset.global_properties["carriers"][asset.out_ports[1].carrier.id] + if "Prim" in asset.out_ports[1].name: + prim_out_id = carrier["id_number_mapping"] + else: + sec_out_id = carrier["id_number_mapping"] + ids = dict( + Primary=dict( + HeatIn=dict(carrier_id=prim_in_id), HeatOut=dict(carrier_id=prim_out_id) + ), + Secondary=dict( + HeatIn=dict(carrier_id=sec_in_id), HeatOut=dict(carrier_id=sec_out_id) + ), + ) + return ids + else: + ids = dict() + try: + carrier = asset.global_properties["carriers"][asset.in_ports[0].carrier.id] + id_number = carrier["id_number_mapping"] + port = f"{carrier['type'].capitalize()}In" + ids[port] = dict(carrier_id=id_number) + carrier = asset.global_properties["carriers"][asset.out_ports[0].carrier.id] + id_number = carrier["id_number_mapping"] + port = f"{carrier['type'].capitalize()}Out" + ids[port] = dict(carrier_id=id_number) + except KeyError: + try: + carrier = asset.global_properties["carriers"][asset.in_ports[0].carrier.id] + id_number = carrier["id_number_mapping"] + port = f"{carrier['type'].capitalize()}In" + ids[port] = dict(carrier_id=id_number) + except KeyError: + carrier = asset.global_properties["carriers"][asset.out_ports[0].carrier.id] + id_number = carrier["id_number_mapping"] + port = f"{carrier['type'].capitalize()}Out" + ids[port] = dict(carrier_id=id_number) + return ids + def get_asset_attribute_value( self, asset: Asset, @@ -309,7 +403,9 @@ def convert_heat_buffer(self, asset: Asset) -> Tuple[Type[HeatBuffer], MODIFIERS **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return HeatBuffer, modifiers @@ -371,7 +467,9 @@ def convert_heat_demand(self, asset: Asset) -> Tuple[Type[HeatDemand], MODIFIERS **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return HeatDemand, modifiers @@ -444,7 +542,9 @@ def convert_cold_demand(self, asset: Asset) -> Tuple[Type[ColdDemand], MODIFIERS **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return ColdDemand, modifiers @@ -507,9 +607,16 @@ def convert_node(self, asset: Asset) -> Tuple[Type[Node], MODIFIERS]: if isinstance(x, esdl.esdl.OutPort): sum_out += len(x.connectedTo) + try: + carrier = asset.global_properties["carriers"][asset.in_ports[0].carrier.id] + except KeyError: + carrier = asset.global_properties["carriers"][asset.out_ports[0].carrier.id] + modifiers = dict( n=sum_in + sum_out, state=self.get_state(asset), + carrier_id=carrier["id_number_mapping"], + **self.get_owner(asset), ) if isinstance(asset.in_ports[0].carrier, esdl.esdl.GasCommodity) or isinstance( @@ -613,7 +720,9 @@ def convert_heat_pipe( GasOut=bounds_nominals, state=self.get_state(asset), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return GasPipe, modifiers @@ -669,10 +778,13 @@ def convert_heat_pipe( **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) if "T_ground" in asset.attributes.keys(): modifiers["T_ground"] = asset.attributes["T_ground"] + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) + return HeatPipe, modifiers def convert_pump(self, asset: Asset) -> Tuple[Type[Pump], MODIFIERS]: @@ -736,7 +848,9 @@ def convert_pump(self, asset: Asset) -> Tuple[Type[Pump], MODIFIERS]: state=self.get_state(asset), **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return Pump, modifiers @@ -891,7 +1005,9 @@ def convert_heat_exchanger(self, asset: Asset) -> Tuple[Type[HeatExchanger], MOD state=self.get_state(asset), **self._get_cost_figure_modifiers(asset), **params, + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return HeatExchanger, modifiers def convert_heat_pump( @@ -1005,7 +1121,9 @@ def convert_heat_pump( state=self.get_state(asset), **self._get_cost_figure_modifiers(asset), **params, + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) if len(asset.in_ports) == 2: return HeatPump, modifiers elif len(asset.in_ports) == 3: @@ -1099,7 +1217,9 @@ def convert_heat_source(self, asset: Asset) -> Tuple[Type[HeatSource], MODIFIERS **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) if asset.asset_type == "GeothermalSource": modifiers["nr_of_doublets"] = asset.attributes["aggregationCount"] @@ -1235,7 +1355,9 @@ def convert_ates(self, asset: Asset) -> Tuple[Type[ATES], MODIFIERS]: **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) # if no maxStorageTemperature is specified we assume a "regular" HT ATES model if ( @@ -1339,7 +1461,9 @@ def convert_control_valve(self, asset: Asset) -> Tuple[Type[ControlValve], MODIF state=self.get_state(asset), **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return ControlValve, modifiers @@ -1394,7 +1518,9 @@ def convert_check_valve(self, asset: Asset) -> Tuple[Type[CheckValve], MODIFIERS HeatOut=dict(Hydraulic_power=dict(nominal=q_nominal * 16.0e5)), **self._supply_return_temperature_modifiers(asset), **self._rho_cp_modifiers, + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return CheckValve, modifiers @@ -1458,7 +1584,9 @@ def convert_electricity_demand(self, asset: Asset) -> Tuple[Type[ElectricityDema V=dict(min=min_voltage, nominal=min_voltage), ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return ElectricityDemand, modifiers @@ -1552,7 +1680,9 @@ def convert_electricity_source(self, asset: Asset) -> Tuple[Type[ElectricitySour Power=dict(min=0.0, max=max_supply, nominal=max_supply / 2.0), ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) if asset.asset_type in ["ElectricityProducer", "Import"]: return ElectricitySource, modifiers @@ -1603,7 +1733,9 @@ def convert_electricity_storage( min=-max_discharge, max=max_charge, nominal=max_charge / 2.0 ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return ElectricityStorage, modifiers @@ -1665,7 +1797,17 @@ def convert_electricity_node(self, asset: Asset) -> Tuple[Type[ElectricityNode], if isinstance(x, esdl.esdl.OutPort): sum_out += len(x.connectedTo) - modifiers = dict(voltage_nominal=nominal_voltage, n=sum_in + sum_out) + try: + carrier = asset.global_properties["carriers"][asset.in_ports[0].carrier.id] + except KeyError: + carrier = asset.global_properties["carriers"][asset.out_ports[0].carrier.id] + + modifiers = dict( + voltage_nominal=nominal_voltage, + carrier_id=carrier["id_number_mapping"], + n=sum_in + sum_out, + **self.get_owner(asset), + ) return ElectricityNode, modifiers @@ -1734,7 +1876,9 @@ def convert_electricity_cable(self, asset: Asset) -> Tuple[Type[ElectricityCable Power=dict(min=-max_power, max=max_power, nominal=max_power / 2.0), ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return ElectricityCable, modifiers def convert_transformer(self, asset: Asset) -> Tuple[Type[Transformer], MODIFIERS]: @@ -1848,7 +1992,9 @@ def convert_gas_demand(self, asset: Asset) -> Tuple[Type[GasDemand], MODIFIERS]: Hydraulic_power=dict(min=0.0, max=0.0, nominal=q_nominal * pressure), ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return GasDemand, modifiers @@ -1917,7 +2063,9 @@ def convert_gas_source(self, asset: Asset) -> Tuple[Type[GasSource], MODIFIERS]: Hydraulic_power=dict(nominal=q_nominal * pressure), ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return GasSource, modifiers @@ -2025,7 +2173,9 @@ def equations(x): V=dict(min=v_min, nominal=v_min), ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return Electrolyzer, modifiers @@ -2092,7 +2242,9 @@ def convert_gas_tank_storage(self, asset: Asset) -> Tuple[Type[GasTankStorage], Hydraulic_power=dict(nominal=q_nominal * pressure), ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return GasTankStorage, modifiers @@ -2158,7 +2310,9 @@ def convert_gas_substation(self, asset: Asset) -> Tuple[Type[GasSubstation], MOD Hydraulic_power=dict(nominal=q_nom_out * pressure_out), ), **self._get_cost_figure_modifiers(asset), + **self.get_owner(asset), ) + modifiers = self.merge_modifiers(modifiers, self.get_carrier_id(asset)) return GasSubstation, modifiers diff --git a/src/mesido/esdl/esdl_mixin.py b/src/mesido/esdl/esdl_mixin.py index f6caf85c0..b43a1b282 100644 --- a/src/mesido/esdl/esdl_mixin.py +++ b/src/mesido/esdl/esdl_mixin.py @@ -148,10 +148,15 @@ def __init__(self, *args, **kwargs) -> None: self._override_gas_pipe_classes = dict() self.override_gas_pipe_classes() + self.__connections = self.__model.get_connections() + self.name_to_esdl_id_map = dict() super().__init__(*args, **kwargs) + def get_connections(self): + return self.__connections + @property def esdl_bytes_string(self) -> bytes: """ diff --git a/src/mesido/esdl/esdl_model_base.py b/src/mesido/esdl/esdl_model_base.py index 3b6e655bc..b4be501c1 100644 --- a/src/mesido/esdl/esdl_model_base.py +++ b/src/mesido/esdl/esdl_model_base.py @@ -24,6 +24,7 @@ class _ESDLModelBase(_Model): primary_port_name_convention = "prim" secondary_port_name_convention = "sec" + __connections = None def _esdl_convert( self, converter: _AssetToComponentBase, assets: Dict, name_to_id_map: Dict, prefix: str @@ -320,6 +321,7 @@ def _esdl_convert( # therefore do the nodes first, and do all remaining connections # after. connections = set() + model_connections = [] for asset in [*node_assets, *bus_assets, *gas_node_assets]: component = getattr(self, asset.name) @@ -353,6 +355,12 @@ def _esdl_convert( == "Pipe" ): self.connect(getattr(component, node_suf)[i], port_map[connected_to.id]) + model_connections.append( + [ + getattr(component, node_suf)[i].name, + port_map[connected_to.id].name, + ] + ) elif connected_to.id not in list(port_map.keys()): # If The asset is not in the # port map means that there is a direct node to node connection with a @@ -374,12 +382,23 @@ def _esdl_convert( getattr(component, node_suf)[i], getattr(getattr(self, connected_node_asset.name), node_suf)[idx], ) + model_connections.append( + [ + getattr(component, node_suf)[i], + getattr(getattr(self, connected_node_asset.name), node_suf)[ + idx + ], + ] + ) else: # If the Connected asset is not of type pipe, there might be # logical link like source to node. self.connect_logical_links( getattr(component, node_suf)[i], port_map[connected_to.id] ) + model_connections.append( + [getattr(component, node_suf)[i], port_map[connected_to.id]] + ) connections.add(conn) i += 1 elif isinstance(port.carrier, esdl.ElectricityCommodity): @@ -394,6 +413,9 @@ def _esdl_convert( self.connect( getattr(component, elec_node_suf)[i], port_map[connected_to.id] ) + model_connections.append( + [getattr(component, elec_node_suf)[i], port_map[connected_to.id]] + ) elif connected_to.id not in list(port_map.keys()): for node in bus_assets: if connected_to.id in [node.in_ports[0].id, node.out_ports[0].id]: @@ -413,11 +435,22 @@ def _esdl_convert( idx ], ) + model_connections.append( + [ + getattr(component, elec_node_suf)[i], + getattr( + getattr(self, connected_node_asset.name), elec_node_suf + )[idx], + ] + ) else: self.connect_logical_links( getattr(component, elec_node_suf)[i], port_map[connected_to.id], ) + model_connections.append( + [getattr(component, elec_node_suf)[i], port_map[connected_to.id]] + ) connections.add(conn) i += 1 elif isinstance(port.carrier, esdl.GasCommodity): @@ -432,6 +465,9 @@ def _esdl_convert( self.connect( getattr(component, gas_node_suf)[i], port_map[connected_to.id] ) + model_connections.append( + [getattr(component, gas_node_suf)[i], port_map[connected_to.id]] + ) elif connected_to.id not in list(port_map.keys()): for node in gas_node_assets: if connected_to.id in [node.in_ports[0].id, node.out_ports[0].id]: @@ -451,10 +487,21 @@ def _esdl_convert( idx ], ) + model_connections.append( + [ + getattr(component, gas_node_suf)[i], + getattr(getattr(self, connected_node_asset.name), gas_node_suf)[ + idx + ], + ] + ) else: self.connect_logical_links( getattr(component, gas_node_suf)[i], port_map[connected_to.id] ) + model_connections.append( + [getattr(component, gas_node_suf)[i], port_map[connected_to.id]] + ) connections.add(conn) i += 1 else: @@ -508,4 +555,12 @@ def _esdl_convert( self.connect(port_map[port.id], port_map[connected_to.id]) else: self.connect_logical_links(port_map[port.id], port_map[connected_to.id]) + model_connections.append( + [port_map[port.id].name, port_map[connected_to.id].name] + ) connections.add(conn) + + self.__connections = model_connections + + def get_connections(self): + return self.__connections diff --git a/src/mesido/esdl/esdl_parser.py b/src/mesido/esdl/esdl_parser.py index 1a1707ee9..d19234af4 100644 --- a/src/mesido/esdl/esdl_parser.py +++ b/src/mesido/esdl/esdl_parser.py @@ -56,7 +56,7 @@ def read_esdl(self) -> None: id=x.id, id_number_mapping=id_to_idnumber_map[x.id], temperature=temperature, - type="milp", + type="heat", ) elif isinstance(x, esdl.esdl.ElectricityCommodity): if x.id not in id_to_idnumber_map: diff --git a/src/mesido/financial_mixin.py b/src/mesido/financial_mixin.py index 81a4cc155..ee64074b6 100644 --- a/src/mesido/financial_mixin.py +++ b/src/mesido/financial_mixin.py @@ -4,6 +4,7 @@ import casadi as ca from mesido.base_component_type_mixin import BaseComponentTypeMixin +from mesido.heat_network_common import NodeConnectionDirection import numpy as np @@ -77,6 +78,7 @@ def __init__(self, *args, **kwargs): # Variable for realized revenue self._asset_revenue_map = {} + self._asset_revenue_variable_port_map = {} self.__asset_revenue_var = {} self.__asset_revenue_nominals = {} self.__asset_revenue_bounds = {} @@ -95,6 +97,163 @@ def pre(self): parameters = self.parameters(0) bounds = self.bounds() + for _owner, assets in self.energy_system_owners.items(): + # All assets already have individual costs that can be used for the business cases + # Here we create revenue variables for the locations where energy is consumed or + # exhanged between owners + + # Create revenue variable for every demand type asset + # The assumption is that if an organizational asset own the demand asset it has the + # single right to sell energy to that demand. + for asset in assets: + if asset in [ + *self.energy_system_components.get("heat_demand", []), + *self.energy_system_components.get("gas_demand", []), + *self.energy_system_components.get("electricity_demand", []), + ]: + # We check if we can find the carrier information used for that demand + # Note that we thus assume that the get_{commodity}_carriers() functions are + # set. In most practical cases this means using the ESDLExtraVarsMixin + var_name = None + nominal_power = None + carrier_name = None + for _id, attr in self.get_electricity_carriers().items(): + # check if we can find the carried id + if ( + attr["id_number_mapping"] + == parameters[f"{asset}.ElectricityIn.carrier_id"] + ): + # Note that we are including the port name at the end as we can have + # multiple revenue variables for a single asset as there can be + # multiple connections with different organizations over the different + # ports. + var_name = f"{asset}__revenue_ElectricityIn" + carrier_name = attr["name"] + nominal_power = self.variable_nominal(f"{asset}.Electricity_demand") + # Here we same the variable which should be used for the energy flow + # between the organizations and thus for computing the revenue. + self._asset_revenue_variable_port_map[var_name] = ( + f"{asset}.ElectricityIn.Power" + ) + for _id, attr in self.get_gas_carriers().items(): + if attr["id_number_mapping"] == parameters[f"{asset}.GasIn.carrier_id"]: + var_name = f"{asset}__revenue_GasIn" + carrier_name = attr["name"] + nominal_power = self.variable_nominal(f"{asset}.Gas_demand") + self._asset_revenue_variable_port_map[var_name] = f"{asset}.GasIn.Q" + for _id, attr in self.get_heat_carriers().items(): + if attr["id_number_mapping"] == parameters[f"{asset}.HeatIn.carrier_id"]: + var_name = f"{asset}__revenue_HeatIn" + carrier_name = attr["name"] + nominal_power = self.variable_nominal(f"{asset}.Heat_demand") + self._asset_revenue_variable_port_map[var_name] = f"{asset}.HeatIn.Heat" + # Only if we could fing the carrier information we create the revenue variable + # and thus the constraints. + if carrier_name is not None: + self._asset_revenue_map[asset] = [var_name] + self.__asset_revenue_var[var_name] = ca.MX.sym(var_name) + try: + price_profile = self.get_timeseries( + f"{carrier_name}.price_profile" + ).values + except KeyError: + price_profile = np.ones(len(self.times())) + avg_price = np.mean(price_profile) if np.mean(price_profile) else 1.0 + self.__asset_revenue_nominals[var_name] = avg_price * nominal_power + self.__asset_revenue_bounds[var_name] = (0.0, np.inf) + + # Find the assets that exchange energy with other owners and create a revenue variable + # for them as well We loop over the assets + for asset in assets: + # This is something extra I created. It is a nested list in the form of: + # [[asset1.HeatIn, asset2.HeatOut], [asset3.HeatIn, asset4.HeatOut], ...] + conn = self.get_connections() + # We loop over the connections (yes this is inefficient) + for connection in conn: + # Here we check if our asset is in the connection hence we find a connected + # asset + if asset == connection[0].split(".")[0] or asset == connection[1].split(".")[0]: + other_asset = ( + connection[0].split(".")[0] + if asset == connection[1].split(".")[0] + else connection[1].split(".")[0] + ) + # This bit of code is used to go from for example: + # HeatExchanger.Secondary.HeatIn -> Secondary.HeatIn + temp = ( + connection[1].split(".")[1:] + if asset == connection[1].split(".")[0] + else connection[0].split(".")[1:] + ) + port = temp[0] + for x in temp[1:]: + port = port + "." + x + + # Here we check that the connected asset is not also owned by the same + # entity. If the other asset belongs to an organization we create a + # revenue variable to represent the cashflow associated with the energy + # exchanged between the two organizations. + if other_asset not in assets: + var_name = f"{asset}__revenue_{port}" + nominal_power = None + carrier_name = None + for _id, attr in self.get_electricity_carriers().items(): + self._asset_revenue_variable_port_map[var_name] = ( + f"{asset}.{port}.Power" + ) + if ( + attr["id_number_mapping"] + == parameters[f"{asset}.{port}.carrier_id"] + ): + nominal_power = self.variable_nominal(f"{asset}.{port}.Power") + carrier_name = attr["name"] + for _id, attr in self.get_gas_carriers().items(): + self._asset_revenue_variable_port_map[var_name] = ( + f"{asset}.{port}.Q" + ) + if ( + attr["id_number_mapping"] + == parameters[f"{asset}.{port}.carrier_id"] + ): + nominal_power = self.variable_nominal(f"{asset}.{port}.Q") + carrier_name = attr["name"] + for _id, attr in self.get_heat_carriers().items(): + self._asset_revenue_variable_port_map[var_name] = ( + f"{asset}.{port}.Heat" + ) + if ( + attr["id_number_mapping"] + == parameters[f"{asset}.{port}.carrier_id"] + ): + nominal_power = self.variable_nominal(f"{asset}.{port}.Heat") + carrier_name = attr["name"] + if carrier_name is not None: + # Please note that we create a dict mapping with a list. As a single + # asset can have multiple revenue variables as it can have + # connections with other organizations on different ports. + try: + self._asset_revenue_map[asset].append(var_name) + except KeyError: + self._asset_revenue_map[asset] = [var_name] + self.__asset_revenue_var[var_name] = ca.MX.sym(var_name) + try: + price_profile = self.get_timeseries( + f"{carrier_name}.price_profile" + ).values + except NameError: + price_profile = np.ones(len(self.times())) + except KeyError: + price_profile = np.ones(len(self.times())) + avg_price = ( + np.mean(price_profile) if np.mean(price_profile) else 1.0 + ) + self.__asset_revenue_nominals[var_name] = ( + avg_price * nominal_power + if nominal_power is not None + else 1.0e2 + ) + self.__asset_revenue_bounds[var_name] = (-np.inf, np.inf) + # Making the cost variables; fixed_operational_cost, variable_operational_cost, # installation_cost and investment_cost for asset_name in [ @@ -391,37 +550,6 @@ def pre(self): else 1.0e2 ) - # Realized revenue - if (asset_name) in [ - *self.energy_system_components.get("electricity_demand", []), - *self.energy_system_components.get("gas_demand", []), - ]: - - carrier_name = None - for _id, attr in self.get_electricity_carriers().items(): - if attr["id_number_mapping"] == parameters[f"{asset_name}.id_mapping_carrier"]: - carrier_name = attr["name"] - for _id, attr in self.get_gas_carriers().items(): - if attr["id_number_mapping"] == parameters[f"{asset_name}.id_mapping_carrier"]: - carrier_name = attr["name"] - if carrier_name is not None: - asset_revenue_var = f"{asset_name}__revenue" - self._asset_revenue_map[asset_name] = asset_revenue_var - self.__asset_revenue_var[asset_revenue_var] = ca.MX.sym(asset_revenue_var) - self.__asset_revenue_bounds[asset_revenue_var] = ( - 0.0, - np.inf, - ) - self.__asset_revenue_nominals[asset_revenue_var] = ( - max( - np.mean(self.get_timeseries(f"{carrier_name}.price_profile").values) - * nominal_fixed_operational, - 1.0e2, - ) - if nominal_fixed_operational is not None - else 1.0e2 - ) - for asset in [ *self.energy_system_components.get("heat_source", []), *self.energy_system_components.get("heat_demand", []), @@ -1305,44 +1433,91 @@ def __revenue_constraints(self, ensemble_member): # TODO: add fixed price default from ESDL in case no price profile is defined. parameters = self.parameters(ensemble_member) - for demand in [ - *self.energy_system_components.get("gas_demand", []), - *self.energy_system_components.get("electricity_demand", []), - ]: - - carrier_name = None - for _id, attr in self.get_electricity_carriers().items(): - if attr["id_number_mapping"] == parameters[f"{demand}.id_mapping_carrier"]: - carrier_name = attr["name"] - cost_multiplier = 1 / 3600.0 # priceprofile electricity is EUR/Wh - for _id, attr in self.get_gas_carriers().items(): - if attr["id_number_mapping"] == parameters[f"{demand}.id_mapping_carrier"]: - carrier_name = attr["name"] - cost_multiplier = 1.0 # priceprofile gas is in EUR/g - if carrier_name is not None: - price_profile = self.get_timeseries(f"{carrier_name}.price_profile").values - - if demand in self.energy_system_components.get("gas_demand", []): - energy_flow = self.__state_vector_scaled( - f"{demand}.Gas_demand_mass_flow", ensemble_member # g/s - ) + # Here we compute the revenue by looping over all the revenue variables. + for asset, variable_revenue_vars in self._asset_revenue_map.items(): + # note we need a nested loop as a asset can have multiple revenue variables + for variable_revenue_var in variable_revenue_vars: + carrier_name = None + port = None + # We find the associated carrier and port (is the full string with asset name + # included) + port_name = self._asset_revenue_variable_port_map[variable_revenue_var] + for _id, attr in self.get_electricity_carriers().items(): + if ( + attr["id_number_mapping"] + == parameters[f"{asset}.{port_name.split('.')[-2]}.carrier_id"] + ): + carrier_name = attr["name"] + port = f"{self._asset_revenue_variable_port_map[variable_revenue_var]}" + for _id, attr in self.get_gas_carriers().items(): + if ( + attr["id_number_mapping"] + == parameters[f"{asset}.{port_name.split('.')[-2]}.carrier_id"] + ): + carrier_name = attr["name"] + port = f"{self._asset_revenue_variable_port_map[variable_revenue_var]}" + for _id, attr in self.get_heat_carriers().items(): + temp = self._asset_revenue_variable_port_map[variable_revenue_var].split(".") + full_port = temp[1] + for x in temp[2:-1]: + full_port = full_port + "." + x + if attr["id_number_mapping"] == parameters[f"{asset}.{full_port}.carrier_id"]: + carrier_name = attr["name"] + port = f"{self._asset_revenue_variable_port_map[variable_revenue_var]}" - elif demand in self.energy_system_components.get("electricity_demand", []): - energy_flow = self.__state_vector_scaled( - f"{demand}.Electricity_demand", ensemble_member + # This if statement should not be needed... + if carrier_name is not None: + try: + price_profile = self.get_timeseries(f"{carrier_name}.price_profile").values + except KeyError: + price_profile = np.ones(len(self.times())) + energy_flow = self.__state_vector_scaled(f"{port}", ensemble_member) + # Here we find the sign of the flow, meaning that energy going into the asset is + # negative revenue. + energy_flow_sign = 1.0 + if "In" in port and asset not in [ + *self.energy_system_components.get("heat_demand", []), + *self.energy_system_components.get("gas_demand", []), + *self.energy_system_components.get("electricity_demand", []), + ]: + energy_flow_sign = -1.0 + # For nodes/busses we do this by checking the type of port in the topology + # object. It is a bit of a hassle as we need to find the index of the port of + # the node for that we find the last digit (which by definition is the port + # index) in the port name string. + if asset in self.energy_system_topology.nodes.keys(): + port_number = int([x for x in port if x.isdigit()][-1]) + if ( + self.energy_system_topology.nodes[f"{asset}"][port_number][1] + == NodeConnectionDirection.IN + ): + energy_flow_sign = -1.0 + if asset in self.energy_system_topology.gas_nodes.keys(): + port_number = int([x for x in port if x.isdigit()][-1]) + if ( + self.energy_system_topology.gas_nodes[f"{asset}"][port_number][1] + == NodeConnectionDirection.IN + ): + energy_flow_sign = -1.0 + if asset in self.energy_system_topology.busses.keys(): + port_number = int([x for x in port if x.isdigit()][-1]) + if ( + self.energy_system_topology.busses[f"{asset}"][port_number][1] + == NodeConnectionDirection.IN + ): + energy_flow_sign = -1.0 + variable_revenue = self.extra_variable(variable_revenue_var, ensemble_member) + nominal = self.variable_nominal(variable_revenue_var) + + sum = 0.0 + timesteps = np.diff(self.times()) / 3600.0 + for i in range(1, len(self.times())): + sum += price_profile[i] * energy_flow[i] * timesteps[i - 1] + + constraints.append( + ((variable_revenue - sum * energy_flow_sign) / (nominal), 0.0, 0.0) ) - variable_revenue_var = self._asset_revenue_map[demand] - variable_revenue = self.extra_variable(variable_revenue_var, ensemble_member) - nominal = self.variable_nominal(variable_revenue_var) - - sum = 0.0 - timesteps = np.diff(self.times()) * cost_multiplier - for i in range(1, len(self.times())): - sum += price_profile[i] * energy_flow[i] * timesteps[i - 1] - - constraints.append(((variable_revenue - sum) / (nominal), 0.0, 0.0)) - return constraints def path_constraints(self, ensemble_member): diff --git a/src/mesido/pycml/component_library/milp/_internal/heat_component.py b/src/mesido/pycml/component_library/milp/_internal/heat_component.py index ed8ef61b3..3e558a673 100644 --- a/src/mesido/pycml/component_library/milp/_internal/heat_component.py +++ b/src/mesido/pycml/component_library/milp/_internal/heat_component.py @@ -22,6 +22,8 @@ def __init__(self, name, **modifiers): super().__init__(name, **modifiers) self.state = 1 + self.owner = "__" + self.variable_operational_cost_coefficient = 0.0 self.fixed_operational_cost_coefficient = 0.0 self.investment_cost_coefficient = 0.0 diff --git a/src/mesido/pycml/component_library/milp/electricity/electricity_base.py b/src/mesido/pycml/component_library/milp/electricity/electricity_base.py index 28c7a6a5f..3675a8b38 100644 --- a/src/mesido/pycml/component_library/milp/electricity/electricity_base.py +++ b/src/mesido/pycml/component_library/milp/electricity/electricity_base.py @@ -12,6 +12,8 @@ class ElectricityPort(ElectricityComponent, Connector): def __init__(self, name, **modifiers): super().__init__(name, **modifiers) + self.carrier_id = -1 + self.add_variable(Variable, "Power") self.add_variable(Variable, "I") self.add_variable(Variable, "V", min=0.0) diff --git a/src/mesido/pycml/component_library/milp/electricity/electricity_node.py b/src/mesido/pycml/component_library/milp/electricity/electricity_node.py index 2aeabf80c..17ec8bc25 100644 --- a/src/mesido/pycml/component_library/milp/electricity/electricity_node.py +++ b/src/mesido/pycml/component_library/milp/electricity/electricity_node.py @@ -24,6 +24,7 @@ def __init__(self, name, **modifiers): self.n = 2 assert self.n >= 2 + self.carrier_id = -1 self.add_variable(ElectricityPort, "ElectricityConn", self.n) self.add_variable(Variable, "V", min=0.0, nominal=self.voltage_nominal) @@ -33,3 +34,6 @@ def __init__(self, name, **modifiers): # Because the orientation of the connected cables are important to setup the energy # conservation, these constraints are added in the mixin. + + for i in range(1, self.n + 1): + self.ElectricityConn[i].carrier_id = self.carrier_id diff --git a/src/mesido/pycml/component_library/milp/gas/gas_base.py b/src/mesido/pycml/component_library/milp/gas/gas_base.py index 3583be82b..3e620ad73 100644 --- a/src/mesido/pycml/component_library/milp/gas/gas_base.py +++ b/src/mesido/pycml/component_library/milp/gas/gas_base.py @@ -14,6 +14,9 @@ def __init__(self, name, **modifiers): super().__init__(name, **modifiers) # TODO: think of more elegant approach for Q_shadow, currently required to ensure that # every port has a unique variable to make the correct port mapping + + self.carrier_id = -1 + self.add_variable(Variable, "Q") # [m3/s] self.add_variable(Variable, "mass_flow") # [g/s] self.add_variable(Variable, "H") # [m] diff --git a/src/mesido/pycml/component_library/milp/gas/gas_node.py b/src/mesido/pycml/component_library/milp/gas/gas_node.py index b68246301..9a8e31e1a 100644 --- a/src/mesido/pycml/component_library/milp/gas/gas_node.py +++ b/src/mesido/pycml/component_library/milp/gas/gas_node.py @@ -20,6 +20,7 @@ def __init__(self, name, **modifiers): self.n = 2 assert self.n >= 2 + self.carrier_id = -1 self.add_variable(GasPort, "GasConn", self.n) self.add_variable(Variable, "H", min=0.0) @@ -29,3 +30,6 @@ def __init__(self, name, **modifiers): # Because the orientation of the connected pipes are important to setup the mass # conservation, these constraints are added in the mixin. + + for i in range(1, self.n + 1): + self.GasConn[i].carrier_id = self.carrier_id diff --git a/src/mesido/pycml/component_library/milp/heat/heat_port.py b/src/mesido/pycml/component_library/milp/heat/heat_port.py index 9ba7ef285..f5a9c61f5 100644 --- a/src/mesido/pycml/component_library/milp/heat/heat_port.py +++ b/src/mesido/pycml/component_library/milp/heat/heat_port.py @@ -12,6 +12,8 @@ class HeatPort(HeatComponent, Connector): def __init__(self, name, **modifiers): super().__init__(name, **modifiers) + self.carrier_id = -1 + self.add_variable(Variable, "Heat") self.add_variable(Variable, "Q") self.add_variable(Variable, "H") diff --git a/src/mesido/pycml/component_library/milp/heat/node.py b/src/mesido/pycml/component_library/milp/heat/node.py index 1c767077f..1b098908d 100644 --- a/src/mesido/pycml/component_library/milp/heat/node.py +++ b/src/mesido/pycml/component_library/milp/heat/node.py @@ -20,6 +20,7 @@ def __init__(self, name, **modifiers): self.n = 2 assert self.n >= 2 + self.carrier_id = -1 self.add_variable(HeatPort, "HeatConn", self.n) self.add_variable(Variable, "H") @@ -31,3 +32,6 @@ def __init__(self, name, **modifiers): for i in range(1, self.n + 1): self.add_equation(self.HeatConn[i].H - self.H) # Q and Heat to be set in the mixin + + for i in range(1, self.n + 1): + self.HeatConn[i].carrier_id = self.carrier_id diff --git a/tests/models/heat_exchange/input/timeseries.csv b/tests/models/heat_exchange/input/timeseries.csv new file mode 100644 index 000000000..fd7ed1290 --- /dev/null +++ b/tests/models/heat_exchange/input/timeseries.csv @@ -0,0 +1,4 @@ +DateTime,HeatingDemand_3322,HeatingDemand_18aa,heat1,heat1_ret,heat2,heat2_ret +01-01-2019 00:00, 350000., 350000.,1.,1.,2.,2. +01-01-2019 01:00, 350000., 350000.,1.,1.,2.,2. +01-01-2019 02:00, 350000., 350000.,1.,1.,2.,2. diff --git a/tests/models/heat_exchange/model/network_with_parties.esdl b/tests/models/heat_exchange/model/network_with_parties.esdl new file mode 100644 index 000000000..7bf18466c --- /dev/null +++ b/tests/models/heat_exchange/model/network_with_parties.esdl @@ -0,0 +1,298 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/models/heat_exchange/src/run_heat_exchanger.py b/tests/models/heat_exchange/src/run_heat_exchanger.py index a1847bcac..d7f22afde 100644 --- a/tests/models/heat_exchange/src/run_heat_exchanger.py +++ b/tests/models/heat_exchange/src/run_heat_exchanger.py @@ -1,8 +1,10 @@ +from mesido.esdl.esdl_additional_vars_mixin import ESDLAdditionalVarsMixin from mesido.esdl.esdl_mixin import ESDLMixin from mesido.esdl.esdl_parser import ESDLFileParser from mesido.esdl.profile_parser import ProfileReaderFromFile from mesido.head_loss_class import HeadLossOption from mesido.physics_mixin import PhysicsMixin +from mesido.techno_economic_mixin import TechnoEconomicMixin from mesido.workflows.io.write_output import ScenarioOutput import numpy as np @@ -70,7 +72,8 @@ def function(self, optimization_problem, ensemble_member): class HeatProblem( ScenarioOutput, _GoalsAndOptions, - PhysicsMixin, + ESDLAdditionalVarsMixin, + TechnoEconomicMixin, LinearizedOrderGoalProgrammingMixin, GoalProgrammingMixin, ESDLMixin, @@ -208,7 +211,7 @@ def constraints(self, ensemble_member): class HeatProblemTvarDisableHEX( _GoalsAndOptions, - PhysicsMixin, + TechnoEconomicMixin, LinearizedOrderGoalProgrammingMixin, GoalProgrammingMixin, ESDLMixin, @@ -315,9 +318,9 @@ def constraints(self, ensemble_member): if __name__ == "__main__": solution = run_optimization_problem( HeatProblem, - esdl_file_name="heat_exchanger.esdl", + esdl_file_name="network_with_parties.esdl", esdl_parser=ESDLFileParser, profile_reader=ProfileReaderFromFile, - input_timeseries_file="timeseries_import.xml", + input_timeseries_file="timeseries.csv", ) results = solution.extract_results() diff --git a/tests/models/heatpump/src/run_heat_pump.py b/tests/models/heatpump/src/run_heat_pump.py index bb5f23e58..a5508b0ee 100644 --- a/tests/models/heatpump/src/run_heat_pump.py +++ b/tests/models/heatpump/src/run_heat_pump.py @@ -127,6 +127,9 @@ def solver_options(self): highs_options["mip_rel_gap"] = 0.001 return options + def times(self, variable=None) -> np.ndarray: + return super().times(variable=variable)[:5] + def temperature_carriers(self): return self.esdl_carriers # geeft terug de carriers met multiple temperature options