<figure><img src="image_1.png" /></figure>

# Operational long-term management of a salt cavern for industry targeted green H₂ production

Robinson Beaucour, Nicolas Lamaison*

Univ. Grenoble Alpes, CEA, Liten, Campus Ines, 73375, Le Bourget du Lac France

## ARTICLE INFO

**Keywords:**

Hydrogen

Electrolyzer

Energy management system

MILP

Long-term

Salt cavern

Rolling-horizon

Carbon constraint

## ABSTRACT

Hydrogen and its derivatives such as ammonia are increasingly important for reducing carbon emissions in sectors with limited alternatives. Large production systems are generally preferred to reduce costs and for those, advanced management is a cornerstone to the project economic viability. The present work explores the development of an Energy Management System for an industrial low-emission ammonia production platform. The latter is connected to the electrical grid and consists of 1.5 GW of renewable production units, around 500 MW of electrolyzers and storage facilities including especially a salt cavern for H₂ storage. The long-term management of the salt cavern is decisive to enhance the profitability while satisfying industry related constraints such as the carbon content of the produced H₂. The literature lacks in addressing both aspects simultaneously. Methodological approaches to tackle the long-term management and industrial constraints are proposed and compared. For the considered platform, it is shown that a strategy based on precalculated cost functions in the long term horizon allows relative gains of 2.5% with respect to a state-of-the-art short term management leading to an annual cost saving of approximately $6.2 million in the platform operational expenses. This translates to a reduction of the gap between the state-of-the-art and the optimal *a posteriori* solution by about 82%.

## 1. Introduction

### 1.1. Context

The utilization of hydrogen stands as a pivotal step towards combatting climate change through the decarbonization of industrial processes, transportation and the seamless integration of variable renewable energy sources into the electricity grid [1]. Presently, there is a significant global upsurge in the demand for hydrogen, highlighting its increasing importance across various sectors [2]. Especially for industry, hydrogen can be converted into ammonia in order to be primarily used for fertilizers.

With this rise in demand, low-emission hydrogen production is expected to grow from 1 to about 20 Mtons/year [2] leading to a constant expansion of the pipeline of projects aiming at fostering low-emission hydrogen production. In order to reduce the specific production costs, gradually larger projects are announced with capacities already well exceeding 300 MW planned in Germany, China and Australia [3].

Concurrently, standards regarding the definition of low-emission hydrogen are being discussed, initiating the need to anticipate new production constraints for the project holders. For instance, the Green Hydrogen Organisation [4] proposes to define a standard requiring that

green hydrogen projects operate at less than 1 ton of CO₂ emitted per ton of H₂ produced taken as an average over a 12-months period.

In pursuit of enhancing the competitiveness and integration of low-emission hydrogen within the broader energy landscape, extensive research has been undertaken. This integration marks a paradigm shift in scientific understanding with the emergence of smart energy systems [5], which employ multiple energy vectors. Hydrogen can be integrated into such systems mainly through electrolysis processes to boost power-to-gas production while the waste heat from these processes can be valorized through district heating for e.g. [6]. The literature presents studies of numerous multi-vector systems. Gabrielli et al. [7] delved into the optimal design of energy systems encompassing both electricity and heat demands, thereby facilitating the integration of a substantial proportion of renewable electricity through hydrogen production and storage mechanisms. Haratyk and Forsberg [8] conducted an in-depth exploration of a nuclear-renewable energy system focusing on both hydrogen and electricity production, addressing particularly the utilization of electrolyzer partial and overload modes and hydrogen storage mechanisms to effectively balance the inherent variability of wind energy.

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: robinson.beaucour@cea.fr (R. Beaucour), nicolas.lamaison@cea.fr (N. Lamaison).</footnote>
Moreover, as planned for the announced projects [3], low-emission hydrogen production is generally associated to storage capacities and particularly large ones. The interest for large storage is also denoted in the literature in the field of power-to-gas-to-power. Gabrielli et al. [9] shed light on the immense potential of extensive hydrogen storage solutions, such as salt caverns, in mitigating CO₂ emissions from energy grids characterized by high variable renewable energy production shares. Clerjon and Perdu [10] uses an innovative signal decomposition approach [11] to show that seasonal electricity storage is particularly difficult to ensure and the only economically viable solution appears to be H₂ storage in salt caverns.

As stated beforehand, the future industrial green H₂ production plant will be large and probably associated to large salt cavern. The latter means that an optimized operational management, i.e. an Energy Management Systems (EMS), should allow for consequent financial savings. However, this EMS will face the following challenges:

* The horizon of prediction of renewable productions does not exceed few days while the system is equipped with a storage which capacity exceeds weeks of production. This time-scale mismatch leads to the necessity of adding long-term management capabilities to the EMS.

* Within the industrial context, the green H₂ must be produced respecting constraints such as CO₂ content [4] or rate of production to match the low flexibility of downstream processes such as Haber-Bosch ammonia production reactor [12].

In the present work, an EMS that accounts for long-term management and industrial constraints is developed in order to evaluate the financial savings obtained on a large-scale plant.

## 1.2. State-of-the-art

The existing literature offers a rich array of methodologies for constructing EMS [13]. As summarized by Alabi et al. [14], the two advanced EMS methodologies that can cope with the dynamic variations of smart energy systems are Model Predictive Control (MPC) and Reinforcement Learning (RL). The latter approaches are illustrated by Stoffel et al. [15], in the context of building energy management systems, who compared traditional reactive expert rules to model-free RL, MPC based on white-box, gray-box and black-box modeling and approximate MPC. While RL seems promising, there are still scientific challenges to overcome such as its lake of explicability, the effective management of specific constraints and the necessary consequent of-line training on precise simulators. Those are probably the reasons why the studies that have probed the domain of EMS tailored specifically to hydrogen production in industrial settings are based on MPC.

Noteworthy among these, Klyapovskiy et al. [16] showcased the considerable potential of EMS in increasing the economic performance of industrial facilities reliant on hydrogen. Their approach involved the utilization of hydrogen and electricity, partially sourced from local renewables, emphasizing the pivotal role of EMS in optimizing such intricate energy systems, all executed within the framework of a Mixed-Integer Linear Programming (MILP) model. Expanding on this, their approach assesses energy management across four representative days, providing a planning horizon with a renewable production forecast of 24 h. Such an approach however lacks the management of long term phenomena.

When addressing specifically long-term management with MPC, the methodological approaches can be classified into 2 categories: (i) the hierarchical approach and (ii) the continuous approach. For the hierarchical approach, the EMS uses 2 different models solved on different time horizons (short/long) and meshes (fine/coarse). It is called hierarchical because the long-term model outputs are used to set constraints and/or objectives to the short-term model. The continuous approach is however based on a single model with different strategies to tackle the long-term phenomena.

Among hierarchical approaches, Petrollese et al. [17] and Weber et al. [18] have explored EMS based on a combination of Optimal Generation Scheduling (OGS) and MPC. The OGS provides reference trajectories for energy storage levels, which the MPC aims to follow with incentives in the objective function. The reference trajectory is updated recursively when the energy system state changes.

In Petrollese et al. [17], the system is a renewable hydrogen-based microgrid. It includes various energy storage systems, and the OGS provides objectives for battery and hydrogen storage levels, taking into account uncertainties in weather and load forecasting. The MPC component handles the real-time management, adjusting actions to meet these objectives while responding to current conditions.

In Weber et al. [18], the EMS focuses on building climate control with seasonal energy storage, particularly a novel thermochemical storage. In this approach, the OGS sets energy storage targets, but rather than merely guiding the MPC, it imposes constraints that the MPC must respect. This distinction results in a different system dynamics where the MPC operates within stricter boundaries, leading to potentially longer computation times and infeasibilities.

Saletti et al. [19] proposed an approach that combines the long-term constraints of a multi-vector system with the short-term operational constraints. Their work introduces a three-module control system. The Distribution module oversees the final energy delivery to the end-user and provides energy demand forecasts using a dynamic programming approach. The Short-term supervisory controls the energy system's actions using a MILP model that takes into account the system state, forecasts, and long-term constraints. The time horizon is sufficiently short to obtain a quick solution. Finally, the Long-term supervisory establishes the long-term constraints through a linear programming problem on a daily basis. This hierarchical approach has also been recently used for multiple power-to-gas systems sharing a long term storage and connected to a district heating network [20].

Regarding the continuous approaches, Abomazid et al. [21] crafted an EMS specifically tailored to tackle challenges arising from non-constant electrolyzer efficiency and seasonal fluctuations in renewable energy production. Through the implementation of a rolling-horizon methodology, their innovative EMS adeptly managed dynamic variables, underscoring its adaptability and efficacy in diverse operational contexts. Notably, their approach introduced a Z-score factor to signal seasonal trends in electricity prices to the EMS, marking an improvement in the integration of long-term strategy within a short-term vision. Darivianakis et al. [22] and Thaler et al. [23] studied the integration of a stored energy value into a traditional EMS. Specifically, Thaler et al. [23] developed a hybrid MPC strategy for a renewable microgrid with seasonal hydrogen storage. This EMS assigns a cost value to the hydrogen stored, serving as an incentive for long-term energy shifts. The study highlighted that even a simple constant value for stored hydrogen significantly enhances EMS performance. Their method determines the optimal hydrogen storage value by evaluating various scenarios, emphasizing the balance between production and consumption throughout the year. Darivianakis et al. [22] introduced a more sophisticated approach, utilizing historical data to dynamically adjust the value of energy stored over time. This method involves constructing bounds around the optimal charging trajectory and developing a piece-wise affine approximation of the storage value function at each time step. The resulting multistage stochastic optimization problem aims to minimize the total energy consumption while ensuring that the system operates close to global optimality. This approach is particularly relevant for systems with significant seasonal variations.

Regarding EMS strategies including long term management but not in the field of hydrogen production, Cuisinier et al. [24,25] presented a rolling horizon approach based on a MILP model that combines the short term, evaluated hourly, and the long term, evaluated at an arbitrarily longer interval. The sole decision variable for the long term becomes the variation in the State of Energy (SOE) of the long-term storage, a thermal pit storage associated to district heating in their
<figure><img src="image_2.png" /><figcaption>Fig. 1. Scheme of the system.</figcaption></figure>

case. This variable is employed by Cost Functions (CF), which represent the operational cost of the system for a given time interval and for a given variation of storage SOE. These CF are precalculated based on numerous MILP models over Representative Periods (RP) of the long term horizon.

Based on the analysis of the state-of-the-art, none of the papers examined cope with industrial constraints such as maximal CO₂ content of the produced H₂ or minimal H₂ rate of production. The most recent and advanced approaches for EMS tailored specifically to hydrogen production in industrial settings are based on hierarchical approaches, that may have difficulties to cope with such constraints. In the field of district heating, the methodological continuous approach of Cuisinier et al. [24] is promising in order to both cope with the long term management of the salt cavern and these industrial constraints.

## 1.3. Main contributions and paper organization

The key contribution of the present paper is to propose two innovative long-term constrained strategies coping with industrial constraints for an EMS of a green hydrogen production chain:

* the *Multiple RP-CF long term* approach which is a stochastic extension of the approach from [24] adapted to a hydrogen-based energy system;

* the *Price signal* approach which uses a modified price of hydrogen calculated based on CF.

These 2 strategies are evaluated using a robust framework allowing to benchmark different EMS performance across diverse scenarios. In this framework, the operation of the system is optimized based on future boundary conditions predicted only for a short term horizon (up to 48 h). The innovative strategies are compared to (i) the long-term approach of Cuisinier et al. [24] and (ii) a state-of-the-art short term approach. Additionally, they are benchmarked against the optimal *a posteriori* solution obtained with a perfect knowledge of the boundary conditions over the entire year of operation.

