Model

The ehubX model is built from sub-models (or modules) that address different aspects of the energy system. Each model adds to (and in some occasions removes from) a Pyomo model object. The models are loosely organized quite similarily to the input file structure:

energy_system_model
|-- stage_model
|-- hub_model
|-- ec_model
|-- times_model
|-- import_model
|-- export_model
|-- demand_model
    |-- load_shedding_model
    |-- load_shifting_model
|-- network_model
|-- tech_model
    |-- stor_tech_model
    |-- conv_tech_model
        |-- solar_tech_model
    |-- ebm_tech_model
    |-- ates_model
    |-- hp_tech_model
|-- self_sufficiency_model
|-- autonomy_model

Each model is responsible for a specific aspect of the energy system and most models are built independently from each other, relying only on models which lie higher up in the hierarchy. This modular structure allows for easy extension of the framework by adding new models or modifying existing ones. The energy system model serves as the overarching bridge that integrates all sub-models into a cohesive whole. We will start here with the individual model although some users may find it more intuitive to start with the energy system model first.

Stage model

This module introduces an index set \(\mathcal{S}_{Stage}\) of stages. They consist of unique string labels (e.g.; \(S1\)).

Hub model

This module introduces an index set \(\mathcal{S}_{Hub}\) of hubs. They consist of unique string labels (e.g.; \(H1\)).

Ec model

This module introduces an index set \(\mathcal{S}_{Ec}\) of ecs (energy carriers). They consist of unique string labels (e.g.; \(E1\)).

Times model

This module introduces two index sets \(\mathcal{S}_{TimeHorizon}\) and \(\mathcal{S}_{Time}\) (energy carriers). They each consist of unique integers \(1, 2, ...\) that act as time indices. In case where reduced-order modeling is employed, \(\mathcal{S}_{Time}\) holds the cluster indices and \(\mathcal{S}_{TimeHorizon}\) holds the full-horizon indices. If no clustering is used, these two sets are identical.

Import model

Import of energy carriers is possible for certain tuples of stages, hubs and ecs, as defined piecewise in imports.yaml. The set of all importable tuples is labeled \(\mathcal{S}_{ImpTuple} \subset \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}\). The variable

\(\mathcal{V}_{Imp}: \mathcal{S}_{ImpTuple} \times \mathcal{S}_{Time} \to \mathbb{R}^+_0\)

denotes the import amounts per time step. These values may be constrained

1. in a piecewise manner by \(min[s, h, e, t] \le \mathcal{V}_{Imp}[s, h, e, t] \le max[s, h, e, t]\) where \(min\) and \(max\) are parameters from imports.yaml.

2. in a summed-up manner by \(sum\_min[s, h, e] \le \sum_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{Imp}[s, h, e, t] \le sum\_max[s, h, e]\) where \(sum\_min\) and \(sum\_max\) are parameters from imports.yaml and \(weight\) is a clustering value.

Certain costs are associated with imports which are defined per import tuple using the parameter \(price\) from imports.yaml by a fixed variable

\begin{align*} \mathcal{V}_{ImpCost}: \mathcal{S}_{ImpTuple} &\to \mathbb{R} \\ \mathcal{V}_{ImpCost}[s, h, e] &= \sum_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot price[s, h, e, t] \cdot \mathcal{V}_{Imp}[s, h, e, t] \end{align*}

Similarily, imports might come with CO2 emissions. These are quantified using the parameter \(co2\) from imports.yaml and tracked in a fixed variable

\begin{align*} \mathcal{V}_{ImpCo2}: \mathcal{S}_{ImpTuple} &\to \mathbb{R} \\ \mathcal{V}_{ImpCo2}[s, h, e] &= \sum_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot co2[s, h, e, t] \cdot \mathcal{V}_{Imp}[s, h, e, t] \end{align*}

Information about costs and CO2 for imports are required in the energy system model. They are bundled in the following fixed variables for convenience:

\begin{align*} \mathcal{V}_{ImpCostTotal} &\in \mathbb{R}^+_0, \quad &\mathcal{V}_{ImpCostTotal} &= \sum_{(s, h, e) \in \mathcal{S}_{ImpTuple}} \mathcal{V}_{ImpCost}[s, h, e] \\ \mathcal{V}_{ImpCo2Total}: \mathcal{S}_{Stage} &\to \mathbb{R}, \quad &\mathcal{V}_{ImpCo2Total}[s] &= \sum_{\substack{(s', h, e) \in \mathcal{S}_{ImpTuple} \\ s'=s}} \mathcal{V}_{ImpCo2}[s', h, e] \end{align*}

Export model

Export of energy carriers is possible for certain tuples of stages, hubs and ecs, as defined piecewise in exports.yaml. The set of all exportable tuples is labeled \(\mathcal{S}_{ExpTuple} \subset \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}\). The variable

\(\mathcal{V}_{Exp}: \mathcal{S}_{ExpTuple} \times \mathcal{S}_{Time} \to \mathbb{R}^+_0\)

denotes the export amounts per time step. These values may be constrained

1. in a piecewise manner by \(min[s, h, e, t] \le \mathcal{V}_{Exp}[s, h, e, t] \le max[s, h, e, t]\) where \(min\) and \(max\) are parameters from exports.yaml.

2. in a summed-up manner by \(sum\_min[s, h, e] \le \sum_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{Exp}[s, h, e, t] \le sum\_max[s, h, e]\) where \(sum\_min\) and \(sum\_max\) are parameters from exports.yaml and \(weight\) is a clustering value.

Certain profits are associated with exports which are defined per export tuple using the parameter \(price\) from exports.yaml by a fixed variable

\begin{align*} \mathcal{V}_{ExpProfit}: \mathcal{S}_{ExpTuple} &\to \mathbb{R} \\ \mathcal{V}_{ExpProfit}[s, h, e] &= \sum_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot price[s, h, e, t] \cdot \mathcal{V}_{Exp}[s, h, e, t] \end{align*}

Similarily, exports might come with reductions in CO2 emissions. These are quantified using the parameter \(co2\) from exports.yaml and tracked in a fixed variable

\begin{align*} \mathcal{V}_{ExpCo2}: \mathcal{S}_{ExpTuple} &\to \mathbb{R} \\ \mathcal{V}_{ExpCo2}[s, h, e] &= \sum_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot co2[s, h, e, t] \cdot \mathcal{V}_{Exp}[s, h, e, t] \end{align*}

Information about profits and CO2 for exports are required in the energy system model. They are bundled in the following fixed variables for convenience:

\begin{align*} \mathcal{V}_{ExpProfitTotal} &\in \mathbb{R}^+_0, \quad &\mathcal{V}_{ExpProfitTotal} &= \sum_{(s, h, e) \in \mathcal{S}_{ExpTuple}} \mathcal{V}_{ExpProfit}[s, h, e] \\ \mathcal{V}_{ExpCo2Total}: \mathcal{S}_{Stage} &\to \mathbb{R}, \quad &\mathcal{V}_{ExpCo2Total}[s] &= \sum_{\substack{(s', h, e) \in \mathcal{S}_{ExpTuple} \\ s'=s}} \mathcal{V}_{ExpCo2}[s', h, e] \end{align*}

Demand model

There are demand values for certain tuples of stages, hubs and ecs, as defined piecewise in demands.yaml. As indicated in this file, there are two types of demands: Demand-profiles (given per time step), and demand-sums (indicating a demand to the total supply over the entire time horizon). Since each of these applies to a tuple of stage, hub, and ec, we define accordingly:

\begin{align*} \mathcal{S}_{DemandProfileTuple}, ~\mathcal{S}_{DemandSumTuple} \subset &\mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}, \quad \mathcal{S}_{DemandProfileTuple} \cap \mathcal{S}_{DemandSumTuple} = \emptyset \end{align*}

We also label the joint set of all demand tuples as

$$ \mathcal{S}_{DemandTuple} = \mathcal{S}_{DemandProfileTuple} \cup \mathcal{S}_{DemandSumTuple} $$

For each of these demand tuples, the system has to provide a supply that meets the demand values.

$$ \mathcal{V}_{DemandSupply}: \mathcal{S}_{DemandTuple} \times \mathcal{S}_{Time} \to \mathbb{R}^+_0 $$

Now these variables need to be properly constrained according to the type of demand. For demand-profiles, this is handled by the overarching energy system model since the supply is not only influenced by the demand values itself, but also by the load shifting and load shedding impacts. For the demand-sums, load shedding and load shifting are disallowed, so we can define the constraint here:

$$ \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{DemandSupply}[s, h, e, t] = demand\_sum[s, h, e]. $$

Additionally, certain formulations in other modules require a so-called big-M parameter which is very common in MILP formulations. Essentially, this is a positive value that is “as large as it needs to be” for the situation at hand. As an example, if a binary variable \(y \in \{0,1\}\) needs to be true for an associated continuous variable \(x \in \mathbb{R}\) to be larger than zero, we would use a constraint of the form \(x \le M \cdot y\) with a value \(M > 0\). However, if \(y = 1\) is the case, the variable \(x\) should still be allowed to reach any value in its permitted range. If \(x\) has an upper limit \(\overline{x}\) from the model formulation already, the case is easy: Just use \(M = \overline{x}\). However, this is often not the case from our model formulations (as an example, a technology’s capacity variable might not have a cap_max setting). For these cases, we require a generic big-M parameter that is larger than any value our variables are likely to assume. This is obviously not an exact science since we are dealing with unrestricted variables but we can use the demand module to make an educated guess:

• For a demand-profile tuple \((s, h, e) \in \mathcal{S}_{DemandProfileTuple}\) and a time \(t \in \mathcal{S}_{Time}\), the value \(demand\_profile[s, h, e, t]\) can serve as an indicator for the magnitude of the system. We can therefore derive a big-M parameter $$ BigM_{DemandProfile}[s, h, e, t] = C_{BigM} \cdot demand\_profile[s, h, e, t] $$ which we postulate to be larger than an an inflow or outflow to an energy balance node (cf. the energy balance equation in the energy system model) from any module. Here, \(C_{BigM} \approx \mathcal{O}(10)\) is reasonable multiplier. The factor of 10 is somewhat arbitrary but should be large enough to cover most situations.

• Similar to the point above, we sometimes require a big-M parameter for the sum of inflows/outflows to energy balance nodes. For this, we use $$ BigM_{DemandProfileSum}[s, h, e] = C_{BigMSum} \cdot \max\limits_{t \in \mathcal{S}_{Time}} demand\_profile[s, h, e, t] $$ with a similar logic than above. Since these values tend to be much larger than the per-timestep values above, we tend to use a smaller multiplier of \(C_{BigMSum} \approx \mathcal{O}(1)\) to keep the system relatively well-scaled.

Note the underlying assumption in this approach that the demand series will be able to quantify the general dimensionality of the system. Especially in systems with small demand profiles (but also in general), it might be prudent to explicitly set values for certain parameters so that the generic big-M parameter does not have to be used. Tailored warning messages are routinely logged for this reason, notifying the user which parameters are missing to calculate a specific big-M parameter (the most common candidate for this is \(cap\_max\) from hubs.yaml and network_links.yaml). Additionally, if no demand series exists for a specific carrier, a fallback option is required. For our case, this comes in the form of the parameter heur_max from ecs.yaml which is therefore mandatory for ecs without a demand series.

Furthermore, in certain formulations demand satisfaction may be relaxed through the introduction of an unmet demand variable. For this purpose, we define

\begin{align*} \mathcal{V}_{DemandUnmet}[s,h,e,t] \in \mathbb{R}_0^+ \quad \forall (s,h,e) \in \mathcal{S}_{DemandTuple},~ t \in \mathcal{S}_{Time} \end{align*}

This variable represents the portion of demand that is not satisfied by the system.

Its use is not part of the core demand formulation, but is instead activated by higher-level modules such as the autonomy model.

If unmet demand is enabled, the demand-sum constraint is extended such that both served and unmet demand contribute to the total:

\begin{align*} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s,t] \cdot \left( \mathcal{V}_{DemandSupply}[s,h,e,t] + \mathcal{V}_{DemandUnmet}[s,h,e,t] \right) = demand\_sum[s,h,e] \quad \forall (s,h,e) \in \mathcal{S}_{DemandSumTuple} \end{align*}

Hence, the prescribed total demand remains fully accounted for, either through actual supply or through unmet demand.

Load shedding model

The load shedding module parses the data model for the set of stage-hub-ec tuples \(\mathcal{S}_{LoadSheddingTuple} \subset \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}\) for which load shedding is enabled (as specified in demands.yaml). Naturally, these are only allowed to be tuples that are also associated with a demand profile, and a data validation procedure ensures this before the model is built. A load shedding variable

$$ \mathcal{V}_{LoadShedding}: \mathcal{S}_{LoadShedding} \times \mathcal{S}_{Time} \to \mathbb{R}^+_0 $$

measures the amount of demand that is chosen not be delivered by the system. This variable is constrained from above in the following way:

\begin{align*} \mathcal{V}_{LoadShedding}[s, h, e, t] \le \min \big( &max\_abs[s, h, e, t], \\ &max\_rel[s, h, e, t] \cdot demand[s, h, e, t] \big) \end{align*}

where \(max\_abs\), \(max\_rel\) and \(demand\) are parameters from demands.yaml. Due to the integration of \(\mathcal{V}_{LoadShedding}\) into the energy system model, the variable is not able to take values larger than \(\mathcal{V}_{Demand}\), thereby ensuring that the model is not able to shed more amounts than the actual demand. Additionally, a warning is logged if \(max\_rel\) is set to a value larger than one to make the user aware of this circumstance.

Shedding loads is associated with certain costs. These are measure by the fixed variable

\begin{align*} \mathcal{V}_{LoadSheddingCost}: \mathcal{S}_{LoadShedding} &\to \mathbb{R}^+_0, \\ \mathcal{V}_{LoadSheddingCost}[s, h, e] &= \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot energy\_cost[s, h, e, t] \cdot \mathcal{V}_{LoadShedding}[s, h, e, t] \end{align*}

Here, \(energy\_cost\) is a parameter from demands.yaml and \(weight\) is a clustering value. Inforomation about these costs are required in the energy system model. They are bundled in the following fixed variable for convenience:

$$ \mathcal{V}_{LoadSheddingCostTotal} \in \mathbb R^+_0, \quad \mathcal{V}_{LoadSheddingCostTotal} = \sum\limits_{(s, h, e) \in \mathcal{S}_{LoadSheddingTuple}} \mathcal{V}_{LoadSheddingCost}[s, h, e] $$

Load shifting model

The load shifting module introduces a set of stage-hub-ec index tuples \(\mathcal{S}_{LoadShiftingTuple} \subset \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}\) for which load shifting is possible. The set is gathered from the load shifting entries defined in demands.yaml.

The variable

$$ \mathcal{V}_{LoadShifting}: \mathcal{S}_{LoadShiftingTuple} \times \mathcal{S}_{Time} \to \mathbb{R} $$

represents the amount of load shifting being performed at each time step, with positive values indicating an over-production of the ec and negative values indicating an under-production. We also use the terminology above/below the demand curve in connection to positive/negative load shifting values. The essential characteristic of any load shifting model is that no energy can be created or destroyed by it, only shifts along the time axis can occur. In ehubX, the way this characteristic is formulated is through load shifting intervals: The time horizon \(\mathcal{S}_{TimeHorizon}\) is segmented into equidistant intervals

\begin{align*} \mathcal{S}_{TimeHorizon} &= \Big \{ \underbrace{1, ... , interval\_length}_{\mathcal{I}_1}, ~\underbrace{interval\_length + 1, ..., 2 \cdot interval\_length}_{\mathcal{I}_2}, ... \Big \} \\ &= \mathcal{I}_1 ~\cup~ \mathcal{I}_2 ~\cup~ ... ~\cup~ \mathcal{I}_N ~\cup~ \mathcal{R} \end{align*}

Here, \(interval\_length\) is a parameter from demands.yaml and \(\mathcal{R}\) holds the remaining time indices at the end of the horizon that did not fit into an interval anymore. We now demand neutrality on each load shifting interval (\(n=1,...,N\)) through the constraints

