A Tri-Level Transaction Method for Microgrid Clusters Considering Uncertainties and Dynamic Hydrogen Prices

Abstract

The advancement of hydrogen technology and rising environmental concerns have shifted research toward renewable energy for green hydrogen production. This study introduces a novel tri-level transaction methodology for microgrid clusters, addressing uncertainties and price fluctuations in hydrogen. We establish a comprehensive microgrid topology with distributed power generation and hydrogen production facilities. A polygonal uncertainty set method quantifies wind and solar energy uncertainties, while an enhanced interval optimization technique refines the model. We integrate a sophisticated demand response model for hydrogen loading, capturing users’ behavior in response to price changes, thereby improving renewable energy utilization and supporting economically viable management practices. Additionally, we propose a tri-level game-theoretic framework for analyzing stakeholder interactions in microgrid clusters, incorporating supply–demand dynamics and a master–slave structure for microgrids and users. A distributed algorithm, “KKT & supply-demand ratio”, solves large-scale optimization problems by integrating Karush–Kuhn–Tucker conditions with a heuristic approach. Our simulations validate the methodology, demonstrating that accounting for uncertainties and dynamic hydrogen prices enhances renewable energy use and economic efficiency, optimizing social welfare for operators and economic benefits for microgrids and users.

1. Introduction

To address resource constraints and environmental challenges while aligning with the United Nations’ UNEP Sustainable Development Goals [1,2], hydrogen energy offers high energy density, environmental sustainability, extensive applications, and abundant resources. It is widely acknowledged as a viable solution for global energy transformation and sustainable development [3,4].

Hydrogen refueling stations are foundational and crucial in advancing the hydrogen energy sector [5]. However, many stations rely on an external supply of high-pressure gaseous hydrogen, raising significant concerns regarding transportation costs and safety [6]. Self-constructed hydrogen refueling stations can streamline hydrogen transportation, achieving both cost-efficiency and zero carbon emissions during operation [7]. Additionally, renewable energy such as wind and solar offers complementary benefits over time scales, enabling the flexible adjustment of hydrogen energy production. Therefore, the reliability of hydrogen supply can be further enhanced through the strategic planning of wind and photovoltaic power generation [8]. Researchers have investigated supply facilities that employ renewable energy for on-site hydrogen production. It is proposed that the hydrogen market integrates renewable energy sources upstream with hydrogen fuel cells downstream [9]. It has been fully demonstrated that multi-energy complementarity offers significant potential for energy consumption. A high-proportion renewable energy multi-energy system is constructed, utilizing large-scale hydrogen production units to balance load demands across multi-energy systems and enhance renewable energy consumption [10]. An optimization scheduling model for private hydrogenation stations is proposed, which facilitates electrolytic water hydrogen production through grid integration and incorporates a dynamic pricing mechanism to ensure profitability [11].

Notably, the existing research predominantly concentrates on optimizing hydrogen energy systems based on energy generation and production, often overlooking the potential impacts on the broader energy landscape [12]. To address uncertainties in energy output, research methods include opportunity constraint planning, scenario analysis, robust optimization, and interval analysis, as shown in Table 1. This paper addresses the solution of large-scale interval uncertainty problems, selecting the interval method based on the analysis of tabular data. Furthermore, to enhance the method’s handling of uncertain interval widths, this paper introduces a convex hull algorithm to establish polygonal uncertainty sets and optimize interval numbers. The improved method significantly enhances the accuracy of traditional interval uncertainty models while reducing the computational complexity.

Table 1. Analysis of methods for dealing with uncertainty problems.

Method	Illustration	Positive	Drawback	Ref.
Constrained planning method	Effectively reflecting the reliability problem of satisfying system constraints under uncertain conditions.	Reduces the complexity of the model; the solution results have the advantages of low cost or high efficiency.	Difficulty in dealing with large-scale scenario problems and violation risk issues in structural constraints.	[13]
Scene analysis method	By predicting and analyzing a series of possible values of uncertain factors in the future environment.	The scene display assists decision-makers in better understanding the problem; the computation scale can be reduced, the efficiency can be improved.	The process is subjective, which affects the accuracy of the results.	[14]
Robust optimization method	The constraints are satisfied within the range of uncertain parameters or in the “worst-case” operating scenario.	Strong reliability, ensuring that constraints are still met even in the harshest scenarios.	The output results are too conservative; high demand for computing resources and complexity of process.	[15]
Interval method	Expands the uncertainty of point variables into interval representations and transforms uncertainty problems into deterministic boundary problems based on interval arithmetic.	The computational complexity is small, suitable for handling large-scale problems; being able to provide a range of uncertain values provides a certain basis for decision making.	The handling of uncertain interval width significantly affects the results; the computational complexity is high when dealing with high-dimensional problems.	[16]

Table 1. Analysis of methods for dealing with uncertainty problems.

Method	Illustration	Positive	Drawback	Ref.
Constrained planning method	Effectively reflecting the reliability problem of satisfying system constraints under uncertain conditions.	Reduces the complexity of the model; the solution results have the advantages of low cost or high efficiency.	Difficulty in dealing with large-scale scenario problems and violation risk issues in structural constraints.	[13]
Scene analysis method	By predicting and analyzing a series of possible values of uncertain factors in the future environment.	The scene display assists decision-makers in better understanding the problem; the computation scale can be reduced, the efficiency can be improved.	The process is subjective, which affects the accuracy of the results.	[14]
Robust optimization method	The constraints are satisfied within the range of uncertain parameters or in the “worst-case” operating scenario.	Strong reliability, ensuring that constraints are still met even in the harshest scenarios.	The output results are too conservative; high demand for computing resources and complexity of process.	[15]
Interval method	Expands the uncertainty of point variables into interval representations and transforms uncertainty problems into deterministic boundary problems based on interval arithmetic.	The computational complexity is small, suitable for handling large-scale problems; being able to provide a range of uncertain values provides a certain basis for decision making.	The handling of uncertain interval width significantly affects the results; the computational complexity is high when dealing with high-dimensional problems.	[16]

