From f6a580b4d84ddddf7c0b84e8b5bff576585f8ddf Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sat, 31 Jan 2026 23:19:50 +0100
Subject: [PATCH 01/16] fix: memory issues due to dense large coeficients
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

1. flixopt/features.py — Added sparse_multiply_sum() function that takes a sparse dict of (group_id, sum_id) -> coefficient instead of a dense DataArray. This avoids ever
  allocating the massive dense array.
  2. flixopt/elements.py — Replaced _coefficients (dense DataArray) and _flow_sign (dense DataArray) with a single _signed_coefficients cached property that returns
  dict[tuple[str, str], float | xr.DataArray] containing only non-zero signed coefficients. Updated create_linear_constraints to use sparse_multiply_sum instead of
  sparse_weighted_sum.

  The dense allocation at line 2385 (np.zeros(n_conv, max_eq, n_flows, *time) ~14.5 GB) is completely eliminated. Memory usage is now proportional to the number of non-zero
  entries (typically 2-3 flows per converter) rather than the full cartesian product.
---
 flixopt/elements.py | 148 ++++++++++++--------------------------------
 flixopt/features.py |  60 ++++++++++++++++++
 2 files changed, 101 insertions(+), 107 deletions(-)

diff --git a/flixopt/elements.py b/flixopt/elements.py
index cbdd23be1..178551edf 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -5,6 +5,7 @@
 from __future__ import annotations
 
 import logging
+from collections import defaultdict
 from functools import cached_property
 from typing import TYPE_CHECKING
 
@@ -15,7 +16,14 @@
 from . import io as fx_io
 from .config import CONFIG
 from .core import PlausibilityError
-from .features import MaskHelpers, StatusBuilder, fast_notnull, sparse_weighted_sum, stack_along_dim
+from .features import (
+    MaskHelpers,
+    StatusBuilder,
+    fast_notnull,
+    sparse_multiply_sum,
+    sparse_weighted_sum,
+    stack_along_dim,
+)
 from .interface import InvestParameters, StatusParameters
 from .modeling import ModelingUtilitiesAbstract
 from .structure import (
@@ -2278,29 +2286,6 @@ def _max_equations(self) -> int:
             return 0
         return max(len(c.conversion_factors) for c in self.converters_with_factors)
 
-    @cached_property
-    def _flow_sign(self) -> xr.DataArray:
-        """(converter, flow) sign: +1 for inputs, -1 for outputs, 0 if not involved."""
-        all_flow_ids = self._flows_model.element_ids
-
-        # Build sign array
-        sign_data = np.zeros((len(self._factor_element_ids), len(all_flow_ids)))
-        for i, conv in enumerate(self.converters_with_factors):
-            for flow in conv.inputs:
-                if flow.label_full in all_flow_ids:
-                    j = all_flow_ids.index(flow.label_full)
-                    sign_data[i, j] = 1.0  # inputs are positive
-            for flow in conv.outputs:
-                if flow.label_full in all_flow_ids:
-                    j = all_flow_ids.index(flow.label_full)
-                    sign_data[i, j] = -1.0  # outputs are negative
-
-        return xr.DataArray(
-            sign_data,
-            dims=['converter', 'flow'],
-            coords={'converter': self._factor_element_ids, 'flow': all_flow_ids},
-        )
-
     @cached_property
     def _equation_mask(self) -> xr.DataArray:
         """(converter, equation_idx) mask: 1 if equation exists, 0 otherwise."""
@@ -2318,95 +2303,45 @@ def _equation_mask(self) -> xr.DataArray:
         )
 
     @cached_property
-    def _coefficients(self) -> xr.DataArray:
-        """(converter, equation_idx, flow, [time, ...]) conversion coefficients.
+    def _signed_coefficients(self) -> dict[tuple[str, str], float | xr.DataArray]:
+        """Sparse (converter_id, flow_id) -> signed coefficient mapping.
 
-        Returns DataArray with dims (converter, equation_idx, flow) for constant coefficients,
-        or (converter, equation_idx, flow, time, ...) for time-varying coefficients.
-        Values are 0 where flow is not involved in equation.
+        Returns a dict where keys are (converter_id, flow_id) tuples and values
+        are the signed coefficients (positive for inputs, negative for outputs).
+        For converters with multiple equations, values are DataArrays with an
+        equation_idx dimension.
         """
         max_eq = self._max_equations
-        all_flow_ids = self._flows_model.element_ids
-        n_conv = len(self._factor_element_ids)
-        n_flows = len(all_flow_ids)
+        all_flow_ids_set = set(self._flows_model.element_ids)
+
+        # Collect signed coefficients per (converter, flow) across equations
+        intermediate: dict[tuple[str, str], list[tuple[int, float | xr.DataArray]]] = defaultdict(list)
 
-        # Build flow_label -> flow_id mapping for each converter
-        conv_flow_maps = []
         for conv in self.converters_with_factors:
             flow_map = {fl.label: fl.label_full for fl in conv.flows.values()}
-            conv_flow_maps.append(flow_map)
-
-        # First pass: collect all coefficients and check for time-varying
-        coeff_values = {}  # (i, eq_idx, j) -> value
-        has_dataarray = False
-        extra_coords = {}
-
-        flow_id_to_idx = {fid: j for j, fid in enumerate(all_flow_ids)}
+            # +1 for inputs, -1 for outputs
+            flow_signs = {f.label_full: 1.0 for f in conv.inputs if f.label_full in all_flow_ids_set}
+            flow_signs.update({f.label_full: -1.0 for f in conv.outputs if f.label_full in all_flow_ids_set})
 
-        for i, (conv, flow_map) in enumerate(zip(self.converters_with_factors, conv_flow_maps, strict=False)):
             for eq_idx, conv_factors in enumerate(conv.conversion_factors):
                 for flow_label, coeff in conv_factors.items():
                     flow_id = flow_map.get(flow_label)
-                    if flow_id and flow_id in flow_id_to_idx:
-                        j = flow_id_to_idx[flow_id]
-                        coeff_values[(i, eq_idx, j)] = coeff
-                        if isinstance(coeff, xr.DataArray) and coeff.ndim > 0:
-                            has_dataarray = True
-                            for d in coeff.dims:
-                                if d not in extra_coords:
-                                    extra_coords[d] = coeff.coords[d].values
-
-        # Build the coefficient array
-        if not has_dataarray:
-            # Fast path: all scalars - use simple numpy array
-            data = np.zeros((n_conv, max_eq, n_flows), dtype=np.float64)
-            for (i, eq_idx, j), val in coeff_values.items():
-                if isinstance(val, xr.DataArray):
-                    data[i, eq_idx, j] = float(val.values)
-                else:
-                    data[i, eq_idx, j] = float(val)
-
-            return xr.DataArray(
-                data,
-                dims=['converter', 'equation_idx', 'flow'],
-                coords={
-                    'converter': self._factor_element_ids,
-                    'equation_idx': list(range(max_eq)),
-                    'flow': all_flow_ids,
-                },
-            )
-        else:
-            # Slow path: some time-varying coefficients - broadcast all to common shape
-            extra_dims = list(extra_coords.keys())
-            extra_shape = [len(c) for c in extra_coords.values()]
-            full_shape = [n_conv, max_eq, n_flows] + extra_shape
-            full_dims = ['converter', 'equation_idx', 'flow'] + extra_dims
-
-            data = np.zeros(full_shape, dtype=np.float64)
-
-            # Create template for broadcasting
-            template = xr.DataArray(coords=extra_coords, dims=extra_dims) if extra_coords else None
-
-            for (i, eq_idx, j), val in coeff_values.items():
-                if isinstance(val, xr.DataArray):
-                    if val.ndim == 0:
-                        data[i, eq_idx, j, ...] = float(val.values)
-                    elif template is not None:
-                        broadcasted = val.broadcast_like(template)
-                        data[i, eq_idx, j, ...] = broadcasted.values
-                    else:
-                        data[i, eq_idx, j, ...] = val.values
-                else:
-                    data[i, eq_idx, j, ...] = float(val)
-
-            full_coords = {
-                'converter': self._factor_element_ids,
-                'equation_idx': list(range(max_eq)),
-                'flow': all_flow_ids,
-            }
-            full_coords.update(extra_coords)
-
-            return xr.DataArray(data, dims=full_dims, coords=full_coords)
+                    sign = flow_signs.get(flow_id, 0.0) if flow_id else 0.0
+                    if sign != 0.0:
+                        intermediate[(conv.label, flow_id)].append((eq_idx, coeff * sign))
+
+        # Stack each (converter, flow) pair's per-equation values into a DataArray
+        result: dict[tuple[str, str], float | xr.DataArray] = {}
+        eq_coords = list(range(max_eq))
+
+        for key, entries in intermediate.items():
+            # Build a list indexed by equation_idx (0.0 where equation doesn't use this flow)
+            per_eq: list[float | xr.DataArray] = [0.0] * max_eq
+            for eq_idx, val in entries:
+                per_eq[eq_idx] = val
+            result[key] = stack_along_dim(per_eq, dim='equation_idx', coords=eq_coords)
+
+        return result
 
     def create_linear_constraints(self) -> None:
         """Create batched linear conversion factor constraints.
@@ -2414,17 +2349,16 @@ def create_linear_constraints(self) -> None:
         For each converter c with equation i:
             sum_f(flow_rate[f] * coefficient[c,i,f] * sign[c,f]) == 0
 
-        Uses sparse_weighted_sum: each converter only touches its own 2-3 flows
-        instead of broadcasting across all flows in the system.
+        Uses sparse_multiply_sum: each converter only touches its own 2-3 flows
+        instead of allocating a dense coefficient array across all flows.
         """
         if not self.converters_with_factors:
             return
 
         flow_rate = self._flows_model[FlowVarName.RATE]
-        signed_coeffs = self._coefficients * self._flow_sign
 
         # Sparse sum: only multiplies non-zero (converter, flow) pairs
-        flow_sum = sparse_weighted_sum(flow_rate, signed_coeffs, sum_dim='flow', group_dim='converter')
+        flow_sum = sparse_multiply_sum(flow_rate, self._signed_coefficients, sum_dim='flow', group_dim='converter')
 
         # Build valid mask: True where converter HAS that equation
         n_equations_per_converter = xr.DataArray(
diff --git a/flixopt/features.py b/flixopt/features.py
index b816bd25f..5b2ba139c 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -24,6 +24,9 @@
 # =============================================================================
 
 
+Numeric = int | float | xr.DataArray
+
+
 def sparse_weighted_sum(var, coeffs: xr.DataArray, sum_dim: str, group_dim: str):
     """Compute (var * coeffs).sum(sum_dim) efficiently using sparse groupby.
 
@@ -100,6 +103,63 @@ def sparse_weighted_sum(var, coeffs: xr.DataArray, sum_dim: str, group_dim: str)
     return result.drop_vars(sum_dim, errors='ignore')
 
 