\begin{align*} \sum\limits_{t \in \mathcal{I}_n} \mathcal{V}_{LoadShifting}[s, h, e, t] = 0 \end{align*}

This constraint ensures that no energy can be created or destroyed within any load shifting interval. To also account for the remainder set \(\mathcal{R}\), we additionally introduce a global neutrality condition

$$ \sum\limits_{t \in \mathcal{S}_{TimeHorizon}} \mathcal{V}_{LoadShifting}[s, h, e, t] = 0 $$

The fact that both over-deliverance and under-deliverance of the demand curve is measured in the variable \(\mathcal{V}_{LoadShifting}\) motivates a splitting of this variable into a positive and negative part. This is done by introducing two nonnegative variables

$$ \mathcal{V}_{LoadShiftingAbove}, \mathcal{V}_{LoadShiftingBelow}: \mathcal{S}_{LoadShiftingTuple} \times \mathcal{S}_{Time} \to \mathbb{R}^+_0 $$

These are connected to the main load shifting variable via the constraint

$$ \mathcal{V}_{LoadShifting}[s, h, e, t] = \mathcal{V}_{LoadShiftingAbove}[s, h, e, t] - \mathcal{V}_{LoadShiftingBelow}[s, h, e, t] $$

In a precise formulation, we would usually be forced to introduce an additional binary variable that tracks whether \(\mathcal{V}_{LoadShifting}\) is positive or negative, allowing only one of the variables \(\mathcal{V}_{LoadShiftingAbove}\) and \(\mathcal{V}_{LoadShiftingBelow}\) to take nonzero values at a time. However, this would drastically increase the complexity of the model by the presence of a time-dependent binary variable. Instead, we rely on the fact that if certain costs (i.e.; the energy costs introduced below) are associated with the the above and below variables, an optimal solution will always try to avoid situations where both variables take nonzero values simulteneously. For example, take the case of \(\mathcal{V}_{LoadShifting}[s, h, e, t] = 3\) which could be comprised of the following two value pairs for the sub-variables:

1. \(\mathcal{V}_{LoadShiftingAbove}[s, h, e, t] = 5\) and \(\mathcal{V}_{LoadShiftingBelow}[s, h, e, t] = 2\)

2. \(\mathcal{V}_{LoadShiftingAbove}[s, h, e, t] = 3\) and \(\mathcal{V}_{LoadShiftingBelow}[s, h, e, t] = 0\)

Without the presence of costs for the sub-variables, these two solutions would be identical for an optimization model. However, if \(\mathcal{V}_{LoadShiftingAbove}[s, h, e, t]\) and \(\mathcal{V}_{LoadShiftingBelow}[s, h, e, t]\) both enter in the problem’s object function which is to be minimized, the optimizer will choose the second option to keep the sub-variable’s values as small as possible. A similar reasoning as the one employed here is used in the storage model for the variables \(\mathcal{V}_{StorTechInflow}\) and \(\mathcal{V}_{StorTechOutflow}\).

Next, let us introduce the fixed variables which quantify the energy costs we mentioned above:

\begin{align*} &\mathcal{V}_{LoadShiftingCostEnergy}: \mathcal{S}_{LoadShiftingTuple} \to \mathbb{R}, \\ &\mathcal{V}_{LoadShiftingCostEnergy}[s, h, e] = \sum\limits_{t \in \mathcal{S}_{LoadShiftingTuple}} \bigg( \\ energy\_cost\_above[s, h, e, t] \cdot \mathcal{V}_{LoadShiftingAbove}[s, h, e, t] \\ &+ energy\_cost\_below[s, h, e, t] \cdot \mathcal{V}_{LoadShiftingBelow}[s, h, e, t] \bigg) \end{align*}

Here, \(energy\_cost\_above\) and \(energy\_cost\_below\) are parameters from demands.yaml. Because of the reasons explained above, a warning is written to the logfile if the energy cost parameters are chosen as zero.

Moving on to the next constraint, the maximum load shifting amounts above and below the demand curve may be bounded as follows:

\begin{align*} \mathcal{V}_{LoadShiftingAbove}[s, h, e, t] \le \min \big( &max\_above\_abs[s, h, e, t], \\ &max\_above\_rel[s, h, e, t] \cdot demand[s, h, e, t] \big) \\ \mathcal{V}_{LoadShiftingBelow}[s, h, e, t] \le \min \big( &max\_below\_abs[s, h, e, t], \\ &max\_below\_rel[s, h, e, t] \cdot demand[s, h, e, t] \big) \end{align*}

The parameters \(max\_above\_abs\), \(max\_above\_rel\), \(max\_below\_abs\) and \(max\_below\_rel\) are input parameters of demands.yaml. A warning is included if \(max\_below\_rel\) is larger than one since we cannot withhold more from the demand side than the actual demand itself.

Another limiting factor that is included in the load shifting module is a concept of available capacity, which determines the maximal amount of load shifting which may occur on a single interval. The optimizer is free to scale this capacity itself by using the variables

\begin{align*} \mathcal{V}_{LoadShiftingCap}, \mathcal{V}_{LoadShiftingCapInstl}: \mathcal{S}_{LoadShiftingTuple} \to \mathbb{R}_0^+ \end{align*}

Here, \(\mathcal{V}_{LoadShiftingCap}\) is the total capacity available for load shifting on a single interval, while \(\mathcal{V}_{LoadShiftingCapInstl}\) is amount of additionally installed capacity by the optimizer. Together with the parameter cap_init from demands.yaml, these two variables are connected by the constraint:

$$ \mathcal{V}_{LoadShiftingCap}[s, h, e] = cap\_init[s, h, e] + \mathcal{V}_{LoadShiftingCapInstl}[s, h, e] $$

The capacity variable limits the amount of load shifting which may occur on each interval:

$$ \sum\limits_{t \in \mathcal{I}_n} \mathcal{V}_{LoadShiftingAbove}[s, h, e, t] \le \mathcal{V}_{LoadShiftingCap}[s, h, e] $$

Note that because of the interval neutrality condition, the above constraint is equivalent to

$$ \sum\limits_{t \in \mathcal{I}_n} \mathcal{V}_{LoadShiftingBelow}[s, h, e, t] \le \mathcal{V}_{LoadShiftingCap}[s, h, e] $$

The installation of capacity is associated with certain costs parametrized by the capex_per_cap from demands.yaml. These costs are measured by the fixed variable

$$ \mathcal{V}_{LoadShiftingCostCapex}: \mathcal{S}_{LoadShiftingTuple} \to \mathbb{R}, \qquad \mathcal{V}_{LoadShiftingCostCapex}[s, h, e] = capex\_per\_cap[s, h, e] \cdot \mathcal{V}_{LoadShiftingCapInstl}[s, h, e] $$

Another kind of cost that can be attributed to the load shifting process is the peak amount of load shifting that ever occurs on the time horizon. In order to quantify this, we first need variables measuring the peak value in each direction:

$$ \mathcal{V}_{LoadShiftingAbovePeak}, \mathcal{V}_{LoadShidftingBelowPeak}: \mathcal{S}_{LoadShiftingTuple} \to \mathbb{R}_0^+ $$

If we want to maintain linearity of the model, we cannot simply set these variables to the maximum of \(\mathcal{V}_{LoadShiftingAbove}\) and \(\mathcal{V}_{LoadShiftingBelow}\). Instead, we repeat the solution we already applied for the energy costs above. First, we set a bound in one direction as follows:

\begin{align*} \mathcal{V}_{LoadShiftingAbove}[s, h, e, t] &\le \mathcal{V}_{LoadShiftingAbovePeak}[s, h, e] \\ \mathcal{V}_{LoadShiftingBelow}[s, h, e, t] &\le \mathcal{V}_{LoadShiftingBelowPeak}[s, h, e] \end{align*}

This will require the peak variables to be larger or equal to the largest shift-variables. Next, if the cost associated with the peak variables is nonzero, the optimizer will try to get them as small as possible, thereby forcing them to exactly the maximum value. We measure this cost by the fixed variable

\begin{align*} \mathcal{V}_{LoadShiftingCostPeak}: \mathcal{S}_{LoadShiftingTuple} \to &~\mathbb{R}, \\ \mathcal{V}_{LoadShiftingCostPeak}[s, h, e, t] = (&peak\_cost\_above[s, h, e] \cdot \mathcal{V}_{LoadShiftingAbovePeak}[s, h, e] \\ &peak\_cost\_below[s, h, e] \cdot \mathcal{V}_{LoadShiftingBelowPeak}[s, h, e]) \end{align*}

Here, \(peak\_cost\_above\) and \(peak\_cost\_below\) are parameters taken from demands.yaml.

The third and last kind of cost that can be associated with load shifting is one that occurs every time any amount of load shifting occurs at all. For this to work, the model requires a binary variable that tracks when this is case:

$$ \mathcal{V}_{YLoadShifting}: \mathcal{S}_{LoadShiftingTuple} \times \mathcal{S}_{Time} \to \{0, 1\} $$

It has to be mentioned that adding this variable to the model might make it exceedingly harder to solve because of the multitude of binary decision variables that are involved. For this reason, ehubX only adds this variable when we need it, i.e.; if the fix cost parameter \(fix\_cost\) from demands.yaml is given as a time series or if its default value is nonzero. In order to force this variable to 1 every time load shifting occurs, we require a big-M parameter \(BigM_{LoadShiftTotal}[s, h, e, t]\) and the constraint

$$ \mathcal{V}_{LoadShiftingAbove}[s, h, e, t] + \mathcal{V}_{LoadShiftingBelow}[s, h, e, t] \le BigM_{LoadShiftTotal}[s, h, e, t] \cdot \mathcal{V}_{YLoadShifting}[s, h, e, t] $$

The big-M parameter needs to be larger than any value the left-hand side could possibly take. For our purposes, we choose the following value:

\begin{align*} BigM_{LoadShiftTotal}[s, h, e, t] = \min \Big( &max\_above\_abs[s, h, e, t], max\_above\_rel[s, h, e, t] \cdot demand[s, h, e, t], \\ &10 \cdot demand[s, h, e, t] \Big) \\ + \min \Big( &max\_below\_abs[s, h, e, t], max\_below\_rel[s, h, e, t] \cdot demand[s, h, e, t], \\ &demand[s, h, e, t] \Big) \end{align*}

With the above constraint, \(\mathcal{V}_{YLoadShifting}\) will always be 1 if \(\mathcal{V}_{LoadShifting}\) is nonzero. Therefore, we can use it to specify the fixed cost variable

\begin{align*} &\mathcal{V}_{LoadShiftingCostFix}: \mathcal{S}_{LoadShifting} \to \mathbb{R}^+_0 \\ &\mathcal{V}_{LoadShiftingCostFix}[s, h, e] = \sum\limits_{t \in \mathcal{S}_{Time}} \cdot fix\_cost[s, h, e, t] \cdot \mathcal{V}_{YLoadShifting}[s, h, e, t] \end{align*}

Since all costs must enter the energy system model, we create a single variable that holds the entire load shifting costs of the model:

\begin{align*} \mathcal{V}_{LoadShiftingCostTotal} \in \mathbb{R}, \quad &\mathcal{V}_{LoadShiftingCostTotal} = \\ \sum\limits_{(s, h, e) \in \mathcal{S}_{LoadShiftingTuple}} \Big( &\mathcal{V}_{LoadShiftingCostCapex}[s, h, e] + \mathcal{V}_{LoadShiftingCostEnergy}[s, h, e] \\ &+ \mathcal{V}_{LoadShiftingCostPeak}[s, h, e] + \mathcal{V}_{LoadShiftingCostFix}[s, h, e] \Big) \end{align*}

Technology model

Technologies are a main functionality of any energy hub model. In ehubX, the set of all technologies is given by \(\mathcal{S}_{Tech}\). A technology may only be eligible for certain stages (due to its technology readiness level, see techs.yaml) or for certain hubs (due to tech_lists, see hubs.yaml). The stage-hub-tech tuples of all allowed combinations is captured in the set \(\mathcal{S}_{TechTuple} \subset \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{Tech}\). It contains all tuples \((s, h, x)\) where

• The technology \(x\) is ready in the stage \(s\), meaning the tech’s TRL at the stage’s start year is larger or equal to the TRL threshold defined in stages.yaml.

• The technology \(x\) occurs in one of the tech lists of hub \(h\).

The key feature of technologies is that they have a certain capacity

\begin{align*} \mathcal{V}_{TechCap}: \mathcal{S}_{TechTuple} \to \mathbb{R}_0^+ \end{align*}

This variable measures in an abstract way the amount of installed assets for each technology in each allowed stage and hub. This is a good moment to mention that the technology module in ehubX operates in a sense like an abstract module, meaning that it has no way of actually contributing to the balance of any energy carrier. Instead, this functionality is performed by the technology submodules which are

• conversion model,

• storage model,

• EBM model,

• ATES model,

• heatpump model.

The way this works internally is that when a technology is added to a model (e.g.; a storage technology), an entry is added to both the technology module and the storage module. The entry in the storage module handles storage-specific behavior whereas the entry in the technology module governs the more general, technology-related functionality described in this section.

Coming back to the capacity variable, each technology submodule has its own physical interpretation of the capacity variable, specifically:

Capacity variable for submodules

Submodule	Physical meaning of capacity	Example unit
Conversion	Output power	kW
Solar	Area of installation	\(m^2\)
Storage	Maximal storable energy	kWh
EBM	Total storage capacity	kWh
ATES	Area of installation	\(m^2\)
Heat pump	Condenser power	kW

Technology capacity is usually comprised of both initial technology (meaning it exists regardless of installation choices) and installed technology (chosen by the optimizer). The amount of installed technology is given by

\begin{align*} \mathcal{V}_{TechCapInstl}: \mathcal{S}_{TechCap} \to \mathbb{R}_0^+ \end{align*}

If a technology is installed in a certain stage, it will count towards the total capacity in this and future stages until its lifetime has run out. The same consideration has to be applied for initial technology whose lifetime will run out at a certain point as well. Therefore, a constraint exists which regulates the total available capacity considering initial and previous installation amounts:

\begin{align*} \mathcal{V}_{TechCap}[s, h, x] = Cap_{init}[s, h, x] + \sum\limits_{\substack{s_{instl} \in \mathcal{S}_{Stage} \\ (s_{instl}, h, x) \in \mathcal{S}_{TechTuple}}} Cap_{instl}[s_{instl}, s, h, x] \end{align*}

The initial capacity that is still operational in stage \(s\) is given by

