Energy Management of Hydrogen Hybrid Electric Vehicles – Online-Capable Control

1. Introduction

A major challenge of today’s society is to decrease the worldwide CO${ }_2$ emissions, of which a major part stems from the transportation sector. In 2022, it was accountable for almost one-fourth of the global CO${ }_2$ emissions [1]. Hence, a major part of the global decarbonization politics includes ever more stringent emissions legislations in this sector [2,3]. This signifies an extension of technological choices for the light duty market in the EU. However, currently, the traditional combustion engine is still the predominant propulsion technology. In 2023, the global electric light-vehicle sales accounted for roughly 16% of new car sales [4]. Therefore, the need to find alternatives to the conventional combustion engine incentivizes investments into advanced propulsion technologies that can comply with the current tank-to-wheel CO${ }_2$ goal and offer the potential for well-to-wheel and life cycle net zero CO${ }_2$ emissions.

1.1. Context and Motivation

Hydrogen is the simplest energy carrier for vehicles, produced from electricity generated by renewable power plants [5]. It can be utilized in two ways: First, hydrogen can be used in fuel cells to generate an electric current, which is used to power an electric motor or stored in a battery. Second, hydrogen can be burned in a hydrogen internal combustion engine (H${ }_2$ICE), which directly provides mechanical power. Both technologies have advantages and disadvantages: The fuel cell is known to emit only water and reaches a high efficiency [6], but requires a high degree of hydrogen purity [7]. The H${ }_2$ICE can burn hydrogen, even of low purity, but produces NO${ }_{\mathrm{x}}$ emissions in the process. In this text, the focus is put on the H${ }_2$ICE and how to tackle the engine-out NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions.

Strategies to operate H${ }_2$ICEs with extremely low NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions usually rely on ultra-lean hydrogen combustion [8,9], which is achieved using adequate turbocharging [10,11]. However, ultra-lean combustion limits the achievable power density of H2ICEs, thus limiting its applicability in mobile application. One possible solution for this issue is the use of a secondary torque source, such as an electric motor, to assist the H${ }_2$ICE during peak load instances. The result is a H${ }_2$-HEV with two onboard energy storage devices, i.e., hydrogen and a battery. This also introduces an additional degree of freedom, which is typically managed by an energy management system (EMS) algorithm. In the following literature review on such algorithms is provided. Particular emphasis is placed on EMS algorithms that explicitly consider NO${ }_{\mathrm{x}}$ emissions.

1.2. Literature Review

The EMS algorithm is concerned with the distribution of the driver’s power request between the vehicle’s different propulsion components. The presented algorithms can be divided into two categories: Offline algorithms, which require large amounts of computational power and are used to evaluate the theoretical potential of a system; and online-capable controllers, which are computationally much leaner and can thus be implemented on a real vehicle. In this text, the focus lies on online-capable EMS controllers.

The earliest EMS controllers reported in the literature are based on simple heuristics (e.g. [12]). Their main goal is to minimize the vehicle’s fuel consumption. However, the power distribution not only affects fuel consumption but also has a widespread impact on a vehicle’s powertrain components. While more and more control objectives can be incorporated into an appropriately designed EMS, this naturally leads to an increase in complexity. Also, the inclusion of predictive data to plan the vehicle’s battery usage for the upcoming driving mission was studied [13]. The integration of thermal effects in the power distribution decision process is investigated in [14] and [15]. The authors of [16] show that, in the case of a Diesel HEV, the decision of the EMS also has a large influence on its NO${ }_{\mathrm{x}}$ emissions. Compared to HEVs powered by traditional fuels, the literature on EMS algorithms for H${ }_2$-HEVs is scarce, especially with a focus on NO${ }_{\mathrm{x}}$ emissions.

The authors of [17] investigate the NO${ }_{\mathrm{x}}$ emissions resulting of simulations of a parallel H${ }_2$-HEV on the world-wide harmonized light duty test cycle. However, they lack the use of an advanced EMS control algorithm as their main focus lies on the investigation of the propulsion architecture. The only study, which discusses the connection between the EMS and the H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off of a H${ }_2$-HEV is presented in [18]. Here, a strong emphasis is laid on the potential of H${ }_2$-HEVs to reduce the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions via dedicated use of the system’s degrees of freedom to satisfy the driver’s power request. The discussed degrees of freedom are managed by the EMS and include, among others, the power distribution between the H${ }_2$ICE and the electric motor, and the driving mode. The presented study shows that, given an optimal use of the EMS algorithm, reduced NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions can only be achieved at an increased amount of H${ }_2$ consumption. Such a trade-off is described using a so-called Pareto curve.

Figure 1 shows a slightly adapted version of a Pareto curve, which is presented in [18]. It depicts the optimal normalized H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off of a mixed H${ }_2$-HEV on a specific driving mission. Each point on this curve represents one optimal EMS calibration. This means that for each point on this curve, accumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions cannot be decreased without increasing H${ }_2$ consumption, and vice-versa. Points above the curve are called Pareto subobtimal, points below the curve cannot be reached. The H${ }_2$-optimal EMS calibration is depicted by the point the farthest to the right on the Pareto curve. The NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-optimal EMS calibration is depicted by the point the farthest to the left on the Pareto curve. As the DP algorithm is used, the results represent the global optimum. However, DP cannot be practically implemented on a real vehicle due to its immense requirements for computational power and computation time.

Three conclusions are drawn in the presented study.

a) 

The potential for NO${ }_{\mathrm{x}}$ reduction is much larger for H${ }_2$-HEVs than for conventional Diesel-powered HEVs, as-well-as standard H${ }_2$ vehicles.

b) 

Although ultra-lean combustion of hydrogen-air mixtures allows H${ }_2$ICEs to emit near-zero NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions, this is a highly delicate process. Small deviations from the chosen OP of the H${ }_2$ICE can increase the instantaneous NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions by over two orders of magnitude.

c) 

The mixed hybrid drivetrain architecture is required to achieve consistent NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ reductions across a wide range of challenging driving missions. However, it is more complex than the standard parallel or series hybrid architectures.

Point b) suggests the use of an online optimization-based control algorithm to balance the delicate H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off. Point c), however, complicates the use of online optimization-based control algorithms, as the mixed hybrid drivetrain architecture introduces integer control variables resulting in a mixed-integer optimization problem. Currently, there is no online-capable controller in the literature that achieves similar performance to what is predicted by the offline computed potential analysis presented in [18].

1.3. Research Statement

This paper aims at closing this gap. A MPC algorithm is used for the EMS algorithm, since the incorporation of conflicting control goals, as for instance the delicate H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off, is straight-forward in optimization problems. A careful separation of the integer and the continuous control variables allows for the circumvention of solving a mixed-integer optimization problem online. Summarizing, the contribution of this publication is twofold:

• To the authors’ best knowledge this publication presents the first online-capable EMS controller for a H${ }_2$-HEV, explicitly accounting for the H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off.
• A case study, using the same mixed H${ }_2$-HEV as discussed in [18], allows for a comparison between the proposed online-capable EMS controller and the full theoretically reachable Pareto front obtained by the DP algorithm. The results show that the proposed online-capable controller reaches close-to-optimal performance on all investigated driving missions, covering a broad range of driving scenarios.

1.4. Paper Structure

The remainder of this paper is organized as follows: In Section 2, the investigated mixed H${ }_2$-HEV is introduced. Moreover, a simplified mathematical model is presented, which is later leveraged to formulate an optimization problem suitable for online-capable control. In Section 3, a method to avoid solving a mixed-integer program is presented. Moreover, a convex optimization problem is presented which can be solved quickly. In Section 4, the proposed EMS controller structure is presented. In Section 5, the EMS algorithm is compared to the DP-optimal solution in a case study encompassing in total four different driving missions. In Section 6, the main contributions are summarized and an outlook on future research is presented.

2. Modeling

In this section, first, the investigated vehicle is outlined and a sketch of the EMS optimization problem including NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ is presented. Second, a motivation for the MPC-based controller structure is given, which then requires some model simplifications to ensure the online-capable use of the MPC’s optimization problem.

2.1. Map-Based Powertrain Model

Figure 2 depicts a schematic representation of the vehicle’s powertrain discussed in this text. The component models are characterized using maps created from testbench measurements, which is why it will be referred to as “map-based" model from here on. The maps for the H${ }_2$ICE and the motor and generator are given in the following, however, a detailed description of the map-based model is provided in [18].