The remainder of the paper proceeds as follows. The system investigated and its modeling are described in Section 2. The specific case study and the comparisons of the EMS performances are detailed in Section 3. Section 4 discusses the obtained results and methodological choices. Section 5 concludes the paper while giving perspectives to potential future work.

# 2. System investigated and methods

## 2.1. System description

The energy platform analyzed in the present paper is illustrated in Fig. 1. It is chosen to be representative of the large projects found in the database of IEA [3]. It consists in producing H₂ for an ammonia

production plant with limited connection to the carbon rich global electricity grid. The energy needed to produce H₂ mainly comes from the installed 1.5 GW of renewable (photovoltaic and windturbine). The local electricity grid supplies electrolysers and compressors. Converters enable renewable installations, electrolysers, and compressors to connect to the local electricity grid by adjusting the voltage. Additionally, the system uses a local H₂ grid. Both grid can be managed using storages (H₂ cavern and pipeline, electrical battery) and external flows (electrical grid and H₂ penalties). The latter is explained in Section 2.2.2. Finally, the system must respect specific industrial constraints of H₂ production rate and CO₂ intensity of the produced H₂ (see Section 2.2.3). As stated in Section 1, the management of the salt cavern and the industrial constraints is specifically investigated in the present paper.

## 2.2. Modeling approach

The platform is modeled using a MILP formulation. The nomenclature of parameters and decision variables are described in Table 1. It uses the index *sys* for each type of system, with the set of indices for each type of system being described in Table 2. Two renewable production systems are considered: photovoltaic (PV) and wind turbine (Wind). Five flow systems are considered: electricity consumed from global grid (ElecCons), electricity injected to the global grid (ElecInj), electricity curtailed from renewables (ElecCurt), H₂ sold from the local grid (H₂sold), H₂ bought or penalties (H₂pen). Three converters are used, namely for renewable (Ren), electrolysers (Ely), compressors (Comp). Three storages are used, namely battery, pipeline, salt cavern (Cavern). In the remainder of the paper, decision variables of the MILP problem are written in bold while binary variables are written in lower case.

The MILP modeling for the renewable production, grid, converters, compressors and storage is rather common and thus described in Appendix A. As detailed in Appendix A, the grid, converters, compressors and storage sets of equations are respectively denoted *F*, *C*, *K* and *S*.

### 2.2.1. Electrolyzer modeling

The Electrolyzer system is composed of numerous electrolysers modules ($N_{ely}$), themselves containing a set of stacks and their own balance of plant (heat exchanger, pumps and separators). Each module uses electricity for auxiliaries consumption and for H₂ conversion. The efficiency of the different modules of electrolysis is not linear. For the sake of simplicity, the wording electrolyzer is used to designate an electrolyzer module (i.e. set of stacks with balance of plant). As shown in Fig. 2, the electrolyzer operating range is considered as a concave linear piecewise function. ($L_i^{(k,elec)}$, $L_i^{(k,H_2)}$) are the coordinates of the vertex of the electrolyzer operating range where *k* is the identification number of each module. Thanks to the operating range concavity and the objective of maximization, additional binary variables are not
<table><caption>**Table 1**

Parameters and decision variables.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td colspan="3"><strong>(E) Renewable production</strong></td></tr><tr><td>E<sub>t</sub><sup>(sys, power)</sup></td><td>Capacity factor at t produced by sys</td><td>[MW/MW<sub>c</sub>]</td></tr><tr><td>E<sub>t</sub><sup>(sys, capacity)</sup></td><td>Capacity installed of sys</td><td>[MW]</td></tr><tr><td>E<sub>t</sub><sup>(sys, CO<sub>2</sub>)</sup></td><td>CO<sub>2</sub> intensity of sys</td><td>[tCO<sub>2</sub>/MWh]</td></tr><tr><td colspan="3"><strong>(F) Grid flow</strong></td></tr><tr><td>F<sub>t</sub><sup>(sys, prices)</sup></td><td>Price of sys flow at t</td><td>[SMWh] or [$/tH<sub>2</sub>]</td></tr><tr><td>F<sub>t</sub><sup>(sys, max/min)</sup></td><td>Max/Min flow of sys</td><td>[MW] or [tH<sub>2</sub>/h]</td></tr><tr><td>F<sub>t</sub><sup>(sys, sens)</sup></td><td>Sens of sys : 1 if input, -1 if output</td><td>[tCO<sub>2</sub>/MWh] or [tCO<sub>2</sub>/tH<sub>2</sub>]</td></tr><tr><td>F<sub>t</sub><sup>(sys, CO<sub>2</sub>)</sup></td><td>CO<sub>2</sub> intensity of sys</td><td>[tCO<sub>2</sub>/MWh] or [tCO<sub>2</sub>/tH<sub>2</sub>]</td></tr><tr><td>F<sub>t</sub><sup>(sys, flow)</sup></td><td>Flow of sys at t</td><td>[MW] or [tH<sub>2</sub>/h]</td></tr><tr><td>F<sub>t</sub><sup>(sys, cost)</sup></td><td>Cost of sys flow at t</td><td>[S]</td></tr><tr><td>F<sub>t</sub><sup>(sys, CO<sub>2</sub>)</sup></td><td>CO<sub>2</sub> emitted by sys flow at t</td><td>[tCO<sub>2</sub>]</td></tr><tr><td colspan="3"><strong>(C) Converter</strong></td></tr><tr><td>C<sup>(sys, eff)</sup></td><td>Conversion efficiency of sys</td><td>[-]</td></tr><tr><td>C<sup>(sys, max/min)</sup></td><td>Max/Min conversion of sys</td><td>[MW]</td></tr><tr><td>C<sub>t</sub><sup>(sys, in)</sup></td><td>Input flow of sys converter at t</td><td>[MW]</td></tr><tr><td>C<sub>t</sub><sup>(sys, out)</sup></td><td>Output flow of sys converter at t</td><td>[MW]</td></tr><tr><td colspan="3"><strong>(K) Compressor</strong></td></tr><tr><td>K<sup>(sys, coe f)</sup></td><td>Coefficient conversion of sys</td><td>[MW/(tH<sub>2</sub>/h)]</td></tr><tr><td>K<sup>(sys, max/min)</sup></td><td>Max/Min flow of sys</td><td>[MW]</td></tr><tr><td>K<sub>t</sub><sup>(sys, elec)</sup></td><td>Electricity of sys compressor at t</td><td>[MW]</td></tr><tr><td>K<sub>t</sub><sup>(sys, H<sub>2</sub>)</sup></td><td>H<sub>2</sub> flow through sys compressor at t</td><td>[MW]</td></tr><tr><td colspan="3"><strong>(S) Storage</strong></td></tr><tr><td>S<sup>(sys, Charge Eff)</sup></td><td>Charge efficiency of sys</td><td>[-]</td></tr><tr><td>S<sup>(sys, Disc. Eff)</sup></td><td>Discharge efficiency of sys</td><td>[-]</td></tr><tr><td>S<sup>(sys, Stat Losses)</sup></td><td>Static losses of sys</td><td>[-]</td></tr><tr><td>S<sup>(sys, Charge Max)</sup></td><td>Max charge of sys</td><td>[MW] or [tH<sub>2</sub>/h]</td></tr><tr><td>S<sup>(sys, Disc. Max)</sup></td><td>Max discharge of sys</td><td>[MW] or [tH<sub>2</sub>/h]</td></tr><tr><td>S<sup>(sys, Max/Min SOE)</sup></td><td>Max/MIN SOE of sys</td><td>[-]</td></tr><tr><td>S<sup>(sys, SOE ini)</sup></td><td>SOE initial of sys</td><td>[-]</td></tr><tr><td>S<sup>(sys, Capacity)</sup></td><td>Capacity of sys</td><td>[tH<sub>2</sub>] or [MWh]</td></tr><tr><td>S<sub>t</sub><sup>(sys, SOE)</sup></td><td>SOE of sys storage at t</td><td>[-]</td></tr><tr><td>S<sub>t</sub><sup>(sys, Charge Flow)</sup></td><td>Charge flow of sys storage at t</td><td>[MW]</td></tr><tr><td>S<sub>t</sub><sup>(sys, Disc. Flow)</sup></td><td>Discharge flow of sys storage at t</td><td>[MW]</td></tr><tr><td>S<sub>t</sub><sup>(sys, Flow)</sup></td><td>Flow of sys storage at t</td><td>[MW]</td></tr><tr><td>S<sub>t</sub><sup>(sys, state)</sup></td><td>state of sys storage at t : 0 if charging else 1</td><td></td></tr><tr><td colspan="3"><strong>(L) Electrolyzer</strong></td></tr><tr><td>L<sup>(k, max/min)</sup></td><td>Max/Min power of the kth electrolyzer</td><td>[MW]</td></tr><tr><td>L<sup>(k, Aux Base)</sup></td><td>Base electricity consumption of the kth electrolyzer</td><td>[MW]</td></tr><tr><td>L<sup>(k, Aux Marg)</sup></td><td>Marginal electricity consumption rate of the kth electrolyzer</td><td></td></tr><tr><td>L<sup>(k, elec)</sup></td><td>Electricity for ith operating point of the kth electrolyzer</td><td>[MW]</td></tr><tr><td>L<sub>i</sub><sup>(k, H<sub>2</sub>)</sup></td><td>H<sub>2</sub> flow for ith operating point of the kth electrolyzer</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>L<sub>i</sub><sup>(k, elec)</sup></td><td>Electricity consumption of the kth electrolyzer at t for H<sub>2</sub></td><td>[MW]</td></tr><tr><td>L<sub>i</sub><sup>(k, H<sub>2</sub>)</sup></td><td>H<sub>2</sub> flow production of the kth electrolyzer at t</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>L<sub>i</sub><sup>(k, total)</sup></td><td>Total electricity consumption of the kth electrolyzer at t</td><td>[MW]</td></tr><tr><td>L<sub>i</sub><sup>(k, state)</sup></td><td>State of the kth electrolyzer at t : 1 if turned on else 0</td><td>[MW]</td></tr><tr><td colspan="3"><strong>Global parameters</strong></td></tr><tr><td>M<sup>(CO<sub>2</sub>)</sup></td><td>Maximum hourly CO<sub>2</sub> intensity of H<sub>2</sub> sold</td><td>[tCO<sub>2</sub>/tH<sub>2</sub>]</td></tr><tr><td>n<sub>i</sub><sup>(H<sub>2</sub>)</sup></td><td>Minimum hourly H<sub>2</sub> sold quantity</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>Δ</td><td>Duration of rolling constraints</td><td>[h]</td></tr><tr><td>δ<sub>i</sub></td><td>Duration of instant t</td><td>[h]</td></tr><tr><td>N<sub>ely</sub></td><td>Number of electrolysers modules</td><td></td></tr></tbody></table>

required thus limiting the computational time. Only one binary $l_t^{(k, state)}$ for each electrolyzer state and each instant is needed in Eq. (L.1).

$$l_t^{(k, state)} \cdot L^{(k, min)} \leq l_t^{(k, elec)} \leq l_t^{(k, state)} \cdot L^{(k, max)} \quad (L.1)$$

For $i \in \{1, n_{ope}\}$, where $n_{ope}$ correspond to the number of vertex in the operation curve, the feasible operating points for each module of electrolysis are described with Eq. (L.2). The latter extend the feasible operating range to the area under the operating curve in Fig. 2. However, because the optimization approach used here tends to maximize electrolyzer efficiency, the chosen operating points should be located on the operational curve of Fig. 2. Appendix B verifies that statement on the obtained results.

