Solving a multi-horizon stochastic facility location problem with capacity expansion

Abstract

Motivated by the real-world problem of locating hydrogen production facilities and assessing different production technologies which differ in the flexibility of operations and costs, we study a multi-stage multi-horizon stochastic facility location problem with capacity expansion. The objective is to minimize the expected sum of investment, production and distribution costs while satisfying customer demand. The multi-horizon formulation allows to capture the effects of both strategic as well as operational uncertainty on the location decisions. Strategic uncertainty is related to uncertain future demand level, while operational uncertainty is related to uncertain future electricity prices resulting in uncertain production costs. We consider multiple production technologies that have different operational characteristics, but can be combined at a location and operated in parallel. To solve the problem, we implement and compare two solution methods: linear relaxation with a restricted MIP approach and Lagrangian relaxation with a two-step restricted MIP approach. We apply the solution approaches to instances based on the real-world problem of locating hydrogen production in Norway. The results show that both approaches can find near-optimal solutions for small instances, but the Lagrangian-based solution method can also find high-quality solutions for larger instances.

1 Introduction

The purpose of facility location problems is to find a trade-off between strategic investment costs in production facilities and operational production and distribution costs (Fernández and Landete 2019). In traditional deterministic approaches, all input parameters are known with certainty and the objective is to minimize cost while satisfying demand. However, when considering the long-term character of investment decisions, several parameters may be uncertain when the investment decisions must be made (Snyder 2006). Uncertainty on a strategic level might be related to the long-term development of future demands, prices, and technology development, among others. There might also be uncertainty on the operational level, such as short-term variations in costs, prices, and demand, that affect the operational decisions, but not necessarily the strategic decisions (Kaut et al. 2014).

A multi-stage multi-horizon stochastic programming model allows adjusting the long-term strategic decisions based on the realization of stochastic parameters, similar to traditional stochastic multi-stage models. But the multi-horizon stochastic scenario tree embeds operational uncertainty in the strategic decision nodes, instead of including it in the multi-stage scenario tree. The multi-horizon approach thus considerably reduces the size of the scenario tree, and of the problem instance (Kaut et al. 2014). The approach therefore allows to embed a more detailed operational problem in each strategic decision node. But this increased level of detail disregards any potential links between the operational decisions embedded in a strategic decision node and the subsequent strategic decision nodes. The links between the strategic decision nodes, however, are preserved.

Hydrogen production is an important step to meet the emission goals set in United Nations (2015). In order to mark hydrogen as a green energy carrier, it must be produced using electrolysis with energy from renewable sources (IRENA 2020a). Currently, there are two commercially available production technologies: Alkaline (AL) electrolyzers and Proton Exchange Membrane (PEM) electrolyzers. As of today, PEM electrolyzers provide higher flexibility in their production range but are more expensive in both investment and operation compared to AL electrolyzers (Andrenacci et al. 2022). However, the costs of producing hydrogen from electrolysis mainly depend on electricity prices (NEL Hydrogen 2018). As electricity prices are volatile throughout the year and uncertain (Nord Pool 2024), a more flexible production technology might be beneficial as production can be reduced during periods with high electricity prices, whereas the installed capacity can be fully utilized during periods with low electricity prices. Choosing a production technology is thus a critical strategic decision, which will affect the cost of hydrogen and the ability of the production system to adapt to short-term uncertainty.

Motivated by the real-world problem of satisfying future hydrogen demand in Norway, we study the problem of locating production facilities for hydrogen in Norway, taking into account different production technologies. Future hydrogen demand is expected to increase over the next years, but annual demand is uncertain (DNV GL 2021). To complicate matters further, electricity prices differ throughout Norway and are highly volatile, adding uncertainty to the operational costs resulting from day-to-day production. Formulating the problem as a multi-stage multi-horizon stochastic facility location problem allows capturing flexibility in operational decisions and thus analyzing the impact of both strategic, long-term uncertainty in demand and operational, short-term uncertainty in production costs on the optimal decisions regarding optimal location, capacity and technology of production facilities.

Our contributions are threefold: First, we formulate a multi-stage multi-horizon stochastic facility location problem with capacity expansion considering uncertainty in demand and electricity prices, representing strategic and operational uncertainty, respectively. Second, we propose two heuristic solution methods and compare their performance: using linear relaxation with a restricted MIP (R-MIP) approach and a solution method based on Lagrangian relaxation together with a two-step R-MIP heuristic to find feasible solutions of high quality for large-scale instances with millions of binary variables. Third, we solve instances based on the real-world problem of locating hydrogen production in Norway and provide managerial insight into how different electricity prices affect the optimal production technology.

The remainder of this paper is structured as follows: we provide a literature review on stochastic facility location problems with capacity expansion in Sect. 2. The mathematical formulation of the problem is then presented in Sect. 3 and the solution method is discussed in Sect. 4. The case study description is provided in Sect. 5. Computational results and conclusions are discussed in Sects. 6 and 7, respectively.

2 Literature review

We structure the literature review in two main parts in order to distinguish between the two modelling choices: one considering the traditional stochastic scenario tree and the second considering the stochastic multi-horizon scenario tree. In Sect. 2.1, we review the literature on traditional stochastic facility location problems. Papers dealing with stochastic multi-stage, multi-horizon problems are discussed in Sect. 2.2.

2.1 Stochastic facility location and capacity expansion problems

Stochastic formulations considering some or all of the input parameters as uncertain are a natural extension of deterministic models. For an overview of deterministic multi-period facility location and capacity expansion problems, see the reviews by Melo et al. (2009), Arabani and Farahani (2012), and Nickel and Saldanha-da Gama (2019). Traditionally, two-stage stochastic facility location problems with recourse have been formulated as single-period problems, see e.g. the reviews by Owen and Daskin (1998), Snyder (2006), Govindan et al. (2017), Correia and Saldanha-da Gama (2019).

Multi-period facility location problems with capacity expansion under uncertainty have only recently started to receive attention in the research literature. In a two-stage formulation, the first-stage decision is usually to decide when, where and with which capacity to open a facility. In the second stage, demand is revealed and capacity expansion decisions as well as production and/or distribution decisions to satisfy customer demand are made. Some of the first examples of two-stage stochastic programming formulations are Correia and Melo (2021), Štádlerová et al. (2022), and Štádlerová et al. (2023).

Multi-stage stochastic problems are required when uncertainty is progressively revealed during the planning horizon and hence, they are inherently multi-period problems (Correia and Saldanha-da Gama 2019). In a multi-stage formulation, the opening and capacity expansion decisions can be made in each stage. However, these problems tend to become extremely large, as the size of the scenario tree usually grows exponentially when increasing the number of stages (Birge and Louveaux 2011). There are only a few studies discussing capacity expansion under uncertainty using a multi-stage stochastic formulation. One of the first studies on multi-stage stochastic capacity expansion in supply chain design is provided by Ahmed et al. (2003). Due to the complexity of the problem, only small instances can be solved. A multi-stage stochastic facility location problem with modular capacity adjustments is studied in Štádlerová et al. (2025). Using Lagrangian relaxation, the authors show that the problem can be solved for relatively large instances with up to 4.6 million binary variables and 30.8 million continuous variables.

2.2 Multi-stage multi-horizon stochastic problems

In multi-stage stochastic problems, the operational activities are usually aggregated which enables solving larger instances. However, highly aggregated operational decisions can have a negative impact on strategic decisions, especially, when the operational conditions are subject to large short-term variations see e.g.,, (Schütz et al., 2009). As capturing both long-term and short-term uncertainty in a single scenario tree will often lead to a large multi-stage stochastic scenario tree, decoupling strategic and operational uncertainty in multi-horizon models has therefore been suggested to reduce the problem size. Multi-horizon models consider operational decisions and subsequent strategic decisions to be independent of each other and thus allow for modelling the operational part of the problem in more detail. This multi-horizon approach was first introduced in Hellemo et al. (2012) and Kaut et al. (2014). The approach has been successfully applied to energy system design problems (ESD) where the objective is to minimize costs while achieving emission reduction targets considering short-term uncertainty in the availability of renewable energy. Su et al. (2015) for example, present a multi-stage multi-horizon stochastic equilibrium model considering strategic uncertainty in the natural gas costs and operational uncertainty in the availability of renewables.