$$\begin{array}{cc}\hfill {\dot{m}}_{\mathrm{f}}^{\mathrm{map}}& =f({\omega }_e,T_e),\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\dot{m}}_{\mathrm{NO}}^{\mathrm{map}}\mathrm{x}& =f({\omega }_e,T_e),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill P_{\mathrm{l}}^{\mathrm{map}}\mathrm{mot}/\mathrm{gen}& =f({\omega }_{mot/gen},T_{mot/gen}).\hfill \end{array}$$

(3)

Here, ${\omega }_e$ denotes the engine speed, $T_e$ denotes the engine torque, and ${\omega }_{mot/gen},T_{mot/gen}$ denote the speed and torque of the motor, or the generator, respectively. The battery state-of-charge is described by the variable $\mathrm{SoC}$ and is modeled using a Thevenin equivalent circuit model.

The vehicle under investigation is a so-called mixed H${ }_2$-HEV. It is equipped with one electric motor, which provides propulsive power and is used for regenerative braking, a H${ }_2$ICE, which is connected to the wheels via a clutch and a fixed transmission ratio, and a generator. The H${ }_2$ICE and the generator together form the engine-generator unit. A mixed hybrid electric vehicle can be operated in three different driving modes, which are introduced in the following.

Driving modes:

The driving modes are denoted by $M\in \{1,2,3\}$ in this work and are listed in Table 1. In series mode, the H${ }_2$ICE is operating, but it is mechanically decoupled from the wheels. Any requested power is fully delivered by the motor and the H${ }_2$ICE is only used to generate electric power through the generator. In parallel mode, the H${ }_2$ICE is operating and is mechanically coupled to the wheels. In this mode, both propulsion systems (H${ }_2$ICE and motor) can provide power to the wheels. In EV mode, the H${ }_2$ICE is off and mechanically decoupled from the wheels. The entire torque request is met by the motor.

Rotational speeds:

Depending on the driving mode, the rotational speeds of the individual powertrain components can be calculated. From the vehicle speed v and the wheel radius $r_w$, the wheel speed calculates as:

$${\omega }_w=\frac{v}{r_w}.$$

(4)

As the motor is always mechanically coupled to the wheels, its speed is directly defined by:

$${\omega }_{\mathrm{mot}}={\omega }_w·{\gamma }_{\mathrm{fd}}·{\gamma }_{\mathrm{mot}},$$

(5)

where ${\gamma }_{\mathrm{fd}}$ and ${\gamma }_{\mathrm{mot}}$ are the fixed gear ratios of the final drive and the motor, respectively. Although, the engine is not always mechanically coupled to the wheels, it is always mechanically coupled to the generator. Therefore, ${\omega }_{\mathrm{e}}$ is defined by the generator speed ${\omega }_{\mathrm{gen}}$:

$${\omega }_{\mathrm{e}}=\frac{{\omega }_{\mathrm{gen}}}{{\gamma }_{\mathrm{gen}}},$$

(6)

where ${\gamma }_{\mathrm{gen}}$ is the fixed gear ratio of the generator.

The generator speed depends on the driving mode. In parallel mode, ${\omega }_{\mathrm{gen}}$ is fixed by the vehicle speed. In EV mode, the generator is not used and its speed is zero. In series mode, the generator is mechanically decoupled from the wheels, but its speed can be chosen freely as it does not depend on the vehicle speed.

$${\omega }_{\mathrm{gen}}\left\{\begin{array}{cc}\in [{\omega }_{\mathrm{gen}}^{\min },{\omega }_{\mathrm{gen}}^{\max }],\hfill & \mathrm{if}M=1,\hfill \\ ={\omega }_w·{\gamma }_{\mathrm{fd}}·{\gamma }_{\mathrm{gen}},\hfill & \mathrm{if}M=2,\hfill \\ =0,\hfill & \mathrm{if}M=3.\hfill \end{array}\right.$$

(7)

Here, ${\omega }_{\mathrm{gen}}^{\min }$ and ${\omega }_{\mathrm{gen}}^{\max }$ represent physical limits imposed by the engine idle speed and the maximum engine speed.

EMS including NO${ }_{\mathrm{x}}^{\mathrm{eo}}$

Using the map-based model and the description of the rotational speeds, the following sketch of the EMS optimization problem including NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ is stated:

$$\begin{array}{cc}\hfill \underset{T_e,T_{\mathrm{mot}},T_{\mathrm{gen}},{\omega }_{\mathrm{gen}},M}{\min }& {\int }_{t_0}^{t_f}{\dot{m}}_f^{\mathrm{map}}dt\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.\frac{d\mathrm{SoC}}{dt}& =f(\mathrm{SoC},I_b)\hfill \end{array}$$

(8b)

$$\begin{array}{cc}\hfill \mathrm{SoC}\left(t_0\right)& ={\mathrm{SoC}}_0\hfill \end{array}$$

(8c)

$$\begin{array}{cc}\hfill \mathrm{SoC}\left(t_f\right)& ={\mathrm{SoC}}_0\hfill \end{array}$$

(8d)

$$\begin{array}{cc}\hfill \frac{d{m}_{\mathrm{NO}}\mathrm{x}}{dt}& ={\dot{m}}_{\mathrm{NO}}^{\mathrm{map}}\mathrm{x}\hfill \end{array}$$

(8e)

$$\begin{array}{cc}\hfill m_{\mathrm{NO}}\mathrm{x}\left(t_0\right)& =0\hfill \end{array}$$

(8f)

$$\begin{array}{cc}\hfill m_{\mathrm{NO}}\mathrm{x}\left(t_f\right)& ={\overline{m}}_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{target}}\hfill \end{array}$$

(8g)

