Net Hydrogen Consumption Minimization of Fuel Cell Hybrid Trains Using a Time-Based Co-Optimization Model

Abstract

With increasing concerns on transportation decarbonization, fuel cell hybrid trains (FCHTs) attract many attentions due to their zero carbon emissions during operation. Since fuel cells alone cannot recover the regenerative braking energy (RBE), energy storage devices (ESDs) are commonly deployed for the recovery of RBE and provide extra traction power to improve the energy efficiency. This paper aims to minimize the net hydrogen consumption (NHC) by co-optimizing both train speed trajectory and onboard energy management using a time-based mixed integer linear programming (MILP) model. In the case with the constraints of speed limits and gradients, the NHC of co-optimization reduces by 6.4% compared to the result obtained by the sequential optimization, which optimizes train control strategies first and then the energy management. Additionally, the relationship between NHC and employed ESD capacity is studied and it is found that with the increase of ESD capacity, the NHC can be reduced by up to 30% in a typical route in urban railway transit. The study shows that ESDs play an important role for FCHTs in reducing NHC, and the proposed time-based co-optimization model can maximize the energy-saving benefits for such emerging traction systems with hybrid energy sources, including both fuel cells and ESD.

1. Introduction

According to the report published by the International Energy Agency (IEA) in 2017, CO${}_2$ emissions of transportation sector from fuel combustion reached 7.74 billion tons in 2015 [1], of which the rail sector emitted 336 million tons, accounting for 4.34% [2]. One major source of carbon in railway transportation is the trains using fossil fuel directly and indirectly. To minimize fossil fuel consumption in rail transportation, hydrogen is regarded as a feasible and effective solution. As a kind of clean and renewable energy source, hydrogen can greatly reduce carbon emission and the use of non-renewable energy when applied to transportation. As a typical study, Chang et al. established an energy consumption model based on traffic data from various parts of China, demonstrating that the use of hydrogen fuel cells in transportation represents one of the future development trends, especially in public transportation [3]. The feasibility of fuel cell hybrid trains (FCHTs) had been verified in many studies [4,5]. Although the technical cost is relatively high at the current stage, FCHTs have some advantages over electric trains and diesel trains on railways with lower traffic density [6]. In this case, optimizing train control and energy management of hybrid power systems is one of the main focuses in the study of cost reduction.

The optimal train speed trajectory with minimum traction energy demand exists when a train is running between two adjacent stations based on Pontryagin’s maximum principle (PMP) [7,8]. Based on PMP, the optimal train control modes, i.e., accelerating with maximum acceleration rate, braking with maximum braking rate, cruising, and coasting, and their optimal sequences, can be obtained for one inter-station journey. The PMP-based method obtains the optimal solution via solving the co-state function and it is referred to as the “indirect method” [9]. Another commonly used method is the direct method, which transforms the optimal control problem into a mathematical programming problem, typically represented by pseudospectral method, mixed integer linear programming (MILP), and other mathematical programming methods. In [10], both pseudospectral and MILP methods were proposed to optimize the train control strategies under fixed arrival time. The MILP method can easily improve the model to add other engineering constraints [11]. In addition, dynamic programming (DP) and various heuristic algorithms are also used to locate the optimal train control strategies. DP simplifies a complex problem by decomposing it into simpler sub-problems in a recursive manner. In [12], two heuristic algorithms, ant colony optimization and genetic algorithm, and DP were applied to search the speed trajectory, and the results show that the optimum train trajectory cannot be solved analytically, and numerical methods must be used.

While the train optimal control of FCHTs shares similar characteristics with other rail vehicles, the optimization of the speed trajectory of FCHTs can generally locate the minimum energy demand from the traction system. However, there are two important traction components of FCHT, i.e., energy storage devices (ESDs) and fuel cell, that deserve further study when the onboard energy management strategies are studied. For FCHTs, the characteristics of ESDs should be considered along with the train optimal control strategies, especially when the ESDs are charging and discharging along with the train control strategies such as braking and motoring. In [13], Snoussi et al. used sliding mode control for energy management, which well absorbed the peak power in fuel cell trains. In the past few years, many studies have focused on the the ESDs applied to electric trains to increase the recovery of regenerative braking energy and to improve the system performances. In [14], ESD is equipped on the electric train to increase the use of regenerative energy and stabilize the voltage, and the DP algorithm-based heuristic is designed to find the optimal speed trajectory. In [15], Wu et al. employed the MILP to locate the optimal train control, and the results indicated that ESD has a significant influence on the speed trajectory of electric train.

On the one hand, the characteristics of ESDs play an important role in train control; on the other hand, the nonlinear characteristics of fuel cell efficiency impose further complexity on the onboard energy management and train controls. The energy management strategy determines the power distribution between fuel cell and ESDs, and has a significant impact on the total consumption of hydrogen fuel. In recent years, there have been many methods to optimize energy management strategies, such as state machine strategy, fuzzy logic control, equivalent consumption minimization strategy, etc. [16,17]. Current studies are conducted with the fixed power–demand [18]. In [19], the convexity of the specific consumption curve is applied to improve the fuel economy. Based on a suggested power–demand curve, a scalable energy management strategy is designed. In [20], an online extremum-seeking method is used to estimate the maximum efficiency and maximum power points of the fuel cell. In [21], a fixed speed trajectory in FCHTs is divided into four states: traction, braking, coasting, and station parking, and then the power distributions between fuel cell and the supercapacitor (SC) are distributed through a multi-mode equivalent energy consumption method. The work in [22] extended to multiple fuel cells, and realized the power splits among multiple fuel cells by means of the equivalent fitting circle method, and optimized the power output of SC by means of equivalent energy consumption method. In [23], Yan et al. optimized the speed trajectory with the aim of achieving the minimum energy consumption, then determined the hybrid system control strategy based on minimum hydrogen consumption. In [24], a rule-based energy management strategy was proposed to maximize regenerative braking energy recovery. In [25], the braking process speed trajectory was obtained based on motor characteristic curve, supercapacitor capacity, maximum acceleration, and other information to ensure that the supercapacitor can obtain more regenerative braking energy, but this method does not consider the efficiency of the fuel cell. In [26], sequential optimization was applied to improve the potential of fuel efficiency, which develops a speed-smoothing strategy first then optimizes the battery charge based on the smoothed speed profile. However, the train control strategies impose direct impacts on the traction energy demand and are not included in the energy management process in the sequential optimization process. This may undermine the effectiveness of the optimization methods.

