Frequency Regulation of an Islanded Microgrid Using Hydrogen Energy Storage Systems: A Data-Driven Control Approach ^†

Abstract

Hydrogen energy storage (HES) systems have recently received attention due to their potential to support real-time power balancing in a power grid. This paper proposes a data-driven model predictive control (MPC) strategy for HES systems in coordination with distributed generators (DGs) in an islanded microgrid (MG). In the proposed strategy, a data-driven model of the HES system is developed to reflect interactive operations of an electrolyzer, hydrogen tank, and fuel cell, and hence the optimal power sharing with DGs is achieved to support real-time grid frequency regulation (FR). The MG-level controller cooperates with a device-level controller of the HES system that overrides the FR support based on the level of hydrogen. Small-signal analysis is used to evaluate the contribution of FR support. Simulation case studies are also carried out to verify the accuracy of the data-driven model and the proposed strategy is effective for improving the real-time MG frequency regulation compared with the conventional PI-based strategy.

1. Introduction

Variable renewable energy (VRE), such as wind and solar, is considered essential for the decarbonization of power grids. Real-time balancing between the power supply and demand is a key challenge on the pathway toward decarbonization [1]. To overcome the challenge, there has been great interest in coupling the electricity and hydrogen energy sectors [2]. Recently, hydrogen energy storage (HES) technologies have advanced to a level where HES systems can support real-time power balancing in power grids [3], and consequently mitigate the reserve requirements imposed on the generation of units and batteries.

A few recent studies have focused on the dynamic responses of HES systems. For example, in [4,5,6], individual models of HES sub-systems, such as electrolyzers, fuel cells (FCs), compressors, and hydrogen tanks, were implemented and combined to analyze the dynamic response of the whole system. However, the HES systems were modeled with a simplified first-order transfer function due to the nonlinearity of the high-order model. Moreover, the nonlinear models included many parameters, most of which need to be extracted using parameter-estimation techniques. The parameter estimations are often time consuming, or even impossible, particularly for HES systems in daily service. To resolve this issue, data-driven modeling has emerged as an alternative to physics-based modeling. Note that data-driven framework approaches are widely adopted for modeling as well as for reducing electricity-market uncertainties [7]. Data-driven modeling aims to accurately represent the input–output dynamics of a system using only statistical records of inputs and corresponding outputs. For example, [8,9] used artificial neural networks to predict the operating voltage of the FC stack, and for data-driven fault diagnosis. However, the internal states of data-driven models have no physical meaning, unlike the actual states of physics-based models. Consequently, the design of state-feedback controllers has rarely been discussed, indicating the need for further studies of stable grid-interactive HES-unit controls.

Several studies have investigated the real-time control of HES systems, based on analysis of their dynamic responses. In particular, for grid power balancing, the power outputs (or inputs) of the HES systems were adjusted using common proportional-integral (PI) controllers [4,10], fuzzy-logic controllers [5], and neural-network-based controllers [6] at the grid level under normal operating conditions. In [11], device-level PI controllers were designed under the condition where line faults occurred at the point of common coupling. However, none of the above studies considered optimal control of the HES system. In [12,13], a model predictive control (MPC) problem was formulated to determine the optimal power of the electrolyzer and HES, but MPC was achieved only for a single HES unit (hence at the device level, rather than the grid level). Further studies are still required for real-time frequency regulation (FR) of an AC grid via optimal power sharing of multiple HES and generating units.

Meanwhile, various research has been conducted on the physics-based and data-driven real-time FR on the power grid without using HES systems. In [14,15], multi-source power generation units (i.e., thermal, hydro, and gas units) are used to enhance system frequency stability with a PID controller and a tilt integral derivative controller. In [16], fuzzy-logic controllers are used for adaptive load–frequency control under different operating conditions of the power grid. For data-driven FR strategies, microgrid (MG) frequency stability is maintained using adaptive controllers [17,18], distributed controllers [19], linear-parameter-varying controllers [20], and MPC controllers [21]. However, to the author’s knowledge, data-driven real-time FR strategies using HES systems have not yet been considered. This provides motivation for studies of the real-time FR of an MG via optimal power sharing of multiple HES and generating units. Table 1 presents a summary of previous studies on the FR control methods of an MG with and without HES systems.

Table 1. Summary of previous studies on the MG frequency regulation.

Ref.	System Types (1)			Control Methods (2)
	OH	WH	Test-Bed	P	D	Controller