$$\begin{array}{cc}\hfill T_{\mathrm{req}}\mathrm{GB}& =\left\{\begin{array}{cc}{T}_{\mathrm{mot}}\mathrm{GB},\hfill & \mathrm{if}M=1,\hfill \\ T_{\mathrm{e}}\mathrm{GB}+T_{\mathrm{mot}}\mathrm{GB}+T_{\mathrm{gen}}\mathrm{GB},\hfill & \mathrm{if}M\in \{2,3\}\hfill \end{array}\right.\hfill \end{array}$$

(8h)

$$\begin{array}{cc}\hfill T_{\mathrm{gen}}\mathrm{GB}& =T_{\mathrm{e}}\mathrm{GB},\mathrm{if}M=1,\hfill \end{array}$$

(8i)

where, $I_b$ denotes the electrical current flowing from the battery, and the subscript ${ }_{\mathrm{GB}}$ denotes the evaluation of the corresponding torques at the level of the gearbox (introduced further down in the text). The goal of this optimization problem is to minimize the cumulated hydrogen consumption over the driving mission, which starts at time instance $t_0$ and ends at time instance $t_f$. Additionally, a cumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ target has to be satisfied at the end of the driving mission. Constraint (8h) ensures that the torque request $T_{\mathrm{req}}$ of the driver must be met by the sum of all the drivetrain’s components. Note that this optimization problem is a multi-objective, nonlinear, mixed-integer program and as such, in general, notoriously hard to solve. One time and resource-consuming solution approach based on the DP algorithm is presented in [18]. However, the aim of this publication is to derive an online-capable controller, the DP algorithm cannot be used for the energy management due to its low computational efficiency.

One control method, which is well-known to cope with multi-objective goals and nonlinear system dynamics, is based on the model predictive control scheme. Therefore, such a controller architecture is aimed for. However, as the technique of model predictive control relies on solving an optimization problem. Therefore, the formulation of a simplified version of the sketch presented in Equation (Section 2.1), which allows a fast and efficient numerical optimization, is needed. One possible ally is the family of convex optimization problems, as there exist highly efficient solvers to compute their solution quickly. In order to formulate a convex optimization problem, two steps must be taken: First, a simplified model needs to be introduced, which allows the formulation of a convex optimization problem. Second, the mixed-integer nature of the optimal control problem in Equation (Section 2.1) needs to be circumvented by dropping the driving mode from the control space.

2.2. Simplified Powertrain Model

The goal with the simplified model structure is to approximate the map-based model with the use of convex functions. Figure 3 depicts the power flows in the powertrain. The vehicle’s differential and the final drive are lumped together forming the gearbox, which has a constant efficiency ${\eta }_{\mathrm{GB}}$. The power request after the gearbox is defined as:

$$P_{\mathrm{req}}\mathrm{GB}=\left\{\begin{array}{cc}{P}_{\mathrm{req}}÷{\eta }_{\mathrm{GB}},\hfill & \mathrm{if}{P}_{\mathrm{req}}\ge 0\hfill \\ P_{\mathrm{req}}·{\eta }_{\mathrm{GB}},\hfill & \mathrm{if}{P}_{\mathrm{req}}<0.\hfill \end{array}\right.$$

(9)

The gearboxes of the motor and the generator have the constant efficiencies ${\eta }_{\mathrm{mot}}$ and ${\eta }_{\mathrm{gen}}$. Based on these, the motor power and the generator power at the gearbox-level are calculated as follows:

$$\begin{array}{cc}\hfill P_{\mathrm{mot}}\mathrm{GB}& =\left\{\begin{array}{cc}{P}_{\mathrm{mot}}·{\eta }_{\mathrm{mot}},\hfill & \mathrm{if}{P}_{\mathrm{mot}}\ge 0,\hfill \\ P_{\mathrm{mot}}÷{\eta }_{\mathrm{mot}},\hfill & \mathrm{if}{P}_{\mathrm{mot}}<0,\hfill \end{array}\right.\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill P_{\mathrm{gen}}\mathrm{GB}& =P_{\mathrm{gen}}·{\eta }_{\mathrm{gen}},\hfill \end{array}$$

(11)

where $P_{\mathrm{mot}}$ and $P_{\mathrm{gen}}$ are the motor’s and generator’s mechanical power. The motor losses and generator losses are denoted by $P_{\mathrm{l}}\mathrm{mot}$ and $P_{\mathrm{l}}\mathrm{gen}$. The sum of the electrical components’ mechanical power and its losses is referred to as their source power:

$$\begin{array}{cc}\hfill P_{{\mathrm{s}}_{\mathrm{mot}}}& =P_{\mathrm{mot}}+P_{\mathrm{l}}\mathrm{mot},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill P_{{\mathrm{s}}_{\mathrm{gen}}}& =P_{\mathrm{gen}}+P_{\mathrm{l}}\mathrm{gen}.\hfill \end{array}$$

(13)

The battery source power $P_{{\mathrm{s}}_b}$ includes the battery losses and is modeled using a piecewise quadratic approximation:

$$P_{\mathrm{s}}\mathrm{b}=\left\{\begin{array}{c}{r}_1^{+}·P_b+r_2^{+}·P_b^2,\mathrm{if}{P}_b\ge 0\hfill \\ r_1^{-}·P_b+r_2^{-}·P_b^2,\mathrm{if}{P}_b<0.\hfill \end{array}\right.$$

(14)

Here, the model parameters $r_1^{+}$, $r_2^{+}$ and $r_1^{-}$, $r_2^{-}$ are used to distinguish between battery charging and discharging losses. The value $P_b$ is defined as the summed up power, which is required by the electric motor, the generator, and the vehicle’s auxiliaries:

$$P_b=P_{s_{\mathrm{mot}}}+P_{s_{\mathrm{gen}}}+P_{\mathrm{aux}}.$$

(15)

The rest of the simplified model is divided into two different modeling structures, depending on the driving mode: In parallel mode, the powertrain components’ rotational speeds are predefined by the vehicle speed, which means that speed-dependent model parameters can be utilized. In series mode, additional speed-independent model parameters need to be identified for the generator and the H${ }_2$ICE.

Parallel mode:

The hydrogen consumption map is approximated using a speed-dependent quadratic Willans approach:

$${\dot{m}}_{\mathrm{f}}=f_0\left({\omega }_{\mathrm{e}}\right)+f_1\left({\omega }_{\mathrm{e}}\right)·P_{\mathrm{e}}\mathrm{GB}+f_2\left({\omega }_{\mathrm{e}}\right)·P_{\mathrm{e}}{\mathrm{GB}}^2.$$

(16)

Figure 4 shows the model fit for the entire hydrogen map.

The motor and generator losses are both approximated using a speed-dependent piece-wise affine model:

$$\begin{array}{cc}\hfill P_{\mathrm{l}}\mathrm{mot}& =\underset{j}{max}\{a_{0,j}\left({\omega }_{\mathrm{mot}}\right)+a_{1,j}\left({\omega }_{\mathrm{mot}}\right)·P_{\mathrm{mot}}\},\mathrm{for}\mathrm{j}\in \{1,...,5\},\hfill \\ \hfill P_{\mathrm{l}}\mathrm{gen}& =\underset{j}{max}\{b_{0,j}\left({\omega }_{\mathrm{gen}}\right)+b_{1,j}\left({\omega }_{\mathrm{gen}}\right)·P_{\mathrm{gen}}\},\mathrm{for}\mathrm{j}\in \{1,...,5\}.\hfill \end{array}$$

(17)

The model accuracy is found to be sufficient using 5 affine functions. Figure 5 shows the model fit for the motor and generator evaluated for one rotational speed of the corresponding component.

The NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions map is approximated using a speed-dependent piece-wise affine model:

$${\dot{m}}_{\mathrm{NO}}\mathrm{x}=\underset{j}{max}\{d_{0,j}\left({\omega }_e\right)+d_{1,j}\left({\omega }_e\right)·P_{\mathrm{e}}\mathrm{GB}\},\mathrm{for}\mathrm{j}\in \{1,...,10\}.$$

(18)

The model accuracy is found to be sufficient using 10 affine functions. Figure 6 shows the model fit for the entire NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions map.

Series mode:

Operating the vehicle in series mode detaches the engine-generator unit from the rest of the drivetrain, which allows to choose its operating point independent of the vehicle speed. For example, to produce a fixed generator power $P_{\mathrm{gen}}$, the selection of $({\omega }_{\mathrm{gen}},T_{\mathrm{gen}})$ can be optimized. From a control perspective, the series mode increases the flexibility of using the individual propulsion systems to satisfy EMS controller goals. From a system modeling perspective, however, the additional degree of freedom impedes the formulation of convex functions. This challenge is tackled by doing a pre-optimization of the engine-generator unit, covering all its operating conditions. The goal is to predefine for each obtainable generator power, which combinations of generator speed and generator torque are optimal candidates with respect to hydrogen consumption and/or NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions. The pre-optimization is based on the following extended cost function:

$$J={\dot{m}}_f+\psi ·{\dot{m}}_{\mathrm{NO}}\mathrm{x}.$$

(19)

Here, the trade-off factor $\psi$ is used as weight factor to reflect the H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off. Algorithm 1 is used to perform the pre-optimization. In the following, it is explained in detail:

• Step 1: The generator power and the trade-off factor $\psi$ are discretized.
• For each realizable $P_{\mathrm{gen}}$, all possible combinations (${\omega }_{\mathrm{e}}$, $T_{\mathrm{e}}$) that result in $P_{\mathrm{gen}}\left(j\right)$ are identified. Using the map-based model, the corresponding hydrogen consumption ${\dot{m}}_f^{\mathrm{map}}$ and the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions ${\dot{m}}_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{map}}$ are calculated (steps 3-5).
• Looping over all $\psi$, Equation 19 is used to formulate the extended cost for all identified pairs of (${\omega }_{\mathrm{e}}$,   $T_{\mathrm{e}}$) and the corresponding trade-off weight $\psi$ (step 7).
• Minimizing the extended cost function over all previously identified operating points (${\omega }_{\mathrm{e}}$,   $T_{\mathrm{e}}$) yields the optimal engine operating point (${\omega }_{\mathrm{e}}^{*}$,   $T_{\mathrm{e}}^{*}$) for the corresponding $\psi$ (step 8).
• Finally, for the generator power $P_{\mathrm{gen}}\left(j\right)$ and the trade-off parameter $\psi$, the following optimal values are stored for later use: Optimal engine power $P_{\mathrm{e}}^{*}$, optimal hydrogen consumption ${\dot{m}}_f^{*}$, and optimal NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions ${\dot{m}}_{\mathrm{NO}}{\mathrm{x}}^{*}$ (steps 9-11).

Algorithm 1 Pre-optimization for series mode

1: ${\overrightarrow{P}}_{\mathrm{s}}\mathrm{gen}\leftarrow \mathrm{Uniformly}\mathrm{distributed}\mathrm{vector}\in [-80kW,0kW]$   2: $\overrightarrow{\psi }\leftarrow \mathrm{Uniformly}\mathrm{distributed}\mathrm{vector}\in [0,150]$   3: ${\overrightarrow{\omega }}_{\mathrm{e}}\leftarrow \mathrm{Uniformly}\mathrm{distributed}\mathrm{vector}\in [{\omega }_e^{\min },{\omega }_e^{\max }]$  4: for$j\in \{1,2,...,|{\overrightarrow{P}}_{\mathrm{s}}\mathrm{gen}\left|\right\}$do  5:     ${\overrightarrow{T}}_{\mathrm{e}}\leftarrow$ for each entry in ${\overrightarrow{\omega }}_{\mathrm{e}}$, find $T_{\mathrm{e}}$ that results in ${\overrightarrow{P}}_{\mathrm{s}}\mathrm{gen}\left[j\right]$ according to eq:lossmap-mot-gen,eq:gengearbox,eq:gen-loss  6:     ${\overrightarrow{\dot{m}}}_f\leftarrow f({\overrightarrow{\omega }}_{\mathrm{e}},{\overrightarrow{T}}_{\mathrm{e}})$                        ▹ eq:fuelmap  7:     ${\overrightarrow{\dot{m}}}_{\mathrm{NO}}\mathrm{x}\leftarrow f({\overrightarrow{\omega }}_{\mathrm{e}},{\overrightarrow{T}}_{\mathrm{e}})$▹ eq:noxmap  8:     for $k\in \{1,2,...,|\overrightarrow{\psi }\left|\right\}$ do  9:         $\overrightarrow{J}\leftarrow {\overrightarrow{\dot{m}}}_f+\psi \left[k\right]·{\overrightarrow{\dot{m}}}_{\mathrm{NO}}\mathrm{x}$▹ eq:costpre-optimization  10:         $i^{*}\leftarrow \arg \min \overrightarrow{J}$▹ index of minimizer in $\overrightarrow{J}$  11:         ${\dot{m}}_f^{*}\leftarrow {\overrightarrow{\dot{m}}}_f\left[i^{*}\right]$  12:         ${\dot{m}}_{\mathrm{NO}}{\mathrm{x}}^{*}\leftarrow {\overrightarrow{\dot{m}}}_{\mathrm{NO}}\mathrm{x}\left[i^{*}\right]$  13:         $P_{\mathrm{e}}{\mathrm{GB}}^{*}\leftarrow {\overrightarrow{\omega }}_{\mathrm{e}}\left[i^{*}\right]·{\overrightarrow{T}}_{\mathrm{e}}\left[i^{*}\right]$  14:         $P_{\mathrm{s}}{\mathrm{gen}}^{*}\leftarrow {\overrightarrow{P}}_{\mathrm{s}}\mathrm{gen}\left[j\right]$  15:         ${\psi }^{*}\leftarrow \psi \left[k\right]$  16:         Store $(P_{\mathrm{e}}{\mathrm{GB}}^{*},{\psi }^{*},{\dot{m}}_f^{*},P_{\mathrm{s}}{\mathrm{gen}}^{*},{\dot{m}}_{\mathrm{NO}}{\mathrm{x}}^{*})$  17:     end for  18: end for  19: Return:${\dot{m}}_f=f(\psi ,P_{\mathrm{e}}\mathrm{GB})$▹ based on stored variables  20: Return:$P_{\mathrm{s}}\mathrm{gen}=f(\psi ,P_{\mathrm{e}}\mathrm{GB})$                        ▹ based on stored variables  21: Return:${\dot{m}}_{\mathrm{NO}}\mathrm{x}=f(\psi ,P_{\mathrm{e}}\mathrm{GB})$                        ▹ based on stored variables

Here, the lower bound ${\underline{P}}_{\mathrm{s}}\mathrm{gen}=-80$ kW is sufficiently large to cover the entire relevant range of engine-generator usage. The upper generator bound of 0 kW is set, as the generator is not used to provide additional mechanical propulsion power at the wheel. Setting the trade-off factor to its lower bound of $\psi =0$ results in an objective solely concerned about hydrogen consumption. Setting $\psi =150$, on the other hand, results in an objective that is solely concerned about NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions. Further increasing $\psi$ does not change the result of step 8. Thus, selecting $\psi \in [0,150]$ represents the complete achievable H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off.

Based on the results of the pre-optimization, the following model parameters are fitted. Notably, compared to the modeling approach for the parallel mode, these mappings have no explicit dependency on rotational speeds.

Figure 7 shows the approximation for the generator source power and the hydrogen consumption, based on above pre-optimized values. The generator source power is approximated using a hyperplane:

$$P_{\mathrm{s}}\mathrm{gen}={\beta }_0+{\beta }_1·P_{\mathrm{e}}\mathrm{GB}+{\beta }_2·\psi .$$

(20)

The hydrogen consumption is approximated using multiple hyperplanes:

$${\dot{m}}_{\mathrm{f}}=\underset{j}{max}\{{\varphi }_{0,j}+{\varphi }_{1{,}_j}·P_{\mathrm{e}}+{\varphi }_{2,j}·\psi \},\mathrm{for}\mathrm{j}\in \{1,...,5\}.$$

(21)

The model accuracy is found to be sufficient using 5 hyperplanes.

Figure 8 depicts the approximation of the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions in series mode, consisting of multiple hyperplanes.

$${\dot{m}}_{\mathrm{NO}}\mathrm{x}=\underset{j}{max}\{{\delta }_{0,j}+{\delta }_{1,j}·P_e+{\delta }_{2,j}·\psi \},\mathrm{for}\mathrm{j}\in \{1,...,17\}.$$

(22)

The fitting process is done using a slightly adapted version of the convex adaptive partitioning algorithm [19].

The convex adaptive partitioning algorithm is free to decide the number of hyperplanes and converges at the number of 17.

2.3. Optimization Parameters

The simplified model provides an approximation of the drivetrain components. Next, to arrive at the formulation of an optimization problem, the set of of dynamic states and the set of control inputs is described.

The dynamic states are the accumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ and the SoC:

$$\mathbf{x}=\left[\begin{array}{c}{m}_{\mathrm{NO}}\mathrm{x}\\ \mathrm{SoC}\end{array}\right].$$

(23)

The differential equation for the accumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ is defined by Equations (18) and (22) and reads as:

$$\frac{\mathrm{d}{m}_{\mathrm{NO}}\mathrm{x}}{\mathrm{d}t}=\left\{\begin{array}{cc}{\delta }_{1,j}·P_e+{\delta }_{2,j}·\psi +{\delta }_{0,j},\hfill & \mathrm{if}M=1\hfill \\ d_{1,j}·P_e+d_{0,j},\hfill & \mathrm{if}M=2.\hfill \end{array}\right.$$

(24)

The differential equation for the SoC is defined as follows:

$$\frac{\mathrm{d}\mathrm{SoC}}{\mathrm{d}t}=-\frac{P_{s_b}}{Q_0},$$

(25)

where $Q_0$ denotes the battery capacity and $P_{s_b}$ denotes the battery source power.

The control inputs for the simplified powertrain model are given in the vector:

$$\mathbf{u}=\left[\begin{array}{c}{P}_{\mathrm{e}}\mathrm{GB}\\ P_{\mathrm{mot}}\mathrm{GB}\\ P_{\mathrm{gen}}\mathrm{GB}\\ \psi \end{array}\right].$$

(26)

A validation of the simplified model description is presented in Figure 9. Both models are simulated using the input trajectories depicted by the black lines. The resulting trajectories for SoC, NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ and the resulting hydrogen consumption are depicted using the red lines and the green lines.

An offset results when the input trajectories are evaluated on the simplified model. The offset of 3 % in the hydrogen consumption and the offset of 2 % in the final SoC are limited. However, the offset of 15.8 % in the cumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions is larger. Nevertheless, the accuracy of the simplified model was found to be good enough for the use as a control-oriented model.

3. Control-Oriented Optimization Problem

In this section, the simplified model description is leveraged to formulate a convex optimization problem. However, in order to omit the formulation of a mixed-integer program (see Equation (Section 2.1)), a time-resolved trajectory of the integer driving mode M has to be found in advance.

3.1. Driving Mode Estimation

The estimation of the optimal driving mode M is based on the work published in [20]. Here, the authors propose a map, which is defined in the ($P_{\mathrm{req}}-v$) plane to estimate the optimal driving mode. If the map is chosen carefully, a close-to-optimal estimation of the driving mode is possible based on the vehicle velocity and the power request. An exemplary map used in this work is depicted in Figure 10.

It is obtained by solving Equation (Section 2.1) using the DP algorithm for a specific driving mission. The parameters $P_{\mathrm{crit}},v_{\mathrm{crit}}$, and $P_{\mathrm{par}}$ describe regions to estimate the optimal driving mode. Assuming predictive knowledge of the vehicle’s upcoming time-resolved velocity and power request, also a time-resolved estimation of the future optimal driving mode can be made.

3.2. Convex Optimization Problem

Based on an estimation of the time-resolved trajectory for the integer driving mode M, the convex optimization problem, which is from here on referred to as the $\mathcal{COP}$, is explained step-by-step in the following.

Cost function and dynamics:

The cost function and dynamic state equations of the $\mathcal{COP}$ read as follows:

$$\begin{array}{cc}\hfill \underset{\mathbf{u}}{\min }{\int }_{t_0}^{t_f}& ({\dot{m}}_f+{\kappa }_{\mathrm{P}}\mathrm{e}·{ϵ}_{\mathrm{P}}\mathrm{e})dt+{\kappa }_{\mathrm{NO}}\mathrm{x}·{ϵ}_{\mathrm{NO}}\mathrm{x}\hfill \end{array}$$

(27a)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.\frac{\mathrm{d}\mathrm{SoC}}{\mathrm{d}t}& =\mathrm{equation}\left(25\right)\hfill \end{array}$$

(27b)

$$\begin{array}{cc}\hfill \mathrm{SoC}\left(t_0\right)& ={\mathrm{SoC}}_0\hfill \end{array}$$

(27c)

$$\begin{array}{cc}\hfill \mathrm{SoC}\left(t_f\right)& \ge {\overline{\mathrm{SoC}}}^{\mathrm{RTG}}\hfill \end{array}$$

(27d)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{m}_{\mathrm{NO}}\mathrm{x}}{\mathrm{d}t}& =\mathrm{equation}\left(24\right)\hfill \end{array}$$

(27e)

$$\begin{array}{cc}\hfill m_{\mathrm{NO}}\mathrm{x}\left(t_0\right)& =m_{\mathrm{NO}}\mathrm{x},0\hfill \end{array}$$

(27f)

$$\begin{array}{cc}\hfill m_{\mathrm{NO}}\mathrm{x}\left(t_f\right)& \le {\overline{m}}_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{RTG}}+{ϵ}_{\mathrm{NO}}\mathrm{x}\hfill \end{array}$$

(27g)

$$\begin{array}{cc}\hfill P_{\mathrm{e}}\mathrm{GB}& \le P_{\mathrm{e}}{\mathrm{GB}}^{\max }+{ϵ}_{\mathrm{P}}\mathrm{e}\hfill \end{array}$$

(27h)

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{P}}\mathrm{e}& \ge 0\hfill \end{array}$$

(27i)

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{NO}}\mathrm{x}& \ge 0\hfill \end{array}$$