In the research area of hydrogen reduction for FCHTs, some recent papers proposed different methods in co-optimizing both train control strategies and onboard energy management concurrently. In [27,28], PMP and DP were combined to solve this problem by applying the Hamiltonian as the objective function of DP, but its application was inevitably constrained by ”the curse of dimensionality” and the boundary-value problem. In [29], Peng et al. proposed dynamic programming for co-optimization of train driving cycles and energy management for fuel cell trains. The parallelization of DP has been proposed to reduce computation time. However, when the dimension of the state variables in the dynamic programming model increases, the calculation time of the model will increase significantly due to the characteristics of the algorithm characteristics. In [30], Jibrin et al. solved the co-optimization of both energy management and speed trajectory by formulating a convex optimization model using convexity relaxation techniques which greatly improved the calculation efficiency, but convex optimization requires that every constraint in the model is convex, which affects the flexibility of the model. The above-mentioned state-of-the-art studies adopted DP and convex programming in dealing with the co-optimization, and they all promoted the advantages of applying co-optimization for FCHTs.

However, it can be seen from the above research that the existing research on co-optimization mainly adopts dynamic programming [27,28,29] and convex optimization [30]. Different from the existing methods of co-optimization, a new method is proposed in this paper to solve the co-optimization problem, and the hydrogen-saving mechanism of a fuel cell train is studied by using the co-optimization model. In addition, the energy-saving potential of regenerative braking process and enhancement of fuel cell efficiency during application by taking the advantage of the hybridization of both fuel cell and onboard ESDs needs further investigation. A detailed overview and comparison of the existing research and this paper are given in

Table 1. Overview of the existing research of FCHT energy management.

Publication	Objective	Solution Approach (es)	Energy Management Method (s)
Snoussi et al. [13]	Minimize hydrogen consumption and preserve battery life	Sliding mode control	Given speed trajectory
Kamal et al. [17]	Minimize hydrogen consumption and maintain the state of charge of the battery	A new fuel cell fuel consumption minimization strategy	Given speed trajectory
Zhang et al. [18]	Minimize hydrogen consumption	Battery state of charge (SOC) balanced strategy and dynamic programming	Given speed trajectory and power–demand curve
Peng et al. [19]	Minimize hydrogen consumption and battery loss	Adaptive Pontryagin’s minimum principle-based strategies	Given speed trajectory
Li et al. [20]	Minimize hydrogen consumption and the degradation of the stack	Online extremum-seeking	Given power–demand curve
Yan et al. [21]	Minimize hydrogen consumption	Multimode equivalent energy consumption	Given power–demand curve
Yan et al. [22]	The lowest system energy consumption	Hierarchical control method	Given power–demand curve
Yan et al. [23]	Minimize hydrogen consumption	Lagrangian algorithm	Sequential optimization
Mendoza et al. [24]	Maximize the energy recovered during braking	A rule-based energy management strategy	Given reference speed trajectory.
Li et al. [25]	Maximize regenerative braking energy	Pontryagin‘s minimum principle	Sequential optimization
Chen et al. [26]	Minimize fuel consumption	Velocity smoothing strategy	Sequential optimization
Uebel et al. [27]	Minimize fuel consumption	Dynamic programming and Pontryagin’s maximum principle	Co-optimization
Kim et al. [28]	Minimize hydrogen consumption	Dynamic programming and Pontryagin’s maximum principle	The possible optimal control state is determined first for co-optimization
Peng et al. [29]	Minimize energy consumption	Dynamic programming	Co-optimization
Jibrin et al. [30]	Minimize hydrogen consumption	Convex optimization	Co-optimization
This paper	Minimize net hydrogen consumption	Mixed Integer Linear Programming	Co-optimization

.

This paper aims to develop a time-based MILP model, which has a relatively short computational time, to find the optimal speed trajectory and energy management strategy simultaneously to minimize the net hydrogen consummations. The contributions of this paper are as follows.

(1)

A time-based co-optimization model based on MILP for FCHT is proposed to tackle the intrinsic nonlinear efficiency characteristics of fuel cell and convert the traditional energy-saving problem into a more realistic hydrogen-saving one. The train control, namely, the train speed trajectory, and energy management between two adjacent stations are optimized simultaneously to guide the autonomous driving.

(2)

The impact of capacity of ESD on the hydrogen consumption has been studied for FCHTs, and the corresponding hydrogen-saving mechanisms are explored with the quantitative simulation results, giving an insightful analysis and discussion on the mutual influence on optimal train operation.

The remainder of the paper is organized as follows: Section 2 proposes the mathematical model in detail, including motion model, energy flow model, piecewise linearization, and objective modeling. Section 3 analyzes the impact of ESD capacity on hydrogen consumption and gives the comparison between co-optimization and sequential optimization under different scenarios. Section 4 concludes the study and discusses the future research directions.

2. Methods

2.1. Motion Model

In this paper, the train is modeled as a particle. In [31], a distance-based MILP model is established to find the optimal speed trajectory of a FCHT, where $\Delta d_i$ are known parameters and $\Delta t_i$ are variables, as shown in Figure 1. However, heavy computational loads are required when calculating the output power from the fuel cells, $P_i=E_i/t_i$. It is noted that both $E_i$ and $t_i$ are both variables and the linearization of the ratio in the MILP model results in computational complexity due to the large magnitude difference and nonlinear relationship between the energy consumption and the corresponding time in the $i^{th}$ distance segment. Alternatively, this paper proposes a time-based MILP model to avoid this nonlinear relationship, where $\Delta t_i$ are known parameters and $\Delta d_i$ are variables. The train speed at each time instance is determined by the MILP model, and it is considered that the train accelerates or decelerates uniformly in each time segment. As shown in Figure 1, the speed trajectory is divided into N time segments, and the train travels through each segment with a fixed period of time, $\Delta t_i$, the length of which can vary according to different case scenarios. The distance traveled in the $i^{th}$ segment is denoted by $\Delta d_i$.

In this study, we assume that the train travels at a constant acceleration or deceleration rate during each time period; thus, $\Delta d_i$ can be calculated by (1).