[4]		✓	Simplified MG	✓		PI
[5]		✓	2-area system	✓		Fuzzy
[6]		✓	2-area system	✓		Neural-network
[9]		✓	Australian test system	✓		PI
[14]	✓		2-area system	✓		PID
[15]	✓		2-area system	✓		Tilt integral derivative
[16]	✓		49-bus test system	✓		Variable universe fuzzy
[17]	✓		Canada distribution system		✓	Adaptive control
[18]	✓		10-bus MG		✓	Adaptive control
[19]	✓		Stand-alone MG		✓	Distributed control
[20]	✓		5-bus MG		✓	Linear parameter varying
[21]	✓		Simplified islanded MG		✓	MPC

⁽¹⁾OH/WH: without/with hydrogen energy storage systems, ⁽²⁾P/D: physics-based/data-driven.

Table 1. Summary of previous studies on the MG frequency regulation.

Ref.	System Types (1)			Control Methods (2)
	OH	WH	Test-Bed	P	D	Controller
[4]		✓	Simplified MG	✓		PI
[5]		✓	2-area system	✓		Fuzzy
[6]		✓	2-area system	✓		Neural-network
[9]		✓	Australian test system	✓		PI
[14]	✓		2-area system	✓		PID
[15]	✓		2-area system	✓		Tilt integral derivative
[16]	✓		49-bus test system	✓		Variable universe fuzzy
[17]	✓		Canada distribution system		✓	Adaptive control
[18]	✓		10-bus MG		✓	Adaptive control
[19]	✓		Stand-alone MG		✓	Distributed control
[20]	✓		5-bus MG		✓	Linear parameter varying
[21]	✓		Simplified islanded MG		✓	MPC

⁽¹⁾OH/WH: without/with hydrogen energy storage systems, ⁽²⁾P/D: physics-based/data-driven.

This paper presents a data-driven MPC strategy for HES systems in coordination with distributed generators (DGs) in an islanded MG. The proposed strategy aims to minimize the MG frequency deviation arising due to the intermittency of renewable power generation and load demand under normal operating conditions. For the proposed strategy, a data-driven model of the HES system is developed using the dynamic mode decomposition with control (DMDc) method. The data-driven model reflects the operation and interaction of the HES components, i.e., an electrolyzer, FC, and hydrogen tank. With the data-driven model, an MPC problem is formulated to design the optimal state-feedback controllers of the HES systems and DGs at the MG level. The MPC problem also involves a device-level controller for the HES system that overrides the MG-level controller when the level of hydrogen (LOH) is too low or high, ensuring the safe operation of the HES system. Time-domain simulations are conducted to demonstrate that the proposed data-driven control strategy is effective for suppressing the MG frequency deviation while still ensuring stable operation of the HES system under power imbalances.

The main contributions of this paper are summarized below:

• To our knowledge, this is the first study to develop a data-driven dynamic model of the HES system for integration into optimal FR support of an MG;
• Physical operating states of the HES systems are reflected in our data-driven model by using the DMDc method;
• MPC based grid-level optimal power sharing and LOH-based device-level coordination of the HES and DG units are designed;
• The effectiveness of the HES systems is carried out in small-signal analysis and simulation case studies under load variation.

2. Dynamic Model of an HES System

2.1. Fundamentals of an HES System

Figure 1 shows a schematic diagram of an HES system that mainly consists of an electrolyzer stack, hydrogen gas tank, and FC stack. Specifically, at a temperature of the electrolyzer stack (T_E), it divides water into hydrogen gas (H_E) and oxygen gas (O_E) using electrical power. An amount of H_T hydrogen gas is contained in the low-pressure gas tank at a temperature of T_T, whereas the oxygen gas is released into the ambient air. The low-pressure tank can reduce the time taken for the transition from hydrogen energy to electrical energy via chemical reactions, thereby improving the dynamic response of the HES system. For long-term energy storage, the hydrogen can be stored as a metal hydride (m_mh). The FC stack at a temperature of T_F then absorbs the hydrogen and oxygen gases from the tank and ambient air in amounts of H_F and O_F, respectively, and combines them to generate electrical power while producing water as a by-product. The HES system can operate in two different modes, which are termed here as a hydrogen gas supply (HGS) and electrical energy supply (EES) modes. In HGS mode, the electrolyzer produces hydrogen gas when the interfaced power grid has surplus power, such as in the case of high renewable power generation and low load demand. By contrast, in EES mode, the FC supplies power when renewable power is low and load demand is high. In general, the electrolyzer and FC operate independently to prevent a decrease in their lifespans and operational safety, thus enabling the HES to operate in HGS and EES modes both simultaneously and separately [12,13].