(27j)

The cost function (27a) consists of three parts. First, the integral of the hydrogen consumption ${\dot{m}}_f$, just as in the original formulation in (8a). Second, the integral of the quantity ${\kappa }_{\mathrm{P}}\mathrm{e}·{ϵ}_{\mathrm{P}}\mathrm{e}$, which, together with equations (27h) and (27i), represents a soft constraint for the engine power. It is used to ensure that the optimization problem remains feasible, even if an unrealistically high power request is predicted. A detailed explanation is provided in [20]. Third, the scalar term ${\kappa }_{\mathrm{NO}}\mathrm{x}·{ϵ}_{\mathrm{NO}}\mathrm{x}$, which, together with equations (27f) and (27j), represents a soft constraint on the accumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$. It is used to ensure that disturbances leading to increased NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ cannot lead to an infeasible optimization problem. It is important to ensure that constraints are only softened if no feasible solution with strict constraint compliance exists. Thus, the weights ${\kappa }_{{\mathrm{P}}_{\mathrm{e}}}$ and ${\kappa }_{\mathrm{NO}}\mathrm{x}$ have to be at least as high as the corresponding optimal Lagrange multipliers for the original problem [21]. Since these Lagrange multipliers are not known a priori, the weighting is chosen very large, i.e., 10${ }^5$ here.