$$L_t^{(k, H_2)} \leq L_t^{(k, elec)} \cdot \alpha_i^{(k)} - L_t^{(k, state)} \cdot \beta_i^{(k)} \quad (L.2)$$

$$\alpha_i^{(k)} = \frac{L_{i+1}^{(k, H_2)} - L_i^{(k, H_2)}}{L_{i+1}^{(k, elec)} - L_i^{(k, elec)}} \quad (L.2.1)$$

$$\beta_i^{(k)} = \frac{L_{i+1}^{(k, H_2)} \cdot L_i^{(k, elec)} - L_i^{(k, H_2)} \cdot L_{i+1}^{(k, elec)}}{L_{i+1}^{(k, elec)} - L_i^{(k, elec)}} \quad (L.2.2)$$

Eq. (L.3) aims to reduce calculation time by setting an order of usage among electrolysers which breaks the symmetries of the problem. These equations are tolerable because all the electrolyzers have the same efficiency curve in the present study. If different aging of the electrolyzers occur, these equations should be removed or adapted.

$$l_t^{(k+1, state)} \leq l_t^{(k, state)} \quad (L.3)$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Electrolyzer operating range.</figcaption></figure>

<table><caption>**Table 2**
System description.</caption><thead><tr><th>Type of system</th><th>System</th><th>Attribute</th></tr></thead><tbody><tr><td>Renewable production (E)</td><td>PV, wind</td><td>Power,capacity, CO<sub>2</sub></td></tr><tr><td>Grid flow (F)</td><td>ElecInj, ElecCons, ElecCurt, H<sub>2</sub>sold, H<sub>2</sub>pen</td><td>Prices, max, min, sens, CO<sub>2</sub> flow, cost</td></tr><tr><td>Converter (C)</td><td>Enr, Ely, Comp</td><td>eff, max, min in, out</td></tr><tr><td>Compressor (K)</td><td>Comp</td><td>coef, max, min elec, H<sub>2</sub></td></tr><tr><td>Storage (S)</td><td>Battery, pipeline, cavern</td><td>ChargeEff, Disc.Eff StatLosses, ChargeMax DiscMax, Max/MinSOE SOEini, Capacity SOE, ChargeFlow Disc.Flow, Flow state</td></tr><tr><td>Electrolyzer (L)</td><td>1, ..., N<sub>ely</sub></td><td>max, min, AuxBase, AuxMarg, elec H<sub>2</sub>, total, state</td></tr></tbody></table>

Eq. (L.4) described additional electricity consumption for the auxiliaries. The total amount of electricity used by the kth electrolyzer at t is described by $L_t^{(k, total)}$.

$$L_t^{(k, \text{total})} = L_t^{(k, \text{elec})} + L_t^{(k, \text{AuxBase})} + L_t^{(k, \text{AuxMarg})} \cdot L_t^{(k, \text{elec})} \quad (L.4)$$

### 2.2.2. Platform modeling

The present Section describes the modeling of the interaction between the components at the scale of the platform. Outputs from renewable production can be sent to the local grid through its converter or be curtailed as shown in Eq. (P.1). The local electricity grid balance is described by Eq. (P.2).

$$\sum_{sys} E_{t}^{(sys,power)} = F_{t}^{(ElecCurt.flow)} + C_{t}^{(enr,in)} \quad (P.1)$$

$$C_{t}^{(enr,out)} + F_{t}^{(ElecCons.flow)} + S_{t}^{(battery.Flow)} = C_{t}^{(ely,in)} + C_{t}^{(comp.in)} + F_{t}^{(ElecInj.flow)} \quad (P.2)$$

The electricity required by the $N_{ely}$ electrolysers go through the assigned converter (Eq. (P.3)). The electricity required by the compressors goes through its converter (Eq. (P.4)), H₂ production from the electrolyzer enters the compressor (Eq. (P.5)). Finally, the local H₂ grid

<figure><img src="image_4.png" /><figcaption>Fig. 3. EMS schematic.</figcaption></figure>

balance is described by Eq. (P.6).

$$\sum_{k=1}^{N_{ely}} L_{t}^{(k, total)} = C_{t}^{(ely,out)} \quad (P.3)$$

$$K_{t}^{(sys,elec)} = C_{t}^{(comp,out)} \quad (P.4)$$

$$\sum_{k=1}^{N_{ely}} L_{t}^{(k,H_{2})} = K_{t}^{(sys,H_{2})} \quad (P.5)$$

$$K_{t}^{(sys,H_{2})} + S_{t}^{(pipeline.Flow)} + S_{t}^{(cavern.Flow)} = F_{t}^{(H_{2}sold,flow)} \quad (P.6)$$

### 2.2.3. Specific constraints modeling

The platform is subjected to additional industrial constraints, which aim to ensure minimal H₂ sell and maximal CO₂ intensity of the produced H₂. The idea of the latter constraints is to set up an energy management that will allow the platform to export green labeled H₂. The specific constraints consist of imposing Δ-rolling average of H₂ sell and CO₂ intensity, where Δ is an arbitrary chosen timelength set to one week in this study. These specific constraints are detailed in Eqs. (SC.1) and (SC.2) where s is an arbitrary instant. Eq. (SC.1) corresponds to the minimal H₂ production constraint while Eq. (SC.2) corresponds to the maximal CO₂ intensity. The parameters $m^{(H₂)}$ and $M^{(CO₂)}$ are respectively the minimum hourly H₂ sold and the maximum hourly CO₂ intensity of H₂.

$$\sum_{s=\Delta \le t < s} F_{t}^{(H_{2}sold,flow)} + F_{t}^{(H_{2}pen,flow)} \geq \Delta \cdot m^{(H_{2})} \quad (SC.1)$$

$$\sum_{s=\Delta \le t < s} \sum_{sys} F_{t}^{(sys,CO_{2})} + \sum_{sys} E_{t}^{(sys,CO_{2})} \le M^{(CO_{2})} \cdot \sum_{s=\Delta \le t < s} F_{t}^{(H_{2}sold,flow)} + F_{t}^{(H_{2}pen,flow)} \quad (SC.2)$$

### 2.3. EMS comparison framework

This subsection describes the framework used to compare different management strategies of the platform. As shown in Fig. 3, an EMS is considered here as an algorithm using a decision tool. The EMS sends adequate information to its decision tool and the decision tool sends back a decision planning. The EMS executes partially or totally the decision planning. The EMS and the decision tool can call secondary tools if needed. A same EMS can use different decision tools and decision tools can be tuned with different secondary tools.

The cosimulation platform PEGASE [26] is used in the present framework. The latter features an efficient cosimulation engine, together with integrated solutions for defining both rule-based and MPC-based control strategies. The principle of rolling-horizon is native in this platform. The latter is illustrated in Fig. 4.

The simulation is composed of cycles operated in sequence at a user defined rate ('timeshift' on Fig. 4), which is set to 1 h for the remainder of the paper. Within each cycle a set of actions are performed
<figure><img src="image_5.png" /><figcaption>Fig. 4. PEGASE cosimulation platform and rolling horizon.</figcaption></figure>

in sequence. Especially, the real system (or more precisely the simulator of the real system in the present study) sends state return to the controller, which then sends back setpoints to the simulator. For each cycle, this is performed within a temporal window illustrated here with the parameter 'future'. For the next cycle, the 'future' window, is moved accordingly to the new current time. For each cycle, boundary conditions (demand, prices, renewable production) are predicted over the future window and fed to the simulator and controller. To ease the interpretation, it is worth mentioning that there is no model mismatch here, i.e. the simulator model is the same as the controller one.

## 2.4. EMS strategies

For the different approaches later presented, each EMS sends an hourly forecast of renewable production to its decision tool, the past decision over a Δ duration interval, the current time τ for the EMS and the current platform status. The length of the forecast corresponds to the length of the short-term horizon (SH).

### 2.4.1. Simple RP-CF long-term approach

This part describes an EMS referred here as Simple RP-CF long-term in which RP and CF stand respectively for Representative Periods and Cost Functions. For this EMS, the decision tool is a MILP problem minimizing cost on short-term (SH) and long-term (LH) horizons, as shown in Fig. 5. LH is composed of temporal blocks aggregating the time steps. The length of the blocks and the number of blocks to consider are arbitrary at this stage but the length is usually one week minimum.

The objective function (Eq. (OF)) of the MILP model is made up of short-term and long-term costs where the cost functions for an instant t, CF<sub>t</sub>, are modeled using SOS2 functions [27]. SH and LH are described with further details in Section 2.4.1.1. Long-term costs are made up using CF<sub>t</sub> which are precalculated affine piecewise function detailed in Sections 2.4.1.2 and 2.4.1.3. S<sub>t</sub><sup>(cavern,SOE)</sup> and S<sub>t</sub><sup>(cavern,Flow)</sup> are the only variables used in LH. The cavern SOE is linked between SH and LH thanks to Eq. (LH.1).

$$ \text{Min} \quad \underbrace{\sum_{t \in SH_{\tau}} \sum_{sys} \mathbf{F}_{t}^{(\text{sys,cost})}}_{\text{Short Term}} + \underbrace{\sum_{t \in LH_{\tau}} \text{CF}_{t} \mathbf{S}_{t}^{(\text{cavern,flow})}}_{\text{Long Term}} \quad (\text{OF}) $$

Subject to:

$$ \forall t \in SH_{\tau}, \text{Eqs. (F, C, K, S, L, P, SC)} \quad (\text{SH}) $$

$$ \forall t \in LH_{\tau}, S_{t}^{(\text{cavern,SOE})} = S_{t-1}^{(\text{cavern,SOE})} + \frac{S_{t}^{(\text{cavern,Flow})}}{S^{(\text{cavern,Capacity})}} \quad (\text{LH.1}) $$

$$ \forall t \in LH_{\tau}, S^{(\text{cavern,MinSOE})} \le S^{(\text{cavern,SOE})} \le S^{(\text{cavern,MaxSOE})} \quad (\text{LH.2}) $$

#### 2.4.1.1. Horizons.
SH<sub>τ</sub> and LH<sub>τ</sub> are the discretized SH and LH after the instant τ which is the current time for the EMS. They can be described as SH<sub>τ</sub> = [τ + 1, θ<sub>τ</sub>] and LH<sub>τ</sub> = [θ<sub>τ</sub> + 1, Θ<sub>τ</sub>]. Therefore, SH<sub>τ</sub> and LH<sub>τ</sub> are sets of instants where each instant has an associated duration δ<sub>t</sub>. The cardinality of the LH can vary, it is even possible that LH is an empty set. In the later case, the EMS strategy is referred as Myopic. Durations of the LH instant are also customizable, with potentially different lengths for each aggregated instant. LH is used to evaluate system cost for a relatively uncertain far future in order to take into account future tendencies in the short-term decisions.

#### 2.4.1.2. Representative periods.
Representative Periods (RP) are a set of time periods with attributed weights aiming to represent a larger time period named split. RP are selected using the algorithm of Poncelet et al. [28]. Given a split [t, t + δ], RP can be assigned to the split as a list of pairs (period, weight). This relates to ({[a<sup>(k)</sup>, a<sup>(k)</sup> + ε], ω<sup>(k)</sup>})<sub>k∈[1,K]</sub> where K is the number of RP, ε the length of RP and a<sup>(k)</sup> is an instant of [t, t + δ]. The number and length of RP for a split are arbitrarily chosen. The sum of the weights is equal to 1 as shown in Eq. (RP.1). The calculation of RP for a given split is further detailed and illustrated in Appendix C.

