diff --git a/examples/00_Minmal/minimal_example_DSM.py b/examples/00_Minmal/minimal_example_DSM.py new file mode 100644 index 000000000..c43e0fda1 --- /dev/null +++ b/examples/00_Minmal/minimal_example_DSM.py @@ -0,0 +1,94 @@ +""" +This script shows how to use the flixopt framework to model a super minimalistic energy system. +""" + +import numpy as np +import pandas as pd +from rich.pretty import pprint + +import flixopt as fx +if __name__ == '__main__': + # --- Define the Flow System, that will hold all elements, and the time steps you want to model --- + timesteps = pd.date_range('2020-01-01', periods=24, freq='h') + flow_system = fx.FlowSystem(timesteps) + + # --- Define Thermal Load Profile --- + # Load profile (e.g., kW) for heating demand over time + thermal_load_profile = np.array([80, 80, 80, 80, 80, 80, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 80, 80, 80, 80, 80, 80]) + #thermal_load_profile = np.array([100, 100, 100, 100, 100, 100, 120, 120, 120, 100, 100, 100, 100, 100, 100, 80, 80, 80, 100, 100, 100, 100, 100, 100]) + #thermal_load_profile = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80]) + + # --- Define Energy Buses --- + # These are balancing nodes (inputs=outputs) and balance the different energy carriers your system + flow_system.add_elements(fx.Bus('District Heating'), fx.Bus('Natural Gas')) + + # --- Define Objective Effect (Cost) --- + # Cost effect representing the optimization objective (minimizing costs) + cost_effect = fx.Effect('costs', '€', 'Cost', is_standard=True, is_objective=True) + + # --- Define Flow System Components --- + # Boiler component with thermal output (heat) and fuel input (gas) + boiler1 = fx.linear_converters.Boiler( + 'Boiler1', + eta=0.5, + Q_th=fx.Flow(label='Thermal Output', bus='District Heating', size=100), + Q_fu=fx.Flow(label='Fuel Input', bus='Natural Gas'), + ) + boiler2 = fx.linear_converters.Boiler( + 'Boiler2', + eta=1/3, + Q_th=fx.Flow(label='Thermal Output', bus='District Heating', size=50), + Q_fu=fx.Flow(label='Fuel Input', bus='Natural Gas'), + ) + + # Heat load component with a fixed thermal demand profile + heat_load = fx.DSMSink( + 'DSM Sink Heat Demand', + sink=fx.Flow(label='Heat Load', bus='District Heating', size=150), + initial_demand=thermal_load_profile, + maximum_cumulated_deficit = -50, + maximum_cumulated_surplus = 50, + maximum_flow_deficit = -20, + maximum_flow_surplus = 20, + relative_loss_per_hour_positive_charge_state = 0.05, + relative_loss_per_hour_negative_charge_state = 0.05, + penalty_costs_positive_charge_states=0, + penalty_costs_negative_charge_states=0.01, + forward_timeshift = 3, + backward_timeshift = 3 + ) + + # Gas source component with cost-effect per flow hour + gas_source = fx.Source( + 'Natural Gas Tariff', + source=fx.Flow(label='Gas Flow', bus='Natural Gas', size=1000, effects_per_flow_hour=0.04), # 0.04 €/kWh + ) + + # --- Build the Flow System --- + # Add all components and effects to the system + flow_system.add_elements(cost_effect, boiler1, boiler2, heat_load, gas_source) + + # --- Define, model and solve a Calculation --- + calculation = fx.FullCalculation('Simulation1', flow_system) + calculation.do_modeling() + #calculation.solve(fx.solvers.HighsSolver(0.01, 60)) + calculation.solve(fx.solvers.GurobiSolver(0.001, 60)) + + # --- Analyze Results --- + # Access the results of an element + df1 = calculation.results['costs'].filter_solution('time').to_dataframe() + + # Original plots + #calculation.results['District Heating'].plot_node_balance_pie() + calculation.results['District Heating'].plot_node_balance() + calculation.results['DSM Sink Heat Demand'].plot_DSM_sink() + + # Save the DSM Sink Heat Demand solution dataset to a CSV file + calculation.results['DSM Sink Heat Demand'].solution.to_dataframe().to_csv('results/DSM_Sink_Heat_Demand_results.csv') + + # Save results to a file + df2 = calculation.results['District Heating'].node_balance().to_dataframe() + #df2.to_csv('results/District Heating.csv') # Save results to csv + + # Print infos to the console. + pprint(calculation.summary) diff --git a/flixopt/__init__.py b/flixopt/__init__.py index b92766449..a647a3982 100644 --- a/flixopt/__init__.py +++ b/flixopt/__init__.py @@ -20,6 +20,7 @@ PiecewiseEffects, SegmentedCalculation, Sink, + DSMSink, Source, SourceAndSink, Storage, diff --git a/flixopt/commons.py b/flixopt/commons.py index 68412d6fe..d1e0801a0 100644 --- a/flixopt/commons.py +++ b/flixopt/commons.py @@ -8,6 +8,7 @@ from .components import ( LinearConverter, Sink, + DSMSink, Source, SourceAndSink, Storage, @@ -29,6 +30,7 @@ 'Effect', 'Source', 'Sink', + 'DSMSink', 'SourceAndSink', 'Storage', 'LinearConverter', diff --git a/flixopt/components.py b/flixopt/components.py index bdac6d2fb..6bcf1ef98 100644 --- a/flixopt/components.py +++ b/flixopt/components.py @@ -9,9 +9,10 @@ import numpy as np from . import utils +from .config import CONFIG from .core import NumericData, NumericDataTS, PlausibilityError, Scalar, TimeSeries from .elements import Component, ComponentModel, Flow -from .features import InvestmentModel, OnOffModel, PiecewiseModel +from .features import InvestmentModel, OnOffModel, PiecewiseModel, StateModel, PreventSimultaneousUsageModel from .interface import InvestParameters, OnOffParameters, PiecewiseConversion from .structure import SystemModel, register_class_for_io @@ -521,10 +522,10 @@ def _initial_and_final_charge_state(self): if utils.is_number(self.element.initial_charge_state): self.add( self._model.add_constraints( - self.charge_state.isel(time=0) == self.element.initial_charge_state, name=name - ), - name_short, - ) + self.charge_state.isel(time=0) == self.element.initial_charge_state, name=name + ), + name_short, + ) elif self.element.initial_charge_state == 'lastValueOfSim': self.add( self._model.add_constraints( @@ -635,3 +636,659 @@ def __init__(self, label: str, sink: Flow, meta_data: Optional[Dict] = None): """ super().__init__(label, inputs=[sink], meta_data=meta_data) self.sink = sink + +@register_class_for_io +class DSMSink(Sink): + """ + Used to model sinks with the ability to perform demand side management. + """ + + def __init__( + self, + label: str, + sink: Flow, + initial_demand: NumericData, + maximum_flow_deficit: NumericData, + maximum_flow_surplus: NumericData, + maximum_cumulated_deficit: NumericData, + maximum_cumulated_surplus: NumericData, + forward_timeshift: Scalar = None, + backward_timeshift: Scalar = None, + relative_loss_per_hour_positive_charge_state: NumericData = 0, + relative_loss_per_hour_negative_charge_state: NumericData = 0, + allow_mixed_charge_states: bool = False, + allow_parallel_charge_and_discharge: bool = False, + penalty_costs_positive_charge_states: NumericData = None, + penalty_costs_negative_charge_states: NumericData = None, + penalty_costs_positive_charge_rates: NumericData = 0, + penalty_costs_negative_charge_rates: NumericData = 0, + meta_data: Optional[Dict] = None + ): + """ + Args: + label: The label of the Element. Used to identify it in the FlowSystem + sink: input-flow of DSM sink after DSM + initial_demand: initial demand of DSM sink before DSM + maximum_flow_deficit: maximum that the supply flow can fall short of the demand + maximum_flow_surplus: maximum that the supply flow can exceed the demand + maximum_cumulated_deficit: maximum cumulated supply deficit in flow hours + maximum_cumulated_surplus: maximum cumulated supply surplus in flow hours + forward_timeshift: Maximum number of hours by which the demand can be shifted forward in time. Default is infinite. + backward_timeshift: Maximum number of hours by which the demand can be shifted backward in time. Default is infinite. + relative_loss_per_hour_positive_charge_state: loss per chargeState-Unit per hour for positive charge states of the virtual storage. The default is 0. + relative_loss_per_hour_negative_charge_state: loss per chargeState-Unit per hour for negative charge states of the virtual storage. The default is 0. + allow_mixed_charge_states: If True, positive and negative charge states can occur simultaneously. + If False, only one type of charge state is allowed at a time. The default is False. + allow_parallel_charge_and_discharge: If True, allows simultaneous charging and discharging in one timestep. + If False, charging and discharging cannot occur simultaneously. The default is False. + penalty_costs_positive_charge_states: penalty costs per flow hour for loss of comfort due to positive charge states of the virtual storage (e.g. increased room temperature). The default is a small epsilon. + penalty_costs_negative_charge_states: penalty costs per flow hour for loss of comfort due to negative charge states of the virtual storage (e.g. decreased room temperature). The default is a small epsilon. + penalty_costs_positive_charge_rates: penalty costs per flow hour for positive charge rates (e.g. costs for increasing supply above demand). The default is 0. + penalty_costs_negative_charge_rates: penalty costs per flow hour for negative charge rates (e.g. costs for decreasing supply below demand). The default is 0. + meta_data: used to store more information about the Element. Is not used internally, but saved in the results. Only use python native types. + """ + super().__init__( + label, + sink, + meta_data + ) + + self.initial_demand: NumericDataTS = initial_demand + self.forward_timeshift = forward_timeshift + self.backward_timeshift = backward_timeshift + self.maximum_flow_deficit: NumericDataTS = maximum_flow_deficit + self.maximum_flow_surplus: NumericDataTS = maximum_flow_surplus + self.maximum_cumulated_deficit: NumericDataTS = maximum_cumulated_deficit + self.maximum_cumulated_surplus: NumericDataTS = maximum_cumulated_surplus + + self.relative_loss_per_hour_positive_charge_state: NumericDataTS = relative_loss_per_hour_positive_charge_state + self.relative_loss_per_hour_negative_charge_state: NumericDataTS = relative_loss_per_hour_negative_charge_state + self.allow_mixed_charge_states = allow_mixed_charge_states + self.allow_parallel_charge_and_discharge = allow_parallel_charge_and_discharge + + self.penalty_costs_positive_charge_states: NumericDataTS = penalty_costs_positive_charge_states if penalty_costs_positive_charge_states is not None else CONFIG.modeling.EPSILON + self.penalty_costs_negative_charge_states: NumericDataTS = penalty_costs_negative_charge_states if penalty_costs_negative_charge_states is not None else CONFIG.modeling.EPSILON + self.penalty_costs_positive_charge_rates: NumericDataTS = penalty_costs_positive_charge_rates + self.penalty_costs_negative_charge_rates: NumericDataTS = penalty_costs_negative_charge_rates + + def create_model(self, model: SystemModel) -> 'DSMSinkModel': + self._plausibility_checks(model) + + # calculate amount of timesteps, by which demand can be shifted forward/backward in time + hours_per_step = model.hours_per_step + self.timesteps_forward = int(self.forward_timeshift/hours_per_step.values[0]) if self.forward_timeshift is not None else None + self.timesteps_backward = int(self.backward_timeshift/hours_per_step.values[0]) if self.backward_timeshift is not None else None + + self.model = DSMSinkModel(model, self) + return self.model + + def transform_data(self, flow_system: 'FlowSystem') -> None: + super().transform_data(flow_system) + self.initial_demand = flow_system.create_time_series( + f'{self.label_full}|initial_demand', + self.initial_demand, + ) + self.maximum_flow_surplus = flow_system.create_time_series( + f'{self.label_full}|maximum_flow_surplus', + self.maximum_flow_surplus, + ) + self.maximum_flow_deficit = flow_system.create_time_series( + f'{self.label_full}|maximum_flow_deficit', + self.maximum_flow_deficit, + ) + self.maximum_cumulated_deficit = flow_system.create_time_series( + f'{self.label_full}|maximum_cumulated_deficit', + self.maximum_cumulated_deficit, + needs_extra_timestep=True, + ) + self.maximum_cumulated_surplus = flow_system.create_time_series( + f'{self.label_full}|maximum_cumulated_surplus', + self.maximum_cumulated_surplus, + needs_extra_timestep=True, + ) + self.relative_loss_per_hour_negative_charge_state = flow_system.create_time_series( + f'{self.label_full}|relative_loss_per_hour_negative_charge_state', + self.relative_loss_per_hour_negative_charge_state, + ) + self.relative_loss_per_hour_positive_charge_state = flow_system.create_time_series( + f'{self.label_full}|relative_loss_per_hour_positive_charge_state', + self.relative_loss_per_hour_positive_charge_state, + ) + self.penalty_costs_negative_charge_states = flow_system.create_time_series( + f'{self.label_full}|penalty_costs_negative_charge_states', + self.penalty_costs_negative_charge_states, + ) + self.penalty_costs_positive_charge_states = flow_system.create_time_series( + f'{self.label_full}|penalty_costs_positive_charge_states', + self.penalty_costs_positive_charge_states, + ) + self.penalty_costs_negative_charge_rates = flow_system.create_time_series( + f'{self.label_full}|penalty_costs_negative_charge_rates', + self.penalty_costs_negative_charge_rates, + ) + self.penalty_costs_positive_charge_rates = flow_system.create_time_series( + f'{self.label_full}|penalty_costs_positive_charge_rates', + self.penalty_costs_positive_charge_rates, + ) + + def _plausibility_checks(self, model: SystemModel): + """ + Check for infeasible or uncommon combinations of parameters + """ + super()._plausibility_checks() + + hours_per_step = model.hours_per_step + + # Check timeshift parameters + if self.forward_timeshift != None or self.backward_timeshift != None: + if any(hours_per_step.values[0]!=hours_per_step.values): + raise ValueError( + f'{self.label_full}:' + f'limits to forward and backward timeshifts can only be used for timesteps of equal length' + ) + + if self.forward_timeshift is not None: + if self.forward_timeshift <= 0: + raise ValueError( + f'{self.label_full}: {self.forward_timeshift=} ' + f'must be positive.' + ) + if self.forward_timeshift%hours_per_step.values[0]!=0: + raise ValueError( + f'{self.label_full}: {self.forward_timeshift=} ' + f'must be a multiple of the timestep length.' + ) + # Check if forward timeshift is specified but maximum flow surplus is zero + if np.all(self.maximum_flow_surplus.active_data == 0): + logger.warning( + f'{self.label_full}: forward_timeshift is specified but maximum_flow_surplus is zero. ' + f'This is contradictory as you cannot shift demand forward if supply cannot exceed initial demand.' + ) + + if self.backward_timeshift is not None: + if self.backward_timeshift <= 0: + raise ValueError( + f'{self.label_full}: {self.backward_timeshift=} ' + f'must be positive.' + ) + if self.backward_timeshift%hours_per_step.values[0]!=0: + raise ValueError( + f'{self.label_full}: {self.backward_timeshift=} ' + f'must be a multiple of the timestep length.' + ) + # Check if backward timeshift is specified but maximum flow deficit is zero + if np.all(self.maximum_flow_deficit.active_data == 0): + logger.warning( + f'{self.label_full}: backward_timeshift is specified but maximum_flow_deficit is zero. ' + f'This is contradictory as you cannot shift demand backward if supply cannot fall short of initial demand.' + ) + + # Check maximum flow bounds + if np.any(self.maximum_flow_deficit.active_data > 0): + raise ValueError( + f'{self.label_full}: maximum_flow_deficit must be non-positive ' + f'(as it represents a deficit).' + ) + if np.any(self.maximum_flow_surplus.active_data < 0): + raise ValueError( + f'{self.label_full}: maximum_flow_surplus must be non-negative ' + f'(as it represents a surplus).' + ) + + # Check maximum cumulated bounds + if np.any(self.maximum_cumulated_deficit.active_data > 0): + raise ValueError( + f'{self.label_full}: maximum_cumulated_deficit must be non-positive ' + f'(as it represents a deficit).' + ) + if np.any(self.maximum_cumulated_surplus.active_data < 0): + raise ValueError( + f'{self.label_full}: maximum_cumulated_surplus must be non-negative ' + f'(as it represents a surplus).' + ) + + # Check loss rates + if np.any(self.relative_loss_per_hour_positive_charge_state.active_data < 0) or \ + np.any(self.relative_loss_per_hour_positive_charge_state.active_data > 1): + raise ValueError( + f'{self.label_full}: relative_loss_per_hour_positive_charge_state must be between 0 and 1 ' + f'(representing 0% to 100% loss).' + ) + if np.any(self.relative_loss_per_hour_negative_charge_state.active_data < 0) or \ + np.any(self.relative_loss_per_hour_negative_charge_state.active_data > 1): + raise ValueError( + f'{self.label_full}: relative_loss_per_hour_negative_charge_state must be between 0 and 1 ' + f'(representing 0% to 100% loss).' + ) + + # Check penalty costs + if np.any(self.penalty_costs_positive_charge_states.active_data < 0): + logger.warning( + f'{self.label_full}: penalty_costs_positive_charge_states should be non-negative.' + ) + if np.any(self.penalty_costs_negative_charge_states.active_data < 0): + logger.warning( + f'{self.label_full}: penalty_costs_negative_charge_states should be non-negative.' + ) + if np.any(self.penalty_costs_positive_charge_rates.active_data < 0): + logger.warning( + f'{self.label_full}: penalty_costs_positive_charge_rates should be non-negative.' + ) + if np.any(self.penalty_costs_negative_charge_rates.active_data < 0): + logger.warning( + f'{self.label_full}: penalty_costs_negative_charge_rates should be non-negative.' + ) + + # Check initial demand + if np.any(self.initial_demand.active_data < 0): + logger.warning( + f'{self.label_full}: initial_demand should be non-negative.' + ) + + # Check for zero bounds that would disable DSM + if np.all(self.maximum_flow_surplus.active_data == 0): + logger.warning( + f'{self.label_full}: maximum_flow_surplus is zero. ' + f'This could effectively disable DSM functionality as supply can never exceed the initial demand.' + ) + if np.all(self.maximum_flow_deficit.active_data == 0): + logger.warning( + f'{self.label_full}: maximum_flow_deficit is zero. ' + f'This could effectively disable DSM functionality as supply can never fall short of the initial demand.' + ) + if np.all(self.maximum_cumulated_surplus.active_data == 0) and np.all(self.maximum_cumulated_deficit.active_data == 0): + logger.warning( + f'{self.label_full}: Both maximum_cumulated_surplus and maximum_cumulated_deficit are zero. ' + f'This effectively disables DSM functionality as no energy can be stored.' + ) + + # Check for parallel charge and discharge + if self.allow_parallel_charge_and_discharge: + logger.warning( + f'{self.label_full}: allow_parallel_charge_and_discharge is True. ' + f'This might lead to unexpected behaviour.' + ) + + # Check for mixed charge states + if self.allow_mixed_charge_states: + logger.warning( + f'{self.label_full}: allow_mixed_charge_states is True. ' + f'This might lead to unexpected behaviour.' + ) + + # Check for zero penalty costs with non-zero bounds + if (np.all(self.penalty_costs_positive_charge_states.active_data == 0) and + not np.all(self.maximum_flow_surplus.active_data == 0) and + not np.any(self.relative_loss_per_hour_positive_charge_state.active_data > 0)): + logger.warning( + f'{self.label_full}: penalty_costs_positive_charge_states is zero and relative_loss_per_hour_positive_charge_state is zero but maximum_flow_surplus is non-zero. ' + f'This might lead to unexpected behavior as there is no cost and no loss associated with exceeding the initial demand.' + ) + if (np.all(self.penalty_costs_negative_charge_states.active_data == 0) and + not np.all(self.maximum_flow_deficit.active_data == 0)): + logger.warning( + f'{self.label_full}: penalty_costs_negative_charge_states is zero but maximum_flow_deficit is non-zero. ' + f'This might lead to unexpected behavior as there is no cost associated with falling short of the initial demand.' + ) + +class DSMSinkModel(ComponentModel): + """Model of DSM Sink""" + + def __init__(self, model: SystemModel, element: DSMSink): + super().__init__(model, element) + self.element: DSMSink = element + self.positive_charge_state: Optional[linopy.Variable] = None + self.negative_charge_state: Optional[linopy.Variable] = None + self.positive_charge_rate: Optional[linopy.Variable] = None + self.negative_charge_rate: Optional[linopy.Variable] = None + + self.is_positive_charge_state: Optional[linopy.Variable] = None + self.is_negative_charge_state: Optional[linopy.Variable] = None + + # Add state models for charge rates + self.is_positive_charge_rate: Optional[StateModel] = None + self.is_negative_charge_rate: Optional[StateModel] = None + self.prevent_simultaneous_charge_rates: Optional[PreventSimultaneousUsageModel] = None + + def do_modeling(self): + super().do_modeling() + + # Add variables for positive and negative charge rates + lb, ub = self.charge_rate_bounds + self.positive_charge_rate = self.add( + self._model.add_variables( + lower=lb, upper=ub, coords=self._model.coords, name=f'{self.label_full}|positive_charge_rate' + ), + 'positive_charge_rate', + ) + + lb, ub = self.discharge_rate_bounds + self.negative_charge_rate = self.add( + self._model.add_variables( + lower=lb, upper=ub, coords=self._model.coords, name=f'{self.label_full}|negative_charge_rate' + ), + 'negative_charge_rate', + ) + + # Add variables for negative charge states + lb, ub = self.negative_charge_state_bounds + self.negative_charge_state = self.add( + self._model.add_variables( + lower=lb, upper=ub, coords=self._model.coords_extra, name=f'{self.label_full}|negative_charge_state' + ), + 'negative_charge_state', + ) + + # Add variables for positive charge states + lb, ub = self.positive_charge_state_bounds + self.positive_charge_state = self.add( + self._model.add_variables( + lower=lb, upper=ub, coords=self._model.coords_extra, name=f'{self.label_full}|positive_charge_state' + ), + 'positive_charge_state' + ) + + positive_charge_state = self.positive_charge_state + negative_charge_state = self.negative_charge_state + etapos = 1 - self.element.relative_loss_per_hour_positive_charge_state.active_data + etaneg = 1 - self.element.relative_loss_per_hour_negative_charge_state.active_data + hours_per_step = self._model.hours_per_step + positive_charge_rate = self.positive_charge_rate + negative_charge_rate = self.negative_charge_rate + initial_demand = self.element.initial_demand.active_data + sink = self.element.sink.model.flow_rate + + # eq: positive_charge_state(t) + negative_charge_state(t) = positive_charge_state(t-1) * etapos + negative_charge_state(t-1) * etaneg + positive_charge_rate(t) + negative_charge_rate(t) + self.add( + self._model.add_constraints( + positive_charge_state.isel(time=slice(1, None)) + + negative_charge_state.isel(time=slice(1,None)) + == positive_charge_state.isel(time=slice(None, -1)) * (etapos ** hours_per_step) + + negative_charge_state.isel(time=slice(None,-1)) * (etaneg ** hours_per_step) + + positive_charge_rate * hours_per_step + + negative_charge_rate * hours_per_step, + name=f'{self.label_full}|charge_state', + ), + 'charge_state', + ) + + # eq: sink(t) = initial_demand(t) + positive_charge_rate(t) + negative_charge_rate(t) + self.add( + self._model.add_constraints( + sink + == positive_charge_rate + + negative_charge_rate + + initial_demand, + name=f'{self.label_full}|resulting_load_profile', + ), + 'resulting_load_profile', + ) + + # Add constraints for preventing mixed charge states if not allowed + if not self.element.allow_mixed_charge_states: + self._add_charge_state_exclusivity_constraints() + + # Add charge rate exclusivity constraints + if not self.element.allow_parallel_charge_and_discharge: + self._add_charge_rate_exclusivity_constraints() + + # Forward and backward timeshift constraints + self._add_timeshift_limits() + + # Initial and final charge state constraints + self._initial_and_final_charge_state() + + # Add penalty costs as effects for positive and negative charge states + penalty_costs_pos = self.element.penalty_costs_positive_charge_states.active_data + penalty_costs_neg = self.element.penalty_costs_negative_charge_states.active_data + + # Add effects for positive charge states + if np.any(penalty_costs_pos != 0): + # Multiply penalty costs with hours_per_step first to get a single coefficient per timestep + penalty_coeff_pos = penalty_costs_pos * hours_per_step + self._model.effects.add_share_to_penalty( + name = self.label_full, + # charge state is shifted backwards in time to apply penalty costs to the charge state at the end of a timestep + expression = (positive_charge_state.shift(time=-1).isel(time=slice(None,-1)) * penalty_coeff_pos).sum() + ) + + # Add effects for negative charge states + if np.any(penalty_costs_neg != 0): + # Multiply penalty costs with hours_per_step first to get a single coefficient per timestep + penalty_coeff_neg = - penalty_costs_neg * hours_per_step + self._model.effects.add_share_to_penalty( + name = self.label_full, + # charge state is shifted backwards in time to apply penalty costs to the charge state at the end of a timestep + expression = (negative_charge_state.shift(time=-1).isel(time=slice(None,-1)) * penalty_coeff_neg).sum() + ) + + # Add penalty costs as effects for positive and negative charge rates + penalty_costs_pos_rate = self.element.penalty_costs_positive_charge_rates.active_data + penalty_costs_neg_rate = self.element.penalty_costs_negative_charge_rates.active_data + + # Add effects for positive charge rates + if np.any(penalty_costs_pos_rate != 0): + # Multiply penalty costs with hours_per_step first to get a single coefficient per timestep + penalty_coeff_pos_rate = penalty_costs_pos_rate * hours_per_step + self._model.effects.add_share_to_penalty( + name = self.label_full, + expression = (positive_charge_rate * penalty_coeff_pos_rate).sum() + ) + + # Add effects for negative charge rates + if np.any(penalty_costs_neg_rate != 0): + # Multiply penalty costs with hours_per_step first to get a single coefficient per timestep + penalty_coeff_neg_rate = - penalty_costs_neg_rate * hours_per_step + self._model.effects.add_share_to_penalty( + name = self.label_full, + expression = (negative_charge_rate * penalty_coeff_neg_rate).sum() + ) + + def _add_charge_state_exclusivity_constraints(self): + """Add constraints to prevent simultaneous positive and negative charge states""" + #TODO: this method currently does not use the implemented statemodel and preventsimultaneous model + #because they do not work with variables using coords_extra + + # Add binary variables to track charge state type + self.is_positive_charge_state = self.add( + self._model.add_variables( + binary=True, + coords=self._model.coords_extra, + name=f'{self.label_full}|is_positive_charge_state' + ), + 'is_positive_charge_state' + ) + + self.is_negative_charge_state = self.add( + self._model.add_variables( + binary=True, + coords=self._model.coords_extra, + name=f'{self.label_full}|is_negative_charge_state' + ), + 'is_negative_charge_state' + ) + + positive_charge_state = self.positive_charge_state + negative_charge_state = self.negative_charge_state + + # If positive_charge_state > 0, then is_positive_charge_state must be 1 + self.add( + self._model.add_constraints( + positive_charge_state <= self.positive_charge_state_bounds[1] * self.is_positive_charge_state, + name=f'{self.label_full}|positive_charge_state_binary_upper' + ), + 'positive_charge_state_binary_upper' + ) + + # If is_positive_charge_state is 1, then positive_charge_state must be > 0 + self.add( + self._model.add_constraints( + positive_charge_state >= CONFIG.modeling.EPSILON * self.positive_charge_state_bounds[1] * self.is_positive_charge_state, # Small epsilon to avoid numerical issues + name=f'{self.label_full}|positive_charge_state_binary_lower' + ), + 'positive_charge_state_binary_lower' + ) + + # If negative_charge_state < 0, then is_negative_charge_state must be 1 + self.add( + self._model.add_constraints( + negative_charge_state >= self.negative_charge_state_bounds[0] * self.is_negative_charge_state, + name=f'{self.label_full}|negative_charge_state_binary_upper' + ), + 'negative_charge_state_binary_upper' + ) + + # If is_negative_charge_state is 1, then negative_charge_state must be < 0 + self.add( + self._model.add_constraints( + negative_charge_state <= CONFIG.modeling.EPSILON * self.negative_charge_state_bounds[0] * self.is_negative_charge_state, # Small epsilon to avoid numerical issues + name=f'{self.label_full}|negative_charge_state_binary_lower' + ), + 'negative_charge_state_binary_lower' + ) + + # Ensure only one type of charge state can be active at a time + self.add( + self._model.add_constraints( + self.is_positive_charge_state + self.is_negative_charge_state <= 1, + name=f'{self.label_full}|mutually_exclusive_charge_states' + ), + 'mutually_exclusive_charge_states' + ) + + def _add_charge_rate_exclusivity_constraints(self): + """Add constraints to prevent simultaneous positive and negative charge rates using StateModel and PreventSimultaneousUsageModel""" + # Create a time series of zeros to be used for both bounds + timeseries_zeros = np.zeros_like(self._model.coords[0], dtype=float) + + # Create StateModel for positive charge rate + self.is_positive_charge_rate = self.add( + StateModel( + model=self._model, + label_of_element=f'{self.label_full}|positive_charge_rate', + defining_variables=[self.positive_charge_rate], + defining_bounds=[(timeseries_zeros, self.charge_rate_bounds[1])], + use_off=False + ) + ) + self.is_positive_charge_rate.do_modeling() + + # Create StateModel for negative charge rate + self.is_negative_charge_rate = self.add( + StateModel( + model=self._model, + label_of_element=f'{self.label_full}|negative_charge_rate', + defining_variables=[-self.negative_charge_rate], # StateModel can only handle positive variables + defining_bounds=[(timeseries_zeros, -self.discharge_rate_bounds[0])], + use_off=False + ) + ) + self.is_negative_charge_rate.do_modeling() + + # Create PreventSimultaneousUsageModel for charge rates + self.prevent_simultaneous_charge_rates = self.add( + PreventSimultaneousUsageModel( + model=self._model, + variables=[self.is_positive_charge_rate.on, self.is_negative_charge_rate.on], + label_of_element=self.label_full, + label='PreventSimultaneousChargeRates' + ) + ) + self.prevent_simultaneous_charge_rates.do_modeling() + + def _add_timeshift_limits(self): + hours_per_step = self._model.hours_per_step + timesteps_forward = self.element.timesteps_forward + timesteps_backward = self.element.timesteps_backward + + etapos = 1 - self.element.relative_loss_per_hour_positive_charge_state.active_data + etaneg = 1 - self.element.relative_loss_per_hour_negative_charge_state.active_data + + positive_charge_state = self.positive_charge_state + negative_charge_state = self.negative_charge_state + + # Add constraints limiting the forward timeshift + if timesteps_forward is not None: + surplus_sum = 0 + for i in range(0, timesteps_forward): + surplus_sum += self.positive_charge_rate.shift(time=i).isel(time=slice(i,None)) * hours_per_step * (etapos ** (hours_per_step * i)) + # eq: positive_charge_state(t) <= sum over n (positive_charge_rate(t-n) * hours_per_step * (etapos ^ (hours_per_step * i)) + # where n ranges from 0 to timesteps_forward + # The positive charge state can't be any higher than the charge that was added during the last x timesteps minus the losses + # x is defined by the timesteps_forward + # This forces the virtual storage to discharge after the maximum timesteps that the demand is allowed to be shifted forward + self.add( + self._model.add_constraints( + positive_charge_state.isel(time = slice(1,None)) + <= surplus_sum, + name=f'{self.label_full}|limit_forward_timeshift' + ), + f'limit_forward_timeshift' + ) + + # Add constraints limiting the backwards timeshift + if timesteps_backward is not None: + deficit_sum = 0 + for i in range(0, timesteps_backward): + deficit_sum += self.negative_charge_rate.shift(time=i).isel(time=slice(i,None)) * hours_per_step * (etaneg ** (hours_per_step * i)) + # eq: -negative_charge_state(t) <= -sum over n (negative_charge_rate(t-n) * hours_per_step * (etaneg ^ (hours_per_step * i)) + # where n ranges from 0 to timesteps_backward + # The negative charge state can't be any higher than the deficit that was accumulated during the last x timesteps minus the losses. + # x is defined by the timesteps_backward. + # This forces the virtual storage to "recharge" (i. e. "discharge" the negative charge state) after the maximum timesteps that the demand is allowed to be shifted backward. + self.add( + self._model.add_constraints( + - negative_charge_state.isel(time=slice(1,None)) + <= - deficit_sum, + name=f'{self.label_full}|limit_backward_timeshift' + ), + f'limit_backward_timeshift' + ) + + def _initial_and_final_charge_state(self): + """Add constraints for initial and final charge states to be zero""" + # Set initial charge state to zero + self.add( + self._model.add_constraints( + self.positive_charge_state.isel(time=0) + self.negative_charge_state.isel(time=0) == 0, + name=f'{self.label_full}|initial_charge_state' + ), + 'initial_charge_state' + ) + + # Set final charge state to zero + self.add( + self._model.add_constraints( + self.positive_charge_state.isel(time=-1) + self.negative_charge_state.isel(time=-1) == 0, + name=f'{self.label_full}|final_charge_state' + ), + 'final_charge_state' + ) + + @property + def positive_charge_state_bounds(self) -> Tuple[NumericData, NumericData]: + return ( + 0, + self.element.maximum_cumulated_surplus.active_data, + ) + + @property + def negative_charge_state_bounds(self) -> Tuple[NumericData, NumericData]: + return ( + self.element.maximum_cumulated_deficit.active_data, + 0, + ) + + @property + def charge_rate_bounds(self) -> Tuple[NumericData, NumericData]: + return( + 0, + self.element.maximum_flow_surplus.active_data, + ) + + @property + def discharge_rate_bounds(self) -> Tuple[NumericData, NumericData]: + return( + self.element.maximum_flow_deficit.active_data, + 0, + ) \ No newline at end of file diff --git a/flixopt/results.py b/flixopt/results.py index d9eb5a654..48e566c24 100644 --- a/flixopt/results.py +++ b/flixopt/results.py @@ -9,6 +9,7 @@ import numpy as np import pandas as pd import plotly +import plotly.express import xarray as xr import yaml @@ -531,6 +532,17 @@ class ComponentResults(_NodeResults): def is_storage(self) -> bool: return self._charge_state in self._variable_names + @property + def is_DSM_sink(self) -> bool: + """Check if this component is a DSM sink by looking for characteristic variables""" + return ( + f'{self.label}|positive_charge_state' in self._variable_names + and f'{self.label}|negative_charge_state' in self._variable_names + and f'{self.label}|positive_charge_rate' in self._variable_names + and f'{self.label}|negative_charge_rate' in self._variable_names + and f'{self.label}|resulting_load_profile' in self._constraint_names + ) + @property def _charge_state(self) -> str: return f'{self.label}|charge_state' @@ -542,6 +554,98 @@ def charge_state(self) -> xr.DataArray: raise ValueError(f'Cant get charge_state. "{self.label}" is not a storage') return self.solution[self._charge_state] + def plot_DSM_sink( + self, + save: Union[bool, pathlib.Path] = False, + show: bool = True, + colors: plotting.ColorType = 'viridis', + engine: plotting.PlottingEngine = 'plotly', + ) -> plotly.graph_objs.Figure: + """ + Plots the operation of a DSM sink, showing the resulting heat load, initial demand, surplus/deficit, and cumulated flow deviation. + + Args: + save: Whether to save the plot or not. If a path is provided, the plot will be saved at that location. + show: Whether to show the plot or not. + colors: The colors to use for the plot. + engine: Plotting engine to use. Only 'plotly' is implemented atm. + + Raises: + ValueError: If the Component is not a DSM sink. + """ + if engine != 'plotly': + raise NotImplementedError( + f'Plotting engine "{engine}" not implemented for ComponentResults.plot_DSM_sink.' + ) + + if not self.is_DSM_sink: + raise ValueError(f'Cannot plot DSM sink. "{self.label}" is not a DSM sink') + + # Get the node balance and the initial demand + node_balance = - self.node_balance(with_last_timestep=True).to_dataframe() + initial_demand = self._calculation_results.flow_system[f'{self.label}|initial_demand'].to_dataframe(name='initial_demand') + + # Get surplus and deficit from the solution + surplus = self.solution[f'{self.label}|positive_charge_rate'].to_dataframe() + deficit = self.solution[f'{self.label}|negative_charge_rate'].to_dataframe() + + # Substract surplus from node blance + node_balance[f'{self.inputs[0]}'] = node_balance[f'{self.inputs[0]}'].values.flatten() - surplus[f'{self.label}|positive_charge_rate'].values.flatten() + + # Merge dataframes into one + data = pd.concat([node_balance, surplus, deficit], axis='columns') + + # Create figure with area plot for node balance + fig = plotting.with_plotly( + data, + mode='area', + colors=colors, + title=f'DSM sink behaviour for {self.label}', + xlabel='Time in h' + ) + + # Add initial demand + fig.add_trace( + plotly.graph_objs.Scatter( + x=initial_demand.index, + y=initial_demand['initial_demand'], + name='Initial Demand', + line=dict(dash='dash', color='black', shape='hv'), # 'hv' for horizontal-vertical steps + mode='lines' + ) + ) + + # Get colors for the cumulated flow + color_samples = plotly.express.colors.sample_colorscale(colors, 8) + cumulated_color = color_samples[3] # Use another color from viridis for cumulated flow + + # Get cumulated flow deviation from the model + cumulated_flow = pd.DataFrame(0, index=surplus.index, columns=['cumulated_flow']) + cumulated_flow['cumulated_flow'] = ( + self.solution[f'{self.label}|positive_charge_state'].values.flatten() + + self.solution[f'{self.label}|negative_charge_state'].values.flatten() + ) + + # Add cumulated flow deviation as diamonds on secondary y-axis + fig.add_trace( + plotly.graph_objs.Scatter( + x=cumulated_flow.index, + y=cumulated_flow.values.flatten(), + name='Virtual Charge State', + line=dict(width=3, color='black'), + mode='lines' + ) + ) + + return plotting.export_figure( + fig, + default_path=self._calculation_results.folder / f'{self.label} (DSM sink)', + default_filetype='.html', + user_path=None if isinstance(save, bool) else pathlib.Path(save), + show=show, + save=True if save else False, + ) + def plot_charge_state( self, save: Union[bool, pathlib.Path] = False,