The state of charge evolution for the battery energy storage system is modeled as:

where S_{ES, t}, S_{ES, t−1} represent the total energy stored in BES at time t and t−1, respectively. P_{c, t}, P_{d, t} represent the charging and discharging power of BES, respectively. η_c, η_d represent the charging and discharging efficiency of BES, respectively. Δt represents the time interval.

The BES model captures the dynamics of energy storage, accounting for conversion losses during both charging and discharging processes [20]. This bidirectional capability is crucial for providing a rapid response to grid fluctuations.

The electrolyzer converts electrical energy into hydrogen through electrolysis, with the following power relationship:

where P_{ele, t}, Q_{ele, t} represent the electrical power consumption and heat generation of the electrolyzer at time t, respectively. ṅ_{H2, t}, η_ele represent the hydrogen production rate and efficiency of the electrolyzer, respectively. H_HHV represents the higher heating value of hydrogen.

The electrolyzer model characterizes the energy conversion from electricity to hydrogen, accounting for the electrochemical process efficiency and associated thermal effects [21, 22]. This process is particularly valuable for utilizing surplus renewable energy that might otherwise be curtailed.

The PEMFC converts hydrogen back to electricity with the following relationship:

where P_{fuel, t}, Q_{fuel, t} represent the electrical power output and heat generation of PEMFC at time t, respectively. ṁ_{H2, t}, η_fuel represent the hydrogen consumption rate and efficiency of PEMFC at time t, respectively.

The PEMFC model captures the reverse conversion process from hydrogen to electricity, completing the energy cycle within the EHESS [23]. This bidirectional conversion capability between electricity and hydrogen provides extended storage duration compared to battery-only systems. The HSS model is similar to the BES model and follows analogous state-of-charge principles, with appropriate adjustments for the physical properties of hydrogen storage.

Demand-side response loads can provide regulatory capacity for renewable energy grids [25]. The specific mathematical model is formulated as:

where, FCL,j,tup, FCL,j,tdown represent the upward and downward regulation capabilities provided by the demand-side load at node j at time t. v_{j, t} represents the state variable indicating whether the demand-side load at node j participates in regulation at time t. P_{j, CL, t} represents the adjustment amount of the demand-side load at node j at time t.

This model accounts for the sequential nature of load control decisions, where the transition between on and off states determines the available flexibility [26]. Demand response resources complement generation-side flexibility by offering rapid load adjustments that can help balance the system.

The flexibility demand is driven by the variability of the net load, which is the load minus renewable output [27]. The upward and downward flexibility demands are defined as the positive ramp rates of the net load [28]:

where FADN,tup, FADN,tdown represent the upward and downward flexibility demands of the renewable energy grid at time t. P_{NL, t} represents the net load output of the renewable energy grid at time t. P_{load, t}, P′RE,t represent the load and forecasted renewable energy output of the grid at time t, respectively. P_{CL, t} represents the demand-side load adjustment of the renewable energy grid at time t.

Based on [29], we decompose flexibility supply as:

where FCL,tup, FCL,tdown represent the upward and downward flexibility regulation capabilities provided by demand-side loads at time t. FES,tup, FES,tdown represent the upward and downward flexibility regulation capabilities provided by EHESS at time t. Ffuel,tup, Fele,tdown represent the upward and downward flexibility regulation capabilities provided by PEMFC and ET at time t, respectively.

A positive margin indicates the system has sufficient flexibility to handle renewable fluctuations. A negative margin implies a flexibility deficit, necessitating load shedding or renewable curtailment [30].

The upper layer planning uses the annual comprehensive cost of the renewable energy grid with EHESS as the optimization objective to optimize the configuration of the electro-hydrogen storage system.

The lower layer planning optimizes the operation of the configured EHESS. While the conceptual goal is to simultaneously minimize operation costs and maximize flexibility margins, we solve this bi-objective problem operationally using the Penalty Function Method.

We introduce a Flexibility Penalty Cost (C_FL) into the objective function. This transforms the maximization of flexibility margins into the minimization of penalty costs associated with flexibility deficits. Thus, the lower-layer problem is formulated as a single-objective minimization:

where C_OP includes grid interaction costs, fuel costs, and maintenance costs (Equation 11), and C_FL penalizes any negative deviation in the flexibility margin (Equation 12). By assigning appropriate penalty coefficients (c_up, c_down), the solver is driven to find an operational schedule that maximizes flexibility (avoiding penalties) while minimizing standard operating expenses.

The AC power flow constraints are modeled as:

where r_ij, x_ij represent the resistance and reactance of branch ij. u(j), v(j) represent the set of nodes with branch j as the terminal node. Pjt, Qjt, Pijt, Qijt represent the active and reactive input power at node j and the active and reactive output power of branch ij. Pj,subt, Qj,subt represent the active and reactive power exchanged between the main grid and the renewable energy grid. Iijt, Uit represent the squared values of branch current and node voltage. Pj,loadt, Pj,CLt, Qj,loadt, Qj,CLt represent the ordinary load and demand-side response load active and reactive power at node j.