Examples of other applications include Escudero and Monge (2018), who present a framework for representing uncertainty in non-symmetric scenario trees considering both strategic and operational uncertainty. The framework is applied in capacity expansion planning, minimizing the expected net present value. The authors study a general integer capacity expansion problem (CEP) in supply chains or networks with focus on risk averse measures, modelling investment decisions as integer variables while capacities are step continuous variables. Alonso-Ayuso et al. (2020) present a multi-stage multi-horizon approach in the context of forestry supply chain planning (FSCP) with strategic uncertainty in timber production and operational uncertainty in the timber price and demand.

Multi-stage multi-horizon problems are difficult to solve. Maggioni et al. (2020) discuss the bounds from traditional stochastic programs in multi-horizon stochastic problems considering an energy-related case study. Castro et al. (2023) introduce a specialized interior-point method for large supply network design (SND) problems with more than 800 million linear variables. Zhang et al. (2024) compares the performance of Benders and Lagrangian decomposition for a energy-related multi-horizon stochastic problem with continuous variables and long-term uncertainty in CO₂ tax and budget as well as in long-term power demand and operational uncertainty in the availability of renewable energy sources. The authors show that a parallel implementation of Lagrangian decomposition method is efficient for large instances even if the method is very sensitive to parameter tuning. Strong methodological contributions regarding heuristics for multi-stage multi-horizon problems are presented in Cadarso et al. (2018) for the rapid transit network design problem (RTND) introducing a decomposition matheuristic based on scenario clustering. Escudero and Monge (2021, 2023) solve different integer capacity expansion problems, where all capacity modules have the same modular capacity and are technologically identical. The investment and expansion decision consists of how many modules to install. Escudero and Monge (2021) present a matheuristic for solving the capacitated multiple allocation hub network problem (HLP), while Escudero and Monge (2023) present a matheuristic for solving a multi-stage multi-scale facility location with capacity expansion problem (FLCE), accounting also for multiple products. The strategic uncertainty is considered in facility building costs and residual values, whereas the operational uncertainty is related to demand, raw material costs, and facility capacity disruptions. The proposed algorithm is tested on instances with up to 930 binary and 310 general integer variables providing good feasible solutions within 7 hours.

2.3 Comparison and contribution

Only a few papers study multi-stage multi-horizon stochastic problems considering facility location with capacity expansion (FLCE) and modular capacities. Table 1 summarizes the multi-stage multi-horizon problems discussed above.

In this paper, we present a multi-stage multi-horizon facility location problem with integer capacity expansion. Capacity expansion decisions consist of opening a new facility and choosing its technology and capacity from a set of available modular capacities. Specific investment and production costs for each modular capacity allows for modelling of economies of scale. To solve the problem, we present two solution methods: linear relaxation with a restricted MIP approach and Lagrangian relaxation with a two-step restricted MIP approach. We further use strengthening constraints to strengthen the LP relaxation bound.

Table 1 Selected works on multi-stage multi-horizon problems

3 Problem description and mathematical formulation

We address the problem of locating production facilities during the planning horizon, considering multiple available technologies and capacity expansion. We take into account uncertainty both related to long-term future demand and in short-term production costs. First, we provide a detailed problem description in Sect. 3.1, followed by an introduction to the modelling approach in Sect. 3.2. Finally, the mathematical model is presented in Sect. 3.3.

3.1 Problem description

We consider a multi-stage multi-horizon stochastic facility location problem with two levels of uncertainty. On the strategic level, we deal with uncertain demand, while on the operational level, production costs are subject to uncertainty. The uncertainty related to demand is revealed in strategic nodes at the beginning of each stage where each node has an embedded operational production model subject to a set of operational scenarios. The detailed structure of the scenario tree is discussed in Sect. 3.2. The objective is to minimize the total expected costs of satisfying the customers’ demand. The total costs consist of investment (both opening and capacity expansion), production, and distribution costs, together with penalty costs for demand shortfall and overproduction.

On the strategic level, the decisions are related to where and when to open new facilities, as well as which capacity and technology to install. Capacity can be expanded by opening new facilities with no restrictions on whether a facility has been opened at that location before. Combining facilities with different technology at one location is also possible. However, at each stage at most one new facility can be opened at a given location. Once a facility has been opened, it must remain open in all subsequent nodes. Consequently, different technologies can be installed at one location and operated simultaneously. The investment and operational costs of a facility depend on its capacity, technology and the time of the investment, i.e. the stage when the facility has been opened.

In our problem, we consider two types of production technologies, representing PEM and AL electrolyzers. The production technologies are distinguished based on their operational characteristics: first, we consider continuous production allowing for arbitrary production quantities between zero (included) and maximum capacity. Second, we consider a semi-continuous production technology. The semi-continuous technology enables continuous production between a strictly positive minimum production requirement and maximum capacity. Still, this technology can be temporarily switched off, but has then to stay offline for a given number of operational time periods before it can be switched on again. Thus, the choice of technology affects to a high degree the possible operational decisions.

On the operational level, we consider uncertainty in production costs. Specifically in our case, uncertain production costs are caused by uncertain electricity prices.

The hydrogen is distributed to customer locations, but not all facility locations can serve all customer locations. Distribution costs are linear in transported volume, but the unit distribution cost depends on the distance between the production facility and the location of the customer. Unsatisfied demand as well as overproduction are penalized.

3.2 Modelling approach

We consider modular capacities where investment costs as well as production costs, depend on both the chosen capacity level and technology (see, e.g., Correia and Captivo 2003). Production costs are modelled as a linear function and depend on technology, capacity and the time of opening the facilities since newer facilities produce at lower cost due to technological development. A common approach in the literature is to model changes in capacity level as a jump between available discrete capacities and modify the existing facility (Jena et al. 2015, 2016, 2017; Sauvey et al. 2020; Štádlerová and Schütz 2021; Štádlerová et al. 2024). However, expanding capacity can also be modelled as opening additional facilities (see, e.g., Shulman (1991); Dias et al. (2007); Ghiani et al. (2002)). In this paper, we model capacity expansion similar to Shulman (1991) by opening a new facility at a given location. Facilities of different size and technology can be combined at any location. This modelling approach accounts for economies of scale in investment and operation as well as technology development over time where facilities installed in later stages can be more efficient than facilities opened in earlier stages.

We model strategic and operational uncertainty by means of a multi-stage multi-horizon stochastic scenario tree. It is composed of stages where strategic investment decisions (both opening and capacity expansion) are made in strategic decision nodes. Strategic decisions are made at the beginning of each stage when strategic uncertainty is revealed. Decisions related to production and distribution are made in the operational periods embedded in the strategic decision nodes as operational nodes. For strategic decisions we use periods with a duration of 5 years in every stage, whereas operational decisions reflect daily production planning on an hourly basis. Figure 1 illustrates such a multi-horizon tree where strategic decision nodes, $n=1\dots 7$, are depicted by squares, while operational nodes are marked with circles. The tree has three stages, and operational subtrees with four periods attached to each strategic node. The operational subtrees consist of operational scenarios, i.e. realizations for the short-term uncertainty in the operational nodes.

Fig. 1

Multi-stage multi-horizon stochastic scenario tree

A strategic decision node can aggregate multiple strategic time periods (i.e. years) characterized by equal costs and demand. The number of aggregated time periods must be equal for all nodes belonging to the same stage, but does not have to be the same across stages. In operational nodes, daily hydrogen production is modelled with hourly production decisions. The operational scenarios for uncertain electricity prices are embedded in the strategic decision nodes, such that each strategic decision node can be interpreted as a two-stage stochastic programming problem, where the first-stage decision is to determine the capacity of the facility and the second-stage decisions cover production and distribution. Kaut et al. (2014) emphasize that this modelling approach is suitable if (1) there is a large difference in time scales for strategic and operational decisions and (2) the operational decision in subsequent strategic nodes do not depend on the operational decisions of the current stage.

As we cannot solve instances with a full year of operational decisions, we use representative periods to estimate the operational production and distribution costs for the strategic investment and expansion decision. A representative period is a realization of hourly electricity prices for one day. An operational scenario can consist of one or more representative periods. For example, we can use four representative periods, where each period represents a different season.

3.3 Mathematical formulation

Let us first introduce the following notation:

3.3.1 Sets

$\mathcal {I}$:

Set of candidate facility locations;

$\mathcal {J}$:

Set of customer locations;

$\mathcal {K}$:

Sorted set of available discrete capacity levels, $\mathcal {K} = \{1,2,...,\overline{K}\}$;

$\mathcal {N}$:

Set of ordered nodes of the strategic tree;

$\mathcal {A}_n$:

Set of ancestor nodes to node n, including node n;

$\mathcal {P}$:

Set of possible production technologies;

$\mathcal {P}^C$:

Set of continuous production technologies, $\mathcal {P}^C \subseteq \mathcal {P}$;

$\mathcal {P}^S$:

Set of semi-continuous production technologies, $\mathcal {P}^S \subseteq \mathcal {P}$;

$\mathcal {R}$:

Set of representative periods;

$\mathcal {S}_{n}$:

Set of operational scenarios at strategic node $n \in \mathcal {N}$;

$\mathcal {T}_r$:

Set of time periods in each representative period $r \in \mathcal {R}$, $\mathcal {T}_r = \{1,2,...,\overline{T}_r\}$.

3.3.2 Parameters and coefficients

$B_{p}$:

Minimum production utilization of technology $p \in \mathcal {P}^S$;

$C_{pk}^{n}$:

Annualized investment costs for a facility opened in node $n \in \mathcal {N}$ with technology $p\in \mathcal {P}$ at capacity level $k\in \mathcal {K}$ acquired in stage $h \in \mathcal {H}$;

$D_{jr}^{n}$:

Demand at customer location $j\in \mathcal {J}$ in representative period $r\in \mathcal {R}$ at strategic node $n \in \mathcal {N}$;

$E_{isrt}^{n}$:

Cost of electricity at location $i \in \mathcal {I}$ in representative period $r \in \mathcal {R}$ and time period $t \in \mathcal {T}$, at strategic node $n \in \mathcal {N}$ in scenario $s\in \mathcal {S}_{n}$;

$F_{p}^{na}$:

Power consumption in node $n \in \mathcal {N}$ for production technology $p \in \mathcal {P}$ acquired in stage $a \in \mathcal {A}_{n}$;

$K_{k}$:

Maximum production quantity at capacity level $k \in \mathcal {K}$;

$L_{ij}$:

1 if demand at location $j\in \mathcal {J}$ can be served from facility $i\in \mathcal {I}$;

$M^u$:

Penalty cost per unit of demand shortfall;

$M^o$:

Penalty cost per unit of overproduction;

$T_{ij}$:

Unit distribution costs from facility location $i\in \mathcal {I}$ to customer $j\in \mathcal {J}$;

$W_{p}$:

Number of operational periods a facility has to remain offline after turning off production for technology $p\in \mathcal {P}$;

$W_{p}^{\prime }$:

Number of operational periods a facility has to remain operational after turning on production for technology $p\in \mathcal {P}$;

$\pi _{n}$:

Probability of strategic node $n \in \mathcal {N}$;

$\tau _{n}$:

Number of strategic periods aggregated in node n;

$\omega _{s}^{n}$:

Probability of an operational scenario $s \in \mathcal {S}_{n}$ related to strategic node $n \in \mathcal {N}$.

3.3.3 Decision variables

The mathematical model uses the following decision variables:

$y_{ipk}^{n}$:

1 if a facility is opened in node $n \in \mathcal {N}$ at location $i \in \mathcal {I}$ with technology $p \in \mathcal {P}$ and capacity level $k \in \mathcal {K}$, 0 otherwise;

$d_{ipksrt}^{na}$:

1 if there is production at location $i\in \mathcal {I}$, using technology $p\in \mathcal {P}^S$ at the capacity level $k\in \mathcal {K}$ in node $n \in \mathcal {N}$, in season $r\in \mathcal {R}$ and time period $t\in \mathcal {T}$ of scenario $s\in \mathcal {S}_{n}$, for facility opened in node $a \in \mathcal {A}_{n}$, 0 otherwise;

$o_{ipksr}^{na}$:

Amount of overproduction at location $i \in \mathcal {I}$ using technology $p \in \mathcal {P}$ and capacity $k \in \mathcal {K}$ in node $n \in \mathcal {N}$ in season $r\in \mathcal {R}$ in scenario $s \in \mathcal {S}_{n}$ for facility opened in node $\mathcal {A}_n$;

$q_{ipksrt}^{na}$:

Production quantity at facility $i\in \mathcal {I}$ in node $n \in \mathcal {N}$ opened in node $a \in \mathcal {A}_{n}$, using technology $p \in \mathcal {P}$ and capacity $k \in \mathcal {K}$ in season $r\in \mathcal {R}$ at time period $t \in \mathcal {T}$ of scenario $s \in \mathcal {S}_{n}$;

$u_{jsr}^{n}$:

Amount of demand shortfall at customer location $j\in \mathcal {J}$ in node $n \in \mathcal {N}$ in season $r\in \mathcal {R}$ in scenario $s\in \mathcal {S}_{n}$;

$x_{ijpksr}^{na}$:

Amount of customer demand at location $j\in \mathcal {J}$ in node $n \in \mathcal {N}$ satisfied from facility $i\in \mathcal {I}$ opened in node $a \in \mathcal {A}_{n}$ using technology $p \in \mathcal {P}$ with installed capacity $k \in \mathcal {K}$, in season $r\in \mathcal {R}$ and scenario $s\in \mathcal {S}_{n}$.

Note that we choose to model investment decisions using binary variables $y_{ipk}^{n}$ which value is one if facility i is opened in node n. For operational variables $\varvec{d,o,q,x}$, we use a pair of indexes na, where $n \in \mathcal {N}$ indicates the current node, while $a \in \mathcal {A}_n$ indicates the node when the facility was first opened.

This structure allows us to capture differences in facility properties and cost functions depending on their opening period. For instance, facilities established in later stages may operate with higher efficiency, while those opened earlier cannot take advantage of technologies that become available only in the future. Since capacity expansion is modeled as the opening of a new facility at the same location, it is possible to have multiple facilities at one location operating simultaneously, distinguished only by their opening period and corresponding efficiency.

The multi-stage multi-horizon stochastic programming model for our problem can then be formulated as:

$$\begin{aligned} & \min \sum _{n \in \mathcal {N}} \pi _{n}\tau _{n} \left( \sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} C_{pk}^{a} y_{ipk}^{a} \right. \nonumber \\ & \left. +\sum _{s \in \mathcal {S}_{n}} \omega ^{n}_{s} \left[ \sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}}\sum _{r\in \mathcal {R}} \sum _{t \in \mathcal {T}_r} E_{isrt}^{n}F_{p}^{na} q_{ipksrt}^{na} \right. \right. \nonumber \\ & \left. \left. +\sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{j \in \mathcal {J}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} \sum _{r\in \mathcal {R}} T_{ij}x_{ijpksr}^{na} \right. \right. \nonumber \\ & \left. \left. +\sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}}\sum _{r\in \mathcal {R}} M^o o_{ipksr}^{na}+\sum _{j \in \mathcal {J}}\sum _{r\in \mathcal {R}} M^u u_{jsr}^{n}\right] \right) , \end{aligned}$$

(1)

subject to:

$$\begin{aligned}&\sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} y_{ipk}^{n} \le 1, \,\,\,\,\,n \in \mathcal {N}, i \in \mathcal {I}, \end{aligned}$$

(2)

$$\begin{aligned}&\sum _{i \in \mathcal {I}} \sum _{a \in \mathcal {A}_{n}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} x_{ijpksr}^{na}+u_{jsr}^{n} = D_{jr}^{n}, \,\,\,\,\,n \in \mathcal {N}, j \in \mathcal {J}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(3)

$$\begin{aligned} & x_{ijpksr}^{na} \le L_{ij} D_{jr}^{n}y_{ipk}^{a}, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, j \in \mathcal {J}, p \in \mathcal {P}, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(4)

$$\begin{aligned} & \sum _{j \in \mathcal {J}} x_{ijpksr}^{na} + o_{ipksr}^{na} = \sum _{t \in \mathcal {T}} q_{ipksrt}^{na}, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(5)

$$\begin{aligned} & q_{ipksrt}^{na} \le K_{k}y_{ipk}^{a}, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^C,k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r, \end{aligned}$$

(6)

$$\begin{aligned} & d_{ipksrt}^{na} \le y_{ipk}^{a}, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r, \end{aligned}$$

(7)

$$\begin{aligned} & q_{ipksrt}^{na} \le K_{k}d_{ipksrt}^{na}, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S,k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r, \end{aligned}$$

(8)

$$\begin{aligned} & q_{ipksrt}^{na} \ge B_{p} K_{k} d_{ipksrt}^{na}, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r, \end{aligned}$$