+def sparse_multiply_sum(
+    var,
+    coefficients: dict[tuple[str, str], Numeric],
+    sum_dim: str,
+    group_dim: str,
+):
+    """Compute weighted sum of var over sum_dim, grouped by group_dim, from sparse coefficients.
+
+    Unlike sparse_weighted_sum (which takes a dense DataArray and finds nonzeros),
+    this function takes an already-sparse dict of coefficients, avoiding the need
+    to ever allocate a dense array.
+
+    Args:
+        var: linopy Variable with sum_dim as a dimension.
+        coefficients: dict mapping (group_id, sum_id) to scalar or DataArray coefficient.
+            Only non-zero entries should be included.
+        sum_dim: Dimension of var to select from and sum over (e.g. 'flow').
+        group_dim: Output dimension name (e.g. 'converter').
+
+    Returns:
+        linopy expression with sum_dim removed, group_dim present.
+    """
+    if not coefficients:
+        raise ValueError('coefficients dict is empty')
+
+    # Unzip the sparse dict into parallel lists
+    group_ids_seen: dict[str, None] = {}
+    pair_group_ids: list[str] = []
+    pair_sum_ids: list[str] = []
+    pair_coeffs_list: list[Numeric] = []
+
+    for (gid, sid), coeff in coefficients.items():
+        group_ids_seen[gid] = None
+        pair_group_ids.append(gid)
+        pair_sum_ids.append(sid)
+        pair_coeffs_list.append(coeff)
+
+    group_ids = list(group_ids_seen)
+
+    # Stack mixed scalar/DataArray coefficients into a single DataArray
+    pair_coords = list(range(len(pair_group_ids)))
+    pair_coeffs = stack_along_dim(pair_coeffs_list, dim='pair', coords=pair_coords)
+
+    # Select var for active pairs, multiply by coefficients, group-sum
+    selected = var.sel({sum_dim: xr.DataArray(pair_sum_ids, dims=['pair'])})
+    weighted = selected * pair_coeffs
+
+    mapping = xr.DataArray(pair_group_ids, dims=['pair'], name=group_dim)
+    result = weighted.groupby(mapping).sum()
+
+    # Reindex to original group order (groupby sorts alphabetically)
+    result = result.sel({group_dim: group_ids})
+
+    # Drop sum_dim coord left by vectorized sel
+    return result.drop_vars(sum_dim, errors='ignore')
+
+
 def fast_notnull(arr: xr.DataArray) -> xr.DataArray:
     """Fast notnull check using numpy (~55x faster than xr.DataArray.notnull()).
 

From 5e6c8be6b7597d2e634e1b95af352d070c194deb Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 00:22:55 +0100
Subject: [PATCH 02/16] fix(effects): avoid massive memory allocation in share
 variable creation

Replace linopy.align(join='outer') with per-contributor accumulation
and linopy.merge(dim='contributor'). The old approach reindexed ALL
dimensions via xr.where(), allocating ~12.7 GB of dense arrays.
Now contributions are split by contributor at registration time and
accumulated via linopy addition (cheap for same-shape expressions),
then merged along the disjoint contributor dimension.
---
 flixopt/effects.py | 63 +++++++++++++++++++++++++++-------------------
 1 file changed, 37 insertions(+), 26 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index b9f446b4e..4dc741fff 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -345,9 +345,9 @@ def __init__(self, model: FlowSystemModel, data):
         self.share_periodic: linopy.Variable | None = None
 
         # Registered contributions from type models (FlowsModel, StoragesModel, etc.)
-        # Each entry: a defining_expr with 'contributor' dim
-        self._temporal_share_defs: list[linopy.LinearExpression] = []
-        self._periodic_share_defs: list[linopy.LinearExpression] = []
+        # Per-contributor accumulation: contributor_id -> LinearExpression (1-element contributor dim)
+        self._temporal_shares: dict[str, linopy.LinearExpression] = {}
+        self._periodic_shares: dict[str, linopy.LinearExpression] = {}
         # Constant (xr.DataArray) contributions with 'contributor' + 'effect' dims
         self._temporal_constant_defs: list[xr.DataArray] = []
         self._periodic_constant_defs: list[xr.DataArray] = []
@@ -374,7 +374,7 @@ def add_temporal_contribution(self, defining_expr, contributor_dim: str = 'contr
         if isinstance(defining_expr, xr.DataArray):
             self._temporal_constant_defs.append(defining_expr)
         else:
-            self._temporal_share_defs.append(defining_expr)
+            self._accumulate_shares(self._temporal_shares, self._as_expression(defining_expr))
 
     def add_periodic_contribution(self, defining_expr, contributor_dim: str = 'contributor') -> None:
         """Register contributors for the share|periodic variable.