$$\begin{array}{c}\hfill \Delta d_i=v_{i,ave}\Delta t\end{array}$$

(1)

where $v_{i,ave}$ is the average speed in $i^{th}$ segment denoted by (2).

When the total running time and travel distance are T and D, respectively, the number of the segments is $N=T/\Delta t$. Therefore, the total distance of the train should meet the constraint (3).

$$\begin{array}{r}\hfill v_{i,ave}=(v_{i+1}+v_i)/2\end{array}$$

(2)

$$\begin{array}{r}\hfill \sum _{i=1}^N\Delta d_i=D.\end{array}$$

(3)

When the train is running, it will be subjected to the resistance from the rotating shaft, rail, and air drag. In the model, the train resistance is modeled by the the Davis equation, Equation (4):

$$\begin{array}{c}\hfill f_{i,drag}=A+B{v}_{i,ave}+C{v}_{i,ave}^2\end{array}$$

(4)

where $f_{i,drag}$ is the drag force in the $i^{th}$ segment. A, B, and C are Davis coefficients.

To ensure safe operation of the train and the comfort of the passenger, the acceleration and deceleration cannot exceed their maximum value. The values for maximum acceleration and deceleration can follow the data from dealing with similar train operation optimization problems [10].

$$\begin{array}{c}\hfill -A_{b,max}\le (v_{i+1}-v_i)/\Delta t=a_i\le A_{a,max}\end{array}$$

(5)

where $A_{a,max}$ and $A_{b,max}$ are the maximum acceleration and deceleration rate, respectively.

The speed of the train should not exceed its speed limit, $v_{i,L}$, at each time instance $t_i$. The definition of $v_{i,L}$ is given by (29).

$$v_i\le v_{i,L}.$$

(6)

2.2. Energy Flow Model

Similar to normal train operations, there are four control strategies of FCHTs: traction, cruising, coasting, and braking. When the train is in traction or cruising, the hydrogen fuel cell and ESD provide the power to the train. During coasting, there is no need to provide energy to the train, and kinetic energy is used to overcome resistance. Because the fuel cells cannot absorb the regenerative braking energy (RBE) when the train is braking, ESD is used to recover the RBE. A simple schematic diagram of a fuel cell hybrid system is shown in Figure 2. This paper will optimize the train energy distribution strategy from the perspective of energy flow.

The conservation laws of energy in traction and braking process are modeled by (7) and (8), respectively. There are four types of energy: the traction energy or braking energy by the electric machines, the change of potential energy and kinetic energy of the train, and the energy converted to heat due to the drag force.

$$\begin{array}{c}\hfill E_{i,seg}{\eta }_m-Mg\Delta h_i-\frac{1}{2}M(v_{i+1}^2-v_i^2)-f_{i,drag}\Delta d_i\ge 0\end{array}$$

(7)

$$\begin{array}{c}\hfill E_{i,seg}/{\eta }_m-Mg\Delta h_i-\frac{1}{2}M(v_{i+1}^2-v_i^2)-f_{i,drag}\Delta d_i\ge 0\end{array}$$

(8)

where M is the total mass of the train, g is the acceleration rate due to gravity, $\Delta h_i$ is the altitude difference of $i^{th}$ segment, $E_{i,seg}$ is the traction or braking energy of the motor in the $i^{th}$ segment, and ${\eta }_m$ is the efficiency of the motor. When $E_{i,seg}$ is greater than zero, it means the train is in traction or cruising mode, and $E_{i,seg}{\eta }_m$ is smaller than $E_{i,seg}/{\eta }_m$. (7) will be applied while (8) is relaxed. When $E_{i,seg}$ is equal to zero, the train is coasting; when $E_{i,seg}$ is less than zero, the train is braking, and $E_{i,seg}{\eta }_m$ is greater than $E_{i,seg}/{\eta }_m$. (8) will be applied while (7) is relaxed. The relaxation corresponding to different train state is listed in

Table 2. Different values of $E_{i,seg}$ and their corresponding train state.

${\mathit{E}}_{\mathit{i},\mathbf{seg}}$	+	−	0
Train state	Traction, cruising	Braking	Coasting
Motoring (7)	Applied	Relaxed	Applied
Braking (8)	Relaxed	Applied	Applied

. In addition, it should be noted that during the braking mode, the mechanical braking and electrical braking can both be applied. In this case, (8) will reach “>” when mechanical braking and regenerative braking happens at the same time, transforming the kinetic energy to heat and electricity, and it will reach “=” when only regenerative braking energy happens, transforming the kinetic energy to electricity.

In order to model the physical characteristic constraints of the motor, the operation of the motor should conform to its traction and braking characteristics. During motoring, the sum of the output power of the fuel cell and ESD should be less than the maximum traction power of the motor. During braking, the braking power should not exceed the maximum braking power of the motor. These related constraints are as follows:

$$\begin{array}{c}\hfill -F_{b,max}\Delta d_i{\eta }_m\le E_{i,seg}\le F_{t,max}\Delta d_i/{\eta }_m\end{array}$$

(9)

$$\begin{array}{c}\hfill -P_{b,max}\Delta t{\eta }_m\le E_{i,seg}\le P_{t,max}\Delta t/{\eta }_m\end{array}$$

(10)

where $P_{t,max}$ and $P_{b,max}$ are the maximum traction power and maximum braking power respectively. $F_{t,max}$ and $F_{b,max}$ are maximum traction force and maximum braking force, respectively.

The energy flow during each time segment can be modeled by using a binary variable ${\lambda }_i$ (11).

$$\begin{array}{c}\hfill {\lambda }_i(E_{i,fc}+E_{i,dis}{\eta }_{ESD})-(1-{\lambda }_i)E_{i,ch}/{\eta }_{ESD}\ge E_{i,seg}\end{array}$$

(11)

$$\begin{array}{r}\hfill E_{trac}=\sum _{i=1}^N(E_{i,fc}+E_{i,dis})\end{array}$$

(12)

$$\begin{array}{r}\hfill E_{fc}=\sum _{i=1}^N{E}_{i,fc}\end{array}$$

(13)

where ${\eta }_{ESD}$ is efficiency of ESD. $E_{trac}$ is the energy required for train traction, which can be calculated by (12). $E_{i,fc}$ is the energy from hydrogen fuel cell, $E_{i,dis}$ is the energy from ESD. $E_{i,ch}$ is the energy recovered by ESD at the $i^{th}$ segment. $E_{fc}$ is the energy provided by the fuel cell throughout the operation. When the train is braking, regenerative energy is generated and charged into ESD.