(9)

$$\begin{aligned} & d_{ipksr(t-1)}^{na} - d_{ipksrt}^{na} \le 1 - \frac{1}{W_{p}} \sum _{t^{\prime } = t+1}^{t+W_p} d_{ipksrt^{\prime }}^{na}, \nonumber \\ &\quad\quad n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r: t \le \overline{T}_r - W_p, \end{aligned}$$

(10)

$$\begin{aligned} & d_{ipksr(t-1)}^{na} - d_{ipksrt}^{na} \le \frac{1}{W_{p}^{\prime }} \sum _{t^{\prime } = t-W_p^{\prime }}^{t-1} d_{ipksrt^{\prime }}^{na}, \nonumber \\ &\\ &\quad\quad n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r: t \ge W_p^{\prime }, \end{aligned}$$

(11)

$$\begin{aligned} & d_{ipksrt}^{na} \in \{0,1\}, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I},p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in {T}_e, \end{aligned}$$

(12)

$$\begin{aligned} & o_{ipksr}^{na} \ge 0, \,\,\,\,\, n \in \mathcal {N}, i \in \mathcal {I}, p \in \mathcal {P}, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(13)

$$\begin{aligned} & q_{ipksrt}^{na} \ge 0, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r, \end{aligned}$$

(14)

$$\begin{aligned}&0\le u_{jsr}^{n} \le D_{jr}^n, \,\,\,\,\, n \in \mathcal {N}, j \in \mathcal {J}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(15)

$$\begin{aligned} & x_{ijpksr}^{na} \ge 0, \,\,\,\,\, \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, j \in \mathcal {J}, p \in \mathcal {P}, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(16)

$$\begin{aligned} & y_{ipk}^{n} \in \{0,1\}, \,\,\,\,\, n \in \mathcal {N}, i \in \mathcal {I}, p \in \mathcal {P}, k \in \mathcal {K}, \end{aligned}$$

(17)

The objective function (1) is the expected sum of annualized investment costs at the strategic level plus the expected sum of operational costs. The operational costs depend on operational scenario $s \in \mathcal {S}$ and consist of production costs, distribution costs and penalty costs for overproduction and demand shortfall, respectively.

Constraints (2) ensure that at most one facility can be opened in each node in each stage. Constraints (3) ensure that demand is either satisfied or penalized as demand shortfall. Constraints (4) specify which facilities can satisfy which customers. Further, they are formulated in the form of strong inequalities to improve the LP bound (see also Jena et al. 2016). Constraints (5) guarantee that the production quantity is either delivered to customers or penalized as overproduction. Note that defining the variable $q_{ipksrt}^{na}$ with a capacity index k further strengthens the formulation. Constraints (6) define the maximum production capacities for technologies with continuous production $\mathcal {P}^C$. Note that the minimum production quantity for continuous production is equal to zero and since the production variable $\varvec{q}$ is defined as non-negative, the minimum production quantity does not need to be defined explicitly. Constraints (7) link the production of facilities with technology $\mathcal {P}^S$ only to existing facilities. Further, constraints (8) and (9) specify the maximum and minimum production quantities for the semi-continuous production technology $\mathcal {P}^S$ with installed capacity k, respectively. Constraints (10) ensure that after turning off semi-continuous production, the facility will be offline for at least $W_p$ consecutive time periods, while Constraints (11) allow the semi-continuous production to be turned off after being on for at least $W_p^{\prime }$ consecutive time periods. Note that production for the continuous technology is defined between zero and the maximum production quantity for given capacity. Therefore, there are no requirement on minimum on or off time. Constraints (12)–(17) are the binary and non-negativity constraints.

We provide an alternative model formulation without strengthening constraints (4) and with variables $\varvec{q}$ without capacity index k in Appendix A. We further provide computational results comparing the LP bound and run times of different model formulation in Table 8. The results show that our model formulation (1)–(17) provides a stronger LP bound and shorter run times than the model without strengthening constraints.

4 Solution approach

The presented problem is a large mixed integer problem with many binary variables. Due to strengthening constraints (4), the LP relaxation provides strong bounds. However, for our large instances, even the LP relaxation cannot be solved to optimality. Therefore, we present and compare two heuristic solution methods to solve the problem. In Sect. 4.1, we present an R-MIP approach based on the optimal solution to the LP relaxation of the original problem. Section 4.2 presents a solution method based on Lagrangian relaxation together with a two-step R-MIP approach to construct an upper bound.

4.1 Linear relaxation and R-MIP

We solve the LP relaxation of the original problem (1)–(17) to obtain a lower bound and a starting point for our heuristic. Our heuristic is based on an R-MIP approach: We define $\Phi$ as a threshold value for binary variables and set variables $y\ge \Phi$ in the LP relaxation equal to one. We then solve the remaining R-MIP to find a feasible solution and upper bound for the original problem.

Choosing a value of $\Phi$ is a critical part of this approach since it can be difficult to find the balance between fixing too many variables and, as a result, finding a solution of poor quality, and fixing too few variables and solving a large problem which remains intractable. Initial tests have shown that $\Phi = 0.75$ is a good compromise in our case (see Appendix B for more details).

4.2 Lagrangian relaxation

Using Lagrangian relaxation, and relaxing the demand constraints in particular, has been a popular choice in facility location problems, (see, e.g., Shulman 1991; Jena et al. 2016, 2017; Štádlerová et al. 2024, 2023). We formulate the relaxed problem in Sect. 4.2.1 before discussing the Lagrangian multipliers in section 4.2.2. Finally, we present our upper bound heuristic in Sect. 4.2.3. The different steps of the Lagrangian-based solution approach are illustrated in Fig. 2.

Fig. 2

Lagrangian relaxation with two-step R-MIP approach

4.2.1 Solving the relaxed problem

We relax the demand constraints (3) and define $\lambda _{jsr}^n$ as the matrix of Lagrangian multipliers. As the demand constraints are the only constraints linking all facility locations together, the resulting relaxed problem becomes separable in facility locations. The Lagrangian relaxation $LR(\varvec{\lambda })$ is given as follows:

$$\begin{aligned} & \min \sum _{n \in \mathcal {N}} \pi _{n}\tau _{n} \left( \sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} C_{pk}^{a} y_{ipk}^{a} \right. \nonumber \\ & \left. \qquad +\sum _{s \in \mathcal {S}_{n}} \omega ^{n}_{s} \left[ \sum _{i \in \mathcal {I}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}}\sum _{r\in \mathcal {R}} \sum _{t \in \mathcal {T}_r} E_{isrt}^{n}F_{p}^{na} q_{ipksrt}^{na} \right. \right. \nonumber \\ & \left. \left. \qquad +\sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{j \in \mathcal {J}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} \sum _{r\in \mathcal {R}} (T_{ij}-\lambda _{jsr}^n)x_{ijpksr}^{na} \right. \right. \nonumber \\ & \left. \left. \qquad +\sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}}\sum _{r\in \mathcal {R}} M^o o_{ipksr}^{na} \right. \right. \nonumber \\ & \left. \left. \qquad +\sum _{j \in \mathcal {J}}\sum _{r\in \mathcal {R}} (M^u - \lambda _{jsr}^n) u_{jsr}^{n} + \sum _{j \in \mathcal {J}} \sum _{r\in \mathcal {R}} \lambda _{jsr}^n D_{jr}^n \right] \right) , \end{aligned}$$

(18)

subject to Constraints (2) and (4)–(17).

For given $\lambda _{jsr}^n$, the expression $\sum _{j \in \mathcal {J}} \sum _{r\in \mathcal {R}} \lambda _{jsr}^n D_{jr}^n$ is constant. For sufficiently large penalties, i.e. $M^u> \lambda _{jsr}^n$, $u_{jsr}^n$ is equal to 0 in the optimal solution to $LR(\varvec{\lambda })$. For $M^u \le \lambda _{jsr}^n$, $u^n_{jsr}$ is equal to $D_{jr}^n$ and the last line of (18) reduces to the constant:

$$\begin{aligned} \sum _{n \in \mathcal {N}}\pi ^n \tau ^n \sum _{s \in \mathcal {S}^n} \omega _s^n M^u D_{jr}^n. \end{aligned}$$

For given multipliers, the last two terms of the objective function (18) reduce to a constant (Štádlerová et al. 2025). The problem becomes separable in facility locations and we can rewrite objective (18) as:

$$\begin{aligned} \sum _{i \in \mathcal {I}}g_i(\varvec{\lambda }) + \sum _{n \in \mathcal {N}} \sum _{s \in \mathcal {S}}\sum _{j \in \mathcal {J}} \sum _{r\in \mathcal {R}} \pi _n \tau _n \omega _s^n \left[ \lambda _{jsr}^n D_{jr}^n + (M^u - \lambda _{jsr}^n) u_{jsr}^{n} \right] , \end{aligned}$$

(19)

where $g_i(\varvec{\lambda })$ is the optimal objective value of the Lagrangian subproblem for a given location $i \in \mathcal {I}$:

$$\begin{aligned} & g_i(\varvec{\lambda }) = \min \sum _{n \in \mathcal {N}} \pi _{n}\tau _{n} \left( \sum _{a \in \mathcal {A}_{n}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} C_{pk}^{a} y_{ipk}^{a} \right. \nonumber \\ & \left. \qquad + \sum _{s \in \mathcal {S}_{n}} \omega ^{n}_{s} \left[ \sum _{a \in \mathcal {A}_{n}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}}\sum _{r\in \mathcal {R}} \sum _{t \in \mathcal {T}_r} E_{isrt}^{n}F_{p}^{na} q_{ipksrt}^{na} \right. \right. \nonumber \\ & \left. \left. \qquad +\sum _{a \in \mathcal {A}_{n}} \sum _{j \in \mathcal {J}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} \sum _{r\in \mathcal {R}} (T_{ij} - \lambda _{jsr}^{na}) x_{ijpksr}^{na} \right. \right. \nonumber \\ & \left. \left. \qquad +\sum _{a \in \mathcal {A}_{n}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}}\sum _{r\in \mathcal {R}} M^o o_{ipksr}^{na}\right] \right) , \end{aligned}$$

(20)

subject to Constraints (2) and (4)–(17) for facility $i \in \mathcal {I}$.

The optimal solution to the Lagrangian subproblem provides the optimal investment and demand allocation decisions in all scenarios for a given facility. Since combining of facilities with different production technologies at the same location is allowed, solving the subproblem using a dynamic programming approach is not computationally tractable for large instances (Shulman 1991). We therefore only solve the LP relaxation of the Lagrangian subproblem. Due to the strengthening constraints (4), the integrality gap for a single facility problem $g_i(\varvec{\lambda })$ is usually below $< 0.05\%$ for our instances. Hence, only solving the linear relaxation of the Lagrangian subproblem still provides a good lower bound.

4.2.2 Lagrangian multipliers

For any given $\varvec{\lambda }$, we obtain a lower bound on (18). In order to obtain the best lower bound, we need to solve the Lagrangian dual problem $LD = \max _{\varvec{\lambda }} LR(\varvec{\lambda })$. We use a cutting plane method with box constraint to iteratively update the multipliers. The method is similar to the one in Štádlerová et al. (2023) and included here for the sake of completeness. We initialize the Lagrangian multipliers as $\lambda _{jsr}^n = \min _{i \in \mathcal {I}}T_{ij}$. In each iteration m, we compute the coordinate of the subgradient $\nabla _{jsr}^{nm}$, where $\nabla _{jsr}^{nm} = D_{jr}^n - \sum _{a \in \mathcal {A}_n} \sum _{i \in \mathcal {I}}\sum _{p \in \mathcal {P}}\sum _{k \in \mathcal {K}}x_{ijpksr}^{nam} - u_{jsr}^{nm}$ and $x_{ijpksr}^{nam}$ is the solution to the Lagrangian subproblem obtained for variables x in iteration m. We then define $L^m = LR(\varvec{\lambda }^m)-\sum _{j \in \mathcal {J}} \sum _{n \in \mathcal {N}} \sum _{s \in \mathcal {S}_n} \sum _{r\in \mathcal {R}} \pi _n \omega _s^n \lambda _{jsr}^{nm} \nabla _{jsr}^{nm}$ and obtain new Lagrangian multipliers by solving the following linear problem:

$$\begin{aligned} \max \phi \end{aligned}$$

(21)

subject to:

$$\begin{aligned} & \phi \le L^g + \sum _{j \in \mathcal {J}} \sum _{n \in \mathcal {N}} \sum _{s \in \mathcal {S}_n} \sum _{e \in \mathcal {R}} \pi _n \omega _s^n \nabla _{jsr}^{ng} \lambda _{jsr}^{n,m+1}, \hfill g = 1, ..., m, \end{aligned}$$

(22)

$$\begin{aligned}&\lambda _{jsr}^{n,m+1} \le \lambda _{jsr}^{nm} + \Delta _{jsr}^{nm}, \hfill j \in \mathcal {J}, n \in \mathcal {N}, s \in \mathcal {S}_n, r\in \mathcal {R}, \end{aligned}$$

(23)

$$\begin{aligned}&\lambda _{jsr}^{n,m+1} \ge \lambda _{jsr}^{nm} - \Delta _{jsr}^{nm}, \hfill j \in \mathcal {J}, n \in \mathcal {N}, s \in \mathcal {S}_n, r\in \mathcal {R}, \end{aligned}$$

(24)

$$\begin{aligned}&\phi \in \mathbb {R}, \lambda _{jsr}^{n,m+1} \in \mathbb {R}. \end{aligned}$$

(25)

Constraints (23) and (24) are the box constraints that limit how much the Lagrangian multipliers can change in each iteration. When the sign of the subgradient element $\nabla _{jse}^{nm}$ changes compared to the previous iteration, we update the box size: $\Delta _{jt}^m = \alpha \Delta _{jt}^{m-1}$, where $0 < \alpha \le 1$. With this step, we reduce the box size in order to speed up the process of finding the optimal multipliers (Marsten et al. 1975).

4.2.3 Calculating the upper bound

Once the Lagrangian relaxation has been solved, we calculate an upper bound for the original problem. To find a feasible solution, we apply a two-step R-MIP approach: We first use the information from solving the Lagrangian relaxation to fix a subset of the location variables y. Let $\alpha _{g}, g = 1,...,m$ be the dual variable associated to constraint (22) for iteration g. The variable $\alpha _g$ can be interpreted as the probability that the solution to the relaxed problem obtained in iteration g is good, (see, e.g., Jena et al. 2017; Štádlerová et al. 2024). We then calculate:

$$\begin{aligned}&\gamma _{ipk}^{n} = \sum _{g = 1}^{m}\alpha _g y_{ipk}^{ng}, \hfill n \in \mathcal {N}, i \in \mathcal {I}, p \in \mathcal {P}, k \in \mathcal {K}, \end{aligned}$$

(26)

where $y_{ipk}^{ng}$ is the solution of the investment decision variable $y_{ipk}^{n}$ in the relaxed problem in iteration g. Note that $y_{ipk}^{ng} \in [0,1]$. We define $\Psi$ as a threshold value and set variables $y_{ipk}^{n}=1$ if $\gamma _{ipk}^{n} \ge \Psi$. Initial tests have shown that $\Psi = 0.85$ works well for our instances, which also corresponds to Jena et al. (2017); Štádlerová et al. (2024). With these variables fixed, we solve a simplified version of the problem that disregards variables $d_{ipksrt}^{na}$ to obtain a feasible solution for all variables $y_{ipk}^{n}$ that have not been fixed. The simplified problem R-MIP1 is given by (1)–(5), (12)–(17), and (27)–(28). We also add the following capacity constraints.

$$\begin{aligned} & q_{ipksrt}^{na} \le K_{k}y_{ipk}^{n}, \hfill \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P},k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r, \end{aligned}$$

(27)

$$\begin{aligned} & q_{ipksrt}^{na} \ge B_p K_{k}y_{ipk}^{n}, \hfill \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S,k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r. \end{aligned}$$

(28)

The formulation of minimum production requirements (28) is, however, stricter for semi-continuous production technology than the original problem using constraints (9). Here, the production level is forced to stay above the minimum production requirements once the facility has been opened. Still, solving this simplified problem results in decisions $y_{ipk}^{n}$ that are feasible for the original problem.

In the second step, we solve R-MIP2, i.e. the original problem (1)–(17) with the variables $y_{ipk}^{n}$ fixed to the values obtained from R-MIP1. With this step, we obtain a feasible solution for all variables. As all $y_{ipk}^{n}$ are now fixed, we can solve the problem independently for each strategic node $n \in \mathcal {N}$, and each operational scenario $s \in \mathcal {S}_n$ using a commercial software.

