diff --git a/.github/workflows/python-app.yaml b/.github/workflows/python-app.yaml index 6f214e88d..068995579 100644 --- a/.github/workflows/python-app.yaml +++ b/.github/workflows/python-app.yaml @@ -342,7 +342,7 @@ jobs: # Wait and retry while PyPI indexes the package INSTALL_SUCCESS=false - for d in 10 20 40 80 120; do + for d in 10 20 40 80 120 240 480; do sleep "$d" echo "Attempting to install $PACKAGE_NAME==$VERSION from PyPI (retry after ${d}s)..." diff --git a/CHANGELOG.md b/CHANGELOG.md index f6097e7a3..56a827078 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -52,7 +52,7 @@ Common use cases are: The weighted sum of the total objective effect of each scenario is used as the objective of the optimization. -#### Improved Data handling: I/O, resampling and more through xarray +### Improved Data handling: I/O, resampling and more through xarray * IO for all Interfaces and the FlowSystem with round-trip serialization support * NetCDF export/import capabilities for all Interface objects and FlowSystem * JSON export for documentation purposes @@ -105,13 +105,13 @@ The weighted sum of the total objective effect of each scenario is used as the o * Better type consistency across all framework components ### Known issues -* IO for single Interfaces/Elements to Datasets might not work properly if the Interface/Element is not part of a fully transformed and connected FlowSystem. This arises from Numeric Data not being stored as xr.DataArray by the user. To avoid this, always use the `to_dataset()` on Elements inside a FlowSystem that's connected and transformed. +* IO for single Interfaces/Elements to Datasets might not work properly if the Interface/Element is not part of a fully transformed and connected FlowSystem. This arises from numeric data not being stored as xr.DataArray by the user. To avoid this, always use `to_dataset()` on Elements inside a FlowSystem that's connected and transformed. ### *Development* * **BREAKING**: Calculation.do_modeling() now returns the Calculation object instead of its linopy.Model * FlowSystem data management simplified - removed `time_series_collection` pattern in favor of direct timestep properties * Change modeling hierarchy to allow for more flexibility in future development. This leads to minimal changes in the access and creation of Submodels and their variables. -* Added new module `.modeling`that contains Modelling primitives and utilities +* Added new module `.modeling` that contains modeling primitives and utilities * Clearer separation between the main Model and "Submodels" * Improved access to the Submodels and their variables, constraints and submodels * Added `__repr__()` for Submodels to easily inspect its content diff --git a/flixopt/elements.py b/flixopt/elements.py index 04bc7a599..8335a00a8 100644 --- a/flixopt/elements.py +++ b/flixopt/elements.py @@ -829,7 +829,7 @@ def _do_modeling(self): ) if self.element.prevent_simultaneous_flows: - # Simultanious Useage --> Only One FLow is On at a time, but needs a Binary for every flow + # Simultaneous Usage --> Only One Flow is On at a time, but needs a Binary for every flow ModelingPrimitives.mutual_exclusivity_constraint( self, binary_variables=[flow.submodel.on_off.on for flow in self.element.prevent_simultaneous_flows], diff --git a/flixopt/features.py b/flixopt/features.py index 51a4832b7..f88e3b987 100644 --- a/flixopt/features.py +++ b/flixopt/features.py @@ -10,8 +10,10 @@ import linopy import numpy as np +import xarray as xr -from .modeling import BoundingPatterns, ModelingPrimitives, ModelingUtilities +from .config import CONFIG +from .modeling import BoundingPatterns, ModelingPrimitives, ModelingUtilities, ModelingUtilitiesAbstract from .structure import FlowSystemModel, Submodel if TYPE_CHECKING: @@ -49,11 +51,14 @@ def __init__( def _do_modeling(self): super()._do_modeling() - self._create_variables_and_constraints() + if self._model.flow_system.years is None: + self._create_variables_and_constraints_without_years() + else: + self._create_variables_and_constraints_with_years() self._add_effects() - def _create_variables_and_constraints(self): - size_min, size_max = (self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size) + def _create_variables_and_constraints_without_years(self): + size_min, size_max = self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size self.add_variables( short_name='size', lower=0 if self.parameters.optional else size_min, @@ -65,29 +70,177 @@ def _create_variables_and_constraints(self): self.add_variables( binary=True, coords=self._model.get_coords(['year', 'scenario']), - short_name='is_invested', + short_name='is_invested', # TODO: rename to invested ) BoundingPatterns.bounds_with_state( self, variable=self.size, variable_state=self.is_invested, - bounds=(self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size), + bounds=(size_min, size_max), ) - if self.parameters.piecewise_effects: - self.piecewise_effects = self.add_submodels( - PiecewiseEffectsModel( - model=self._model, - label_of_element=self.label_of_element, - label_of_model=f'{self.label_of_element}|PiecewiseEffects', - piecewise_origin=(self.size.name, self.parameters.piecewise_effects.piecewise_origin), - piecewise_shares=self.parameters.piecewise_effects.piecewise_shares, - zero_point=self.is_invested, - ), - short_name='segments', + ######################################################################## + # Global investment decision variables (same as with_years version) + previous_size = self.parameters.previous_size if self.parameters.previous_size is not None else xr.DataArray(0) + + self.add_variables( + binary=True, + coords=self._model.get_coords(['year', 'scenario']), + short_name='size|investment_occurs', + ) + self.add_variables( + binary=True, + coords=self._model.get_coords(['year', 'scenario']), + short_name='size|decommissioning_occurs', + ) + + BoundingPatterns.state_transition_bounds( + self, + state_variable=self.is_invested, + switch_on=self.investment_occurs, + switch_off=self.decommissioning_occurs, + name=self.is_invested.name, + previous_state=ModelingUtilitiesAbstract.to_binary(values=previous_size, epsilon=CONFIG.modeling.EPSILON), + coord='year', + ) + + ######################################################################## + # Size change variables (same as with_years version) + self.add_variables( + coords=self._model.get_coords(['year', 'scenario']), + short_name='size|increase', + lower=0, + upper=size_max, + ) + self.add_variables( + coords=self._model.get_coords(['year', 'scenario']), + short_name='size|decrease', + lower=0, + upper=size_max, + ) + + BoundingPatterns.link_changes_to_level_with_binaries( + self, + level_variable=self.size, + increase_variable=self.size_increase, + decrease_variable=self.size_decrease, + increase_binary=self.investment_occurs, + decrease_binary=self.decommissioning_occurs, + name=f'{self.label_of_element}|size|changes', + max_change=size_max, + initial_level=previous_size, + coord='year', + ) + + def _create_variables_and_constraints_with_years(self): + size_min, size_max = self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size + + # Always create size variable + self.add_variables( + short_name='size', + lower=0, + upper=size_max, + coords=self._model.get_coords(['year', 'scenario']), + ) + + # Always create is_invested variable (yearly) + if self.parameters.optional: + self.add_variables( + binary=True, + coords=self._model.get_coords(['year', 'scenario']), + short_name='is_invested', + ) + else: + # When not optional, fix is_invested to 1 for all years + self.add_variables( + coords=self._model.get_coords(['year', 'scenario']), + short_name='is_invested', + lower=1, + upper=1, + ) + + # Apply bounds with state for both cases + BoundingPatterns.bounds_with_state( + self, + variable=self.size, + variable_state=self.is_invested, + bounds=(size_min, size_max), + ) + + ######################################################################## + # Global investment decision variables + previous_size = self.parameters.previous_size if self.parameters.previous_size is not None else xr.DataArray(0) + + self.add_variables( + binary=True, + coords=self._model.get_coords(['year', 'scenario']), + short_name='size|investment_occurs', + ) + self.add_variables( + binary=True, + coords=self._model.get_coords(['year', 'scenario']), + short_name='size|decommissioning_occurs', + ) + + BoundingPatterns.state_transition_bounds( + self, + state_variable=self.is_invested, + switch_on=self.investment_occurs, + switch_off=self.decommissioning_occurs, + name=self.is_invested.name, + previous_state=ModelingUtilitiesAbstract.to_binary(values=previous_size, epsilon=CONFIG.modeling.EPSILON), + coord='year', + ) + + # Global investment constraints + if self.parameters.optional: + # Can invest at most once + self.add_constraints( + self.investment_occurs.sum('year') <= 1, + short_name='investment_occurs|once', + ) + else: + # Must invest exactly once + self.add_constraints( + self.investment_occurs.sum('year') == 1, + short_name='investment_occurs|once', ) + # Can decommission at most once + self.add_constraints( + self.decommissioning_occurs.sum('year') <= 1, + short_name='decommissioning_occurs|once', + ) + + ######################################################################## + # Size change variables + self.add_variables( + coords=self._model.get_coords(['year', 'scenario']), + short_name='size|increase', + lower=0, + upper=size_max, + ) + self.add_variables( + coords=self._model.get_coords(['year', 'scenario']), + short_name='size|decrease', + lower=0, + upper=size_max, + ) + + BoundingPatterns.link_changes_to_level_with_binaries( + self, + level_variable=self.size, + increase_variable=self.size_increase, + decrease_variable=self.size_decrease, + increase_binary=self.investment_occurs, + decrease_binary=self.decommissioning_occurs, + name=f'{self.label_of_element}|size|changes', + max_change=size_max, + initial_level=previous_size, + coord='year', + ) + def _add_effects(self): """Add investment effects""" if self.parameters.fix_effects: @@ -117,6 +270,19 @@ def _add_effects(self): target='invest', ) + if self.parameters.piecewise_effects: + self.piecewise_effects = self.add_submodels( + PiecewiseEffectsModel( + model=self._model, + label_of_element=self.label_of_element, + label_of_model=f'{self.label_of_element}|PiecewiseEffects', + piecewise_origin=(self.size.name, self.parameters.piecewise_effects.piecewise_origin), + piecewise_shares=self.parameters.piecewise_effects.piecewise_shares, + zero_point=self.is_invested, + ), + short_name='segments', + ) + @property def size(self) -> linopy.Variable: """Investment size variable""" @@ -129,6 +295,26 @@ def is_invested(self) -> linopy.Variable | None: return None return self._variables['is_invested'] + @property + def investment_occurs(self) -> linopy.Variable: + """Binary increase decision variable""" + return self._variables['size|investment_occurs'] + + @property + def decommissioning_occurs(self) -> linopy.Variable: + """Binary decrease decision variable""" + return self._variables['size|decommissioning_occurs'] + + @property + def size_decrease(self) -> linopy.Variable: + """Binary decrease decision variable""" + return self._variables['size|decrease'] + + @property + def size_increase(self) -> linopy.Variable: + """Binary increase decision variable""" + return self._variables['size|increase'] + class OnOffModel(Submodel): """OnOff model using factory patterns""" diff --git a/flixopt/interface.py b/flixopt/interface.py index da456efee..318ee6ba3 100644 --- a/flixopt/interface.py +++ b/flixopt/interface.py @@ -8,6 +8,8 @@ import logging from typing import TYPE_CHECKING, Literal, Optional +import xarray as xr + from .config import CONFIG from .structure import Interface, register_class_for_io @@ -705,6 +707,7 @@ class InvestParameters(Interface): divest_effects: Costs incurred if the investment is NOT made, such as demolition of existing equipment, contractual penalties, or lost opportunities. Dictionary mapping effect names to values. + previous_size: Initial size of the investment. Only relevant in multi-period investments. Cost Annualization Requirements: All cost values must be properly weighted to match the optimization model's time horizon. @@ -861,6 +864,7 @@ def __init__( specific_effects: NonTemporalEffectsUser | None = None, # costs per Flow-Unit/Storage-Size/... piecewise_effects: PiecewiseEffects | None = None, divest_effects: NonTemporalEffectsUser | None = None, + previous_size: NonTemporalDataUser | None = None, ): self.fix_effects: NonTemporalEffectsUser = fix_effects or {} self.divest_effects: NonTemporalEffectsUser = divest_effects or {} @@ -870,6 +874,7 @@ def __init__( self.piecewise_effects = piecewise_effects self.minimum_size = minimum_size if minimum_size is not None else CONFIG.modeling.EPSILON self.maximum_size = maximum_size if maximum_size is not None else CONFIG.modeling.BIG # default maximum + self.previous_size = previous_size def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None: self.fix_effects = flow_system.fit_effects_to_model_coords( @@ -904,6 +909,9 @@ def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None self.fixed_size = flow_system.fit_to_model_coords( f'{name_prefix}|fixed_size', self.fixed_size, dims=['year', 'scenario'] ) + self.previous_size = flow_system.fit_to_model_coords( + f'{name_prefix}|previous_size', self.previous_size, dims=['scenario'] + ) @property def minimum_or_fixed_size(self) -> NonTemporalData: diff --git a/tests/test_models.py b/tests/test_models.py new file mode 100644 index 000000000..3fc38c1f1 --- /dev/null +++ b/tests/test_models.py @@ -0,0 +1,231 @@ +from typing import Union + +import linopy +import numpy as np +import pandas as pd +import pytest +import xarray as xr + +import flixopt as fx + +from .conftest import ( + Buses, + Effects, + LoadProfiles, + Sinks, + Sources, + assert_conequal, + assert_sets_equal, + assert_var_equal, + create_linopy_model, +) + + +def calculate_annual_payment(total_cost: float, remaining_years: int, discount_rate: float) -> float: + """Calculate annualized payment for given remaining years. + + Args: + total_cost: Total cost to be annualized. + remaining_years: Number of remaining years. + discount_rate: Discount rate for annualization. + + Returns: + Annual payment amount. + """ + if remaining_years == 1: + return total_cost + + return ( + total_cost + * (discount_rate * (1 + discount_rate) ** remaining_years) + / ((1 + discount_rate) ** remaining_years - 1) + ) + + +def create_annualized_effects( + year_of_investments: range | list | pd.Index, + all_years: range | list | pd.Index, + total_cost: float, + discount_rate: float, + horizon_end: int, + extra_dim: str = 'year_of_investment', +) -> xr.DataArray: + """Create a 2D effects array for annualized costs. + + Creates an array where investing in year Y results in annualized costs + applied to years Y through horizon_end. + + Args: + year_of_investments: Years when investment decisions can be made. + all_years: All years in the model (for the 'year' dimension). + total_cost: Total upfront cost to be annualized. + discount_rate: Discount rate for annualization calculation. + horizon_end: Last year when effects apply. + extra_dim: Name for the investment year dimension. + + Returns: + xr.DataArray with dimensions [extra_dim, 'year'] containing annualized costs. + """ + + # Convert to lists for easier iteration + year_of_investments_list = list(year_of_investments) + all_years_list = list(all_years) + + # Initialize cost matrix + cost_matrix = np.zeros((len(year_of_investments_list), len(all_years_list))) + + # Fill matrix with annualized costs + for i, year_of_investment in enumerate(year_of_investments_list): + remaining_years = horizon_end - year_of_investment + 1 + if remaining_years > 0: + annual_cost = calculate_annual_payment(total_cost, remaining_years, discount_rate) + + # Apply cost to years from year_of_investment through horizon_end + for j, cost_year in enumerate(all_years_list): + if year_of_investment <= cost_year <= horizon_end: + cost_matrix[i, j] = annual_cost + + return xr.DataArray( + cost_matrix, coords={extra_dim: year_of_investments_list, 'year': all_years_list}, dims=[extra_dim, 'year'] + ) + + +@pytest.fixture +def flow_system() -> fx.FlowSystem: + """Create basic elements for component testing with coordinate parametrization.""" + years = pd.Index([2020, 2021, 2022, 2023, 2024, 2026, 2028, 2030], name='year') + timesteps = pd.date_range('2020-01-01', periods=24, freq='h', name='time') + flow_system = fx.FlowSystem(timesteps=timesteps, years=years) + + thermal_load = LoadProfiles.random_thermal(len(timesteps)) + p_el = LoadProfiles.random_electrical(len(timesteps)) + + costs = Effects.costs() + heat_load = Sinks.heat_load(thermal_load) + gas_source = Sources.gas_with_costs() + electricity_sink = Sinks.electricity_feed_in(p_el) + + flow_system.add_elements(*Buses.defaults()) + flow_system.buses['Fernwärme'].excess_penalty_per_flow_hour = 0 + flow_system.add_elements(costs, heat_load, gas_source, electricity_sink) + + return flow_system + + +class TestYearAwareInvestParameters: + """Test the YearAwareInvestParameters interface.""" + + def test_basic_initialization(self): + """Test basic parameter initialization.""" + params = fx.YearAwareInvestParameters( + minimum_size=10, + maximum_size=100, + ) + + assert params.minimum_size == 10 + assert params.maximum_size == 100 + assert params.fixed_size is None + assert not params.allow_divestment + assert params.fixed_year_of_investment is None + assert params.fixed_year_of_decommissioning is None + assert params.fixed_duration is None + + def test_fixed_size_initialization(self): + """Test initialization with fixed size.""" + params = fx.YearAwareInvestParameters(fixed_size=50) + + assert params.minimum_or_fixed_size == 50 + assert params.maximum_or_fixed_size == 50 + assert params.is_fixed_size + + def test_timing_constraints_initialization(self): + """Test initialization with various timing constraints.""" + params = fx.YearAwareInvestParameters( + fixed_year_of_investment=2, + minimum_duration=3, + maximum_duration=5, + earliest_year_of_decommissioning=4, + ) + + assert params.fixed_year_of_investment == 2 + assert params.minimum_duration == 3 + assert params.maximum_duration == 5 + assert params.earliest_year_of_decommissioning == 4 + + def test_effects_initialization(self): + """Test initialization with effects.""" + params = fx.YearAwareInvestParameters( + effects_of_investment={'costs': 1000}, + effects_of_investment_per_size={'costs': 100}, + allow_divestment=True, + effects_of_divestment={'costs': 500}, + effects_of_divestment_per_size={'costs': 50}, + ) + + assert params.effects_of_investment == {'costs': 1000} + assert params.effects_of_investment_per_size == {'costs': 100} + assert params.allow_divestment + assert params.effects_of_divestment == {'costs': 500} + assert params.effects_of_divestment_per_size == {'costs': 50} + + def test_property_methods(self): + """Test property methods.""" + # Test with fixed size + params_fixed = fx.YearAwareInvestParameters(fixed_size=50) + assert params_fixed.minimum_or_fixed_size == 50 + assert params_fixed.maximum_or_fixed_size == 50 + assert params_fixed.is_fixed_size + + # Test with min/max size + params_range = fx.YearAwareInvestParameters(minimum_size=10, maximum_size=100) + assert params_range.minimum_or_fixed_size == 10 + assert params_range.maximum_or_fixed_size == 100 + assert not params_range.is_fixed_size + + +class TestYearAwareInvestmentModelDirect: + """Test the YearAwareInvestmentModel class directly with linopy.""" + + def test_flow_invest_new(self, flow_system): + da = xr.DataArray( + [10] * 8, + coords=(flow_system.years_of_investment,), + ).expand_dims(year=flow_system.years) + da = da.where(da.year == da.year_of_investment).fillna(0) + + flow = fx.Flow( + 'Wärme', + bus='Fernwärme', + size=fx.InvestTimingParameters( + # year_of_decommissioning=2030, + minimum_lifetime=2, + maximum_lifetime=3, + minimum_size=9, + maximum_size=10, + specific_effects=xr.DataArray( + [25, 30, 35, 40, 45, 50, 55, 60], + coords=(flow_system.years,), + ) + * -0, + # fix_effects=-2e3, + specific_effects_by_investment_year=-1 * da, + ), + relative_maximum=np.linspace(0.5, 1, flow_system.timesteps.size), + ) + + flow_system.add_elements(fx.Source('Source', source=flow)) + calculation = fx.FullCalculation('GenericName', flow_system) + calculation.do_modeling() + # calculation.model.add_constraints(calculation.model['Source(Wärme)|is_invested'].sel(year=2022) == 1) + calculation.solve(fx.solvers.GurobiSolver(0, 60)) + + ds = calculation.results['Source'].solution + filtered_ds_year = ds[[v for v in ds.data_vars if ds[v].dims == ('year',)]] + print(filtered_ds_year.round(0).to_pandas().T) + + filtered_ds_scalar = ds[[v for v in ds.data_vars if ds[v].dims == tuple()]] + print(filtered_ds_scalar.round(0).to_pandas().T) + + print(calculation.results.solution['costs(invest)|total'].to_pandas()) + + print('##') diff --git a/tests/test_storage.py b/tests/test_storage.py index 5b84214b8..526fa04a9 100644 --- a/tests/test_storage.py +++ b/tests/test_storage.py @@ -481,6 +481,6 @@ def test_investment_parameters( if 'InvestStorage|is_invested' in model.variables: var = model.variables['InvestStorage|is_invested'] # Check if the lower and upper bounds are both 1 - assert var.upper == 1 and var.lower == 1, ( + assert np.all(var.upper == 1) and np.all(var.lower == 1), ( 'is_invested variable should be fixed to 1 when optional=False' )