${\lambda }_i$ is a binary variable; it guarantees that the charging and discharging of ESD will not occur simultaneously, and the fuel cell will not charge ESD directly, as constrained by (14)–(16). L is a sufficiently large number that it is greater than the maximum energy that the fuel cell and ESD can provide in $\Delta t$. When ${\lambda }_i=1$, ESD discharges and the fuel cell operates, and the constraints (15) and (16) are relaxed, and $E_{i,ch}$ is equal to 0. When ${\lambda }_i=0$, ESD is charged and the fuel cell stops operating, and the constraint (14) is relaxed.

$$\begin{array}{c}\hfill 0\le E_{i,ch}\le (1-{\lambda }_i)L\end{array}$$

(14)

$$\begin{array}{c}\hfill 0\le E_{i,dis}\le {\lambda }_{i}L\end{array}$$

(15)

$$\begin{array}{c}\hfill 0\le E_{i,fc}\le {\lambda }_{i}L.\end{array}$$

(16)

In this study, the state of energy (SOE) is used to denote the energy state of ESD; a similar modeling was adopted in [32], and its value is the current available energy of ESD divided by its total capacity (17). To ensure that the total charged and discharged energy does not exceed the capacity of ESD, (18) should be satisfied.

$$\begin{array}{c}\hfill SO{E}_i=\frac{E_{ini}+{\sum }_{j=1}^i{E}_{j,ch}-{\sum }_{j=1}^i{E}_{j,dis}}{E_{cap}}\end{array}$$

(17)

$$\begin{array}{r}\hfill 0\le SO{E}_i\le 1\end{array}$$

(18)

where $E_{ini}$ is the initial available energy in ESD. $E_{cap}$ is the capacity of ESD.

For the output power of the fuel cell and ESD, they should be less than the rated power, and the charging power of ESD should be less than the rated charging power. These constraints are given by

$$\begin{array}{c}\hfill E_{i,fc}\le P_{fc,max}\Delta t\end{array}$$

(19)

$$\begin{array}{c}\hfill E_{i,ch}\le P_{d,max}\Delta t\end{array}$$

(20)

$$\begin{array}{c}\hfill E_{i,dis}\le P_{c,max}\Delta t\end{array}$$

(21)

where $P_{fc,max}$ is the maximum output power of the fuel cell, and $P_{d,max}$ and $P_{c,max}$ are the maximum discharging and charging power of ESD, respectively.

2.3. Piecewise Linearization Using Special Ordered Set Type 2 (SOS2)

2.3.1. Speed-Related Variables

In the MILP model, a series of variables related to speed are included, and the relationship between them is nonlinear. Piecewise linearization (PWL) is applied to handle nonlinear relationships, which can represent a nonlinear function with a series of nonnegative variables of the special ordered set type 2 (SOS2), among which only two adjacent ones can be greater than 0, with the total sum of all variables equal to 1 [33]. To linearize the nonlinear constraints in (4) and (7), a series of ascending numbers, $C_1,C_2,\cdots ,C_K$, are defined to represent any speed within the range from $v_0$ to $v_N$. As a result, the decision variables $v_i^2$, the approximation of the speed $v_i^{{}^{\prime }}$, $v_{i,ave}^{{}^{\prime }}$, and $v_{i,ave}^{{}^{\prime }2}$, can be expressed by

$$\begin{array}{c}\hfill v_i^2=\sum _{j=1}^K{C}_j^2{\alpha }_{i,j}\end{array}$$

(22)

$$\begin{array}{c}\hfill v_i^{{}^{\prime }}=\sum _{j=1}^K{C}_j{\alpha }_{i,j}\end{array}$$

(23)

$$\begin{array}{c}\hfill v_{i,ave}^{{}^{\prime }}=\frac{v_i^{{}^{\prime }}+v_{i+1}^{{}^{\prime }}}{2}=\sum _{j=1}^K{C}_j{\beta }_{i,j}\end{array}$$

(24)

$$\begin{array}{c}\hfill v_{i,ave}^{{}^{\prime }2}=\sum _{j=1}^K{C}_j^2{\beta }_{i,j}\end{array}$$

(25)

where ${\alpha }_{i,j}$ and ${\beta }_{i,j}$ are variables of SOS2. The two variable of SOS2 should satisfy the constraints (26) and (27).

$$\begin{array}{cc}\hfill \sum _{j=1}^K{\alpha }_{i,j}=1,& \sum _{j=1}^K{\beta }_{i,j}=1\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill 0\le {\alpha }_{i,j}\le 1,& 0\le {\beta }_{i,j}\le 1.\hfill \end{array}$$

(27)

Using the preset points, e.g., $C_1,C_2,\cdots ,C_K$ and $C_1^2,C_2^2,\cdots ,C_K^2$, along any nonlinear 2D curve to present the relationship between two variables, it is possible to apply a variable of SOS2 to approximate the PWL relationship in between. In this section, we present the PWL modeling technique to approximate the quadratic relationship between the instant speed $v_i$ and its square $v_i^2$. In the following sections, our discussion will be focused on identifying the preset points, and the detailed modeling technique can be referred to the one using SOS2 variables discussed in this section.

2.3.2. Speed Limit and Altitude

For speed limits, as shown in Figure 3, the speed at each point cannot exceed the speed limit of the red curve, but speed limit generally is based on distance, so $d_i$ is used to indicate the position of the train (28) where $i=2,3,\cdots ,N$, and $d_1=\Delta d_1$.

$$\begin{array}{c}\hfill d_i=\sum _{j=1}^i\Delta d_j.\end{array}$$

(28)

By introducing PWL, the maximum speed limit at each point, $v_{i,L}$, can be formulated by

$$\begin{array}{cc}\hfill v_{i,L}=a_j^l{d}_i+b_j^l& ,d_i\in [D_j^L,D_{j+1}^L)\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill a_j^l=\frac{V_{j+1}-V_j}{D_{j+1}^L-D_j^L},& b_j^l=V_j-\frac{V_{j+1}-V_j}{D_{j+1}^L-D_j^L}{D}_j^L\hfill \end{array}$$