Given the structural and operational fundamentals, an HES system is modeled with state variables (i.e., X_H = [T_E, V_Fa, T_F, H_F, O_F, H_T, m_mh, T_T]^T) that significantly affect its dynamic response, as discussed in Section 2.2 and Section 2.3.

2.2. HES System Modeling

The electrolyzer stack includes unit cells connected in series, a heat exchanger, and a water pump. Considering their thermal dynamics according to cell voltage variation, the dynamic model of the electrolyzer stack is given by:

$${\dot{X}}_E(t)=f_E\left(X_E(t), U_E(t)\right),$$

(1)

where X_E = T_E and U_E = I_E. Please refer to [22,23] for details. By assigning T_E and I_E as the state and input variables, Equation (1) reflects the aggregated thermal (and hence gas) dynamics of the electrolyzer stack, rather than the individual cells, thereby facilitating grid-level application and analysis of the HES system.

Similarly, a dynamic model of the FC stack can be implemented considering the thermal dynamics and mass flow rates of hydrogen and oxygen gases with the cell voltages, as:

$${\dot{X}}_F(t)=f_F\left(X_F(t), U_F(t)\right),$$

(2)

where X_F = [V_Fa, T_F, H_F, O_F]^T, and U_F = I_F. Please see [24,25] for details. In Equation (2), a sudden change in the input current I_F leads to first-order dynamic variation in the activation voltage V_Fa due to the double-layer charging capacitance.

Moreover, the dynamic responses of the hydrogen gas and metal-hydride tanks are determined by the combined effects of the energy balance and mass and heat transfers, as:

$${\dot{X}}_T(t)=f_T\left(X_T(t), U_T(t)\right),$$

(3)

where X_T = [H_T, m_mh, T_T]^T and U_T = [H_E, H_F]^T. A detailed procedure for Equation (3) is provided in [26]. In Equation (3), H_E is greater than zero when the HES system operates in HGS mode: i.e., H_E > 0 when U_E = I_E > 0 in (1). Similarly, H_F > 0 is valid for EES mode, i.e., H_F > 0 for U_F = I_F > 0 in Equation (2).

The complete dynamic model of the HES system can be represented using the coupled nonlinear differential Equations (1)–(3), or, in an aggregated form as:

$${\dot{X}}_H(t)=f_H\left(X_H(t), U_H(t)\right),$$

(4)

where $X_H={\left[X_E^T, X_F^T, X_T^T\right]}^T={\left[T_E, V_{Fa}, T_F, H_F, O_F, H_T, m_{mh}, T_T\right]}^T$ and $U_H(t)={\left[U_E^T, U_F^T\right]}^T.$

2.3. Data-Driven Linearized Model

The physics-based nonlinear model Equation (4) includes numerous HES system parameters, rendering FR support of the HES system difficult in practice. This implies the need to develop a data-driven dynamic model of the HES system. This paper uses a DMDc method [27] to construct the underlying dynamics of a complex system from a limited number of temporal measurements. It gives the best-fit of a linearized relationship of X_H(t + Δt) with X_H(t) and U_H(t), thereby finding the best approximation of a linearized state-space representation of Equation (4).

Briefly, the states and inputs of the HES system (i.e., X_H and U_H in Equation (4)) are measured over the time period t − m·Δt to t (or, equivalently, from k − m to k in discrete time). Given the sequences of the measured data, the following matrices are constructed, as shown in Figure 2:

$${\Theta }_{\mathrm{H}}={\left[X_{H,k-m}^T, X_{H,k-(m-1)}^T,\cdots , X_{H,k-1}^T\right]}^T,$$

(5)

$$\Theta {\prime }_{\mathrm{H}}={\left[X_{H,k-(m-1)}^T, X_{H,k-(m-2)}^T,\cdots , X_{H,k}^T\right]}^T,$$

(6)

$${\Gamma }_H={\left[U_{H,k-m}^T, U_{H,k-(m-1)}^T,\cdots , U_{H,k-1}^T\right]}^T.$$

(7)

Using Equations (5)–(7), Equation (4) can then be rewritten in matrix form, as:

$${{\Theta }^{\prime }}_{\mathrm{H}}\approx { \tilde{A}}_{\mathrm{H}}{\Theta }_{\mathrm{H}}+{ \tilde{B}}_{\mathrm{H}}{\Gamma }_{\mathrm{H}}=\left[{ \tilde{A}}_{\mathrm{H}}{ \tilde{B}}_{\mathrm{H}}\right]\left[\begin{array}{l}{\Theta }_{\mathrm{H}}\\ {\Gamma }_{\mathrm{H}}\end{array}\right]=G\Omega ,$$

(8)