These constraints ensure that the power system operates within physical limits while satisfying Kirchhoff's laws.

Based on the distributed robust optimization method [32], the planning model considering flexibility for electro-hydrogen energy storage in renewable energy grids can be written as:

where x, X represent the decision variables and feasible region of Stage 1. g(x), ξ represent the objective function of Stage 1 and uncertainty variables. P represents the probability distribution function of ξ. y represents the decision variables of Stage 2. Ω(x, ξ) represents the feasible region of y. f(y) represents the objective function of Stage 2. 𝔼_P(·) represents the expectation function under distribution function P.

This formulation enables decision-making under uncertainty by considering a range of possible scenarios rather than relying on point forecasts.

Let W(P₁, P₂) represent the Wasserstein distance between distributions P₁ and P₂, defined as:

Where Π represents the joint distribution of P₁ and P₂. S_{ES, t}, S_{ES, t−1} represents the distance between ξ₁ andξ₂.

The ambiguity set P for the EHESS planning model can be written as:

where M(ϖ) represents the set of all distributions on the polyhedron. ϖ={ξ∈RW:Hξ≤h} is the polyhedron.

The Wasserstein distance-based ambiguity set allows us to account for distributional uncertainty in renewable energy outputs, providing robustness against forecasting errors and variability.

For ease of description, the electro-hydrogen energy storage planning model can be rewritten in the following abstract form:

subject to:

where □ represents the feasible region of the first-stage decision variables. A_j(Y), b_j(x) represent the corresponding matrix and vector. j∈J≜{GG, GD}, with GG and GD representing unit and network constraints, respectively.

The feasible region of the second-stage objective function in the distributed robust planning is:

This linearization facilitates the solution process while maintaining the key characteristics of the original model.

Applying duality theory, we transform the maximization problem in the distributed robust model into a minimization problem, converting (Equation 19) into:

where λ₀, s₀, γi0 represent auxiliary variables in duality theory. N represents the number of training samples. t represents the dual norm of t−1.

This transformation produces a computationally tractable formulation of the robust optimization problem.

We decompose matrix A(Y) and vector b(x) as:

For computational efficiency, we use the Bonferroni inequality to decompose the original constraints into $m$ more manageable constraints, approximating the feasible region Ω_CC as Ω_B:

According to the CVaR theorem in [33], the feasible region Ω_B can be written as:

Based on Theorem 2 in [34], (Equation 28) can be transformed into:

The proposed distributed robust optimization model for electro-hydrogen energy storage planning is implemented in MATLAB R2016b (MathWorks, Natick, MA, United States) and solved using the CPLEX solver, providing an efficient computational framework for practical applications.

This study uses a modified IEEE 33-node system based on a renewable energy grid in Liaoning Province, China, as the test case. The network topology is illustrated in Figure 2, which includes wind turbines (WT), photovoltaic systems (PV), and controllable loads (CL) at various nodes. The system integrates both renewable generation sources and flexible demand-side resources, representing a typical modern distribution network with high renewable penetration. To represent the uncertainty of renewable generation, we utilized historical data from the Liaoning Province grid. We applied K-means clustering to generate typical scenarios. To avoid arbitrary scenario selection, we employed the Elbow Method to determine the optimal number of clusters. The analysis of the sum of squared errors (SSE) indicated that K = 4 provides the optimal balance between scenario representativeness and computational efficiency for the robust optimization model.

The economic parameters used in this study are based on market data and industry standards. The demand-side load cost is set at 0.1 yuan/(kW·h). For Battery Energy Storage (BES), the investment cost is 1,130 yuan/kW, and the operation cost is 0.74 yuan/kW. The investment costs for the Electrolyzer (ET) and Proton Exchange Membrane Fuel Cell (PEMFC) are 2,340 yuan/kW and 1,410 yuan/kW, respectively. These parameters reflect the current market conditions and technological maturity of each component in the electro-hydrogen energy storage system.

The renewable generation and load data utilized in this study originate from actual operational records of a regional power grid in Liaoning Province, China, spanning the year 2022 with hourly temporal resolution. The dataset comprises 8,760 observations for each variable, including wind power output from wind farms with a combined capacity of 45 MW, photovoltaic generation from distributed installations totaling 12 MW, and aggregated load demand from mixed residential-commercial consumers. All data were collected through the grid operator's SCADA system (the SCADA system refers to the grid operator's existing supervisory control and data acquisition infrastructure used for data collection—it is not a specific commercial product purchased by the authors, so manufacturer details are not applicable). The wind power data exhibit characteristic patterns typical of northern China, with higher output during winter and nighttime hours, while photovoltaic profiles follow expected solar irradiance patterns. The load profiles display a double-peak pattern with morning and evening peaks, alongside seasonal variations driven by heating and cooling demands.