$$ \sum_{k=1}^{K} \omega_{k} = 1 \quad (\text{RP.1}) $$

In this paper, splits are set to have a constant length δ. The associated RP are calculated for several years. The split set for each year can be written {[t<sub>i</sub>, t<sub>i</sub> + δ]|t<sub>i</sub> = δ · i, 0 ≤ δ · i, δ · (i + 1) < 8760, i ∈ N}. Given that the split length is constant, the split set is represented as a set of instants T = {t<sub>i</sub>|t<sub>i</sub> = δ · i, 0 ≤ δ · i, δ · (i + 1) < 8760, i ∈ N}. Therefore, the RP for each split t<sub>i</sub> in each year Y are written as a triplet {(a<sup>(k)</sup>(t<sub>i</sub>, Y), ε, ω<sup>(k)</sup>(t<sub>i</sub>, Y))}k∈[1,K].

#### 2.4.1.3. Cost functions.
Cost Functions (CF) objective is to obtain a representative estimation of the system costs for different splits of the year and for different long-term storage behavior in order to give a long-term vision to the EMS. They are affine piecewise functions used to describe aggregate operational cost over a period. These CF are defined using a set of point denoted CF<sub>t</sub><sup>(δ,Y)</sup>(Λ) with t ∈ T = {t<sub>1</sub>, ..., t<sub>n</sub>} an ordered set of instants in a year, δ ∈ D a set of split durations, Y ∈ Y a set of years and Λ ∈ V = {Λ<sub>0</sub>, ..., Λ<sub>L</sub>} an ordered set of long term storage SOE variations (between −1 and 1), which cardinality is L − 1.

The CF<sub>t</sub><sup>(δ,Y)</sup>(Λ) are defined in Eq. (CF.1) as the weighted average result of the MILP problem described in the following. To ease the reading, a<sup>(k)</sup>(t, Y) and a<sup>(k)</sup>(t, Y) + ε have respectively been replaced by a<sup>(k)</sup> and b<sup>(k)</sup> in the equations. The MILP model used for the present EMS strategy is similar to the model of Section 2.2 (see Eq. (CF.2)). The specific constraints Equations (SC) are rewritten in Eqs. (CF.3) and (CF.4) to adapt to the length ε of the RP. The constraint on cavern SOE is adapted relatively to RP and split length in Eq. (CF.5).

$$ \dot{\text{CF}}_{t}^{(\delta, Y)}(\Lambda) = \frac{\delta}{\epsilon} \sum_{k=1}^{K} \omega^{(k)}(t, Y) \cdot \text{Min} \sum_{a^{(k)} \le t < b^{(k)}} \sum_{sys} \mathbf{F}_{t}^{(\text{sys,cost})} \quad (\text{CF.1}) $$
<figure><img src="image_6.png" /><figcaption>Fig. 5. Planning horizons.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Illustrations of obtained CF.</figcaption></figure>

Subject to:

$$ \forall k \in [1, K], a^{(k)} \le t < b^{(k)}, \text{Eqs. (F, C, K, S, L, P)} \quad (\text{CF.2}) $$

$$ \forall k \in [1, K], \quad \sum_{a(k) \le t < b(k)} (\mathbf{F}_{t}^{(\text{H}_2\text{sold,flow})} + \mathbf{F}_{t}^{(\text{H}_2\text{pen,flow})}) \ge \epsilon \cdot m^{(\text{H}_2)} \quad (\text{CF.3}) $$

$$ \begin{aligned} \forall k \in [1, K], \quad & \sum_{a(k) \le t < b(k)} \left( \sum_{sys} \mathbf{F}_{t}^{(\text{sys,CO}_2)} + \sum_{sys} E_{t}^{(\text{sys,CO}_2)} \right) \\ & \le M^{(\text{CO}_2)} \cdot \sum_{a(k) \le t < b(k)} (\mathbf{F}_{t}^{(\text{H}_2\text{sold,flow})} + \mathbf{F}_{t}^{(\text{H}_2\text{pen,flow})}) \end{aligned} \quad (\text{CF.4}) $$

$$ S_{b^{(k)}-1}^{(\text{cavern,SOE})} = S_{a^{(k)}}^{(\text{cavern,SOE})} + \frac{\epsilon}{\delta} \cdot \Lambda \quad (\text{CF.5}) $$

Fig. 6 illustrates the results obtained when the CF are calculated for 24 different splits of 2 weeks and all $\Lambda \in \mathcal{V} = \{\Lambda_0, \dots, \Lambda_L\}$ (storage variations) for a given year. With that Figure, the estimation of the system costs for different splits of the year and for different long-term storage behavior is obtained. Appendix D details the set of results for each point of the curves of Fig. 6. For a more comprehensive grasp, additional visual representations and analysis of cost function derivatives can be found in Appendix E.

The $\bar{C}_t^{(\delta,Y)}(\Lambda)$ are now defined for every $t \in \mathcal{T}$ and for every SOE variation $\Lambda \in \mathcal{V}$. First, it can be extended for every instant $t$ of the year using the interpolation described in Eq. (CF.6). Second, it can

be extended for all values of SOE variations using the interpolation described in Eq. (CF.7)

$$ t_i \le t \le t_{i+1}, \quad \bar{C}_t^{(\delta,Y)}(\Lambda) = \frac{t_{i+1} - t}{t_{i+1} - t_i} \cdot \bar{C}_t^{(\delta,Y)}(\Lambda) + \frac{t - t_i}{t_{i+1} - t_i} \cdot \bar{C}_t^{(\delta,Y)}(\Lambda) \quad (\text{CF.6}) $$

$$ \begin{aligned} 0 \le k < L - 1, \Lambda \in [\Lambda_k, \Lambda_{k+1}], & \quad \bar{C}_t^{(\delta,Y)}(\Lambda) = \frac{\Lambda_{k+1} - \Lambda}{\Lambda_{k+1} - \Lambda_k} \cdot \bar{C}_t^{(\delta,Y)}(\Lambda_k) \\ & + \frac{\Lambda - \Lambda_k}{\Lambda_{k+1} - \Lambda_k} \cdot \bar{C}_t^{(\delta,Y)}(\Lambda_{k+1}) \end{aligned} \quad (\text{CF.7}) $$

### 2.4.2. Multiple RP-CF long-term approach

This section proposes a strategy to extend the simple RP-CF long-term approach of Section 2.4.1 by using an equivalent deterministic model of a stochastic approach [29]. It is referred here as the Multiple RP-CF long-term approach. The decision tool of the Multiple RP-CF long-term approach is a similar MILP problem than for the single approach but with additional scenarios on the long horizon (LH). The MILP is described as follows.

$$ \text{Min} \quad \sum_{t \in SH_r} \sum_{sys} \mathbf{F}_{t}^{(\text{sys, cost})} + \sum_{Y \in \mathcal{Y}} p_Y \sum_{t \in LH_r} \mathbf{C}_t^{(Y)} S_{t,Y}^{(\text{cavern,flow})} \quad (\text{M-OF}) $$

Subject to:

$$ \forall t \in SH_r, \text{Eqs. (F, C, K, S, L, P, SC)} \quad (\text{M-SH}) $$

$$ \forall Y \in \mathcal{Y}, \forall t \in LH_r, S_{t,Y}^{(\text{cavern,SOE})} = S_{t-1,Y}^{(\text{cavern,SOE})} + \frac{S_{t,Y}^{(\text{cavern,Flow})}}{S_{t,Y}^{(\text{cavern,Capacity})}} \quad (\text{M-LH.1}) $$

$$ \forall Y \in \mathcal{Y}, \forall t \in LH_r, S^{(\text{cavern,MinSOE})} \le S_{t,Y}^{(\text{cavern,SOE})} \le S^{(\text{cavern,MaxSOE})} \quad (\text{M-LH.2}) $$

In the objective function of Eq. (M-OF), the probability $p_Y$ can be interpreted as the weight of the cost function of year $Y$ in the total cost function calculation. For the remainder of the paper, these probabilities are equally distributed. On the LH, the salt cavern storage is thus managed over many scenarios with $S_{t,Y}^{(\text{cavern,SOE})}$ and $S_{t,Y}^{(\text{cavern,Flow})}$ representing salt cavern management for each year $Y$. For all scenarios, LH is continuous with SH because $\forall Y \in \mathcal{Y}, S_{\theta_r, Y}^{(\text{cavern,SOE})} = S_{\theta_r}^{(\text{cavern,SOE})}$ in Eq. (M-LH.1). Eq. (M-SH) is completely identical to Eq. (SH).

In order to illustrate the different approaches described beforehand, Fig. 7 highlights the EMS vision at $\tau = 2020-01-15 00:00$ (vertical red line) on the cavern management for the Myopic (top), Single RP-CF long term (middle) and Multiple RP-CF long term (bottom) approaches. The part to the left of the red line corresponds to past actions over a week. Note that these are taken into account for future optimization. The part to the right of the red line corresponds to actions planned on future horizons. In Fig. 7(a), since the Myopic approach is used, the salt cavern
<figure><img src="image_8.png" /><figcaption>Fig. 7. Example of salt cavern management for different approaches, left part of the red line showing past actions, right part showing future actions depending of the approach. (a) Myopic approach with a short horizon size of 48 h. (b) Single RP-CF long-term approach with a short horizon size of 48 h and 3 times 2 weeks on the long horizon. (c) Multiple RP-CF long-term approach with a short horizon size of 48 h and multiple planning on the long horizon of 3 times 2 weeks.</figcaption></figure>

is managed only on SH. Consequently the right part of the red line is the same length as SH. In Fig. 7(b), the Single RP-CF long-term approach is used. Therefore, the salt cavern is hourly managed on SH then managed with aggregated time steps on LH. Consequently, the right part of the red line exhibits hourly changes on SH and longer changes on LH. In Fig. 7(c), the Multiple RP-CF long-term is used. Therefore, the salt cavern is hourly managed on SH but managed with multiple scenarios on LH. Consequently, the right part of the red line exhibits hourly changes on SH and multiple changes on LH (4 different future scenarios).

### 2.4.3. RP-CF price signal approach

This Section describes an innovative approach to build an EMS referred as RP-CF price signal. This approach consists in using a conventional Myopic approach (i.e. Simple RP-CF long-term with no long horizon) but with a modified selling price for H₂, referred as a proxy price. The CF are used to calculate this proxy price signal in order to account for a more distant future.

To do so, the CF derivatives $\frac{dCF_t^{(\delta,Y)}}{d\Lambda}$ are defined in Eq. (PS.1). In addition, the affine transformation function $\phi_\tau$ translating from $[\tau + 1, \theta_\tau + 1]$ to $[-1, 1]$ is also defined in Eq. (PS.2). With these new functions, the proxy H₂ prices are eventually defined using Eq. (PS.3) where $F_t^{(H_2\text{sell,prices})}$ is the standard prediction of the H₂ selling price. A maximum is used in Eq. (PS.3) in order to prevent too low prices preventing H₂ production.