Additionally, uncertainties in the output and load of new energy sources lead to imbalances in hydrogen energy supply and demand, which could impact the economic efficiency and reliability of hydrogenation stations [17]. Currently, hydrogen fuel prices are predominantly fixed and fail to adjust dynamically with hydrogen energy supply and demand. This rigidity hampers the optimization of profits for new energy hydrogen production stations and the efficient consumption of new energy [18]. Therefore, implementing a dynamic hydrogen pricing strategy is proposed. The hydrogen price would be adjusted dynamically based on fluctuations in energy output and hydrogen storage. This would enable the hydrogen load to respond in real time, facilitating effective interaction between supply and demand and achieving equilibrium in hydrogen energy supply. The recent literature has explored dynamic hydrogen pricing. It is suggested that when new energy generation exceeds load requirements, electricity prices should be reduced to encourage timely electric vehicle charging, thereby promoting the consumption of new energy [19]. In the power distribution system, a solution for optimal current conditions yields marginal distribution prices for various hydrogen refueling stations. A user demand model based on traffic data determines hydrogen demand, and a 24 h optimal hydrogen pricing scheme is developed, integrating power distribution costs and hydrogen demand [20]. However, this method does not account for the impact of hydrogen storage levels on pricing decisions. By implementing dynamic adjustment strategies, prices can be real-time adjusted based on various factors such as market supply and demand, cost fluctuations, and policy directions. This real-time responsiveness and flexibility allow for a more accurate reflection of market conditions, thereby enhancing the effectiveness of the pricing mechanism.

The introduction of microgrid clusters (MGs) enhances the flexibility and economic efficiency of regional energy systems by optimizing energy scheduling to achieve overall regional energy savings and efficiency gains [21]. However, increased system degrees of freedom and complex stakeholder interactions make optimizing the entire energy network challenging. In [22,23], the authors propose a principal–agent game model involving a microgrid cluster operator (MGCO) and microgrid (MG), providing a solution for the equilibrium by using particle swarm optimization. However, this method is prone to local optima, compromising solution reliability. Ref. [24] establish a bi-level optimization model of a principal–agent game between energy service providers and power market operators, solved using Karush–Kuhn–Tucker (KKT) conditions. However, this method suffers from slow convergence rates. Ref. [25] introduce a tri-level market with a multi-leader–multi-follower principal–agent game model, obtaining equilibrium trading strategies through iterative methods. Overall, existing research has primarily focused on the operational mechanisms and business models without delving into the interests and interactions among multiple stakeholders in microgrid groups that incorporate new energy into the Internet. Ref. [26] suggest the use of model-free deep reinforcement learning to address energy management issues in the MG sector of the power industry, achieving significant success in optimizing decision making. However, this approach requires a substantial amount of training data, which considerably increase the complexity of the problem. Furthermore, the algorithm is trained on data from a single location, which limit its scalability. Therefore, this study constructs a tri-level trading system comprising MGCOs, MGs, and DUs and employs distributed algorithms to solve equilibrium strategies among the three stakeholder levels in microgrid clusters incorporating renewable energy. The constructed model encompasses the interactions among multiple stakeholders, complicating the solution process. Currently, no efficient method exists to tackle these optimization challenges. Traditional bi-level optimization typically employs two primary approaches: one utilizes the KKT conditions to convert the bi-level problem into a single-level format, solvable through mixed-integer linear programming; effective linear transformation techniques are essential here. The other approach involves designing distributed iterative algorithms based on analytical or heuristic methods. Although the KKT method offers rapid computation, it requires direct information flow, raising concerns about privacy. In contrast, distributed iterative algorithms can become complex in large-scale problems, while heuristic algorithms often fall short in computational efficiency. To effectively address these challenges, we propose a “KKT & Supply-Demand Ratio” distributed iterative strategy aimed at solving the proposed three-level optimization model.

The main contributions are as follows: (1) Considering the uncertainty from renewable energy generation, an improved polygonal uncertainty set is proposed to represent intervals. (2) A dynamic hydrogen pricing strategy is proposed, and based on optimized hydrogen prices, a hydrogen demand response mechanism is established to enhance stability and economic efficiency. (3) A tri-level game structure is developed, enabling energy cooperation and mutual benefit among the three stakeholders through game theory.

2. Optimization Model of the System

2.1. Topology of Microgrid Cluster

Figure 1 illustrates the structure of the system. The renewable energy hydrogenation station primarily generates wind and solar energy. It employs electrolysis technology to produce hydrogen, stores it locally in hydrogen storage tanks (HTs), or supplies it to hydrogen loads, achieving the station’s economic operation. To meet energy demands, the system is equipped with wind turbines (WTs) and photovoltaic panels (PVs) as energy production devices. Electrolyzers (ELs) and hydrogen fuel cells (FCs) serve as devices for electricity-to-hydrogen energy conversion; HTs act as energy storage devices, facilitating the absorption of wind and solar energy while providing a reserve for hydrogen demand.

The hydrogenation station’s electrical energy primarily comes from WTs and PVs, with a portion supplied by the external grid. The station operates under two distinct modes: (1) WTs-PVs output power is insufficient. The station purchases electricity from the external grid to ensure daily operation; (2) WTs-PVs output exceeds demand. The station converts a portion of the excess hydrogen into electricity using FCs and sells it to the external grid, thereby generating economic returns.

2.2. Interval Uncertainty Model

Given the topology mentioned above, it is essential to account for the potential negative impacts of WTs-PVs power uncertainty on the system. The study employs an interval number model to depict the uncertainty and subsequently refines the interval width based on this model. Then, an interval linear optimization approach is utilized to optimize the system.

2.2.1. Traditional Uncertainty Interval Model