The definition of the dynamic states requires three ingredients. First, both dynamic states are initialized with the values ${\mathrm{SoC}}_0$, $m_{\mathrm{NO}}\mathrm{x},0$ (Equations (27c) and (27f)). Second, the evolution of the dynamic states is defined by Equations (27b) and (27e). Third, the terminal conditions for both dynamic states are given by Equations (27d) and (27g).

Power split:

The goal of the $\mathcal{COP}$ is to satisfy the time-resolved power request trajectory optimally using the optimal control inputs $\mathbf{u}$. In the following, the power balance is described depending on the driving mode.

$$\begin{array}{cc}\hfill \underline{M=1:}& \left\{\begin{array}{cc}{P}_{\mathrm{req}}\mathrm{GB}\hfill & =P_{\mathrm{mot}}\mathrm{GB}\hfill \\ P_{\mathrm{s}}\mathrm{gen}\hfill & ={\beta }_1·P_{\mathrm{e}}\mathrm{GB}+{\beta }_2·\psi +{\beta }_0\hfill \end{array}\right.\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \underline{M=2:}& \left\{\begin{array}{cc}{P}_{\mathrm{req}}\mathrm{GB}\hfill & =P_{\mathrm{gen}}\mathrm{GB}+P_{\mathrm{mot}}\mathrm{GB}+P_{\mathrm{e}}\mathrm{GB}\hfill \end{array}\right.\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \underline{M=3:}& \left\{\begin{array}{cc}{P}_{\mathrm{req}}\mathrm{GB}\hfill & =P_{\mathrm{mot}}\mathrm{GB}\hfill \\ P_{\mathrm{e}}\mathrm{GB}\hfill & =0\hfill \\ P_{\mathrm{gen}}\mathrm{GB}\hfill & =0\hfill \end{array}\right.\hfill \end{array}$$

(30)

Based on the optimal power distribution, the hydrogen consumption, NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions, and electric losses result, which are formulated using constraint relaxations.

Constraint relaxations:

In the following, the equality constraints (16), (17), (18), (21), and (22) are relaxed, which is necessary to obtain a convex feasible set for the optimization parameters:

$$\begin{array}{cc}\hfill {\dot{m}}_f& \ge \left\{\begin{array}{c}{\alpha }_i{P}_{\mathrm{e}}\mathrm{GB}+{\beta }_{i}s+{\gamma }_i,i\in \{1,..,5\},\mathrm{if}M=1\hfill \\ f_0+f_1{P}_{\mathrm{e}}\mathrm{GB}+f_2{P}_{\mathrm{e}}{\mathrm{GB}}^2,\mathrm{if}M=2\hfill \\ 0,\mathrm{if}M=3,\hfill \end{array}\right.\hfill \end{array}$$

(31a)

$$\begin{array}{cc}\hfill {\dot{m}}_{\mathrm{NO}}\mathrm{x}& \ge \left\{\begin{array}{c}{\alpha }_i{P}_{\mathrm{e}}\mathrm{GB}+{\beta }_{i}s+{\gamma }_i,i\in \{1,..,17\},\mathrm{if}M=1\hfill \\ {\alpha }_i{P}_{\mathrm{e}}\mathrm{GB}+{\gamma }_i,i\in \{1,..,10\},\mathrm{if}M=2\hfill \\ 0,\mathrm{if}M=3,\hfill \end{array}\right.\hfill \end{array}$$

(31b)

$$\begin{array}{cc}\hfill P_{\mathrm{l}}\mathrm{gen}& \ge \left\{\begin{array}{c}0,\mathrm{if}M\in \{1,3\}\hfill \\ {\alpha }_i{P}_{\mathrm{gen}}+{\gamma }_i,i\in \{1,..,5\},\mathrm{if}M=2,\hfill \end{array}\right.\hfill \end{array}$$

(31c)

$$\begin{array}{cc}\hfill P_{\mathrm{l}}\mathrm{mot}& \ge {\alpha }_i{P}_{\mathrm{mot}}+{\gamma }_i,i\in \{1,..,5\}.\hfill \end{array}$$

(31d)

An in-depth explanation of constraint relaxations is presented in [22]. However, intuitively, each parameter subject to the relaxed Equations (31a), (31b), (31c) and (31d) should be as small as possible to accumulate as little hydrogen consumption, NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions, or electrical losses as possible. Therefore, the equations hold with equality if used within an optimization problem.

Battery:

The following set of constraints reformulates the convex battery model introduced by Equations (14) and (15):

$$\begin{array}{cc}\hfill P_b^{\mathrm{tot}}& =P_{\mathrm{s}}\mathrm{mot}+P_{\mathrm{s}}\mathrm{gen}+P_{\mathrm{aux}}\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill P_b^{\mathrm{tot}}& =P_b^{+}+P_b^{-}\hfill \end{array}$$

(32b)

$$\begin{array}{cc}\hfill P_{\mathrm{s}}\mathrm{b},1& \ge r_1^{+}{P}_b^{+}+r_2^{+}{\left(P_b^{+}\right)}^2\hfill \end{array}$$

(32c)

$$\begin{array}{cc}\hfill P_{\mathrm{s}}\mathrm{b},2& \ge r_1^{-}{P}_b^{-}+r_2^{-}{\left(P_b^{-}\right)}^2\hfill \end{array}$$

(32d)

$$\begin{array}{cc}\hfill P_{\mathrm{s}}\mathrm{b}& \ge P_{\mathrm{s}}\mathrm{b},1+P_{\mathrm{s}}\mathrm{b},2\hfill \end{array}$$

(32e)

$$\begin{array}{cc}\hfill P_b^{+}& \ge 0\hfill \end{array}$$

(32f)

$$\begin{array}{cc}\hfill P_b^{-}& \le 0.\hfill \end{array}$$

(32g)

The auxiliary variables $P_b^{\mathrm{tot}}$, $P_b^{+}$, $P_b^{-}$, $P_{\mathrm{s}}\mathrm{b},1$, and $P_{\mathrm{s}}\mathrm{b},2$ are introduced to formulate the piece-wise quadratic battery model using convex set representations. Equations (32a), (32b) and (32c) are relaxed equality constraints.

Input and state domains:

The feasible control input domain and state domain are given in the following:

$$\begin{array}{c}\hfill \mathbf{u}\in \mathcal{U}\end{array}$$

(33a)

$$\begin{array}{c}\hfill \mathbf{x}\in \mathcal{X},\end{array}$$

(33b)

where $\mathcal{U}$ includes the physical limitations of the propulsion systems. The feasible set $\mathcal{X}$ consists of constant upper and lower battery bounds, i.e., SoC ∈ [0.5, 0.9]. It is important to note that the $\mathcal{COP}$ is only a convex optimization problem, if the driving mode M is known in advance. The $\mathcal{COP}$ is parsed using the software package YALMIP [23] and solved using the solver MOSEK [24]. Euler forward integration is used and the time discretization is set to 1 s.

4. Controller Structure

In this section, the proposed control architecture to calculate the vehicle control inputs is presented. The overall controller architecture is presented in Figure 11 and is based on the work presented in [20]. The controller consists of three controller levels: the lower-level controller, the so-called MPC, and the reference trajectory generator, which is from here on referred to as the RTG. Each controller level receives a power request signal, i.e., $P_{\mathrm{req}}$, ${\widehat{P}}_{\mathrm{req}}$ and ${\overline{P}}_{\mathrm{req}}$, respectively. The corresponding power request signals are either determined by the driver or calculated using the elevation profile $\Gamma$ of the driving mission and a corresponding predicted velocity signal ($\overline{v}$ or $\widehat{v}$). As predictive data is likely to be imperfect and of lower accuracy for longterm predictions than for short-term predictions, the accents $\overline{\left\{\right\}}$ and $\widehat{\left\{\right\}}$ are used to denote signals of differing prediction accuracy. The corresponding power request signals are calculated in the block labeled PR using a set of equations of the longitudinal vehicle dynamics [25].

In the following, a bottom-up approach is used to explain the overall controller structure in detail. The focus is on the communication between the individual controller levels.

4.1. Lower-Level Controller

The lower-level controller calculates the powertrain control inputs $T_{\mathrm{e}}$, $T_{\mathrm{mot}}$, $T_{\mathrm{gen}}$, ${\omega }_{\mathrm{gen}}$ and M in a real-time-feasible fashion. It is a non-predictive controller and obtains the current value $P_{\mathrm{req}}$. The control inputs are passed on to the vehicle plant model, which can be interpreted as a digital twin of the mixed H${ }_2$-HEV. The battery SoC and the accumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions can be accessed from the vehicle plant model in a causal fashion. They are denoted with SoC${ }^{\mathrm{meas}}$ and $m_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{meas}}$. The lower-level controller is based on the concept of Pontryagin’s Minimum Principle (PMP) and requires knowledge of the so-called costates, which can be interpreted as equivalent cost terms in a local Hamiltonian minimization problem [26]. The theory of PMP states that for an optimization problem, there exists one costate variable for each dynamic state variable. Therefore, two costates are required for the EMS at hand: ${\lambda }_{\mathrm{SoC}}$ and ${\lambda }_{N{O}_{\mathrm{x}}}$. Given an estimation of the costates, the calculation of controller inputs based on the application of PMP has already been shown multiple times (see e.g. [27]) and is therefore not further discussed in this text.

4.2. MPC

The MPC is used to estimate the value of the optimal costates, which are required by the lower-level controller. It solves the $\mathcal{COP}$ based on medium-term time-resolved predictive information for the power request of the upcoming driving mission. In Equation Section 3.2, the parameter $t_f$ is replaced by $t_0+H_p$, where $H_p$ denotes the length of the predictive information and is called the prediction horizon. The MPC also consists of a map, denoted in Figure 11 by the block $\widehat{\mathrm{Map}}$, which is defined in the (${\widehat{P}}_{\mathrm{req}}-\widehat{v}$) plane, and the $\mathcal{COP}$. This map is not static, but is adapted over time using a learning algorithm, which is explained in detail in [20]. The controller is updated every $T_s^{\mathrm{MPC}}$ seconds, where each update includes the following steps: The map is used to estimate the time-resolved optimal driving mode $\widehat{M}$ based on the predicted power request. The dynamic states of the $\mathcal{COP}$ (Equations (27c) and (27f)) are initialized using the signals SoC${ }^{\mathrm{meas}}$ and $m_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{meas}}$. The terminal conditions for the dynamic states (Equations (27d) and (27g)) are set. This is is done using interpolation of the spatially defined vectors SoC${ }^{\mathrm{RTG}}\left(s\right)$ and $m_{{\mathrm{NO}}_{\mathrm{x}}}^{\mathrm{RTG}}\left(s\right)$ and the predicted vehicle position $s(t_0+H_p)$ at the end of the prediction horizon. The two costates ${\lambda }_{\mathrm{SoC}}$ and ${\lambda }_{N{O}_{\mathrm{x}}}$ are obtained by solving the $\mathcal{COP}$ [20].

4.3. Reference Trajectory Generator

The reference trajectory generator is used to calculate suitable reference state trajectories for the MPC. The use of reference trajectories for online-capable EMS algorithms is a well-documented technique (see e.g. [28]). Typically, these approaches only encompass a reference for the battery. In order to ensure the possibility to aim at a user-defined long-term NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ target ${\overline{m}}_{\mathrm{NO}}{\mathrm{x}}^t$, also a NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trajectory needs to be included in this work. The RTG is a predictive algorithm, receiving long-term time-resolved predictive information for the power request of the entire driving mission. It consists of a map, denoted in Figure 11 by the block $\overline{\mathrm{Map}}$, which is defined in the (${\overline{P}}_{\mathrm{req}}-\overline{v}$) plane, and the $\mathcal{COP}$. The map is used to estimate the time-resolved optimal driving mode $\overline{M}$ based on the predicted power request. The $\mathcal{COP}$ in the RTG is set up only once and solved prior to departure. The initialization of the dynamic states is done using the current battery SoC and $m_{\mathrm{NO}}\mathrm{x},0=0$ (Equations (27c) and (27f)). The terminal condition for the battery is set to be charge sustaining (Equation (27d)) and the terminal condition for the cumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions is set to the user-defined long-term target ${\overline{m}}_{\mathrm{NO}}{\mathrm{x}}^t$.

Algorithm 2 explains the procedure of obtaining reference trajectories. First, the algorithm is initialized with a user-defined NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ target and a critical power request $P_{\mathrm{crit}}^{\mathrm{RTG}}$ (see Figure 10). After calculating the predicted power request for the entire driving mission, the $\mathcal{COP}$ is solved iteratively using the estimated optimal driving mode, resulting from the current $P_{\mathrm{crit}}^{\mathrm{RTG}}$. If the optimal control problem can be solved while respecting the user-defined NO${ }_{\mathrm{x}}^{\mathrm{eo}}$, i.e., ${ϵ}_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{RTG}}=0$, the resulting optimal state trajectories are used as reference trajectories for the MPC. In most cases only one iteration suffices. However, Equation (27g) is checked for slacking first. If this constraint requires slacking, then the resulting NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trajectory doesn’t satisfy the user-defined accumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ target. If this trajectory would be used as a reference for the MPC, the controller would track a NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ target, which is different from the user-defined target by design. If slacking is needed, $P_{\mathrm{crit}}^{\mathrm{RTG}}$ is decreased for the next iteration. This allows to use the series mode more often, allowing more ultra-low NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ operating points of the H${ }_2$ICE. During simulation, it was seen that decreasing $P_{\mathrm{crit}}^{\mathrm{RTG}}$ below $-40$ kW does not lead to noticeable changes in the solution of the $\mathcal{COP}$. It was therefore used as a termination criterion. If the termination criterion is met but slacking is still required, the resulting reference state trajectories are used. However, the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ state trajectory is linearly scaled to ${\overline{m}}_{\mathrm{NO}}{\mathrm{x}}^t$. Although such a post-processing results in an infeasible NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ state trajectory for the RTG, it can be feasible for the MPC, which calculates optimal solutions based on different predictive data.

Algorithm 2RTG Iterations

1.55 1: $P_{\mathrm{crit}}^{\mathrm{RTG}}=12000\mathrm{W}$ 2: ${ϵ}_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{RTG}}=1$ 3: ${\overline{P}}_{\mathrm{req}}^{\mathrm{RTG}}\leftarrow \overline{\Gamma },\overline{v}$ 4: while${ϵ}_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{RTG}}\ne 0|P_{\mathrm{crit}}^{\mathrm{RTG}}>-40\mathrm{kW}$do 5:     $M\leftarrow \mathrm{map}\{{\overline{P}}_{\mathrm{req}},\overline{v}\}$ 6:     {${\mathrm{SoC}}^{\mathrm{RTG}},m_{{\mathrm{NO}}_{\mathrm{x}}}^{\mathrm{RTG}},{ϵ}_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{RTG}}\}\leftarrow \mathrm{solve}\mathcal{COP}({\overline{P}}_{\mathrm{req}},M)$ 7:     if ${ϵ}_{\mathrm{NO}}{\mathrm{x}}^{\mathrm{RTG}}\ne 0$ then 8:          $P_{\mathrm{crit}}^{\mathrm{RTG}}$ = $P_{\mathrm{crit}}^{\mathrm{RTG}}$ - 500 W 9:     end if 10: end while

