From c1fa5d6d0d534b9fa3c26f6f1318669b179aa23b Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 18:12:40 +0100 Subject: [PATCH 01/19] Updated CHangelog and migration guide --- CHANGELOG.md | 2 +- docs/user-guide/migration-guide-v5.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d104dc819..7d556bf4a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -393,7 +393,7 @@ results = calc.results # Returns SegmentedResults ``` **Planned for future release:** -- `flow_system.optimize.rolling(solver, ...)` - Rolling horizon optimization to replace `SegmentedOptimization` +- `flow_system.optimize.rolling_horizon(solver, horizon, overlap, ...)` - Rolling horizon optimization to replace `SegmentedOptimization` ### Migration Checklist diff --git a/docs/user-guide/migration-guide-v5.md b/docs/user-guide/migration-guide-v5.md index 5c19761e0..044f26897 100644 --- a/docs/user-guide/migration-guide-v5.md +++ b/docs/user-guide/migration-guide-v5.md @@ -327,7 +327,7 @@ Clustered optimization uses the new transform accessor: ### Segmented Optimization (Not Yet Migrated) -Segmented optimization still uses the class-based API. A new `optimize.rolling()` method is planned for a future release. +Segmented optimization still uses the class-based API. A new `optimize.rolling_horizon()` method is planned for a future release. ```python # Still use the class-based API (unchanged from v4.x) From 57f5ea1535094dcfe8104957d6d2e077cb08f4e1 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 18:23:43 +0100 Subject: [PATCH 02/19] The rolling_horizon method is now implemented on OptimizeAccessor: API: segments = flow_system.optimize.rolling_horizon( solver, horizon=100, # Timesteps per segment overlap=0, # Lookahead timesteps nr_of_previous_values=1, # State transfer depth ) flow_system.solution # Combined solution (overlaps trimmed) segments # List[FlowSystem] with individual segment solutions Key features implemented: 1. Segment calculation - Divides timesteps into overlapping windows 2. State transfer - Transfers storage charge states and flow rates between segments 3. Solution combining - Trims overlaps and concatenates time-indexed variables 4. Effect recalculation - Correctly sums costs from combined per-timestep values 5. Investment validation - Raises error if InvestParameters are used 6. Progress bar - Shows solving progress with tqdm Location: flixopt/optimize_accessor.py:91-399 --- flixopt/optimize_accessor.py | 319 +++++++++++++++++++++++++++++++++-- 1 file changed, 309 insertions(+), 10 deletions(-) diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py index 5428cd855..7f27172e2 100644 --- a/flixopt/optimize_accessor.py +++ b/flixopt/optimize_accessor.py @@ -7,12 +7,21 @@ from __future__ import annotations +import logging +import sys from typing import TYPE_CHECKING +import xarray as xr +from tqdm import tqdm + +from .config import CONFIG + if TYPE_CHECKING: from .flow_system import FlowSystem from .solvers import _Solver +logger = logging.getLogger('flixopt') + class OptimizeAccessor: """ @@ -79,13 +88,303 @@ def __call__(self, solver: _Solver, normalize_weights: bool = True) -> FlowSyste self._fs.solve(solver) return self._fs - # Future methods can be added here: - # - # def clustered(self, solver: _Solver, aggregation: AggregationParameters, - # normalize_weights: bool = True) -> FlowSystem: - # """Clustered optimization with time aggregation.""" - # ... - # - # def mga(self, solver: _Solver, alternatives: int = 5) -> FlowSystem: - # """Modeling to Generate Alternatives.""" - # ... + def rolling_horizon( + self, + solver: _Solver, + horizon: int = 100, + overlap: int = 0, + nr_of_previous_values: int = 1, + ) -> list[FlowSystem]: + """ + Solve the optimization using a rolling horizon approach. + + Divides the time horizon into overlapping segments that are solved sequentially. + Each segment uses final values from the previous segment as initial conditions, + ensuring dynamic continuity across the solution. The combined solution is stored + on the original FlowSystem. + + This approach is useful for: + - Large-scale problems that exceed memory limits + - Annual planning with seasonal variations + - Operational planning with limited foresight + + Args: + solver: The solver to use (e.g., HighsSolver, GurobiSolver). + horizon: Number of timesteps in each segment (excluding overlap). + Must be > 2. Larger values provide better optimization at the cost + of memory and computation time. Default: 100. + overlap: Number of additional timesteps added to each segment for lookahead. + Improves storage optimization by providing foresight. Higher values + improve solution quality but increase computational cost. Default: 0. + nr_of_previous_values: Number of previous timestep values to transfer between + segments for initialization (e.g., for uptime/downtime tracking). Default: 1. + + Returns: + List of segment FlowSystems, each with their individual solution. + The combined solution (with overlaps trimmed) is stored on the original FlowSystem. + + Raises: + ValueError: If horizon <= 2 or overlap < 0. + ValueError: If horizon + overlap > total timesteps. + ValueError: If InvestParameters are used (not supported in rolling horizon). + + Examples: + Basic rolling horizon optimization: + + >>> segments = flow_system.optimize.rolling_horizon( + ... solver, + ... horizon=168, # Weekly segments + ... overlap=24, # 1-day lookahead + ... ) + >>> print(flow_system.solution) # Combined result + + Inspect individual segments: + + >>> for i, seg in enumerate(segments): + ... print(f'Segment {i}: {seg.solution["costs(total)"].item():.2f}') + + Note: + - InvestParameters are not supported as investment decisions require + full-horizon optimization. + - Global constraints (flow_hours_max, etc.) may produce suboptimal results + as they cannot be enforced globally across segments. + - Storage optimization may be suboptimal compared to full-horizon solutions + due to limited foresight in each segment. + """ + + # Validation + if horizon <= 2: + raise ValueError('horizon must be greater than 2 to avoid internal side effects.') + if overlap < 0: + raise ValueError('overlap must be non-negative.') + if nr_of_previous_values < 0: + raise ValueError('nr_of_previous_values must be non-negative.') + if nr_of_previous_values > horizon: + raise ValueError('nr_of_previous_values cannot exceed horizon.') + + total_timesteps = len(self._fs.timesteps) + horizon_with_overlap = horizon + overlap + + if horizon_with_overlap > total_timesteps: + raise ValueError( + f'horizon + overlap ({horizon_with_overlap}) cannot exceed total timesteps ({total_timesteps}).' + ) + + # Ensure flow system is connected + if not self._fs.connected_and_transformed: + self._fs.connect_and_transform() + + # Calculate segment indices + segment_indices = self._calculate_segment_indices(total_timesteps, horizon, overlap) + n_segments = len(segment_indices) + + logger.info(f'{"":#^80}') + logger.info(f'{" Rolling Horizon Optimization ":#^80}') + logger.info(f'Segments: {n_segments}, Horizon: {horizon}, Overlap: {overlap}') + + # Store original initial values for restoration and state transfer + original_initial_values = self._store_initial_values() + + # Create and solve segments + segment_flow_systems: list[FlowSystem] = [] + + progress_bar = tqdm( + enumerate(segment_indices), + total=n_segments, + desc='Solving segments', + unit='segment', + file=sys.stdout, + disable=not CONFIG.Solving.log_to_console, + ) + + try: + for i, (start_idx, end_idx) in progress_bar: + progress_bar.set_description(f'Segment {i + 1}/{n_segments} (timesteps {start_idx}-{end_idx})') + + # Create segment FlowSystem + segment_fs = self._fs.transform.isel(time=slice(start_idx, end_idx)) + + # Transfer state from previous segment + if i > 0 and nr_of_previous_values > 0: + self._transfer_state( + source_fs=segment_flow_systems[i - 1], + target_fs=segment_fs, + horizon=horizon, + nr_of_previous_values=nr_of_previous_values, + ) + + # Build and solve + segment_fs.build_model() + + # Check for investments (only on first segment) + if i == 0: + self._check_no_investments(segment_fs) + + segment_fs.solve(solver) + segment_flow_systems.append(segment_fs) + + finally: + progress_bar.close() + + # Combine solutions and store on original FlowSystem + combined_solution = self._combine_solutions(segment_flow_systems, horizon) + self._fs._solution = combined_solution + + # Restore original initial values + self._restore_initial_values(original_initial_values) + + logger.info(f'Rolling horizon optimization completed: {n_segments} segments solved.') + + return segment_flow_systems + + def _calculate_segment_indices(self, total_timesteps: int, horizon: int, overlap: int) -> list[tuple[int, int]]: + """Calculate start and end indices for each segment.""" + segments = [] + start = 0 + while start < total_timesteps: + end = min(start + horizon + overlap, total_timesteps) + segments.append((start, end)) + start += horizon # Move by horizon (not horizon + overlap) + if end == total_timesteps: + break + return segments + + def _store_initial_values(self) -> dict: + """Store original initial values for later restoration.""" + from .components import Storage + + values = {} + for flow in self._fs.flows.values(): + values[f'flow|{flow.label_full}'] = flow.previous_flow_rate + + for comp in self._fs.components.values(): + if isinstance(comp, Storage): + values[f'storage|{comp.label_full}'] = comp.initial_charge_state + + return values + + def _restore_initial_values(self, values: dict) -> None: + """Restore original initial values after rolling horizon.""" + from .components import Storage + + for flow in self._fs.flows.values(): + key = f'flow|{flow.label_full}' + if key in values: + flow.previous_flow_rate = values[key] + + for comp in self._fs.components.values(): + if isinstance(comp, Storage): + key = f'storage|{comp.label_full}' + if key in values: + comp.initial_charge_state = values[key] + + def _transfer_state( + self, + source_fs: FlowSystem, + target_fs: FlowSystem, + horizon: int, + nr_of_previous_values: int, + ) -> None: + """Transfer final state from source segment to target segment.""" + + from .components import Storage + + # Transfer flow rates (for uptime/downtime tracking) + for source_flow in source_fs.flows.values(): + target_flow = target_fs.flows.get(source_flow.label_full) + if target_flow is None: + continue + + # Get last nr_of_previous_values from source solution + flow_rate_var = f'{source_flow.label_full}|flow_rate' + if flow_rate_var in source_fs.solution: + # Select from the non-overlap portion (first 'horizon' timesteps) + values = ( + source_fs.solution[flow_rate_var].isel(time=slice(horizon - nr_of_previous_values, horizon)).values + ) + target_flow.previous_flow_rate = values if len(values) > 1 else values.item() + + # Transfer storage charge states + for source_comp in source_fs.components.values(): + if not isinstance(source_comp, Storage): + continue + + target_comp = target_fs.components.get(source_comp.label_full) + if target_comp is None or not isinstance(target_comp, Storage): + continue + + charge_var = f'{source_comp.label_full}|charge_state' + if charge_var in source_fs.solution: + # Get charge state at the end of the non-overlap portion + # Use horizon index (0-indexed, so horizon-1 is last non-overlap) + charge_state = source_fs.solution[charge_var].isel(time=horizon - 1).values.item() + target_comp.initial_charge_state = charge_state + + def _check_no_investments(self, segment_fs: FlowSystem) -> None: + """Check that no InvestParameters are used (not supported in rolling horizon).""" + from .features import InvestmentModel + + invest_elements = [] + for component in segment_fs.components.values(): + for model in component.submodel.all_submodels: + if isinstance(model, InvestmentModel): + invest_elements.append(model.label_full) + + if invest_elements: + raise ValueError( + f'InvestParameters are not supported in rolling horizon optimization. ' + f'Found InvestmentModels: {invest_elements}. ' + f'Use standard optimize() for problems with investments.' + ) + + def _combine_solutions(self, segment_flow_systems: list[FlowSystem], horizon: int) -> xr.Dataset: + """Combine segment solutions, trimming overlaps.""" + if not segment_flow_systems: + raise ValueError('No segments to combine.') + + # Get all variable names from first segment + var_names = list(segment_flow_systems[0].solution.data_vars) + + # Identify effect names for special handling + effect_labels = {e.label for e in self._fs.effects.values()} + + combined_vars = {} + for var_name in var_names: + arrays = [] + for i, seg_fs in enumerate(segment_flow_systems): + da = seg_fs.solution[var_name] + + # Check if this variable has a time dimension + if 'time' in da.dims: + if i < len(segment_flow_systems) - 1: + # Not the last segment: trim to horizon (exclude overlap) + da = da.isel(time=slice(None, horizon)) + # Last segment: keep all timesteps + arrays.append(da) + + # Concatenate along time if time dimension exists + if 'time' in segment_flow_systems[0].solution[var_name].dims: + combined_vars[var_name] = xr.concat(arrays, dim='time') + else: + # For non-time scalars, just take the last value for now + # Effect totals will be recalculated below + combined_vars[var_name] = arrays[-1] + + # Recalculate effect totals from combined per-timestep data + for effect_label in effect_labels: + per_timestep_key = f'{effect_label}(temporal)|per_timestep' + temporal_key = f'{effect_label}(temporal)' + total_key = effect_label + + if per_timestep_key in combined_vars: + # Recalculate temporal sum from per-timestep values + combined_vars[temporal_key] = combined_vars[per_timestep_key].sum() + + # Recalculate total (temporal + periodic) + periodic_key = f'{effect_label}(periodic)' + if periodic_key in combined_vars: + combined_vars[total_key] = combined_vars[temporal_key] + combined_vars[periodic_key] + else: + combined_vars[total_key] = combined_vars[temporal_key] + + return xr.Dataset(combined_vars) From 58cce188c6608915b0b77c1a11a7e879bdbeac6a Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 18:59:02 +0100 Subject: [PATCH 03/19] Summary: Rolling Horizon Optimization The optimize.rolling_horizon() method has been implemented in flixopt/optimize_accessor.py with the following features: API segments = flow_system.optimize.rolling_horizon( solver, horizon=100, # Timesteps per segment (excluding overlap) overlap=0, # Additional lookahead timesteps nr_of_previous_values=1, # For uptime/downtime tracking ) Key Features 1. Segment creation: Divides time horizon into overlapping segments using transform.isel() 2. State transfer: Transfers storage charge states and previous flow rates between segments 3. Solution combining: Trims overlaps from segment solutions and concatenates 4. Effect recomputation: Computes total effect values from combined per-timestep values (no re-solve needed) 5. Investment check: Raises error if InvestParameters are detected (not supported) Implementation Details - Uses _calculate_segment_indices() to compute start/end indices - Uses _transfer_state() to pass storage charge states and flow rates to next segment - Uses _combine_solutions() to merge segment solutions and recompute effect totals - Progress bar shows segment progress via tqdm --- flixopt/optimize_accessor.py | 153 +++++++++++++++++------------------ 1 file changed, 72 insertions(+), 81 deletions(-) diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py index 7f27172e2..e59c07df5 100644 --- a/flixopt/optimize_accessor.py +++ b/flixopt/optimize_accessor.py @@ -37,10 +37,10 @@ class OptimizeAccessor: >>> flow_system.optimize(solver) >>> print(flow_system.solution) - Future specialized modes: + Rolling horizon optimization: - >>> flow_system.optimize.clustered(solver, aggregation=params) - >>> flow_system.optimize.mga(solver, alternatives=5) + >>> segments = flow_system.optimize.rolling_horizon(solver, horizon=168) + >>> print(flow_system.solution) # Combined result """ def __init__(self, flow_system: FlowSystem) -> None: @@ -141,7 +141,7 @@ def rolling_horizon( Inspect individual segments: >>> for i, seg in enumerate(segments): - ... print(f'Segment {i}: {seg.solution["costs(total)"].item():.2f}') + ... print(f'Segment {i}: {seg.solution["costs"].item():.2f}') Note: - InvestParameters are not supported as investment decisions require @@ -182,9 +182,6 @@ def rolling_horizon( logger.info(f'{" Rolling Horizon Optimization ":#^80}') logger.info(f'Segments: {n_segments}, Horizon: {horizon}, Overlap: {overlap}') - # Store original initial values for restoration and state transfer - original_initial_values = self._store_initial_values() - # Create and solve segments segment_flow_systems: list[FlowSystem] = [] @@ -226,12 +223,9 @@ def rolling_horizon( finally: progress_bar.close() - # Combine solutions and store on original FlowSystem - combined_solution = self._combine_solutions(segment_flow_systems, horizon) - self._fs._solution = combined_solution - - # Restore original initial values - self._restore_initial_values(original_initial_values) + # Combine segment solutions + logger.info('Combining segment solutions...') + self._finalize_solution(segment_flow_systems, horizon) logger.info(f'Rolling horizon optimization completed: {n_segments} segments solved.') @@ -249,35 +243,6 @@ def _calculate_segment_indices(self, total_timesteps: int, horizon: int, overlap break return segments - def _store_initial_values(self) -> dict: - """Store original initial values for later restoration.""" - from .components import Storage - - values = {} - for flow in self._fs.flows.values(): - values[f'flow|{flow.label_full}'] = flow.previous_flow_rate - - for comp in self._fs.components.values(): - if isinstance(comp, Storage): - values[f'storage|{comp.label_full}'] = comp.initial_charge_state - - return values - - def _restore_initial_values(self, values: dict) -> None: - """Restore original initial values after rolling horizon.""" - from .components import Storage - - for flow in self._fs.flows.values(): - key = f'flow|{flow.label_full}' - if key in values: - flow.previous_flow_rate = values[key] - - for comp in self._fs.components.values(): - if isinstance(comp, Storage): - key = f'storage|{comp.label_full}' - if key in values: - comp.initial_charge_state = values[key] - def _transfer_state( self, source_fs: FlowSystem, @@ -286,7 +251,6 @@ def _transfer_state( nr_of_previous_values: int, ) -> None: """Transfer final state from source segment to target segment.""" - from .components import Storage # Transfer flow rates (for uptime/downtime tracking) @@ -316,7 +280,6 @@ def _transfer_state( charge_var = f'{source_comp.label_full}|charge_state' if charge_var in source_fs.solution: # Get charge state at the end of the non-overlap portion - # Use horizon index (0-indexed, so horizon-1 is last non-overlap) charge_state = source_fs.solution[charge_var].isel(time=horizon - 1).values.item() target_comp.initial_charge_state = charge_state @@ -337,54 +300,82 @@ def _check_no_investments(self, segment_fs: FlowSystem) -> None: f'Use standard optimize() for problems with investments.' ) - def _combine_solutions(self, segment_flow_systems: list[FlowSystem], horizon: int) -> xr.Dataset: - """Combine segment solutions, trimming overlaps.""" + def _finalize_solution( + self, + segment_flow_systems: list[FlowSystem], + horizon: int, + ) -> None: + """Combine segment solutions and compute derived values directly (no re-solve).""" + # Combine all solution variables from segments + combined_solution = self._combine_solutions(segment_flow_systems, horizon) + + # Assign combined solution to the original FlowSystem + self._fs._solution = combined_solution + + def _combine_solutions( + self, + segment_flow_systems: list[FlowSystem], + horizon: int, + ) -> xr.Dataset: + """Combine segment solutions into a single Dataset, recomputing effect totals.""" if not segment_flow_systems: raise ValueError('No segments to combine.') - # Get all variable names from first segment + combined_vars: dict[str, xr.DataArray] = {} var_names = list(segment_flow_systems[0].solution.data_vars) - # Identify effect names for special handling - effect_labels = {e.label for e in self._fs.effects.values()} + # Identify effect-related variables for later recomputation + effect_totals = {} # effect_name -> will be recomputed + temporal_effects = {} # effect_name(temporal) -> will be recomputed - combined_vars = {} for var_name in var_names: - arrays = [] - for i, seg_fs in enumerate(segment_flow_systems): - da = seg_fs.solution[var_name] - - # Check if this variable has a time dimension - if 'time' in da.dims: + first_var = segment_flow_systems[0].solution[var_name] + + if 'time' in first_var.dims: + # Time-dependent variable: concatenate segments, trimming overlaps + arrays = [] + for i, seg_fs in enumerate(segment_flow_systems): + da = seg_fs.solution[var_name] + # Trim overlap for all segments except the last if i < len(segment_flow_systems) - 1: - # Not the last segment: trim to horizon (exclude overlap) da = da.isel(time=slice(None, horizon)) - # Last segment: keep all timesteps - arrays.append(da) - - # Concatenate along time if time dimension exists - if 'time' in segment_flow_systems[0].solution[var_name].dims: + arrays.append(da) combined_vars[var_name] = xr.concat(arrays, dim='time') else: - # For non-time scalars, just take the last value for now - # Effect totals will be recalculated below - combined_vars[var_name] = arrays[-1] - - # Recalculate effect totals from combined per-timestep data - for effect_label in effect_labels: - per_timestep_key = f'{effect_label}(temporal)|per_timestep' - temporal_key = f'{effect_label}(temporal)' - total_key = effect_label - - if per_timestep_key in combined_vars: - # Recalculate temporal sum from per-timestep values - combined_vars[temporal_key] = combined_vars[per_timestep_key].sum() - - # Recalculate total (temporal + periodic) - periodic_key = f'{effect_label}(periodic)' - if periodic_key in combined_vars: - combined_vars[total_key] = combined_vars[temporal_key] + combined_vars[periodic_key] + # Scalar variable: check if it's an effect total that needs recomputation + if var_name.endswith('(temporal)'): + # Will recompute from per_timestep values + effect_name = var_name.replace('(temporal)', '') + temporal_effects[effect_name] = var_name + elif var_name.endswith('(periodic)'): + # Sum periodic costs from all segments (no overlap trimming for periodic) + combined_vars[var_name] = sum(seg_fs.solution[var_name] for seg_fs in segment_flow_systems) + elif '(' not in var_name and any(f'{var_name}(temporal)|per_timestep' in v for v in var_names): + # This is an effect total - will recompute after we have temporal + periodic + effect_totals[var_name] = True else: - combined_vars[total_key] = combined_vars[temporal_key] + # Other scalar variables: take from last segment or sum as appropriate + # For most scalars, the last segment's value is representative + combined_vars[var_name] = segment_flow_systems[-1].solution[var_name] + + # Recompute temporal effect totals from per-timestep values + for effect_name, temporal_var_name in temporal_effects.items(): + per_timestep_name = f'{effect_name}(temporal)|per_timestep' + if per_timestep_name in combined_vars: + # Sum per-timestep values, ignoring NaN (final time point) + per_timestep = combined_vars[per_timestep_name] + temporal_total = per_timestep.sum(dim='time', skipna=True) + combined_vars[temporal_var_name] = temporal_total + + # Recompute total effect values (temporal + periodic) + for effect_name in effect_totals: + temporal_var = f'{effect_name}(temporal)' + periodic_var = f'{effect_name}(periodic)' + total = xr.DataArray(0.0) + if temporal_var in combined_vars: + total = total + combined_vars[temporal_var] + if periodic_var in combined_vars: + total = total + combined_vars[periodic_var] + combined_vars[effect_name] = total return xr.Dataset(combined_vars) From ac1384f33d5f2aaeb6c04da7e30aa06368f36edc Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:01:37 +0100 Subject: [PATCH 04/19] Add output supression --- flixopt/optimize_accessor.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py index e59c07df5..9549eb7d6 100644 --- a/flixopt/optimize_accessor.py +++ b/flixopt/optimize_accessor.py @@ -15,6 +15,7 @@ from tqdm import tqdm from .config import CONFIG +from .io import suppress_output if TYPE_CHECKING: from .flow_system import FlowSystem @@ -210,14 +211,19 @@ def rolling_horizon( nr_of_previous_values=nr_of_previous_values, ) - # Build and solve - segment_fs.build_model() - - # Check for investments (only on first segment) - if i == 0: - self._check_no_investments(segment_fs) + # Build and solve (suppress output when progress bar is shown) + if CONFIG.Solving.log_to_console: + with suppress_output(): + segment_fs.build_model() + if i == 0: + self._check_no_investments(segment_fs) + segment_fs.solve(solver) + else: + segment_fs.build_model() + if i == 0: + self._check_no_investments(segment_fs) + segment_fs.solve(solver) - segment_fs.solve(solver) segment_flow_systems.append(segment_fs) finally: From b46c6be9880d1093cde6dbb6dc8a2677fd093462 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:03:20 +0100 Subject: [PATCH 05/19] Improve output supression --- flixopt/optimize_accessor.py | 46 ++++++++++++++++++++++-------------- 1 file changed, 28 insertions(+), 18 deletions(-) diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py index 9549eb7d6..33ad4e1f5 100644 --- a/flixopt/optimize_accessor.py +++ b/flixopt/optimize_accessor.py @@ -199,26 +199,36 @@ def rolling_horizon( for i, (start_idx, end_idx) in progress_bar: progress_bar.set_description(f'Segment {i + 1}/{n_segments} (timesteps {start_idx}-{end_idx})') - # Create segment FlowSystem - segment_fs = self._fs.transform.isel(time=slice(start_idx, end_idx)) - - # Transfer state from previous segment - if i > 0 and nr_of_previous_values > 0: - self._transfer_state( - source_fs=segment_flow_systems[i - 1], - target_fs=segment_fs, - horizon=horizon, - nr_of_previous_values=nr_of_previous_values, - ) - - # Build and solve (suppress output when progress bar is shown) + # Suppress output when progress bar is shown (including logger and solver) if CONFIG.Solving.log_to_console: - with suppress_output(): - segment_fs.build_model() - if i == 0: - self._check_no_investments(segment_fs) - segment_fs.solve(solver) + # Temporarily raise logger level to suppress INFO messages + original_level = logger.level + logger.setLevel(logging.WARNING) + try: + with suppress_output(): + segment_fs = self._fs.transform.isel(time=slice(start_idx, end_idx)) + if i > 0 and nr_of_previous_values > 0: + self._transfer_state( + source_fs=segment_flow_systems[i - 1], + target_fs=segment_fs, + horizon=horizon, + nr_of_previous_values=nr_of_previous_values, + ) + segment_fs.build_model() + if i == 0: + self._check_no_investments(segment_fs) + segment_fs.solve(solver) + finally: + logger.setLevel(original_level) else: + segment_fs = self._fs.transform.isel(time=slice(start_idx, end_idx)) + if i > 0 and nr_of_previous_values > 0: + self._transfer_state( + source_fs=segment_flow_systems[i - 1], + target_fs=segment_fs, + horizon=horizon, + nr_of_previous_values=nr_of_previous_values, + ) segment_fs.build_model() if i == 0: self._check_no_investments(segment_fs) From dfb44d3e9cdb0699c01c4327ce317b60c6d7097f Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:05:37 +0100 Subject: [PATCH 06/19] Improve log message --- flixopt/optimize_accessor.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py index 33ad4e1f5..ded6d2506 100644 --- a/flixopt/optimize_accessor.py +++ b/flixopt/optimize_accessor.py @@ -178,10 +178,9 @@ def rolling_horizon( # Calculate segment indices segment_indices = self._calculate_segment_indices(total_timesteps, horizon, overlap) n_segments = len(segment_indices) - - logger.info(f'{"":#^80}') - logger.info(f'{" Rolling Horizon Optimization ":#^80}') - logger.info(f'Segments: {n_segments}, Horizon: {horizon}, Overlap: {overlap}') + logger.info( + f'Starting Rolling Horizon Optimization - Segments: {n_segments}, Horizon: {horizon}, Overlap: {overlap}' + ) # Create and solve segments segment_flow_systems: list[FlowSystem] = [] From 29ed8babf74cabaa09272ad6d96b1d2ca0f85d5c Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:10:16 +0100 Subject: [PATCH 07/19] Improve solution concatenation --- flixopt/optimize_accessor.py | 153 ++++++++++++++++++----------------- 1 file changed, 79 insertions(+), 74 deletions(-) diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py index ded6d2506..e126caf85 100644 --- a/flixopt/optimize_accessor.py +++ b/flixopt/optimize_accessor.py @@ -265,38 +265,30 @@ def _transfer_state( horizon: int, nr_of_previous_values: int, ) -> None: - """Transfer final state from source segment to target segment.""" + """Transfer final state from source segment to target segment. + + Transfers: + - Flow previous_flow_rate: Last nr_of_previous_values from non-overlap portion + - Storage initial_charge_state: Charge state at end of non-overlap portion + """ from .components import Storage + solution = source_fs.solution + time_slice = slice(horizon - nr_of_previous_values, horizon) + # Transfer flow rates (for uptime/downtime tracking) - for source_flow in source_fs.flows.values(): - target_flow = target_fs.flows.get(source_flow.label_full) - if target_flow is None: - continue - - # Get last nr_of_previous_values from source solution - flow_rate_var = f'{source_flow.label_full}|flow_rate' - if flow_rate_var in source_fs.solution: - # Select from the non-overlap portion (first 'horizon' timesteps) - values = ( - source_fs.solution[flow_rate_var].isel(time=slice(horizon - nr_of_previous_values, horizon)).values - ) - target_flow.previous_flow_rate = values if len(values) > 1 else values.item() + for label, target_flow in target_fs.flows.items(): + var_name = f'{label}|flow_rate' + if var_name in solution: + values = solution[var_name].isel(time=time_slice).values + target_flow.previous_flow_rate = values.item() if values.size == 1 else values # Transfer storage charge states - for source_comp in source_fs.components.values(): - if not isinstance(source_comp, Storage): - continue - - target_comp = target_fs.components.get(source_comp.label_full) - if target_comp is None or not isinstance(target_comp, Storage): - continue - - charge_var = f'{source_comp.label_full}|charge_state' - if charge_var in source_fs.solution: - # Get charge state at the end of the non-overlap portion - charge_state = source_fs.solution[charge_var].isel(time=horizon - 1).values.item() - target_comp.initial_charge_state = charge_state + for label, target_comp in target_fs.components.items(): + if isinstance(target_comp, Storage): + var_name = f'{label}|charge_state' + if var_name in solution: + target_comp.initial_charge_state = solution[var_name].isel(time=horizon - 1).item() def _check_no_investments(self, segment_fs: FlowSystem) -> None: """Check that no InvestParameters are used (not supported in rolling horizon).""" @@ -332,65 +324,78 @@ def _combine_solutions( segment_flow_systems: list[FlowSystem], horizon: int, ) -> xr.Dataset: - """Combine segment solutions into a single Dataset, recomputing effect totals.""" + """Combine segment solutions into a single Dataset, recomputing effect totals. + + Effect variables are identified from the FlowSystem's effect labels: + - {effect}: Total = temporal + periodic (recomputed) + - {effect}(temporal): Sum of per_timestep values (recomputed) + - {effect}(periodic): Sum across segments (no overlap) + - {effect}(temporal)|per_timestep: Time-dependent (concatenated with overlap trimming) + """ if not segment_flow_systems: raise ValueError('No segments to combine.') - combined_vars: dict[str, xr.DataArray] = {} - var_names = list(segment_flow_systems[0].solution.data_vars) + # Get effect labels from the original FlowSystem (includes Penalty) + effect_labels = set(self._fs.effects.keys()) - # Identify effect-related variables for later recomputation - effect_totals = {} # effect_name -> will be recomputed - temporal_effects = {} # effect_name(temporal) -> will be recomputed + combined_vars: dict[str, xr.DataArray] = {} + first_solution = segment_flow_systems[0].solution - for var_name in var_names: - first_var = segment_flow_systems[0].solution[var_name] + for var_name in first_solution.data_vars: + first_var = first_solution[var_name] if 'time' in first_var.dims: - # Time-dependent variable: concatenate segments, trimming overlaps - arrays = [] - for i, seg_fs in enumerate(segment_flow_systems): - da = seg_fs.solution[var_name] - # Trim overlap for all segments except the last - if i < len(segment_flow_systems) - 1: - da = da.isel(time=slice(None, horizon)) - arrays.append(da) - combined_vars[var_name] = xr.concat(arrays, dim='time') + # Time-dependent: concatenate with overlap trimming + combined_vars[var_name] = self._concat_time_variable(segment_flow_systems, var_name, horizon) + + elif var_name in effect_labels: + # Effect total: skip for now, will recompute from temporal + periodic + pass + + elif var_name.endswith('(temporal)') and var_name[:-10] in effect_labels: + # Temporal effect total: skip for now, will recompute from per_timestep + pass + + elif var_name.endswith('(periodic)') and var_name[:-10] in effect_labels: + # Periodic effect: sum across all segments (no overlap issue) + combined_vars[var_name] = sum(seg.solution[var_name] for seg in segment_flow_systems) + else: - # Scalar variable: check if it's an effect total that needs recomputation - if var_name.endswith('(temporal)'): - # Will recompute from per_timestep values - effect_name = var_name.replace('(temporal)', '') - temporal_effects[effect_name] = var_name - elif var_name.endswith('(periodic)'): - # Sum periodic costs from all segments (no overlap trimming for periodic) - combined_vars[var_name] = sum(seg_fs.solution[var_name] for seg_fs in segment_flow_systems) - elif '(' not in var_name and any(f'{var_name}(temporal)|per_timestep' in v for v in var_names): - # This is an effect total - will recompute after we have temporal + periodic - effect_totals[var_name] = True - else: - # Other scalar variables: take from last segment or sum as appropriate - # For most scalars, the last segment's value is representative - combined_vars[var_name] = segment_flow_systems[-1].solution[var_name] - - # Recompute temporal effect totals from per-timestep values - for effect_name, temporal_var_name in temporal_effects.items(): - per_timestep_name = f'{effect_name}(temporal)|per_timestep' - if per_timestep_name in combined_vars: - # Sum per-timestep values, ignoring NaN (final time point) - per_timestep = combined_vars[per_timestep_name] - temporal_total = per_timestep.sum(dim='time', skipna=True) - combined_vars[temporal_var_name] = temporal_total - - # Recompute total effect values (temporal + periodic) - for effect_name in effect_totals: - temporal_var = f'{effect_name}(temporal)' - periodic_var = f'{effect_name}(periodic)' + # Other scalar: use last segment's value + combined_vars[var_name] = segment_flow_systems[-1].solution[var_name] + + # Recompute effect totals from combined per-timestep values + for effect_label in effect_labels: + per_timestep_var = f'{effect_label}(temporal)|per_timestep' + temporal_var = f'{effect_label}(temporal)' + periodic_var = f'{effect_label}(periodic)' + + # Temporal total = sum of per_timestep + if per_timestep_var in combined_vars: + combined_vars[temporal_var] = combined_vars[per_timestep_var].sum(dim='time', skipna=True) + + # Total = temporal + periodic total = xr.DataArray(0.0) if temporal_var in combined_vars: total = total + combined_vars[temporal_var] if periodic_var in combined_vars: total = total + combined_vars[periodic_var] - combined_vars[effect_name] = total + combined_vars[effect_label] = total return xr.Dataset(combined_vars) + + def _concat_time_variable( + self, + segment_flow_systems: list[FlowSystem], + var_name: str, + horizon: int, + ) -> xr.DataArray: + """Concatenate a time-dependent variable from segments, trimming overlaps.""" + arrays = [] + for i, seg_fs in enumerate(segment_flow_systems): + da = seg_fs.solution[var_name] + # Trim overlap for all segments except the last + if i < len(segment_flow_systems) - 1: + da = da.isel(time=slice(None, horizon)) + arrays.append(da) + return xr.concat(arrays, dim='time') From e5c16e4fe2bb07984886df197e3fa8436ae51e11 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:12:50 +0100 Subject: [PATCH 08/19] Improve solution concatenation --- flixopt/optimize_accessor.py | 86 ++++++++++++------------------------ 1 file changed, 28 insertions(+), 58 deletions(-) diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py index e126caf85..7e9a4bec0 100644 --- a/flixopt/optimize_accessor.py +++ b/flixopt/optimize_accessor.py @@ -324,78 +324,48 @@ def _combine_solutions( segment_flow_systems: list[FlowSystem], horizon: int, ) -> xr.Dataset: - """Combine segment solutions into a single Dataset, recomputing effect totals. + """Combine segment solutions into a single Dataset. - Effect variables are identified from the FlowSystem's effect labels: - - {effect}: Total = temporal + periodic (recomputed) - - {effect}(temporal): Sum of per_timestep values (recomputed) - - {effect}(periodic): Sum across segments (no overlap) - - {effect}(temporal)|per_timestep: Time-dependent (concatenated with overlap trimming) + - Time-dependent variables: concatenated with overlap trimming + - Effect totals: recomputed from per-timestep values + - Other scalars: set to NaN (unknown how to combine) """ if not segment_flow_systems: raise ValueError('No segments to combine.') - # Get effect labels from the original FlowSystem (includes Penalty) effect_labels = set(self._fs.effects.keys()) - combined_vars: dict[str, xr.DataArray] = {} first_solution = segment_flow_systems[0].solution - for var_name in first_solution.data_vars: - first_var = first_solution[var_name] - + # Step 1: Concatenate all time-dependent variables + for var_name, first_var in first_solution.data_vars.items(): if 'time' in first_var.dims: - # Time-dependent: concatenate with overlap trimming - combined_vars[var_name] = self._concat_time_variable(segment_flow_systems, var_name, horizon) - - elif var_name in effect_labels: - # Effect total: skip for now, will recompute from temporal + periodic - pass - - elif var_name.endswith('(temporal)') and var_name[:-10] in effect_labels: - # Temporal effect total: skip for now, will recompute from per_timestep - pass - - elif var_name.endswith('(periodic)') and var_name[:-10] in effect_labels: - # Periodic effect: sum across all segments (no overlap issue) - combined_vars[var_name] = sum(seg.solution[var_name] for seg in segment_flow_systems) - + arrays = [ + seg.solution[var_name].isel( + time=slice(None, horizon if i < len(segment_flow_systems) - 1 else None) + ) + for i, seg in enumerate(segment_flow_systems) + ] + combined_vars[var_name] = xr.concat(arrays, dim='time') else: - # Other scalar: use last segment's value - combined_vars[var_name] = segment_flow_systems[-1].solution[var_name] + # Scalar: NaN placeholder (will recompute effects below) + combined_vars[var_name] = xr.DataArray(float('nan')) - # Recompute effect totals from combined per-timestep values - for effect_label in effect_labels: - per_timestep_var = f'{effect_label}(temporal)|per_timestep' - temporal_var = f'{effect_label}(temporal)' - periodic_var = f'{effect_label}(periodic)' + # Step 2: Recompute effect totals + for effect in effect_labels: + per_ts = f'{effect}(temporal)|per_timestep' + temporal = f'{effect}(temporal)' + periodic = f'{effect}(periodic)' - # Temporal total = sum of per_timestep - if per_timestep_var in combined_vars: - combined_vars[temporal_var] = combined_vars[per_timestep_var].sum(dim='time', skipna=True) + # Temporal = sum of per_timestep + if per_ts in combined_vars: + combined_vars[temporal] = combined_vars[per_ts].sum(dim='time', skipna=True) + + # Periodic = sum across segments (no overlap issue for periodic costs) + if periodic in first_solution: + combined_vars[periodic] = sum(seg.solution[periodic] for seg in segment_flow_systems) # Total = temporal + periodic - total = xr.DataArray(0.0) - if temporal_var in combined_vars: - total = total + combined_vars[temporal_var] - if periodic_var in combined_vars: - total = total + combined_vars[periodic_var] - combined_vars[effect_label] = total + combined_vars[effect] = combined_vars.get(temporal, 0.0) + combined_vars.get(periodic, 0.0) return xr.Dataset(combined_vars) - - def _concat_time_variable( - self, - segment_flow_systems: list[FlowSystem], - var_name: str, - horizon: int, - ) -> xr.DataArray: - """Concatenate a time-dependent variable from segments, trimming overlaps.""" - arrays = [] - for i, seg_fs in enumerate(segment_flow_systems): - da = seg_fs.solution[var_name] - # Trim overlap for all segments except the last - if i < len(segment_flow_systems) - 1: - da = da.isel(time=slice(None, horizon)) - arrays.append(da) - return xr.concat(arrays, dim='time') From c4ed2f3474a651f4b3746e2d1698726bdd109f82 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:21:24 +0100 Subject: [PATCH 09/19] Improve solution concatenation --- flixopt/optimize_accessor.py | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py index 7e9a4bec0..f88cdf982 100644 --- a/flixopt/optimize_accessor.py +++ b/flixopt/optimize_accessor.py @@ -288,7 +288,7 @@ def _transfer_state( if isinstance(target_comp, Storage): var_name = f'{label}|charge_state' if var_name in solution: - target_comp.initial_charge_state = solution[var_name].isel(time=horizon - 1).item() + target_comp.initial_charge_state = solution[var_name].isel(time=horizon).item() def _check_no_investments(self, segment_fs: FlowSystem) -> None: """Check that no InvestParameters are used (not supported in rolling horizon).""" @@ -327,8 +327,8 @@ def _combine_solutions( """Combine segment solutions into a single Dataset. - Time-dependent variables: concatenated with overlap trimming - - Effect totals: recomputed from per-timestep values - - Other scalars: set to NaN (unknown how to combine) + - Effect temporal/total: recomputed from per-timestep values + - Other scalars (including periodic): NaN (not meaningful for rolling horizon) """ if not segment_flow_systems: raise ValueError('No segments to combine.') @@ -337,7 +337,7 @@ def _combine_solutions( combined_vars: dict[str, xr.DataArray] = {} first_solution = segment_flow_systems[0].solution - # Step 1: Concatenate all time-dependent variables + # Step 1: Time-dependent → concatenate; Scalars → NaN for var_name, first_var in first_solution.data_vars.items(): if 'time' in first_var.dims: arrays = [ @@ -348,24 +348,14 @@ def _combine_solutions( ] combined_vars[var_name] = xr.concat(arrays, dim='time') else: - # Scalar: NaN placeholder (will recompute effects below) combined_vars[var_name] = xr.DataArray(float('nan')) - # Step 2: Recompute effect totals + # Step 2: Recompute effect totals from per-timestep values for effect in effect_labels: per_ts = f'{effect}(temporal)|per_timestep' - temporal = f'{effect}(temporal)' - periodic = f'{effect}(periodic)' - - # Temporal = sum of per_timestep if per_ts in combined_vars: - combined_vars[temporal] = combined_vars[per_ts].sum(dim='time', skipna=True) - - # Periodic = sum across segments (no overlap issue for periodic costs) - if periodic in first_solution: - combined_vars[periodic] = sum(seg.solution[periodic] for seg in segment_flow_systems) - - # Total = temporal + periodic - combined_vars[effect] = combined_vars.get(temporal, 0.0) + combined_vars.get(periodic, 0.0) + temporal_sum = combined_vars[per_ts].sum(dim='time', skipna=True) + combined_vars[f'{effect}(temporal)'] = temporal_sum + combined_vars[effect] = temporal_sum # Total = temporal (periodic is NaN/unsupported) return xr.Dataset(combined_vars) From 47799494aac7d3c83bc00019ca2f2d8235ac9266 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:28:09 +0100 Subject: [PATCH 10/19] Adding notebook for rolling --- docs/notebooks/08b-rolling-horizon.ipynb | 445 +++++++++++++++++++++++ mkdocs.yml | 1 + 2 files changed, 446 insertions(+) create mode 100644 docs/notebooks/08b-rolling-horizon.ipynb diff --git a/docs/notebooks/08b-rolling-horizon.ipynb b/docs/notebooks/08b-rolling-horizon.ipynb new file mode 100644 index 000000000..90efe99bb --- /dev/null +++ b/docs/notebooks/08b-rolling-horizon.ipynb @@ -0,0 +1,445 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Rolling Horizon\n", + "\n", + "Solve large operational problems by decomposing the time horizon into sequential segments.\n", + "\n", + "This notebook introduces:\n", + "\n", + "- **Rolling horizon optimization**: Divide time into overlapping segments\n", + "- **State transfer**: Pass storage states and flow history between segments\n", + "- **When to use**: Memory limits, operational planning with limited foresight" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Setup" + ] + }, + { + "cell_type": "code", + "id": "2", + "metadata": {}, + "source": [ + "import timeit\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "import flixopt as fx\n", + "\n", + "fx.CONFIG.notebook()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## Create a Test System\n", + "\n", + "We'll use a simple heating system with storage to demonstrate the rolling horizon approach:" + ] + }, + { + "cell_type": "code", + "id": "4", + "metadata": {}, + "source": [ + "# One week at hourly resolution\n", + "timesteps = pd.date_range('2024-01-01', periods=168, freq='h') # 7 days\n", + "hours = np.arange(len(timesteps))\n", + "hour_of_day = hours % 24\n", + "\n", + "np.random.seed(42)\n", + "\n", + "# Heat demand: daily pattern\n", + "daily_pattern = np.select(\n", + " [\n", + " (hour_of_day >= 6) & (hour_of_day < 9),\n", + " (hour_of_day >= 9) & (hour_of_day < 17),\n", + " (hour_of_day >= 17) & (hour_of_day < 22),\n", + " ],\n", + " [80, 50, 70],\n", + " default=30,\n", + ").astype(float)\n", + "\n", + "heat_demand = daily_pattern + np.random.normal(0, 5, len(timesteps))\n", + "heat_demand = np.clip(heat_demand, 20, 100)\n", + "\n", + "print(f'Timesteps: {len(timesteps)} hours ({len(timesteps) / 24:.0f} days)')\n", + "print(f'Heat demand: {heat_demand.min():.0f} - {heat_demand.max():.0f} kW')" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "id": "5", + "metadata": {}, + "source": [ + "def build_system(timesteps, heat_demand):\n", + " \"\"\"Build a simple heating system with storage.\"\"\"\n", + " fs = fx.FlowSystem(timesteps)\n", + " fs.add_elements(\n", + " # Buses\n", + " fx.Bus('Gas'),\n", + " fx.Bus('Heat'),\n", + " # Effects\n", + " fx.Effect('costs', '€', 'Total Costs', is_standard=True, is_objective=True),\n", + " # Gas Supply\n", + " fx.Source(\n", + " 'GasGrid',\n", + " outputs=[fx.Flow('Gas', bus='Gas', size=500, effects_per_flow_hour=0.05)],\n", + " ),\n", + " # Gas Boiler\n", + " fx.linear_converters.Boiler(\n", + " 'Boiler',\n", + " thermal_efficiency=0.92,\n", + " thermal_flow=fx.Flow('Q_th', bus='Heat', size=150),\n", + " fuel_flow=fx.Flow('Q_fuel', bus='Gas'),\n", + " ),\n", + " # Thermal Storage\n", + " fx.Storage(\n", + " 'Storage',\n", + " capacity_in_flow_hours=100,\n", + " initial_charge_state=0,\n", + " eta_charge=0.95,\n", + " eta_discharge=0.95,\n", + " relative_loss_per_hour=0.01,\n", + " charging=fx.Flow('Charge', bus='Heat', size=50),\n", + " discharging=fx.Flow('Discharge', bus='Heat', size=50),\n", + " ),\n", + " # Heat Demand\n", + " fx.Sink(\n", + " 'HeatDemand',\n", + " inputs=[fx.Flow('Q_th', bus='Heat', size=1, fixed_relative_profile=heat_demand)],\n", + " ),\n", + " )\n", + " return fs\n", + "\n", + "\n", + "flow_system = build_system(timesteps, heat_demand)\n", + "print(f'System: {len(timesteps)} timesteps')" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## Full Optimization (Baseline)\n", + "\n", + "First, solve the full problem as a baseline:" + ] + }, + { + "cell_type": "code", + "id": "7", + "metadata": {}, + "source": [ + "solver = fx.solvers.HighsSolver()\n", + "\n", + "start = timeit.default_timer()\n", + "fs_full = flow_system.copy()\n", + "fs_full.optimize(solver)\n", + "time_full = timeit.default_timer() - start\n", + "\n", + "print(f'Full optimization: {time_full:.2f} seconds')\n", + "print(f'Cost: {fs_full.solution[\"costs\"].item():.2f} €')" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Rolling Horizon Optimization\n", + "\n", + "The `optimize.rolling_horizon()` method divides the time horizon into segments that are solved sequentially:\n", + "\n", + "```\n", + "Full horizon: |-------- 168 hours ---------------------------------------...|\n", + " \n", + "Segment 1: |==== 24h ====|-- overlap --|\n", + "Segment 2: |==== 24h ====|-- overlap --|\n", + "Segment 3: |==== 24h ====|-- overlap --|\n", + "... \n", + "```\n", + "\n", + "Key parameters:\n", + "- **horizon**: Timesteps per segment (excluding overlap)\n", + "- **overlap**: Additional lookahead timesteps (improves storage optimization)\n", + "- **nr_of_previous_values**: Flow history transferred between segments" + ] + }, + { + "cell_type": "code", + "id": "9", + "metadata": {}, + "source": [ + "start = timeit.default_timer()\n", + "fs_rolling = flow_system.copy()\n", + "segments = fs_rolling.optimize.rolling_horizon(\n", + " solver,\n", + " horizon=24, # Daily segments\n", + " overlap=6, # 6-hour lookahead\n", + ")\n", + "time_rolling = timeit.default_timer() - start\n", + "\n", + "print(f'Rolling horizon: {time_rolling:.2f} seconds ({len(segments)} segments)')\n", + "print(f'Cost: {fs_rolling.solution[\"costs\"].item():.2f} €')" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## Compare Results" + ] + }, + { + "cell_type": "code", + "id": "11", + "metadata": {}, + "source": [ + "cost_full = fs_full.solution['costs'].item()\n", + "cost_rolling = fs_rolling.solution['costs'].item()\n", + "cost_gap = (cost_rolling - cost_full) / cost_full * 100\n", + "\n", + "results = pd.DataFrame(\n", + " {\n", + " 'Method': ['Full optimization', 'Rolling horizon'],\n", + " 'Time [s]': [time_full, time_rolling],\n", + " 'Cost [€]': [cost_full, cost_rolling],\n", + " 'Cost Gap [%]': [0.0, cost_gap],\n", + " }\n", + ").set_index('Method')\n", + "\n", + "results.round(2)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## Visualize: Heat Balance Comparison" + ] + }, + { + "cell_type": "code", + "id": "13", + "metadata": {}, + "source": [ + "fs_full.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Full)')" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "id": "14", + "metadata": {}, + "source": [ + "fs_rolling.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Rolling)')" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "## Storage State Continuity\n", + "\n", + "Rolling horizon transfers storage charge states between segments to ensure continuity:" + ] + }, + { + "cell_type": "code", + "id": "16", + "metadata": {}, + "source": [ + "import plotly.graph_objects as go\n", + "from plotly.subplots import make_subplots\n", + "\n", + "fig = make_subplots(\n", + " rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.1, subplot_titles=['Full Optimization', 'Rolling Horizon']\n", + ")\n", + "\n", + "# Full optimization\n", + "charge_full = fs_full.solution['Storage|charge_state'].values[:-1] # Drop final NaN\n", + "fig.add_trace(go.Scatter(x=timesteps, y=charge_full, name='Full', line=dict(color='blue')), row=1, col=1)\n", + "\n", + "# Rolling horizon\n", + "charge_rolling = fs_rolling.solution['Storage|charge_state'].values[:-1]\n", + "fig.add_trace(go.Scatter(x=timesteps, y=charge_rolling, name='Rolling', line=dict(color='orange')), row=2, col=1)\n", + "\n", + "fig.update_yaxes(title_text='Charge State [kWh]', row=1, col=1)\n", + "fig.update_yaxes(title_text='Charge State [kWh]', row=2, col=1)\n", + "fig.update_layout(height=400, showlegend=False)\n", + "fig.show()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Inspect Individual Segments\n", + "\n", + "The method returns the individual segment FlowSystems, which can be inspected:" + ] + }, + { + "cell_type": "code", + "id": "18", + "metadata": {}, + "source": [ + "print(f'Number of segments: {len(segments)}')\n", + "print()\n", + "for i, seg in enumerate(segments):\n", + " start_time = seg.timesteps[0]\n", + " end_time = seg.timesteps[-1]\n", + " cost = seg.solution['costs'].item()\n", + " print(f'Segment {i + 1}: {start_time} → {end_time} | Cost: {cost:.2f} €')" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "## Effect of Overlap\n", + "\n", + "The overlap parameter provides lookahead for storage optimization. Let's compare different overlap values:" + ] + }, + { + "cell_type": "code", + "id": "20", + "metadata": {}, + "source": [ + "overlaps = [0, 3, 6, 12, 24]\n", + "overlap_results = []\n", + "\n", + "for overlap in overlaps:\n", + " fs = flow_system.copy()\n", + " start = timeit.default_timer()\n", + " fs.optimize.rolling_horizon(solver, horizon=24, overlap=overlap)\n", + " elapsed = timeit.default_timer() - start\n", + " cost = fs.solution['costs'].item()\n", + " gap = (cost - cost_full) / cost_full * 100\n", + " overlap_results.append({'Overlap [h]': overlap, 'Time [s]': elapsed, 'Cost [€]': cost, 'Gap [%]': gap})\n", + "\n", + "pd.DataFrame(overlap_results).round(2)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "## When to Use Rolling Horizon\n", + "\n", + "| Use Case | Recommendation |\n", + "|----------|----------------|\n", + "| **Memory limits** | Large problems that exceed available memory |\n", + "| **Operational planning** | When limited foresight is realistic |\n", + "| **Quick approximate solutions** | Faster than full optimization |\n", + "| **Investment decisions** | Use full optimization instead |\n", + "\n", + "### Limitations\n", + "\n", + "- **No investments**: `InvestParameters` are not supported (raises error)\n", + "- **Suboptimal storage**: Limited foresight may miss long-term storage opportunities\n", + "- **Global constraints**: `flow_hours_max` etc. cannot be enforced globally" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "## API Reference\n", + "\n", + "```python\n", + "segments = flow_system.optimize.rolling_horizon(\n", + " solver, # Solver instance\n", + " horizon=100, # Timesteps per segment\n", + " overlap=0, # Additional lookahead timesteps\n", + " nr_of_previous_values=1, # Flow history for uptime/downtime tracking\n", + ")\n", + "\n", + "# Combined solution on original FlowSystem\n", + "flow_system.solution['costs'].item()\n", + "\n", + "# Individual segment solutions\n", + "for seg in segments:\n", + " print(seg.solution['costs'].item())\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "You learned how to:\n", + "\n", + "- Use **`optimize.rolling_horizon()`** to decompose large problems\n", + "- Choose **horizon** and **overlap** parameters\n", + "- Understand the **trade-offs** vs. full optimization\n", + "\n", + "### Key Takeaways\n", + "\n", + "1. **Rolling horizon** is useful for memory-limited or operational planning problems\n", + "2. **Overlap** improves solution quality at the cost of computation time\n", + "3. **Storage states** are automatically transferred between segments\n", + "4. Use **full optimization** for investment decisions" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/mkdocs.yml b/mkdocs.yml index 186e109fd..bdc21fabc 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -70,6 +70,7 @@ nav: - Scaling: - Scenarios: notebooks/07-scenarios-and-periods.ipynb - Large-Scale: notebooks/08-large-scale-optimization.ipynb + - Rolling Horizon: notebooks/08b-rolling-horizon.ipynb - Results: - Plotting: notebooks/09-plotting-and-data-access.ipynb From 47bb27e250f443140d00c5ded30fd7790a223e78 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:28:30 +0100 Subject: [PATCH 11/19] Update migration guide --- docs/user-guide/migration-guide-v5.md | 28 ++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/docs/user-guide/migration-guide-v5.md b/docs/user-guide/migration-guide-v5.md index 044f26897..0c43e18f0 100644 --- a/docs/user-guide/migration-guide-v5.md +++ b/docs/user-guide/migration-guide-v5.md @@ -325,17 +325,27 @@ Clustered optimization uses the new transform accessor: # Results in clustered_fs.solution ``` -### Segmented Optimization (Not Yet Migrated) +### Segmented / Rolling Horizon Optimization -Segmented optimization still uses the class-based API. A new `optimize.rolling_horizon()` method is planned for a future release. +=== "v4.x (Old)" + ```python + calc = fx.SegmentedOptimization('model', flow_system, + timesteps_per_segment=96) + calc.do_modeling_and_solve(solver) + results = calc.results # Returns SegmentedResults + ``` -```python -# Still use the class-based API (unchanged from v4.x) -calc = fx.SegmentedOptimization('model', flow_system, - timesteps_per_segment=96) -calc.do_modeling_and_solve(solver) -results = calc.results # Returns SegmentedResults -``` +=== "v5.0.0 (New)" + ```python + # Use optimize.rolling_horizon() method + segments = flow_system.optimize.rolling_horizon( + solver, + horizon=96, # Timesteps per segment + overlap=12, # Lookahead for storage optimization + ) + # Combined solution on original FlowSystem + flow_system.solution['costs'].item() + ``` --- From ae143bff95674dde84af396ff3c8402f7293edb1 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:36:12 +0100 Subject: [PATCH 12/19] Improve notebook --- docs/notebooks/08b-rolling-horizon.ipynb | 321 +++++++++++++---------- 1 file changed, 184 insertions(+), 137 deletions(-) diff --git a/docs/notebooks/08b-rolling-horizon.ipynb b/docs/notebooks/08b-rolling-horizon.ipynb index 90efe99bb..fae589ff6 100644 --- a/docs/notebooks/08b-rolling-horizon.ipynb +++ b/docs/notebooks/08b-rolling-horizon.ipynb @@ -4,17 +4,7 @@ "cell_type": "markdown", "id": "0", "metadata": {}, - "source": [ - "# Rolling Horizon\n", - "\n", - "Solve large operational problems by decomposing the time horizon into sequential segments.\n", - "\n", - "This notebook introduces:\n", - "\n", - "- **Rolling horizon optimization**: Divide time into overlapping segments\n", - "- **State transfer**: Pass storage states and flow history between segments\n", - "- **When to use**: Memory limits, operational planning with limited foresight" - ] + "source": "# Rolling Horizon\n\nSolve large operational problems by decomposing the time horizon into sequential segments.\n\nThis notebook introduces:\n\n- **Rolling horizon optimization**: Divide time into overlapping segments\n- **State transfer**: Pass storage states and flow history between segments\n- **When to use**: Memory limits, operational planning with limited foresight\n\nWe use a realistic district heating system with CHP, boiler, and storage to demonstrate the approach." }, { "cell_type": "markdown", @@ -26,114 +16,188 @@ }, { "cell_type": "code", + "execution_count": 3, "id": "2", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:35:42.432171Z", + "start_time": "2025-12-13T18:35:42.279884Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "flixopt.config.CONFIG" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "import timeit\n", "\n", - "import numpy as np\n", "import pandas as pd\n", + "import plotly.graph_objects as go\n", + "from plotly.subplots import make_subplots\n", "\n", "import flixopt as fx\n", "\n", "fx.CONFIG.notebook()" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, - "source": [ - "## Create a Test System\n", - "\n", - "We'll use a simple heating system with storage to demonstrate the rolling horizon approach:" - ] + "source": "## Load Time Series Data\n\nWe use real-world district heating data at 15-minute resolution (two weeks):" }, { "cell_type": "code", + "execution_count": 4, "id": "4", - "metadata": {}, - "source": [ - "# One week at hourly resolution\n", - "timesteps = pd.date_range('2024-01-01', periods=168, freq='h') # 7 days\n", - "hours = np.arange(len(timesteps))\n", - "hour_of_day = hours % 24\n", - "\n", - "np.random.seed(42)\n", - "\n", - "# Heat demand: daily pattern\n", - "daily_pattern = np.select(\n", - " [\n", - " (hour_of_day >= 6) & (hour_of_day < 9),\n", - " (hour_of_day >= 9) & (hour_of_day < 17),\n", - " (hour_of_day >= 17) & (hour_of_day < 22),\n", - " ],\n", - " [80, 50, 70],\n", - " default=30,\n", - ").astype(float)\n", - "\n", - "heat_demand = daily_pattern + np.random.normal(0, 5, len(timesteps))\n", - "heat_demand = np.clip(heat_demand, 20, 100)\n", - "\n", - "print(f'Timesteps: {len(timesteps)} hours ({len(timesteps) / 24:.0f} days)')\n", - "print(f'Heat demand: {heat_demand.min():.0f} - {heat_demand.max():.0f} kW')" + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:35:43.695556Z", + "start_time": "2025-12-13T18:35:42.878212Z" + } + }, + "outputs": [ + { + "ename": "KeyError", + "evalue": "'Value based partial slicing on non-monotonic DatetimeIndexes with non-existing keys is not allowed.'", + "output_type": "error", + "traceback": [ + "\u001b[31m---------------------------------------------------------------------------\u001b[39m", + "\u001b[31mKeyError\u001b[39m Traceback (most recent call last)", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[4]\u001b[39m\u001b[32m, line 9\u001b[39m\n\u001b[32m 6\u001b[39m notebook_dir = pathlib.Path(\u001b[34m__file__\u001b[39m).parent / \u001b[33m'\u001b[39m\u001b[33mdata\u001b[39m\u001b[33m'\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[33m'\u001b[39m\u001b[33m__file__\u001b[39m\u001b[33m'\u001b[39m \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mdir\u001b[39m() \u001b[38;5;28;01melse\u001b[39;00m pathlib.Path(\u001b[33m'\u001b[39m\u001b[33mdocs/notebooks/data\u001b[39m\u001b[33m'\u001b[39m)\n\u001b[32m 8\u001b[39m data = pd.read_csv(notebook_dir / \u001b[33m'\u001b[39m\u001b[33mZeitreihen2020.csv\u001b[39m\u001b[33m'\u001b[39m, index_col=\u001b[32m0\u001b[39m, parse_dates=\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[32m----> \u001b[39m\u001b[32m9\u001b[39m data = \u001b[43mdata\u001b[49m\u001b[43m[\u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43m2020-01-01\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m:\u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43m2020-01-14 23:45:00\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m]\u001b[49m \u001b[38;5;66;03m# Two weeks\u001b[39;00m\n\u001b[32m 11\u001b[39m timesteps = data.index\n\u001b[32m 13\u001b[39m \u001b[38;5;66;03m# Extract profiles\u001b[39;00m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/PycharmProjects/flixopt_719231/.venv/lib/python3.11/site-packages/pandas/core/frame.py:4096\u001b[39m, in \u001b[36mDataFrame.__getitem__\u001b[39m\u001b[34m(self, key)\u001b[39m\n\u001b[32m 4094\u001b[39m \u001b[38;5;66;03m# Do we have a slicer (on rows)?\u001b[39;00m\n\u001b[32m 4095\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(key, \u001b[38;5;28mslice\u001b[39m):\n\u001b[32m-> \u001b[39m\u001b[32m4096\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_getitem_slice\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 4098\u001b[39m \u001b[38;5;66;03m# Do we have a (boolean) DataFrame?\u001b[39;00m\n\u001b[32m 4099\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(key, DataFrame):\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/PycharmProjects/flixopt_719231/.venv/lib/python3.11/site-packages/pandas/core/generic.py:4371\u001b[39m, in \u001b[36mNDFrame._getitem_slice\u001b[39m\u001b[34m(self, key)\u001b[39m\n\u001b[32m 4366\u001b[39m \u001b[38;5;250m\u001b[39m\u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 4367\u001b[39m \u001b[33;03m__getitem__ for the case where the key is a slice object.\u001b[39;00m\n\u001b[32m 4368\u001b[39m \u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 4369\u001b[39m \u001b[38;5;66;03m# _convert_slice_indexer to determine if this slice is positional\u001b[39;00m\n\u001b[32m 4370\u001b[39m \u001b[38;5;66;03m# or label based, and if the latter, convert to positional\u001b[39;00m\n\u001b[32m-> \u001b[39m\u001b[32m4371\u001b[39m slobj = \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mindex\u001b[49m\u001b[43m.\u001b[49m\u001b[43m_convert_slice_indexer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkind\u001b[49m\u001b[43m=\u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mgetitem\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[32m 4372\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(slobj, np.ndarray):\n\u001b[32m 4373\u001b[39m \u001b[38;5;66;03m# reachable with DatetimeIndex\u001b[39;00m\n\u001b[32m 4374\u001b[39m indexer = lib.maybe_indices_to_slice(\n\u001b[32m 4375\u001b[39m slobj.astype(np.intp, copy=\u001b[38;5;28;01mFalse\u001b[39;00m), \u001b[38;5;28mlen\u001b[39m(\u001b[38;5;28mself\u001b[39m)\n\u001b[32m 4376\u001b[39m )\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/PycharmProjects/flixopt_719231/.venv/lib/python3.11/site-packages/pandas/core/indexes/base.py:4288\u001b[39m, in \u001b[36mIndex._convert_slice_indexer\u001b[39m\u001b[34m(self, key, kind)\u001b[39m\n\u001b[32m 4286\u001b[39m indexer = key\n\u001b[32m 4287\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m-> \u001b[39m\u001b[32m4288\u001b[39m indexer = \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mslice_indexer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstart\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstop\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstep\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 4290\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m indexer\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/PycharmProjects/flixopt_719231/.venv/lib/python3.11/site-packages/pandas/core/indexes/datetimes.py:697\u001b[39m, in \u001b[36mDatetimeIndex.slice_indexer\u001b[39m\u001b[34m(self, start, end, step)\u001b[39m\n\u001b[32m 694\u001b[39m in_index &= (end_casted == \u001b[38;5;28mself\u001b[39m).any()\n\u001b[32m 696\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m in_index:\n\u001b[32m--> \u001b[39m\u001b[32m697\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m(\n\u001b[32m 698\u001b[39m \u001b[33m\"\u001b[39m\u001b[33mValue based partial slicing on non-monotonic DatetimeIndexes \u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 699\u001b[39m \u001b[33m\"\u001b[39m\u001b[33mwith non-existing keys is not allowed.\u001b[39m\u001b[33m\"\u001b[39m,\n\u001b[32m 700\u001b[39m )\n\u001b[32m 701\u001b[39m indexer = mask.nonzero()[\u001b[32m0\u001b[39m][::step]\n\u001b[32m 702\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(indexer) == \u001b[38;5;28mlen\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n", + "\u001b[31mKeyError\u001b[39m: 'Value based partial slicing on non-monotonic DatetimeIndexes with non-existing keys is not allowed.'" + ] + } ], - "outputs": [], - "execution_count": null + "source": [ + "import pathlib\n", + "\n", + "# Load time series data (15-min resolution)\n", + "notebook_dir = pathlib.Path('data')\n", + "if not notebook_dir.exists():\n", + " notebook_dir = (\n", + " pathlib.Path(__file__).parent / 'data' if '__file__' in dir() else pathlib.Path('docs/notebooks/data')\n", + " )\n", + "\n", + "data = pd.read_csv(notebook_dir / 'Zeitreihen2020.csv', index_col=0, parse_dates=True).sort_index()\n", + "data = data['2020-01-01':'2020-01-14 23:45:00'] # Two weeks\n", + "\n", + "timesteps = data.index\n", + "\n", + "# Extract profiles\n", + "electricity_demand = data['P_Netz/MW'].to_numpy()\n", + "heat_demand = data['Q_Netz/MW'].to_numpy()\n", + "electricity_price = data['Strompr.€/MWh'].to_numpy()\n", + "gas_price = data['Gaspr.€/MWh'].to_numpy()\n", + "\n", + "print(f'Timesteps: {len(timesteps)} ({len(timesteps) / 96:.0f} days at 15-min resolution)')\n", + "print(f'Heat demand: {heat_demand.min():.1f} - {heat_demand.max():.1f} MW')\n", + "print(f'Electricity price: {electricity_price.min():.1f} - {electricity_price.max():.1f} €/MWh')" + ] }, { "cell_type": "code", + "execution_count": null, "id": "5", "metadata": {}, + "outputs": [], "source": [ - "def build_system(timesteps, heat_demand):\n", - " \"\"\"Build a simple heating system with storage.\"\"\"\n", + "def build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price):\n", + " \"\"\"Build a district heating system with CHP, boiler, and storage.\"\"\"\n", " fs = fx.FlowSystem(timesteps)\n", + "\n", + " # Effects\n", + "\n", + " # Buses\n", " fs.add_elements(\n", - " # Buses\n", - " fx.Bus('Gas'),\n", + " fx.Bus('Electricity'),\n", " fx.Bus('Heat'),\n", - " # Effects\n", + " fx.Bus('Gas'),\n", + " fx.Bus('Coal'),\n", " fx.Effect('costs', '€', 'Total Costs', is_standard=True, is_objective=True),\n", - " # Gas Supply\n", - " fx.Source(\n", - " 'GasGrid',\n", - " outputs=[fx.Flow('Gas', bus='Gas', size=500, effects_per_flow_hour=0.05)],\n", + " fx.Effect('CO2', 'kg', 'CO2 Emissions'),\n", + " fx.linear_converters.CHP(\n", + " 'CHP',\n", + " thermal_efficiency=0.58,\n", + " electrical_efficiency=0.22,\n", + " status_parameters=fx.StatusParameters(effects_per_startup=24000),\n", + " electrical_flow=fx.Flow('P_el', bus='Electricity', size=200),\n", + " thermal_flow=fx.Flow('Q_th', bus='Heat', size=200),\n", + " fuel_flow=fx.Flow('Q_fu', bus='Coal', size=288, relative_minimum=87 / 288, previous_flow_rate=100),\n", " ),\n", - " # Gas Boiler\n", " fx.linear_converters.Boiler(\n", " 'Boiler',\n", - " thermal_efficiency=0.92,\n", - " thermal_flow=fx.Flow('Q_th', bus='Heat', size=150),\n", - " fuel_flow=fx.Flow('Q_fuel', bus='Gas'),\n", + " thermal_efficiency=0.85,\n", + " thermal_flow=fx.Flow('Q_th', bus='Heat'),\n", + " fuel_flow=fx.Flow(\n", + " 'Q_fu',\n", + " bus='Gas',\n", + " size=95,\n", + " relative_minimum=12 / 95,\n", + " previous_flow_rate=20,\n", + " status_parameters=fx.StatusParameters(effects_per_startup=1000),\n", + " ),\n", " ),\n", - " # Thermal Storage\n", " fx.Storage(\n", " 'Storage',\n", - " capacity_in_flow_hours=100,\n", - " initial_charge_state=0,\n", - " eta_charge=0.95,\n", - " eta_discharge=0.95,\n", - " relative_loss_per_hour=0.01,\n", - " charging=fx.Flow('Charge', bus='Heat', size=50),\n", - " discharging=fx.Flow('Discharge', bus='Heat', size=50),\n", + " capacity_in_flow_hours=684,\n", + " initial_charge_state=137,\n", + " minimal_final_charge_state=137,\n", + " maximal_final_charge_state=158,\n", + " eta_charge=1,\n", + " eta_discharge=1,\n", + " relative_loss_per_hour=0.001,\n", + " prevent_simultaneous_charge_and_discharge=True,\n", + " charging=fx.Flow('Charge', size=137, bus='Heat'),\n", + " discharging=fx.Flow('Discharge', size=158, bus='Heat'),\n", + " ),\n", + " fx.Source(\n", + " 'GasGrid',\n", + " outputs=[fx.Flow('Q_Gas', bus='Gas', size=1000, effects_per_flow_hour={'costs': gas_price, 'CO2': 0.3})],\n", + " ),\n", + " fx.Source(\n", + " 'CoalSupply',\n", + " outputs=[fx.Flow('Q_Coal', bus='Coal', size=1000, effects_per_flow_hour={'costs': 4.6, 'CO2': 0.3})],\n", + " ),\n", + " fx.Source(\n", + " 'GridBuy',\n", + " outputs=[\n", + " fx.Flow(\n", + " 'P_el',\n", + " bus='Electricity',\n", + " size=1000,\n", + " effects_per_flow_hour={'costs': electricity_price + 0.5, 'CO2': 0.3},\n", + " )\n", + " ],\n", + " ),\n", + " fx.Sink(\n", + " 'GridSell',\n", + " inputs=[fx.Flow('P_el', bus='Electricity', size=1000, effects_per_flow_hour=-(electricity_price - 0.5))],\n", " ),\n", - " # Heat Demand\n", + " fx.Sink('HeatDemand', inputs=[fx.Flow('Q_th', bus='Heat', size=1, fixed_relative_profile=heat_demand)]),\n", " fx.Sink(\n", - " 'HeatDemand',\n", - " inputs=[fx.Flow('Q_th', bus='Heat', size=1, fixed_relative_profile=heat_demand)],\n", + " 'ElecDemand', inputs=[fx.Flow('P_el', bus='Electricity', size=1, fixed_relative_profile=electricity_demand)]\n", " ),\n", " )\n", + "\n", " return fs\n", "\n", "\n", - "flow_system = build_system(timesteps, heat_demand)\n", + "flow_system = build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price)\n", "print(f'System: {len(timesteps)} timesteps')" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -147,8 +211,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "7", "metadata": {}, + "outputs": [], "source": [ "solver = fx.solvers.HighsSolver()\n", "\n", @@ -158,54 +224,34 @@ "time_full = timeit.default_timer() - start\n", "\n", "print(f'Full optimization: {time_full:.2f} seconds')\n", - "print(f'Cost: {fs_full.solution[\"costs\"].item():.2f} €')" - ], - "outputs": [], - "execution_count": null + "print(f'Cost: {fs_full.solution[\"costs\"].item():,.0f} €')" + ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, - "source": [ - "## Rolling Horizon Optimization\n", - "\n", - "The `optimize.rolling_horizon()` method divides the time horizon into segments that are solved sequentially:\n", - "\n", - "```\n", - "Full horizon: |-------- 168 hours ---------------------------------------...|\n", - " \n", - "Segment 1: |==== 24h ====|-- overlap --|\n", - "Segment 2: |==== 24h ====|-- overlap --|\n", - "Segment 3: |==== 24h ====|-- overlap --|\n", - "... \n", - "```\n", - "\n", - "Key parameters:\n", - "- **horizon**: Timesteps per segment (excluding overlap)\n", - "- **overlap**: Additional lookahead timesteps (improves storage optimization)\n", - "- **nr_of_previous_values**: Flow history transferred between segments" - ] + "source": "## Rolling Horizon Optimization\n\nThe `optimize.rolling_horizon()` method divides the time horizon into segments that are solved sequentially:\n\n```\nFull horizon: |---------- 1344 timesteps (14 days) ----------|\n \nSegment 1: |==== 192 (2 days) ====|-- overlap --|\nSegment 2: |==== 192 (2 days) ====|-- overlap --|\nSegment 3: |==== 192 (2 days) ====|-- overlap --|\n... \n```\n\nKey parameters:\n- **horizon**: Timesteps per segment (excluding overlap)\n- **overlap**: Additional lookahead timesteps (improves storage optimization)\n- **nr_of_previous_values**: Flow history transferred between segments" }, { "cell_type": "code", + "execution_count": null, "id": "9", "metadata": {}, + "outputs": [], "source": [ "start = timeit.default_timer()\n", "fs_rolling = flow_system.copy()\n", "segments = fs_rolling.optimize.rolling_horizon(\n", " solver,\n", - " horizon=24, # Daily segments\n", - " overlap=6, # 6-hour lookahead\n", + " horizon=192, # 2-day segments (192 timesteps at 15-min resolution)\n", + " overlap=48, # 12-hour lookahead\n", ")\n", "time_rolling = timeit.default_timer() - start\n", "\n", "print(f'Rolling horizon: {time_rolling:.2f} seconds ({len(segments)} segments)')\n", - "print(f'Cost: {fs_rolling.solution[\"costs\"].item():.2f} €')" - ], - "outputs": [], - "execution_count": null + "print(f'Cost: {fs_rolling.solution[\"costs\"].item():,.0f} €')" + ] }, { "cell_type": "markdown", @@ -217,8 +263,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "11", "metadata": {}, + "outputs": [], "source": [ "cost_full = fs_full.solution['costs'].item()\n", "cost_rolling = fs_rolling.solution['costs'].item()\n", @@ -233,10 +281,8 @@ " }\n", ").set_index('Method')\n", "\n", - "results.round(2)" - ], - "outputs": [], - "execution_count": null + "results.style.format({'Time [s]': '{:.2f}', 'Cost [€]': '{:,.0f}', 'Cost Gap [%]': '{:.2f}'})" + ] }, { "cell_type": "markdown", @@ -248,23 +294,23 @@ }, { "cell_type": "code", + "execution_count": null, "id": "13", "metadata": {}, + "outputs": [], "source": [ "fs_full.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Full)')" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "id": "14", "metadata": {}, + "outputs": [], "source": [ "fs_rolling.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Rolling)')" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -278,31 +324,28 @@ }, { "cell_type": "code", + "execution_count": null, "id": "16", "metadata": {}, + "outputs": [], "source": [ - "import plotly.graph_objects as go\n", - "from plotly.subplots import make_subplots\n", - "\n", "fig = make_subplots(\n", " rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.1, subplot_titles=['Full Optimization', 'Rolling Horizon']\n", ")\n", "\n", "# Full optimization\n", - "charge_full = fs_full.solution['Storage|charge_state'].values[:-1] # Drop final NaN\n", + "charge_full = fs_full.solution['Storage|charge_state'].values[:-1] # Drop final value\n", "fig.add_trace(go.Scatter(x=timesteps, y=charge_full, name='Full', line=dict(color='blue')), row=1, col=1)\n", "\n", "# Rolling horizon\n", "charge_rolling = fs_rolling.solution['Storage|charge_state'].values[:-1]\n", "fig.add_trace(go.Scatter(x=timesteps, y=charge_rolling, name='Rolling', line=dict(color='orange')), row=2, col=1)\n", "\n", - "fig.update_yaxes(title_text='Charge State [kWh]', row=1, col=1)\n", - "fig.update_yaxes(title_text='Charge State [kWh]', row=2, col=1)\n", + "fig.update_yaxes(title_text='Charge State [MWh]', row=1, col=1)\n", + "fig.update_yaxes(title_text='Charge State [MWh]', row=2, col=1)\n", "fig.update_layout(height=400, showlegend=False)\n", "fig.show()" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -316,8 +359,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "18", "metadata": {}, + "outputs": [], "source": [ "print(f'Number of segments: {len(segments)}')\n", "print()\n", @@ -325,10 +370,10 @@ " start_time = seg.timesteps[0]\n", " end_time = seg.timesteps[-1]\n", " cost = seg.solution['costs'].item()\n", - " print(f'Segment {i + 1}: {start_time} → {end_time} | Cost: {cost:.2f} €')" - ], - "outputs": [], - "execution_count": null + " print(\n", + " f'Segment {i + 1}: {start_time.strftime(\"%Y-%m-%d %H:%M\")} → {end_time.strftime(\"%Y-%m-%d %H:%M\")} | Cost: {cost:,.0f} €'\n", + " )" + ] }, { "cell_type": "markdown", @@ -342,25 +387,27 @@ }, { "cell_type": "code", + "execution_count": null, "id": "20", "metadata": {}, + "outputs": [], "source": [ - "overlaps = [0, 3, 6, 12, 24]\n", + "overlaps = [0, 24, 48, 96] # 0, 6h, 12h, 24h lookahead\n", "overlap_results = []\n", "\n", "for overlap in overlaps:\n", " fs = flow_system.copy()\n", " start = timeit.default_timer()\n", - " fs.optimize.rolling_horizon(solver, horizon=24, overlap=overlap)\n", + " fs.optimize.rolling_horizon(solver, horizon=192, overlap=overlap)\n", " elapsed = timeit.default_timer() - start\n", " cost = fs.solution['costs'].item()\n", " gap = (cost - cost_full) / cost_full * 100\n", - " overlap_results.append({'Overlap [h]': overlap, 'Time [s]': elapsed, 'Cost [€]': cost, 'Gap [%]': gap})\n", + " overlap_results.append(\n", + " {'Overlap': f'{overlap} ({overlap * 15 / 60:.0f}h)', 'Time [s]': elapsed, 'Cost [€]': cost, 'Gap [%]': gap}\n", + " )\n", "\n", - "pd.DataFrame(overlap_results).round(2)" - ], - "outputs": [], - "execution_count": null + "pd.DataFrame(overlap_results).style.format({'Time [s]': '{:.2f}', 'Cost [€]': '{:,.0f}', 'Gap [%]': '{:.2f}'})" + ] }, { "cell_type": "markdown", @@ -393,8 +440,8 @@ "```python\n", "segments = flow_system.optimize.rolling_horizon(\n", " solver, # Solver instance\n", - " horizon=100, # Timesteps per segment\n", - " overlap=0, # Additional lookahead timesteps\n", + " horizon=192, # Timesteps per segment (e.g., 2 days at 15-min resolution)\n", + " overlap=48, # Additional lookahead timesteps (e.g., 12 hours)\n", " nr_of_previous_values=1, # Flow history for uptime/downtime tracking\n", ")\n", "\n", From a2288d8f7eb489adfdd4d1b0b663902fca4d6b9d Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:41:28 +0100 Subject: [PATCH 13/19] Speed up notebook --- docs/notebooks/08b-rolling-horizon.ipynb | 152 +++++++++++++---------- 1 file changed, 87 insertions(+), 65 deletions(-) diff --git a/docs/notebooks/08b-rolling-horizon.ipynb b/docs/notebooks/08b-rolling-horizon.ipynb index fae589ff6..cb99fd7cf 100644 --- a/docs/notebooks/08b-rolling-horizon.ipynb +++ b/docs/notebooks/08b-rolling-horizon.ipynb @@ -16,14 +16,24 @@ }, { "cell_type": "code", - "execution_count": 3, "id": "2", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:35:42.432171Z", - "start_time": "2025-12-13T18:35:42.279884Z" + "end_time": "2025-12-13T18:42:50.619667Z", + "start_time": "2025-12-13T18:42:36.115178Z" } }, + "source": [ + "import timeit\n", + "\n", + "import pandas as pd\n", + "import plotly.graph_objects as go\n", + "from plotly.subplots import make_subplots\n", + "\n", + "import flixopt as fx\n", + "\n", + "fx.CONFIG.notebook()" + ], "outputs": [ { "data": { @@ -31,22 +41,12 @@ "flixopt.config.CONFIG" ] }, - "execution_count": 3, + "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], - "source": [ - "import timeit\n", - "\n", - "import pandas as pd\n", - "import plotly.graph_objects as go\n", - "from plotly.subplots import make_subplots\n", - "\n", - "import flixopt as fx\n", - "\n", - "fx.CONFIG.notebook()" - ] + "execution_count": 1 }, { "cell_type": "markdown", @@ -56,31 +56,13 @@ }, { "cell_type": "code", - "execution_count": 4, "id": "4", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:35:43.695556Z", - "start_time": "2025-12-13T18:35:42.878212Z" + "end_time": "2025-12-13T18:42:51.475411Z", + "start_time": "2025-12-13T18:42:51.081258Z" } }, - "outputs": [ - { - "ename": "KeyError", - "evalue": "'Value based partial slicing on non-monotonic DatetimeIndexes with non-existing keys is not allowed.'", - "output_type": "error", - "traceback": [ - "\u001b[31m---------------------------------------------------------------------------\u001b[39m", - "\u001b[31mKeyError\u001b[39m Traceback (most recent call last)", - "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[4]\u001b[39m\u001b[32m, line 9\u001b[39m\n\u001b[32m 6\u001b[39m notebook_dir = pathlib.Path(\u001b[34m__file__\u001b[39m).parent / \u001b[33m'\u001b[39m\u001b[33mdata\u001b[39m\u001b[33m'\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[33m'\u001b[39m\u001b[33m__file__\u001b[39m\u001b[33m'\u001b[39m \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mdir\u001b[39m() \u001b[38;5;28;01melse\u001b[39;00m pathlib.Path(\u001b[33m'\u001b[39m\u001b[33mdocs/notebooks/data\u001b[39m\u001b[33m'\u001b[39m)\n\u001b[32m 8\u001b[39m data = pd.read_csv(notebook_dir / \u001b[33m'\u001b[39m\u001b[33mZeitreihen2020.csv\u001b[39m\u001b[33m'\u001b[39m, index_col=\u001b[32m0\u001b[39m, parse_dates=\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[32m----> \u001b[39m\u001b[32m9\u001b[39m data = \u001b[43mdata\u001b[49m\u001b[43m[\u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43m2020-01-01\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m:\u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43m2020-01-14 23:45:00\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m]\u001b[49m \u001b[38;5;66;03m# Two weeks\u001b[39;00m\n\u001b[32m 11\u001b[39m timesteps = data.index\n\u001b[32m 13\u001b[39m \u001b[38;5;66;03m# Extract profiles\u001b[39;00m\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/PycharmProjects/flixopt_719231/.venv/lib/python3.11/site-packages/pandas/core/frame.py:4096\u001b[39m, in \u001b[36mDataFrame.__getitem__\u001b[39m\u001b[34m(self, key)\u001b[39m\n\u001b[32m 4094\u001b[39m \u001b[38;5;66;03m# Do we have a slicer (on rows)?\u001b[39;00m\n\u001b[32m 4095\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(key, \u001b[38;5;28mslice\u001b[39m):\n\u001b[32m-> \u001b[39m\u001b[32m4096\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_getitem_slice\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 4098\u001b[39m \u001b[38;5;66;03m# Do we have a (boolean) DataFrame?\u001b[39;00m\n\u001b[32m 4099\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(key, DataFrame):\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/PycharmProjects/flixopt_719231/.venv/lib/python3.11/site-packages/pandas/core/generic.py:4371\u001b[39m, in \u001b[36mNDFrame._getitem_slice\u001b[39m\u001b[34m(self, key)\u001b[39m\n\u001b[32m 4366\u001b[39m \u001b[38;5;250m\u001b[39m\u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 4367\u001b[39m \u001b[33;03m__getitem__ for the case where the key is a slice object.\u001b[39;00m\n\u001b[32m 4368\u001b[39m \u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 4369\u001b[39m \u001b[38;5;66;03m# _convert_slice_indexer to determine if this slice is positional\u001b[39;00m\n\u001b[32m 4370\u001b[39m \u001b[38;5;66;03m# or label based, and if the latter, convert to positional\u001b[39;00m\n\u001b[32m-> \u001b[39m\u001b[32m4371\u001b[39m slobj = \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mindex\u001b[49m\u001b[43m.\u001b[49m\u001b[43m_convert_slice_indexer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkind\u001b[49m\u001b[43m=\u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mgetitem\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[32m 4372\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(slobj, np.ndarray):\n\u001b[32m 4373\u001b[39m \u001b[38;5;66;03m# reachable with DatetimeIndex\u001b[39;00m\n\u001b[32m 4374\u001b[39m indexer = lib.maybe_indices_to_slice(\n\u001b[32m 4375\u001b[39m slobj.astype(np.intp, copy=\u001b[38;5;28;01mFalse\u001b[39;00m), \u001b[38;5;28mlen\u001b[39m(\u001b[38;5;28mself\u001b[39m)\n\u001b[32m 4376\u001b[39m )\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/PycharmProjects/flixopt_719231/.venv/lib/python3.11/site-packages/pandas/core/indexes/base.py:4288\u001b[39m, in \u001b[36mIndex._convert_slice_indexer\u001b[39m\u001b[34m(self, key, kind)\u001b[39m\n\u001b[32m 4286\u001b[39m indexer = key\n\u001b[32m 4287\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m-> \u001b[39m\u001b[32m4288\u001b[39m indexer = \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mslice_indexer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstart\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstop\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstep\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 4290\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m indexer\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/PycharmProjects/flixopt_719231/.venv/lib/python3.11/site-packages/pandas/core/indexes/datetimes.py:697\u001b[39m, in \u001b[36mDatetimeIndex.slice_indexer\u001b[39m\u001b[34m(self, start, end, step)\u001b[39m\n\u001b[32m 694\u001b[39m in_index &= (end_casted == \u001b[38;5;28mself\u001b[39m).any()\n\u001b[32m 696\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m in_index:\n\u001b[32m--> \u001b[39m\u001b[32m697\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m(\n\u001b[32m 698\u001b[39m \u001b[33m\"\u001b[39m\u001b[33mValue based partial slicing on non-monotonic DatetimeIndexes \u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 699\u001b[39m \u001b[33m\"\u001b[39m\u001b[33mwith non-existing keys is not allowed.\u001b[39m\u001b[33m\"\u001b[39m,\n\u001b[32m 700\u001b[39m )\n\u001b[32m 701\u001b[39m indexer = mask.nonzero()[\u001b[32m0\u001b[39m][::step]\n\u001b[32m 702\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(indexer) == \u001b[38;5;28mlen\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n", - "\u001b[31mKeyError\u001b[39m: 'Value based partial slicing on non-monotonic DatetimeIndexes with non-existing keys is not allowed.'" - ] - } - ], "source": [ "import pathlib\n", "\n", @@ -105,14 +87,29 @@ "print(f'Timesteps: {len(timesteps)} ({len(timesteps) / 96:.0f} days at 15-min resolution)')\n", "print(f'Heat demand: {heat_demand.min():.1f} - {heat_demand.max():.1f} MW')\n", "print(f'Electricity price: {electricity_price.min():.1f} - {electricity_price.max():.1f} €/MWh')" - ] + ], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Timesteps: 1344 (14 days at 15-min resolution)\n", + "Heat demand: 122.2 - 254.3 MW\n", + "Electricity price: -3.3 - 72.6 €/MWh\n" + ] + } + ], + "execution_count": 2 }, { "cell_type": "code", - "execution_count": null, "id": "5", - "metadata": {}, - "outputs": [], + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:42:52.825012Z", + "start_time": "2025-12-13T18:42:52.669724Z" + } + }, "source": [ "def build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price):\n", " \"\"\"Build a district heating system with CHP, boiler, and storage.\"\"\"\n", @@ -197,7 +194,17 @@ "\n", "flow_system = build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price)\n", "print(f'System: {len(timesteps)} timesteps')" - ] + ], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "System: 1344 timesteps\n" + ] + } + ], + "execution_count": 3 }, { "cell_type": "markdown", @@ -211,10 +218,15 @@ }, { "cell_type": "code", - "execution_count": null, "id": "7", - "metadata": {}, - "outputs": [], + "metadata": { + "jupyter": { + "is_executing": true + }, + "ExecuteTime": { + "start_time": "2025-12-13T18:42:53.822926Z" + } + }, "source": [ "solver = fx.solvers.HighsSolver()\n", "\n", @@ -225,7 +237,17 @@ "\n", "print(f'Full optimization: {time_full:.2f} seconds')\n", "print(f'Cost: {fs_full.solution[\"costs\"].item():,.0f} €')" - ] + ], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001B[2m2025-12-13 19:42:53.853\u001B[0m \u001B[33mWARNING \u001B[0m │ FlowSystem is not connected_and_transformed. Connecting and transforming data now.\n" + ] + } + ], + "execution_count": null }, { "cell_type": "markdown", @@ -235,10 +257,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "9", "metadata": {}, - "outputs": [], "source": [ "start = timeit.default_timer()\n", "fs_rolling = flow_system.copy()\n", @@ -251,7 +271,9 @@ "\n", "print(f'Rolling horizon: {time_rolling:.2f} seconds ({len(segments)} segments)')\n", "print(f'Cost: {fs_rolling.solution[\"costs\"].item():,.0f} €')" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -263,10 +285,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "11", "metadata": {}, - "outputs": [], "source": [ "cost_full = fs_full.solution['costs'].item()\n", "cost_rolling = fs_rolling.solution['costs'].item()\n", @@ -282,7 +302,9 @@ ").set_index('Method')\n", "\n", "results.style.format({'Time [s]': '{:.2f}', 'Cost [€]': '{:,.0f}', 'Cost Gap [%]': '{:.2f}'})" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -294,23 +316,23 @@ }, { "cell_type": "code", - "execution_count": null, "id": "13", "metadata": {}, - "outputs": [], "source": [ "fs_full.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Full)')" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "code", - "execution_count": null, "id": "14", "metadata": {}, - "outputs": [], "source": [ "fs_rolling.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Rolling)')" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -324,10 +346,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "16", "metadata": {}, - "outputs": [], "source": [ "fig = make_subplots(\n", " rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.1, subplot_titles=['Full Optimization', 'Rolling Horizon']\n", @@ -345,7 +365,9 @@ "fig.update_yaxes(title_text='Charge State [MWh]', row=2, col=1)\n", "fig.update_layout(height=400, showlegend=False)\n", "fig.show()" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -359,10 +381,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "18", "metadata": {}, - "outputs": [], "source": [ "print(f'Number of segments: {len(segments)}')\n", "print()\n", @@ -373,7 +393,9 @@ " print(\n", " f'Segment {i + 1}: {start_time.strftime(\"%Y-%m-%d %H:%M\")} → {end_time.strftime(\"%Y-%m-%d %H:%M\")} | Cost: {cost:,.0f} €'\n", " )" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -387,10 +409,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "20", "metadata": {}, - "outputs": [], "source": [ "overlaps = [0, 24, 48, 96] # 0, 6h, 12h, 24h lookahead\n", "overlap_results = []\n", @@ -407,7 +427,9 @@ " )\n", "\n", "pd.DataFrame(overlap_results).style.format({'Time [s]': '{:.2f}', 'Cost [€]': '{:,.0f}', 'Gap [%]': '{:.2f}'})" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", From 5b5a4aea4d44283edc1fcbd99766586b1d1a6654 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 19:59:04 +0100 Subject: [PATCH 14/19] Improve notebook --- docs/notebooks/08b-rolling-horizon.ipynb | 5091 +++++++++++++++++++++- 1 file changed, 4993 insertions(+), 98 deletions(-) diff --git a/docs/notebooks/08b-rolling-horizon.ipynb b/docs/notebooks/08b-rolling-horizon.ipynb index cb99fd7cf..4d1d019e9 100644 --- a/docs/notebooks/08b-rolling-horizon.ipynb +++ b/docs/notebooks/08b-rolling-horizon.ipynb @@ -16,24 +16,14 @@ }, { "cell_type": "code", + "execution_count": 1, "id": "2", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:42:50.619667Z", - "start_time": "2025-12-13T18:42:36.115178Z" + "end_time": "2025-12-13T18:50:32.137761Z", + "start_time": "2025-12-13T18:50:18.627145Z" } }, - "source": [ - "import timeit\n", - "\n", - "import pandas as pd\n", - "import plotly.graph_objects as go\n", - "from plotly.subplots import make_subplots\n", - "\n", - "import flixopt as fx\n", - "\n", - "fx.CONFIG.notebook()" - ], "outputs": [ { "data": { @@ -46,7 +36,19 @@ "output_type": "execute_result" } ], - "execution_count": 1 + "source": [ + "import timeit\n", + "\n", + "import pandas as pd\n", + "import plotly.express as px\n", + "import plotly.graph_objects as go\n", + "import xarray as xr\n", + "from plotly.subplots import make_subplots\n", + "\n", + "import flixopt as fx\n", + "\n", + "fx.CONFIG.notebook()" + ] }, { "cell_type": "markdown", @@ -56,13 +58,25 @@ }, { "cell_type": "code", + "execution_count": 2, "id": "4", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:42:51.475411Z", - "start_time": "2025-12-13T18:42:51.081258Z" + "end_time": "2025-12-13T18:50:32.508210Z", + "start_time": "2025-12-13T18:50:32.280206Z" } }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Timesteps: 1344 (14 days at 15-min resolution)\n", + "Heat demand: 122.2 - 254.3 MW\n", + "Electricity price: -3.3 - 72.6 €/MWh\n" + ] + } + ], "source": [ "import pathlib\n", "\n", @@ -87,29 +101,27 @@ "print(f'Timesteps: {len(timesteps)} ({len(timesteps) / 96:.0f} days at 15-min resolution)')\n", "print(f'Heat demand: {heat_demand.min():.1f} - {heat_demand.max():.1f} MW')\n", "print(f'Electricity price: {electricity_price.min():.1f} - {electricity_price.max():.1f} €/MWh')" - ], - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Timesteps: 1344 (14 days at 15-min resolution)\n", - "Heat demand: 122.2 - 254.3 MW\n", - "Electricity price: -3.3 - 72.6 €/MWh\n" - ] - } - ], - "execution_count": 2 + ] }, { "cell_type": "code", + "execution_count": 3, "id": "5", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:42:52.825012Z", - "start_time": "2025-12-13T18:42:52.669724Z" + "end_time": "2025-12-13T18:50:32.824703Z", + "start_time": "2025-12-13T18:50:32.746993Z" } }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "System: 1344 timesteps\n" + ] + } + ], "source": [ "def build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price):\n", " \"\"\"Build a district heating system with CHP, boiler, and storage.\"\"\"\n", @@ -194,17 +206,7 @@ "\n", "flow_system = build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price)\n", "print(f'System: {len(timesteps)} timesteps')" - ], - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "System: 1344 timesteps\n" - ] - } - ], - "execution_count": 3 + ] }, { "cell_type": "markdown", @@ -218,15 +220,112 @@ }, { "cell_type": "code", + "execution_count": 4, "id": "7", "metadata": { - "jupyter": { - "is_executing": true - }, "ExecuteTime": { - "start_time": "2025-12-13T18:42:53.822926Z" + "end_time": "2025-12-13T18:50:46.537647Z", + "start_time": "2025-12-13T18:50:32.897621Z" } }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001B[2m2025-12-13 19:50:32.902\u001B[0m \u001B[33mWARNING \u001B[0m │ FlowSystem is not connected_and_transformed. Connecting and transforming data now.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Writing constraints.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 74/74 [00:01<00:00, 43.82it/s]\n", + "Writing continuous variables.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 56/56 [00:00<00:00, 361.47it/s]\n", + "Writing binary variables.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 11/11 [00:00<00:00, 177.80it/s]" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Set parameter Username\n", + "Academic license - for non-commercial use only - expires 2026-11-11\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Read LP format model from file /private/var/folders/2s/46_0tgfd5gq5kkpx7k42xr_w0000gn/T/linopy-problem-vur951cb.lp\n", + "Reading time = 0.17 seconds\n", + "obj: 53792 rows, 51102 columns, 168036 nonzeros\n", + "Set parameter MIPGap to value 0.01\n", + "Set parameter TimeLimit to value 300\n", + "Set parameter LogToConsole to value 1\n", + "Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (mac64[arm] - Darwin 25.0.0 25A362)\n", + "\n", + "CPU model: Apple M1\n", + "Thread count: 8 physical cores, 8 logical processors, using up to 8 threads\n", + "\n", + "Non-default parameters:\n", + "TimeLimit 300\n", + "MIPGap 0.01\n", + "\n", + "Optimize a model with 53792 rows, 51102 columns and 168036 nonzeros\n", + "Model fingerprint: 0x6fccf55d\n", + "Variable types: 36318 continuous, 14784 integer (14784 binary)\n", + "Coefficient statistics:\n", + " Matrix range [1e-05, 2e+04]\n", + " Objective range [1e+00, 1e+00]\n", + " Bounds range [1e+00, 1e+03]\n", + " RHS range [1e-05, 2e+02]\n", + "Presolve removed 35869 rows and 35378 columns\n", + "Presolve time: 0.77s\n", + "Presolved: 17923 rows, 15724 columns, 51909 nonzeros\n", + "Variable types: 5461 continuous, 10263 integer (10263 binary)\n", + "Found heuristic solution: objective 1665125.5429\n", + "\n", + "Root relaxation: objective 1.540085e+06, 10522 iterations, 0.56 seconds (0.17 work units)\n", + "\n", + " Nodes | Current Node | Objective Bounds | Work\n", + " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", + "\n", + " 0 0 1540084.73 0 2000 1665125.54 1540084.73 7.51% - 3s\n", + "H 0 0 1541048.7119 1540084.73 0.06% - 3s\n", + "\n", + "Explored 1 nodes (13246 simplex iterations) in 3.69 seconds (1.51 work units)\n", + "Thread count was 8 (of 8 available processors)\n", + "\n", + "Solution count 2: 1.54105e+06 1.66513e+06 \n", + "\n", + "Optimal solution found (tolerance 1.00e-02)\n", + "Best objective 1.541048711885e+06, best bound 1.540084733469e+06, gap 0.0626%\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Dual values of MILP couldn't be parsed\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Full optimization: 13.62 seconds\n", + "Cost: 1,541,049 €\n" + ] + } + ], "source": [ "solver = fx.solvers.HighsSolver()\n", "\n", @@ -237,17 +336,7 @@ "\n", "print(f'Full optimization: {time_full:.2f} seconds')\n", "print(f'Cost: {fs_full.solution[\"costs\"].item():,.0f} €')" - ], - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\u001B[2m2025-12-13 19:42:53.853\u001B[0m \u001B[33mWARNING \u001B[0m │ FlowSystem is not connected_and_transformed. Connecting and transforming data now.\n" - ] - } - ], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -258,22 +347,648 @@ { "cell_type": "code", "id": "9", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:53:48.870902Z", + "start_time": "2025-12-13T18:53:19.360210Z" + } + }, "source": [ "start = timeit.default_timer()\n", "fs_rolling = flow_system.copy()\n", "segments = fs_rolling.optimize.rolling_horizon(\n", " solver,\n", " horizon=192, # 2-day segments (192 timesteps at 15-min resolution)\n", - " overlap=48, # 12-hour lookahead\n", + " overlap=96, # 1-day lookahead\n", ")\n", "time_rolling = timeit.default_timer() - start\n", "\n", "print(f'Rolling horizon: {time_rolling:.2f} seconds ({len(segments)} segments)')\n", "print(f'Cost: {fs_rolling.solution[\"costs\"].item():,.0f} €')" ], - "outputs": [], - "execution_count": null + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Segment 1/6 (timesteps 0-384): 0%| | 0/6 [00:00" + ], + "text/html": [ + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
 Time [s]Cost [€]Cost Gap [%]