which leads to G = [Ã_H, ${\tilde{B}}_H$] = Θ’_H∙Ω^–1. The inverse of Ω can be calculated readily using the singular value decomposition of Ω, as shown in Equations (9) and (10):

$$\Omega =T\mathrm{\Sigma }{V}^{*}=\left[{ \tilde{T} \tilde{T}}_{\mathrm{rem}}\right]\left[\begin{array}{l}\tilde{\mathrm{\Sigma }} 0\\ 0{\mathrm{\Sigma }}_{\mathrm{rem}}\end{array}\right]\left[\begin{array}{l}{ \tilde{V}}^{*}\\ { \tilde{V}}_{\mathrm{rem}}^{*}\end{array}\right]\approx \tilde{T}\tilde{\mathrm{\Sigma }}{ \tilde{V}}^{*},$$

(9)

$$\mathrm{G}=\left[{ \tilde{A}}_{\mathrm{H}}{ \tilde{B}}_{\mathrm{H}}\right]={{\Theta }^{\prime }}_{\mathrm{H}}\left( \tilde{V}{\tilde{\mathrm{\Sigma }}}^{-1}{ \tilde{T}}^{*}\right),$$

(10)

where $T\in ℝ$^(n+q)×(n+q), $\Sigma \in ℝ$^{(n+q)×(m−1)}, $V^{*}\in ℝ$^{(m−1)×(m−1)}, $\tilde{T}\in ℝ$^(n+q)×r, $\tilde{\mathrm{\Sigma }}\in ℝ$^r×r, ${\tilde{V}}^{*}\in ℝ$^r×(m−1), $G\in ℝ$^r×(n+q), ${\tilde{A}}_H\in ℝ$^n×n, and ${\tilde{B}}_H\in ℝ$^n×q for the truncation value r of the singular values. Given ${\tilde{A}}_H$ and ${\tilde{B}}_H$ in Equation (10), the linearized dynamics of the HES system are represented in discrete time as:

$$X_{H,k+1}={\tilde{A}}_H\cdot X_{H,k}+{\tilde{B}}_H\cdot U_{H,k},$$

(11)

and, consequently, in continuous time as:

$${\dot{X}}_H(t)=A_H\cdot X_H(t)+B_H\cdot U_H(t),$$

(12)

$$\mathrm{where} A_H=\ln \left({ \tilde{A}}_H\right)/T_s \mathrm{and} B_H=A_H{\left({ \tilde{A}}_H-I\right)}^{-1}{ \tilde{B}}_H.$$

(13)

In Equations (12) and (13), A_H and B_H are calculated using only data obtained from the HES system under normal operating conditions, unlike the case of Equation (4). This facilitates FR support of HES systems of different types and sizes. Moreover, in Equation (12), X_H and U_H still represent the physical properties of the state and input variables, thereby ensuring safe operation of the HES system when used for FR support.

3. Optimal Real-Time Control of the HES System

Figure 3 shows a schematic of the optimal control of the HES systems, in coordination with DGs, to support real-time FR in an islanded MG via primary and secondary controllers (i.e., PFR and SFR). Such hierarchical control is consistent with common practice [28], ensuring the wide applicability of the proposed FR strategy. The PFR signals are produced locally when individual HES systems and DGs are connected to the MG, enabling any frequency deviations to be counteracted quickly. PFR is commonly achieved using droop controllers [29]. The signals for SFR are generated centrally and distributed via communications links so that the power outputs of all of HESs and DGs are adjusted to bring the MG frequency back to its nominal value, following an imbalance between total power generation and consumption. Considering the droop controllers, the optimal state-feedback controller for SFR is designed to enhance the coordination of the HESs and DGs, thereby improving real-time FR compared to a conventional PI-based controller. This is achieved using an MPC method with the current and predicted data obtained from the HESs and DGs as discussed in Section 3.2.

Without loss of generality, the DGs are assumed to be a type with a synchronous machine. The HES systems and DGs operate as grid-supporting and -forming units, respectively; MG voltage is regulated mainly by the DGs. Note that the proposed modeling and control methods can be readily applied to the case of an HES operating in grid-forming mode.

3.1. Extension of a Data-Driven HES Model

Specifically, the data-driven dynamic model of the HES system in (12) and (13) is extended for the case with N_h HESs, each of which is combined with models of power electronics converters serving as the interfaces with the MG (see Figure 1). Considering the PFR support with droop constants 1/R_hnj, the extended dynamic models of the HESs can be represented in an aggregated form, as:

$${\dot{X}}_{HE}(t)={\mathrm{A}}_{\mathrm{HE}}\cdot X_{\mathrm{HE}}(t)+{\mathrm{B}}_{\mathrm{HE}}\cdot U_{\mathrm{HE}}(t)+{\mathrm{B}}_{\mathrm{FH}}\cdot \mathrm{\Delta }f(t)$$

(14)

$$P_{HE}(t)={\mathrm{C}}_{\mathrm{HE}}\cdot X_{\mathrm{HE}}(t),$$

(15)

$$\mathrm{where} X_{HE}={\left[X_{HE,1}^T, X_{HE,2}^T, \cdots , X_{HE,Nh}^T\right]}^T,$$

(16)

$$X_{HE,j}={\left[X_{H,j}^T, I_{Ej}(t), {\mathrm{V}}_{EDj}(t), I_{Fj}(t), {\mathrm{V}}_{FDj}(t), I_{LCL,j}^T,{ \mathrm{V}}_{LCL,j}^T\right]}^T,$$

(17)

$$A_{HE}=di\mathrm{A}g\left({\mathrm{A}}_{\mathrm{HE},1},{ \mathrm{A}}_{\mathrm{HE},2}, \dots ,{ \mathrm{A}}_{\mathrm{HE},Nh}\right),$$

(18)

$$B_{HE}=di\mathrm{A}g\left({\mathrm{B}}_{\mathrm{HE},1},{ \mathrm{B}}_{\mathrm{HE},2}, \dots ,{ \mathrm{B}}_{\mathrm{HE},Nh}\right),$$

(19)

$$B_{FH}=di\mathrm{A}g\left({\mathrm{B}}_{\mathrm{FH},1},{ \mathrm{B}}_{\mathrm{FH},2}, \dots ,{ \mathrm{B}}_{\mathrm{FH},Nh}\right),$$

(20)

$$C_{HE}=di\mathrm{A}g\left({\mathrm{C}}_{\mathrm{HE},1},{ \mathrm{C}}_{\mathrm{HE},2}, \dots ,{ \mathrm{C}}_{\mathrm{HE},Nh}\right),$$

(21)

$$A_{HE,j}=di\mathrm{A}g\left({\mathrm{A}}_{\mathrm{H},j},{ \mathrm{A}}_{\mathrm{C},j}\right), B_{HE,j}=di\mathrm{A}g\left({\mathrm{B}}_{\mathrm{H},j},{ \mathrm{B}}_{\mathrm{C},j}\right),B_{\mathrm{FH}}=-R_{hnj}^{-1}{B}_{\mathrm{HE},j}, {\mathrm{C}}_{\mathrm{HE},j}=\left[\begin{array}{cc}{\mathrm{O}}_{8\times 1}& {\mathrm{C}}_{\mathrm{C},j}\end{array}\right].$$

(22)

Moreover, the extended HES systems Equations (14)–(22) are combined with the dynamic models of N_g DGs with droop constants R_gni and give the dynamic response model of the MG frequency [30], as:

$$M\cdot \mathrm{\Delta }\dot{f}(t)+D\cdot \mathrm{\Delta }f(t)={\displaystyle \sum _{j=1}^{N_h}{P}_{HE,j}(t)}+{\displaystyle \sum _{i=1}^{N_g}{P}_{G,i}(t)}-P_L(t),$$

(23)

where M is the moment of inertia; D is the load-damping constant; Δf is the MG frequency deviation; P_HE,j and P_G,i are the power outputs of HES unit j and DG unit i; N_h and N_g are the numbers of HESs and DGs; and P_L is the total load demand. This leads to:

$${\dot{X}}_N(t)=A_{\mathrm{N}}\cdot X_{\mathrm{N}}(t)+B_{\mathrm{N}}\cdot U_{\mathrm{N}}(t)+B_{\mathrm{L}}\cdot P_L(t),$$

(24)

$$\mathrm{where} X_N={\left[\mathrm{\Delta }f(t), {\displaystyle \int \mathrm{\Delta }f(t)}dt, X_G, X_{HE}\right]}^T,{ \mathrm{U}}_N={\left[U_G,{ \mathrm{U}}_{HE}\right]}^T,$$

(25)

$${\mathrm{A}}_{\mathrm{N}}=\left[\begin{array}{cccc}-D{M}^{-1}& 0& -{\mathrm{C}}_{\mathrm{G}}{M}^{-1}& {\mathrm{C}}_{\mathrm{HE}}{M}^{-1}\\ 0& 1& \mathrm{O}& \mathrm{O}\\ B_{\mathrm{FG}}& \mathrm{O}& {\mathrm{A}}_{\mathrm{G}}& \mathrm{O}\\ B_{\mathrm{FH}}& \mathrm{O}& \mathrm{O}& {\mathrm{A}}_{\mathrm{HE}}\end{array}\right], B_{\mathrm{N}}=\left[\begin{array}{cc}\mathrm{O}& \mathrm{O}\\ \mathrm{O}& \mathrm{O}\\ B_{\mathrm{G}}& \mathrm{O}\\ \mathrm{O}& B_{\mathrm{HE}}\end{array}\right],B_{\mathrm{L}}=\left[\begin{array}{c}\raisebox{1ex}{$-1$}\!\left/ \!\raisebox{-1ex}{$M$}\right.\\ 0\\ \mathrm{O}\\ \mathrm{O}\end{array}\right].$$

(26)

For the optimal SFR, an average observer is incorporated into Equations (24)–(26) to estimate the integral of the average frequency deviation, so that in Equation (24) X_N(t) includes the term ∫Δf(t)dt. For a sampling time of Δt = T_s, a discrete-time model of Equations (24)–(26) with the MG-level model output Y_N(k) is represented as:

$$X_N(k+1)={\mathrm{A}}_{\mathrm{Nd}}\cdot X_{\mathrm{N}}(k)+{\mathrm{B}}_{\mathrm{Nd}}\cdot U_{\mathrm{N}}(k),$$

(27)

$$Y_N(k)={\left[\mathrm{\Delta }f(k), {\displaystyle \int \mathrm{\Delta }f(k)}dk\right]}^T={ \mathrm{C}}_{\mathrm{Nd}}\cdot X_{\mathrm{N}}(k),$$

(28)

$$\mathrm{where} A_{Nd}=e^{A_{N}Ts}, B_{Nd}={\displaystyle {\int }_0^{Ts}{e}^{A_N\tau }{{\rm B}}_{{\rm N}}d\tau }, C_{Nd}=\left[\begin{array}{ccc}1& 1 & O\end{array}\right].$$

(29)

3.2. MPC-Based FR Support

Given Equations (27)–(29), the optimal SFR is achieved by solving the discrete-time quadratic optimization problem:

$$\underset{{\mathrm{U}}_{\mathrm{N}}}{\min }{\displaystyle \sum _{k=0}^{N_p-1}\left({\mathrm{Y}}_{\mathrm{N}}{(k+1)}^T{\mathrm{QY}}_{\mathrm{N}}(k+1)+{\mathrm{U}}_{\mathrm{N}}{(k)}^T{\mathrm{RU}}_{\mathrm{N}}(k)\right)}$$

(30)

$$\mathrm{subject} \mathrm{to} {\underset{\_}{U}}_N\le U_N(k)\le {\overline{U}}_N,$$

(31)

$$\underset{\_}{Y}\le C_P\cdot X_N(k)\le 0,\forall O_{hj}=1,$$

(32)

$$0\le C_P\cdot X_N(k)\le \overline{Y},\forall O_{hj}=2,$$

(33)

$$\mathrm{\Delta }\underset{\_}{Y}\le C_R\cdot \mathrm{\Delta }{X}_N(k)\le \mathrm{\Delta }\overline{Y}.$$

(34)

In Equation (30), Q is the diagonal matrix with weighting factors, each of which is multiplied by the square of each output variable in Y_N(k + 1). Similarly, R is a diagonal matrix with weighting factors for the squares of the input variables in U_N(k). Moreover, N_p represents the number of sequences over the time period from k to k + N_p where the control inputs U_N(k) are determined based on one-step-ahead prediction of Y_N(k) at each time step k. After the first control sequences are distributed to the HESs and DGs, the optimization problem is solved for the next time step k + 1. U̲_N and ${\overline{U}}_N$ are defined as the minimum and maximum variation in the SFR reference signals, respectively; Y̲ and $\overline{Y}$ are the respective minimum and maximum output powers of the HESs and DGs; and ∆Y̲ and $\mathrm{\Delta }\overline{Y}$ are the minimum and maximum ramp rate limits of the HESs and DGs, respectively. Moreover, C_P and C_R are the coefficient matrices used to extract the power outputs and ramp rates of the HESs and DGs. Note that, to develop a predictive model that captures the dynamics of line power flows in response to the power output of HESs and DGs, the proposed strategy can consider the line-limit constraints [31]. However, considering line congestion is out of the scope in this paper and has been left for future work.

Constraints (Equations (32) and (33)) reflect the device-level controllers that override the FR supports of the electrolyzer and FC stacks. In HGS mode, the power inputs of the electrolyzers are controlled to mitigate any increase in the MG frequency, for example, due to high renewable power generation and low load demand, whereas the power outputs of the FCs are set to zero. Similarly, in EES mode, the FC compensates for a decrease in MG frequency occurring due to low renewable power and high load demand, while the electrolyzer remains idle. This implies that an optimal LOH scheduling method needs to be integrated with real-time FR support to overcome the energy-storage capacity limit and reduce the occurrence of the override control, although this is not investigated further here.

3.3. Small-Signal Analysis of HES Systems

Using Equations (24)–(26), a small-signal analysis was conducted with the model parameters specified in Section 4 (see

Table 2. Parameters for the case studies.

Devices	Parameters	Values
DGs	Sn (MVA), Vn (V)	1.5, 380
	M (s), D	3, 0.1
	Tg (s)	0.1
	Tt (s)	2
	Rgni	0.1
HESs	Sn (MVA), Vn (V)	1, 380
	CE, CF (kJ/K)	162, 35
	E0 (V)	1.229
	NE, NF	60, 35
	Rhnj	0.15
Loads	SL (MVA)	2.0 + j0.8

). Figure 4a shows the eigenvalues for the case where the PFR support (i.e., increasing 1/R_gni) is provided by the HES units. All eigenvalues are placed in the left half plane (LHP), confirming MG frequency stability under normal conditions. The dominant complex-conjugate eigenvalues in Figure 4a are located farther away from the imaginary axis than those for the case with no grid-interactive HES, as shown in Figure 4b. The comparisons verify that the HES can effectively attenuate MG frequency deviation with a faster settling time and smaller oscillations. The improvement is mainly attributed to the fact that the electrolyzers and FCs can adjust their respective power inputs and outputs rapidly in response to the MG power imbalance, thereby assisting the DGs in compensating for any remaining imbalance. This could not be achieved by increasing the DG droop control gains for the case with no grid-interactive HES. Moreover, an excessive increase in 1/R_gni leads to MG frequency instability, as shown in Figure 4b.

The small-signal analysis was further conducted considering variations in the HES parameters. The HES droop control gain 1/R_hnj affected the system dynamic responses. To simplify the analysis, the electrolyzer and FC were assumed to have same values of 1/R_hnj. In Figure 5, as 1/R_hnj increases, the complex-conjugate eigenvalues move away from the imaginary axis, indicating an improvement in FR performance. When 1/R_hnj increases further, the eigenvalues change direction, implying that an excessive increase in HES control inputs is likely to act as an additional disturbance. However, the eigenvalues still remain on the LHP, ensuring frequency stability.

4. Case Studies and Simulation Results

4.1. Test System and Simulation Conditions

The proposed FR strategy was tested on the islanded MG shown in Figure 6, which was implemented using the IEEE 37-node test feeder with modifications based on [32,33,34]. A comprehensive numerical simulator of the test MG was implemented using Matlab/Simulink. Note that due to the lack of experimental data, a data-driven HES system model and data-driven MPC were developed using data measured from the physics-based HES model in the test MG. Specifically, the test MG contained three DGs and two HES systems, with total power ratings of 1.5 and 1 MVA, respectively. It also included three photovoltaic (PV) arrays with the total power rating of 1 MVA. Given the power ratings, the total load demand was set to 2.0 + j0.8 MVA. For simplicity, the respective ZIP load coefficients of all nodes were set to 1.4, –2.0, and 1.6 [35], and 3-ph balanced lines were adopted with impedances determined based on the average value over the three phases for each line configuration [33]. The intermittency of the PV arrays and loads caused an instantaneous imbalance between the power supply and demand. Table 2 lists the MG model parameters and grid- and device-level controller gains for the proposed and conventional FR strategies.

Table 3. MPC hyper-parameters for the case studies.

State Variables (Q)		Input Variables (R)
Δf (k)	∫Δf (k)dk	UG(k)	UHE1(k)	UHE2(k)
14 × 102	86 × 102	0.1	0.7	0.6

shows the weighting factors for the input and output variables in the proposed data-driven MPC strategy. The initial LOH values for HES₁ and HES₂ were set to 35% and 30%, respectively. Consequently, the weighting factors for the input variables of HES₁ and HES₂ were assumed to be 0.7 and 0.6. Different values can be assigned to the weighting factors without affecting the proposed strategy.

Table 4. Features of the FR strategies.

FR Strategies		SFR and PFR Gains
Proposed	Case 1	Set as Table 3.
Conventional	Case 2	KP = 2 and KI = 6 for SFR gains
	Case 3	KP = 15 and KI = 60 for SFR gains

shows the main features of the proposed (Case 1) and conventional (Cases 2 and 3) strategies. Cases 1 and 2 were compared to examine the effect of the optimal MPC on the real-time FR. For a fair comparison, Case 3 was also considered to compare the performance of the proposed strategy to that of the conventional strategy with an increase in SFR gain.