$$0 \le k < L-1, \Lambda \in [\Lambda_k, \Lambda_{k+1}[, \frac{dCF_t^{(\delta,Y)}}{d\Lambda}(\Lambda) = \frac{CF_t^{(\delta,Y)}(\Lambda_{k+1}) - CF_t^{(\delta,Y)}(\Lambda_k)}{\Lambda_{k+1} - \Lambda_k} \quad (\text{PS.1})$$

$$\forall t \in [\tau + 1, \theta_{\tau} + 1], \phi_{\tau}(t) = \frac{t - (\tau + 1)}{\theta_{\tau} - \tau} - \frac{t - (\theta_{\tau} + 1)}{\tau - \theta_{\tau}} \quad (\text{PS.2})$$

$$\forall t \in SH_{\tau}, F_t^{(H_2\text{sell,prices})} = \max \left( F_{\text{standard}}^{(H_2\text{sell,prices})}, \frac{\sum_{Y \in \mathcal{Y}} p_Y \cdot \frac{dCF_t^{(\delta,Y)}}{d\Lambda} \phi_{\tau}(t)}{S^{(\text{Cavern, Capacity})}} \right) \quad (\text{PS.3})$$

With $F_t^{(H_2\text{sell,prices})}$ defined in Eq. (PS.3), the decision tool is returning an optimal solution to the following MILP problem.

$$\text{Min} \sum_{t \in SH_{\tau}} \sum_{\text{sys}} F_t^{(\text{sys,cost})} \quad (\text{PS-OF})$$

Subject to:

$$\forall t \in SH_{\tau}, \text{Eqs. (F, C, K, S, L, P, SC)} \quad (\text{PS-SH})$$

## 3. Results

### 3.1. Case study presentation

For the case study, the platform described in Section 2.1 is considered to be installed in Western Australia. Renewable production data from 2022 are used to test and compare different EMS approaches. Years 2018 to 2021 are used to precalculate CF (see 2.4.1). It is crucial to use different data for precalculation and testing to prevent bias. More details on the used time series are given in Appendix F.
<table><caption>**Table 3**
EMS tested.</caption><thead><tr><th>Name</th><th>Description</th><th>EMS strategy</th></tr></thead><tbody><tr><td>Myopic</td><td>No LH</td><td>Myopic</td></tr><tr><td>1 × 2w</td><td>LH of one 2-weeks long split</td><td>Multiple RP-CF long-term</td></tr><tr><td>2 × 2w</td><td>LH of two 2-weeks long splits</td><td>Multiple RP-CF long-term</td></tr><tr><td>6 × 2w</td><td>LH of six 2-weeks long splits</td><td>Multiple RP-CF long-term</td></tr><tr><td>9 × 2w</td><td>LH of nine 2-weeks long splits</td><td>Multiple RP-CF long-term</td></tr><tr><td>12 × 2w</td><td>LH of twelve 2-weeks long splits</td><td>Multiple RP-CF long-term</td></tr><tr><td>Price Signal</td><td>Price signal based on 2 weeks CF</td><td>Price Signal</td></tr><tr><td>A posteriori</td><td>Optimal a posteriori management</td><td>a posteriori</td></tr></tbody></table>

RP and CF are calculated using a set of instants $T = \{t_i|t_i = \delta \cdot i, \delta \cdot (i+1) < 8760, i \in \mathbb{N}\}$ for which $\delta$ is equal to two weeks and a set of variations $V = \{\Lambda_k|\Lambda_k = k \cdot \frac{2}{L-1} - 1, 0 \le k < L, k \in \mathbb{N}\}$ where $L$ is arbitrarily set to 15.

Seven different EMS are tested and compared with the optimal a posteriori management. The latter is equivalent to a Myopic approach with an horizon of 8760 h for which the boundary conditions are perfectly known. The reason it is called a posteriori is because it can only be performed once the year of operation is over. The different EMS are listed in Table 3.

EMS are tested with a perfect forecast on the short horizon (SH) and with the following identical initial conditions: initial storage SOE is set to 50%, past week H₂ hourly sell is set to 4.5 tH₂/h, past hourly CO₂ intensity is set to 3 tCO₂/tH₂. The latter values are also used respectively for m(H₂) and M(CO₂). The different parameters used in the case study are summarized in Table 4. The electrolyzer operating curve (see Fig. 2) was obtained based on a physical model of a commercial electrolyzer module. The latter physical model was performed using the Modelica-based library of Nepveu et al. [30].

The EMS performance is based on the final economic assessment. The latter is calculated using the amount of electricity and H₂ exchanged. Earnings are set to 4 k$/tH₂ for H₂ sell, 18.6 $/MWh for electricity sell. Costs are set to 12 k$/tH₂ for H₂ penalties, 18.6 $/MWh for electricity bought. Differential H₂ in salt cavern at the end of the year relatively to initial state is valorized at 4 k$/tH₂.

## 3.2. Comparison

Fig. 8 and Table 5 summarize the results obtained for the different EMS. Strategies based on RP-CF are significantly effective in reducing the penalty H₂pen. As shown in Fig. 8(a), Price signal used 9 times less H₂pen than Myopic and Multiple RP-CF long-term EMS used almost no H₂pen for long enough horizon. Nevertheless, once the penalties are avoided, there is no advantage in lengthening the horizon. As shown in Fig. 8(b), H₂ production does not improve for LH longer than 6 weeks leading no additional economical benefits. As shown in Table 5, economic performance is strongly correlated with negative electricity balance once penalties are avoided highlighting the economic benefits associated with consuming electricity to generate and subsequently sell H₂ while still satisfying the CO₂ constraints. As expected, the highest economical output is given by a posteriori, based on an higher H₂ production and a more significant electricity reinjection into the grid. For the considered platform, it is shown that the strategy based on Multiple RP-CF long-term allows relative gains of 2.5% with respect to the Myopic management leading to an annual cost saving of approximately $6.2 million in the platform's operational expenses. This translates to a reduction of the gap between the Myopic and the optimal a posteriori solution by 82.3%.

Strategies lead to different salt cavern management. As shown in Fig. 9, Myopic uses less cavern capacity storage with lower SOE than other strategies at every instant. Price signal and Myopic have similar cavern management during January, February, November and December. Even though Price signal does not lead to 100% SOE, this strategy

<table><caption>**Table 4**
Parameters value.</caption><thead><tr><th>Parameter</th><th>System</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>E<sup>t</sup><sub>t</sub>(PV power)</td><td>PV</td><td>a</td><td>[MW/MW<sub>c</sub>]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(W ind. power)</td><td>Wind</td><td>a</td><td>[MW/MW<sub>c</sub>]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(PV capacity)</td><td>PV</td><td>733</td><td>[MW]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(W ind. capacity)</td><td>Wind</td><td>843</td><td>[MW]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(PV CO<sub>2</sub>)</td><td>PV</td><td>0.029</td><td>[tCO<sub>2</sub>/MWh]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(W ind. CO<sub>2</sub>)</td><td>Wind</td><td>0.013</td><td>[tCO<sub>2</sub>/MWh]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecCons, prices)</td><td>ElecCons</td><td>18.6</td><td>[SMWh]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecInj, prices)</td><td>ElecInj</td><td>18.6</td><td>[SMWh]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecCurt, prices)</td><td>ElecCurt</td><td>0.0</td><td>[SMWh]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(H<sub>2</sub>sold, prices)</td><td>H<sub>2</sub>sold</td><td>4000</td><td>[$/tH<sub>2</sub>]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(H<sub>2</sub>pen, prices)</td><td>H<sub>2</sub>pen</td><td>12 000</td><td>[$/tH<sub>2</sub>]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(sys, min)</td><td>All system</td><td>0</td><td>[MW]/[tH<sub>2</sub>/h]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecCons, max)</td><td>ElecCons</td><td>250</td><td>[MW]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecInj, max)</td><td>ElecInj</td><td>250</td><td>[MW]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecCurt, max)</td><td>ElecCurt</td><td>1e8</td><td>[MW]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(H<sub>2</sub>sold, max)</td><td>H<sub>2</sub>sold</td><td>9</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(H<sub>2</sub>pen, max)</td><td>H<sub>2</sub>pen</td><td>1000</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecCons, sens)</td><td>ElecCons</td><td>1</td><td></td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecInj, sens)</td><td>ElecInj</td><td>-1</td><td></td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecCurt, sens)</td><td>ElecCurt</td><td>-1</td><td></td></tr><tr><td>E<sup>t</sup><sub>t</sub>(H<sub>2</sub>sold, sens)</td><td>H<sub>2</sub>sold</td><td>-1</td><td></td></tr><tr><td>E<sup>t</sup><sub>t</sub>(H<sub>2</sub>pen, sens)</td><td>H<sub>2</sub>pen</td><td>1</td><td></td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecCons, CO<sub>2</sub>)</td><td>ElecCons</td><td>0.565</td><td>[tCO<sub>2</sub>/MWh]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecInj, CO<sub>2</sub>)</td><td>ElecInj</td><td>0</td><td>[tCO<sub>2</sub>/MWh]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(ElecCurt, CO<sub>2</sub>)</td><td>ElecCurt</td><td>0</td><td>[tCO<sub>2</sub>/MWh]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(H<sub>2</sub>sold, CO<sub>2</sub>)</td><td>H<sub>2</sub>sold</td><td>0</td><td>[tCO<sub>2</sub>/tH<sub>2</sub>]</td></tr><tr><td>E<sup>t</sup><sub>t</sub>(H<sub>2</sub>pen, CO<sub>2</sub>)</td><td>H<sub>2</sub>pen</td><td>1</td><td>[tCO<sub>2</sub>/tH<sub>2</sub>]</td></tr><tr><td>C<sup>t</sup>(sys, eff)</td><td>All system</td><td>0.9</td><td>[-]</td></tr><tr><td>C<sup>t</sup>(sys, min)</td><td>All system</td><td>0</td><td>[MW]</td></tr><tr><td>C<sup>t</sup>(sys, max)</td><td>All system</td><td>10<sup>8</sup></td><td>[MW]</td></tr><tr><td>K<sup>t</sup>(sys, eff)</td><td>All system</td><td>0.663</td><td>[MW]/[tH<sub>2</sub>/h]</td></tr><tr><td>K<sup>t</sup>(sys, min)</td><td>All system</td><td>0</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>K<sup>t</sup>(sys, max)</td><td>All system</td><td>9</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>S<sup>t</sup>(Battery, Charge, E f f)</td><td>Battery</td><td>0.837</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Pipeline, Charge, E f f)</td><td>Pipeline</td><td>1</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Cavern, Charge, E f f)</td><td>Cavern</td><td>1</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Battery, Disc, E f f)</td><td>Battery</td><td>0.837</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Pipeline, Disc, E f f)</td><td>Pipeline</td><td>1</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Cavern, Disc, E f f)</td><td>Cavern</td><td>1</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Battery, Stat Losses)</td><td>Battery</td><td>2.5<sup>-3</sup></td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Pipeline, Stat Losses)</td><td>Pipeline</td><td>0</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Cavern, Stat Losses)</td><td>Cavern</td><td>0</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Battery, Charge Max)</td><td>Battery</td><td>500</td><td>[MW]</td></tr><tr><td>S<sup>t</sup>(Pipeline, Charge Max)</td><td>Pipeline</td><td>9</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>S<sup>t</sup>(Cavern, Charge Max)</td><td>Cavern</td><td>9</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>S<sup>t</sup>(Battery, Disc, Max)</td><td>Battery</td><td>500</td><td>[MW]</td></tr><tr><td>S<sup>t</sup>(Pipeline, Disc, Max)</td><td>Pipeline</td><td>9</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>S<sup>t</sup>(Cavern, Disc, Max)</td><td>Cavern</td><td>9</td><td>[tH<sub>2</sub>/h]</td></tr><tr><td>S<sup>t</sup>(sys, Max, SOE)</td><td>All system</td><td>1</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(sys, Min, SOE)</td><td>All system</td><td>0</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(sys, SOEini)</td><td>All system</td><td>b</td><td>[-]</td></tr><tr><td>S<sup>t</sup>(Battery, Capacity)</td><td>Battery</td><td>500</td><td>[MWh]</td></tr><tr><td>S<sup>t</sup>(Pipeline, Capacity)</td><td>Pipeline</td><td>18</td><td>[tH<sub>2</sub>]</td></tr><tr><td>S<sup>t</sup>(Cavern, Capacity)</td><td>Cavern</td><td>378</td><td>[tH<sub>2</sub>]</td></tr><tr><td>L<sup>t</sup>(k, min)</td><td>All system</td><td>5.902</td><td>[MW]</td></tr><tr><td>L<sup>t</sup>(k, max)</td><td>All system</td><td>59.02</td><td>[MW]</td></tr><tr><td>L<sup>t</sup>(k, Aux, Base)</td><td>All system</td><td>0.5</td><td>[MW]</td></tr><tr><td>L<sup>t</sup>(k, Aux, Marg)</td><td>All system</td><td>0.085</td><td></td></tr><tr><td>Δ</td><td></td><td>168</td><td>[b]</td></tr><tr><td>M<sup>t</sup>(CO<sub>2</sub>)</td><td></td><td>3</td><td>[tCO<sub>2</sub>/tH<sub>2</sub>]</td></tr><tr><td>m<sup>t</sup>(H<sub>2</sub>)</td><td></td><td>4.5</td><td>[tH<sub>2</sub>/h]</td></tr></tbody></table>