Method   
Full optimization13.621,541,0490.00
Rolling horizon29.491,541,2020.01
\n" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "execution_count": 15 }, { "cell_type": "markdown", @@ -317,22 +1082,4004 @@ { "cell_type": "code", "id": "13", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:51:39.718316Z", + "start_time": "2025-12-13T18:51:39.273832Z" + } + }, "source": [ "fs_full.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Full)')" ], - "outputs": [], - "execution_count": null + "outputs": [ + { + "data": { + "text/html": [ + " \n", + " \n", + " " + ] + }, + "metadata": {}, + "output_type": "display_data", + "jetTransient": { + "display_id": null + } + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data", + "jetTransient": { + "display_id": null + } + } + ], + "execution_count": 7 }, { "cell_type": "code", "id": "14", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:51:43.807912Z", + "start_time": "2025-12-13T18:51:43.428483Z" + } + }, "source": [ "fs_rolling.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Rolling)')" ], - "outputs": [], - "execution_count": null + "outputs": [ + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data", + "jetTransient": { + "display_id": null + } + } + ], + "execution_count": 8 }, { "cell_type": "markdown", @@ -347,7 +5094,12 @@ { "cell_type": "code", "id": "16", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:53:58.854957Z", + "start_time": "2025-12-13T18:53:58.669213Z" + } + }, "source": [ "fig = make_subplots(\n", " rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.1, subplot_titles=['Full Optimization', 'Rolling Horizon']\n", @@ -366,8 +5118,45 @@ "fig.update_layout(height=400, showlegend=False)\n", "fig.show()" ], - "outputs": [], - "execution_count": null + "outputs": [ + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data", + "jetTransient": { + "display_id": null + } + } + ], + "execution_count": 16 }, { "cell_type": "markdown", @@ -382,7 +5171,12 @@ { "cell_type": "code", "id": "18", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:54:02.881690Z", + "start_time": "2025-12-13T18:54:02.848506Z" + } + }, "source": [ "print(f'Number of segments: {len(segments)}')\n", "print()\n", @@ -394,42 +5188,143 @@ " f'Segment {i + 1}: {start_time.strftime(\"%Y-%m-%d %H:%M\")} → {end_time.strftime(\"%Y-%m-%d %H:%M\")} | Cost: {cost:,.0f} €'\n", " )" ], - "outputs": [], - "execution_count": null + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of segments: 6\n", + "\n", + "Segment 1: 2020-01-01 00:00 → 2020-01-04 23:45 | Cost: 396,333 €\n", + "Segment 2: 2020-01-03 00:00 → 2020-01-06 23:45 | Cost: 407,184 €\n", + "Segment 3: 2020-01-05 00:00 → 2020-01-08 23:45 | Cost: 466,360 €\n", + "Segment 4: 2020-01-07 00:00 → 2020-01-10 23:45 | Cost: 544,281 €\n", + "Segment 5: 2020-01-09 00:00 → 2020-01-12 23:45 | Cost: 420,676 €\n", + "Segment 6: 2020-01-11 00:00 → 2020-01-14 23:45 | Cost: 412,611 €\n" + ] + } + ], + "execution_count": 17 }, { "cell_type": "markdown", "id": "19", "metadata": {}, - "source": [ - "## Effect of Overlap\n", - "\n", - "The overlap parameter provides lookahead for storage optimization. Let's compare different overlap values:" - ] + "source": "## Visualize Segment Overlaps\n\nUnderstanding how segments overlap is key to tuning rolling horizon. Let's visualize the flow rates from each segment including their overlap regions:" }, { "cell_type": "code", "id": "20", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:55:25.578028Z", + "start_time": "2025-12-13T18:55:25.350305Z" + } + }, "source": [ - "overlaps = [0, 24, 48, 96] # 0, 6h, 12h, 24h lookahead\n", - "overlap_results = []\n", - "\n", - "for overlap in overlaps:\n", - " fs = flow_system.copy()\n", - " start = timeit.default_timer()\n", - " fs.optimize.rolling_horizon(solver, horizon=192, overlap=overlap)\n", - " elapsed = timeit.default_timer() - start\n", - " cost = fs.solution['costs'].item()\n", - " gap = (cost - cost_full) / cost_full * 100\n", - " overlap_results.append(\n", - " {'Overlap': f'{overlap} ({overlap * 15 / 60:.0f}h)', 'Time [s]': elapsed, 'Cost [€]': cost, 'Gap [%]': gap}\n", - " )\n", + "# Concatenate all segment solutions into one dataset (including overlaps)\n", + "ds = xr.concat([seg.solution for seg in segments], dim=pd.RangeIndex(len(segments), name='segment'), join='outer')\n", "\n", - "pd.DataFrame(overlap_results).style.format({'Time [s]': '{:.2f}', 'Cost [€]': '{:,.0f}', 'Gap [%]': '{:.2f}'})" + "# Plot CHP thermal flow across all segments - each segment as a separate line\n", + "px.line(\n", + " ds['Boiler(Q_th)|flow_rate'].to_pandas().T,\n", + " labels={'value': 'Boiler Thermal Output [MW]', 'index': 'Timestep'},\n", + ")" + ], + "outputs": [ + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data", + "jetTransient": { + "display_id": null + } + } + ], + "execution_count": 21 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T18:58:34.250178Z", + "start_time": "2025-12-13T18:58:34.172449Z" + } + }, + "cell_type": "code", + "source": [ + "px.line(\n", + " ds['Storage|charge_state'].to_pandas().T,\n", + " labels={'value': 'Storage Charge State [MW]', 'index': 'Timestep'},\n", + ")" + ], + "id": "d7c660381f2190e0", + "outputs": [ + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data", + "jetTransient": { + "display_id": null + } + } ], - "outputs": [], - "execution_count": null + "execution_count": 27 }, { "cell_type": "markdown", From ba3a05855b458fa9906a10f34e4dc27e605c495d Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Sat, 13 Dec 2025 20:04:42 +0100 Subject: [PATCH 15/19] Improve notebook --- docs/notebooks/08b-rolling-horizon.ipynb | 918 ++++++----------------- 1 file changed, 217 insertions(+), 701 deletions(-) diff --git a/docs/notebooks/08b-rolling-horizon.ipynb b/docs/notebooks/08b-rolling-horizon.ipynb index 4d1d019e9..191f4bc7c 100644 --- a/docs/notebooks/08b-rolling-horizon.ipynb +++ b/docs/notebooks/08b-rolling-horizon.ipynb @@ -16,26 +16,13 @@ }, { "cell_type": "code", - "execution_count": 1, "id": "2", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:50:32.137761Z", - "start_time": "2025-12-13T18:50:18.627145Z" + "end_time": "2025-12-13T19:01:44.873555Z", + "start_time": "2025-12-13T19:01:40.936227Z" } }, - "outputs": [ - { - "data": { - "text/plain": [ - "flixopt.config.CONFIG" - ] - }, - "execution_count": 1, - "metadata": {}, - "output_type": "execute_result" - } - ], "source": [ "import timeit\n", "\n", @@ -48,7 +35,20 @@ "import flixopt as fx\n", "\n", "fx.CONFIG.notebook()" - ] + ], + "outputs": [ + { + "data": { + "text/plain": [ + "flixopt.config.CONFIG" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "execution_count": 1 }, { "cell_type": "markdown", @@ -58,25 +58,13 @@ }, { "cell_type": "code", - "execution_count": 2, "id": "4", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:50:32.508210Z", - "start_time": "2025-12-13T18:50:32.280206Z" + "end_time": "2025-12-13T19:01:45.078418Z", + "start_time": "2025-12-13T19:01:44.973157Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Timesteps: 1344 (14 days at 15-min resolution)\n", - "Heat demand: 122.2 - 254.3 MW\n", - "Electricity price: -3.3 - 72.6 €/MWh\n" - ] - } - ], "source": [ "import pathlib\n", "\n", @@ -101,27 +89,29 @@ "print(f'Timesteps: {len(timesteps)} ({len(timesteps) / 96:.0f} days at 15-min resolution)')\n", "print(f'Heat demand: {heat_demand.min():.1f} - {heat_demand.max():.1f} MW')\n", "print(f'Electricity price: {electricity_price.min():.1f} - {electricity_price.max():.1f} €/MWh')" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "5", - "metadata": { - "ExecuteTime": { - "end_time": "2025-12-13T18:50:32.824703Z", - "start_time": "2025-12-13T18:50:32.746993Z" - } - }, + ], "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "System: 1344 timesteps\n" + "Timesteps: 1344 (14 days at 15-min resolution)\n", + "Heat demand: 122.2 - 254.3 MW\n", + "Electricity price: -3.3 - 72.6 €/MWh\n" ] } ], + "execution_count": 2 + }, + { + "cell_type": "code", + "id": "5", + "metadata": { + "ExecuteTime": { + "end_time": "2025-12-13T19:01:45.204918Z", + "start_time": "2025-12-13T19:01:45.183230Z" + } + }, "source": [ "def build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price):\n", " \"\"\"Build a district heating system with CHP, boiler, and storage.\"\"\"\n", @@ -206,7 +196,17 @@ "\n", "flow_system = build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price)\n", "print(f'System: {len(timesteps)} timesteps')" - ] + ], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "System: 1344 timesteps\n" + ] + } + ], + "execution_count": 3 }, { "cell_type": "markdown", @@ -220,123 +220,109 @@ }, { "cell_type": "code", - "execution_count": 4, "id": "7", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:50:46.537647Z", - "start_time": "2025-12-13T18:50:32.897621Z" + "end_time": "2025-12-13T19:02:56.367270Z", + "start_time": "2025-12-13T19:01:45.486690Z" } }, + "source": [ + "solver = fx.solvers.HighsSolver()\n", + "\n", + "start = timeit.default_timer()\n", + "fs_full = flow_system.copy()\n", + "fs_full.optimize(solver)\n", + "time_full = timeit.default_timer() - start\n", + "\n", + "print(f'Full optimization: {time_full:.2f} seconds')\n", + "print(f'Cost: {fs_full.solution[\"costs\"].item():,.0f} €')" + ], "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "\u001B[2m2025-12-13 19:50:32.902\u001B[0m \u001B[33mWARNING \u001B[0m │ FlowSystem is not connected_and_transformed. Connecting and transforming data now.\n" + "\u001B[2m2025-12-13 20:01:45.496\u001B[0m \u001B[33mWARNING \u001B[0m │ FlowSystem is not connected_and_transformed. Connecting and transforming data now.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "Writing constraints.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 74/74 [00:01<00:00, 43.82it/s]\n", - "Writing continuous variables.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 56/56 [00:00<00:00, 361.47it/s]\n", - "Writing binary variables.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 11/11 [00:00<00:00, 177.80it/s]" + "Writing constraints.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 74/74 [00:00<00:00, 152.15it/s]\n", + "Writing continuous variables.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 56/56 [00:00<00:00, 378.63it/s]\n", + "Writing binary variables.: 100%|\u001B[38;2;128;191;255m██████████\u001B[0m| 11/11 [00:00<00:00, 335.35it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "Set parameter Username\n", - "Academic license - for non-commercial use only - expires 2026-11-11\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Read LP format model from file /private/var/folders/2s/46_0tgfd5gq5kkpx7k42xr_w0000gn/T/linopy-problem-vur951cb.lp\n", - "Reading time = 0.17 seconds\n", - "obj: 53792 rows, 51102 columns, 168036 nonzeros\n", - "Set parameter MIPGap to value 0.01\n", - "Set parameter TimeLimit to value 300\n", - "Set parameter LogToConsole to value 1\n", - "Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (mac64[arm] - Darwin 25.0.0 25A362)\n", + "Running HiGHS 1.12.0 (git hash: 755a8e0): Copyright (c) 2025 HiGHS under MIT licence terms\n", + "MIP linopy-problem-lqmu4n6b has 53792 rows; 51102 cols; 168036 nonzeros; 14784 integer variables (14784 binary)\n", + "Coefficient ranges:\n", + " Matrix [1e-05, 2e+04]\n", + " Cost [1e+00, 1e+00]\n", + " Bound [1e+00, 1e+03]\n", + " RHS [1e-05, 2e+02]\n", + "WARNING: Problem has some excessively small row bounds\n", + "Presolving model\n", + "29568 rows, 24168 cols, 76581 nonzeros 0s\n", + "25322 rows, 18981 cols, 72124 nonzeros 0s\n", + "24294 rows, 18173 cols, 69378 nonzeros 0s\n", + "Presolve reductions: rows 24294(-29498); columns 18173(-32929); nonzeros 69378(-98658) \n", "\n", - "CPU model: Apple M1\n", - "Thread count: 8 physical cores, 8 logical processors, using up to 8 threads\n", + "Solving MIP model with:\n", + " 24294 rows\n", + " 18173 cols (13978 binary, 0 integer, 0 implied int., 4195 continuous, 0 domain fixed)\n", + " 69378 nonzeros\n", "\n", - "Non-default parameters:\n", - "TimeLimit 300\n", - "MIPGap 0.01\n", + "Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic;\n", + " I => Shifting; J => Feasibility jump; L => Sub-MIP; P => Empty MIP; R => Randomized rounding;\n", + " S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution; Y => HiGHS solution;\n", + " Z => ZI Round; l => Trivial lower; p => Trivial point; u => Trivial upper; z => Trivial zero\n", "\n", - "Optimize a model with 53792 rows, 51102 columns and 168036 nonzeros\n", - "Model fingerprint: 0x6fccf55d\n", - "Variable types: 36318 continuous, 14784 integer (14784 binary)\n", - "Coefficient statistics:\n", - " Matrix range [1e-05, 2e+04]\n", - " Objective range [1e+00, 1e+00]\n", - " Bounds range [1e+00, 1e+03]\n", - " RHS range [1e-05, 2e+02]\n", - "Presolve removed 35869 rows and 35378 columns\n", - "Presolve time: 0.77s\n", - "Presolved: 17923 rows, 15724 columns, 51909 nonzeros\n", - "Variable types: 5461 continuous, 10263 integer (10263 binary)\n", - "Found heuristic solution: objective 1665125.5429\n", + " Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work \n", + "Src Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time\n", "\n", - "Root relaxation: objective 1.540085e+06, 10522 iterations, 0.56 seconds (0.17 work units)\n", + " 0 0 0 0.00% 1030491.76973 inf inf 0 0 0 0 0.5s\n", + " 0 0 0 0.00% 1540084.733469 inf inf 0 0 0 9609 1.0s\n", + " C 0 0 0 0.00% 1540120.790012 1662404.794591 7.36% 10533 2312 0 13317 3.5s\n", + " 0 0 0 0.00% 1540129.709896 1662404.794591 7.36% 10993 2339 0 14788 8.6s\n", + " 0 0 0 0.00% 1540135.248328 1662404.794591 7.35% 10379 2761 0 17251 13.8s\n", + " 0 0 0 0.00% 1540139.906068 1662404.794591 7.35% 10794 2560 0 18445 19.2s\n", + " 0 0 0 0.00% 1540142.257624 1662404.794591 7.35% 10289 2512 0 19682 24.4s\n", + " L 0 0 0 0.00% 1540142.371061 1591562.49514 3.23% 10213 2579 0 19922 31.8s\n", + " L 0 0 0 0.00% 1540142.371061 1540200.034522 0.00% 10213 2579 0 23652 68.0s\n", + " 1 0 1 100.00% 1540142.371061 1540200.034522 0.00% 9659 2579 0 31373 68.0s\n", "\n", - " Nodes | Current Node | Objective Bounds | Work\n", - " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", - "\n", - " 0 0 1540084.73 0 2000 1665125.54 1540084.73 7.51% - 3s\n", - "H 0 0 1541048.7119 1540084.73 0.06% - 3s\n", - "\n", - "Explored 1 nodes (13246 simplex iterations) in 3.69 seconds (1.51 work units)\n", - "Thread count was 8 (of 8 available processors)\n", - "\n", - "Solution count 2: 1.54105e+06 1.66513e+06 \n", - "\n", - "Optimal solution found (tolerance 1.00e-02)\n", - "Best objective 1.541048711885e+06, best bound 1.540084733469e+06, gap 0.0626%\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Dual values of MILP couldn't be parsed\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Full optimization: 13.62 seconds\n", - "Cost: 1,541,049 €\n" + "Solving report\n", + " Model linopy-problem-lqmu4n6b\n", + " Status Optimal\n", + " Primal bound 1540200.03452\n", + " Dual bound 1540142.37106\n", + " Gap 0.00374% (tolerance: 1%)\n", + " P-D integral 3.24885572414\n", + " Solution status feasible\n", + " 1540200.03452 (objective)\n", + " 0 (bound viol.)\n", + " 6.93576991062e-07 (int. viol.)\n", + " 0 (row viol.)\n", + " Timing 68.01\n", + " Max sub-MIP depth 2\n", + " Nodes 1\n", + " Repair LPs 0\n", + " LP iterations 31373\n", + " 0 (strong br.)\n", + " 10313 (separation)\n", + " 11450 (heuristics)\n", + "Full optimization: 70.87 seconds\n", + "Cost: 1,540,200 €\n" ] } ], - "source": [ - "solver = fx.solvers.HighsSolver()\n", - "\n", - "start = timeit.default_timer()\n", - "fs_full = flow_system.copy()\n", - "fs_full.optimize(solver)\n", - "time_full = timeit.default_timer() - start\n", - "\n", - "print(f'Full optimization: {time_full:.2f} seconds')\n", - "print(f'Cost: {fs_full.solution[\"costs\"].item():,.0f} €')" - ] + "execution_count": 4 }, { "cell_type": "markdown", @@ -349,8 +335,8 @@ "id": "9", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:53:48.870902Z", - "start_time": "2025-12-13T18:53:19.360210Z" + "end_time": "2025-12-13T19:03:27.454194Z", + "start_time": "2025-12-13T19:02:56.525964Z" } }, "source": [ @@ -371,7 +357,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Segment 1/6 (timesteps 0-384): 0%| | 0/6 [00:00" + "" ], "text/html": [ "\n", - "\n", + "
\n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", @@ -1049,27 +564,27 @@ " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", "
 Time [s]Cost [€]Cost Gap [%]Time [s]Cost [€]Cost Gap [%]
Method
Full optimization13.621,541,0490.00Full optimization70.871,540,2000.00
Rolling horizon29.491,541,2020.01Rolling horizon30.921,540,6760.03
\n" ] }, - "execution_count": 15, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], - "execution_count": 15 + "execution_count": 6 }, { "cell_type": "markdown", @@ -1084,8 +599,8 @@ "id": "13", "metadata": { "ExecuteTime": { - "end_time": "2025-12-13T18:51:39.718316Z", - "start_time": "2025-12-13T18:51:39.273832Z" + "end_time": "2025-12-13T19:03:28.570509Z", + "start_time": "2025-12-13T19:03:27.661Z" } }, "source": [ @@ -4993,9 +4508,9 @@ { "data": { "text/html": [ - "