@@ -389,7 +389,21 @@ def add_periodic_contribution(self, defining_expr, contributor_dim: str = 'contr
         if isinstance(defining_expr, xr.DataArray):
             self._periodic_constant_defs.append(defining_expr)
         else:
-            self._periodic_share_defs.append(defining_expr)
+            self._accumulate_shares(self._periodic_shares, self._as_expression(defining_expr))
+
+    @staticmethod
+    def _accumulate_shares(
+        accum: dict[str, linopy.LinearExpression],
+        expr: linopy.LinearExpression,
+    ) -> None:
+        """Split expression by contributor and accumulate per-contributor."""
+        for cid in expr.data.coords['contributor'].values:
+            cid_str = str(cid)
+            single = expr.sel(contributor=[cid])
+            if cid_str in accum:
+                accum[cid_str] = accum[cid_str] + single
+            else:
+                accum[cid_str] = single
 
     def create_variables(self) -> None:
         """Create batched effect variables with 'effect' dimension."""
@@ -543,8 +557,8 @@ def finalize_shares(self) -> None:
             sm.add_effect_contributions(self)
 
         # === Create share|temporal variable ===
-        if self._temporal_share_defs:
-            self.share_temporal = self._create_share_var(self._temporal_share_defs, 'share|temporal', temporal=True)
+        if self._temporal_shares:
+            self.share_temporal = self._create_share_var(self._temporal_shares, 'share|temporal', temporal=True)
             self._eq_per_timestep.lhs -= self.share_temporal.sum('contributor')
 
         # === Apply temporal constants directly ===
@@ -552,8 +566,8 @@ def finalize_shares(self) -> None:
             self._eq_per_timestep.lhs -= const.sum('contributor').reindex({'effect': self.data.effect_index})
 
         # === Create share|periodic variable ===
-        if self._periodic_share_defs:
-            self.share_periodic = self._create_share_var(self._periodic_share_defs, 'share|periodic', temporal=False)
+        if self._periodic_shares:
+            self.share_periodic = self._create_share_var(self._periodic_shares, 'share|periodic', temporal=False)
             self._eq_periodic.lhs -= self.share_periodic.sum('contributor').reindex({'effect': self.data.effect_index})
 
         # === Apply periodic constants directly ===
@@ -573,35 +587,32 @@ def _share_coords(self, element_dim: str, element_index, temporal: bool = True)
 
     def _create_share_var(
         self,
-        share_defs: list[linopy.LinearExpression],
+        shares: dict[str, linopy.LinearExpression],
         name: str,
         temporal: bool,
     ) -> linopy.Variable:
-        """Create a share variable from registered contributor definitions.
+        """Create a share variable from per-contributor expressions.
 
-        Aligns all contributor expressions (outer join on contributor dimension),
-        then sums them to produce a single expression with the full contributor dimension.
+        Each entry in shares is a single-contributor LinearExpression (already
+        accumulated). Uses linopy.merge along 'contributor' to combine them.
         """
         import pandas as pd
 
-        # Ensure all share defs have canonical effect order before alignment.
-        # linopy merge uses join="override" when shapes match, which aligns by
-        # position not label — mismatched effect order silently shuffles coefficients.
         effect_index = self.data.effect_index
-        normalized = []
-        for expr in share_defs:
+
+        # Normalize effect order
+        exprs: list[linopy.LinearExpression] = []
+        for expr in shares.values():
             if 'effect' in expr.dims:
-                expr_effects = list(expr.data.coords['effect'].values)
+                expr_effects = list(expr.coords['effect'].values)
                 if expr_effects != list(effect_index):
-                    expr = linopy.LinearExpression(expr.data.reindex(effect=effect_index), expr.model)
-            normalized.append(expr)
+                    expr = expr.reindex(effect=effect_index)
+            exprs.append(expr)
 
-        aligned = linopy.align(*normalized, join='outer', fill_value=0)
-        combined_expr = sum(aligned[1:], start=aligned[0])
+        combined_expr = linopy.merge(exprs, dim='contributor')
 
-        # Extract contributor IDs from the combined expression
-        all_ids = [str(cid) for cid in combined_expr.data.coords['contributor'].values]
-        contributor_index = pd.Index(all_ids, name='contributor')
+        # Create variable and constraint
+        contributor_index = pd.Index(list(shares.keys()), name='contributor')
         coords = self._share_coords('contributor', contributor_index, temporal=temporal)
         var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=name)
 

From 032d4eadcd87e5dc0a312321a1d240a6b7bca8c6 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 00:35:15 +0100
Subject: [PATCH 03/16] Switch to per contributor constraints to solve memmory
 issues

---
 flixopt/effects.py | 37 +++++++++++++++++++++++--------------
 tests/test_flow.py |  4 +++-
 2 files changed, 26 insertions(+), 15 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index 4dc741fff..f7d03b2f8 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -593,30 +593,39 @@ def _create_share_var(
     ) -> linopy.Variable:
         """Create a share variable from per-contributor expressions.
 
-        Each entry in shares is a single-contributor LinearExpression (already
-        accumulated). Uses linopy.merge along 'contributor' to combine them.
+        Groups contributors by nterm (number of terms) so that expressions
+        with the same shape can be merged cheaply via linopy.merge (no padding
+        needed when _term sizes match). This avoids both the massive memory
+        cost of full alignment and the overhead of per-contributor constraints.
         """
         import pandas as pd
 
         effect_index = self.data.effect_index
 
-        # Normalize effect order
-        exprs: list[linopy.LinearExpression] = []
-        for expr in shares.values():
+        contributor_index = pd.Index(list(shares.keys()), name='contributor')
+        coords = self._share_coords('contributor', contributor_index, temporal=temporal)
+        var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=name)
+
+        # Group contributors by nterm for efficient batching
+        by_nterm: dict[int, list[str]] = {}
+        normalized: dict[str, linopy.LinearExpression] = {}
+        for cid, expr in shares.items():
             if 'effect' in expr.dims:
                 expr_effects = list(expr.coords['effect'].values)
                 if expr_effects != list(effect_index):
                     expr = expr.reindex(effect=effect_index)
-            exprs.append(expr)
-
-        combined_expr = linopy.merge(exprs, dim='contributor')
-
-        # Create variable and constraint
-        contributor_index = pd.Index(list(shares.keys()), name='contributor')
-        coords = self._share_coords('contributor', contributor_index, temporal=temporal)
-        var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=name)
+            normalized[cid] = expr
+            nt = expr.nterm
+            by_nterm.setdefault(nt, []).append(cid)
+
+        for i, (_nt, cids) in enumerate(by_nterm.items()):
+            group_exprs = [normalized[cid].sel(contributor=cid) for cid in cids]
+            combined = linopy.merge(group_exprs, dim='contributor')
+            self.model.add_constraints(
+                var.sel(contributor=cids) == combined,
+                name=f'{name}|g{i}',
+            )
 
-        self.model.add_constraints(var == combined_expr, name=name)
         return var
 
     def get_periodic(self, effect_id: str) -> linopy.Variable:
diff --git a/tests/test_flow.py b/tests/test_flow.py
index 64b9bf84d..72a96a202 100644
--- a/tests/test_flow.py
+++ b/tests/test_flow.py
@@ -102,7 +102,9 @@ def test_effects_per_flow_hour(self, basic_flow_system_linopy_coords, coords_con
         model = create_linopy_model(flow_system)
 
         # Batched temporal shares are managed by the EffectsModel
-        assert 'share|temporal' in model.constraints, 'Batched temporal share constraint should exist'
+        assert any(c.startswith('share|temporal') for c in model.constraints), (
+            'Batched temporal share constraint(s) should exist'
+        )
 
         # Check batched effect variables exist
         assert 'effect|per_timestep' in model.variables, 'Batched effect per_timestep should exist'

From a9c2e22b7ae036f3ab19d1284a34d6d6dff7e94d Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 00:36:25 +0100
Subject: [PATCH 04/16] fix(effects): avoid massive memory allocation in share
 variable creation

Replace linopy.align(join='outer') with per-contributor accumulation
and individual constraints. The old approach reindexed ALL dimensions
via xr.where(), allocating ~12.7 GB of dense arrays.

Now contributions are split by contributor at registration time and
accumulated via linopy addition (cheap for same-shape expressions).
Each contributor gets its own constraint, avoiding any cross-contributor
alignment. Reduces effects expression memory from 1.2 GB to 5 MB.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 flixopt/effects.py | 20 ++++----------------
 1 file changed, 4 insertions(+), 16 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index f7d03b2f8..b0216b2f3 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -593,10 +593,8 @@ def _create_share_var(
     ) -> linopy.Variable:
         """Create a share variable from per-contributor expressions.
 
-        Groups contributors by nterm (number of terms) so that expressions
-        with the same shape can be merged cheaply via linopy.merge (no padding
-        needed when _term sizes match). This avoids both the massive memory
-        cost of full alignment and the overhead of per-contributor constraints.
+        Each entry is a single-contributor LinearExpression. Creates one
+        constraint per contributor to avoid expensive alignment/merge.
         """
         import pandas as pd
 
@@ -606,24 +604,14 @@ def _create_share_var(
         coords = self._share_coords('contributor', contributor_index, temporal=temporal)
         var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=name)
 
-        # Group contributors by nterm for efficient batching
-        by_nterm: dict[int, list[str]] = {}
-        normalized: dict[str, linopy.LinearExpression] = {}
         for cid, expr in shares.items():
             if 'effect' in expr.dims:
                 expr_effects = list(expr.coords['effect'].values)
                 if expr_effects != list(effect_index):
                     expr = expr.reindex(effect=effect_index)
-            normalized[cid] = expr
-            nt = expr.nterm
-            by_nterm.setdefault(nt, []).append(cid)
-
-        for i, (_nt, cids) in enumerate(by_nterm.items()):
-            group_exprs = [normalized[cid].sel(contributor=cid) for cid in cids]
-            combined = linopy.merge(group_exprs, dim='contributor')
             self.model.add_constraints(
-                var.sel(contributor=cids) == combined,
-                name=f'{name}|g{i}',
+                var.sel(contributor=cid) == expr.sel(contributor=cid),
+                name=f'{name}|{cid}',
             )
 
         return var

From cdc2da97010fb2814391394a3ea85d7140b196fe Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 00:47:20 +0100
Subject: [PATCH 05/16] Switch to per contributor constraints to solve memmory
 issues

---
 flixopt/effects.py | 23 +++++++++++++++++++++--
 1 file changed, 21 insertions(+), 2 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index b0216b2f3..7692725cb 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -559,7 +559,8 @@ def finalize_shares(self) -> None:
         # === Create share|temporal variable ===
         if self._temporal_shares:
             self.share_temporal = self._create_share_var(self._temporal_shares, 'share|temporal', temporal=True)
-            self._eq_per_timestep.lhs -= self.share_temporal.sum('contributor')
+            # Sum expressions directly instead of summing the variable (avoids NaN-padded where())
+            self._eq_per_timestep.lhs -= self._sum_share_expressions(self._temporal_shares)
 
         # === Apply temporal constants directly ===
         for const in self._temporal_constant_defs:
@@ -568,7 +569,7 @@ def finalize_shares(self) -> None:
         # === Create share|periodic variable ===
         if self._periodic_shares:
             self.share_periodic = self._create_share_var(self._periodic_shares, 'share|periodic', temporal=False)
-            self._eq_periodic.lhs -= self.share_periodic.sum('contributor').reindex({'effect': self.data.effect_index})
+            self._eq_periodic.lhs -= self._sum_share_expressions(self._periodic_shares)
 
         # === Apply periodic constants directly ===
         for const in self._periodic_constant_defs:
@@ -585,6 +586,24 @@ def _share_coords(self, element_dim: str, element_index, temporal: bool = True)
             }
         )
 
+    def _sum_share_expressions(self, shares: dict[str, linopy.LinearExpression]) -> linopy.LinearExpression:
+        """Sum per-contributor expressions, dropping the contributor dim.
+
+        Uses linopy.merge along _term (default), which is cheap when all
+        expressions have the same non-helper dimensions. This avoids
+        Variable.sum('contributor') which triggers expensive nansum/where().
+        """
+        effect_index = self.data.effect_index
+        flat = []
+        for cid, expr in shares.items():
+            e = expr.sel(contributor=cid)  # drop contributor dim
+            if 'effect' in e.dims:
+                expr_effects = list(e.coords['effect'].values)
+                if expr_effects != list(effect_index):
+                    e = e.reindex(effect=effect_index)
+            flat.append(e)
+        return linopy.merge(flat)  # default dim=_term = addition
+
     def _create_share_var(
         self,
         shares: dict[str, linopy.LinearExpression],

From a1d57a23bd606926ff6c5ec5c4d63fed6a7df04a Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 00:48:28 +0100
Subject: [PATCH 06/16] perf: improve bus balance to be more memmory efficient

---
 flixopt/elements.py | 22 ++++++++--------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/flixopt/elements.py b/flixopt/elements.py
index 178551edf..8f64394d5 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -1641,26 +1641,20 @@ def create_constraints(self) -> None:
         flow_dim = self._flows_model.dim_name  # 'flow'
         bus_dim = self.dim_name  # 'bus'
 
-        # Get ordered lists for coefficient matrix
         bus_ids = list(self.elements.keys())
-        flow_ids = list(flow_rate.coords[flow_dim].values)
-
-        if not bus_ids or not flow_ids:
-            logger.debug('BusesModel: no buses or flows, skipping balance constraints')
+        if not bus_ids:
+            logger.debug('BusesModel: no buses, skipping balance constraints')
             return
 
-        # Build coefficient matrix: +1 for inputs, -1 for outputs, 0 otherwise
-        coeffs = np.zeros((len(bus_ids), len(flow_ids)), dtype=np.float64)
-        for i, bus in enumerate(self.elements.values()):
+        # Build sparse coefficients: +1 for inputs, -1 for outputs
+        coefficients: dict[tuple[str, str], float] = {}
+        for bus in self.elements.values():
             for f in bus.inputs:
-                coeffs[i, flow_ids.index(f.label_full)] = 1.0
+                coefficients[(bus.label_full, f.label_full)] = 1.0
             for f in bus.outputs:
-                coeffs[i, flow_ids.index(f.label_full)] = -1.0
-
-        coeffs_da = xr.DataArray(coeffs, dims=[bus_dim, flow_dim], coords={bus_dim: bus_ids, flow_dim: flow_ids})
+                coefficients[(bus.label_full, f.label_full)] = -1.0
 
-        # Balance = sum(inputs) - sum(outputs)
-        balance = sparse_weighted_sum(flow_rate, coeffs_da, sum_dim=flow_dim, group_dim=bus_dim)
+        balance = sparse_multiply_sum(flow_rate, coefficients, sum_dim=flow_dim, group_dim=bus_dim)
 
         if self.buses_with_imbalance:
             imbalance_ids = [b.label_full for b in self.buses_with_imbalance]

From 76442c71c160e04f42f5500ce626d6bdfde0c88f Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 01:27:05 +0100
Subject: [PATCH 07/16] Switch to per effect shares

---
 flixopt/effects.py             | 114 ++++++++++++++++-----------------
 flixopt/results.py             |  10 ++-
 flixopt/statistics_accessor.py |  32 +++++----
 flixopt/transform_accessor.py  |   2 +-
 tests/test_effect.py           |   2 +-
 tests/test_flow.py             |   2 +-
 tests/test_functional.py       |   2 +-
 7 files changed, 80 insertions(+), 84 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index 7692725cb..46205ff25 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -345,9 +345,9 @@ def __init__(self, model: FlowSystemModel, data):
         self.share_periodic: linopy.Variable | None = None
 
         # Registered contributions from type models (FlowsModel, StoragesModel, etc.)
-        # Per-contributor accumulation: contributor_id -> LinearExpression (1-element contributor dim)
-        self._temporal_shares: dict[str, linopy.LinearExpression] = {}
-        self._periodic_shares: dict[str, linopy.LinearExpression] = {}
+        # Per-effect, per-contributor accumulation: effect_id -> {contributor_id -> expr (no effect dim)}
+        self._temporal_shares: dict[str, dict[str, linopy.LinearExpression]] = {}
+        self._periodic_shares: dict[str, dict[str, linopy.LinearExpression]] = {}
         # Constant (xr.DataArray) contributions with 'contributor' + 'effect' dims
         self._temporal_constant_defs: list[xr.DataArray] = []
         self._periodic_constant_defs: list[xr.DataArray] = []
@@ -393,17 +393,24 @@ def add_periodic_contribution(self, defining_expr, contributor_dim: str = 'contr
 
     @staticmethod
     def _accumulate_shares(
-        accum: dict[str, linopy.LinearExpression],
+        accum: dict[str, dict[str, linopy.LinearExpression]],
         expr: linopy.LinearExpression,
     ) -> None:
-        """Split expression by contributor and accumulate per-contributor."""
-        for cid in expr.data.coords['contributor'].values:
-            cid_str = str(cid)
-            single = expr.sel(contributor=[cid])
-            if cid_str in accum:
-                accum[cid_str] = accum[cid_str] + single
-            else:
-                accum[cid_str] = single
+        """Split expression by effect and contributor, accumulate per pair."""
+        effects = list(expr.data.coords['effect'].values) if 'effect' in expr.dims else [None]
+        for eid in effects:
+            eid_str = str(eid) if eid is not None else None
+            for cid in expr.data.coords['contributor'].values:
+                cid_str = str(cid)
+                single = expr.sel(contributor=[cid])
+                if eid is not None:
+                    single = single.sel(effect=eid, drop=True)
+                if eid_str not in accum:
+                    accum[eid_str] = {}
+                if cid_str in accum[eid_str]:
+                    accum[eid_str][cid_str] = accum[eid_str][cid_str] + single
+                else:
+                    accum[eid_str][cid_str] = single
 
     def create_variables(self) -> None:
         """Create batched effect variables with 'effect' dimension."""
@@ -556,20 +563,23 @@ def finalize_shares(self) -> None:
         if (sm := self.model._storages_model) is not None:
             sm.add_effect_contributions(self)
 
-        # === Create share|temporal variable ===
+        # === Create share|temporal variables (one per effect) ===
         if self._temporal_shares:
-            self.share_temporal = self._create_share_var(self._temporal_shares, 'share|temporal', temporal=True)
-            # Sum expressions directly instead of summing the variable (avoids NaN-padded where())
-            self._eq_per_timestep.lhs -= self._sum_share_expressions(self._temporal_shares)
+            self.share_temporal = self._create_share_vars(self._temporal_shares, 'share|temporal', temporal=True)
+            for eid, var in self.share_temporal.items():
+                summed = var.sum('contributor').expand_dims(effect=[eid])
+                self._eq_per_timestep.lhs -= summed.reindex({'effect': self.data.effect_index}, fill_value=0)
 
         # === Apply temporal constants directly ===
         for const in self._temporal_constant_defs:
             self._eq_per_timestep.lhs -= const.sum('contributor').reindex({'effect': self.data.effect_index})
 
-        # === Create share|periodic variable ===
+        # === Create share|periodic variables (one per effect) ===
         if self._periodic_shares:
-            self.share_periodic = self._create_share_var(self._periodic_shares, 'share|periodic', temporal=False)
-            self._eq_periodic.lhs -= self._sum_share_expressions(self._periodic_shares)
+            self.share_periodic = self._create_share_vars(self._periodic_shares, 'share|periodic', temporal=False)
+            for eid, var in self.share_periodic.items():
+                summed = var.sum('contributor').expand_dims(effect=[eid])
+                self._eq_periodic.lhs -= summed.reindex({'effect': self.data.effect_index}, fill_value=0)
 
         # === Apply periodic constants directly ===
         for const in self._periodic_constant_defs:
@@ -586,54 +596,44 @@ def _share_coords(self, element_dim: str, element_index, temporal: bool = True)
             }
         )
 
-    def _sum_share_expressions(self, shares: dict[str, linopy.LinearExpression]) -> linopy.LinearExpression:
-        """Sum per-contributor expressions, dropping the contributor dim.
-
-        Uses linopy.merge along _term (default), which is cheap when all
-        expressions have the same non-helper dimensions. This avoids
-        Variable.sum('contributor') which triggers expensive nansum/where().
-        """
-        effect_index = self.data.effect_index
-        flat = []
-        for cid, expr in shares.items():
-            e = expr.sel(contributor=cid)  # drop contributor dim
-            if 'effect' in e.dims:
-                expr_effects = list(e.coords['effect'].values)
-                if expr_effects != list(effect_index):
-                    e = e.reindex(effect=effect_index)
-            flat.append(e)
-        return linopy.merge(flat)  # default dim=_term = addition
-
-    def _create_share_var(
+    def _create_share_vars(
         self,
-        shares: dict[str, linopy.LinearExpression],
+        shares_by_effect: dict[str, dict[str, linopy.LinearExpression]],
         name: str,
         temporal: bool,
-    ) -> linopy.Variable:
-        """Create a share variable from per-contributor expressions.
+    ) -> dict[str, linopy.Variable]:
+        """Create one share variable per effect, each with only its contributors.
+
+        Within a single effect, all expressions have the same non-helper dims
+        (no effect dim), so linopy.merge along 'contributor' is cheap.
 
-        Each entry is a single-contributor LinearExpression. Creates one
-        constraint per contributor to avoid expensive alignment/merge.
+        Returns:
+            dict mapping effect_id -> linopy.Variable with dims (contributor, time/period).
         """
         import pandas as pd
 