uses considerably more storage capacity than Myopic. The Multiple RP-CF long-term EMS conducts to use salt cavern at 100% SOE many times.

## 4. Discussion

In the following discussion, we explore the significance of various parameters and configurations in the EMS optimization process, with a

<footnote>a Time series depending on scenario.</footnote>

<footnote>b Case dependent.</footnote>
<caption>**Table 5**</caption>

Results. The relative and absolute gains are calculated with respect to the *Myopic* results. The gap reduction is calculated with respect to the maximum potential improvement, i.e. the difference between *Myopic* and a *posteriori* results.

<table><thead><tr><th></th><th>H<sub>2</sub> sold (kt)</th><th>H<sub>2</sub> produced (kt)</th><th>Electricity balance (GWh)</th><th>H<sub>2</sub> penalties (t)</th><th>Gain (M$)</th><th>Relative gain (%)</th><th>Absolute gain (M$)</th><th>Gap reduction (%)</th></tr></thead><tbody><tr><td>Myopic</td><td>63.77</td><td>63.58</td><td>-18.29</td><td>459.5</td><td>248.48</td><td>0.0</td><td>0.0</td><td>0.0</td></tr><tr><td>1 x 2w</td><td>63.95</td><td>63.79</td><td>-30.64</td><td>13.0</td><td>254.43</td><td>+2.39</td><td>+5.95</td><td>+78.5</td></tr><tr><td>2 x 2w</td><td>64.01</td><td>63.83</td><td>-33.93</td><td>0.0</td><td>254.70</td><td>+2.50</td><td>+6.22</td><td>+82.0</td></tr><tr><td>6 x 2w</td><td>64.04</td><td>63.85</td><td>-36.14</td><td>0.9</td><td>254.72</td><td>+2.51</td><td>+6.24</td><td>+82.3</td></tr><tr><td>9 x 2w</td><td>63.71</td><td>63.84</td><td>-34.36</td><td>0.0</td><td>254.71</td><td>+2.51</td><td>+6.24</td><td>+82.3</td></tr><tr><td>12 x 2w</td><td>63.73</td><td>63.84</td><td>-34.18</td><td>1.1</td><td>254.71</td><td>+2.51</td><td>+6.24</td><td>+82.3</td></tr><tr><td>Price signal</td><td>63.85</td><td>63.67</td><td>-24.83</td><td>74.4</td><td>253.31</td><td>+1.94</td><td>+4.83</td><td>+63.7</td></tr><tr><td>A posteriori</td><td>64.50</td><td>64.31</td><td>-62.82</td><td>0.0</td><td>256.06</td><td>+3.05</td><td>+7.58</td><td>+100.0</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison outputs from EMS experiments. (a) Comparison of $H_{2pen}$ quantities used, with a lower value indicating superior performance. (b) Comparison of $H_2$ production, different from selling, where higher values denote enhanced outcomes.</figcaption></figure>

specific focus on the utilization of RP and CF. RP might play a crucial role in capturing essential characteristics of the energy system's operation, enabling efficient decision-making. We examine the advantages of using RP, including their potential to improve EMS performance. Moreover, we assess the impact of parameterizing the long-term horizon and split configurations on EMS optimization, highlighting the importance of striking a balance to achieve optimal performance. Additionally, we contrast the multiple scenario approach with the single scenario approach. Furthermore, we discuss the advantages and disadvantages of the *Price signal* approach. Lastly, we touch upon the process of EMS

performance validation and its significance in ensuring the reliability and robustness of the optimization approaches.

### 4.1. Representative periods and cost functions

The CF calculation could be obtained by solving an optimization on the entire split. However, depending on the complexity of the system and the number of estimation required, the calculation time of an entire set of CF could become prohibitive. The use of RP permits to limit the time calculation of a set of CF by limiting the problem to solve in duration. This shorter time calculation could be useful for operational purposes. Many reasons can lead to the need of recalculating the CFs and thus the need to do it with a limited computational effort:

* A change in system nature because of long periods of maintenance, expansion of the system or components aging.

* After one year of operation, we may want to add this new scenario into the *Multiple RP-CF approach* and thus recalculate the CF for this new year.

* The exploration of new system architectures/sizing and their control for R&D purposes.

RP customization options, such as the selection method, length, and number of periods, provide flexibility in adapting to various energy system configurations and maintain performance. Although the EMS performance might not be highly sensitive to these choices, employing RP allows for the inclusion of representative time series data.

### 4.2. Long-term horizon

The long-term horizon in EMS optimization emerges as a main factor influencing system performance. Its parametrization, involving both the length of the horizon and the split size, plays a crucial role in achieving good results. Particularly, the length of the long-term horizon exhibits a substantial impact on EMS performance. The latter is tightly linked to the time constants associated to the system (storage size, time-integral constraints, etc.). In the present study, the relative low size of the long term storage leads to a relative small length of the LH. Configurations with longer horizons tend to lead to improved outcomes up to a limit, as they allow for better anticipation and adaptability to system dynamics. Moreover, at equivalent horizon lengths, the nature of the split – the arrangement of data subsets within the horizon – may have minimal influence on EMS performance. This observation indicates that the split arrangement can be flexibly customized without compromising the overall effectiveness of the optimization process.

### 4.3. Multiple scenarios

The multiple scenarios approach presents an intriguing trade-off in EMS optimization. While it is possible that this approach yields slightly inferior results compared to the single scenario approach for a fixed test year, it should exhibit significant potential over a larger set of test years, yielding to a better robustness. This highlights the
<figure><img src="image_10.png" /><figcaption>(b) monthly statistics using boxplot (top and bottom of the box representing 1st and 3rd quartiles while the horizontal line is the median, points are outliers)</figcaption></figure>

<caption>Fig. 9. Salt cavern management.</caption>

importance of considering the performance of the multiple scenarios approach over a diverse range of scenarios, which better captures the system complexities and uncertainties. Additionally, a higher number of scenarios may lead to a long-term signal that is too constant, potentially limiting the approach's adaptability to dynamic changes in the energy system. To address this, emphasizing extreme scenarios in the long-term signal could be beneficial. By strategically highlighting scenarios with extreme conditions, the approach becomes more responsive to unforeseen events and perturbations, contributing to improved long-term decision-making and enhancing overall EMS performance.

## 4.4. *Price signal approach*

The *Price signal* approach, while demonstrating more limited results compared to the RP-CF multiple scenarios approach, still represents a significant improvement over the *Myopic* approach. The *Price signal* approach offers a valuable alternative for achieving enhanced EMS performance relatively simple to set up. Although it may not match the performance of the multiple scenarios approach with RP-CF in all cases, the *Price signal* approach proves to be efficient and effective, especially in situations where a rapid and resource-efficient EMS is essential. Moreover, depending on different short-term horizon configurations, the *Price signal* approach has the potential to exhibit greater efficiency gains. Its ability to optimize EMS operations using minimal resources makes it an attractive option for real-time decision-making and situations where computational complexities need to be reduced.

## 4.5. *Comparison framework*

The evaluation of EMS performance is a critical aspect that requires careful consideration to ensure unbiased and reliable results. One potential bias lies in the limited number of years of data used for evaluation. To obtain a comprehensive understanding of EMS performance,

it is advisable to employ data from multiple test years. By doing so, the variance in EMS performance can be better assessed, allowing for a more robust evaluation of the optimization approaches. Additionally, it is essential to avoid using the same years for both secondary tool calculations (see Fig. 3) and performance testing to prevent bias in the evaluation. By using distorted data or different sets of years for these purposes, any potential bias can be mitigated. The relative performance gaps calculated in this study demonstrate significant improvements with the RP-CF and *Price signal* approaches, validating their potential to enhance EMS performance. Overall, the EMS performance evaluation framework should strive to ensure unbiased assessments and draw on multiple years of data to capture a comprehensive view of the optimization approaches effectiveness.

# 5. Conclusion and perspectives

This study has presented a comprehensive analysis of an industry-targeted green H₂ production plant through the development of a MILP optimization model. The model incorporates a variable efficiency for multiple independent electrolyzer modules and introduces an innovative long-term H₂ storage mechanism. Two pioneering EMS strategies tailored for industry-focused H₂ production systems were proposed, with a particular emphasis on a novel extension of the Representative Periods-Cost Functions (RP-CF) strategy of Cuisinier et al. [24,25] to a stochastic approach, integrating historical data for enhanced accuracy. The main findings of this study reveal a substantial advancement in the RP-CF strategy, demonstrating relative gains of 2.5% with respect to the *Myopic* management leading to an annual cost saving of approximately $6.2 million in the platform's operational expenses. This translates to a reduction of the gap between the *Myopic* and the optimal *a posteriori* solution by 82.3%. The innovative *Price signal* approach also leads to consequent savings ($4.8 million) while keeping the EMS setup simple.
This highlights the effectiveness of the RP-CF strategy in optimizing the operational performance of the H₂ production plant with large storage, particularly in the context of long-term management. However, certain limitations and avenues for future research have emerged from this study.

Firstly, the absence of predefined rules for determining the optimal long-term horizon poses a challenge, and future investigations could explore methodologies for establishing such guidelines. The resilience of the proposed strategies to variations in long-term storage parameters, including capacity, loss and other behavioral characteristics, remains unexplored and warrants further examination. Additionally, the adaptability of the EMS strategies to different energy system configurations is an area that could benefit from in-depth exploration. Understanding how these strategies perform in diverse contexts will enhance their applicability across a range of industrial settings. Furthermore, the study acknowledges that more operational constraints, such as annual maintenance activities, electrical grid day-ahead constraints or other annual constraints, have not been investigated. Future research efforts should address these constraints to provide a more comprehensive understanding of the strategies practical implementation and operational robustness. Another field for future research is the extension of the RP-CF methodology to multiple long term storages possibly with different energy vectors (e.g. heat and H₂ for Power-to-Gas system including heat recovery).

### CRediT authorship contribution statement

**Robinson Beaucour:** Writing – original draft, Visualization, Validation, Methodology. **Nicolas Lamaison:** Writing – review & editing, Validation, Supervision, Methodology, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

The authors do not have permission to share data.

### Declaration of Generative AI and AI-assisted technologies in the writing process

During the preparation of this work the authors used ChatGPT 3.5 in order to improve the language and readability with the appropriate disclosure. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

### Acknowledgments