5. Case Study

This section provides a case study on the performance of the MPC. To quantify the performance of controllers, the following relative performance criteria are used:

$$\begin{array}{cc}\hfill e_{{\mathrm{NO}}_{\mathrm{x}}}& =1-\frac{m_{\mathrm{NO}}\mathrm{x}\left(t_f\right)}{{\overline{m}}_{\mathrm{NO}}{\mathrm{x}}^t},\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill e_{\mathrm{SoC}}& =1-\frac{\mathrm{SoC}\left(t_f\right)}{\mathrm{SoC}\left(t_0\right)},\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill e_{{\mathrm{H}}_2}& =\left(m_{\mathrm{H}}2\left(t_f\right)+\frac{Q_0·(\mathrm{SoC}\left(t_f\right)-\mathrm{SoC}\left(t_0\right))}{{\eta }_{\mathrm{H}}2{\mathrm{ICE}}^{\mathrm{const}}·{\eta }_{\mathrm{gen}}^{\mathrm{const}}·H_l}\right)/m_{\mathrm{H}}{2}^{\mathrm{DP}}-1.\hfill \end{array}$$

(36)

The first two, i.e., the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-error $e_{{\mathrm{NO}}_{\mathrm{x}}}$ and the SoC-error $e_{\mathrm{SoC}}$ quantify the violation of the user-defined NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target and the charge-sustainibility constraint at the end of the driving mission, respectively. The corrected hydrogen error $e_{{\mathrm{H}}_2}$ denotes the additional fuel spent with respect to the DP-optimal result. To account for the additional energy stored in the battery at the end of the driving mission, the fuel consumption is corrected using the lower heating value of hydrogen $H_l$ and a constant engine efficiency ${\eta }_{\mathrm{H}}2{\mathrm{ICE}}^{\mathrm{const}}=35$ % and generator efficiency ${\eta }_{\mathrm{gen}}^{\mathrm{const}}=90$ %.