Traditional interval number models use upper and lower bounds to characterize uncertain variables, lacking explicit requirements for their statistical properties and membership functions. Typically, these models construct uncertainty sets using historical data to derive the interval bounds. The intervals for wind and solar power variables and the total output of renewable energy are represented by Equations (1) and (2), respectively, as follows:

$$\{\begin{array}{c}{P}_{\mathit{WT}}^{\pm }=\left[\begin{array}{ccc}{P}_{\mathit{WT}}^{-}& ,& P_{\mathit{WT}}^{+}\end{array}\right]\\ P_{\mathit{PV}}^{\pm }=\left[\begin{array}{ccc}{P}_{\mathit{PV}}^{-}& ,& P_{\mathit{PV}}^{+}\end{array}\right]\end{array}$$

(1)

$$P_{\mathit{REG}}^{\pm }=\left[\begin{array}{ccc}{P}_{\mathit{WD}}^{-}+P_{\mathit{PV}}^{-}& ,& P_{\mathit{WD}}^{+}+P_{\mathit{PV}}^{+}\end{array}\right]$$

(2)

where P_WT^± and P_PV^± represent the positive and negative intervals of WT and PV power, respectively; P_REG^± is the output power for new energy. In practical wind and solar power plants, the probability of encountering extreme scenarios with either maximum or minimum output is very low or may not occur. Traditional uncertain mathematical models often use excessively broad intervals for renewable energy output, which can exacerbate errors in the optimization process.

2.2.2. Improved Interval Model

This study employs convex hull algorithms to construct polygonal uncertainty sets for interval optimization, enhancing the accuracy of traditional interval uncertainty models and reducing the computational complexity. Polygonal uncertainty sets offer advantages in optimizing two-dimensional variables due to their simple formulation and high accuracy. The convex hull algorithm identifies and connects all boundary points to construct the boundary of the polygonal uncertainty set, thereby narrowing the range of joint wind and solar power output intervals. Each boundary of the polygon can be expressed as a linear function, as shown in Equation (3):

$$MA{X}_{1\le i\le n}{a}_{i}x+b_{i}y+c_i=0$$

(3)

where n is the number of sides of the convex bag; i is the convex boundary; a_i, b_i, and c_i are coefficients of a linear equation calculated based on the coordinates of the two endpoints of each boundary. The uncertainty set modeling is carried out using box boundary and polygon boundary, as shown in Figure 2.

As shown in Figure 2, the upper limit of the interval uncertainty set is often on the boundary of the polygon uncertainty set. The interval width of the polygon uncertainty set is significantly smaller than that of the box uncertainty set. This proves that research methods can significantly enhance the accuracy and precision of system optimization to achieve better economic results in interval optimization.

3. Dynamic Hydrogen Pricing and Demand Response

In the electro-hydrogen system, hydrogen prices are prone to fluctuations due to factors such as renewable energy output and supply–demand dynamics. Therefore, this section proposes a dynamic hydrogen pricing scheme that considers the power output and storage levels of the station. When the P_REG(t) is high, and the E_HT(t) is sufficient, the C_H(t) is adjusted downward to attract and increase hydrogen demand, thereby promoting the absorption of renewable energy. Conversely, when P_REG(t) is low and E_HT(t) is insufficient, C_H(t) is adjusted upward to reduce demand and enhance the economic operation. The model is given by Equations (4) and (5).

$$C_H(t)=C_{H\text{-}\mathit{MAX}}-C_{H\text{-}\mathit{TRANS}}$$

(4)

$$C_{H\text{-}\mathit{TRANS}}={\lambda }_1{P}_{\mathit{REG}}(t)/(P_{\mathit{WT}\text{-}\mathit{MAX}}+P_{\mathit{PV}\text{-}\mathit{MAX}})+{\lambda }_2[E_{\mathit{HT}}(t)-E_{\mathit{HT}}(t-1)]/E_{\mathit{HT}\text{-}\mathit{MAX}}$$

(5)

where at time t, C_H(t) is the hydrogen price; C_H-TRANS(t) is the hydrogen adjustment price; P_REG(t) is the WTs-PVs power output; E_HT(t) is the HTs’ capacity. C_H-MAX and C_H-MIN are the maximum and minimum dynamic hydrogen prices. P_WT-MAX, P_PV-MAX, and E_HT-MAX are the maximum power output of WT, PV, and HT maximum power output. λ₁ and λ₂ are the ratio coefficients of new energy output and hydrogen storage, respectively, representing the weight of promoting the consumption of new energy and the weight of improving the utilization rate of hydrogen storage tanks. Specifically, λ₁ is set to 18 and λ₂ to 0.15.

Based on the aforementioned dynamic hydrogen pricing scheme, hydrogen demand adjusts according to fluctuations in hydrogen prices. Additionally, considering the temporal scale, this approach assists hydrogen stations in adapting to the intermittency and variability of renewable energy output. A hydrogen energy demand response model based on different periods is constructed, as illustrated in Equation (6):

$$\{\begin{array}{l}{\epsilon }_{ii}=\Delta E_i{M}_i/E_i\Delta M_i\\ {\epsilon }_{ij}=\Delta E_i{M}_j/E_i\Delta M_j\end{array}$$

(6)

where ε_ii and ε_ij represent the own-elasticity and the cross-elasticity coefficient. ε_ii represents the demand response of users in time period i to the hydrogen price in time period i; ε_ij represents the demand response of users in time period i to hydrogen prices in other time periods j. ΔE denotes the change in hydrogen load, and E represents the forecasted hydrogen load. ΔM indicates the change in hydrogen price, while M refers to the original hydrogen price.

One day’s scheduling period is divided into n intervals. The hydrogenation demand based on dynamic hydrogen prices across these n intervals collectively forms the elasticity matrix, as shown in Equation (7).

$$\left[\begin{array}{c}{\displaystyle \frac{\Delta E_1}{E_1}}\\ {\displaystyle \frac{\Delta E_2}{E_2}}\\ ⋮\\ {\displaystyle \frac{\Delta E_n}{E_n}}\end{array}\right]=\left[\begin{array}{cccc}{\epsilon }_{11}& {\epsilon }_{12}& \cdots & {\epsilon }_{1n}\\ {\epsilon }_{21}& {\epsilon }_{22}& \cdots & {\epsilon }_{2n}\\ ⋮& ⋮& & ⋮\\ {\epsilon }_{n1}& {\epsilon }_{n2}& \cdots & {\epsilon }_{nn}\end{array}\right]\left[\begin{array}{c}{\displaystyle \frac{\Delta M_1}{M_1}}\\ {\displaystyle \frac{\Delta M_2}{M_2}}\\ ⋮\\ {\displaystyle \frac{\Delta M_n}{M_n}}\end{array}\right]$$

(7)

The hydrogen load at a renewable hydrogenation station is reducible, with the total daily hydrogen load varying in response to price changes. The hydrogen load after responding to hydrogen demand, denoted as E_HRES-n, is given by Equation (8):

$$E_{\mathit{HRES}\text{-}n}=E_n+\Delta E_n$$

(8)

4. Tri-Level Transaction Model

Figure 3 illustrates the framework of the system trading model. The MGCO, positioned at the upper layer of the system, serves as the central hub connecting the internal and external microgrids and distribution networks. It coordinates the supply and demand of microgrid energy, acting as the information collection platform, trading center, and energy flow allocator within the system. MGs, located at the middle layer of the system, connect microgrid operators with end users, facilitating power and information exchanges while independently conducting financial transactions with the MGCO. EUs, situated at the lower layer of the system, connect with MGs to purchase electrical energy.

As shown in Figure 3, the MGCO engages in transactions with the distribution network to establish the electricity trading prices between MGs and each microgrid. MGs set the electricity purchasing prices for EUs based on the trading prices determined by the MGCO and the load demands of EUs, and subsequently provide this information to the MGCO. EUs exhibit autonomous capability, enabling them to adjust their actual energy consumption in response to the trading prices set by MGs, thereby balancing expenditure with energy usage experience. They can also perform demand response and relay feedback to MGs. MGs and EUs conform to a leader–follower game model, with MGs acting as the leaders and EUs as the followers.

4.1. Multi-Agent Transaction Model

4.1.1. MGCO Model

The MGCO participates in microgrid cluster optimization to maximize system social welfare, with a benefit of 0. The constructed benefit objective function is shown in Equation (9):

$$\begin{array}{l}{F}_{UP\text{-}MGCO}={\displaystyle \sum _{t=1}^T\left[k_{SU}^{MGs}(t)E_{SU}^{MGs}(t)-k_{BU}^{MGs}(t)E_{BU}^{MGs}(t)\right]}\\ +{\displaystyle \sum _{t=1}^T\left[{\displaystyle \sum _{n=1}^N\left(k_{BMGs}^{MG}(t)E_{BMGs}^{MG}(t)-k_{SMGs}^{MG}(t)E_{SMGs}^{MG}(t)\right)}\right]}=0\end{array}$$

(9)

where F_UP-MGCO is the benefit, and n is the number of microgrids. k^MGs_SU(t) and k^MGs_BU(t), and E^MGs_SU(t), E^MGs_BU(t) denote the sales, purchase price, and quantities of electricity from the MGCO to the utility grid at time t. k^MG_BMGs(t), k^MG_SMGs(t), and E^MG_BMGs(t), E^MG_SMGs(t) represent the purchase, sale price, and quantities of electricity from the MGCO to MGs at time t.

The electricity transaction discrepancy ΔE(t) between the MGCO and the utility grid is defined as shown in Equation (10):

$$\Delta E\left(t\right)={\displaystyle \sum _{t=1}^T[{\displaystyle \sum _{n=1}^N{E}_{SU}^{MGs}\left(t\right)-E_{BU}^{MGs}\left(t\right)}}]$$

(10)

When ΔE(t) ≥ 0, the MGCO is selling electricity; when ΔE(t) < 0, the MGCO is purchasing electricity.

Based on the price and energy supply–demand relationship, the supply–demand ratio of the MGs can be expressed as shown in Equation (11):

$$K(t)={\displaystyle \sum _{n=1}^N{E}_{SMGs}^{MG}(t)}/{\displaystyle \sum _{n=1}^N{E}_{BMGs}^{MG}(t)}$$

(11)

where K(t) represents the ratio of energy supply to demand in the MGs during time t.

Equations (12) and (13) illustrate the relationship between the supply–demand ratio and price.

$$k_{SMGs}^{MG}(t)=\{\begin{array}{l}{\displaystyle \frac{k_{SU}^{MGs}(t)\left[k_{SU}^{MGs}(t)+k_{BU}^{MGs}(t)\right]}{k_{SU}^{MGs}(t)\left[1+K\left(t\right)\right]+k_{BU}^{MGs}(t)\left[1-K\left(t\right)\right]}},0⩽K\left(t\right)⩽1\hfill \\ k_{BMGs}^{MG}(t)G\left(t\right)+[1-K\left(t\right)]k_{SU}^{MGs}(t),1<K\left(t\right)\hfill \end{array}$$

(12)

$$k_{BMGs}^{MG}(t)=\{\begin{array}{l}{k}_{SMGs}^{MG}(t)K\left(t\right)+k_{BU}^{MGs}(t)[1-K\left(t\right)],0⩽K\left(t\right)⩽1\hfill \\ {\displaystyle \frac{k_{SU}^{MGs}(t)[k_{SU}^{MGs}(t)+k_{BU}^{MGs}(t)]}{k_{SU}^{MGs}(t)[1+G\left(t\right)]+k_{BU}^{MGs}(t)[1-G\left(t\right)]}},1<K\left(t\right)\hfill \end{array}$$

(13)

The MGCO only acts as either a sole supplier or seller when K(t) or G(t) = 0. Currently, the internal electricity trading price is set equal to the price established by the distribution network.

4.1.2. MG Model

(1)

Objective function

The MG derives its revenue primarily from the energy payments made by users and transaction fees with the MGCO, minus the costs associated with energy storage system maintenance, as represented by Equation (14):

$$F_{MID-MG}=max\left[\begin{array}{l}{\displaystyle \sum _{t=1}^T\left({\displaystyle \sum _{i=1}^I{k}_{BMG}^{LOADn}(t)}{E}_{BMG-i}^{LOADn}\right)}+{\displaystyle \sum _{t=1}^T\left(k_{SMGs}^{MG}(t)E_{SMGs}^{MG}(t)-k_{BMGs}^{MG}(t)E_{BMGs}^{MG}(t)\right)}\\ -{\displaystyle \sum _{t=1}^T{C}_{ES}}\left(P_{EL}(t)+P_{FUEL}(t)\right)\end{array}\right]$$

(14)

where F_MID-MG represents the benefit of the nth microgrid; i denotes the ith user within the MG, where i = 1,2, …, I; K^LOADn_BMG(t) and E^LOADn_BMG-i are the electricity purchase price and the load of the ith EUs at time t. C_ES is the operational cost of the energy storage system; P_EL(t) and P_FUEL(t) are the power of the ELs and FCs at time t.

(2)

Constraints

1. Power balance constraints: the system power balance during operation, as shown in Equation (15):

$$P_{PV}(t)+P_{WD}(t)+P_{BMGs}^{MG}(t)+P_{FUEL}(t)=P^{LOAD}(t)+P_{SMGs}^{MG}(t)+P_{EL}(t)$$

(15)

where P^LOAD(t) represents the load power at time t;

2. Energy storage System constraints: it is necessary to model each component of the energy storage system and impose capacity and power constraints, as shown in Equations (16) to (21):

The production rate of the Els is dependent on the input power and operational efficiency, as shown in Equation (16):

$$E_{EL}(t)={\eta }_{EL}{P}_{EL}(t)\Delta t/L_1$$

(16)

$$0\le P_{EL}(t)\le P_{EL,max}$$

(17)

where E_EL(t) and P_EL(t) are the hydrogen capacity and power of the ELs at time t; η_EL is the efficiency; L₁ is the conversion coefficient; P_EL,max is the upper limit of the ELs power output.

In renewable hydrogenation stations, produced hydrogen is stored directly in HTs and simultaneously dispensed to hydrogen loads, as shown in Equations (18) and (19):

$$E_{HT}(t)=E_{HT}(t-1)+E_{EL}(t-1)-E_{HLOAD}(t-1)-E_{FUEL}(t-1)$$

(18)

$$E_{HT,min}\le E_{HT}(t)\le E_{HT,max}$$

(19)

where E_HT,min and E_HT,max represent the minimum and maximum capacity limits of the HTs, respectively.

The FCs utilize proton exchange membrane technology to facilitate the coupling between the hydrogen energy and electrical energy, as shown in Equations (20) and (21):

$$P_{FC}(t)={\eta }_{FC}{E}_{FC}(t)L_2/\Delta t$$

(20)

$$0\le P_{FC}(t)\le E_{FC,max}$$

(21)

where L₂ is the hydrogen-to-electricity conversion coefficient of the FCs.

3. Operational constraints: constraints on the energy storage system’s charging and discharging behaviors, and the MG’s electricity purchasing and selling activities, as shown in Equations (22) and (23):

$$P_{HESS}(t)=a{P}_{EL}(t)+(1-a)P_{FUEL}(t)$$

(22)

$$P_{MG}(t)=b{P}_{MGB}(t)+(1-b)P_{MGS}(t)$$

(23)

where P_HESS(t) is the operational power of the energy storage system at time t. a represents the operational status indicator for the energy storage system, where a = 1, system charging. And a = 0 indicates system discharging. P_MG(t), P_MGB(t), and P_MGS(t) are the operational, purchasing, and selling power of the MGs. b is the operational status indicator; where b = 1, the system is purchasing electricity, and where b = 0, the system is selling electricity.

4.1.3. EUs Model

(1)

Objective function

The objective function of the EUs is defined as the difference between the user’s utility and the cost of electricity procurement. The quadratic function is used as depicted in Equation (24).

$$F_{LOW-EUsn}=max{\displaystyle \sum _{t=1}^T\left[a{E}_{BMGs}^{MG}(t)-\frac{b}{2}{\left(E_{BMGs}^{MG}(t)\right)}^2-k_{SU}^{MGs}(t)E_{BMGs}^{MG}(t)\right]}$$

(24)

where F_LOW-EUsn represents the benefit of user i in the nth microgrid; a and b denote the preference constants for energy consumption.

(2)

Constraints

The constraints of the user model are presented in Equations (25) and (26).

$$\{\begin{array}{l}{E}_{BMGs}^{MG}(t)=\overline{E_{BMGs}^{MG}(t)}+T_{BMGs}^{MG}(t)\\ {\displaystyle \sum _{t=1}^T{T}_{BMGs}^{MG}(t)\Delta t=N_{BMGs}^{MG}}\end{array}$$

(25)

$$0\le T_{BMGs}^{MG}(t)\le T_{BMGs,max}^{MG}$$

(26)

where $\overline{E_{BMGs}^{MG}(t)}$ and $T_{BMGs}^{MG}(t)$ denote the base electrical load and the transferable electrical load of user i in the nth microgrid at time t. $T_{BMGs,max}^{MG}$ is the upper limit of the user’s transferable load; and $N_{BMGs}^{MG}$ indicates the total amount of transferable load over the Δt.

5. Establish and Solve Model

This section develops and solves the optimization model. Firstly, it addresses modeling for interval uncertainty. Subsequently, a tri-level optimization model for multi-stakeholder game theory is constructed based on this foundation. Finally, a distributed method is employed for solving the model, as illustrated in Figure 4.

To address the uncertainty in wind and solar power outputs within defined intervals, we developed a mathematical model. This model integrates uncertainty sets based on spectral box technology to create a polygonal uncertainty framework. Building upon this foundation, we incorporated dynamic hydrogen pricing adjustments and a demand response mechanism, leading to the design of a corresponding system operation strategy. Furthermore, we constructed a tri-level optimization game model involving multiple stakeholders. At the Middle–Lower Level, the MG serves as the leader, while EUs act as followers, establishing a principal–agent game structure. To solve this game problem, we employed the KKT conditions alongside convex optimization techniques, reformulating it as a bi-level mixed-integer quadratic programming model for efficient resolution. To further analyze and optimize system performance, we introduced the concept of supply–demand ratios, subsequently constructing an Upper–Middle Level model. Through iterative calculations, this model produces equilibrium strategy solutions under various conditions, providing a scientific basis for decision making.

5.1. Solve Uncertainty Model

In the interval uncertainty modeling, the decision variables are P_EL(t) and P_FC(t), while the interval variable represents the range of WTs-PVs power output. The interval model is transformed into two sub-models, namely the optimistic and pessimistic models, which are solved using distinct procedures. The optimistic model corresponds to the scenario of maximum WTs-PVs power output, yielding the optimistic optimal value Z⁻, as shown in Equation (27). Conversely, the pessimistic model represents the scenario of minimum output, resulting in the pessimistic optimal value Z⁺, as depicted in Equation (28).

$$\{\begin{array}{l}min{Z}^{-}=A^{-}X+B^{-}\\ \begin{array}{cc}s.t.& C^{-}X\le w^{+}\end{array}\\ \begin{array}{cc}& D^{-}X\le k^{+}\end{array}\\ \begin{array}{cc}& D^{+}X\le k^{-}\end{array}\\ \begin{array}{cc}& L^{-}X\le U\end{array}\end{array}$$

(27)

$$\{\begin{array}{l}min{Z}^{+}=A^{+}X+B^{+}\\ \begin{array}{cc}s.t.& C^{+}X\le w^{-}\end{array}\\ \begin{array}{cc}& D^{-}X\le k^{+}\end{array}\\ \begin{array}{cc}& or{D}^{+}X\le k^{-}\end{array}\\ \begin{array}{cc}& L\le X\le U\end{array}\end{array}$$

(28)

where A∈R^1×Ng, X∈R^Ng×1, B∈R^Ng×1, C∈R^M×Ng, and D∈R^T×Ng; R is the real number matrix; N_g is the number of decision variables; M is the number of interval inequality; T is the number of interval equations; U and L are the set of upper and lower boundaries. Optimistic solution Z⁻ and pessimistic solution Z⁺ constitute the optimal value interval [Z⁻, Z⁺].

5.2. Game Model

The game is specifically characterized as follows: (1) The trading prices set by the MGCO influence the energy trading strategies of the MG, while the real-time prices offered by the MG also affect the load demand of the EUs. (2) The energy purchasing strategies of the EUs, following demand response, impact the MG’s pricing strategy and its energy procurement from the MGCO, thereby indirectly influencing the MGCO’s pricing strategy.

In the iterative interaction process, stakeholders continuously adjust their strategies, ultimately achieving a state of equilibrium. Specifically, the trading strategies between the MG and EUs adhere to a leader–follower game model. By defining the participants, decisions, and payoff functions of G, the optimization problem can be expressed as Equation (29).

$$G=\{S;\{{\gamma }_1,{\gamma }_2\};\{F_1,F_2\left\}\right\}$$

(29)

(1) Participants: S denotes the participants, which include the MG and EUs.

(2) Decision Set: The strategy of the MG is derived from the MGCO’s transaction plan and the real-time electricity purchasing prices set for users. Among those, γ₁ = $[P_B^{MGs}(t),P_S^{MGs}(t),P_{BMG}^{LOADn}(t)]$. The strategy of the EUs, which conforms to demand response, is denoted as γ₂ = $[T_{BMGs}^{MG}(t)]$.

(3) Utility Functions: {F₁, F₂} represents the set of utility functions, which denotes the objective functions of the stakeholders.

5.3. KKT and Convex Optimization

To address the bi-level optimization game model involving the MG and EUs, this paper transforms the lower-level EU model into constraint conditions for the MG optimization model using KKT conditions. Consequently, the bi-level optimization problem is reduced to a single-level optimization problem, as expressed in Equations (30) and (31):

$${\displaystyle \frac{\partial L_{BMGs}^{MG}(t)}{\partial T_{BMGs}^{MG}(t)}}=a-b[\overline{E_{BMGs}^{MG}(t)}+T_{BMGs}^{MG}(t)]-k_{BMG}^{LOADn}(t)+{\alpha }_1-{\alpha }_2+\mu =0$$

(30)

$$\{\begin{array}{l}0\le T_{BMGs}^{MG}(t)\perp {\alpha }_1\ge 0\\ 0\le [T_{BMGs,max}^{MG}-T_{BMGs}^{MG}(t)]\perp {\alpha }_2\ge 0\end{array}$$

(31)

where $L_{BMGs}^{MG}$ is the Lagrange function constructed by the EUs, with Formula (28) as a partial derivative and set to 0;

α₁, α_2, and μ are dual variables introduced in the KKT process. Among those, 0 ≤ A⊥B ≥ 0, means a ≥ 0, b ≥ 0 and a·b = 0.

There are nonlinear items in the constraint condition, so it needs to be convex optimized. In the converted single-layer model, the constraint Equation (31) is a nonlinear constraint. And $k_{SU}^{MGs}(t)E_{BMGs}^{MG}(t)$ is a nonlinear constraint in Equation (24). Thus, we adopt the large M method, by attracting several 0–1 variables, and the original nonlinear constraint Equation (31) is converted into a mixed-integer linear constraint, as shown in Equation (32):

$$\{\begin{array}{l}0\le T_{BMGs}^{MG}(t)\le {\beta }_{1}M\\ 0\le {\alpha }_1\le (1-{\beta }_1)M\\ 0\le [T_{BMGs,max}^{MG}-T_{BMGs}^{MG}(t)]\le {\beta }_{2}M\\ 0\le {\alpha }_2\le (1-{\beta }_2)M\end{array}$$

(32)

where β₁ and β₂ are the binary variables; M is a constant large number.

5.4. Distributed Method Solution Process

This article proposes a distributed algorithm to solve the master–slave game model of “MGCO-MG-EUs”. The specific solution process is shown in Figure 5.

(1) Solving the Interval Optimization Model

The initial step involves the initialization of system parameters, followed by the input of uncertainty ranges within the WTs-PVs framework, aimed at refining the interval range. Subsequently, the interval optimization objective function, denoted as min[f], is established. The optimistic and pessimistic models are then solved to yield the corresponding optimal solutions, min f⁻ and min f⁺, along with the dynamic hydrogen prices. These dynamic prices, informed by demand response mechanisms, are generated for both scenarios and subsequently integrated into the tri-level optimization model.

(2) Solving the Tri-Level Stakeholders Trading Model

The trading prices of the MGCO are initially defined as the foundational input conditions for the optimization process. Following this, the KKT conditions and convex optimization techniques are employed to delineate the principal–agent game strategies between the MG and its users. The optimized energy purchase plans from each microgrid are subsequently relayed back to the MGCO. In its capacity as a public welfare entity, the MGCO engages with microgrids to enhance social welfare, iteratively adjusting the trading prices until convergence is reached and an equilibrium solution is established. Throughout each iteration, participants exchange only essential information, such as energy trading prices and purchase strategies, thereby preventing information leakage and ensuring the confidentiality of all parties involved.

6. Case Analysis

6.1. Example Setting

The case study constructs a model of a WTs-PVs AC-DC grid-connected microgrid as the smallest unit for simulation experiments. Wind and solar output forecasting data are based on historical wind and solar data from typical summer and winter days (June and December) at a specific location. Hydrogen load data for the refueling station are shown in Figure 6. Time-of-use electricity prices are presented in

Table 2. Time-of-use electricity price.

Period Type	Period	Electricity Purchase Price ($/kWh)	Electricity Sales Price ($/kWh)
Low valley	23:00–07:00	0.055	0.042
Daily	11:00–14:00, 18:00–23:00	0.105	0.042
Peak	07:00–11:00, 14:00–18:00	0.162	0.042

. Elasticity coefficients are adjusted for different periods, with two distinct coefficients set to increase the hydrogen load during lower price periods. The specific settings are detailed in

Table 3. Elasticity coefficient setting.

Classification	Period	Value
Sub elasticity coefficient Ε1	00:00–24:00	−1.5
Mutual elasticity coefficient Ε2	06:00–10:00, 18:00–24:00	0.03
Mutual elasticity coefficient Ε3	11:00–17:00	0.05
Mutual elasticity coefficient Ε4	01:00–05:00	0

. Multiple scenarios are constructed to validate the model, as shown in

Table 4. Scenario condition setting.

Scenarios	Multilateral Uncertainty Set	Dynamic Hydrogen Price	Demand Response
1	×	×	×
2	√	×	×
3	√	√	×
4	√	√	√

. The case study is modeled using the MATLAB 2021b platform and solved with the Gurobi solver, with simulation performed on hardware equipped with an Intel Core i7-12700H CPU @ 2.30 GHz.

6.2. Results Analysis

6.2.1. Uncertainty Results Analysis

This section aims to validate the effectiveness of the proposed method for handling uncertainties by analyzing the operation of the hydrogenation station under different scenarios, as shown in

Table 5. Operational performance of hydrogenation stations.

Scenarios	Typical Days	The Price of Hydrogen ($/kg)	Total Hydrogen Load (kg)	Revenue ($)
1	Summer	5.601	289.6	[885.54, 1661.88]
	Winter	5.601	289.6	[1050.64, 1817.41]
2	Summer	5.601	289.6	[1070.51, 1687.47]
	Winter	5.601	289.6	[1053.94, 1820]
3	Summer	[4.901, 6.301]	289.6	[1133.91, 1534.37]
	Winter	[4.901, 6.301]	289.6	[1104.84, 1719.43]
4	Summer	[4.901, 6.301]	[273.1, 301.4]	[1206.34, 1741.3]
	Winter	[4.901, 6.301]	[273.1, 301.4]	[1103.27, 1666.54]

. Furthermore, Figure 7 compares the output ranges for Scenario 1 and Scenario 2 during typical days.

Analysis of Table 5 and Figure 7 reveals that in summer, the increase in the minimum wind and solar output results in a reduction in the electricity purchase costs for the hydrogen refueling station and an improvement in the minimum revenue. In winter, the reduced sunlight leads to a decrease in the peak renewable energy output, and the difference between the optimized polygonal uncertainty set interval and the original box-based interval is minimal. Therefore, it can be concluded that when wind and solar output is excessive in summer and insufficient in winter, optimizing the wind and solar data using the polygonal uncertainty set results in improved revenue intervals for the hydrogen refueling station.

Further analysis of Scenarios 1 and 2 in Figure 7 indicates that the original data align with the energy output boundary defined by the box-type uncertainty set. In contrast, the optimized data conform to the same boundary defined by the polygonal uncertainty set. As demonstrated in Figure 7, applying polygonal uncertainty in power output optimization effectively reduces some extreme scenarios, leading to a decrease in the upper boundary and an increase in the lower boundary. Consequently, the interval derived from the polygonal uncertainty set offers a more precise representation than that obtained from the box-type uncertainty set.

6.2.2. Dynamic Hydrogen Price and Demand Response Results Analysis

Scenario 3 is constructed to assess the validity of the dynamic hydrogen pricing. As illustrated in Figure 8, the dynamic variations in hydrogen prices are demonstrated.

Analysis combining Figure 8 and Table 5 reveals that during the 11:00–14:00 period, when the total power output is high, hydrogen prices significantly decrease. Conversely, in the evening, when hydrogen demand is high, and hydrogen storage fluctuates more frequently, hydrogen prices notably increase. Considering the dynamic hydrogen prices influenced by the renewable energy output and hydrogen storage levels, these prices fluctuate due to the variations in the power output and hydrogen storage. Furthermore, in the optimization process, the hydrogen production volume primarily depends on hydrogen demand, resulting in no significant changes in the hydrogen production volume at these stations. Fluctuations in hydrogen prices directly affect the revenue range of refueling stations. Specifically, the minimum revenue increases while the maximum revenue decreases, which adversely impacts the profitability of these stations. Therefore, adopting a demand response model based on dynamic hydrogen prices, which adjusts hydrogen demand according to these prices, is a crucial strategy for enhancing the economic efficiency of the system.

Furthermore, the hydrogen load curve is illustrated in Figure 9. An analysis of Figure 8 and Figure 9 reveals that on a typical winter day, when power output is at its lowest, the minimum revenue in Scenario 4 is slightly lower than in Scenario 3 but still higher than in Scenario 2 due to the reduction in hydrogen load. This demonstrates that the introduction of dynamic hydrogen pricing strategies and demand response mechanisms can enhance the revenue of hydrogen refueling stations under the worst-case scenarios. Comparative analysis of typical day results shows that during seasons with higher overall power output, dynamic hydrogen pricing strategies and demand response mechanisms can significantly increase the revenue of hydrogen refueling stations. The results thoroughly validate the proposed dynamic hydrogen pricing strategy and hydrogen load demand response mechanism, demonstrating effectiveness in achieving interaction and enhancing the revenue of hydrogenation stations during periods of both surplus and deficit in energy power output.

6.2.3. Transaction Result Analysis

• MG Operation Result Analysis

Figure 10 presents the optimized electrical output results for each microgrid within the cluster. MG2 typically functions as the electricity demand entity, while MG1 and MG3 primarily serve as the electricity supply sources.

It is noteworthy that there is a significant disparity in the energy supply and demand at the peak electrical output of MG1 and MG3. Consequently, the MGCO integrates resources such that surplus electricity within the cluster is preferentially sold to other microgrids experiencing a shortage, facilitating energy exchange between microgrids.

• EUs Optimization Result analysis

Table 6. User economic optimization results.

User Type	Optimize the Goal ($)			Cost Function ($)
	Before Optimization	After Optimization	Optimization Rate	Before Optimization	After Optimization	Optimization Rate
EU1	410.0294	490.1148	19.53%	798.9328	722.6499	9.54%
EU2	380.6569	467.5588	22.83%	796.3768	719.056	9.71%
EU3	166.7465	227.3978	36.37%	443.8627	403.5896	9.07%

presents the results before and after user optimization. The study employs a benefit function to quantify user satisfaction, indicating that higher benefits are associated with greater user satisfaction.

As shown in Table 6, optimizing EU1, EU2, and EU3 results in reduced costs and increased overall benefits. Specifically, the costs for EU1, EU2, and EU3 have been reduced by 9.54%, 9.71%, and 9.07%, respectively. Therefore, the proposed model, which utilizes the KKT method to solve the bi-level game problem for microgrids and their users, yields results that facilitate demand response while enhancing user satisfaction.

• MGs Result analysis

This part establishes two comparative scenarios. Figure 11 illustrates the comparison of energy transactions under the two scenarios. Scenario 1: independent optimization, where each microgrid is optimized independently as a standalone entity without the involvement of the MGCO, and each microgrid directly engages in transactions with the external distribution network. Scenario 2: joint optimization, wherein the method proposed in this study is used to set intra-group trading prices. This allows the MGCO to centrally manage and coordinate energy transactions within the group.

From the perspective of the overall energy trading volume, the joint optimization of the microgrid cluster allowed the three proposed microgrids to engage in significant internal trading via the MGCO, resulting in a 24.78% increase in traded energy and notably decreasing the cluster’s dependence on the distribution network. This demonstrates that post-optimization, the microgrid cluster effectively balanced excess and deficit energy within the group, further reflecting and validating the experimental results shown in Figure 10. Consequently, the approach prioritized the internal absorption of surplus energy, enhancing the safety and reliability of the microgrid cluster.

7. Conclusions

This study considers the fluctuations in power generation and hydrogen load demand, proposing a game-theoretic framework method for the tri-level transaction of microgrid clusters. The conclusions are as follows:

(1) The interval number model based on polygonal uncertainty sets, which incorporates wind–solar correlation, effectively eliminates extreme scenarios and enhances the rationality of subsequent design phases and the economic viability of returns.

(2) Implementing dynamic hydrogen pricing strategies and incorporating demand response adjustments effectively ensure the minimum economic returns for hydrogen stations, broaden the interaction between energy sources and loads, and achieve proactive energy consumption of hydrogenation stations.

(3) The tri-layer method incentivizes microgrids to participate in local energy trading markets. This approach significantly improves the economic efficiency of multiple system stakeholders and reduces microgrid clusters’ reliance on the distribution network.

This paper primarily investigates a single new energy hydrogenation station. Future research will extend to examining operational scenarios for multiple hydrogenation stations at the urban scale and addressing various hydrogen load demands, including industrial applications.