5 Case study

In this section, we present the input data for our case study regarding the optimal location of hydrogen production facilities in Norway. We first discuss the planning horizon (see Sect. 5.1), while Sect. 5.2 presents the demand and the generation of strategic scenarios. Section 5.3 provides the investment, production, distribution, and penalty costs. Scenarios for operational uncertainty are discussed together with production costs. Finally, Sect. 5.4 presents the candidate facility locations and customer locations.

5.1 Planning horizon

We consider a planning horizon of 15 years, split into 3 stages. Each stage consists of 5 years. On the operational level, we consider up to 4 representative periods representing four seasons (spring, summer, autumn, winter) with a scenario specifying the 24 hourly electricity prices of a representative day within each season.

5.2 Demand

Future hydrogen demand represents strategic uncertainty in our problem. However, hydrogen demand is expected to increase in the coming years. Figure 3 illustrates the minimum and maximum potential demand level during our 15 years planning horizon. The minimum demand is derived from demand in the maritime transportation sector. This demand is based on the estimated energy needs for providing public transportation services (such as ferries and high speed passenger vessels) using hydrogen as fuel (Ocean Hyway Cluster 2020). The maximum demand level combines demand from the maritime sector with estimates for demand from land-based transportation (DNV GL 2019). In the first stage, demand is deterministic and consists exclusively of maritime demand. In the second and third stage, the share of land-based transportation demand is uncertain.

Fig. 3

Development of hydrogen demand during the planning horizon

Daily demand in our scenarios is always between the minimum demand level, $D_h^m$, and the maximum demand level, $D_h^u$, of each stage h. Note that we assume demand to be constant within a stage. To ensure that demand is increasing from one stage to another, we use a growth factor $g>1$. We first sample demand in stage $h = 1$ from the interval $[D_1^m,D_1^u]$. We then sample growth factors $g_2^n$ from the interval $(1,\frac{D_2^u}{D_1}]$ and calculate demand in each node of the second stage as $D_2^n=D_1\cdot g_2^n$. To generate scenarios for the third stage, we repeat this procedure: for every second stage node m, we sample growth factors $g_3^n$ from the interval $(1,\frac{D_3^u}{D_2^m}]$ and calculate demand in the third stage as $D_3^n=D_2^m\cdot g_3^n$. The daily demand is equal in all representative periods.

For our instances, we draw the same number of samples for the growth factors in the strategic nodes of stages 1 and 2. We generate 3, 4, 5, and 8 scenarios, resulting in up to 64 scenarios in stage 3.

5.3 Costs

We assume that investment costs are independent of the location, but both investment costs and production costs change over the planning horizon, e.g. due to technology development. However, the production costs depend on the location.

5.3.1 Investment and maintenance costs

We consider two production technologies: AL and PEM electrolyzers. AL electrolysis is modelled as semi-continuous production technology with a minimum production requirement, whereas PEM electrolysis is modelled as continuous production technology. The investment costs for both technologies are taken from Andrenacci et al. (2022). The development of these costs over time is adjusted based on IRENA (2020b). We consider 9 discrete capacities and use costs from Refsdal and Sindre (2023). Table 2 shows annualized investment costs at each stage assuming a lifetime of 15 years, and Table 3 shows the yearly maintenance costs for both technologies.

Table 2 Investment costs in [mill. €] (Refsdal and Sindre 2023)

Table 3 Maintenance costs [mill. €/year] (Refsdal and Sindre 2023)

5.3.2 Production costs

Electricity costs are a critical part of the production costs when using electrolysis to produce hydrogen. In general, electricity prices in the northern part of Norway are considerably lower and less volatile than in the southern part. We therefore consider specific electricity prices for each facility location depending on its regional electricity market (bidding zone). Figure 4 illustrates the electricity prices for each of the Norwegian bidding zones (NO1–NO5) during a 5-day period (120 hours) in December 2022. Note that prices for NO1, NO2, and NO5 are equal.

Fig. 4

Hourly electricity prices in Norway starting 10^th December 2022

We collect historical day-ahead electricity prices from Nord Pool (2024) for the years 2020, 2021, and 2022. We further group electricity prices based on season (spring, summer, autumn, and winter). Operational scenarios are generated by sampling electricity prices for each bidding zone from the historical data. We generate operational scenarios with one representative period per scenario and operational scenarios with four representative periods per scenario. For instances with only one representative period, we sample electricity price scenarios for each bidding zone directly from the historical data without considering which season the scenario belongs to. For instances with four representative periods per scenario, we sample one electricity price realization from each season.

We generate instances with up to 150 operational scenarios corresponding to 150 representative periods. Note that instances with one representative period and 100 operational scenarios and instances with four representative periods and 25 operational scenarios are equal in size.

We assume the production costs to be linear for both technologies. AL electrolyzers have a production range between $20 - 100\%$. They can be switched off, but have an average start-up time of 2 hours (Refsdal and Sindre 2023). The operational range of a PEM electrolyzer is $0-100 \%$. Power consumption data for both technologies are taken from Andrenacci et al. (2022). Table 4 shows the power consumption per produced kg of hydrogen at each stage and production range.

Table 4 Power consumption and specifications of PEM and AL (Refsdal and Sindre 2023)

5.3.3 Distribution costs

In this paper, we focus on hydrogen distribution by truck. Distribution costs depend solely on the distance and the amount of hydrogen transported. Further, there is a distribution limit of a maximum of 1000km. The distribution cost per kilogram and kilometer are given in Table 5 for different distance intervals. Note that parameter $T_{ij}$ is the unit distribution cost per kg given the distance between facility i and customer j.

Table 5 Distribution costs per kilometer and kilogram (Refsdal and Sindre 2023)

5.3.4 Penalty costs

In order to focus on independent domestic production, the penalty costs are set to a high value, $10^6$ € per unit, to provide a strong incentive to avoid demand shortfall and overproduction.

5.4 Facility and customer locations

In the base case, we consider 17 candidate locations. Figure 5 shows the candidate facility locations based on Ocean Hyway Cluster (2020) and the bidding zones they belong to. Keep in mind that the bidding zone impacts the production costs of the facilities, see Sect. 5.3.2.

Fig. 5

Candidate facility locations and their electricity bidding zones

We start by using the 71 customer locations from Štádlerová et al. (2024). To reduce the complexity of the problem, we aggregate them into 30 customer locations using K-means clustering based on geographical location (Refsdal and Sindre 2023).

6 Computational study

Our algorithm is implemented in Julia 1.8.2. using OpenMPI 4.1.6 for distributed computations. We also enable parallelization in Julia with up to 32 threads. Gurobi 11.0 is used as the commercial solver. All calculations have been carried out on a cluster with 9 nodes running CentOS 7 (kernel 3.10), where each node has two 3.6 GHz Intel Xeon Gold 6244 CPU (8 cores) processors and 384 GB RAM.

We denote our instances with the number of strategic scenarios (N), operational scenarios (S) and representative periods per operational scenario (R), e.g. instance F17D30N25S3R4 consists of 17 candidate locations, 30 demand locations, 25 strategic scenarios, and 3 operational scenarios in each of the 4 representative periods. To give some insights in the size of our problem instances, we report more detailed information in Appendix C in Table 10.

We discuss the solution quality of the proposed solution methods in Sect. 6.1. In Sect. 6.2, we analyze the solution structure.

6.1 Solution quality

We compare our two solution approaches, the linear relaxation with R-MIP (LP + R-MIP) and the Lagrangian relaxation with the two-step R-MIP (LR + 2R-MIP), to Gurobi (Gur). In particular, we compare the run time, the best feasible solution, and the optimality gap. We set a time limit of 48 hours for each method which are divided equally between calculating the lower bound and the upper bound. We set a time limit of 24 hours for solving each of the following: the linear relaxation of the original problem (LP), the R-MIP based on the LP relaxation (R-MIP), and calculating the upper bound with a two-step R-MIP based on the Lagrangian relaxation (UB). The limit for solving the Lagrangian dual is set to 24 hours or 100 iterations (one iteration consists of solving the relaxed problem and updating the Lagrangian multipliers). In Table 6, we report the optimality gaps of Gurobi after 48 hours, the LP + R-MIP approach (Gap LP), and the LR + 2R-MIP approach (Gap LR). We also report the best feasible solution (Best FS) for all tested solution approaches. Further, we show the run times for solving the linear relaxation and the R-MIP as well as solving the Lagrangian dual and the two-step R-MIP. Due to the strengthening constraints, the LP bound provides a strong lower bound for the original problem. Initial tests have shown that, on average, the integrality gap is $0.5\%$ for the original problem formulation (1)–(17). Note that the LP bound is on average $0.2 \%$ higher than the LR bound. However, the LP bound is not available for larger instances.