For all results, the MPC’s prediction horizon is set to 450 s. The resulting $\mathcal{COP}$ takes below 0.5 s to solve using a standard laptop computer. The MPC’s update time is $T_s=2$ s. The map-based model is used as the vehicle plant.

Overall, after introducing the driving missions, the subsequent case study is built up in three steps, which together form an inductive approach to analyze the proposed controller’s performance. First, the controller’s ability to adhere to a single NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target is investigated on one specific driving mission. Second, on the same driving mission, the full range of achievable NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target is investigated. Third, the controller’s performance is analyzed for a broad range of different realistic driving missions.

5.1. Driving Missions

Figure 12 depicts the four driving missions that are used in the remainder of this text to test the performance of the suggested EMS algorithm. The first driving mission is extracted from the open-source software SUMO [29] and is referred to as urban cycle. The second driving mission is based on recorded sensor data of a real vehicle and features a mix of urban driving, rural driving, and highway driving and is referred to as real driving cycle. The third and fourth driving missions are generated from combinations of measured and simulated velocity profiles. They are used to represent extreme driving conditions and are referred to as mountain driving cycle and highway driving cycle.

The driving missions are visualized in Figure 12 over the travel distance s. The visualization includes three velocity trajectories for each mission in the upper plot and the elevation profile in the lower plot. The blue lines depict the actual vehicle speed v, which is interpreted as the driver’s request and is not known in advance by the MPC and the RTG. The red lines depict a predicted velocity $\widehat{v}$, which is available to the MPC. The dashed black lines depict a crude velocity estimation $\overline{v}$, which could represent, for instance, the legal speed limits.

5.2. Single NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-Target Adherence

First of all, the MPC’s ability to comply with a user-defined NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target is analyzed. Here, two sources of disturbance are investigated. First, the model mismatch between the the map-based model and the simplified model is always present. Second, the MPC relies on inaccurate predictive data. To distinguish between the performance loss resulting from either of the two disturbances, a Prescient MPC is introduced, which has access to the entire trajectory of the actual vehicle speed for the entire driving mission. Together with the full elevation profile, the Prescient MPC therefore has perfect knowledge of the entire upcoming driving mission, which means that the only disturbance to this controller remains the model mismatch.

Figure 13 shows a comparison between the results obtained by the DP algorithm, the Prescient MPC and the MPC on the real driving cycle. For all three methods, the same NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ target is chosen. The upper plot shows the accumulated NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions, normalized by the amount achieved by the DP algorithm. Both MPC controllers are able to satisfy the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ target. When looking at the three trajectories, it can be observed that all three methods have a similar way of distributing the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions across the driving mission: In time window A, the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions are comparably high. In time window B, all algorithms choose to emit less NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions than in the first part of the driving mission. Finally, in time window C, all three algorithms show increased NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions. The Prescient MPC yields a NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-error of $0.01$% and the MPC yields a NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-error of $0.4$%. It is therefore concluded that realistic velocity mispredictions have a very low impact on satisfying the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target on this driving mission and for this specific NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target.