To account for the stochastic nature of renewable energy resources and load variations, we employed K-means clustering on historical operational data from the Liaoning Province grid. The raw dataset comprising 365 daily profiles was normalized to eliminate scale differences among variables. The optimal number of clusters was determined using the elbow method combined with silhouette coefficient analysis, indicating that four representative scenarios provide an appropriate balance between computational tractability and statistical fidelity. Each cluster centroid represents a typical daily pattern, with probability weights corresponding to the proportion of historical days in each cluster. The four scenarios capture: (i) high-solar, moderate-wind days in spring and autumn (probability 0.28); (ii) low-solar, high-wind winter days (probability 0.31); (iii) moderate renewable output with average loads (probability 0.25); and (iv) peak-load days with below-average renewable generation (probability 0.16). This approach preserves the statistical characteristics of the original dataset while reducing computational burden.

As shown in Figure 3, photovoltaic output follows a bell-shaped curve with peak generation occurring around midday, reaching approximately 1,750 kW in the highest scenario. The intermittent nature of PV generation is evident, with zero output during nighttime hours. Figure 4 illustrates the wind power profiles, which exhibit less predictable patterns with generation capacity ranging from 650 kW to 2,900 kW across different scenarios. Wind power tends to be higher during nighttime and early morning hours in this region. Load profiles in Figure 5 display the characteristic double-peak pattern of residential and commercial consumption, with morning and evening peaks reaching up to 3,200 kW in high-demand scenarios.

To validate the effectiveness and economic efficiency of the proposed electro-hydrogen energy storage planning model for renewable energy grids, we established three comparative schemes:

As the energy transition progresses, increasing penetration of wind and solar power places growing demands on system flexibility. To evaluate this impact, three renewable penetration scenarios are considered: a Baseline level, a Medium Growth level with a 25% increase in wind and solar output, and a High Growth level with a 50% increase. The optimized economic performance of the three schemes under these scenarios is summarized in Table 3.

As shown in Table 3, total system costs increase with higher renewable penetration for all schemes; however, the magnitude of this increase differs substantially, reflecting their varying abilities to accommodate rising flexibility requirements. Scheme 1, which operates without energy storage, exhibits the most pronounced cost escalation. Its total cost rises from 3,636.1 × 104 yuan in the baseline scenario to 8,159.3 × 104 yuan under high growth, indicating limited capability to cope with the increased variability introduced by renewable generation.

Scheme 2 introduces EHESS but lacks coordinated flexibility-oriented planning. As a result, its total cost remains consistently high across all scenarios, increasing from 5,879.1 × 104 yuan at baseline to 9,461.2 × 104 yuan at the high growth level. This suggests that storage investment alone, without strategic adaptation to flexibility demand, is insufficient to effectively mitigate the economic impacts of higher renewable penetration.

In contrast, Scheme 3 demonstrates superior adaptability and scalability. Its total cost increases more moderately, from 3,249.4 × 104 yuan in the baseline scenario to 4,886.5 × 104 yuan under high renewable penetration, significantly lower than the corresponding costs of Schemes 1 and 2. This indicates that the proposed method can proactively respond to increasing flexibility requirements, enabling more efficient system operation and restraining excessive cost growth. These results confirm the economic robustness of Scheme 3 and its suitability for future power systems with high shares of renewable energy.

Figure 10 illustrates the 24-h coordinated operation of EHESS components in Scheme 3. Electrolyzers operate primarily during midday, absorbing surplus photovoltaic generation. Battery storage manages short-term fluctuations, charging during low-demand periods and discharging during evening peaks. Fuel cells provide sustained generation during periods of low renewable output, enabling long-duration energy shifting and maintaining positive flexibility margins.

Figure 11 further reveals that operational costs are dominated by EHESS operation and maintenance, while grid power purchases are significantly reduced. This confirms that active EHESS utilization effectively offsets external energy procurement costs.

Table 4 presents the sensitivity of Scheme 3 to ±20% variations in equipment costs. Reductions in capital cost lead to increased EHESS capacity deployment and an 8.1% decrease in total system cost. Conversely, higher costs result in moderated investment and a controlled increase in operational and penalty costs. Notably, the proportional mix of EHESS components remains stable across scenarios, indicating a structurally robust configuration insensitive to moderate cost fluctuations.

The computational performance of the proposed framework merits discussion for assessing its applicability to larger systems. For the IEEE 33-node system, the two-layer optimization model comprises approximately 2,400 continuous variables, 180 binary variables, and 5,600 constraints when considering four scenarios with 24 h periods. Using CPLEX 12.8 on a desktop computer with an Intel Core i7-8700 processor (Intel Corporation, Santa Clara, CA, United States) processor and 16 GB RAM, the average computation time was 847 s with a 0.1% optimality gap tolerance.

The computational complexity scales primarily with three factors: candidate EHESS nodes, uncertainty scenarios, and temporal resolution. The Wasserstein distance-based uncertainty quantification maintains linear scaling with training samples, offering advantages over moment-based methods. For extension to larger systems, several strategies can be employed: decomposition algorithms such as Benders decomposition can exploit the two-stage structure for parallel subproblem solutions; network reduction techniques can decrease dimensionality while preserving essential characteristics; and warm-starting strategies can accelerate convergence. Preliminary experiments on the IEEE 118-node system indicate that the framework with Benders decomposition achieves solutions within 5% of optimality in approximately 3 h, demonstrating reasonable scalability for practical applications.