-        effect_index = self.data.effect_index
-
-        contributor_index = pd.Index(list(shares.keys()), name='contributor')
-        coords = self._share_coords('contributor', contributor_index, temporal=temporal)
-        var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=name)
+        result: dict[str, linopy.Variable] = {}
 
-        for cid, expr in shares.items():
-            if 'effect' in expr.dims:
-                expr_effects = list(expr.coords['effect'].values)
-                if expr_effects != list(effect_index):
-                    expr = expr.reindex(effect=effect_index)
-            self.model.add_constraints(
-                var.sel(contributor=cid) == expr.sel(contributor=cid),
-                name=f'{name}|{cid}',
+        for eid, contributors in shares_by_effect.items():
+            contributor_index = pd.Index(list(contributors.keys()), name='contributor')
+            base_dims = None if temporal else ['period', 'scenario']
+            coords = xr.Coordinates(
+                {
+                    'contributor': contributor_index,
+                    **{k: v for k, v in (self.model.get_coords(base_dims) or {}).items()},
+                }
             )
+            var_name = f'{name}({eid})'
+            var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=var_name)
+
+            # Merge all contributor expressions into one constraint
+            exprs = [expr.sel(contributor=cid) for cid, expr in contributors.items()]
+            combined = linopy.merge(exprs, dim='contributor')
+            self.model.add_constraints(var == combined, name=var_name)
+
+            result[eid] = var
 
-        return var
+        return result
 
     def get_periodic(self, effect_id: str) -> linopy.Variable:
         """Get periodic variable for a specific effect."""