(30)

Let ${\alpha }_{i,j}^v$ be the SOS2 variable sets to formulate the piecewise linear relationship of the speed limit for each $d_i$; the calculation can be conducted as shown in (31).

$$v_{i,L}=\frac{V_{j+1}-V_j}{D_{j+1}^L-D_j^L}{d}_i+V_j-\frac{V_{j+1}-V_j}{D_{j+1}^L-D_j^L}{D}_j^L=\sum _{j=1}^{M-1}{D}_j{\alpha }_{i,j}^v,\sum _{j=1}^{M-1}{\alpha }_{i,j}^v=1$$

(31)

where $j=0,1,2,\cdots ,M-1$, $D_0^L=0,D_M^L=D$. The speed at the distance of D is limited by

$$\begin{array}{c}\hfill v_{N,L}=V_M.\end{array}$$

(32)

Similarly, the piecewise linear function of altitude, $f_h$, can be obtained. It is shown in Figure 3, and the gradient of $d_i$ can be calculated by

$$\begin{array}{cc}\hfill h_i=a_k^h{d}_i+b_k^h,& d_i\in [D_k^h,D_{k+1}^h)\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill a_k^h=\frac{H_{k+1}-H_k}{D_{k+1}^h-D_k^h},& b_k^h=H_k-\frac{H_{k+1}-H_k}{D_{k+1}^h-D_h^h}{D}_k^h\hfill \end{array}$$

(34)

Let ${\alpha }_{i,k}^h$ be the SOS2 variable sets to formulate the piecewise linear relationship of the altitude for each $d_i$; the calculation can be conducted as shown in (35).

$$h_i=\frac{H_{k+1}-H_k}{D_{k+1}^h-D_k^h}{d}_i+H_k-\frac{H_{k+1}-H_k}{D_{k+1}^h-D_h^h}{D}_k^h=\sum _{k=1}^{O-1}{D}_k{\alpha }_{i,k}^h,\sum _{k=1}^{O-1}{\alpha }_{i,k}^h=1$$

(35)

where $k=0,1,2,\cdots ,O-1$, $D_0^h=0,D_O^h=D$. The altitude at distance D is given as follows.

$$\begin{array}{c}\hfill h_N=H_O\end{array}$$

(36)

Then, the gradient difference can be obtained by

$$\begin{array}{c}\hfill \Delta h_i=h_{i+1}-h_i.\end{array}$$

(37)

2.3.3. Hydrogen Efficiency

In this model, hydrogen consumption rate is used to calculate the final consumption of the hybrid system. Based on the fuel-cell efficiency, the hydrogen consumption rate can be calculated by

$$\begin{array}{c}\hfill C{r}_{ifc}=\frac{P_{i,fc}}{H·{\eta }_{i,fc}}\end{array}$$

(38)

where $P_{i,fc}$ and ${\eta }_{i,fc}$ are the power output and efficiency of the fuel cell in $i^{th}$ segment, H is the combustion heat value of hydrogen, and $C{r}_{i,fc}$ is the corresponding hydrogen consumption rate whose modeling is further illustrated in Figure 4.

The hydrogen consumption rate for fuel cell is closely related to its output power, and it is normally represented with respect to the normalized power [29,34]; thus, in this proposed method, the hydrogen consumption rate can be modeled as follows:

$$\begin{array}{cc}\hfill C{r}_{i,fc}=a_r^{fc}{\theta }_{i,fc}+b_r^{fc},& {\theta }_{i,fc}\in [X_r,X_{r+1})\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill a_r^{fc}=\frac{Z_{r+1}-Z_r}{X_{r+1}-X_r},& b_r^{fc}=Z_r-\frac{Z_{r+1}-Z_r}{X_{r+1}-X_r}{X}_r\hfill \end{array}$$

(40)

Let ${\alpha }_{i,r}^C$ be the SOS2 variable sets to formulate the piecewise linear relationship of the hydrogen consumption rate; the calculation can be conducted as shown in (41).

$$C{r}_{i,fc}=\frac{Z_{r+1}-Z_r}{X_{r+1}-X_r}{\theta }_{i,fc}+Z_r-\frac{Z_{r+1}-Z_r}{X_{r+1}-X_r}{X}_r=\sum _{r=1}^{S-1}{X}_r{\alpha }_{i,r}^C,\sum _{r=1}^{S-1}{\alpha }_{i,r}^C=1$$

(41)

where $r=0,1,2,\cdots ,S-1$, $X_0=0,X_S=100$, ${\theta }_{i,fc}$ is the fuel cell normalized power which can be calculated by

$$\begin{array}{c}\hfill {\theta }_{i,fc}=\frac{P_{i,fc}}{P_{fc,max}}·100\%.\end{array}$$

(42)

When the normalized power is $X_S$, the hydrogen consumption rate is given by

$$\begin{array}{c}\hfill C{r}_{N,fc}=Z_S.\end{array}$$

(43)

In the time-based MILP model, the power of the fuel cell is obtained by

$$\begin{array}{c}\hfill P_{i,fc}=E_{i,fc}/\Delta t.\end{array}$$

(44)

2.4. Co-Optimization and Sequential Optimization Model

For a hydrogen fuel cell, the consumption of the hydrogen fuel in $i^{th}$ segment, $m_i$, can be calculated by (45). The total consumption of hydrogen fuel, $m_{fc}$, can be calculated by (46).

$$\begin{array}{cc}\hfill m_{i,fc}& =C{r}_{i,fc}\Delta t\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill m_{fc}& =\sum _{i=1}^N{m}_{i,fc}.\hfill \end{array}$$

(46)

To evaluate the value of hydrogen consumption, we apply an equivalent conversion of the energy discharged from or charged into ESD to hydrogen consumption, i.e., the hydrogen consumed by the fuel cell operating at the fixed efficiency to produce the same amount of energy discharged from ESD. In this paper, we take the maximum efficiency of the fuel cell, ${\eta }_{fc,max}$.

$$\begin{array}{c}\hfill m_{i,ESD}=\frac{P_{i,ESD}·\Delta t}{H{\eta }_{fc,max}}\end{array}$$

(47)

where $m_{i,ESD}$ is the equivalent hydrogen consumption of ESD in the $i^{th}$ segment. $P_{i,ESD}$ is the power of ESD in the $i^{th}$ segment. The value of $P_{i,ESD}$ is calculated by

$$\begin{array}{c}\hfill P_{i,ESD}=\frac{E_{i,dis}-E_{i,ch}}{\Delta t}.\end{array}$$

(48)

In the whole operation process, the RBE recovered by ESD can be calculated by (49), and its equivalent hydrogen can be calculated by (50).

$$\begin{array}{cc}\hfill E_{ch}& =\sum _{i=1}^N{E}_{i,ch}\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill m_{ch}& =\frac{E_{ch}}{H{\eta }_{fc,max}}.\hfill \end{array}$$

(50)

2.4.1. Co-Optimization

The objective function of co-optimization is applied in the co-optimization model, as shown in (51), which is to minimize the net hydrogen consumption (NHC) in the whole operation process.

$$\begin{array}{c}\hfill \begin{array}{c}\hfill NHC=\underset{v_i,P_{i,fc},P_{i,ESD}}{min}\sum _{i=1}^N(m_{i,ESD}(v_i,P_{i,ESD})\\ \hfill +m_{i,fc}(v_i,P_{i,fc}))\end{array}\end{array}$$

(51)

where $v_i$, $P_{i,fc}$, and $P_{i,ESD}$ are the decision variables. The demanded power curve is determined by $v_i$, and the energy management strategies are determined by $P_{i,fc}$ and $P_{i,ESD}$, both of which help the train operate at the minimum hydrogen consumption are obtained simultaneously.

2.4.2. Sequential Optimization

On the other hand, the objective function of the sequential optimization is introduced as follows. Compared with the control strategy of electric or fossil-fuel vehicles, hybrid electric vehicles (HEV) are challenging, since there are more control states resulted by the nonlinear characteristics of fuel cells and the charge and discharge states of ESD. Sequential optimization can reduce the computing load and memory usage. The objective equations of sequential optimization are as follows:

• Step 1: Speed trajectory optimization
  The first step is to obtain the speed trajectory with the aim of minimizing net energy consumption (NEC) of the motor. Namely, we take (52) as the objective function. Thus, it needs to be subject to constraints related to the solving speed trajectory, (1) ∼ (10) and (22) ∼ (35).

  $$\begin{array}{c}\hfill \begin{array}{c}\hfill NEC=\underset{v_i}{min}\sum _{i=1}^N{E}_{i,seg}\left(v_i\right)\end{array}\end{array}$$

  (52)
  Subject to: (1)–(10) and (22)–(35) where $v_i$ is the decision variable.
• Step 2: Energy management optimization
  The second step is to optimize the energy management with the aim of minimizing the net hydrogen consumption (NHC). In this step, the speed trajectory is obtained during the first step, determining the energy demand of the motor, $E_{i,seg}$. Thus, the second step can be described as (53).

  $$\begin{array}{c}\hfill \begin{array}{c}\hfill NHC=\underset{P_{i,fc},P_{i,ESD}}{min}\sum _{i=1}^N(m_{i,ESD}(P_{i,ESD})\\ \hfill +m_{i,fc}(P_{i,fc}))\\ \hfill \mathrm{Subject}\mathrm{to}:(11)-(21) \mathrm{and} (36)-(47)\end{array}\end{array}$$

  (53)

  where $P_{i,fc}$ and $P_{i,ESD}$ are the decision variables.

3. Results

The modeling parameters of FCHT in this paper are listed in

Table 3. Key parameters used in the FCHT model.

Symbol	Description	Value
$M\left(\mathrm{t}\right)$	Total mass of the train	72.2
$D\left(\mathrm{km}\right)$	Total travel distance	10
$T\left(\mathrm{s}\right)$	Total travel time	450
$\Delta t\left(\mathrm{s}\right)$	The time of each interval	9
$A\left(\mathrm{kN}\right)$	Davis coefficient	1.5
$B\left(\frac{\mathrm{kN}}{\mathrm{m}/\mathrm{s}}\right)$	Davis coefficient	0.006
$C\left(\frac{\mathrm{kN}}{{\mathrm{m}}^2/{\mathrm{s}}^2}\right)$	Davis coefficient	0.0067
$A_{a,max}(\mathrm{m}/{\mathrm{s}}^2)$	Maximum acceleration	1
$A_{b,max}(\mathrm{m}/{\mathrm{s}}^2)$	Maximum deceleration	−1
${\eta }_m(\%)$	Motor efficiency	0.9
${\eta }_{ESD}(\%)$	ESD efficiency	0.95
$H(\mathrm{MJ}/\mathrm{kg})$	Combustion heat value of hydrogen	140
${\eta }_{fc,max}$	Maximum efficiency of fuel cell	60%
$P_{fc,max}\left(\mathrm{kW}\right)$	Maximum power of fuel cell	250
$P_{ESD,max}\left(\mathrm{kW}\right)$	Maximum power of ESD	400
$P_{m,max}\left(\mathrm{kW}\right)$	Maximum traction power	600
$P_{b,max}\left(\mathrm{kW}\right)$	Maximum braking power	445
$F_{m,max}\left(\mathrm{kN}\right)$	Maximum traction effort	80

with references to the design proposed in [29]. The global optimizer used in this study to solve the MILP model is GUROBI9.0 [35]. In this section, the optimization results of the time-based MILP model are analyzed. In this study, the supercapacitor is applied as ESD for FCHT. The energy density of the supercapacitor is in the range of 2.5∼15 kWh/t [36], and it also has a relatively high efficiency, generally around 95% [37]. In this paper, the energy density is selected to be 10 kWh/t, i.e., 36 MJ/kg. At the same time, the converters (such as GTO and IGBT), motors, and gear system also have energy losses. As mentioned in [38], the efficiency of them are in the range of 98.5∼99.5%, 90∼94%, and 96∼98%, respectively. In this paper, the efficiency of the converters and gear system is incorporated in the motor efficiency, set as 90%.

Proton exchange membrane fuel cells (PEMFCs) are widely used for their stability and efficiency. A power–efficiency curve of a PEMFC stack is shown in Figure 5. The power–efficiency model does not consider the characteristics of current and auxiliary system in detail, but only considers the relationship between fuel cell efficiency and output power. When the output power of the fuel cell is low, the efficiency is low due to the factors related to the operation of the air compressor, cooling, and humidifying system. The efficiency decreases due to the increase of current at high load demand. It is noted that since the power–efficiency relationship shows that the lowest efficiency of the fuel cell is 20% with the output power of it being 0, when the efficiency is 20%, the fuel cell can be regarded to be in the non-working status.

3.1. The Impact of ESD Capacity on Hydrogen Consumption

In this section, the effect of ESD on the equivalent hydrogen consumption of a hybrid power system is studied based on the co-optimization model. Two case studies, i.e., Case 1 and Case 2, are proposed to model the FCHT with and without ESDs, respectively. In Case 1, the capacity of ESD is assumed to be 0 MJ to indicate that the FCHT is not equipped with ESD, and ESD capacity in Case 2 is set to be 40 MJ. Comparison of results between the two cases is shown in

Table 4. Key parameter settings and optimization results for Case 1 and Case 2.

	Case 1	Case 2
$Capacity\left(\mathrm{MJ}\right)$	0	40
$E_{trac}\left(\mathrm{MJ}\right)$	78.25	75.78
$E_{fc}\left(\mathrm{MJ}\right)$	78.25	36.27
$m_{fc}\left(\mathrm{g}\right)$	1067.66	418.54
$E_{ch}\left(\mathrm{MJ}\right)$	/	11.51
$m_{ch}\left(\mathrm{g}\right)$	/	137.02
$NHC\left(\mathrm{g}\right)$	1067.66	751.87

. The initial SOE for ESD in Case 2 is 1.

In Figure 6a, the optimized speed trajectories are demonstrated. It can be seen that the speed trajectory varies when ESD is introduced. The acceleration of Case 1 without ESD is smaller, so a longer acceleration phase and higher coasting speed are required to meet the time constraint and this leads to more traction energy needed for Case 1.

Power splits between fuel cell and ESD are shown in Figure 6b, where the solid line represents the optimization results when ESD capacity is 40 MJ, and the dashed line represents the results without ESD. It can be seen from this figure that the fuel cell in Case 2 outputs $100\mathrm{kW}$ for all working hours except for 150 kW for 10 to 20 s. ESD provides extra power so that fuel cell in Case 2 can work with a higher efficiency. In addition, in Case 1, the RBE is consumed directly as heat while the RBE is recovered by 11.51 MJ, which allows 12.83% equivalent hydrogen consumption to be saved in Case 2.

It can be seen that the net hydrogen consumption of the two is significantly different. The net hydrogen consumption of Case 2 can obtain a saving of 29.58% compared with Case 1. This comparison shows that the introduction of ESD can reduce net hydrogen consumption from two aspects. On the one hand, ESD assists fuel cells to avoid high-power but low-efficient working states; on the other hand, RBE can be recovered to improve the energy efficiency.

To show the detailed mechanism on how the fuel cell and ESD influence each other, the relationship between fuel cell and ESD capacity is studied for this hybrid power system. Optimization results with the varied ESD capacities by using the proposed model are shown in Figure 7a,b. Note that the initial SOE for ESD is 1.

As can be seen from Figure 7a, as ESD capacity increases, net hydrogen consumption decreases, and this trend is particularly evident in the low capacity, from 0 to 13 MJ. At this stage, the equivalent hydrogen consumption of ESD is 0. It can be seen from Figure 7b that the reason is that the energy provided by ESD is equal to the energy recovered. This indicates that at low capacity, the increase in recovered regenerative braking energy can significantly reduce net hydrogen consumption.

When ESD capacity increases from 13 to 40 MJ, the net hydrogen consumption of the train decreases slowly, while the equivalent hydrogen consumption of ESD increases and the hydrogen consumption of the fuel cell decreases. From the perspective of energy flow, the hybrid power system is more inclined to use the energy in ESD due to the higher energy efficiency of ESD. This is reflected in the fact that with the increase of ESD capacity, ESD provides more energy for traction, which allows the fuel cell to work at a higher efficiency and obtain a higher average efficiency. In addition, although the net hydrogen consumption decreases slowly, it can be seen from Figure 7b that the demanded traction energy of the train increases slowly, while the RBE also slightly decreases. This also reflects the higher energy efficiency of the hybrid power system.

It is noted that when ESD capacity is greater than 40 MJ, the net hydrogen consumption remains unchanged, and it can be seen from Figure 7a that the fuel cell average efficiency is 60%, which is the maximum value. Since the efficiency of ESD in this model is a constant, and in this model, the maximum efficiency of the fuel cell, which is 60%, is used to equate the energy provided by ESD to hydrogen consumption, the energy conversion rate of the system reaches the highest, and the net hydrogen consumption of the system is the minimum. Therefore, when the objective is to minimize hydrogen consumption, 40 MJ can be considered as the optimal ESD capacity of the system.

3.2. Case Study on Flat Track

Section 3.2 and Section 3.3 are the comparison between co-optimization and sequential optimization. In this section, a simple case is conducted to explore the difference between the two methods on a flat track to avoid the effect of gradient and speed limit on traction or braking power. The capacity of ESD in this case is selected as 40 MJ and its initial SOE is 1, and the travel distance and time are 10 km and 450 s. The performance comparison between the two is listed in

Table 5. The performance comparison between co-optimization and sequential optimization without speed limit and gradient.

	Case 2	Case 3
Computation time(s)	6.65	step 1	1.94
		step 2	0.19
$NEC\left(\mathrm{MJ}\right)$	64.27	62.28
$E_{fc}\left(\mathrm{MJ}\right)$	36.27	34.58
$m_{fc}\left(\mathrm{g}\right)$	418.54	446.11
$NHC\left(\mathrm{g}\right)$	751.87	775.87

, where Case 2 and Case 3 are based on co-optimization and sequential optimization, respectively.

The speed trajectory and motor power are shown in Figure 8a. The solid line and the dashed line represent the results of co-optimization and sequential optimization, respectively. It can be seen that the speed trajectory of sequential optimization can be divided into four stages, acceleration with maximum acceleration, cruising, coasting, and deceleration with maximum deceleration, which meets the optimal control theory, and it should be noted that the train of Case 3 which is based on sequential optimization accelerates to the cruise speed with the maximum traction power, while Case 2 does not. Figure 8b represents the power output of the fuel cell and ESD. From this figure, it can be seen that ESD of Case 2 and Case 3 have reached the maximum power output in the acceleration stage, and the power difference at this stage is only reflected in the output power of the fuel cell.