Table 6 Computational results general formulation

Gurobi is not able to find the optimal solution for any of our instances. The LP relaxation can be solved for 12 instances. Gurobi further finds feasible solutions for these 12 instances within 48 hours. The best solutions from Gurobi are comparable with solutions from LP + R-MIP, and slightly better than LR + 2R-MIP solutions.

In general, the LP + R-MIP approach provides high-quality solutions for the smaller instances with an average optimality gap of $0.56 \%$. However, the R-MIP is not always solved to optimality. If the optimal solution to the R-MIP is not found within 86400 seconds, we report the best found solution.

The LR + 2R-MIP approach also provides high-quality solutions but with a slightly larger average optimality gap of $1.06 \%$. This is mainly due to the fact that the Lagrangian bound is weaker than the LP bound. Still, the best feasible solutions provided by the LP + R-MIP approach are on average only $0.34 \%$ better than solutions provided by the LR + 2R-MIP approach. The run time of LR + 2R-MIP is considerably shorter than LP + R-MIP for all but one instance with 17 candidate locations. The main advantage of the approach based on Lagrangian relaxation over the linear relaxation approach however, is that it also finds valid lower and upper bounds for the large instances with more than 100 operational scenarios. Note that the largest instances with 16 strategic scenarios and 100 operational scenarios have approximately 25 million binary variables.

6.2 Solution structure

In this section, we analyze the technology choice of our proposed model formulation given by (1)–(17), where different technologies can be combined at one location and operated in parallel without restrictions. Below, we focus on the results for the instances with more than 100 operational scenarios since the strategic technology decisions seem to stabilize for these instances. Please note that we do not test for in-sample or out-of-sample stability here, as available RAM limits the size (i.e. the number of scenarios) of instances.

We discuss the best known solution for instance F17D30N9S100R1 as this is the instance for which we obtain the lowest optimality gap among the large problem instances. The technology choices for this instance are identical to the other instances with more than 100 operational scenarios, but the capacity decisions may differ. In Fig. 6, we plot the location and technology decisions as these are identical in all scenarios, but omit the installed capacity. Each of the black underlines represents one stage, and the square indicates the technology installed in that stage at the location. Blue and white colour indicates AL technology, while orange colour is used for PEM.

As shown in Fig. 6, the best-known solution chooses to open only AL facilities in the first stage as they are considerably cheaper in investment and production. In the second stage, AL technology is still cheaper in investment and operation, but we see that PEM technology is being installed in southern Norway. Electricity prices in this area can be much higher than in northern Norway, such that the additional flexibility of PEM allows it to outperform production using AL technology. Hellesylt is an exception as the installed capacity is so small that it operates at capacity limit in all scenarios (and thus does not require flexibility). In the third stage, the solution chooses to install only PEM to benefit from the flexibility the technology provides. Despite being approximately $2\%$ more expensive in investment, PEM technology becomes cheaper in the third stage due to lower yearly maintenance costs. Note that PEM provides cheaper maintenance in all stages. Nevertheless, AL still provides a lower sum of yearly investment and maintenance costs in the first and second stage compared to PEM.

Note that most of the production capacity is located in zone NO4 in northern Norway due to lower and less volatile electricity prices. Even though approximately $80\%$ of the demand is located in the area south of Trondheim, there are no facilities south of Florø in any of the strategic scenarios. On average, installed capacity in Trondheim and Florø corresponds to more than $50\%$ of the overall demand. This indicates that electricity costs, and hence production costs, are too high in zones NO1 and NO2 to install production capacity. These results further imply that it is cheaper to pay for distribution costs instead of producing with high electricity costs.

Fig. 6

Location and technology decisions for F17D30N9S100R1

Figure 6 shows that the solution chooses to combine technologies at all but one location in Glomfjord. In order to better understand and analyze the advantages of combining technologies at a location, we also consider the special case where mixing technologies at a single location is not allowed (single technology formulation) by imposing additional Constraints (29).

$$\begin{aligned} & \sum _{k \in \mathcal {K}}y_{ipk}^{n} + \sum _{p^{\prime } \in \mathcal {P}\setminus {\{p\}}}\sum _{k \in \mathcal {K}} y_{ipk}^{n^{\prime }} \le 1, \hfill \nonumber n \in \mathcal {N}, n^{\prime } \in \mathcal {A}_n, i \in \mathcal {I}, p \in \mathcal {P}. \end{aligned}$$

(29)

Using the Lagrangian relaxation with the two-step R-MIP, we obtain high-quality solutions for all instances with optimality gaps below $1.5\%$ for the single technology formulation as well. Similar to the general model formulation (1)–(17), the technology decisions are the same for all strategic scenarios and differ only in the installed capacity level.

When allowing only one technology per location, all selected locations south of Glomfjord use PEM technology due to its flexibility, see Fig. 7. Since PEM has equal operational costs with AL in the third stage, and lower maintenance costs, all locations, which aim to install large capacity (above capacity level 5) in the last stage, also use PEM. Only two locations, Finnsnes and Svolvær, with relatively small facilities (below capacity level 3) use AL technology.

Fig. 7

Location and technology decisions for F17D30N9S100R1 for the single technology model formulation

The objective function value of the general model solution is about $1\%$ lower than for the single technology solution despite installing more capacity. The expected operational costs show that the general model enables cheaper hydrogen production in the first stage since the solution installs the AL technology at all locations. In addition, the general model solution installs less capacity in the first stage. A similar situation can be observed in the second stage where, however, the general model solution already installs more capacity, see Table 7. On the other hand, the single technology solution results in slightly lower operational costs in the third stage since there is flexible PEM technology at almost all locations. Further, Table 7 shows that also in the third stage, the installed capacity is higher for the general model formulation.

Table 7 Expected capacity utilization

To compare operations in the different solutions, we show the expected capacity utilization over the operational scenarios in the third stage for the strategic scenario with the lowest demand in Fig. 8. In this scenario, the differences in the capacity utilization are more apparent than in other scenarios. In general, the utilization of newer facilities, i.e. facilities opened in the third stage, is higher than the capacity utilization of older facilities due to technology development leading to more efficient facilities.

For the general model formulation, Fig. 8 a shows that higher installed capacity in the third stage allows to not use the oldest AL facilities in Florø and Trondheim (highlighted with red circles). While these facilities were beneficial in the first stage, they are relatively expensive in operation in the third stage and less flexible than PEM technology. Figure 8b shows the expected capacity utilization for the single technology formulation, where all facilities in Trondheim, Hellesylt and Florø are approaching $100\%$ utilization independent of their opening stage. As a result of the capacity utilization, the option not to produce during peak hours is limited. The advantage of flexible PEM technology is exploited in Storekorsnes and Berlevåg in the north, which are both quite large facilities with relatively low capacity utilization.

Fig. 8

Expected capacity utilization over operational scenarios for F17D30N9S100R1

7 Conclusion

We have presented a multi-horizon stochastic facility location problem formulation motivated by the real-world problem of locating hydrogen production facilities in Norway. Since the problem becomes easily very large, we present the Lagrangian relaxation method to solve the problem. With our solution method, we can solve instances with up to 150 operational scenarios with a gap below $1.5 \%$. For instances with more than 100 operational scenarios, Gurobi is not able to solve the linear relaxation.

When solving the problem of locating hydrogen production facilities in Norway, the results indicate that the technology choice is affected not only by investment and production costs but also by technology flexibility. In areas with high and volatile electricity costs, the solution chooses to open the PEM technology due to its flexibility despite its higher costs. In later stages, the solution further prefers newer facilities due to their lower operational costs, while the utilization of facilities installed in the first stage is relatively low.

In future work, the model can be further extended by considering facility closing. Uncertainty in technology development can be captured and analyzed by also considering uncertainty in investment and distribution costs. If lower run times are needed, further reformulation and decomposition by scenario could be considered. In addition, our approach should be combined with other methods to solve instances with a higher number of scenarios, which would also allow for stability testing.

Appendices

Appendix A Alternative model formulation

We provide an alternative model formulation to (1)–(17) presented in Sect. 3 and compare the LP bound and Gurobi run times. In the alternative model formulation (A1)–(A17), variables $\varvec{x}$ are defined without the technology and capacity index. Hence, constraints () only specify which customers can be served from facility i, but do not serve as strengthening constraints. Further, we define variables $\varvec{q}$ without the capacity index. The alternative model is given as follows:

$$\begin{aligned} & \min \sum _{n \in \mathcal {N}} \pi _{n}\tau _{n} \left( \sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} C_{pk}^{a} y_{ipk}^{a} \right. \nonumber \\ & \left. \qquad +\sum _{s \in \mathcal {S}_{n}} \omega ^{n}_{s} \left[ \sum _{a \in \mathcal {A}_{n}} \sum _{i \in \mathcal {I}} \sum _{p \in \mathcal {P}} \sum _{r\in \mathcal {R}} \sum _{t \in \mathcal {T}} E_{isrt}^{n}F_{p}^{na} q_{ipsrt}^{na} \right. \right. \nonumber \\ & \left. \left. \qquad +\sum _{i \in \mathcal {I}} \sum _{j \in \mathcal {J}} \sum _{r\in \mathcal {R}} T_{ij}x_{ijsr}^{n} \right. \right. \nonumber \\ & \left. \left. \qquad +\sum _{j \in \mathcal {J}}\sum _{r\in \mathcal {R}} M^u u_{jsr}^{n} + \sum _{i \in \mathcal {I}} \sum _{r\in \mathcal {R}} M^o o_{isr}^{n}\right] \right) , \end{aligned}$$

(A1)

subject to:

$$\begin{aligned}&\sum _{p \in \mathcal {P}} \sum _{k \in \mathcal {K}} y_{ipk}^{n} \le 1, \hfill n \in \mathcal {N}, i \in \mathcal {I}, \end{aligned}$$

(A2)

$$\begin{aligned}&\sum _{i \in \mathcal {I}} x_{ijsr}^{n}+u_{jsr}^{n} = D_{jr}^{n}, \hfill n \in \mathcal {N}, j \in \mathcal {J}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(A3)

$$\begin{aligned}&x_{ijsr}^{n} \le L_{ij} D_{jr}^{n}, \hfill n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, j \in \mathcal {J}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(A4)

$$\begin{aligned} \sum _{j \in \mathcal {J}} x_{ijsr}^{n} + o_{isr}^{n} = \sum _{a \in \mathcal {A}_{n}} \sum _{p \in \mathcal {P}} \sum _{t \in \mathcal {T}} q_{ipsrt}^{na}, & \hfill n \in \mathcal {N}, i \in \mathcal {I}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(A5)

$$\begin{aligned} & q_{ipsrt}^{na} \le \sum _{k \in \mathcal {K}} K_{k}y_{ipk}^{a}, \hfill \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^C, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}, \end{aligned}$$

(A6)

$$\begin{aligned} & d_{ipksrt}^{na} \le y_{ipk}^{a}, \hfill \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}, \end{aligned}$$

(A7)

$$\begin{aligned} & q_{ipsrt}^{na} \le \sum _{k \in \mathcal {K}} K_{k}d_{ipksrt}^{na}, \hfill\nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S,s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}, \end{aligned}$$

(A8)

$$\begin{aligned} & q_{ipsrt}^{na} \ge \sum _{k \in \mathcal {K}} B_{p} K_{k} d_{ipksrt}^{na}, \hfill \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}, \end{aligned}$$

(A9)

$$\begin{aligned} & d_{ipksr(t-1)}^{na} - d_{ipksrt}^{na} \le 1 - \frac{1}{W_{p}} \sum _{t^{\prime } = t+1}^{t+W_p} d_{ipksrt^{\prime }}^{na}, \hfill \nonumber \\ &\quad n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r: t \le \overline{T}_r - W_p, \end{aligned}$$

(A10)

$$\begin{aligned} & d_{ipksr(t-1)}^{na} - d_{ipksrt}^{na} \le \frac{1}{W_{p}^{\prime }} \sum _{t^{\prime } = t-W_p^{\prime }}^{t-1} d_{ipksrt^{\prime }}^{na}, \hfill\nonumber \\ &\quad n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}_r: t \ge W_p^{\prime }, \end{aligned}$$

(A11)

$$\begin{aligned} & d_{ipksrt}^{na} \in \{0,1\}, \hfill \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I},p \in \mathcal {P}^S, k \in \mathcal {K}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in {T}, \end{aligned}$$

(A12)

$$\begin{aligned}&o_{isr}^{n} \ge 0, \hfill n \in \mathcal {N}, i \in \mathcal {I}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(A13)

$$\begin{aligned} & q_{ipsrt}^{na} \ge 0, \hfill \nonumber n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, t \in \mathcal {T}, \end{aligned}$$

(A14)

$$\begin{aligned}&u_{jsr}^{n} \ge 0, \hfill n \in \mathcal {N}, j \in \mathcal {J}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(A15)

$$\begin{aligned} & x_{ijsr}^{n} \ge 0, \hfill \nonumber n \in \mathcal {N}, i \in \mathcal {I}, j \in \mathcal {J}, s \in \mathcal {S}_{n}, r\in \mathcal {R}, \end{aligned}$$

(A16)

$$\begin{aligned}&y_{ipk}^{n} \in \{0,1\}, \hfill n \in \mathcal {N}, a \in \mathcal {A}_{n}, i \in \mathcal {I}, p \in \mathcal {P}, k \in \mathcal {K}, \end{aligned}$$

(A17)

We compare the computational results for different alternative model formulations to illustrate the effect of strengthening constraints presented in our model (1)–(17) (M1). Model (A1)–(A17) is referred to as M2. We further present results for model M3, where variables $\varvec{x}$ are again defined with the technology and capacity index (as in (1)–(17)), while variables $\varvec{q}$ remain without the capacity index. We also include strengthening constraints (4) in M3.

We compare the LP relaxation bound (LP bound), the best found objective (Best FS), the optimality gap (Gap) and the run time (Time). The time limit is set 24 hours for all instances. Table 8 shows that M1 provides the highest LP relaxation bound, which is on average $0.58\%$ and $1.12\%$ higher than M2 and M3, respectively. Further, M1 outperforms M2 and M3 in run time as well as in solution quality, providing always Best FS with the lowest objective. Note that the two largest instances F17D30N9S5R4 and F17D30N25S5R4 can be solved only using M1. Using M2 and M3 causes an out-of-memory error.

Table 8 Computational results for different model formulations

Appendix B Results from initial tests for LP+R-MIP heuristic

We provide initial results for different values of parameter $\Phi$ when using the LP+R-MIP approach. Parameter $\Phi$ serves as threshold value and variables $y \ge \Phi$ in the LP relaxation are set equal to one when solving the R-MIP.

For small instances, high values of $\Phi$ can be favourable because it leads to high-quality solutions, and longer run times are not an issue. However, for larger instances, the problem might be too big if the parameter $\Phi$ is chosen too high and only a few variables are fixed. Table 9 shows the optimality gap of the LP+R-MIP approach together with R-MIP run times for different $\Phi$ values. The run time limit for solving the R-MIP is set to 24 hours.

We choose $\Phi = 0.75$, which provides a good compromise between solution quality and run time.

Table 9 LP+R-MIP results for different $\Phi$ values

Appendix C Dimensions and complexity of problem instances

In Table 10, we report the number of constraints (m), number of binary variables (n01), number of continuous variables (nc), constraint matrix nonzero elements (nel), LP solution value (Z_LP), incumbent solution ($\overline{\text{ Z}_\textrm{BB}}$) and time to find this solution ($\overline{\text{ t}_\textrm{BB}}$), as well as the gap of the incumbent solution vs the LP bound (GAP_LP), and the gap of the incumbent solution vs the best lower bound from branch and bound (GAP_BB). All times reported by Gurobi are wall clock times. The optimality gap is calculated by Gurobi as: $\frac{UB-LB}{UB}$, where UB represents the upper bound and LB the lower bound. We have set a run time limit for Gurobi of 48 hours.

Note that Gurobi is not able to solve the LP relaxation for the three largest instances F17D30N9S100R1, F17D30N9S150R1, and F17D30N16S100R1 with more than 100 operational scenarios within the run time limit. The three instances have approximately 13.9, 20.8, and 23.3 million binary variables, respectively.

Table 10 Complexity of representative problem instances