diff --git a/flixopt/results.py b/flixopt/results.py
index 66f738465..74d1e71ae 100644
--- a/flixopt/results.py
+++ b/flixopt/results.py
@@ -794,20 +794,18 @@ def get_effect_shares(
 
         ds = xr.Dataset()
 
-        share_var_name = f'share|{mode}'
+        share_var_name = f'share|{mode}({effect})'
         if share_var_name in self.solution:
             share_var = self.solution[share_var_name]
-            # Find the contributor dimension
             contributor_dim = None
             for dim in ['contributor', 'flow', 'storage', 'component', 'source']:
                 if dim in share_var.dims:
                     contributor_dim = dim
                     break
             if contributor_dim is not None and element in share_var.coords[contributor_dim].values:
-                if effect in share_var.coords['effect'].values:
-                    selected = share_var.sel({contributor_dim: element, 'effect': effect}, drop=True)
-                    label = f'{element}->{effect}({mode})'
-                    ds = xr.Dataset({label: selected})
+                selected = share_var.sel({contributor_dim: element}, drop=True)
+                label = f'{element}->{effect}({mode})'
+                ds = xr.Dataset({label: selected})
 
         if include_flows:
             if element not in self.components:
diff --git a/flixopt/statistics_accessor.py b/flixopt/statistics_accessor.py
index 50471e2d8..c1f05b210 100644
--- a/flixopt/statistics_accessor.py
+++ b/flixopt/statistics_accessor.py
@@ -805,16 +805,19 @@ def _create_effects_dataset(self, mode: Literal['temporal', 'periodic', 'total']
 
         # Determine modes to process
         modes_to_process = ['temporal', 'periodic'] if mode == 'total' else [mode]
-        share_var_map = {'temporal': 'share|temporal', 'periodic': 'share|periodic'}
-
-        # Detect contributors from batched share variables
+        # Detect contributors from per-effect share variables (share|temporal|{effect}, share|periodic|{effect})
         detected_contributors: set[str] = set()
+        share_var_names: dict[str, list[str]] = {'temporal': [], 'periodic': []}
+        for var_name in solution:
+            for mode_key in ('temporal', 'periodic'):
+                prefix = f'share|{mode_key}('
+                if str(var_name).startswith(prefix):
+                    share_var_names[mode_key].append(str(var_name))
+
         for current_mode in modes_to_process:
-            share_name = share_var_map[current_mode]
-            if share_name in solution:
+            for share_name in share_var_names[current_mode]:
                 share_da = solution[share_name]
                 for c in share_da.coords['contributor'].values:
-                    # Exclude effect-to-effect shares
                     base_name = str(c).split('(')[0] if '(' in str(c) else str(c)
                     if base_name not in effect_labels:
                         detected_contributors.add(str(c))
@@ -851,15 +854,6 @@ def get_contributor_type(contributor: str) -> str:
                 share_total: xr.DataArray | None = None
 
                 for current_mode in modes_to_process:
-                    share_name = share_var_map[current_mode]
-                    if share_name not in solution:
-                        continue
-                    share_da = solution[share_name]
-
-                    # Check if this contributor exists in the share variable
-                    if contributor not in share_da.coords['contributor'].values:
-                        continue
-
                     # Get conversion factors: which source effects contribute to this target effect
                     conversion_factors = {
                         key[0]: value
@@ -869,9 +863,13 @@ def get_contributor_type(contributor: str) -> str:
                     conversion_factors[effect] = 1  # Direct contribution
 
                     for source_effect, factor in conversion_factors.items():
-                        if source_effect not in share_da.coords['effect'].values:
+                        share_name = f'share|{current_mode}({source_effect})'
+                        if share_name not in solution:
+                            continue
+                        share_da = solution[share_name]
+                        if contributor not in share_da.coords['contributor'].values:
                             continue
-                        da = share_da.sel(contributor=contributor, effect=source_effect) * factor
+                        da = share_da.sel(contributor=contributor) * factor
                         # For total mode, sum temporal over time (apply cluster_weight for proper weighting)
                         if mode == 'total' and current_mode == 'temporal' and 'time' in da.dims:
                             weighted = da * self._fs.weights.get('cluster', 1.0)
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index b8951bf56..93c99491e 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -471,7 +471,7 @@ def _cluster_indices_per_timestep(self) -> xr.DataArray:
     def _get_mode(var_name: str) -> ExpansionMode:
         """Look up expansion mode for a variable name via suffix matching."""
         for suffix, mode in NAME_TO_EXPANSION.items():
-            if var_name == suffix or var_name.endswith(suffix):
+            if var_name == suffix or var_name.endswith(suffix) or var_name.startswith(suffix):
                 return mode
         return ExpansionMode.REPEAT
 
diff --git a/tests/test_effect.py b/tests/test_effect.py
index f26ad3438..af3ad9dc4 100644
--- a/tests/test_effect.py
+++ b/tests/test_effect.py
@@ -164,7 +164,7 @@ def test_shares(self, basic_flow_system_linopy_coords, coords_config):
 
         # Check share allocation variables exist (e.g., share|temporal_from_effect for effect-to-effect shares)
         # These are managed by the EffectsModel
-        assert 'share|temporal' in model.variables, 'Temporal share variable should exist'
+        assert any(v.startswith('share|temporal') for v in model.variables), 'Temporal share variable(s) should exist'
 
         # Access individual effect variables using batched model + sel
         _effect2_periodic = model.variables['effect|periodic'].sel(effect='Effect2')
diff --git a/tests/test_flow.py b/tests/test_flow.py
index 72a96a202..2aafb08a1 100644
--- a/tests/test_flow.py
+++ b/tests/test_flow.py
@@ -357,7 +357,7 @@ def test_flow_invest_divest_effects(self, basic_flow_system_linopy_coords, coord
         assert 'effect|periodic' in model.variables
 
         # Check that temporal share exists for the flow's effects
-        assert 'share|temporal' in model.variables
+        assert any(v.startswith('share|temporal') for v in model.variables)
 
 
 class TestFlowOnModel:
diff --git a/tests/test_functional.py b/tests/test_functional.py
index 509be5d04..341bc9ee8 100644
--- a/tests/test_functional.py
+++ b/tests/test_functional.py
@@ -129,7 +129,7 @@ def test_minimal_model(solver_fixture, time_steps_fixture):
     )
 
     assert_allclose(
-        flow_system.solution['share|temporal'].sel(contributor='Gastarif(Gas)', effect='costs').values[:-1],
+        flow_system.solution['share|temporal(costs)'].sel(contributor='Gastarif(Gas)').values[:-1],
         [-0.0, 20.0, 40.0, -0.0, 20.0],
         rtol=1e-5,
         atol=1e-10,

From 00f92c159bbd28ee956c5bb9540e2a82f6893438 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 02:43:19 +0100
Subject: [PATCH 08/16] Firs succesfull drop to 10 GB

---
 flixopt/components.py |  56 +++++++++++---
 flixopt/effects.py    | 172 ++++++++++++++++++++++++++++++++----------
 flixopt/elements.py   |  94 +++++++++++++++++++----
 3 files changed, 257 insertions(+), 65 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index e25eaf757..9dd837ace 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -917,25 +917,59 @@ def add_effect_contributions(self, effects_model) -> None:
         if inv.effects_per_size is not None:
             factors = inv.effects_per_size
             size = self.size.sel({dim: factors.coords[dim].values})
-            effects_model.add_periodic_contribution(size * factors, contributor_dim=dim)
+            for eid in factors.coords['effect'].values:
+                f_single = factors.sel(effect=eid, drop=True)
+                if (f_single == 0).all():
+                    continue
+                effects_model.add_periodic_contribution(size * f_single, contributor_dim=dim, effect=str(eid))
 
             # Investment/retirement effects
             invested = self.invested
             if invested is not None:
-                if (f := inv.effects_of_investment) is not None:
-                    effects_model.add_periodic_contribution(
-                        invested.sel({dim: f.coords[dim].values}) * f, contributor_dim=dim
-                    )
-                if (f := inv.effects_of_retirement) is not None:
-                    effects_model.add_periodic_contribution(
-                        invested.sel({dim: f.coords[dim].values}) * (-f), contributor_dim=dim
-                    )
+                if (ff := inv.effects_of_investment) is not None:
+                    for eid in ff.coords['effect'].values:
+                        f_single = ff.sel(effect=eid, drop=True)
+                        if (f_single == 0).all():
+                            continue
+                        effects_model.add_periodic_contribution(
+                            invested.sel({dim: f_single.coords[dim].values}) * f_single,
+                            contributor_dim=dim,
+                            effect=str(eid),
+                        )
+                if (ff := inv.effects_of_retirement) is not None:
+                    for eid in ff.coords['effect'].values:
+                        f_single = ff.sel(effect=eid, drop=True)
+                        if (f_single == 0).all():
+                            continue
+                        effects_model.add_periodic_contribution(
+                            invested.sel({dim: f_single.coords[dim].values}) * (-f_single),
+                            contributor_dim=dim,
+                            effect=str(eid),
+                        )
 
         # === Constants: mandatory fixed + retirement ===
         if inv.effects_of_investment_mandatory is not None:
-            effects_model.add_periodic_contribution(inv.effects_of_investment_mandatory, contributor_dim=dim)
+            mandatory = inv.effects_of_investment_mandatory
+            if 'effect' in mandatory.dims:
+                for eid in mandatory.coords['effect'].values:
+                    effects_model.add_periodic_contribution(
+                        mandatory.sel(effect=eid, drop=True),
+                        contributor_dim=dim,
+                        effect=str(eid),
+                    )
+            else:
+                effects_model.add_periodic_contribution(mandatory, contributor_dim=dim)
         if inv.effects_of_retirement_constant is not None:
-            effects_model.add_periodic_contribution(inv.effects_of_retirement_constant, contributor_dim=dim)
+            ret_const = inv.effects_of_retirement_constant
+            if 'effect' in ret_const.dims:
+                for eid in ret_const.coords['effect'].values:
+                    effects_model.add_periodic_contribution(
+                        ret_const.sel(effect=eid, drop=True),
+                        contributor_dim=dim,
+                        effect=str(eid),
+                    )
+            else:
+                effects_model.add_periodic_contribution(ret_const, contributor_dim=dim)
 
     # --- Investment Cached Properties ---
 
diff --git a/flixopt/effects.py b/flixopt/effects.py
index 46205ff25..255a33912 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -361,56 +361,68 @@ def effect_index(self):
         """Public access to the effect index for type models."""
         return self.data.effect_index
 
-    def add_temporal_contribution(self, defining_expr, contributor_dim: str = 'contributor') -> None:
+    def add_temporal_contribution(
+        self,
+        defining_expr,
+        contributor_dim: str = 'contributor',
+        effect: str | None = None,
+    ) -> None:
         """Register contributors for the share|temporal variable.
 
         Args:
-            defining_expr: Expression with a contributor dimension.
-                Accepts linopy LinearExpression/Variable or plain xr.DataArray (constants).
+            defining_expr: Expression with a contributor dimension (no effect dim if effect is given).
             contributor_dim: Name of the element dimension to rename to 'contributor'.
+            effect: If provided, the expression is for this specific effect (no effect dim needed).
         """
         if contributor_dim != 'contributor':
             defining_expr = defining_expr.rename({contributor_dim: 'contributor'})
         if isinstance(defining_expr, xr.DataArray):
+            if effect is not None:
+                defining_expr = defining_expr.expand_dims(effect=[effect])
             self._temporal_constant_defs.append(defining_expr)
         else:
-            self._accumulate_shares(self._temporal_shares, self._as_expression(defining_expr))
+            self._accumulate_shares(self._temporal_shares, self._as_expression(defining_expr), effect)
 
-    def add_periodic_contribution(self, defining_expr, contributor_dim: str = 'contributor') -> None:
+    def add_periodic_contribution(
+        self,
+        defining_expr,
+        contributor_dim: str = 'contributor',
+        effect: str | None = None,
+    ) -> None:
         """Register contributors for the share|periodic variable.
 
         Args:
-            defining_expr: Expression with a contributor dimension.
-                Accepts linopy LinearExpression/Variable or plain xr.DataArray (constants).
+            defining_expr: Expression with a contributor dimension (no effect dim if effect is given).
             contributor_dim: Name of the element dimension to rename to 'contributor'.
+            effect: If provided, the expression is for this specific effect (no effect dim needed).
         """
         if contributor_dim != 'contributor':
             defining_expr = defining_expr.rename({contributor_dim: 'contributor'})
         if isinstance(defining_expr, xr.DataArray):
+            if effect is not None:
+                defining_expr = defining_expr.expand_dims(effect=[effect])
             self._periodic_constant_defs.append(defining_expr)
         else:
-            self._accumulate_shares(self._periodic_shares, self._as_expression(defining_expr))
+            self._accumulate_shares(self._periodic_shares, self._as_expression(defining_expr), effect)
 
     @staticmethod
     def _accumulate_shares(
-        accum: dict[str, dict[str, linopy.LinearExpression]],
+        accum: dict[str, list],
         expr: linopy.LinearExpression,
+        effect: str | None = None,
     ) -> None:
-        """Split expression by effect and contributor, accumulate per pair."""
-        effects = list(expr.data.coords['effect'].values) if 'effect' in expr.dims else [None]
-        for eid in effects:
-            eid_str = str(eid) if eid is not None else None
-            for cid in expr.data.coords['contributor'].values:
-                cid_str = str(cid)
-                single = expr.sel(contributor=[cid])
-                if eid is not None:
-                    single = single.sel(effect=eid, drop=True)
-                if eid_str not in accum:
-                    accum[eid_str] = {}
-                if cid_str in accum[eid_str]:
-                    accum[eid_str][cid_str] = accum[eid_str][cid_str] + single
-                else:
-                    accum[eid_str][cid_str] = single
+        """Append expression to per-effect list."""
+        # accum structure: {effect_id: [(expr, contributor_ids), ...]}
+        if effect is not None:
+            # Expression has no effect dim — tagged with specific effect
+            accum.setdefault(effect, []).append(expr)
+        elif 'effect' in expr.dims:
+            # Expression has effect dim — split per effect (DataArray sel is cheap)
+            for eid in expr.data.coords['effect'].values:
+                eid_str = str(eid)
+                accum.setdefault(eid_str, []).append(expr.sel(effect=eid, drop=True))
+        else:
+            raise ValueError('Expression must have effect dim or effect parameter must be given')
 
     def create_variables(self) -> None:
         """Create batched effect variables with 'effect' dimension."""
@@ -598,41 +610,121 @@ def _share_coords(self, element_dim: str, element_index, temporal: bool = True)
 
     def _create_share_vars(
         self,
-        shares_by_effect: dict[str, dict[str, linopy.LinearExpression]],
+        accum: dict[str, list[linopy.LinearExpression]],
         name: str,
         temporal: bool,
     ) -> dict[str, linopy.Variable]:
-        """Create one share variable per effect, each with only its contributors.
+        """Create one share variable per effect from accumulated per-effect expression lists.
+
+        accum structure: {effect_id: [expr1, expr2, ...]} where each expr has
+        (contributor, ...other_dims) dims — no effect dim.
 
-        Within a single effect, all expressions have the same non-helper dims
-        (no effect dim), so linopy.merge along 'contributor' is cheap.
+        Builds combined constraint arrays with numpy to avoid expensive linopy.merge.
 
         Returns:
             dict mapping effect_id -> linopy.Variable with dims (contributor, time/period).
         """
+        from collections import defaultdict
+
         import pandas as pd
 
+        if not accum:
+            return {}
+
         result: dict[str, linopy.Variable] = {}
 
-        for eid, contributors in shares_by_effect.items():
-            contributor_index = pd.Index(list(contributors.keys()), name='contributor')
+        for eid, expr_list in accum.items():
+            # Determine canonical dim order from variable coords
             base_dims = None if temporal else ['period', 'scenario']
-            coords = xr.Coordinates(
-                {
-                    'contributor': contributor_index,
-                    **{k: v for k, v in (self.model.get_coords(base_dims) or {}).items()},
-                }
-            )
+            model_coords = self.model.get_coords(base_dims) or {}
+
+            # Phase 1: Extract per-contributor term slices using xarray (handles dim order)
+            # contributor_id -> list of (coeffs_np, vars_np, const_np) in canonical order
+            contrib_terms: dict[str, list[tuple]] = defaultdict(list)
+
+            # Canonical other_dim order: from model coords (excludes contributor)
+            canonical_dims = list(model_coords.keys())
+
+            for expr in expr_list:
+                ds = expr.data
+                for cid in ds.coords['contributor'].values:
+                    cid_str = str(cid)
+                    # sel contributor, then transpose to canonical order
+                    c_da = ds['coeffs'].sel(contributor=cid)
+                    v_da = ds['vars'].sel(contributor=cid)
+                    k_da = ds['const'].sel(contributor=cid)
+                    # Transpose to canonical: (...canonical_dims, _term) for c/v, (...canonical_dims) for k
+                    c_order = [d for d in canonical_dims if d in c_da.dims] + ['_term']
+                    k_order = [d for d in canonical_dims if d in k_da.dims]
+                    contrib_terms[cid_str].append(
+                        (
+                            c_da.transpose(*c_order).values,
+                            v_da.transpose(*c_order).values,
+                            k_da.transpose(*k_order).values if k_order else k_da.values,
+                        )
+                    )
+
+            # Phase 2: Build combined numpy arrays
+            contributor_ids = list(contrib_terms.keys())
+            contributor_index = pd.Index(contributor_ids, name='contributor')
+            coords = xr.Coordinates({'contributor': contributor_index, **{k: v for k, v in model_coords.items()}})
             var_name = f'{name}({eid})'
             var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=var_name)
 
-            # Merge all contributor expressions into one constraint
-            exprs = [expr.sel(contributor=cid) for cid, expr in contributors.items()]
-            combined = linopy.merge(exprs, dim='contributor')
-            self.model.add_constraints(var == combined, name=var_name)
+            # Compute max total _term across contributors
+            max_nterm = 0
+            for terms in contrib_terms.values():
+                total = sum(t[0].shape[-1] for t in terms)
+                if total > max_nterm:
+                    max_nterm = total
+
+            total_terms = 1 + max_nterm
+
+            # other_shape from canonical dims
+            sample_const = next(iter(contrib_terms.values()))[0][2]
+            other_shape = sample_const.shape
+
+            n_contrib = len(contributor_ids)
+            all_coeffs = np.zeros((n_contrib, *other_shape, total_terms))
+            all_vars = np.full((n_contrib, *other_shape, total_terms), -1, dtype=np.int64)
+            all_const = np.zeros((n_contrib, *other_shape))
+
+            var_data = var.data
+            for i, cid in enumerate(contributor_ids):
+                var_label_da = var_data['labels'].sel(contributor=cid)
+                var_label = var_label_da.transpose(*canonical_dims).values
+                all_coeffs[i, ..., 0] = 1.0
+                all_vars[i, ..., 0] = var_label
+
+                offset = 1
+                for c_slice, v_slice, k_slice in contrib_terms[cid]:
+                    nt = c_slice.shape[-1]
+                    all_coeffs[i, ..., offset : offset + nt] = -c_slice
+                    all_vars[i, ..., offset : offset + nt] = v_slice
+                    all_const[i, ...] -= k_slice
+                    offset += nt
+
+            all_dims = ['contributor'] + canonical_dims + ['_term']
+            const_dims = ['contributor'] + canonical_dims
+
+            other_coords = {d: model_coords[d] for d in canonical_dims}
+            ds_coords = {'contributor': contributor_ids, **other_coords}
+
+            con_ds = xr.Dataset(
+                {
+                    'coeffs': (all_dims, all_coeffs),
+                    'vars': (all_dims, all_vars),
+                    'const': (const_dims, all_const),
+                },
+                coords=ds_coords,
+            )
+
+            lhs = linopy.LinearExpression(con_ds, self.model)
+            self.model.add_constraints(lhs == 0, name=var_name)
 
             result[eid] = var
 
+        accum.clear()
         return result
 
     def get_periodic(self, effect_id: str) -> linopy.Variable:
diff --git a/flixopt/elements.py b/flixopt/elements.py
index 8f64394d5..5ad9be90e 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -1221,7 +1221,17 @@ def add_effect_contributions(self, effects_model) -> None:
         factors = self.data.effects_per_flow_hour
         if factors is not None:
             rate = self.rate.sel({dim: factors.coords[dim].values})
-            effects_model.add_temporal_contribution(rate * (factors * dt), contributor_dim=dim)
+            for eid in factors.coords['effect'].values:
+                f_single = factors.sel(effect=eid, drop=True)  # (flow,) or (flow, time) — pure DataArray, cheap
+                # Only include flows with nonzero factor
+                nonzero = f_single != 0
+                if not nonzero.any():
+                    continue
+                effects_model.add_temporal_contribution(
+                    rate * (f_single * dt),
+                    contributor_dim=dim,
+                    effect=str(eid),
+                )
 
         # === Temporal: status effects ===
         if self.status is not None:
@@ -1229,38 +1239,94 @@ def add_effect_contributions(self, effects_model) -> None:
             if factor is not None:
                 flow_ids = factor.coords[dim].values
                 status_subset = self.status.sel({dim: flow_ids})
-                effects_model.add_temporal_contribution(status_subset * (factor * dt), contributor_dim=dim)
+                for eid in factor.coords['effect'].values:
+                    f_single = factor.sel(effect=eid, drop=True)
+                    nonzero = f_single != 0
+                    if not nonzero.any():
+                        continue
+                    effects_model.add_temporal_contribution(
+                        status_subset * (f_single * dt),
+                        contributor_dim=dim,
+                        effect=str(eid),
+                    )
 
             factor = self.data.effects_per_startup
             if self.startup is not None and factor is not None:
                 flow_ids = factor.coords[dim].values
                 startup_subset = self.startup.sel({dim: flow_ids})
-                effects_model.add_temporal_contribution(startup_subset * factor, contributor_dim=dim)
+                for eid in factor.coords['effect'].values:
+                    f_single = factor.sel(effect=eid, drop=True)
+                    nonzero = f_single != 0
+                    if not nonzero.any():
+                        continue
+                    effects_model.add_temporal_contribution(
+                        startup_subset * f_single,
+                        contributor_dim=dim,
+                        effect=str(eid),
+                    )
 
         # === Periodic: size * effects_per_size ===
         inv = self.data._investment_data
         if inv is not None and inv.effects_per_size is not None:
             factors = inv.effects_per_size
             size = self.size.sel({dim: factors.coords[dim].values})
-            effects_model.add_periodic_contribution(size * factors, contributor_dim=dim)
+            for eid in factors.coords['effect'].values:
+                f_single = factors.sel(effect=eid, drop=True)
+                nonzero = f_single != 0
+                if not nonzero.any():
+                    continue
+                effects_model.add_periodic_contribution(size * f_single, contributor_dim=dim, effect=str(eid))
 
             # Investment/retirement effects
             if self.invested is not None:
-                if (f := inv.effects_of_investment) is not None:
-                    effects_model.add_periodic_contribution(
-                        self.invested.sel({dim: f.coords[dim].values}) * f, contributor_dim=dim
-                    )
-                if (f := inv.effects_of_retirement) is not None:
-                    effects_model.add_periodic_contribution(
-                        self.invested.sel({dim: f.coords[dim].values}) * (-f), contributor_dim=dim
-                    )
+                if (ff := inv.effects_of_investment) is not None:
+                    for eid in ff.coords['effect'].values:
+                        f_single = ff.sel(effect=eid, drop=True)
+                        nonzero = f_single != 0
+                        if not nonzero.any():
+                            continue
+                        effects_model.add_periodic_contribution(
+                            self.invested.sel({dim: f_single.coords[dim].values}) * f_single,
+                            contributor_dim=dim,
+                            effect=str(eid),
+                        )
+                if (ff := inv.effects_of_retirement) is not None:
+                    for eid in ff.coords['effect'].values:
+                        f_single = ff.sel(effect=eid, drop=True)
+                        nonzero = f_single != 0
+                        if not nonzero.any():
+                            continue
+                        effects_model.add_periodic_contribution(
+                            self.invested.sel({dim: f_single.coords[dim].values}) * (-f_single),
+                            contributor_dim=dim,
+                            effect=str(eid),
+                        )
 
         # === Constants: mandatory fixed + retirement ===
         if inv is not None:
             if inv.effects_of_investment_mandatory is not None:
-                effects_model.add_periodic_contribution(inv.effects_of_investment_mandatory, contributor_dim=dim)
+                # These already have effect dim — split per effect
+                mandatory = inv.effects_of_investment_mandatory
+                if 'effect' in mandatory.dims:
+                    for eid in mandatory.coords['effect'].values:
+                        effects_model.add_periodic_contribution(
+                            mandatory.sel(effect=eid, drop=True),
+                            contributor_dim=dim,
+                            effect=str(eid),
+                        )
+                else:
+                    effects_model.add_periodic_contribution(mandatory, contributor_dim=dim)
             if inv.effects_of_retirement_constant is not None:
-                effects_model.add_periodic_contribution(inv.effects_of_retirement_constant, contributor_dim=dim)
+                ret_const = inv.effects_of_retirement_constant
+                if 'effect' in ret_const.dims:
+                    for eid in ret_const.coords['effect'].values:
+                        effects_model.add_periodic_contribution(
+                            ret_const.sel(effect=eid, drop=True),
+                            contributor_dim=dim,
+                            effect=str(eid),
+                        )
+                else:
+                    effects_model.add_periodic_contribution(ret_const, contributor_dim=dim)
 
     # === Status Variables (cached_property) ===
 

From 3a61f5bab40fdb4b168b2987268c70f06be70477 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 02:52:58 +0100
Subject: [PATCH 09/16] Make more readable

---
 flixopt/effects.py | 100 ++++++++-------------------------------------
 1 file changed, 17 insertions(+), 83 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index 255a33912..e6d91b589 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -619,13 +619,12 @@ def _create_share_vars(
         accum structure: {effect_id: [expr1, expr2, ...]} where each expr has
         (contributor, ...other_dims) dims — no effect dim.
 
-        Builds combined constraint arrays with numpy to avoid expensive linopy.merge.
+        Within a single effect, all expressions share the same non-helper dims,
+        so linopy.merge along 'contributor' is cheap (no cross-dim alignment).
 
         Returns:
             dict mapping effect_id -> linopy.Variable with dims (contributor, time/period).
         """
-        from collections import defaultdict
-
         import pandas as pd
 
         if not accum:
@@ -634,93 +633,28 @@ def _create_share_vars(
         result: dict[str, linopy.Variable] = {}
 
         for eid, expr_list in accum.items():
-            # Determine canonical dim order from variable coords
-            base_dims = None if temporal else ['period', 'scenario']
-            model_coords = self.model.get_coords(base_dims) or {}
-
-            # Phase 1: Extract per-contributor term slices using xarray (handles dim order)
-            # contributor_id -> list of (coeffs_np, vars_np, const_np) in canonical order
-            contrib_terms: dict[str, list[tuple]] = defaultdict(list)
-
-            # Canonical other_dim order: from model coords (excludes contributor)
-            canonical_dims = list(model_coords.keys())
-
+            # Sum expressions that share the same contributors (e.g. rate + status + startup)
+            combined: dict[str, linopy.LinearExpression] = {}
             for expr in expr_list:
-                ds = expr.data
-                for cid in ds.coords['contributor'].values:
+                for cid in expr.data.coords['contributor'].values:
                     cid_str = str(cid)
-                    # sel contributor, then transpose to canonical order
-                    c_da = ds['coeffs'].sel(contributor=cid)
-                    v_da = ds['vars'].sel(contributor=cid)
-                    k_da = ds['const'].sel(contributor=cid)
-                    # Transpose to canonical: (...canonical_dims, _term) for c/v, (...canonical_dims) for k
-                    c_order = [d for d in canonical_dims if d in c_da.dims] + ['_term']
-                    k_order = [d for d in canonical_dims if d in k_da.dims]
-                    contrib_terms[cid_str].append(
-                        (
-                            c_da.transpose(*c_order).values,
-                            v_da.transpose(*c_order).values,
-                            k_da.transpose(*k_order).values if k_order else k_da.values,
-                        )
-                    )
-
-            # Phase 2: Build combined numpy arrays
-            contributor_ids = list(contrib_terms.keys())
+                    single = expr.sel(contributor=[cid])
+                    if cid_str in combined:
+                        combined[cid_str] = linopy.merge([combined[cid_str], single])
+                    else:
+                        combined[cid_str] = single
+
+            contributor_ids = list(combined.keys())
             contributor_index = pd.Index(contributor_ids, name='contributor')
+            base_dims = None if temporal else ['period', 'scenario']
+            model_coords = self.model.get_coords(base_dims) or {}
             coords = xr.Coordinates({'contributor': contributor_index, **{k: v for k, v in model_coords.items()}})
             var_name = f'{name}({eid})'
             var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=var_name)
 
-            # Compute max total _term across contributors
-            max_nterm = 0
-            for terms in contrib_terms.values():
-                total = sum(t[0].shape[-1] for t in terms)
-                if total > max_nterm:
-                    max_nterm = total
-
-            total_terms = 1 + max_nterm
-
-            # other_shape from canonical dims
-            sample_const = next(iter(contrib_terms.values()))[0][2]
-            other_shape = sample_const.shape
-
-            n_contrib = len(contributor_ids)
-            all_coeffs = np.zeros((n_contrib, *other_shape, total_terms))
-            all_vars = np.full((n_contrib, *other_shape, total_terms), -1, dtype=np.int64)
-            all_const = np.zeros((n_contrib, *other_shape))
-
-            var_data = var.data
-            for i, cid in enumerate(contributor_ids):
-                var_label_da = var_data['labels'].sel(contributor=cid)
-                var_label = var_label_da.transpose(*canonical_dims).values
-                all_coeffs[i, ..., 0] = 1.0
-                all_vars[i, ..., 0] = var_label
-
-                offset = 1
-                for c_slice, v_slice, k_slice in contrib_terms[cid]:
-                    nt = c_slice.shape[-1]
-                    all_coeffs[i, ..., offset : offset + nt] = -c_slice
-                    all_vars[i, ..., offset : offset + nt] = v_slice
-                    all_const[i, ...] -= k_slice
-                    offset += nt
-
-            all_dims = ['contributor'] + canonical_dims + ['_term']
-            const_dims = ['contributor'] + canonical_dims
-
-            other_coords = {d: model_coords[d] for d in canonical_dims}
-            ds_coords = {'contributor': contributor_ids, **other_coords}
-
-            con_ds = xr.Dataset(
-                {
-                    'coeffs': (all_dims, all_coeffs),
-                    'vars': (all_dims, all_vars),
-                    'const': (const_dims, all_const),
-                },
-                coords=ds_coords,
-            )
-
-            lhs = linopy.LinearExpression(con_ds, self.model)
-            self.model.add_constraints(lhs == 0, name=var_name)
+            # Merge all per-contributor expressions along 'contributor'
+            merged = linopy.merge(list(combined.values()), dim='contributor')
+            self.model.add_constraints(var == merged, name=var_name)
 
             result[eid] = var
 

From 0ef2a02f0052ecc92045eeb17bade792ef5cbfce Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 03:15:48 +0100
Subject: [PATCH 10/16] Go back to one variable for all shares

---
 flixopt/effects.py             | 76 ++++++++++++++++++++--------------
 flixopt/statistics_accessor.py | 29 ++++++-------
 tests/test_effect.py           |  2 +-
 tests/test_flow.py             |  4 +-
 tests/test_functional.py       |  2 +-
 5 files changed, 62 insertions(+), 51 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index e6d91b589..bf5bfa8f5 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -575,23 +575,19 @@ def finalize_shares(self) -> None:
         if (sm := self.model._storages_model) is not None:
             sm.add_effect_contributions(self)
 
-        # === Create share|temporal variables (one per effect) ===
+        # === Create share|temporal variable (one combined with contributor × effect dims) ===
         if self._temporal_shares:
-            self.share_temporal = self._create_share_vars(self._temporal_shares, 'share|temporal', temporal=True)
-            for eid, var in self.share_temporal.items():
-                summed = var.sum('contributor').expand_dims(effect=[eid])
-                self._eq_per_timestep.lhs -= summed.reindex({'effect': self.data.effect_index}, fill_value=0)
+            self.share_temporal = self._create_share_var(self._temporal_shares, 'share|temporal', temporal=True)
+            self._eq_per_timestep.lhs -= self.share_temporal.sum('contributor')
 
         # === Apply temporal constants directly ===
         for const in self._temporal_constant_defs:
             self._eq_per_timestep.lhs -= const.sum('contributor').reindex({'effect': self.data.effect_index})
 
-        # === Create share|periodic variables (one per effect) ===
+        # === Create share|periodic variable (one combined with contributor × effect dims) ===
         if self._periodic_shares:
-            self.share_periodic = self._create_share_vars(self._periodic_shares, 'share|periodic', temporal=False)
-            for eid, var in self.share_periodic.items():
-                summed = var.sum('contributor').expand_dims(effect=[eid])
-                self._eq_periodic.lhs -= summed.reindex({'effect': self.data.effect_index}, fill_value=0)
+            self.share_periodic = self._create_share_var(self._periodic_shares, 'share|periodic', temporal=False)
+            self._eq_periodic.lhs -= self.share_periodic.sum('contributor')
 
         # === Apply periodic constants directly ===
         for const in self._periodic_constant_defs:
@@ -608,30 +604,43 @@ def _share_coords(self, element_dim: str, element_index, temporal: bool = True)
             }
         )
 
-    def _create_share_vars(
+    def _create_share_var(
         self,
         accum: dict[str, list[linopy.LinearExpression]],
         name: str,
         temporal: bool,
-    ) -> dict[str, linopy.Variable]:
-        """Create one share variable per effect from accumulated per-effect expression lists.
+    ) -> linopy.Variable:
+        """Create one share variable with (contributor, effect, ...) dims.
 
         accum structure: {effect_id: [expr1, expr2, ...]} where each expr has
         (contributor, ...other_dims) dims — no effect dim.
 
-        Within a single effect, all expressions share the same non-helper dims,
-        so linopy.merge along 'contributor' is cheap (no cross-dim alignment).
+        Constraints are added per-effect: var.sel(effect=eid) == merged_for_eid,
+        which avoids cross-effect alignment.
 
         Returns:
-            dict mapping effect_id -> linopy.Variable with dims (contributor, time/period).
+            linopy.Variable with dims (contributor, effect, time/period).
         """
         import pandas as pd
 
         if not accum:
-            return {}
+            return None
+
+        # Collect all contributor IDs across all effects
+        all_contributor_ids: set[str] = set()
+        for expr_list in accum.values():
+            for expr in expr_list:
+                all_contributor_ids.update(str(c) for c in expr.data.coords['contributor'].values)
+
+        contributor_index = pd.Index(sorted(all_contributor_ids), name='contributor')
+        effect_index = self.data.effect_index
+        coords = self._share_coords('contributor', contributor_index, temporal=temporal)
+        var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=name)
 
-        result: dict[str, linopy.Variable] = {}
+        # Track which (effect, contributor) combos are constrained
+        constrained: dict[str, set[str]] = {}
 
+        # Add per-effect constraints
         for eid, expr_list in accum.items():
             # Sum expressions that share the same contributors (e.g. rate + status + startup)
             combined: dict[str, linopy.LinearExpression] = {}
@@ -644,22 +653,27 @@ def _create_share_vars(
                     else:
                         combined[cid_str] = single
 
-            contributor_ids = list(combined.keys())
-            contributor_index = pd.Index(contributor_ids, name='contributor')
-            base_dims = None if temporal else ['period', 'scenario']
-            model_coords = self.model.get_coords(base_dims) or {}
-            coords = xr.Coordinates({'contributor': contributor_index, **{k: v for k, v in model_coords.items()}})
-            var_name = f'{name}({eid})'
-            var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=var_name)
-
-            # Merge all per-contributor expressions along 'contributor'
+            # Merge per-contributor expressions along 'contributor'
             merged = linopy.merge(list(combined.values()), dim='contributor')
-            self.model.add_constraints(var == merged, name=var_name)
-
-            result[eid] = var
+            # Constrain: var.sel(effect=eid, contributor=contributors) == merged
+            var_slice = var.sel(effect=eid, contributor=list(combined.keys()))
+            self.model.add_constraints(var_slice == merged, name=f'{name}({eid})')
+            constrained[eid] = set(combined.keys())
+
+        # Constrain uncovered (effect, contributor) combos to 0 to avoid unbounded vars
+        all_effects = [str(e) for e in effect_index]
+        all_contributors = [str(c) for c in contributor_index]
+        for eid in all_effects:
+            covered = constrained.get(eid, set())
+            uncovered = [c for c in all_contributors if c not in covered]
+            if uncovered:
+                self.model.add_constraints(
+                    var.sel(effect=eid, contributor=uncovered) == 0,
+                    name=f'{name}({eid})|zero',
+                )
 
         accum.clear()
-        return result
+        return var
 
     def get_periodic(self, effect_id: str) -> linopy.Variable:
         """Get periodic variable for a specific effect."""
diff --git a/flixopt/statistics_accessor.py b/flixopt/statistics_accessor.py
index c1f05b210..bd454e48e 100644
--- a/flixopt/statistics_accessor.py
+++ b/flixopt/statistics_accessor.py
@@ -805,22 +805,17 @@ def _create_effects_dataset(self, mode: Literal['temporal', 'periodic', 'total']
 
         # Determine modes to process
         modes_to_process = ['temporal', 'periodic'] if mode == 'total' else [mode]
-        # Detect contributors from per-effect share variables (share|temporal|{effect}, share|periodic|{effect})
+        # Detect contributors from combined share variables (share|temporal, share|periodic)
         detected_contributors: set[str] = set()
-        share_var_names: dict[str, list[str]] = {'temporal': [], 'periodic': []}
-        for var_name in solution:
-            for mode_key in ('temporal', 'periodic'):
-                prefix = f'share|{mode_key}('
-                if str(var_name).startswith(prefix):
-                    share_var_names[mode_key].append(str(var_name))
-
         for current_mode in modes_to_process:
-            for share_name in share_var_names[current_mode]:
-                share_da = solution[share_name]
-                for c in share_da.coords['contributor'].values:
-                    base_name = str(c).split('(')[0] if '(' in str(c) else str(c)
-                    if base_name not in effect_labels:
-                        detected_contributors.add(str(c))
+            share_name = f'share|{current_mode}'
+            if share_name not in solution:
+                continue
+            share_da = solution[share_name]
+            for c in share_da.coords['contributor'].values:
+                base_name = str(c).split('(')[0] if '(' in str(c) else str(c)
+                if base_name not in effect_labels:
+                    detected_contributors.add(str(c))
 
         contributors = sorted(detected_contributors)
 
@@ -863,13 +858,15 @@ def get_contributor_type(contributor: str) -> str:
                     conversion_factors[effect] = 1  # Direct contribution
 
                     for source_effect, factor in conversion_factors.items():
-                        share_name = f'share|{current_mode}({source_effect})'
+                        share_name = f'share|{current_mode}'
                         if share_name not in solution:
                             continue
                         share_da = solution[share_name]
+                        if source_effect not in share_da.coords['effect'].values:
+                            continue
                         if contributor not in share_da.coords['contributor'].values:
                             continue
-                        da = share_da.sel(contributor=contributor) * factor
+                        da = share_da.sel(effect=source_effect, contributor=contributor, drop=True) * factor
                         # For total mode, sum temporal over time (apply cluster_weight for proper weighting)
                         if mode == 'total' and current_mode == 'temporal' and 'time' in da.dims:
                             weighted = da * self._fs.weights.get('cluster', 1.0)
diff --git a/tests/test_effect.py b/tests/test_effect.py
index af3ad9dc4..f26ad3438 100644
--- a/tests/test_effect.py
+++ b/tests/test_effect.py
@@ -164,7 +164,7 @@ def test_shares(self, basic_flow_system_linopy_coords, coords_config):
 
         # Check share allocation variables exist (e.g., share|temporal_from_effect for effect-to-effect shares)
         # These are managed by the EffectsModel
-        assert any(v.startswith('share|temporal') for v in model.variables), 'Temporal share variable(s) should exist'
+        assert 'share|temporal' in model.variables, 'Temporal share variable should exist'
 
         # Access individual effect variables using batched model + sel
         _effect2_periodic = model.variables['effect|periodic'].sel(effect='Effect2')
diff --git a/tests/test_flow.py b/tests/test_flow.py
index 2aafb08a1..2feabf39e 100644
--- a/tests/test_flow.py
+++ b/tests/test_flow.py
@@ -103,7 +103,7 @@ def test_effects_per_flow_hour(self, basic_flow_system_linopy_coords, coords_con
 
         # Batched temporal shares are managed by the EffectsModel
         assert any(c.startswith('share|temporal') for c in model.constraints), (
-            'Batched temporal share constraint(s) should exist'
+            'Temporal share constraint(s) should exist'
         )
 
         # Check batched effect variables exist
@@ -357,7 +357,7 @@ def test_flow_invest_divest_effects(self, basic_flow_system_linopy_coords, coord
         assert 'effect|periodic' in model.variables
 
         # Check that temporal share exists for the flow's effects
-        assert any(v.startswith('share|temporal') for v in model.variables)
+        assert 'share|temporal' in model.variables
 
 
 class TestFlowOnModel:
diff --git a/tests/test_functional.py b/tests/test_functional.py
index 341bc9ee8..f309b52de 100644
--- a/tests/test_functional.py
+++ b/tests/test_functional.py
@@ -129,7 +129,7 @@ def test_minimal_model(solver_fixture, time_steps_fixture):
     )
 
     assert_allclose(
-        flow_system.solution['share|temporal(costs)'].sel(contributor='Gastarif(Gas)').values[:-1],
+        flow_system.solution['share|temporal'].sel(effect='costs', contributor='Gastarif(Gas)').values[:-1],
         [-0.0, 20.0, 40.0, -0.0, 20.0],
         rtol=1e-5,
         atol=1e-10,

From e640104f849ec27c08e73968f8c90e45ba68d02d Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 03:22:57 +0100
Subject: [PATCH 11/16] =?UTF-8?q?=E2=8F=BA=20Instead=20of=20adding=20zero-?=
 =?UTF-8?q?constraints=20for=20uncovered=20combos,=20we=20should=20just=20?=
 =?UTF-8?q?set=20lower=3D0,=20upper=3D0=20on=20those=20entries=20(fix=20th?=
 =?UTF-8?q?e=20bounds),=20or=20better=20yet=20=E2=80=94=20use=20a=20mask?=
 =?UTF-8?q?=20on=20the=20per-effect=20=20=20constraints=20and=20set=20the?=
 =?UTF-8?q?=20variable=20bounds=20to=200=20for=20uncovered=20combos.=20The?=
 =?UTF-8?q?=20simplest=20fix:=20create=20the=20variable=20with=20lower=3D0?=
 =?UTF-8?q?,=20upper=3D0=20by=20default,=20then=20only=20the=20covered=20e?=
 =?UTF-8?q?ntries=20need=20=20=20constraints.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 flixopt/effects.py | 36 +++++++++++++++++-------------------
 1 file changed, 17 insertions(+), 19 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index bf5bfa8f5..57bbdd1d5 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -633,14 +633,25 @@ def _create_share_var(
                 all_contributor_ids.update(str(c) for c in expr.data.coords['contributor'].values)
 
         contributor_index = pd.Index(sorted(all_contributor_ids), name='contributor')
-        effect_index = self.data.effect_index
         coords = self._share_coords('contributor', contributor_index, temporal=temporal)
-        var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=name)
 
-        # Track which (effect, contributor) combos are constrained
-        constrained: dict[str, set[str]] = {}
+        # Start with bounds fixed at 0; widen to -inf/+inf for covered (effect, contributor) combos
+        lower = xr.DataArray(0.0, coords=coords)
+        upper = xr.DataArray(0.0, coords=coords)
 
-        # Add per-effect constraints
+        # Pre-scan to find covered combos and widen bounds
+        covered_map: dict[str, list[str]] = {}
+        for eid, expr_list in accum.items():
+            cids = set()
+            for expr in expr_list:
+                cids.update(str(c) for c in expr.data.coords['contributor'].values)
+            covered_map[eid] = sorted(cids)
+            lower.loc[dict(effect=eid, contributor=covered_map[eid])] = -np.inf
+            upper.loc[dict(effect=eid, contributor=covered_map[eid])] = np.inf
+
+        var = self.model.add_variables(lower=lower, upper=upper, coords=coords, name=name)
+
+        # Add per-effect constraints (only for covered combos)
         for eid, expr_list in accum.items():
             # Sum expressions that share the same contributors (e.g. rate + status + startup)
             combined: dict[str, linopy.LinearExpression] = {}
@@ -656,21 +667,8 @@ def _create_share_var(
             # Merge per-contributor expressions along 'contributor'
             merged = linopy.merge(list(combined.values()), dim='contributor')
             # Constrain: var.sel(effect=eid, contributor=contributors) == merged
-            var_slice = var.sel(effect=eid, contributor=list(combined.keys()))
+            var_slice = var.sel(effect=eid, contributor=covered_map[eid])
             self.model.add_constraints(var_slice == merged, name=f'{name}({eid})')
-            constrained[eid] = set(combined.keys())
-
-        # Constrain uncovered (effect, contributor) combos to 0 to avoid unbounded vars
-        all_effects = [str(e) for e in effect_index]
-        all_contributors = [str(c) for c in contributor_index]
-        for eid in all_effects:
-            covered = constrained.get(eid, set())
-            uncovered = [c for c in all_contributors if c not in covered]
-            if uncovered:
-                self.model.add_constraints(
-                    var.sel(effect=eid, contributor=uncovered) == 0,
-                    name=f'{name}({eid})|zero',
-                )
 
         accum.clear()
         return var

From 6e933af0820b89c023db5c0666970bf372dfd988 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 03:26:39 +0100
Subject: [PATCH 12/16] Only create variables needed

---
 flixopt/effects.py | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index 57bbdd1d5..67a2850d4 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -633,23 +633,24 @@ def _create_share_var(
                 all_contributor_ids.update(str(c) for c in expr.data.coords['contributor'].values)
 
         contributor_index = pd.Index(sorted(all_contributor_ids), name='contributor')
+        effect_index = self.data.effect_index
         coords = self._share_coords('contributor', contributor_index, temporal=temporal)
 
-        # Start with bounds fixed at 0; widen to -inf/+inf for covered (effect, contributor) combos
-        lower = xr.DataArray(0.0, coords=coords)
-        upper = xr.DataArray(0.0, coords=coords)
-
-        # Pre-scan to find covered combos and widen bounds
+        # Build mask: only create variables for (effect, contributor) combos that have expressions
+        mask = xr.DataArray(
+            np.zeros((len(contributor_index), len(effect_index)), dtype=bool),
+            dims=['contributor', 'effect'],
+            coords={'contributor': contributor_index, 'effect': effect_index},
+        )
         covered_map: dict[str, list[str]] = {}
         for eid, expr_list in accum.items():
             cids = set()
             for expr in expr_list:
                 cids.update(str(c) for c in expr.data.coords['contributor'].values)
             covered_map[eid] = sorted(cids)
-            lower.loc[dict(effect=eid, contributor=covered_map[eid])] = -np.inf
-            upper.loc[dict(effect=eid, contributor=covered_map[eid])] = np.inf
+            mask.loc[dict(effect=eid, contributor=covered_map[eid])] = True
 
-        var = self.model.add_variables(lower=lower, upper=upper, coords=coords, name=name)
+        var = self.model.add_variables(lower=-np.inf, upper=np.inf, coords=coords, name=name, mask=mask)
 
         # Add per-effect constraints (only for covered combos)
         for eid, expr_list in accum.items():

From 327de3ed9c02382623813629c05073095f6c4db3 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 03:39:00 +0100
Subject: [PATCH 13/16] =?UTF-8?q?=5Fcreate=5Fshare=5Fvar=20went=20from=201?=
 =?UTF-8?q?,674ms=20=E2=86=92=20116ms=20=E2=80=94=20a=2014x=20speedup!=20T?=
 =?UTF-8?q?he=20reindex=20+=20+=20approach=20is=20much=20faster=20than=20p?=
 =?UTF-8?q?er-contributor=20sel=20+=20merge?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 flixopt/effects.py | 25 ++++++++++---------------
 1 file changed, 10 insertions(+), 15 deletions(-)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index 67a2850d4..6053d2220 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -654,21 +654,16 @@ def _create_share_var(
 
         # Add per-effect constraints (only for covered combos)
         for eid, expr_list in accum.items():
-            # Sum expressions that share the same contributors (e.g. rate + status + startup)
-            combined: dict[str, linopy.LinearExpression] = {}
-            for expr in expr_list:
-                for cid in expr.data.coords['contributor'].values:
-                    cid_str = str(cid)
-                    single = expr.sel(contributor=[cid])
-                    if cid_str in combined:
-                        combined[cid_str] = linopy.merge([combined[cid_str], single])
-                    else:
-                        combined[cid_str] = single
-
-            # Merge per-contributor expressions along 'contributor'
-            merged = linopy.merge(list(combined.values()), dim='contributor')
-            # Constrain: var.sel(effect=eid, contributor=contributors) == merged
-            var_slice = var.sel(effect=eid, contributor=covered_map[eid])
+            contributors = covered_map[eid]
+            if len(expr_list) == 1:
+                merged = expr_list[0].reindex(contributor=contributors)
+            else:
+                # Reindex all to common contributor set, then sum via linopy.merge (_term addition)
+                aligned = [e.reindex(contributor=contributors) for e in expr_list]
+                merged = aligned[0]
+                for a in aligned[1:]:
+                    merged = merged + a
+            var_slice = var.sel(effect=eid, contributor=contributors)
             self.model.add_constraints(var_slice == merged, name=f'{name}({eid})')
 
         accum.clear()

From 47711f843169c86d4de9793d836fd6fa661cf980 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 15:15:38 +0100
Subject: [PATCH 14/16] Revert

---
 flixopt/results.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/flixopt/results.py b/flixopt/results.py
index 74d1e71ae..3d95357eb 100644
--- a/flixopt/results.py
+++ b/flixopt/results.py
@@ -794,7 +794,7 @@ def get_effect_shares(
 
         ds = xr.Dataset()
 
-        share_var_name = f'share|{mode}({effect})'
+        share_var_name = f'share|{mode}'
         if share_var_name in self.solution:
             share_var = self.solution[share_var_name]
             contributor_dim = None
@@ -803,9 +803,10 @@ def get_effect_shares(
                     contributor_dim = dim
                     break
             if contributor_dim is not None and element in share_var.coords[contributor_dim].values:
-                selected = share_var.sel({contributor_dim: element}, drop=True)
-                label = f'{element}->{effect}({mode})'
-                ds = xr.Dataset({label: selected})
+                if effect in share_var.coords['effect'].values:
+                    selected = share_var.sel({contributor_dim: element, 'effect': effect}, drop=True)
+                    label = f'{element}->{effect}({mode})'
+                    ds = xr.Dataset({label: selected})
 
         if include_flows:
             if element not in self.components:

From b30f8c38f548923716494ab69409f3fd48fa1565 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 15:20:13 +0100
Subject: [PATCH 15/16] Revert

---
 flixopt/transform_accessor.py | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 93c99491e..1123b9f66 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -469,11 +469,8 @@ def _cluster_indices_per_timestep(self) -> xr.DataArray:
 
     @staticmethod
     def _get_mode(var_name: str) -> ExpansionMode:
-        """Look up expansion mode for a variable name via suffix matching."""
-        for suffix, mode in NAME_TO_EXPANSION.items():
-            if var_name == suffix or var_name.endswith(suffix) or var_name.startswith(suffix):
-                return mode
-        return ExpansionMode.REPEAT
+        """Look up expansion mode for a variable name."""
+        return NAME_TO_EXPANSION.get(var_name, ExpansionMode.REPEAT)
 
     def _append_final_state(self, expanded: xr.DataArray, da: xr.DataArray) -> xr.DataArray:
         """Append final state value from original data to expanded data."""

From 35af3d2388c3c3d2a61c8cb1c8958fbfcd38c248 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 1 Feb 2026 15:23:28 +0100
Subject: [PATCH 16/16]   1. effects.py: add_temporal_contribution and
 add_periodic_contribution now raise ValueError if a DataArray has no effect
 dimension and no effect= argument is provided.   2. statistics_accessor.py:
 Early return with empty xr.Dataset() when no contributors are detected,
 preventing xr.concat from failing on an empty list.

---
 flixopt/effects.py             | 8 ++++++++
 flixopt/statistics_accessor.py | 3 +++
 2 files changed, 11 insertions(+)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index 6053d2220..87ed65776 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -379,6 +379,10 @@ def add_temporal_contribution(
         if isinstance(defining_expr, xr.DataArray):
             if effect is not None:
                 defining_expr = defining_expr.expand_dims(effect=[effect])
+            elif 'effect' not in defining_expr.dims:
+                raise ValueError(
+                    "DataArray contribution must have an 'effect' dimension or an explicit effect= argument."
+                )
             self._temporal_constant_defs.append(defining_expr)
         else:
             self._accumulate_shares(self._temporal_shares, self._as_expression(defining_expr), effect)
@@ -401,6 +405,10 @@ def add_periodic_contribution(
         if isinstance(defining_expr, xr.DataArray):
             if effect is not None:
                 defining_expr = defining_expr.expand_dims(effect=[effect])
+            elif 'effect' not in defining_expr.dims:
+                raise ValueError(
+                    "DataArray contribution must have an 'effect' dimension or an explicit effect= argument."
+                )
             self._periodic_constant_defs.append(defining_expr)
         else:
             self._accumulate_shares(self._periodic_shares, self._as_expression(defining_expr), effect)
diff --git a/flixopt/statistics_accessor.py b/flixopt/statistics_accessor.py
index bd454e48e..9b088bad8 100644
--- a/flixopt/statistics_accessor.py
+++ b/flixopt/statistics_accessor.py
@@ -819,6 +819,9 @@ def _create_effects_dataset(self, mode: Literal['temporal', 'periodic', 'total']
 
         contributors = sorted(detected_contributors)
 
+        if not contributors:
+            return xr.Dataset()
+
         # Build metadata for each contributor
         def get_parent_component(contributor: str) -> str:
             if contributor in self._fs.flows: