Add Scenarios to flixOpt

docs/user-guide/Mathematical Notation/Investment.md

@@ -0,0 +1,115 @@
+# Investments
+
+## Current state
+$$
+\beta_{\text{invest}} \cdot \text{max}(\epsilon, \text V^{\text L}) \leq V \leq  \beta_{\text{invest}} \cdot \text V^{\text U}
+$$
+With:
+- $V$ = size
+- $V^{\text L}$ = minimum size
+- $V^{\text U}$ = maximum size
+- $\epsilon$ = epsilon, a small number (such as $1e^{-5}$)
+- $\beta_{invest} \in {0,1}$ = wether the size is invested or not
+
+_Please edit the use cases as needed_
+## Quickfix 1: Optimize the single best size overall
+### Single variable
+This is already possible and should be, as this is a needed use case
+An additional factor to when the size is actually available might me practical (Which indicates the (fixed) time of investment)
+## Math
+$$
+V(p) = V * a(p)
+$$
+with:
+- $V$ = size
+- $a(p)$ = factor for availlability per period
+
+Factor $a(p)$ is simply multiplied with relative minimum or maximum(t). This is already possible by doing this yourself.
+Effectively, the relative minimum or maximum are altered before using the same constraiints as before.
+THis might lead to some issues regariding minimum_load factor, or others, as the size is not 0 in a scenario where the component cant produce.
+**Therefore this might not be the best choice. See (#Variable per Scenario)
+
+## Variable per Scenario
+- **size** and **invest** as a variable per period $V(s)$ and $\beta_{invest}(s)$
+- with scenario $s \in S$
+
+### Usecase 1: Optimize the size for each Scenario independently
+Restrictions are seperatly for each scenario
+No changes needed. This could be the default behaviour.
+
+### Usecase 2: Optimize ONE size for ALL scenarios
+The size is the same globally, but not a scalar, but a variable per scenario $V(s)$
+#### 2a: The same size in all scenarios
+$$
+V(s) = V(s') \quad \forall s,s' \in S
+$$
+
+With:
+- $V(s)$ and $V(s')$ = size
+- $S$ = set of scenarios
+
+#### 2b: The same size, but can be 0 prior to the first increment
+- Find the Optimal time of investment.
+- Force an investment in a certain scenario (parameter optional as a list/array ob booleans)
+- Combine optional and minimum/maximum size to force an investment inside a range if scenarios
+
+$$
+\beta_{\text{invest}}(s) \leq \beta_{\text{invest}}(s+1) \quad \forall s \in \{1,2,\ldots,S-1\}
+$$
+
+$$
+V(s') - V(s) \leq M \cdot (2 - \beta_{\text{invest}}(s) - \beta_{\text{invest}}(s')) \quad \forall s, s' \in S
+$$
+$$
+V(s') - V(s) \geq M \cdot (2 - \beta_{\text{invest}}(s) - \beta_{\text{invest}}(s')) \quad \forall s, s' \in S
+$$
+
+This could be the default behaviour. (which would be consistent with other variables)
+
+
+### Switch
+
+$$
+\begin{aligned}
+& \text{SWITCH}_s \in \{0,1\} \quad \forall s \in \{1,2,\ldots,S\} \\
+& \sum_{s=1}^{S} \text{SWITCH}_s = 1 \\
+& \beta_{\text{invest}}(s) = \sum_{s'=1}^{s} \text{SWITCH}_{s'} \quad \forall s \in \{1,2,\ldots,S\} \\
+\end{aligned}
+$$
+
+$$
+\begin{aligned}
+& V(s) \leq V_{\text{actual}} \quad \forall s \in \{1,2,\ldots,S\} \\
+& V(s) \geq V_{\text{actual}} - M \cdot (1 - \beta_{\text{invest}}(s)) \quad \forall s \in \{1,2,\ldots,S\}
+\end{aligned}
+$$
+
+
+
+
+### Usecase 3: Find the best scenario to increment the size (Timing of the investment)
+The size can only increment once (based on a starting point). This allows to optimize the timing of an investment.
+#### Math
+Treat $\beta_{invest}$ like an ON/OFF variable, and introduce a SwitchOn, that can only be active once.
+
+*Thoughts:*
+- Treating $\beta_{invest}$ like an ON/OFF variable suggest using the already presentconstraints linked to On/OffModel
+- The timing could be constraint to be first in scenario x, or last in scenario y
+- Restrict the number of consecutive scenarios
+THis might needs the OnOffModel to be more generic (HOURS). Further, the span between scenarios needs to be weighted (like dt_in_hours), or the scenarios need to be measureable (integers)
+
+
+### Others
+
+#### Usecase 4: Only increase/decrease the size
+Start from a certain size. For each scenario, the size can increase, but never decrease. (Or the other way around).
+This would mean that a size expansion is possible,
+
+#### Usecase 5: Restrict the increment in size per scenario
+Restrict how much the size can increase/decrease for in scenario, based on the prior scenario.
+
+
+
+
+
+Many more are possible

examples/04_Scenarios/scenario_example.py

@@ -0,0 +1,123 @@
+"""
+This script shows how to use the flixopt framework to model a simple energy system.
+"""
+
+import numpy as np
+import pandas as pd
+from rich.pretty import pprint  # Used for pretty printing
+
+import flixopt as fx
+
+if __name__ == '__main__':
+    # Create datetime array starting from '2020-01-01' for the given time period
+    timesteps = pd.date_range('2020-01-01', periods=9, freq='h')
+    scenarios = pd.Index(['Base Case', 'High Demand'])
+
+    # --- Create Time Series Data ---
+    # Heat demand profile (e.g., kW) over time and corresponding power prices
+    heat_demand_per_h = pd.DataFrame({'Base Case':[30, 0, 90, 110, 110, 20, 20, 20, 20],
+                                      'High Demand':[30, 0, 100, 118, 125, 20, 20, 20, 20]}, index=timesteps)
+    power_prices = np.array([0.08, 0.09])
+
+    flow_system = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios)
+
+    # --- Define Energy Buses ---
+    # These represent nodes, where the used medias are balanced (electricity, heat, and gas)
+    flow_system.add_elements(fx.Bus(label='Strom'), fx.Bus(label='Fernwärme'), fx.Bus(label='Gas'))
+
+    # --- Define Effects (Objective and CO2 Emissions) ---
+    # Cost effect: used as the optimization objective --> minimizing costs
+    costs = fx.Effect(
+        label='costs',
+        unit='€',
+        description='Kosten',
+        is_standard=True,  # standard effect: no explicit value needed for costs
+        is_objective=True,  # Minimizing costs as the optimization objective
+    )
+
+    # CO2 emissions effect with an associated cost impact
+    CO2 = fx.Effect(
+        label='CO2',
+        unit='kg',
+        description='CO2_e-Emissionen',
+        specific_share_to_other_effects_operation={costs.label: 0.2},
+        maximum_operation_per_hour=1000,  # Max CO2 emissions per hour
+    )
+
+    # --- Define Flow System Components ---
+    # Boiler: Converts fuel (gas) into thermal energy (heat)
+    boiler = fx.linear_converters.Boiler(
+        label='Boiler',
+        eta=0.5,
+        Q_th=fx.Flow(label='Q_th', bus='Fernwärme', size=50, relative_minimum=0.1, relative_maximum=1),
+        Q_fu=fx.Flow(label='Q_fu', bus='Gas'),
+    )
+
+    # Combined Heat and Power (CHP): Generates both electricity and heat from fuel
+    chp = fx.linear_converters.CHP(
+        label='CHP',
+        eta_th=0.5,
+        eta_el=0.4,
+        P_el=fx.Flow('P_el', bus='Strom', size=60, relative_minimum=5 / 60),
+        Q_th=fx.Flow('Q_th', bus='Fernwärme'),
+        Q_fu=fx.Flow('Q_fu', bus='Gas'),
+    )
+
+    # Storage: Energy storage system with charging and discharging capabilities
+    storage = fx.Storage(
+        label='Storage',
+        charging=fx.Flow('Q_th_load', bus='Fernwärme', size=1000),
+        discharging=fx.Flow('Q_th_unload', bus='Fernwärme', size=1000),
+        capacity_in_flow_hours=fx.InvestParameters(fix_effects=20, fixed_size=30, optional=False),
+        initial_charge_state=0,  # Initial storage state: empty
+        relative_maximum_charge_state=np.array([80, 70, 80, 80, 80, 80, 80, 80, 80, 80]) * 0.01,
+        eta_charge=0.9,
+        eta_discharge=1,  # Efficiency factors for charging/discharging
+        relative_loss_per_hour=0.08,  # 8% loss per hour. Absolute loss depends on current charge state
+        prevent_simultaneous_charge_and_discharge=True,  # Prevent charging and discharging at the same time
+    )
+
+    # Heat Demand Sink: Represents a fixed heat demand profile
+    heat_sink = fx.Sink(
+        label='Heat Demand',
+        sink=fx.Flow(label='Q_th_Last', bus='Fernwärme', size=1, fixed_relative_profile=heat_demand_per_h),
+    )
+
+    # Gas Source: Gas tariff source with associated costs and CO2 emissions
+    gas_source = fx.Source(
+        label='Gastarif',
+        source=fx.Flow(label='Q_Gas', bus='Gas', size=1000, effects_per_flow_hour={costs.label: 0.04, CO2.label: 0.3}),
+    )
+
+    # Power Sink: Represents the export of electricity to the grid
+    power_sink = fx.Sink(
+        label='Einspeisung', sink=fx.Flow(label='P_el', bus='Strom', effects_per_flow_hour=-1 * power_prices)
+    )
+
+    # --- Build the Flow System ---
+    # Add all defined components and effects to the flow system
+    flow_system.add_elements(costs, CO2, boiler, storage, chp, heat_sink, gas_source, power_sink)
+
+    # Visualize the flow system for validation purposes
+    flow_system.plot_network(show=True)
+
+    # --- Define and Run Calculation ---
+    # Create a calculation object to model the Flow System
+    calculation = fx.FullCalculation(name='Sim1', flow_system=flow_system)
+    calculation.do_modeling()  # Translate the model to a solvable form, creating equations and Variables
+
+    # --- Solve the Calculation and Save Results ---
+    calculation.solve(fx.solvers.HighsSolver(mip_gap=0, time_limit_seconds=30))
+
+    # --- Analyze Results ---
+    calculation.results['Fernwärme'].plot_node_balance_pie()
+    calculation.results['Fernwärme'].plot_node_balance()
+    calculation.results['Storage'].plot_node_balance()
+    calculation.results.plot_heatmap('CHP(Q_th)|flow_rate')
+
+    # Convert the results for the storage component to a dataframe and display
+    df = calculation.results['Storage'].node_balance_with_charge_state()
+    print(df)
+
+    # Save results to file for later usage
+    calculation.results.to_file()

flixopt/calculation.py

@@ -44,19 +44,22 @@ def __init__(
         name: str,
         flow_system: FlowSystem,
         active_timesteps: Optional[pd.DatetimeIndex] = None,
+        selected_scenarios: Optional[pd.Index] = None,
         folder: Optional[pathlib.Path] = None,
     ):
         """
         Args:
             name: name of calculation
             flow_system: flow_system which should be calculated
-            active_timesteps: list with indices, which should be used for calculation. If None, then all timesteps are used.
+            active_timesteps: timesteps which should be used for calculation. If None, then all timesteps are used.
+            selected_scenarios: scenarios which should be used for calculation. If None, then all scenarios are used.
             folder: folder where results should be saved. If None, then the current working directory is used.
         """
         self.name = name
         self.flow_system = flow_system
         self.model: Optional[SystemModel] = None
         self.active_timesteps = active_timesteps
+        self.selected_scenarios = selected_scenarios
 
         self.durations = {'modeling': 0.0, 'solving': 0.0, 'saving': 0.0}
         self.folder = pathlib.Path.cwd() / 'results' if folder is None else pathlib.Path(folder)
@@ -74,47 +77,49 @@ def __init__(
     def main_results(self) -> Dict[str, Union[Scalar, Dict]]:
         from flixopt.features import InvestmentModel
 
-        return {
+        main_results = {
             'Objective': self.model.objective.value,
-            'Penalty': float(self.model.effects.penalty.total.solution.values),
+            'Penalty': self.model.effects.penalty.total.solution.values,
             'Effects': {
                 f'{effect.label} [{effect.unit}]': {
-                    'operation': float(effect.model.operation.total.solution.values),
-                    'invest': float(effect.model.invest.total.solution.values),
-                    'total': float(effect.model.total.solution.values),
+                    'operation': effect.model.operation.total.solution.values,
+                    'invest': effect.model.invest.total.solution.values,
+                    'total': effect.model.total.solution.values,
                 }
                 for effect in self.flow_system.effects
             },
             'Invest-Decisions': {
                 'Invested': {
-                    model.label_of_element: float(model.size.solution)
+                    model.label_of_element: model.size.solution
                     for component in self.flow_system.components.values()
                     for model in component.model.all_sub_models
-                    if isinstance(model, InvestmentModel) and float(model.size.solution) >= CONFIG.modeling.EPSILON
+                    if isinstance(model, InvestmentModel) and model.size.solution.max() >= CONFIG.modeling.EPSILON
                 },
                 'Not invested': {
-                    model.label_of_element: float(model.size.solution)
+                    model.label_of_element: model.size.solution
                     for component in self.flow_system.components.values()
                     for model in component.model.all_sub_models
-                    if isinstance(model, InvestmentModel) and float(model.size.solution) < CONFIG.modeling.EPSILON
+                    if isinstance(model, InvestmentModel) and model.size.solution.max() < CONFIG.modeling.EPSILON
                 },
             },
             'Buses with excess': [
                 {
                     bus.label_full: {
-                        'input': float(np.sum(bus.model.excess_input.solution.values)),
-                        'output': float(np.sum(bus.model.excess_output.solution.values)),
+                        'input': bus.model.excess_input.solution.sum('time'),
+                        'output': bus.model.excess_output.solution.sum('time'),
                     }
                 }
                 for bus in self.flow_system.buses.values()
                 if bus.with_excess
                 and (
-                    float(np.sum(bus.model.excess_input.solution.values)) > 1e-3
-                    or float(np.sum(bus.model.excess_output.solution.values)) > 1e-3
+                    bus.model.excess_input.solution.sum() > 1e-3
+                    or bus.model.excess_output.solution.sum() > 1e-3
                 )
             ],
         }
 
+        return utils.round_floats(main_results)
+
     @property
     def summary(self):
         return {
@@ -183,8 +188,8 @@ def solve(self, solver: _Solver, log_file: Optional[pathlib.Path] = None, log_ma
 
     def _activate_time_series(self):
         self.flow_system.transform_data()
-        self.flow_system.time_series_collection.activate_timesteps(
-            active_timesteps=self.active_timesteps,
+        self.flow_system.time_series_collection.set_selection(
+            timesteps=self.active_timesteps, scenarios=self.selected_scenarios
         )
 
 
@@ -217,6 +222,8 @@ def __init__(
                 list with indices, which should be used for calculation. If None, then all timesteps are used.
             folder: folder where results should be saved. If None, then the current working directory is used.
         """
+        if flow_system.time_series_collection.scenarios is not None:
+            raise ValueError('Aggregation is not supported for scenarios yet. Please use FullCalculation instead.')
         super().__init__(name, flow_system, active_timesteps, folder=folder)
         self.aggregation_parameters = aggregation_parameters
         self.components_to_clusterize = components_to_clusterize
@@ -272,9 +279,9 @@ def _perform_aggregation(self):
 
         # Aggregation - creation of aggregated timeseries:
         self.aggregation = Aggregation(
-            original_data=self.flow_system.time_series_collection.to_dataframe(
-                include_extra_timestep=False
-            ),  # Exclude last row (NaN)
+            original_data=self.flow_system.time_series_collection.as_dataset(
+                with_extra_timestep=False, with_constants=False
+            ).to_dataframe(),
             hours_per_time_step=float(dt_min),
             hours_per_period=self.aggregation_parameters.hours_per_period,
             nr_of_periods=self.aggregation_parameters.nr_of_periods,
@@ -286,9 +293,11 @@ def _perform_aggregation(self):
         self.aggregation.cluster()
         self.aggregation.plot(show=True, save=self.folder / 'aggregation.html')
         if self.aggregation_parameters.aggregate_data_and_fix_non_binary_vars:
-            self.flow_system.time_series_collection.insert_new_data(
-                self.aggregation.aggregated_data, include_extra_timestep=False
-            )
+            for col in self.aggregation.aggregated_data.columns:
+                data = self.aggregation.aggregated_data[col].values
+                if col in self.flow_system.time_series_collection._has_extra_timestep:
+                    data = np.append(data, data[-1])
+                self.flow_system.time_series_collection.update_time_series(col, data)
         self.durations['aggregation'] = round(timeit.default_timer() - t_start_agg, 2)
 
 
@@ -327,8 +336,8 @@ def __init__(
         self.nr_of_previous_values = nr_of_previous_values
         self.sub_calculations: List[FullCalculation] = []
 
-        self.all_timesteps = self.flow_system.time_series_collection.all_timesteps
-        self.all_timesteps_extra = self.flow_system.time_series_collection.all_timesteps_extra
+        self.all_timesteps = self.flow_system.time_series_collection._full_timesteps
+        self.all_timesteps_extra = self.flow_system.time_series_collection._full_timesteps_extra
 
         self.segment_names = [
             f'Segment_{i + 1}' for i in range(math.ceil(len(self.all_timesteps) / self.timesteps_per_segment))