\begin{align*} Cap_{init}[s, h, x] = \left \{ \begin{array}{rl} cap\_init[h, x], &\text{if } start\_year[s] - init\_year < lifetime[x] - age\_init[h, x] \\ 0, &\text{else} \end{array} \right . \end{align*}

Here, \(cap\_init\) and \(age\_init\) are initial capacity parameters from hubs.yaml, \(start\_year\) is the first year of a stage from stages.yaml and \(init\_year\) is the start year of the first stage (earliest \(start\_year\) in stages.yaml). The capacity installed in stage \(s_{instl}\) that is still operational in stage \(s\) is given by

\begin{align*} Cap_{instl}[s_{instl}, s, h, x] = \left \{ \begin{array}{rl} \mathcal{V}_{TechCapInstl}[s_{instl}, h, x], &\text{if } 0 \le start\_year[s] - start\_year[s_{instl}] < lifetime[x] \\ 0, &\text{else} \end{array} \right . \end{align*}

For the upcoming variable \(\mathcal{V}_{TechCostCapex}\), we require a binary variable that measures whether any amount of technology was installed at all. This variable is given by

\begin{align*} \mathcal{V}_{YTechCapInstl}: \mathcal{S}_{TechTuple} \to \{0, 1\} \end{align*}

This is achieved by a bigM constraint

\begin{align*} \mathcal{V}_{TechCapInstl}[s, h, x] \le BigM_{TechCap}[s, h, x] \cdot \mathcal{V}_{YTechCapInstl}[s, h, x] \end{align*}

where the big-M parameter has to provide an upper bound to the technology capacity and is chosen as

\begin{align*} BigM_{TechCap}[s, h, x] = \left \{ \begin{array}{rl} cap\_max[s, h, x] + 10^{-5}, &\text{if } cap\_max[s, h, x] < \infty \\ \mathcal{P}_{BigMGeneric}, &\text{else} \end{array} \right . \end{align*}

Here, \(cap\_max\) is the maximal capacity parameter from hubs.yaml and \({P}_{BigMGeneric}\) is the demand-based default big-M parameter introduced in the demand model.

Similar to \(\mathcal{V}_{YTechCapInstl}\), we require another binary variable to monitor whether a certain technology was used at all, later used for the cost variable \(\mathcal{V}_{TechCostOpexCap}\). This variable is declared in the technology module as

\begin{align*} \mathcal{V}_{YTechUsed}: \mathcal{S}_{TechTuple} \to \{0, 1 \} \end{align*}

Since the concept of usage may mean something different depending on the type of technology, each submodule declares its own constraint that ties this variable to a specific behavior pattern.

Technology capacity is possibly subject to some constraint based on the input parameters, starting with \(unit\_cap\_min\) from techs.yaml:

\begin{align*} \mathcal{V}_{TechCapInstl}[s, h, x] \ge unit\_cap\_min[s, x] \cdot \mathcal{V}_{YTechCapInstl}[s, h, x] \end{align*}

Another set of parameters that limits the capacity are \(cap\_min\) and \(cap\_max\) from hubs.yaml:

\begin{align*} cap\_min[s, h, x] \le \mathcal{V}_{TechCap}[s, h, x] \le cap\_max[s, h, x] \end{align*}

Furthermore, there is the parameter \(last_inst_year\) from techs.yaml which may forbid any installation after that year:

\begin{align*} \mathcal{V}_{YTechCapInstl}[s, h, x] = 0 \text{ if } start\_year[s] > last\_instl\_year[x] \end{align*}

For technologies, we allow the concept of coupled technologies (based on the section coupling_params in techs.yaml). All technologies with this field are collected in a subset \(\mathcal{S}_{SubTech} \subset \mathcal{S}_{Tech}\) and every \(x \in \mathcal{S}_{SubTech}\) is assigned a \(main\_tech\_id[x] \in \mathcal{S}_{Tech}\) with their associated main technology, and a \(cap\_factor[x] \in \mathbb{R}^+_0\) describing the factor relating the capacities between this tech and their main tech:

\begin{align*} \mathcal{V}_{TechCap}[s, h, x] = cap\_factor[x] \cdot \mathcal{V}_{TechCap}[s, h, main\_tech\_id[x]] \end{align*}

Certain costs are connected to the technology model, such as CAPEX costs:

\begin{align*} \mathcal{V}_{TechCostCapex}: \mathcal{S}_{TechTuple} &\to \mathbb{R}, \\ \mathcal{V}_{TechCostCapex}[s, h, x] &= CRF(interest\_rate[x], lifetime[x]) \cdot \sum\limits_{\substack{s_{instl} \in \mathcal{S}_{Stage} \\ (s_{instl}, h, x) \in \mathcal{S}_{TechTuple}}} Capex_{Stage}[s_{instl}, s, h, x] \end{align*}

Here \(CRF\) is a standard capital recovery factor calculated from the technology’s lifetime and interest rate as follows:

\begin{align*} CRF(i, N) = \frac{i \cdot (i+1)^N}{(i+1)^N-1} \end{align*}

The summand takes the following form:

\begin{align*} Capex_{Stage}[s_{instl}, s, h, x] ~=~ \left \{ \begin{array}{rl} capex\_per\_cap[s, x] \cdot \mathcal{V}_{TechCapInstl}[s_{instl}, h, x] & \\ +~ one\_time\_capex[s, x] \cdot \mathcal{V}_{YTechCapInstl}[s_{instl}, h, x], &\text{if } 0 \le start\_year[s] - start\_year[s_{instl}] < lifetime[x] \\ 0, &\text{else} \end{array} \right . \end{align*}

Here, \(capex\_per\_cap\) and \(one\_time\_capex\) are parameters from techs.yaml.

In addition to these installation-related CAPEX costs, operating and maintaining the technologies is associated with the OPEX costs parametrized by \(opex\_per\_cap\) and \(one\_time\_opex\) from techs.yaml:

\begin{align*} \mathcal{V}_{TechCostOpexCap}: \mathcal{S}_{TechTuple} &\to \mathbb{R}, \\ \mathcal{V}_{TechCostOpexCap}[s, h, x] &= opex\_per\_cap[s, x] \cdot \mathcal{V}_{TechCap}[s, h, x] + one\_time\_opex[s, x] \cdot \mathcal{V}_{YTechUsed}[s, h, x] \end{align*}

As can be seen, both costs have a part that arises per amount of installed capacity while another enters any time any amount of technology is installed or used at all. Adding up these costs gives us the following fixed variable:

\begin{align*} \mathcal{V}_{TechCostTotal} &\in \mathbb{R}, \\ \mathcal{V}_{TechCostTotal} &= \sum\limits_{(s, h, x) \in \mathcal{S}_{TechTuple}} \big( \mathcal{V}_{TechCostCapex}[s, h, x] + \mathcal{V}_{TechCostOpexCap}[s, h, x] \big) \end{align*}

Similar to the costs associated with tech installation, certain CO2-embodied emissions are associated with this process as well. We collect these in the variable

\begin{align*} \mathcal{V}_{TechCo2Instl}: \mathcal{S}_{TechTuple} &\to \mathbb{R} \\ \mathcal{V}_{TechCo2Instl}[s, h, x] &= \sum\limits_{\substack{s_{instl} \in \mathcal{S}_{Stage} \\ (s_{instl}, h, x) \in \mathcal{S}_{TechTuple} \\ 0 \le start\_year[s] - start\_year[s_{instl}] \le lifetime[x]}} \frac{co2\_per\_cap[s, x]}{lifetime[x]} \cdot \mathcal{V}_{TechCapInstl}[s_{instl}, h, x] \end{align*}

The parameter \(co\_per\_cap\) comes from the techs.yaml file. As with the costs, we gather the total amount of CO2 emissions into a fixed variable for convenience:

\begin{align*} \mathcal{V}_{TechCo2Total}: \mathcal{S} \to \mathbb{R}, \qquad \mathcal{V}_{TechCo2Total}[s] = \sum\limits_{\substack{(s', h, x) \in \mathcal{S}_{TechTuple} \\ s = s'}} \mathcal{V}_{TechCo2Instl}[s', h, x] \end{align*}

Conversion model

The conversion model defines its own set of tech tuples

\begin{align*} \mathcal{S}_{ConvTechTuple} = \big \{ (s, h, x) \in \mathcal{S}_{TechTuple}: ~ x \text{is a conversion technology } \big \} \end{align*}

In addition, a conversion technology may have multiple input and output ecs, as defined in techs.yaml and illustrated below for an example with three inputs and two outputs:

The model defines sets of input and output combinations

\begin{align*} \mathcal{S}_{ConvTechIn} &= \big \{ (s, h, x, e): (s, h, x) \in \mathcal{S}_{ConvTechTuple}, e \text{ is in\_ec of x } \big \} \\ \mathcal{S}_{ConvTechOut} &= \big \{ (s, h, x, e): (s, h, x) \in \mathcal{S}_{ConvTechTuple}, e \text{ is out\_ec of x } \big \} \end{align*}

Based on these sets, the following variables track intakes and outputs of conversion technologies at every time step:

\begin{align*} \mathcal{V}_{ConvTechIn}: \mathcal{S}_{ConvTechIn} \times \mathcal{S}_{Time} &\to \mathbb{R}_0^+ \\ \mathcal{V}_{ConvTechOut}: \mathcal{S}_{ConvTechOut} \times \mathcal{S}_{Time} &\to \mathbb{R}_0^+ \end{align*}

On the intake side, the ec amounts stand in a fixed relation to each other, expressed by the constraint

\begin{align*} \frac{\mathcal{V}_{ConvTechIn}[s, h, x, e, t]}{in\_part[s, x, e]} = \frac{\mathcal{V}_{ConvTechIn}[s, h, x, main\_in\_ec[x], t]}{in\_part[s, x, main\_in\_ec[x]]} \end{align*}

where \(in\_part\) is a parameter defined in techs.yaml. Similarily, the output amounts of each output ec are specified based on the efficiency parameter \(out\_eff\) from techs.yaml which relates to the input of the main input ec \(main\_in\_ec\):

\begin{align*} \mathcal{V}_{ConvTechOut}[s, h, x, e, t] = out\_eff[s, x, e, t] \cdot \mathcal{V}_{ConvTechIn}[s, h, x, main\_in\_ec[x], t] \end{align*}

From the tech model, we still have to define a constraint that determines the usage of a conversion technology.

\begin{align*} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{ConvTechOut}[s, h, x, main\_out\_ec[x], t] \le BigM_{ConvMainOutSum}[s, h, x] \cdot \mathcal{V}_{YTechUsed}[s, h, x] \end{align*}

Here, \(weight\) is a clustering parameter from reduced-order modeling and \(BigM_{ConvMainOutSum}[s, h, x]\) is a big-M parameter which has to provide an upper bound to the summed-up output of the conversion technology’s main output ec over the time horizon. It is chosen as

\begin{align*} BigM_{ConvMainOutSum}[s, h, x] = \left \{ \begin{array}{rl} out\_sum\_max[s, h, x], &\text{if } out\_sum\_max[s, h, x] < \infty \\ cap\_max[s, h, x] \cdot |\mathcal{S}_{TimeHorizon}| + 10^{-5}, &\text{elif } cap\_max[s, h, x] < \infty \\ \mathcal{P}_{BigMGeneric}, &\text{else} \end{array} \right . \end{align*}

We recall that \({P}_{BigMGeneric}\) is the demand-based default big-M parameter introduced in the demand model. Furthermore, \(out\_sum\_max\) and :math:`cap_max are parameters from hubs.yaml.

The amount of output a conversion technology can generate is not only limited by the technology model’s capacity*` but also by the *availability parameter from hubs.yaml. Since availability works as a percentage value, the constraint for this restriction looks as follows:

\begin{align*} \mathcal{V}_{ConvTechOut}[s, h, x, main\_out\_ec, t] \le availability[s, h, x, t] \cdot \mathcal{V}_{TechCap}[s, h, x] \end{align*}

Note that restricting the \(main\_out\_ec\) (see techs.yaml), all other output ecs are restricted accordingly due to the input-output dynamics defined above.

The techs.yaml offers the parameter \(out\_sum\_min\) and \(out\_sum\_max\) to limit the output amounts of the main output ec over the entire time horizon. Using the parameter \(weight\) from the clustering logic, this is formulated as follows:

\begin{align*} out\_sum\_min[s, h, x] \le \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{ConvTechOut}[s, h, x, out\_ec\_main[x], t] \le out\_sum\_max[s, h, x] \end{align*}

Lastly, the conversion technology model has its own cost concept in addition to the costs defined in the tech model, namely OPEX costs that arise proportionally to the amount of energy that is being output by the conversion technologies. The calculation of these costs uses the cost parameter \(opex\_per\_energy\) from techs.yaml:

\begin{align*} \mathcal{V}_{ConvTechCostOpexOut}&: \mathcal{S}_{ConvTechTuple} \to \mathbb{R}, \\ \mathcal{V}_{ConvTechCostOpexOut}[s, h, x] &= \sum\limits_{\substack{e \in \mathcal{S}_{Ec} \\ (s, h, x, e) \in \mathcal{S}_{ConvTechOut}}} opex\_per\_energy[s, h, x] \cdot \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{ConvTechOut}[s, h, x, e, t] \end{align*}

Finally, all costs in the conversion model (i.e.; only the energy-related OPEX costs above) are bundled into a single variable for convenient use outside of the conversion module:

\begin{align*} \mathcal{V}_{ConvTechCostTotal} &\in \mathbb{R}, \\ \mathcal{V}_{ConvTechCostTotal} &= \sum\limits_{(s, h, x) \in \mathcal{S}_{ConvTechTuple}} \mathcal{V}_{ConvTechCostOpexOut}[s, h, x] \end{align*}

Solar model

Solar technologies are a special case of conversion technologies, converting solar irradiation into electrical power. The conversion model defines its own set of tech tuples

\begin{align*} \mathcal{S}_{SolarTechTuple} = \big \{ (s, h, x) \in \mathcal{S}_{ConvTechTuple}: x \text{ is a solar technology} \big \} \end{align*}

Additionally, solar technologies may only have a single energy carrier which is the solar irradition. Within this section, we therefore simply write \(e_{sol}[x]\) for this energy carrier. If a conversion technology with more than one input ec is designated as solar, ehubX will throw an exception.

The main input variable associated with solar technologies is the solar tech incident which relates to the total amount of solar irradiation that is collected through the installed technology, as follows:

\begin{align*} \mathcal{V}_{SolarTechIncident}[s, h, x, t] = solar\_irradiation[s, e_{sol}[x], t] \cdot \mathcal{V}_{TechCap}[s, h, x] \end{align*}

The parameter \(solar\_irradiation\) comes from the file solar_irradiation.csv which contains an irradiation amount per \(m^2\) which is why the capacity variable \(\mathcal{V}_{TechCap}\) for solar technologies is understood in \(m^2\) as well. In order for this to work, the solar model disables the conversion-specific capacity limitation constraint (there understood in \(kW\)) and defines its own version of this as follows:

\begin{align*} \mathcal{V}_{ConvTechIn}[s, h, x, e_{sol}[x], t] \le availability[s, h, x, t] \cdot \mathcal{V}_{SolarTechIncident}[s, h, x, t] \end{align*}

In contrast to the conversion model, it is the single input ec which is limited proportionally by the availability. Another difference lies in the fact that production of solar electricity can not always be freely chosen between zero and the upper limit from the previous constraint. Instead, curtailment must be employed to down-regulate the production of PV panels which is subject to the parameter \(curtail\_max\_rel\) from techs.yaml:

\begin{align*} \mathcal{V}_{ConvTechIn}[s, h, x, e_{sol}[x], t] \ge (1 - curtail\_max\_rel[s, x]) \cdot availability[s, h, x, t] \cdot \mathcal{V}_{SolarTechIncident}[s, h, x, t] \end{align*}

Finally, we need to limit installed solar capacity (i.e.; area) by the total available area as defined in the file solar_areas.csv. This file is parsed into a parameter \(solar\_area\) for which the following constraint has to hold:

\begin{align*} \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x) \in \mathcal{S}_{SolarTechTuple} \\ e = e_{sol}[x]}} \mathcal{V}_{TechCap}[s, h, x] \le solar\_area[s, h, e] \end{align*}

Storage model

The storage model defines its own set of technologies

\begin{align*} \mathcal{S}_{StorTech} = \big \{ x \in \mathcal{S}_{Tech}: ~ x \text{is a storage technology } \big \} \end{align*}

and storage tech tuples

\begin{align*} \mathcal{S}_{StorTechTuple} = \big \{ (s, h, x) \in \mathcal{S}_{TechTuple}: ~ x \in \mathcal{S}_{StorTech} \big \} \end{align*}

In addition, a storage technology has exactly one ec that can be stored, defined as the parameter

\begin{align*} \mathcal{P}_{StorTechEc}: \mathcal{S}_{StorTech} \to \mathcal{S}_{Ec} \end{align*}

This mapping from storage technology to ec works using the parameter \(ec\) from the section \(storage_params\) in techs.yaml. During every time step, a storage technology is able to charge and discharge its stored ec using the variables

\begin{align*} \mathcal{V}_{StorTechInflow}, \mathcal{V}_{StorTechOutflow}: \mathcal{S}_{StorTechTuple} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \end{align*}

The extent to which this is possible can be limited by the parameters \(charge\_max\) and \(discharge\_max\) from techs.yaml:

\begin{align*} \mathcal{V}_{StorTechInflow}[s, h, x, t] &\le charge\_max[s, x] \cdot \mathcal{V}_{TechCap}[s, h, x] \\ \mathcal{V}_{StorTechOutflow}[s, h, x, t] &\le discharge\_max[s, x] \cdot \mathcal{V}_{TechCap}[s, h, x] \end{align*}

Recall that \(\mathcal{V}_{TechCap}\) is the total installed capacity from the tech model which, in case of the storage model, takes on the meaning of storage capacity, i.e.; total storable energy. Therefore, setting e.g. the parameter \(charge\_max\) to 0.1 means that the storage technology would be able to charge by 10% of its total storage capacity at every time step, therefore it would take at least 10 time step to fully charge the storage.

The charging and discharging amounts can also be used to determine whether a storage technology was used in the sense of the tech model by adding the constraint

\begin{align*} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] (\mathcal{V}_{StorTechInflow}[s, h, x, t] + \mathcal{V}_{StorTechOutflow}[s, h, x, t]) \le BigM_{StorTotalFlow}[s, h, x] \cdot \mathcal{V}_{YTechUsed}[s, h, x] \end{align*}

Ideally, we would like to make use of the \(charge\_max\) and \(discharge\_max\) parameters for \(BigM_{StorTotalFlow}[s, h, x]\) and therefore choose the following form:

\begin{align*} BigM_{StorTotalFlow}[s, h, x] = \left \{ \begin{array}{rl} 10^{-5} + cap\_max[s, h, x] \cdot \big( \min(charge\_max[s, x], ~1) + &\\ \min(discharge\_max[s, x], ~1) \big) \cdot |\mathcal{S}_{TimeHorizon}|, &\text{if } cap\_max[s, h, x] < \infty \\ \mathcal{P}_{BigMGeneric}, &\text{else} \end{array} \right . \end{align*}

Next, we introduce a storage energy variable to the model which keeps track of how much energy is currently stored in a storage technology:

\begin{align*} \mathcal{V}_{StorTechEnergy}: \mathcal{S}_{StorTechTuple} \times \mathcal{S}_{TimeHorizon} \to \mathbb{R}_0^+ \end{align*}

Note that this variable is not defined on the (possibly clustered) time domain \(\mathcal{S}_{Time}\) but on the full time horizon \(\mathcal{S}_{TimeHorizon}\). This is done because the dynamics imposed on the variable will rely on comparing adjacent time steps, something that is not possible for clustered time steps. The storage energy’s evolution is subject to a standby loss and the charging and discharging process:

\begin{align*} \mathcal{V}_{StorTechEnergy}[s, h, x, t_h^+] = &(1 - standby\_loss[s, x]) \cdot \mathcal{V}_{StorTechEnergy}[s, h, x, t_h] \\ &+ in\_eff[s, x] \cdot \mathcal{V}_{StorTechInflow}[s, h, x, t\_cl[t_h]] - out\_eff[s, x]^{-1} \cdot \mathcal{V}_{StorTechOutflow}[s, h, x, t\_cl[t_h]] \end{align*}

The parameters \(standby\_loss\), \(in\_eff\) and \(out_eff\) come from techs.yaml. It is important to note that the above constraint is formulated for every horizon time step \(t_h \in \mathcal{S}_{TimeHorizon}\). We have used the parameter \(t\_cl[t_h]\) from the clustering procedure which assigns every horizon time step \(t_h\) to its clustering time step. There is also the concept of a next time step which we define as

\begin{align*} t_h^+ = \left \{ \begin{array}{rl} \text{first horizon time step}, &\text{if } t_h \text{ is the last horizon time step } \\ \text{horizon time step after } t_h, &\text{else} \end{array} \right . \end{align*}

This way, we use a circular time model so no storage energy can be procued or destroyed over the time horizon. This assumption is fairly common in energy system models where no long-term energy storage effects have to be taken into account. For example, if the time horizon \(\mathcal{S}_{TimeHorizon}\) encompasses a single year, storage technologies will be able to shift around energy throughout the year but not over multiple years.

The storage energy levels are additionally constrained by the parameters \(soc\_min\) and \(soc\_max\) from \(techs_yaml\):

\begin{align*} soc\_min[s, x] \cdot \mathcal{V}_{TechCap}[s, h, x] \le \mathcal{V}_{StorTechEnergy}[s, h, x, t_h] \le \min(soc\_max[s, x], 1) \cdot \mathcal{V}_{TechCap}[s, h, x] \end{align*}

For a final consideration, the energy storage level is currently able to ‘start’ at any initial value on the time horizon but has to follow the prescribed dynamics afterwards. Depending on the model at hand, this behavior may be desirable or we may want to set an initial value for the storage level. ehubX handles this by considering the parameter \(soc\_init\) from hubs.yaml: If this parameter is infinite (which is its default value), no constraint will be set. Otherwise, the parameter will be used to fix the initial energy storage value:

\begin{align*} \mathcal{V}_{StorTechEnergy}[s, h, x, t_{h,0}] = \left \{ \begin{array}{rl} \mathcal{V}_{StorTechEnergy}[s_0, h, x, t_{h,0}], &\text{if } s \neq s_0 \\ soc\_init[h, x] \cdot \mathcal{V}_{TechCap}[s, h, x], &\text{if } s = s_0 \text{ and } soc\_init[h, x] < \infty \end{array} \right . \end{align*}

Here, \(s_0 \in \mathcal{S}_{Stage}\) is the first stage in the model, i.e.; the one with the smallest \(start\_year\), and \(t_{h,0} \in \mathcal{S}_{TimeHorizon}\) is the initial horizon time step. Note that regardless of whether an initial energy level is set in the first stage, the model has to maintain the same initial storage level across all stages.

It has to be mentioned at this point that the chosen storage formulation would allow for simultaneous inflow and outflow from the storage technology (i.e.; \(mathcal{V}_{StorTechInflow}[s, h, x, t] > 0\) and \(mathcal{V}_{StorTechOutflow}[s, h, x, t] > 0\)). It is generally possible to avoid this by introducing a time-dependent binary variable which monitors whether a charging or discharing process is taking place at the current time step. However, this drastically increases the complexity of the resulting MILP problem which is why ehubX decides to avoid this. Instead, we rely on the fact that if certain costs are associated with the flow variables, an optimal solution will always try to avoid situations where both variables take nonzero values simulteneously. For example, take the case of \(\mathcal{V}_{StorTechOutflow} - \mathcal{V}_{StorTechInflow} = 3\) which could be comprised of the following two value pairs for the sub-variables:

1. \(\mathcal{V}_{StorTechOutflow}[s, h, e, t] = 5\) and \(\mathcal{V}_{StorTechInflow}[s, h, e, t] = 2\)

2. \(\mathcal{V}_{StorTechOutflow}[s, h, e, t] = 3\) and \(\mathcal{V}_{StorTechInflow}[s, h, e, t] = 0\)

If obtaining an amount of the ec \(e\) is associated with certain cost within the energy system (e.g.; through costs in the import model) and either the parameters \(in\_eff[s, x]\) or \(out\_eff[s, x]\) are smaller than 1, it will be undesirable from a mathematical perspective to operate the storage technology more than necessary, meaning the case b) above would be preferred to a). A similar reasoning as the one employed here is used in the load shifting model for the variables \(\mathcal{V}_{LoadShiftingAbove}\) and \(\mathcal{V}_{LoadShiftingBelow}\).

EBM model

First, we note that the electricity-based mobility (EBM) technology shares a lot of design overlap with the storage model since electric vehicles can essentially be thought of as batteries with an added consumption behavior.

The EBM model defines its own set of technologies

\begin{align*} \mathcal{S}_{EbmTech} = \big \{ x \in \mathcal{S}_{Tech}: ~ x \text{is an EBM technology } \big \} \end{align*}

and EBM tech tuples

\begin{align*} \mathcal{S}_{EbmTechTuple} = \big \{ (s, h, x) \in \mathcal{S}_{TechTuple}: ~ x \in \mathcal{S}_{EbmTech} \big \} \end{align*}

The first deviation from the storage technology comes from the fact that the capacity of the EBM fleet is a fixed value that cannot be changed by the optimizer, originating from the parameters num_vehicles (from hubs.yaml) and storage_cap (from techs.yaml):

\begin{align*} \mathcal{V}_{TechCap}[s, h, x] = num\_vehicles[s, h, x] \cdot storage\_cap[s, x] \end{align*}

All capacity-related parameters from the technology model can still be specified and will work as intended yet their values will have no outcome on the optimization results. Note specifically that if no cap_init is specified in hubs.yaml, ehubX will install the necessary capacity itself, thereby modifying the total system cost.

Moving on, an EBM technology has exactly one ec that operates it, defined as the parameter

\begin{align*} \mathcal{P}_{EbmTechEc}: \mathcal{S}_{EbmTech} \to \mathcal{S}_{Ec} \end{align*}

This mapping from EBM technology to ec works using the parameter \(ec\) from the section \(ebm_params\) in techs.yaml. During every time step, an EBM technology can charge and discharge its stored ec using the variables

\begin{align*} \mathcal{V}_{EbmTechInflow}, \mathcal{V}_{EbmTechOutflow}: \mathcal{S}_{EbmTechTuple} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \end{align*}

The extent to which this is possible is firstly limited by the parameters \(charge\_max\) and \(discharge\_max\) from techs.yaml. Second, the time-dependent parameter availability from hubs.yaml specifies the percentage of the fleet that is available for charging and discharging at any given point. Third and final, the parameter discharge_controllability from techs.yaml further dampens the fleet’s ability to discharge at will. This results in the following constraints:

\begin{align*} \mathcal{V}_{EbmTechInflow}[s, h, x, t] &\le num\_vehicles[s, h, x] \cdot availability[s, h, x] \cdot charge\_max[s, x] \\ \mathcal{V}_{StorTechOutflow}[s, h, x, t] &\le num\_vehicles[s, h, x] \cdot availability[s, h, x] \cdot discharge\_controllability[s, x] \cdot discharge\_max[s, x] \end{align*}

The charging and discharging amounts can also be used to determine whether an EBM technology was used in the sense of the tech model by adding the constraint

\begin{align*} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] (\mathcal{V}_{EbmTechInflow}[s, h, x, t] + \mathcal{V}_{EbmTechOutflow}[s, h, x, t]) \le BigM_{EbmTotalFlow}[s, h, x] \cdot \mathcal{V}_{YTechUsed}[s, h, x] \end{align*}

Due to the constraints above this, the big-M parameter can be chosen as

\begin{align*} BigM_{EbmTotalFlow}[s, h, x] = &10^{-5} + num\_vehicles[s, h, x] \cdot \Big( \min(charge\_max[s, x], storage\_cap[s, x]) \\ &+ \min(discharge\_max[s, x], storage\_cap[s, x]) \Big) \cdot |\mathcal{S}_{TimeHorizon}| \end{align*}

To avoid the case where \(charge\_max\) or \(discharge\_max\) are infinite, we have made additional use of the fact that an electric vehicle can never charge or discharge more than its storage capacity \(storage\_cap\) in a single timestep due to model limitations.

Next, we introduce an energy variable to the model which keeps track of how much energy is currently stored in the entire EBM fleet:

\begin{align*} \mathcal{V}_{EbmTechEnergy}: \mathcal{S}_{EbmTechTuple} \times \mathcal{S}_{TimeHorizon} \to \mathbb{R}_0^+ \end{align*}

Note that this variable is not defined on the (possibly clustered) time domain \(\mathcal{S}_{Time}\) but on the full time horizon \(\mathcal{S}_{TimeHorizon}\). This is done because the dynamics imposed on the variable will rely on comparing adjacent time steps, something that is not possible for clustered time steps. The storage energy’s evolution is subject to a standby loss, the charging and discharging process and the fleet’s energy consumption due to the demand:

\begin{align*} \mathcal{V}_{EbmTechEnergy}[s, h, x, t_h^+] = &(1 - standby\_loss[s, x]) \cdot \mathcal{V}_{EbmTechEnergy}[s, h, x, t_h] \\ &+ in\_eff[s, x] \cdot \mathcal{V}_{EbmTechInflow}[s, h, x, t\_cl[t_h]] - out\_eff[s, x]^{-1} \cdot \mathcal{V}_{EbmTechOutflow}[s, h, x, t\_cl[t_h]] \\ &- num\_vehicles[s, h, x] \cdot demand\_modifier[s, h, x] \cdot demand\_nominal[s, h, x, t] \end{align*}

The parameters \(standby\_loss\), \(in\_eff\) and \(out_eff\) come from techs.yaml whereas \(num\_vehicles\), \(demand_modifier\) and \(demand_nominal\) can be specified in hubs.yaml. It is important to note that the above constraint is formulated for every horizon time step \(t_h \in \mathcal{S}_{TimeHorizon}\). We have used the parameter \(t\_cl[t_h]\) from the clustering procedure which assigns every horizon time step \(t_h\) to its clustering time step. There is also the concept of a next time step which we define as

\begin{align*} t_h^+ = \left \{ \begin{array}{rl} \text{first horizon time step}, &\text{if } t_h \text{ is the last horizon time step } \\ \text{horizon time step after } t_h, &\text{else} \end{array} \right . \end{align*}

This way, we use a circular time model so no EBM energy can be procued or destroyed over the time horizon. This assumption is fairly common in energy system models where no long-term energy storage effects have to be taken into account.

The EBM energy levels are additionally constrained by the parameters \(soc\_min\) and \(soc\_max\) from \(techs_yaml\):

\begin{align*} \mathcal{V}_{EbmTechEnergy}[s, h, x, t_h] &\ge soc\_min[s, x] \cdot num\_vehicles[s, h, x] \cdot storage\_cap[s, x] \\ \mathcal{V}_{EbmTechEnergy}[s, h, x, t_h] &\le soc\_max[s, x] \cdot num\_vehicles[s, h, x] \cdot storage\_cap[s, x] \end{align*}

For a final consideration, the energy level is currently able to ‘start’ at any initial value on the time horizon but has to follow the prescribed dynamics afterwards. Depending on the model at hand, this behavior may be desirable or we may want to set an initial value for the energy level. ehubX handles this by considering the parameter \(soc\_init\) from hubs.yaml: If this parameter is infinite (which is its default value), no constraint will be set. Otherwise, the parameter will be used to fix the initial energy storage value:

\begin{align*} \mathcal{V}_{EbmTechEnergy}[s, h, x, t_{h,0}] = \left \{ \begin{array}{rl} \mathcal{V}_{EbmTechEnergy}[s_0, h, x, t_{h,0}], &\text{if } s \neq s_0 \\ soc\_init[h, x] \cdot num\_vehicles[s, h, x] \cdot storage\_cap[s, x], &\text{if } s = s_0 \text{ and } soc\_init[h, x] < \infty \end{array} \right . \end{align*}

Here, \(s_0 \in \mathcal{S}_{Stage}\) is the first stage in the model, i.e.; the one with the smallest \(start\_year\), and \(t_{h,0} \in \mathcal{S}_{TimeHorizon}\) is the initial horizon time step. Note that regardless of whether an initial energy level is set in the first stage, the model has to maintain the same initial storage level across all stages.

For a final point, the current model allows for a simultaneous charging and discharging process of EBM vehicles (\(\mathcal{V}_{EbmTechInflow}[s, h, x, t] > 0\) and \(\mathcal{V}_{EbmTechOutflow}[s, h, x, t] > 0\)). We would like to refer to the final section of the storage model because the considerations expressed there hold true in the exact same way for the EBM model.

ATES model

We would like to precede this chapter on ATES technology by mentioning a publication where the mathematical model has been first introduced by the authors of ehubX in [Ates25]:

[Ates25]

4. Beermann, S.N. Sørensen, M. Tønder, C. Doughty, P. Blum, K. Menberg, M. Wetter, M. Sulzer, R. Mutschler, “Integration of Aquifer Thermal Energy Storage (ATES) to linear energy system models”, submitted to Applied Energy, 2026.

The ATES model defines its own set of tech tuples

\begin{align*} \mathcal{S}_{AtesTechTuple} = \big \{ (s, h, x) \in \mathcal{S}_{TechTuple}: ~ x \text{ is an ATES technology } \big \} \end{align*}

A key characteristic of an ATES is that any given year will contain two time periods: One during which the fluid is pumped from warm wells to cold wells (warm-to-cold phase) and another where the direction is reversed (cold-to-warm phase). We call this division of the horizon an ATES schedule. The decision when to switch from one phase to the other cannot be fully optimized by ehubX while still maintaining a linear model. Instead, we offer a set of schedules per hub (cf. schedules in hubs.yaml), collected in the set \(\mathcal{S}_{AtesSchedule}\). This leads to the set of ATES technologies per schedule:

\begin{align*} \mathcal{S}_{AtesTechTupleSchedule} = \big \{ (s, h, x, i) \in \mathcal{S}_{AtesTechTuple} \times \mathcal{S}_{AtesSchedule}: i \text{ is an ATES schedule of } h \big \} \end{align*}

The technology capacity \(\mathcal{V}_{TechCap}\) for ATES technologies will be split into the different schedules by ehubX, thereby allowing for a limited amount of time schedule optimization. This is done via an introduction of the variable

\begin{align*} \mathcal{V}_{AtesTechCapSchedule}: \mathcal{S}_{AtesTechTupleSchedule} \to \mathbb{R}_0^+ \end{align*}

The capacity subdivison itself is organized in a straightforward manner by the following constraint for all \((s, h, x) \in \mathcal{S}_{AtesTechTuple}\):

\begin{align*} \mathcal{V}_{TechCap}[s, h, x] = \sum\limits_{\substack{i \in \mathcal{S}_{AtesSchedule} \\ (s, h, x, i) \in \mathcal{S}_{AtesTechTupleSchedule}}} \mathcal{V}_{AtesTechCapSchedule}[s, h, x, i] \end{align*}

Similar to the solar model, the ATES model interprets the technology capacity \(\mathcal{V}_{TechCap}\) (and \(\mathcal{V}_{AtesTechCapSchedule}\)) as the amount of area that is covered by installed ATES technologies. Therefore, we also need to respect the total available area which is quantified by the parameter available_area from hubs.yaml:

\begin{align*} \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x) \in \mathcal{S}_{AtesTechTuple}}} \mathcal{V}_{TechCap}[s, h, x] \le available\_area[s, h]. \end{align*}

Next, we now turn our attention towards the main dynamics and functionality of an ATES system. In ehubX, the ATES model assumes that all wells are organized in pairs of warm and cold wells. It has to be noted that this is not the only way to design an ATES system, and additional configuration options may be added in the future. Additionally, the optimizer might choose to install not one but a certain number of well pairs which we will detail further below. For now, let us consider a single well pair to illustrate the model:

In cooling mode during the cold-to-warm phase, a heat exchanger transmits heating energy to the cold fluid to heat it up, after which it is sent to the warm well. The thermal power at which this happens is denoted in the graphic by \(P_{Cool}\) and will be handled in the model by the variable \(\mathcal{V}_{AtesTechCoolSchedule}\) (see further below). In heating mode during the warm-to-cold phase, the heat exchanger extracts energy from the warm fluid to cool it down, after which it is sent to the cold well. In this case, the thermal power is denoted in the graphic by \(P_{Heat}\) and will be handled in the model by the variable \(\mathcal{V}_{AtesTechHeatSchedule}\) (see futher below). In both cases, the electricity consumed by the pumps and the heat exchanger for the entire process is refered to as \(P_{Elec}\) in the graphic.

While the amounts \(P_{Cool}\), \(P_{Heat}\) and \(P_{Elec}\) belong to a single pair of wells, the ATES model in ehubX considers a group of well pairs whose size might be chosen by the optimizer. Therefore, we introduce the following variables:

\begin{align*} \mathcal{V}_{AtesTechCoolSchedule}:~ &\mathcal{S}_{AtesTechTupleSchedule} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \\ \mathcal{V}_{AtesTechHeatSchedule}:~ &\mathcal{S}_{AtesTechTupleSchedule} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \\ \mathcal{V}_{AtesTechElecSchedule}:~ &\mathcal{S}_{AtesTechTupleSchedule} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \end{align*}

These variables quantify how much power is produced (for heating and cooling) or consumed (for electricity) for every combination of stage, hub, ATES technology, ATES schedule and time step. Obviously, heating is only permitted in warm-to-cold phases and cooling only in cold-to-warm phases which leads to the constraints

\begin{align*} \mathcal{V}_{AtesTechHeatSchedule}[s, h, x, i, t] = 0 \text{ if } t \text{ is not in } i\text{'s warm-to-cold phase} \\ \mathcal{V}_{AtesTechCoolSchedule}[s, h, x, i, t] = 0 \text{ if } t \text{ is not in } i\text{'s cold-to-warm phase} \end{align*}

Furthermore, the relation between heating/cooling outputs and the electricity consumption is defined by the parameters elec_per_energy_heat and elec_per_energy_cool from the ATES technology data (cf. hubs.yaml):

\begin{align*} \mathcal{V}_{AtesTechElecSchedule}[s, h, x, i, t] = &elec\_per\_energy\_heat[s, h, x] \cdot \mathcal{V}_{AtesTechHeatSchedule}[s, h, x, i, t] \\ +~ &elec\_per\_energy\_cool[s, h, x] \cdot \mathcal{V}_{AtesTechCoolSchedule}[s, h, x, i, t] \end{align*}

It might be desireable to keep an ATES system in balance between its warm-to-cold and cold-to-warm periods, as not to completely deplete one of the wells in favor of the other. To achieve this, we use the parameters max_heat_over_cool and max_cool_over_heat from hubs.yaml to formulate the following constraints:

\begin{align*} \sum_{t \in \mathcal{S}_{Time}} \mathcal{V}_{AtesTechHeatSchedule}[s, h, x, i, t] \le max\_heat\_over\_cool[s, h, x, i] \cdot \sum_{t \in \mathcal{S}_{Time}} \mathcal{V}_{AtesTechCoolSchedule}[s, h, x, i, t] \\ \sum_{t \in \mathcal{S}_{Time}} \mathcal{V}_{AtesTechCoolSchedule}[s, h, x, i, t] \le max\_cool\_over\_heat[s, h, x, i] \cdot \sum_{t \in \mathcal{S}_{Time}} \mathcal{V}_{AtesTechHeatSchedule}[s, h, x, i, t] \end{align*}

Continuing with the next ATES property, the produced thermal power(s) need to be limited by the amount of installed capacity (i.e.; covered aquifer area) and the availability (an ATES parameter from hubs.yaml) at every timestep. The following constraints will relate the generated thermal powers to the concept of area capacity, respecting the availability:

\begin{align*} \mathcal{V}_{AtesTechHeatSchedule}[s, h, x, i, t] \le &~ max\_power\_density\_heat[s, h, x, i] \cdot availability[s, h, x, i, t] \\ &* \mathcal{V}_{AtesTechCapSchedule}[s, h, x, i] \\ \mathcal{V}_{AtesTechCoolSchedule}[s, h, x, i, t] \le &~ max\_power\_density\_cool[s, h, x, i] \cdot availability[s, h, x, i, t] \\ &* \mathcal{V}_{AtesTechCapSchedule}[s, h, x, i] \end{align*}

These equations contain maximal power densities for heating and cooling which are calculated as

\begin{align*} max\_power\_density\_heat[s, h, x, i] = \frac{max\_heat\_power\_per\_well\_pair[s, h, x, i]}{area\_per\_well\_pair[s, h, x, i]} \end{align*}

and similarly for cooling (in general, the following equations are calculated identically for heating and cooling so we will only write down the heating version here to keep things concise). Starting with the nominator, the maximal heating power that a well pair can produce is calculated as

\begin{align*} max\_heat\_power\_per\_well\_pair[s, h, x, i] = &~density\_fluid[x] \cdot specific\_heat\_capacity\_fluid[x] \\ &\cdot max\_pump\_rate\_per\_warm\_well[s, h, x, i] \cdot max\_temperature\_spread\_warm[h] \end{align*}

The parameters \(density\_fluid\), \(specific\_heat\_capacity\_fluid\), and \(max\_temperature\_spread\_warm\) are input parameters. The parameter \(max\_pump\_rate\_per\_warm\_well\) is also a potential input parameter. However, if it is not specified, ehubX will try to alternatively calculate a maximal pumping rate based on the Cooper-Jacobs approximation of the Theis equation, using the parameters \(well\_radius\), \(well\_distance\), \(hydraulic\_conductivity\_aquifer\), \(thickness\_aquifer\), and \(max\_drawdown\). The details of this calculation are carried out in [Ates25].

Returning to the maximal power density equation above, the denominator is given by the thermally affected area of a well pair. ehubX will start the calculation of this property by considering the circular areas of each well that are thermally affected. The radius of these circles is called the thermal radius of a well. This thermal radius can either be given as an input parameter (thermal_radius_per_warm_well and thermal_radius_per_cold_well). If it is not specified, ehubX will try to calculate it based on the input parameters specific_heat_capacity_fluid*`, max_pump_rate_per_warm_well, specific_heat_capacity_aquifer, thickness_aquifer, and darcy_velocity. For the details of this calculation, we refer again to [Ates25].

Given the now calculated areas of the two wells, ehubX offers different options for how these will combine to the total area of a well pair. The calculation method can be set by the input parameter well_pair_area_calculation_method, leading to:

\begin{align*} area\_per\_well\_pair[s, h, x, i] = \left \{ \begin{array}{rl} \pi \cdot (r_{th,warm}[s, h, x, i]^2 + r_{th,cold}[s, h, x, i]^2), &\text{if method = 'two\_circles'} \\ \max(r_{th,warm}, r_{th,cold}) \cdot (r_{th,warm} + r_{th,cold}), &\text{if method = 'smallest\_rectangle'} \end{array} \right . \end{align*}

We can see from the previous equations that the area for a well pair is generally not easily known at the input level, and instead calculated by ehubX based on a multitude of other parameters. Therefore, it will be hard for a user to correctly specify certain parameters from the tech model in the conventional way. For example, the parameter capex_per_cap would in this case need to be given in the units \(CHF/m^2\) which is fairly unconventional for energy system designers. Instead, a system designer would rather specify the price for each pair of wells. For this reason, we deactivate the conventional CAPEX, OPEX and CO2 calculations from the tech model and offer the parameters capex_per_well_pair, opex_per_well_pair, and co2_per_well_pair instead. These parameters correspond to the technology parameters capex_per_cap, opex_per_cap, and co2_per_cap which are approximated dynamically as, e.g.;

\begin{align*} capex\_per\_cap[s, h, x] \approx capex\_per\_well\_pair[s, x] / area\_per\_well\_pair\_max[s, h, x], \end{align*}

and accordingly for OPEX and CO2. The parameter \(area\_per\_well\_pair\_max[s, h, x]\) comes from an area calculation which is very similar to \(area\_per\_well\_pair\) above, using maximal-area estimations. For full details of this calculation, we refer again to [Ates25]. Additionally, we directly incorporate the number of well pairs per schedule in the module by introducing the fixed variable

\begin{align*} \mathcal{V}_{AtesTechNumWellPairs}:~ \mathcal{S}_{AtesTechTupleSchedule} &\to \mathbb{R}_0^+, \\ \mathcal{V}_{AtesTechNumWellPairs}[s, h, x, i] &= \frac{\mathcal{V}_{AtesTechCapSchedule}[s, h, x, i]}{area\_per\_well\_pair[s, h, x, i]} \end{align*}

Using the parameters well_pairs_min and well_pairs_max from hubs.yaml, we can then add the constraints

\begin{align*} well\_pairs\_min[s, h, x, i] \le \mathcal{V}_{AtesTechNumWellPairs}[s, h, x, i] \le well\_pairs\_max[s, h, x, i] \end{align*}

Finally, the variables \(\mathcal{V}_{AtesTechCoolSchedule}\), \(\mathcal{V}_{AtesTechHeatSchedule}\) and \(\mathcal{V}_{AtesTechElecSchedule}\) need to be capsulated into central input and output variables for the energy system model. These are indexed by the sets

\begin{align*} \mathcal{S}_{AtesTechOut} &= \big \{ (s, h, x, e) \in \mathcal{S}_{AtesTechTuple} \times \mathcal{S}_{Ec}: e \in \{ec\_heat[x], ec\_cool[x] \} \big \} \\ \mathcal{S}_{AtesTechIn} &= \big \{ (s, h, x, e) \in \mathcal{S}_{AtesTechTuple} \times \mathcal{S}_{Ec}: e = ec\_elec[x] \big \} \end{align*}

where ec_heat, ec_cool and ec_elec are the ecs specified in the ates_params|ecs block of the ATES technology x in techs.yaml. Based on these set, the central variables are defined as

\begin{align*} \mathcal{V}_{AtesTechOut}: \mathcal{S}_{AtesTechOut} \times \mathcal{S}_{Time} &\to \mathbb{R}_0^+, \\ \mathcal{V}_{AtesTechOut}[s, h, x, ec\_heat[x], t] &= \sum\limits_{\substack{i \in \mathcal{S}_{AtesSchedule} \\ (s, h, x, i) \in \mathcal{S}_{AtesTechHeatSchedule}}} \mathcal{V}_{AtesTechHeatSchedule}[s, h, x, i, t] \\ \mathcal{V}_{AtesTechOut}[s, h, x, ec\_cool[x], t] &= \sum\limits_{\substack{i \in \mathcal{S}_{AtesSchedule} \\ (s, h, x, i) \in \mathcal{S}_{AtesTechCoolSchedule}}} \mathcal{V}_{AtesTechCoolSchedule}[s, h, x, i, t] \end{align*}

Similarly, the central input variable is defined as

\begin{align*} \mathcal{V}_{AtesTechIn}: \mathcal{S}_{AtesTechIn} \times \mathcal{S}_{Time} &\to \mathbb{R}_0^+, \\ \mathcal{V}_{AtesTechIn}[s, h, x, ec\_elec[x], t] &= \sum\limits_{\substack{i \in \mathcal{S}_{AtesSchedule} \\ (s, h, x, i) \in \mathcal{S}_{AtesTechElecSchedule}}} \mathcal{V}_{AtesTechElecSchedule}[s, h, x, i, t] \end{align*}

Heat pump model

The heat pump model first defines its own set of tech tuples

\begin{align*} \mathcal{S}_{HpTechTuple} = \big \{ (s, h, x) \in \mathcal{S}_{TechTuple}: ~ x \text{ is a heat pump technology } \big \} \end{align*}

Next, we need to understand that a heat pump has multiple inputs and outputs, as illustrated by the figure below:

The heat pump can operate either in heating or cooling mode but the interior workings are always comparable: The evaporator is connected to an input medium at temperature \(T_{Evap,in}\) which is cooled down to \(T_{Evap,out}\). This absorbed heating power \(P_{Evap}\) is compressed using an input power \(P_{Elec}\). On the condenser side, the added-up power of these two processes results in a heating power \(P_Cond\) which increases the temperature \(T_{Cond,in}\) of an input medium to an output temperature of \(T_{Cond,out}\).

• In heating mode, the evaporator is connected to a source medium (air, water, …) and releases the condenser energy to a radiator or similar device.

• In cooling mode, the heating intake on the evaporator side is used to operate a chiller and the condenser side takes the cold medium from a source.

Since ehubX works based on energy carriers, we need to assign one energy carrier from \(\mathcal{S}_{Ec}\) to each power node of this schematic. These are specified in techs.yaml and have the following interpretation:

Energy carrier of heat pump technology \(x\)

Mode	Power node	Energy carrier from \(\mathcal{S}_{Ec}\)
Heating	Electricity consumption	\(ec\_el[x]\)
Heating	Evaporator power	\(ec\_ht\_in[x]\)
Heating	Condenser power	\(ec\_ht\_out[x]\)
Cooling	Electricity consumption	\(ec\_el[x]\)
Cooling	Evaporator power	\(ec\_co\_out[x]\)
Cooling	Condenser power	\(ec\_co\_in[x]\)

As can be seen, we label the energy carriers from the perspective of the energy usage. It is also essential that the five energy carriers do not contain any duplicates since this would break the heat pump model below, leading to potential infeasibilities or unbounded solutions. Checks are included in the data validation phase to ensure that this is not the case. Using the above ecs, we define sets of input/output tuples for the heat pump:

\begin{align*} \mathcal{S}_{HpTechIn} &= \Big \{ (s, h, x, e): (s, h, x) \in \mathcal{S}_{HpTechTuple}, ~e \in \{ec\_el[x], ec\_ht\_in[x], ec\_co\_in[x] \} \Big \} \\ \mathcal{S}_{HpTechOut} &= \Big \{ (s, h, x, e): (s, h, x) \in \mathcal{S}_{HpTechTuple}, ~e \in \{ec\_ht\_out[x], ec\_co\_out[x] \} \Big \} \end{align*}

Based on these sets, we can now define heat pump input and output variables as follows:

\begin{align*} \mathcal{V}_{HpTechIn}: \mathcal{S}_{HpTechIn} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \\ \mathcal{V}_{HpTechOut}: \mathcal{S}_{HpTechOut} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \end{align*}

Note that this formulation so far suggests that heating and cooling modes may occur simultaneously, e.g.; if \(\mathcal{V}_{HpTechOut}[s, h, x, ec\_ht\_out, t] > 0\) and \(\mathcal{V}_{HpTechOut}[s, h, x, ec\_co\_out, t] > 0\). Indeed, we allow for this to be the case because ehubX never considers a single asset but a collection of them. As such, the heat pump capacity in a single is allowed to consist of multiple real-life heat pumps where some operate in heating and others in cooling mode. However, we do need to respect the distinction between different modes in the model itself, which is why we define two additional power consumption variable for heating and cooling mode:

\begin{align*} \mathcal{V}_{HpTechElecHt}, \mathcal{V}_{HpTechElecCo}: \mathcal{S}_{HpTechTuple} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \end{align*}

The sum of these two variables needs to be identical to the total power consumption which is handled in the following constraint:

\begin{align*} \mathcal{V}_{HpTechIn}[s, h, x, ec\_el[x], t] = \mathcal{V}_{HpTechElecHt}[s, h, x, t] + \mathcal{V}_{HpTechElecCo}[s, h, x, t] \end{align*}

The heat pump needs to satisfy a power balance equation. More specifically, the condenser power is comprised of the sum of compressor consumption and evaporator power. Formulating this for both heating and cooling modes reads:

\begin{align*} \mathcal{V}_{HpTechOut}[s, h, x, ec\_ht\_out[x], t] &= \mathcal{V}_{HpTechElecHt}[s, h, x, t] + \mathcal{V}_{HpTechIn}[s, h, x, ec\_ht\_in[x], t] \\ \mathcal{V}_{HpTechIn}[s, h, x, ec\_co\_in[x], t] &= \mathcal{V}_{HpTechElecCo}[s, h, x, t] + \mathcal{V}_{HpTechOut}[s, h, x, ec\_co\_out[x], t]\end{align*}

An essential indicator for heat pumps is the Coefficient of Performance (COP). In ehubX, this is considered as a potentially stage-, hub- and time-dependent parameter \(cop[s, h, x, t]\) that acts as the quotient between condenser power and electricity consumption. Translated into heating and cooling modes, this yields the constraint

\begin{align*} \mathcal{V}_{HpTechOut}[s, h, x, ec\_ht\_out[x], t] &= cop[s, h, x, t] \cdot \mathcal{V}_{HpTechElecHt}[s, h, x, t] \\ \mathcal{V}_{HpTechIn}[s, h, x, ec\_co\_in[x], t] &= cop[s, h, x, t] \cdot \mathcal{V}_{HpTechElecCo}[s, h, x, t] \end{align*}

As an input, \(cop\) can either be specified directly in hubs.yaml. If it is not given as a default or profile value, ehubX will try to calculate it from the parameters \(temp\_heat\_in\), \(temp\_heat\_out\) and \(cop\_factor\) (also from hubs.yaml) using a dampened Carnot efficiency model as follows:

\begin{align*} cop[s, h, x, t] = cop\_factor[s, x] \cdot \frac{temp\_heat\_out[s, h, x, t]}{temp\_heat\_out[s, h, x, t] - temp\_heat\_in[s, h, x, t]} \end{align*}

The heat pump module still needs to be linked to two quantities from the tech model, namely installed capacity and technology usage. For the capacity, it is common to measure heat pump capacity in terms of condenser power. For this model, this results in the following constraint, dampened by the factor availability from hubs.yaml:

\begin{align*} \mathcal{V}_{HpTechOut}[s, h, x, ec\_ht\_out, t] + \mathcal{V}_{HpTechIn}[s, h, x, ec\_co\_in, t] \le availability[s, h, x, t] \cdot \mathcal{V}_{TechCap}[s, h, x] \end{align*}

Finally, detecting heat pump usage is also formulated using the condenser power:

\begin{align*} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \big( \mathcal{V}_{HpTechOut}[s, h, x, ec\_ht\_out, t] + \mathcal{V}_{HpTechIn}[s, h, x, ec\_co\_in, t] \big) \le BigM_{HpCondSum}[s, h, x] \cdot \mathcal{V}_{YTechUsed}[s, h, x] \end{align*}

Here, \(weight\) is a clustering value and \(BigM_{HpCondSum}[s, h, x]\) is a big-M parameter which has to provide an upper bound to the summed-up power of the heat pump condenser over the time horizon. It is chosen as

\begin{align*} BigM_{HpTechCap}[s, h, x] = \left \{ \begin{array}{rl} cap\_max[s, h, x] \cdot |\mathcal{S}_{TimeHorizon}| + 10^{-5}, &\text{if } cap\_max[s, h, x] < \infty \\ \mathcal{P}_{BigMGeneric}, &\text{else} \end{array} \right . \end{align*}

We recall that \({P}_{BigMGeneric}\) is the demand-based default big-M parameter introduced in the demand model. Furthermore, \(cap\_max\) is a parameter from hubs.yaml.

Network model

The network model makes it possible to transfer energy between different hubs of the model. In network_links.yaml, a set of link is defined running between hubs. We gather these in a set

\begin{align*} \mathcal{S}_{NetLink} = \big \{ \ell: \ell \text{ is a network link } \big \} \end{align*}

In addition, network_techs.yaml defines a set of network technologies that can be installed and used for the actual transfer, here collected in the set

\begin{align*} \mathcal{S}_{NetTech} = \big \{ n: n \text{ is a network technology } \big \} \end{align*}

network_links.yaml additionally defines a connectivity structure between hubs and ecs defined on the links. This structure can be expressed for all inputs to the network by the following set:

\begin{align*} \mathcal{S}_{NetLinkIn} = &\big \{ (h, \ell, e): \ell \in \mathcal{S}_{NetLink}, ~ e \text{ is ec of } \ell, ~ h \text{ is } start\_hub \text{ of } \ell \big \} \\ &\cup \big \{ (h, \ell, e): \ell \in \mathcal{S}_{NetLink}, ~ \ell \text{ is bidirectional}, ~e \text{ is ec of } \ell, ~ h \text{ is } end\_hub \text{ of } \ell \big \} \end{align*}

This clarifies also what we understand by bidirectionality: Every link is able to transfer an ec from its \(start\_hub\) to its \(end\_hub\) but only bidirectional links can also transfer energy the other way. The same logic allows to define all outputs of the network as

\begin{align*} \mathcal{S}_{NetLinkOut} = &\big \{ (h, \ell, e): \ell \in \mathcal{S}_{NetLink}, ~ e \text{ is ec of } \ell, ~ h \text{ is } end\_hub \text{ of } \ell \big \} \\ &\cup \big \{ (h, \ell, e): \ell \in \mathcal{S}_{NetLink}, ~ \ell \text{ is bidirectional}, ~e \text{ is ec of } \ell, ~ h \text{ is } start\_hub \text{ of } \ell \big \} \end{align*}

On these two sets, we have variables quantifying the amount of inputs and outputs to the network transfer system:

\begin{align*} \mathcal{V}_{NetLinkIn} &: \mathcal{S}_{Stage} \times \mathcal{S}_{NetLinkIn} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \\ \mathcal{V}_{NetLinkOut} &: \mathcal{S}_{Stage} \times \mathcal{S}_{NetLinkOut} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \end{align*}

From the perspective of the rest of the energy system, links are not really an aspect that is considered. So it stands to reason that for an energy balance perspective, the only relevant quantities are the inputs and outputs at each hub, aggregated over all possible links. We prepare this procedure by defining the hub inputs/output sets by

\begin{align*} \mathcal{S}_{NetHubIn} &= \big \{ (h, e) \in \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}: (h, \ell, e) \in \mathcal{S}_{NetLinkIn} \text{ for any } \ell \in \mathcal{S}_{NetLink} \big \} \\ \mathcal{S}_{NetHubOut} &= \big \{ (h, e) \in \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}: (h, \ell, e) \in \mathcal{S}_{NetLinkOut} \text{ for any } \ell \in \mathcal{S}_{NetLink} \big \} \end{align*}

Accordingly, we define the hub inputs/outputs on these sets:

\begin{align*} \mathcal{V}_{NetHubIn} &: \mathcal{S}_{Stage} \times \mathcal{S}_{NetHubIn} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \\ \mathcal{V}_{NetHubOut} &: \mathcal{S}_{Stage} \times \mathcal{S}_{NetHubOut} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \end{align*}

The connection between link inputs/outputs and hub inputs/outputs is naturally given by

\begin{align*} \mathcal{V}_{NetHubIn}[s, h, e, t] &= \sum\limits_{\substack{\ell \in \mathcal{S}_{NetLink} \\ (h, \ell, e) \in \mathcal{S}_{NetLinkIn}}} \mathcal{V}_{NetLinkIn}[s, h, e, t] \\ \mathcal{V}_{NetHubOut}[s, h, e, t] &= \sum\limits_{\substack{\ell \in \mathcal{S}_{NetLink} \\ (h, \ell, e) \in \mathcal{S}_{NetLinkOut}}} \mathcal{V}_{NetLinkOut}[s, h, e, t] \end{align*}

In a next step, we consider that network transfer is possible through network technologies. Similar to the technology model, these technologies may only be eligible for certain stages (due to its technology readiness level, see network_techs.yaml) or certain links (due to allowed_net_tech_lists, see network_links.yaml). The stage-link-tech tuples of all allowed combinations is captured in the set \(\mathcal{S}_{NetTechTuple} \subset \mathcal{S}_{Stage} \times \mathcal{S}_{NetLink} \times \mathcal{S}_{NetTech}\). Based on these set, we can define the set of all possible network technologies that could contribute to a network transfer along a specific link in the sets

\begin{align*} \mathcal{S}_{NetTechIn} &= \big \{ (s, h, \ell, n) \in \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{NetLink} \times \mathcal{S}_{NetTech}: (s, \ell, n) \in \mathcal{S}_{NetTechTuple}, (h, \ell, ec[n]) \in \mathcal{S}_{NetLinkIn} \big \} \\ \mathcal{S}_{NetTechOut} &= \big \{ (s, h, \ell, n) \in \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{NetLink} \times \mathcal{S}_{NetTech}: (s, \ell, n) \in \mathcal{S}_{NetTechTuple}, (h, \ell, ec[n]) \in \mathcal{S}_{NetLinkOut} \big \} \end{align*}

Here we have used the parameter \(ec[n]\) from network_techs.yaml which maps each network technology to its unique ec. As for the other sets, we assign the amount of network link transfer handled by each network technology to the variables

\begin{align*} \mathcal{V}_{NetTechIn}&: \mathcal{S}_{NetTechIn} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \\ \mathcal{V}_{NetTechOut}&: \mathcal{S}_{NetTechOut} \times \mathcal{S}_{Time} \to \mathbb{R}_0^+ \end{align*}

Naturally, the link inputs/outputs are then simply given be the sum over all network technologies:

\begin{align*} \mathcal{V}_{NetLinkIn}[s, h, \ell, e, t] &= \sum\limits_{\substack{n \in \mathcal{S}_{NetTech} \\ (s, h, \ell, n) \in \mathcal{S}_{NetTechIn} \\ e = ec[n]}} \mathcal{V}_{NetTechIn}[s, h, \ell, n, t] \\ \mathcal{V}_{NetLinkOut}[s, h, \ell, e, t] &= \sum\limits_{\substack{n \in \mathcal{S}_{NetTech} \\ (s, h, \ell, n) \in \mathcal{S}_{NetTechOut} \\ e = ec[n]}} \mathcal{V}_{NetTechOut}[s, h, \ell, n, t] \end{align*}

This concludes the long establishment process for all variables. To summarize, we have defined three layers:

1. Hub inputs/outputs, set \(\mathcal{S}_{NetHubIn/Out} \subset \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}\).

2. Link inputs/outputs, set \(\mathcal{S}_{NetLinkIn/Out} \subset \mathcal{S}_{Hub} \times \mathcal{S}_{NetLink} \times \mathcal{S}_{Ec}\)

3. Tech inputs/outputs, set \(\mathcal{S}_{NetTechIn/Out} \subset \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{NetLink} \times \mathcal{S}_{NetTech}\)

where each layer is defined as the aggregated sum of the layer below it.

In a next step, we connect the input and output amounts for the network technologies by the following constraint:

\begin{align*} \mathcal{V}_{NetTechOut}[s, h_{out}, \ell, n, t] = exp \big( -trans\_decay[s,n] \cdot length[\ell] \big) \cdot \mathcal{V}_{NetTechIn}[s, h_{in}[\ell], \ell, n, t] \end{align*}

where \(trans\_decay\) is the transmission decay parameter from network_techs.yaml, \(length\) is the link length from network_links.yaml and \(h_{in}\) is given by

\begin{align*} h_{in}[\ell] = \left \{ \begin{array}{rl} hub\_start[\ell], &\text{if } h_{out} = hub\_end[\ell] \\ hub\_end[\ell], &\text{if } h_{out} = hub\_start[\ell] \end{array} \right . \end{align*}

In network_links.yaml, the parameters \(trans\_sum\_min\_forward\), \(trans\_sum\_min\_backward\), \(trans\_sum\_max\_forward\) and \(trans\_sum\_max\_backward\) denote maximal and minimal transmission amounts for an ec along each link in forward and backward directions. They result in the constraints

\begin{align*} trans\_sum\_min[s, h_{out}, \ell, e] ~\le~ \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{NetLinkOut}[s, h_{out}, \ell, e, t] ~\le~ trans\_sum\_max[s, h_{out}, \ell, e] \end{align*}

where we have used the parameters

\begin{align*} trans\_sum\_min[s, h_{out}, \ell, e] &= \left \{ \begin{array}{rl} trans\_sum\_min\_forward[s, \ell, e], &\text{if } h_{out} = hub\_end[\ell] \\ trans\_sum\_min\_backard[s, \ell, e], &\text{if } h_{out} = hub\_start[\ell] \end{array} \right . \\ trans\_sum\_max[s, h_{out}, \ell, e] &= \left \{ \begin{array}{rl} trans\_sum\_max\_forward[s, \ell, e], &\text{if } h_{out} = hub\_end[\ell] \\ trans\_sum\_max\_backard[s, \ell, e], &\text{if } h_{out} = hub\_start[\ell] \end{array} \right . \end{align*}

As it was done for the technology model, ehubX needs to be able to have an understanding of capacity for network technologies, and be able to make its own decisions about the optimal capacity along eacah link. We prepare this by introducing the variable

\begin{align*} \mathcal{V}_{NetTechCap}: \mathcal{S}_{Stage} \times \mathcal{S}_{NetLink} \times \mathcal{S}_{NetTech} \to \mathbb{R} \end{align*}

These capacities are to be understood as the maximal amount of input power that the technologies can handle. In contrast to the technology model, we do not specify minimal or maximal capacities for each technology, but rather each ec along a link. To formulate this mathematically, we require the set

\begin{align*} \mathcal{S}_{NetLinkAndEc} = \big \{ (\ell, e) \in \mathcal{S}_{NetLink} \times \mathcal{S}_{Ec}: e \in ecs[\ell] \big \} \end{align*}

The capacity restrictions then take the following form for all \(s \in \mathcal{S}_{Stage}\) and all \((\ell, e) \in \mathcal{S}_{NetLinkAndEc}\):

\begin{align*} cap\_min[s, \ell, e] \le \sum\limits_{\substack{n \in \mathcal{S}_{NetTech} \\ e = ec[n]}} \mathcal{V}_{NetTechCap}[s, \ell, n] \le cap\_max[s, \ell, e] \end{align*}

The total capacity of a network technology along a link is composed of both the remnants of the initial capacity (\(cap\_init\) in network_links.yaml) and the installed capacity from previous stages, which ehubX can decide upon via the variable

\begin{align*} \mathcal{V}_{NetTechCapInstl}: \mathcal{S}_{Stage} \times \mathcal{S}_{NetLink} \times \mathcal{S}_{NetTech} \to \mathbb{R}_0^+ \end{align*}

To formulate the sum of network technology capacity, we add the constraint

\begin{align*} \mathcal{V}_{NetTechCap}[s, \ell, n] = Cap_{init}[s, \ell, n] + \sum\limits_{s_{instl} \in \mathcal{S}_{Stage}} Cap_{instl}[s_{instl}, s, \ell, n] \end{align*}

The initial capacity that is still operational in stage \(s\) is given by

\begin{align*} Cap_{init}[s, \ell, n] = \left \{ \begin{array}{rl} cap\_init[\ell, n], &\text{if } start\_year[s] - init\_year < lifetime[n] - age\_init[\ell, n] \\ 0, &\text{else} \end{array} \right . \end{align*}

Here, \(cap\_init\) and \(age\_init\) are initial capacity parameters from network_links.yaml, \(start\_year\) is the first year of a stage from stages.yaml, and \(init\_year\) is the start year of the first stage (earliest \(start\_year\) in stages.yaml). The capacity installed in stage \(s_{instl}\) that is still operational in stage \(s\) is given by

\begin{align*} Cap_{instl}[s_{instl}, s, \ell, n] = \left \{ \begin{array}{rl} \mathcal{V}_{NetTechCapInstl}[s_{instl}, \ell, n], &\text{if } 0 \le start\_year[s] - start\_year[s_{instl}] < lifetime[n] \\ 0, &\text{else} \end{array} \right . \end{align*}

\(lifetime\) is the lifetime of a network technology parametrized in network_techs.yaml.

Now that network technology capacity can be installed by the optimizer, we need to limit the actual network inputs by this capacity, as follows:

\begin{align*} \sum\limits_{\substack{h \in \mathcal{S}_{Hub} \\ (s, h, \ell, n) \in \mathcal{S}_{NetTechIn}}} \mathcal{V}_{NetTechIn}[s, h, \ell, n, t] \le \mathcal{V}_{NetTechCap}[s, \ell, n] \end{align*}

The sum over the hubs in this formulation only accounts for the fact that the link \(\ell\) might have a single input hub (just \(hub\_start[\ell]\) if \(\ell\) is unidirectional) or two input hubs (\(hub\_start[\ell]\) and \(hub\_end[\ell]\) if \(\ell\) is bidirectional).

As a next step, we need to respect the parameter \(unit\_cap\_min\) from network_techs.yaml which declares the minimal amount of capacity that has to be installed if any installation takes place at all. To implement this, we first require a binary variable monitoring whether any installation takes place at all, given by

\begin{align*} \mathcal{V}_{YNetTechCapInstl}: \mathcal{S}_{Stage} \times \mathcal{S}_{NetLink} \times \mathcal{S}_{NetTech} \to \{0, 1 \} \end{align*}

This variable is forced to activate if \(\mathcal{V}_{NetTechCapInstl}\) is nonnegative by the big-M constraint

\begin{align*} \mathcal{V}_{NetTechCapInstl[s, \ell, n]} \le BigM_{NetTechCap}[s, \ell, n] \cdot \mathcal{V}_{YNetTechCapInstl}[s, \ell, n] \end{align*}

where the big-M parameter can be parametrized with \(cap\_max\) from network_links.yaml and the generic big-M value \(\mathcal{P}_{BigMGeneric}\) from the demand model as

\begin{align*} BigM_{NetTechCap}[s, \ell, n] = \left \{ \begin{array}{rl} 10^{-5} + cap\_max[s, \ell, ec[n]], &\text{if } cap\_max[s, \ell, ec[n]] < \infty \\ \mathcal{P}_{BigMGeneric}, &\text{else} \end{array} \right . \end{align*}

Now that we have the binary variable \(\mathcal{V}_{YNetTechCapInstl}\), we can include the functionality of \(unit\_cap\_min\) from network_techs.yaml by demanding that

\begin{align*} \mathcal{V}_{NetTechCapInstl}[s, \ell, n] \ge unit\_cap\_min[s, n] \cdot \mathcal{V}_{YNetTechCapInstl}[s, \ell, n] \end{align*}

We also need to track if a network technology is used at all on a specific stage-link tuple which is monitored by a binary variable

\begin{align*} \mathcal{V}_{YNetTechUsed}: \mathcal{S}_{NetTechTuple} \to \{0, 1 \} \end{align*}

This relationship is enforced by the a big-M constraint

\begin{align*} \sum\limits_{\substack{h \in \mathcal{S}_{Hub} \\ (s, h, \ell, n) \in \mathcal{S}_{NetTechIn}}} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{NetTechIn}[s, h, \ell, n, t] \le BigM_{NetTechLinkIn}[s, \ell, n] \cdot \mathcal{V}_{YNetTechUsed}[s, \ell, n] \end{align*}

Here, the big-M parameter \(BigM_{NetTechLinkIn}[s, \ell, n]\) is a tight upper limit to the maximal amount of network input that is possible, given by

\begin{align*} BigM_{NetTechLinkIn}[s, \ell, n] = \left \{ \begin{array}{rl} 10^{-5} + cap\_max[s, \ell, e] \cdot |\mathcal{S}_{TimeHorizon}|, &\text{if } cap\_max[s, \ell, e] < \infty \\ \mathcal{P}_{BigMGeneric}, &\text{else} \end{array} \right .\end{align*}

\(cap\_max\) is the maximal capacity parameter from network_links.yaml and \(\mathcal{P}_{BigMGeneric}\) is the demand-based default big-M parameter introduced in the demand model.

At this point, we can consider the costs that arise from network technology installation and maintenance. Starting with CAPEX costs, we introduce the fixed variable

\begin{align*} \mathcal{V}_{NetTechCostCapex}: \mathcal{S}_{NetTechTuple} &\to \mathbb{R}, \\ \mathcal{V}_{NetTechCostCapex}[s, \ell, n] &= length[\ell] \cdot CRF(interest\_rate[n], lifetime[n]) \cdot \sum\limits_{\substack{s_{instl} \in \mathcal{S}_{Stage} \\ (s_{instl}, \ell, n) \in \mathcal{S}_{NetTechTuple}}} \cdot Capex_{Stage}[s_{instl}, s, \ell, n] \end{align*}

where we have used the same formulation as in the technology model but added the link length parameter \(length\) from network_links.yaml. \(CRF\) is a standard capital recovery factor calculated from the technology’s lifetime and interest rate as follows:

\begin{align*} CRF(i, N) = \frac{i \cdot (i+1)^N}{(i+1)^{N-1}} \end{align*}

The summand takes the following form:

\begin{align*} Capex_{Stage}[s_{instl}, s, \ell, n] ~=~ \left \{ \begin{array}{rl} capex\_per\_cap[s, n] \cdot \mathcal{V}_{NetTechCapInstl}[s_{instl}, \ell, n] & \\ +~ one\_time\_capex[s, n] \cdot \mathcal{V}_{YNetTechCapInstl}[s_{instl}, \ell, n], &\text{if } 0 \le start\_year[s] - start\_year[s_{instl}] < lifetime[n] \\ 0, &\text{else} \end{array} \right . \end{align*}

Here, \(capex\_per\_cap\) and \(one\_time\_capex\) are parameters from network_techs.yaml.

In addition to these installation-related CAPEX costs, operating and maintaining the network technologies is associated with the OPEX costs parametrized by \(opex\_per\_cap\) and \(one\_time\_opex\) from network_techs.yaml:

\begin{align*} \mathcal{V}_{NetTechCostOpexCap}: \mathcal{S}_{NetTechTuple} &\to \mathbb{R}, \\ \mathcal{V}_{NetTechCostOpexCap}[s, \ell, n] &= length[ell] \cdot \big( opex\_per\_cap[s, n] \cdot \mathcal{V}_{NetTechCap}[s, \ell, n] + one\_time\_opex[s, n] \cdot \mathcal{V}_{YNetTechUsed}[s, \ell, n] \big) \end{align*}

As can be seen, both costs have a part that arises per amount of installed capacity while another enters any time any amount of network technology is installed or used at all. On top of these capacity-related costs, the network model contains a cost factor that arises from the parameter \(opex\_per\_energy\) in network_techs.yaml, given for the amount of transmitted energy by the variable

\begin{align*} \mathcal{V}_{NetTechCostOpexTrans}: \mathcal{S}_{NetTechTuple} &\to \mathbb{R}, \\ \mathcal{V}_{NetTechCostOpexTrans}[s, \ell, n] &= length[\ell] \cdot opex\_per\_energy[s, n] \cdot \sum\limits_{\substack{h \in \mathcal{S}_{Hub} \\ (s, h, \ell, n) \in \mathcal{S}_{NetTechIn}}} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{NetTechIn}[s, h, \ell, n, t] \end{align*}

Adding up these costs gives us the following fixed variable:

\begin{align*} \mathcal{V}_{NetTechCostTotal} &\in \mathbb{R}, \\ \mathcal{V}_{NetTechCostTotal} &= \sum\limits_{(s, \ell, n) \in \mathcal{S}_{NetTechTuple}} \big( \mathcal{V}_{NetTechCostCapex}[s, \ell, n] + \mathcal{V}_{NetTechCostOpexCap}[s, \ell, n] + \mathcal{V}_{NetTechCostOpexTrans}[s, \ell, n] \big) \end{align*}

Similar to the costs associated with network technologies, certain CO2-embodied emissions are associated with this process as well. We collect these in the variable

\begin{align*} \mathcal{V}_{NetTechCo2Instl}: \mathcal{S}_{NetTechTuple} &\to \mathbb{R} \\ \mathcal{V}_{NetTechCo2Instl}[s, \ell, n] &= length[\ell] \cdot \sum\limits_{\substack{s_{instl} \in \mathcal{S}_{Stage} \\ (s_{instl}, \ell, n) \in \mathcal{S}_{NetTechTuple} \\ 0 \le start\_year[s] - start\_year[s_{instl}] \le lifetime[n]}} \frac{co2\_per\_cap[s, n]}{lifetime[n]} \cdot \mathcal{V}_{NetTechCapInstl}[s_{instl}, \ell, n] \end{align*}

The parameter \(co\_per\_cap\) comes from the network_techs.yaml file.

In addition to installation-related CO2 emissions, they may arise relative to the total amount of transmitted energy, as parametrized by \(co2\_per\_energy\) from network_techs.yaml and tracked by the fixed variable

\begin{align*} \mathcal{V}_{NetTechCo2Trans}: \mathcal{S}_{NetTechTuple} &\to \mathbb{R} \\ \mathcal{V}_{NetTechCo2Trans}[s, \ell, n] &= length[\ell] \cdot co2\_per\_energy[s, n] \cdot \sum\limits_{\substack{h \in \mathcal{S}_{Hub} \\ (s, h, \ell, n) \in \mathcal{S}_{NetTechIn}}} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, n] \cdot \mathcal{V}_{NetTechIn}[s, h, \ell, n, t] \end{align*}

As with the costs, we gather the total amount of CO2 emissions into a fixed variable for convenience:

\begin{align*} \mathcal{V}_{NetTechCo2Total}: \mathcal{S} \to \mathbb{R}, \qquad \mathcal{V}_{NetTechCo2Total}[s] = \sum\limits_{\substack{(s', \ell, n) \in \mathcal{S}_{NetTechTuple} \\ s = s'}} \big( \mathcal{V}_{NetTechCo2Instl}[s', \ell, n] + \mathcal{V}_{NetTechCo2Trans}[s', \ell, n] \big) \end{align*}

Self-sufficiency model

In ehubX, the concept of self-sufficiency is calculated by two kinds of energy imports. The first is called cross-border imports, measured in a fixed scalar variable

\begin{align*} \mathcal{V}_{SelfSufficiencyImpCross} &\in \mathbb{R}_0^+, \\ \mathcal{V}_{SelfSufficiencyImpCross} &= \sum\limits_{\substack{(s, h, e) \in \mathcal{S}_{ImpTuple} \\ is\_energy[e] = True \\ imp\_exp\_type[e] = cross}} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{Imp}[s, h, e, t] \end{align*}

The ec-specific parameters \(is\_energy\) and \(imp\_exp\_type\) can be specified in ecs.yaml, and \(weight\) is a clustering parameter. The idea behind this variable is that it should measure all imports into the system that would be unavailable if neighboring energy systems were not able to provide them. We think of energy carriers like electricity or gas that are purchased from different countries or regions.

In contrast the these cross-import, the ehubX model also allows for “imports” that originate in natural resources and would remain be available if the energy system became completely isolated. These are e.g.; solar energy or thermal resources. The total amount of these imports is measured in the fixed scalar variable

\begin{align*} \mathcal{V}_{SelfSufficiencyImpInternal} &\in \mathbb{R}_0^+, \\ \mathcal{V}_{SelfSufficiencyImpInternal} &= \sum\limits_{\substack{(s, h, e) \in \mathcal{S}_{ImpTuple} \\ is\_energy[e] = True \\ imp\_exp\_type[e] = internal}} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{Imp}[s, h, e, t] + \sum\limits_{\substack{(s, h, x, e) \in \mathcal{S}_{ConvTechOut} \\ x \in \mathcal{S}_{ConvTechInternal}}} \sum\limits_{t \in \mathcal{S}_{Time}} weight[s, t] \cdot \mathcal{V}_{ConvTechOut}[s, h, x, e, t] \end{align*}

The first summand in the above equation is a similar type of value as in the definition of \(\mathcal{V}_{SelfSufficiencyImpCross}\). The second term considers the internal import of ecs that are not directly measured in terms of energy, like solar irradiation. These ecs are characterized by the fact that they have \(is\_energy[e] = False\) and \(imp\_exp\_type[e] = internal\). In order to accurately include these carriers into the internal imports, we have considered a special subset of conversion technologies, given by

\begin{align*} \mathcal{S}_{ConvTechInternal} = \big \{ x \in \mathcal{S}_{ConvTech}: ~ x \text{ has } &\text{a single input } ec\_in \text{ with } is\_energy[ec\_in] = False \text{ and } imp\_exp\_type[ec\_in] = internal \\ \text{and } &\text{a single output } ec\_out \text{ with } is\_energy[ec\_out] = True \big \} \end{align*}

These technologies transform a single non-energy internal input ec into a single energy output ec, and the outputs of the conversion technology are included in the internal imports.

Having calculated both internal imports and cross-border imports, the self-sufficiency measure is defined as

\begin{align*} SelfSufficiency = \frac{\mathcal{V}_{SelfSufficiencyImpInternal}}{\mathcal{V}_{SelfSufficiencyImpInternal} + \mathcal{V}_{SelfSufficiencyImpCross}} \end{align*}

The actual variable we use to track the self-sufficiency measure is given by

\begin{align*} \mathcal{V}_{SelfSufficiency} \in \mathbb{R}_0^+ \end{align*}

There are now two options, based on the self-sufficiency calculation method selected in stages.yaml. First, if the method ‘quadratic’ is chosen, this nonlinear constraint is directly included in the model, making it a quadratically constrained system:

\begin{align*} \mathcal{V}_{SelfSufficiency} \cdot \big(\mathcal{V}_{SelfSufficiencyImpInternal} + \mathcal{V}_{SelfSufficiencyImpCross}) = \mathcal{V}_{SelfSufficiencyImpInternal} \end{align*}

Second, if the method ‘linearized’ is chosen, the original quadratic constraint relating \(\mathcal{V}_{SelfSufficiency}\) to \(\mathcal{V}_{SelfSufficiencyImpInternal}\) and \(\mathcal{V}_{SelfSufficiencyImpCross}\) is replaced by a linearization to ensure compatibility with MILP solvers. ehubX will do this by computing a boxed domain

\begin{align*} \Omega &= [\mathcal{V}_{\mathrm{SelfSufficiencyImpInternal, min}}, \mathcal{V}_{\mathrm{SelfSufficiencyImpInternal, max}}] \times [\mathcal{V}_{\mathrm{SelfSufficiencyImpCross, min}}, \mathcal{V}_{\mathrm{SelfSufficiencyImpCross, max}}] \\[8pt] \text{where} \quad &0 \leq \mathcal{V}_{\mathrm{SelfSufficiencyImpInternal, min}} < \mathcal{V}_{\mathrm{SelfSufficiencyImpInternal, max}}, \quad 0 \leq \mathcal{V}_{\mathrm{SelfSufficiencyImpCross, min}} < \mathcal{V}_{\mathrm{SelfSufficiencyImpCross, max}}, \quad 0 \notin \Omega. \end{align*}

Furthermore, \(\mathcal{V}_{SelfSufficiencyImpInternal}\) and \(\mathcal{V}_{SelfSufficiencyImpCross}\) represent input variables. To facilitate the linearization process, the rectangular domain is segmented into simplexes (triangles). A grid is imposed on each dimension as follows:

\begin{align*} \mathcal{V}_{\mathrm{SelfSufficiencyImpInternal, min}} &= \hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpInternal, 1}} < \dots < \hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpInternal, n}} = \mathcal{V}_{\mathrm{SelfSufficiencyImpInternal, max}}, \quad n \in \mathbb{N}. \end{align*}

\begin{align*} \mathcal{V}_{\mathrm{SelfSufficiencyImpCross, min}} &= \hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpCross, 1}} < \dots < \hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpCross, m}} = \mathcal{V}_{\mathrm{SelfSufficiencyImpCross, max}}, \quad m \in \mathbb{N}. \end{align*}

This results in \(n \cdot m\) grid points, labeled as:

\begin{align*} \hat{x}_{ij} = (\hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpInternal, i}}, \hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpCross, j}})^\top, \quad i = 1,\ldots,n,\quad j = 1,\ldots,m. \end{align*}

A key aspect of the formulation is that any point \(x = (\mathcal{V}_{SelfSufficiencyImpInternal},\mathcal{V}_{SelfSufficiencyImpCross}) \in \Omega\) can be expressed as a convex combination of these grid points:

\begin{align*} x &= \sum_{i=1}^n \sum_{j=1}^m \lambda_{ij}\,\hat{x}_{ij}, \\[8pt] 1 &= \sum_{i=1}^n \sum_{j=1}^m \lambda_{ij}, \\[8pt] \lambda_{ij} &\geq 0, \quad \forall\,i = 1,\dots,n, \quad \forall\,j = 1,\dots,m. \end{align*}

Alternatively, splitting by dimension gives the equivalent formulation:

\begin{align*} \mathcal{V}_{SelfSufficiencyImpInternal} &= \sum_{i=1}^n \sum_{j=1}^m \lambda_{ij}\,\hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpInternal, i}}, \\[4pt] \mathcal{V}_{SelfSufficiencyImpCross} &= \sum_{i=1}^n \sum_{j=1}^m \lambda_{ij}\,\hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpCross, j}}. \end{align*}

The representation above for \(x\) is not unique, thus it is important that \(x\) is assigned to exactly one triangle. Initially, there are \(L := 2 \cdot (n-1) \cdot (m-1)\) triangles in total. A binary variable is introduced as

\begin{align*} z_\ell \in \{0,1\}, \quad \forall\,\ell = 1,\ldots,L. \end{align*}

The goal is to arrive at an equation system where \(x\) is assigned to exactly one simplex \(\ell\), i.e., \(z_\ell = 1\) for one \(\ell\) and zero otherwise:

\begin{align*} \sum_{\ell=1}^L z_\ell = 1. \end{align*}

Lastly \(x\) must be connected to the correct triangle. In the expression for \(x\) above, the convex multipliers \(\lambda_{ij}\) were used to express \(x\) via the grid points \(\hat{x}_{ij}\). Therefore, what remains is to establish the connection between these multipliers and the triangle indicator variables \(z_\ell\). For each triangle \(1 \le \ell \le L\), let \(\mathcal{E}_\ell\) denote the corner points of this triangle. An additional restriction is then imposed to ensure that only the active triangle’s grid points may contribute:

\begin{align*} \lambda_{ij} \;\le\; \sum_{\{\ell=1,\dots,L \mid \hat{x}_{ij}\,\in\,\mathcal{E}_\ell\}} z_\ell. \end{align*}

Once the indexing of the triangles and the corresponding assignment of grid points (via the corner sets) is set, the constraints above are fully specified.

With the constraints in place, the linearized value for \(\mathcal{V}_{\mathrm{SelfSufficiency}}\bigl(\mathcal{V}_{\mathrm{SelfSufficiencyImpInternal}}, \mathcal{V}_{\mathrm{SelfSufficiencyImpCross}}\bigr) = \mathcal{V}_{\mathrm{SelfSufficiency}}(x)\) is given by

\begin{align*} \mathcal{V}_{\mathrm{SelfSufficiency}} = \sum_{i=1}^n \sum_{j=1}^m \lambda_{ij}\,\hat{\mathcal{V}}_{\mathrm{SelfSufficiency, i, j}}, \end{align*}

where the precomputed values of \(\mathcal{V}_{SelfSufficiency}\) at the grid points have been inserted

\begin{align*} \hat{\mathcal{V}}_{\mathrm{SelfSufficiency, i, j}} &= SelfSufficiency\bigl(\hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpInternal, i, j}},\hat{\mathcal{V}}_{\mathrm{SelfSufficiencyImpCross, i, j}}\bigr), \quad \forall\,i = 1,\ldots,n,\quad \forall\,j = 1,\ldots,m. \end{align*}

Finally, independent of the method that is used for the calculation of the variable \(\mathcal{V}_{SelfSufficiency}\), two additional constraints based on the parameters self_sufficiency_min and self_sufficiency_max from stages.yaml are added:

\begin{align*} self\_sufficiency\_min ~\le~ \mathcal{V}_{SelfSufficiency} ~\le~ self\_sufficiency\_max \end{align*}

Autonomy model

In ehubX, the concept of autonomy measures how long an energy system can operate without relying on cross-border imports. The autonomy is defined as the length of the longest consecutive time period during which demand can be satisfied using only internal resources.

The core of the formulation is a binary variable indicating whether the system is still autonomous (“alive”) at a given stage and time step:

\begin{align*} \mathcal{V}_{AutAlive}[s,t] \in \{0,1\} \end{align*}

where \(\mathcal{V}_{AutAlive}[s,t] = 1\) indicates that the system is autonomous at stage \(s\) and time \(t\), and \(\mathcal{V}_{AutAlive}[s,t] = 0\) indicates that autonomy has been lost.

To ensure that autonomy represents a continuous operation until the first failure, a prefix structure is enforced:

\begin{align*} \mathcal{V}_{AutAlive}[s,t] \ge \mathcal{V}_{AutAlive}[s,t+1] \quad \forall s \in \mathcal{S}_{Stage},~ t \in \mathcal{S}_{TimeHorizon} \end{align*}

The total autonomy duration per stage is defined as

\begin{align*} \mathcal{V}_{Autonomy}[s] \in \mathbb{R}_0^+, \quad \mathcal{V}_{Autonomy}[s] = \sum\limits_{t \in \mathcal{S}_{TimeHorizon}} \mathcal{V}_{AutAlive}[s,t] \end{align*}

This variable measures the number of consecutive time steps during which the system remains autonomous.

While the system is autonomous, cross-border imports are prohibited. Let \(imp\_exp\_type[e] = cross\) denote energy carriers that represent external imports. Then:

\begin{align*} \mathcal{V}_{Imp}[s,h,e,t] \le M[s,h,e,t] \cdot (1 - \mathcal{V}_{AutAlive}[s,t]) \end{align*}

for all \((s,h,e) \in \mathcal{S}_{ImpTuple}\) with \(imp\_exp\_type[e] = cross\) and \(t \in \mathcal{S}_{TimeHorizon}\).

Cross-border imports are therefore permitted only once autonomy is no longer maintained.

The autonomy formulation interacts with the unmet demand variable

\begin{align*} \mathcal{V}_{DemandUnmet}[s,h,e,t] \in \mathbb{R}_0^+ \end{align*}

which is defined in the demand model.

While the system is autonomous, unmet demand is forbidden. This requirement is enforced by the following constraint:

\begin{align*} \mathcal{V}_{DemandUnmet}[s,h,e,t] \le M[s,h,e,t] \cdot (1 - \mathcal{V}_{AutAlive}[s,t]) \end{align*}

As a result, all demand must be satisfied throughout autonomous operation.

Two different configurations are supported for the autonomy calculation, depending on whether unmet demand is allowed.

1. Strict autonomy (unmet demand disabled)

In this case, unmet demand is globally forbidden. The system must fully satisfy demand at all times using only available ressources and technologies. Autonomy therefore measures the longest period during which demand can be met without cross-border imports and without any slack.

This formulation is more restrictive and may lead to infeasibility if the system cannot supply all demand

2. Relaxed autonomy (unmet demand enabled)

In this configuration, unmet demand is allowed after autonomy failure. During autonomous operation (\(\mathcal{V}_{AutAlive}[s,t] = 1\)), unmet demand remains forbidden, but once autonomy is lost (\(\mathcal{V}_{AutAlive}[s,t] = 0\)), unmet demand may be used to maintain feasibility.

Consequently, the model can explore trade-offs between system cost and autonomy while avoiding infeasibility through controlled demand shortfalls after autonomy ends.

Energy system model

The energy system model connects the variables from the submodels and defines the objective function for the optimization model. The most important task of the energy system model is to define energy balance equations which ensure that no energy can be created or destroyed within the model. This balance equation is given for all stages \(s \in \mathcal{S}_{Stage}\), hubs \(h \in \mathcal{S}_{Hub}\), ecs \(e \in \mathcal{S}_{Ec}\), and time steps \(t \in \mathcal{S}_{Time}\) by

\begin{align*} 0 =~ & \mathcal{V}_{Imp}^*[s, h, e, t] - \mathcal{V}_{Exp}^*[s, h, e, t] + \mathcal{V}_{DemandSupply}^*[s, h, e, t] \\ +&\sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x) \in \mathcal{S}_{StorTechTuple}}} \big( \mathcal{V}_{StorTechOutflow}[s, h, x, t] - \mathcal{V}_{StorTechInflow}[s, h, x, t] \big) \\ +& \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x) \in \mathcal{S}_{EbmTechTuple}}} \big( \mathcal{V}_{EbmTechOutflow}[s, h, x, t] - \mathcal{V}_{EbmTechInflow}[s, h, x, t] \big) \\ +& \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x, e) \in \mathcal{S}_{ConvTechOut}}} \mathcal{V}_{ConvTechOut}[s, h, x, e, t] - \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x, e) \in \mathcal{S}_{ConvTechIn}}} \mathcal{V}_{ConvTechIn}[s, h, x, e, t] \\ +& \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x, e) \in \mathcal{S}_{HpTechOut}}} \mathcal{V}_{HpTechOut}[s, h, x, e, t] - \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x, e) \in \mathcal{S}_{HpTechIn}}} \mathcal{V}_{HpTechIn}[s, h, x, e, t] \\ +& \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x, e) \in \mathcal{S}_{AtesTechOut}}} \mathcal{V}_{AtesTechOut}[s, h, x, e, t] - \sum\limits_{\substack{x \in \mathcal{S}_{Tech} \\ (s, h, x, e) \in \mathcal{S}_{AtesTechIn}}} \mathcal{V}_{AtesTechIn}[s, h, x, e, t] \\ +& \mathcal{V}_{NetHubOut}^*[s, h, e, t] - \mathcal{V}_{NetHubIn}^*[s, h, e, t] \end{align*}

In this constraint, a star symbol has been added to simplify the notation. Take for example a tuple \((s, h, e) \in \mathcal{S}_{Stage} \times \mathcal{S}_{Hub} \times \mathcal{S}_{Ec}\). Then we define

\begin{align*} \mathcal{V}_{Imp}^*[s, h, e, t] = \left \{ \begin{array}{rl} \mathcal{V}_{Imp}[s, h, e, t], &\text{if } (s, h, e) \in \mathcal{S}_{ImpTuple} \\ 0, &\text{else} \end{array} \right . \end{align*}

and similarily for other star-notated variables where the tuple \((s, h, e)\) might not belong to the respective definition set.

We also have to keep in mind that the supply variable \(\mathcal{V}_{DemandSupply}\) has already been constrained by the demand model in case of a sum-type demand. For profile demands, the constraint needs to be set here in order to include load shedding and load shifting contributions:

\begin{align*} \mathcal{V}_{DemandSupply}&[s, h, e, t] = demand\_profile[s, h, e, t] - \mathcal{V}_{LoadShedding}^*[s, h, e, t] + \mathcal{V}_{LoadShifting}^*[s, h, e, t] \qquad \text{if } (s, h, e) \in \mathcal{S}_{DemandProfileTuple} \end{align*}

Here, \(demand\_profile\) is the demand-profile from demands.yaml and \(\mathcal{V}_{LoadShedding}\), \(\mathcal{V}_{LoadShifting}\) come from the load shedding model and load shifting model respectively.

In addition to energy balancing, the energy system model has the job of defining the overall CO2 emissions which are gathered in a global fixed variable

\begin{align*} \mathcal{V}_{SystemCo2}: \mathcal{S}_{Stage} &\to \mathbb{R}, \\ \mathcal{V}_{SystemCo2}[s] &= \mathcal{V}_{TechCo2Total}[s] + \mathcal{V}_{ImpCo2Total}[s] - \mathcal{V}_{ExpCo2Total}[s] + \mathcal{V}_{NetTechCo2Total}[s] \end{align*}

Every stage has the potential to include CO2 penalization costs via the optional parameter \(co2\_price\) from stages.yaml which are defined as:

\begin{align*} \mathcal{V}_{SystemCostCo2Penalty} \in \mathbb{R}, \quad \mathcal{V}_{SystemCostCo2Penalty} = \sum\limits_{s \in \mathcal{S}{Stage}} co2\_price[s] \cdot \mathcal{V}_{SystemCo2}[s] \end{align*}

In addition, stages.yaml contains the parameters \(co2\_min\) and \(co2\_max\) which result in the following straightforward constraints for the CO2 emission amounts:

\begin{align*} co2\_min[s] \le \mathcal{V}_{SystemCo2}[s] \le co2\_max[s] \end{align*}

Finally, when CO2 emissions are to be chosen as an objective function, they need to be gathered in a single variable which is given by

\begin{align*} \mathcal{V}_{SystemCo2Total} \in \mathbb{R}, \quad \mathcal{V}_{SystemCo2Total} = \sum\limits_{s \in \mathcal{S}{Stage}} \mathcal{V}_{SystemCo2}[s] \end{align*}

Moving from system CO2 emissions to system costs, these are collected in a similar global variable

\begin{align*} \mathcal{V}_{SystemCost} \in~ &\mathbb{R}, \\ \mathcal{V}_{SystemCost} =~ &\mathcal{V}_{TechCostTotal} + \mathcal{V}_{ConvTechCostTotal} + \mathcal{V}_{ImpCostTotal} - \mathcal{V}_{ExpProfitTotal} \\ &~ + \mathcal{V}_{LoadSheddingCostTotal} + \mathcal{V}_{LoadShiftingCostTotal} + \mathcal{V}_{NetTechCostTotal} + \mathcal{V}_{SystemCostCo2Penalty} \end{align*}

This variable is a possible objective function of the system alongside the CO2 total emissions.