4.2. Validating the Data-Driven HES System

Figure 7 compares dynamic responses of the three HES models, i.e., the experimental [23,25,26] and nonlinear physics-based models discussed in Section 2.2 and the linearized data-driven model discussed in Section 2.3. Specifically, Figure 7a–c show the operating profiles of the electrolyzer, hydrogen tank, and FC, respectively, for the three different models under the test scenarios discussed in [23,25,26]. This shows that the responses of the three models were very similar for each operating profile. In particular, the profiles of the data-driven HES model are almost the same as those of the comprehensive nonlinear physics-based model. Figure 7d,e also shows the comparisons of the HES power outputs between the data-driven and physics-based models. This confirms the good accuracy of the data-driven dynamic model with respect to not only the thermal dynamics and gas reactions of the individual components but also the electrical operation of the whole system. It hence validates the accuracy of the case study results obtained in Section 4.3.

4.3. Effects of Stepwise Load Demand Variation

Figure 8a–d show f, P_{HE1, 2}, Σ_i P_G,i, and LOH_{1, 2} for the step response of the test MG according to the load variation ΔS_L = 0.1 pu at t = 3 s.

Table 5. Effects of stepwise load demand variation.

Comparison Factors	Proposed	Conventional
	(Case 1)	(Case 2)	(Case 3)
Δfpk (HZ)	0.0629	0.0823 (23.6%)	0.0835 (24.7%)
Δfrms (z)	0.0069	0.0133 (48.1%)	0.0108 (36.1%)
ΔTset (s)	4.2331	6.4674 (34.6%)	6.9254 (38.9%)

lists the corresponding numerical results of the peak-to-peak variation Δf_pk and root-mean-square variation Δf_rms. In Figure 8a, Case 1 had a significantly decreased Δf at t = 3 s, compared to Cases 2 and 3. In Case 1, Δf_pk was 23.6% and 24.7% smaller than in Cases 2 and 3, respectively. Case 1 showed a reduction in settling time T_set of 34.6% and 38.9%, respectively, compared with Cases 2 and 3. Moreover, in Case 1, f was restored to the nominal value more rapidly, with almost no overshoot, than in Cases 2 and 3. Case 1 showed reductions in Δf_rms of 48.1% and 36.1%, compared to Cases 2 and 3, respectively. This verifies that the proposed FR strategy with the data-driven HES models effectively suppressed MG frequency deviation via the fast dynamic responses of the HESs and optimal coordinated control achieved with the DGs. Figure 8b,c show that, in Case 1, P_HE increased more sharply and was maintained at a higher level than in Cases 2 and 3. In other words, the proposed strategy led to optimal power sharing between the HESs and DGs, thereby reducing Δf_pk and Δf_rms. Note that, after t = 3s sum of P_{HE1, 2} and Σ_i P_G,i were equal to 0.1 pu, which was the value of the load variation ΔS_L. This result demonstrated that the proposed method successfully updated the optimal control references of the HES and DG units, thereby improving the optimal power sharing among HES and DG units, and real-time FR over the entire period.

In all cases, the LOH₂ was reduced from 30 to 10 at approximately t = 1245 s, thereby changing the operating mode of HES₂ from FR support mode to override mode. This led to an abrupt change in the HES power output: i.e., increasing P_HE1 and decreasing P_HE2. The proposed FR strategy is still more effective at mitigating the effect of the HES mode change. This confirms that the integration of device-level control into the MPC enables the MPC problem formulated in Equations (30)–(34) to reflect the mode change successfully, thereby ensuring stable HES operation and reliable FR support. Note that the case studies considered a large increase in load demand and small initial LOH values to explicitly analyze the performance of the proposed and conventional strategies in FR support and override modes, as discussed in Section 3.

5. Conclusions

This paper proposed a data-driven MPC strategy for HES systems in coordination with DGs in an islanded MG. For the proposed MPC strategy, a data-driven model of the HES system was developed using a DMDc method, and then integrated into a linearized MG network when designing the optimal FR strategy. The accuracy of the data-driven HES model was verified by comparison with a nonlinear physics-based model and experimental dataset. Small-signal analysis was conducted to analyze the effects of the grid-interactive HES units. Given the data-driven model, the HES systems were optimally coordinated with the DGs and the proposed data-driven control strategy was effective in suppressing the MG frequency deviation, while ensuring stable operation of HES systems. The simulation results revealed that the proposed FR strategy regulates the average frequency more effectively than the conventional PI-based strategy.