The authors express their gratitude to the researchers at the Centre for Applied Mathematics of Mines Paristech: Welington Luis de Oliveira, Sophie Demassey, Nadia Maizi, and Gilles Guerassimoff for their valuable support. The authors are grateful to Electricity Maps for providing power mix data including wind and PV production in Western Australia for the years 2018 to 2022 at an hourly time step.

## Appendix A. MILP models

The present Appendix describes the MILP modeling for the renewable production, grid, converters, compressors and storage.

### A.1. Renewable production

PV and wind production are modeled with capacity factor time series considered as inputs for the optimization problem. Thus, no decision variables are related to renewable production.

### A.2. Grid flow

Local electricity grid can consume (related to ElecCons) or inject (related to ElecInj) electricity from/to the global electricity grid. The local H₂ grid can sell H₂ (related to H₂sold). If the platform does not satisfy production or CO₂ constraints (see Section 2.2.3), it can “buy” expensive external green H₂ (related to H₂pen). H₂pen can be interpreted as a penalty to prevent unfeasibility during the platform operation. Eqs. (F.1) and (F.2) set hourly cost and CO₂ emissions for each flow, system and instant. Eq. (F.3) sets flow constraints related to infrastructures limitation.

$$F_{t}^{(\text{sys,cost})} = F^{(\text{sys,sens})} \cdot F_{t}^{(\text{sys,prices})} \cdot F_{t}^{(\text{sys,flow})} \quad (\text{F.1})$$

$$F_{t}^{(\text{sys,CO2})} = F^{(\text{sys,CO2})} \cdot F_{t}^{(\text{sys,flow})} \quad (\text{F.2})$$

$$F^{(\text{sys,min})} \le F_{t}^{(\text{sys,flow})} \le F^{(\text{sys,max})} \quad (\text{F.3})$$

### A.3. Converters

In order to use local electricity grid, entities must change voltage through converter. Converters imply some losses as described in Eq. (C.1). Infrastructure limitations are modeled with Eq. (C.2).

$$C_{t}^{(\text{sys,out})} = C^{(\text{sys,eff})} \cdot C_{t}^{(\text{sys,in})} \quad (\text{C.1})$$

$$C^{(\text{sys,min})} \le C_{t}^{(\text{sys,in})} \le C^{(\text{sys,max})} \quad (\text{C.2})$$

### A.4. Compressors

H₂ production from electrolyzer is compressed before storing and/or selling. Compressors uses electricity from local grid proportionally to H₂ flow (Eq. (K.1)) with infrastructure limitation (Eq. (K.2)).

$$K_{t}^{(\text{sys,elec})} = K^{(\text{sys,coef})} \cdot K_{t}^{(\text{sys,H}_2)} \quad (\text{K.1})$$

$$K^{(\text{sys,min})} \le K_{t}^{(\text{sys,H}_2)} \le K^{(\text{sys,max})} \quad (\text{K.2})$$

### A.5. Storage

The platform includes electrical storage (battery) and H₂ storage (pipeline and cavern). A binary variable is introduced for each storage to prevent simultaneously charging and discharging. An initial state of energy (SOE) is set depending on the case study. Storages are modeled with the constraints in Eqs. (S.1) to (S.6) (see Box I).

## Appendix B. MILP tolerance

In this Appendix, the impact of the modeling of the electrolysers with SOS2 variables is evaluated. The objective is to check that the operating points of the electrolysers on the obtained results do not deviate from the operating curve of Fig. 2. The gap between operating point given by the EMS and operating curve is calculated with Eq. (L.2). Fig. B.10 provides a visual representation of the gap, demonstrating that only a few operating points are significantly outside the allowed range.

$$\text{Gap}_t = L_t^{(\text{k,elec})} \cdot \alpha_t^{(k)} - L_t^{(\text{k,state})} \cdot \beta_t^{(k)} - L_t^{(\text{k,H}_2)}$$
$$ \forall t > 1, S_t^{(\text{sys,SOE})} - S_{t-1}^{(\text{sys,SOE})} = S^{(\text{sys,StatLosses})} \cdot S_t^{(\text{sys,SOE})} + \frac{S_t^{(\text{sys,ChargeFlow})} \cdot S^{(\text{sys,ChargeEff})} - \frac{S_t^{(\text{sys,Disc.Flow})}}{S^{(\text{sys,Disc.Eff})}}}{S^{(\text{sys,Capacity})}} \quad (S.1) $$

$$ \text{for } t = 1, S_t^{(\text{sys,SOE})} - S^{(\text{sys,SOEini})} = S^{(\text{sys,StatLosses})} \cdot S_t^{(\text{sys,SOE})} + \frac{S_t^{(\text{sys,ChargeFlow})} \cdot S^{(\text{sys,ChargeEff})} - \frac{S_t^{(\text{sys,Disc.Flow})}}{S^{(\text{sys,Disc.Eff})}}}{S^{(\text{sys,Capacity})}} \quad (S.2) $$

$$ 0 \le S_t^{(\text{sys,ChargeFlow})} \le S^{(\text{sys,ChargeMax})} \cdot S_t^{(\text{sys,state})} \quad (S.3) $$

$$ 0 \le S_t^{(\text{sys,Disc.Flow})} \le S^{(\text{sys,Disc.Max})} \cdot (1 - S_t^{(\text{sys,state})}) \quad (S.4) $$

$$ S_t^{(\text{sys.Flow})} = S_t^{(\text{sys.ChargeFlow})} - S_t^{(\text{sys.Disc.Flow})} \quad (S.5) $$

$$ S^{(\text{sys,MinSOE})} \le S_t^{(\text{sys,SOE})} \le S^{(\text{sys,MaxSOE})} \quad (S.6) $$

Box I.

<figure><img src="image_11.png" /><figcaption>Fig. B.10. Absolute gap between operation range of electrolyzer and operating point of the EMS.</figcaption></figure>

Despite minor deviations, the solutions remain close to the operating curve, attesting to the effectiveness and reliability of the MILP-based optimization approaches in achieving near-optimal EMS operation while adhering to constraints.

It is worth mentioning here that for all the results obtained, the solver used is CPLEX 22.1.0 and the relative gap between the solution returned and the optimal solution is notably small, with a value of 0.01 for both the Myopic and Price Signal approaches and an even smaller value of 0.001 for the Multiple RP-CF approach. The reduced tolerance for RP-CF can be attributed to the consideration of long-term cost, requiring higher precision in the optimization process.

## Appendix C. Representatives periods calculation

Fig. C.11 provides insights into the calculation process of RP and the steps involved in their determination. The RP calculation is based on the work of Poncelet et al. [28]. The first step encompasses the extraction of candidate RP from a segmented time series. Fig. C.12(a) illustrates this procedure, exemplified by the division of a time series spanning from 2020-02-12 to 2020-02-26 into seven distinct periods, each spanning two days. Subsequently, the second step involves the selection of RP, accompanied by the assignment of individual weights

to each RP. This process is visually depicted in Fig. C.12(a), where the chosen RP are highlighted in green, and the corresponding weight for each RP is indicated. The selection of RP is strategically based on the replication of the original duration curve through the new RP, as shown in Fig. C.12(b). This intricate methodology ensures that the chosen RP correctly capture the essential characteristics of the original time series.

## Appendix D. Cost function focus

Within this Appendix, the cost function computation process is detailed, providing a comprehensive understanding of its underlying mechanisms. To illustrate this, we compare the calculation outcomes over the same segmented time split, albeit with distinct storage variations, i.e. one with a +14% storage charge setting and the other with a -6% storage discharge setting. Figs. D.13 and D.14 showcase these calculations through a series of subfigures.

Figs. D.13(a) and D.14(a) graphically present the hourly SOE within each RP, highlighting the temporal fluctuations corresponding to the imposed variations of +14% and -6%. It is worth noting here that both RP are subjected to this imposed variation. Compared to a variant (not shown here) where the imposed variation is shared in a weighted fashion by the different RP, the present one allows enhancing differences in costs for the different storage variation settings leading to better results of the EMS. Figs. D.13(b) and D.14(b) show insights into the hourly H₂ production of each electrolyzer within each RP. Notably, distinct patterns emerge, revealing variations in production rates between the subfigures. An interesting observation is the synchronization of electrolyzers to optimize yield, emphasizing the intricate interplay between energy generation and storage.

Figs. D.13(c) and D.14(c) also provide a view of H₂ sales within each RP, with Fig. 2.c notably demonstrating the requirement for H₂pen to adhere to the minimum H₂ sales constraint. In a similar vein, Figs. D.13(d) and D.14(d) show the CO₂ content of the H₂ sold across each RP. These subfigures underscore the utilization of a cost-effective, albeit CO₂-intense, electricity grid to facilitate H₂ production while adhering to stringent CO₂ content constraints.

Lastly, Figs. D.13(e) and D.14(e) present the gains achieved through H₂ and electricity transactions within each RP. Gains are higher on Fig. D.13(e) than on Fig. D.14(e) as a results of different cavern storage constraint. The weighted mean of these gains, utilizing RP weights, is the base metric for the calculation of the cost function value.
<figure><img src="image_12.png" /><figcaption>Fig. C.11. Schematic of RP selection principles.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. C.12. Process to select RP. (a) Original time series sliced in RP candidates, RP candidates selected to be RP for the period are highlighted in green. (b) Original DC and DC from RP (DC approx), RP are selected to maximize similarity of DC from RP with original DC.</figcaption></figure>

## Appendix E. Cost function results and derivative

In this present Appendix, a nuanced analysis of the CF values and their derivatives serves as a pivotal aspect of comprehending the optimization landscape. Derivatives are calculated relatively to variation in tH₂ and not in percent of SOE, as described below:

$$ \frac{dCF_t^{(Y)}}{d\Lambda}(\Lambda) = \frac{1}{S^{(Cavern, Capacity)}} \frac{CF_t^{(Y)}(\Lambda_{k+1}) - CF_t^{(Y)}(\Lambda_k)}{\Lambda_{k+1} - \Lambda_k} $$

Splitting each year into 24 two-week segments, a total of 24 CF are generated for each year. Illustrated in Figs. E.15(a) and E.15(b), the distribution of CF values presents intriguing variations between the two distinct years, offering insights into the dynamic energy system behaviors.

Figs. E.15(c) and E.15(d) provide an in-depth exploration of the distribution of derivatives, shedding light on the complex relationship

between imposed storage variations and resulting CF values. While the CF exhibit convex-like behaviors, the derivatives unveil instances of non-convexity. Notably, some CF experience transition from decreasing to increasing trend.

A particularly informative pattern emerges from the derivatives, as most values tend to cluster around -4k $/tH₂. This observation signifies that, under specific conditions, a reduction of 1 tH₂ in the imposed storage variation leads to the potential sale of 1 tH₂ for 4k$. A derivative falling within the range of -4k to -16k $/tH₂ suggests that a storage variation reduction of 1 tH₂ enables the avoidance of H₂ penalties and the sale of H₂ in a given RP. On the other hand, a derivative ranging from -4k to -2k $/tH₂ indicates that a similar reduction in storage variation does not necessarily facilitate increased H₂ sales, but rather encourages the sale of more electricity, a less profitable path.
<figure><img src="image_14.png" /><figcaption>**Fig. D.13.** Illustration of the outputs of the calculation of a CF value where the storage variation is imposed by +14% in each RP. (a) Hourly cavern SOE variation rising in each RP to respect constraint. (b) Hourly H₂ production, due to large renewable output on RP1 H₂ production is almost maximum. (c) H₂ sell, H₂pen is needed on RP0 to respect minimum sell constraint. (d) CO₂ constraints. (e) Gains on each RP, on RP0 H₂pen implies an import cost, on RP1 high renewable production conducts to electricity injection.</figcaption></figure>

## Appendix F. Data