In Figure 8b, it is shown that in the acceleration stage, the fuel cell efficiency of Case 3 is 55%, but for Case 2, except for the fuel cell output power of 11 s to 22 s, which is 150 kW when the efficiency of fuel cell is 58%, the output power of the other acceleration time is 100 kW when the efficiency of fuel cell is 60%. That is to say, the fuel cell in Case 2 worked at the point of higher efficiency during the whole acceleration process. As shown in Table 5, even though the fuel cell of Case 3 provides less traction energy, the hydrogen consumption of the fuel cell is still higher, which means that the average efficiency of the fuel cell in Case 3 is lower than the one in Case 2.

It can be seen from Table 5 that the calculation time of sequential optimization is short and it has an advantage in terms of net energy consumption, i.e., 62.28 MJ. In the meantime, the result by co-optimization consumes 64.27 MJ of net energy. However, since the efficiency of the fuel cell is improved during the operation obtained by co-optimization, the net hydrogen consumption in the co-optimization results is lower, saving 3.1% compared to the one in Case 3. The sequential optimization first optimizes the speed trajectory with the minimum energy consumption, so the NEC is less than co-optimization. In terms of hydrogen consumption, the nonlinear characteristics of fuel cell efficiency are considered during optimization in the co-optimization model proposed in this paper. This method allows energy management to actively change the speed trajectory accordingly, thereby changing the corresponding power–demand curve to further improve the energy conversion efficiency of the fuel cell. As a result, the fuel cell efficiency in the operation obtained by the co-optimization result is higher than the one in the results of sequential optimization.

3.3. Case Study with Speed Limits and Gradients

The above cases are carried out on the assumption that there is no speed limit and gradient in the journey. However, in the real-world railway operations, there are speed limits and gradients along the journey. In this case, ESD capacity, travel distance, and time were 40 MJ, 10 km, and 450 s, respectively. The initial SOE of ESD is 1. Case 4 and Case 5 indicate co-optimization and sequential optimization, respectively. The geographic information is shown in Figure 9a; a speed limit section, an uphill area, and a downhill area are added.

With extra constraints of gradients and speed limits introduced, from the speed trajectory profiles presented in Figure 9b, it can be seen that the train cruises at 72 km/h in 196 s to 450 s in the speed limit area, and the train of sequential optimization is still accelerating with maximum acceleration, which requires higher traction power. There are still some differences in the two acceleration phases; the traction power of the train in Case 5, from 0 to 62 s and 252 s to 285 s, is 600 kW, which means that the train has the maximum acceleration.

As shown in Figure 9c, the efficiency in Case 4 in general is higher than the one in Case 5 due to a better energy management on the train load demand to enable a higher working efficiency of the fuel cell. As can be seen from

Table 6. The performance comparison between co-optimization (Case 4) and sequential optimization (Case 5) with speed limit and gradients.

	Case 4	Case 5
Computation time(s)	12.71	step 1	4.57
		step 2	0.06
$NEC\left(\mathrm{MJ}\right)$	66.57	64.52
$E_{fc}\left(\mathrm{MJ}\right)$	39.28	37.15
$m_{fc}\left(\mathrm{g}\right)$	459.47	512.99
$NHC\left(\mathrm{g}\right)$	784.48	838.85

, the net hydrogen consumption in Case 4, i.e., co-optimization, is reduced by 6.48% compared with that in Case 5, i.e., sequential optimization.

Different from the case study on the flat railway track, the speed trajectories of the two methods are different in the cruising phase. It can be seen from the discharge of ESD of Case 4 in Figure 9d that ESD discharges all the available energy during the acceleration process of the train, so the train in Case 4 is not running at a constant cruise speed to avoid the high fuel cell power. From the energy-efficiency point of view, the train needs higher power during acceleration, so ESD is more inclined to discharge during the acceleration process to increase the train’s speed, and this will support the fuel cell to work in states with higher efficiency.

It can be seen from the above case studies that the power split strategy of both the fuel cell and ESD is determined by the algorithm itself but not preset by rules. In the model, only the total energy demand or generation of the train is considered, and by conducting the optimization, the power and energy split of the fuel cell and the supercapacitor can be configured by the model automatically to achieve the objective of minimizing the hydrogen consumption.

4. Conclusions and Future Work

In this paper, we establish a time-based MILP co-optimization model to optimize the train control and energy management simultaneously to minimize the net hydrogen consumption. Based on this model, two problems are discussed. The first is the relationship between ESD capacity and net hydrogen consumption in a hybrid power system. The results show that ESD can reduce net hydrogen consumption by up to 29.58% in the proposed case study. It provides extra power to enable the fuel cell to work at the point of higher efficiency, leading to a significant reduction of net hydrogen consumption.

Another problem is the difference between co-optimization and sequential optimization. Two specific scenarios are simulated to compare the net hydrogen consumption, and the results show that co-optimization can save 6.48% net hydrogen consumption in the proposed case study. The advantage of co-optimization is that it tends to avoid the working states of high power and low efficiency for a fuel cell, which mainly exist during the acceleration phase, and the power–demand curve can be optimized simultaneously to find a global optimal solution.

In this paper, the optimal ESD capacity is obtained based on the minimum hydrogen consumption, without considering the device initial cost and other engineering factors, and the results may not fully reflect the economic indicators of the railway line operation. In addition, the whole life-cycle optimization of fuel cell and ESD is not considered in the proposed model, and the multi-station operation, considering optimal running time allocation, can also be explored in the future. Therefore, in future, the cost and detailed characteristics of ESD and the whole life-cycle cost can be further studied in the model to achieve a more realistic optimization result, and flexible running time can be considered. It also should be noted that the linearization technique we used in the paper is to transfer the nonlinear optimization problem into an MILP problem to achieve a more efficient method of finding the global optimal solution. However, the approximation applied in the modeling process brings the errors, and the evaluation of the error and the comparison between the linearized and nonlinear method on the same problem is also the potential future direction from the perspective of methodology. The size of voltage of the energy system is also an important variable in the train control, which is of great significance to implement the energy management strategy of the train and ensure that the train runs at the optimal speed trajectory. Therefore, in the future work, the voltage will be considered in the model.