The middle plot depicts the SoC trajectories. In time window D, the MPC chooses to divert from the SOC trajectory chosen by the DP algorithm and the Prescient MPC. This is a result of the mispredictions and its reduced prediction horizon. Even though in time window E, both model predictive controllers choose a less aggressive use of the battery than the DP algorithm, both algorithms are able to operate the vehicle in charge sustaining mode.

The bottom plot depicts the corrected hydrogen error of the three algorithms. Time windows D and E highlight the reflection of the algorithms’ different choices of battery usage in their corresponding hydrogen allocation strategy: If an algorithm chooses to charge the battery more than the DP algorithm, its hydrogen consumption is increased as a result of transferring energy from the hydrogen storage into the battery. Ultimately, both model predictive controllers result in a close-to-optimal hydrogen consumption: The Prescient MPC reaches $e_{{\mathrm{H}}_2}$ = 1.3% and the MPC reaches $e_{{\mathrm{H}}_2}$ = 2.3%. The slightly lower hydrogen consumption of the Prescient MPC is due to its access to perfect predictive data, which is leveraged to calculate a better estimation of the costates.

In this subsection, it was shown that the MPC is able to comply with one user-defined NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ target and satisfy charge-sustaining vehicle operation on the real driving cycle, while reaching close-to-optimal hydrogen consumption. Moreover, the additional loss of optimality by using realistic predictive data instead of perfect predictive data is found to be marginal. The next step is to show that the proposed MPC is able to satisfy all physically possible choices of NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets, without the need to adjust or re-tune the controller.

5.3. NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-Target Expansion

To investigate the ability of the MPC to reach different NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets, a simulation study is performed on the real driving cycle for varying NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets. Figure 14 includes a total of three plots, which are explained in detail in the following paragraphs.

The top plot shows the entire physically achievable H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ Pareto front in gray. This is calculated offline, using the DP algorithm with perfect predictive data and the full map-based model. The gray arrow depicts the maximum amount of NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ reduction over the H${ }_2$ optimal solution calculated by the DP algorithm for this driving mission. It is defined as 100%. The blue curve depicts the performance achieved by the online-capable MPC, which is subject to model mismatch and mispredictions, as mentioned above. This curve will from now on be referred to as the MPC’s performance front. Its results are obtained by running twelve simulations on the real driving cycle with twelve different NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets spanning the entire Pareto front. The MPC is able to reproduce a large part of the offline achievable Pareto front, including extremely low NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets, recovering more than 99% of the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ reduction potential over the H${ }_2$-optimal solution. The blue dashed line is a horizontal extension of the point that is marked by a blue X. It serves as a visual aid and depicts a part of the Pareto front that cannot be reached by the MPC, even if the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target is set to a higher value. This is a result of the model mismatch which is introduced by using the simplified emission models, which underestimate the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$. However, as indicated by the DP-optimal Pareto front, there is an extremely flat trade-off between H${ }_2$ and NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ in this region. Consequently, the benefits in H${ }_2$ efficiency are offset by a significant surplus of NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions.

The middle plot depicts the MPC’s NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-error and SoC-error for all twelve iterations. The blue line is very steep approaching the lowest NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets of the theoretically reachable Pareto front. This shows that more NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ is emitted by the MPC and thus these targets cannot be achieved by the online-capable controller. However, already for the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target which is only 1.4% of the maximum NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions on this driving mission, the MPC’s NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-error is 5% and well contained below 1% for higher NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets. The region of very high NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets is here visualized by the steeply declining error towards negative values, highlighting that the MPC cannot achieve extremely high NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets. This part corresponds to the region represented by the dashed blue line in the top plot.

The bottom plot shows the MPC’s corrected hydrogen consumption error for the different NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ targets. It shows that for the entire achievable Pareto front, the MPC is able to operate the vehicle with close-to-optimal hydrogen consumption.

In this subsection, it was shown that the MPC is able to achieve close-to-optimal hydrogen consumption consistently for a wide variety of NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets. Furthermore, the MPC can achieve almost the entire possible NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ reduction potential that was predicted by DP. This shows that the MPC is able to reproduce almost the entire H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ Pareto front on the real driving cycle. To conclude the case study, it remains to be shown whether these results can be generalized to other driving missions.

5.4. Driving Mission Generalization

In the following, the MPC’s ability to reproduce the H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ Pareto front is presented for a broad range of driving missions. Figure 15 depicts the comparison between the offline calculated DP-optimal Pareto fronts and the MPC’s online calculated performance fronts obtained for the urban driving mission, the mountain driving mission, and the highway driving mission. For the sake of comparison, the results obtained for the real driving mission are also included. The gray lines depict the entire, physically achievable H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ Pareto front, calculated by the DP algorithm with perfect predictive data and the full map-based model. The blue curves represent the performance fronts achieved by the online-capable MPC, which is subject to model mismatch as well as mispredictions. Its results are obtained by running multiple simulations on the corresponding drive mission, but with different NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets.

On all four driving missions, the MPC is able to reach extremely low NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions, recovering more than 99% of the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ reduction potential over the H${ }_2$-optimal solution. It cannot capture the theoretically possible EMS calibrations to the far right hand side of the Pareto front, which result in the lowest possible H${ }_2$ consumption. However, as indicated by the DP-optimal Pareto fronts, there is an extremely flat trade-off between H${ }_2$ and NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ for all investigated driving missions in this region. Consequently, the benefits in H${ }_2$ efficiency are offset by a significant surplus of NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions.

Table 2 summarizes the MPC’s performance fronts for all investigated driving missions. The two columns $e_{\mathrm{H}}{2}^{\mathrm{mean}}$ and $e_{\mathrm{H}}{2}^{\max }$ represent the mean and the maximum values, respectively, of all the $e_{\mathrm{H}}2$ resulting from the simulation on the corresponding driving mission and the individual NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-targets. Two things can be noted here: First, in all four driving missions, the maximum hydrogen error is close to the mean hydrogen error. This suggests that variance in the controller performance can be characterized as low. Second, across all four driving missions, the mean hydrogen error is similar. This suggests that the controller performance is similarly consistent across different driving missions. Overall, these results show that the proposed MPC is able to provide predictable, close-to-optimal hydrogen consumption for a wide variety of driving missions and for almost the entire physically reachable range of H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off.

6. Conclusion

To the authors’ best knowledge, this paper provides the first results of an online-capable control algorithm that for the EMS of a H${ }_2$-HEV, explicitly including the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions. The algorithm is based on a multi-level control approach and features a supervisory MPC algorithm, which solves a convex optimization problem. The proposed controller is evaluated in a case study, comprising three levels of complexity: First, the controller is shown to be able to satisfy a user-defined NO${ }_{\mathrm{x}}^{\mathrm{eo}}$-target, while operating the vehicle in charge-sustaining mode and achieving close-to-optimal hydrogen consumption. Second, the controller is shown to be able to extend its close-to-optimal hydrogen consumption to a wide range of different NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ targets. Third, it is shown that the proposed controller exhibits predictable performance for a wide variety of driving missions. In all tested scenarios, at most an increase of 4.66% in hydrogen consumption is observed, compared to the globally optimal solution. More importantly, for each driving mission, more than 99% of the possible NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ reduction potential can be reached by the online-capable MPC controller. In summary, this means that the proposed controller can almost fully exploit the theoretical potential of an H${ }_2$-HEV regarding the H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off.

It is important to emphasize that the obtained NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ stem from a map, which is derived from steady-state measurements. While these maps can predict the NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions of an H${ }_2$ICE well during steady-state operation, higher NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions may occur during transient operating point switches of the H${ }_2$ICE. Therefore, no statements regarding compliance with legislative limits should be made.

In future research, the development and inclusion of an additional dynamic NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ emissions model and its inclusion in the MPC could be looked into. Also, the possibility to adjust the H${ }_2$-NO${ }_{\mathrm{x}}^{\mathrm{eo}}$ trade-off online, could lead to changed requirements for the exhaust-gas aftertreatment system. This could create an avenue for exploring the optimization of EMS calibration in conjunction with the design and operation of the exhaust gas aftertreatment system.

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