In this Appendix, detailed information about the data used in this paper, sourced from ElectricityMaps [31], are provided. The dataset includes capacity factor time series for photovoltaic and wind production in the region of Western Australia. Fig. F.16(a) displays the duration curve of the time series for each of the 5 years, offering insights into the

distribution and variability of PV and wind production. Furthermore, Fig. F.16(b) illustrates the rolling 2-week average of the time series for each year, allowing for a visual representation of the trends and patterns in PV and wind production over time. The data provided by ElectricityMap serves as a valuable foundation for the analysis and evaluation of various EMS optimization approaches presented in the core of the paper.
<figure><img src="image_15.png" /><figcaption>(a) Cavern management</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) $H_2$ constraints</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Electrolyser</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) $CO_2$ constraints</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(e) Gain</figcaption></figure>

**Fig. D.14.** Illustration of the outputs of the calculation of a CF value where the storage variation is imposed by -6% in each RP. (a) Hourly cavern SOE variation falling in each RP to respect constraint. (b) Hourly $H_2$ production, production on RP1 is not maximum due to use of cavern $H_2$ (c) $H_2$ sell, No more $H_{2pen}$ is needed on RP0 as a result of different constraints on cavern storage. (d) $CO_2$ constraints. (e) Gains on each RP.
<figure><img src="image_20.png" /><figcaption>Fig. E.15. Cost functions.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. F.16. Renewable production time series. (a) Duration curves showcase the temporal distribution of renewable energy generation for the years spanning 2018 to 2022. (b) The rolling two-weeks mean plot of PV and wind production from 2018 to 2022 offers a dynamic perspective on the average energy generation trends over time.</figcaption></figure>

## References

[1] Shukla PR, Skea J, Slade R, Al Khourdajie A, Van Diemen R, McCollum D, Pathak M, Some S, Vyas P, Fradera R, et al. Climate change 2022: Mitigation of climate change. In: Contribution of working group III to the sixth assessment report of the intergovernmental panel on climate change. vol. 10, Cambridge University Press Cambridge, UK; 2022, 9781009157926.

[2] IEA. Global hydrogen review 2022. Technical report, Paris: IEA; 2022, URL: https://www.iea.org/reports/global-hydrogen-review-2022.

[3] IEA. Hydrogen production and infrastructure projects database - Data product. 2023, URL: https://www.iea.org/data-and-statistics/data-product/hydrogen-production-and-infrastructure-projects-database.

[4] Green Hydrogen Organisation. Green hydrogen standard - The global standard for green hydrogen and green hydrogen derivatives including green ammonia. 2023, URL: https://gh2.org/sites/default/files/2023-01/GH2_Standard_A5_JAN%202023_1.pdf.

[5] Lund H, Østergaard PA, Connolly D, Mathiesen BV. Smart energy and smart energy systems. Energy 2017;137:556-65. http://dx.doi.org/10.1016/j.energy.2017.05.123, URL: https://www.sciencedirect.com/science/article/pii/S0360544217308812.

[6] Guelpa E, Bischi A, Verda V, Chertkov M, Lund H. Towards future infrastructures for sustainable multi-energy systems: A review. Energy 2019;184:2-21. http://dx.doi.org/10.1016/j.energy.2019.05.057, URL: https://linkinghub.elsevier.com/retrieve/pii/S0360544219309260.

[7] Gabrielli P, Gazzani M, Martelli E, Mazzotti M. Optimal design of multi-energy systems with seasonal storage. Appl Energy 2018;219:408-24. http://dx.doi.org/10.1016/j.apenergy.2017.07.142, URL: https://linkinghub.elsevier.com/retrieve/pii/S03606261917310139.

[8] Haratyk G, Forsberg CW. Nuclear-renewables energy system for hydrogen and electricity production. Nucl Technol 2012;178(1):66-82. http://dx.doi.org/10.13182/NT12-A13548, arXiv:https://doi.org/10.13182/NT12-A13548.

[9] Gabrielli P, Poluzzi A, Kramer GJ, Spiers C, Mazzotti M, Gazzani M. Seasonal energy storage for zero-emissions multi-energy systems via underground hydrogen storage. Renew Sustain Energy Rev 2020;121:109629. http://dx.doi.org/10.1016/j.rser.2019.109629, URL: https://linkinghub.elsevier.com/retrieve/pii/S1364032119308366.

[10] Clerjon A, Perdu F. Matching intermittency and electricity storage characteristics through time scale analysis: an energy return on investment comparison. Energy Environ Sci 2019;12(2):693-705. http://dx.doi.org/10.1039/C8EE01940A, URL: https://pubs.rsc.org/en/content/article/landing/2019/ee/c8ee01940a, Publisher: The Royal Society of Chemistry.

[11] Clerjon A, Perdu F. Matching intermittent electricity supply and demand with electricity storage - An optimization based on a time scale analysis. Energy 2022;241:122799. http://dx.doi.org/10.1016/j.energy.2021.122799, URL: https://www.sciencedirect.com/science/article/pii/S0360544221030486.

[12] Shamiri A, Aliabadi N. Modeling and performance improvement of an industrial ammonia synthesis reactor. Chem Eng J Adv 2021;8:100177. http://dx.doi.org/10.1016/j.ceja.2021.100177, URL: https://linkinghub.elsevier.com/retrieve/pii/S2666821121000922.

[13] Weitzel T, Glock CH. Energy management for stationary electric energy storage systems: A systematic literature review. European J Oper Res 2018;264(2):582-606. http://dx.doi.org/10.1016/j.ejor.2017.06.052, URL: https://linkinghub.elsevier.com/retrieve/pii/S0377221717305933.

[14] Alabi TM, Agbajor FD, Yang Z, Lu L, Ogungbile AJ. Strategic potential of multi-energy system towards carbon neutrality: A forward-looking overview. Energy Built Environ 2023;4(6):689-708. http://dx.doi.org/10.1016/j.enbenv.2022.06.007, URL: https://linkinghub.elsevier.com/retrieve/pii/S2666123322000435.

[15] Stoffel P, Maier L, Kümpel A, Schreiber T, Müller D. Evaluation of advanced control strategies for building energy systems. Energy Build 2023;280:112709. http://dx.doi.org/10.1016/j.enbuild.2022.112709, URL: https://www.sciencedirect.com/science/article/pii/S0378778822008805.

[16] Klyapovskiy S, Zheng Y, You S, Bindner HW. Optimal operation of the hydrogen-based energy management system with P2X demand response and ammonia plant. Appl Energy 2021;304:117559. http://dx.doi.org/10.13182/NT12-A13548, arXiv:https://doi.org/10.13182/NT12-A13548.
1016/j.apenergy.2021.117559, URL: https://linkinghub.elsevier.com/retrieve/pii/S0306261921009375.

[17] Petrollese M, Valverde L, Cocco D, Cau G, Guerra J. Real-time integration of optimal generation scheduling with MPC for the energy management of a renewable hydrogen-based microgrid. Appl Energy 2016;166:96-106. http://dx.doi.org/10.1016/j.apenergy.2016.01.014, URL: https://linkinghub.elsevier.com/retrieve/pii/S0306261916000337.

[18] Weber SO, Oei M, Linder M, Böhm M, Leistner P, Sawodny O. Model predictive approaches for cost-efficient building climate control with seasonal energy storage. Energy Build 2022;270:112285. http://dx.doi.org/10.1016/j.enbuild.2022.112285, URL: https://linkinghub.elsevier.com/retrieve/pii/S037877882200456X.

[19] Saletti C, Morini M, Gambarotta A. Smart management of integrated energy systems through co-optimization with long and short horizons. Energy 2022;250:123748. http://dx.doi.org/10.1016/j.energy.2022.123748, URL: https://www.sciencedirect.com/science/article/pii/S036054422200651X.

[20] Marzi E, Morini M, Saletti C, Gambarotta A. Coordinating multiple power-to-gas plants for optimal management of e-fuel seasonal storage. Smart Energy 2024;14:100143. http://dx.doi.org/10.1016/j.segy.2024.100143, URL: https://www.sciencedirect.com/science/article/pii/S2666955224000133.

[21] Abomazid AM, El-Taweel NA, Farag HEZ. Optimal energy management of hydrogen energy facility using integrated battery energy storage and solar photovoltaic systems. IEEE Trans Sustain Energy 2022;13(3):1457-68. http://dx.doi.org/10.1109/TSTE.2022.3161891, URL: https://ieeexplore.ieee.org/document/9740444/.

[22] Darivianakis G, Eichler A, Smith RS, Lygeros J. A data-driven stochastic optimization approach to the seasonal storage energy management. IEEE Control Syst Lett 2017;1(2):394-9. http://dx.doi.org/10.1109/LCSYS.2017.2714426, URL: https://ieeexplore.ieee.org/abstract/document/7945264 Conference Name: IEEE Control Systems Letters.

[23] Thaler B, Posch S, Wimmer A, Pirker G. Hybrid model predictive control of renewable microgrids and seasonal hydrogen storage. Int J Hydrog Energy 2023;48(97):38125-42. http://dx.doi.org/10.1016/j.ijhydene.2023.06.067, URL: https://linkinghub.elsevier.com/retrieve/pii/S0360319923029269.

[24] Cuisinier E, Lemaire P, Penz B, Ruby A, Bourasseau C. New rolling horizon optimization approaches to balance short-term and long-term decisions: An application to energy planning. Energy 2022;245:122773. http://dx.doi.org/10.1016/j.energy.2021.122773, URL: https://www.sciencedirect.com/science/article/pii/S036054422103022X.

[25] Cuisinier E, Lemaire P, Ruby A, Bourasseau C, Penz B. Impact of operational modelling choices on techno-economic modelling of local energy systems. Energy 2023;276:127599. http://dx.doi.org/10.1016/j.energy.2023.127599, URL: https://www.sciencedirect.com/science/article/pii/S0360544223009933.

[26] Vallée M, Bavière R, Seguin V, Vuillerme V, Lamaison N, Descamps M, Aurousseau A. An efficient co-simulation and control approach to tackle complex multi-domain energetic systems: Concepts and applications of the PEGASE platform. In: ECOS 2019 - proceedings of the 32nd international conference on efficiency, cost, optimization, simulation and environmental impact of energy systems. 2019, p. 2853-64, URL: https://www.scopus.com/inward/record.uri?eid=2-s2.0-85079592390&partnerID=40&md5=79cc82ddedfb7bf976fdc399f3583027, Type: Conference paper.

[27] Vielma JP, Nemhauser GL. Modeling disjunctive constraints with a logarithmic number of binary variables and constraints. Math Program 2011;128(1):49-72. http://dx.doi.org/10.1007/s10107-009-0295-4.

[28] Poncelet K, Hoschle H, Delarue E, Virag A, Drhaeseleer W. Selecting representative days for capturing the implications of integrating intermittent renewables in generation expansion planning problems. IEEE Trans Power Syst 2017;32(3):1936-48. http://dx.doi.org/10.1109/TPWRS.2016.2596803, URL: http://ieeexplore.ieee.org/document/7527691/.

[29] Wets RJ-B. Stochastic programs with fixed recourse: The equivalent deterministic program. SIAM Rev 1974;16(3):309-39. http://dx.doi.org/10.1137/1016053.

[30] Nepveu F, Mathonniere S, Lamaison N, Enee G. Presentation, validation and application of the EnergyProcess library. 2023, p. 311-20. http://dx.doi.org/10.3384/ecp204311, URL: https://ecp.ep.liu.se/index.php/modelica/article/view/939.

[31] ElectricityMaps. Electricitymaps. 2023, URL: http://electricitymap tmrow.co.