Time-Optimal Low-Level Control and Gearshift Strategies for the Formula 1 Hybrid Electric Powertrain

Abstract

Today, Formula 1 race cars are equipped with complex hybrid electric powertrains that display significant cross-couplings between the internal combustion engine and the electrical energy recovery system. Given that a large number of these phenomena are strongly engine-speed dependent, not only the energy management but also the gearshift strategy significantly influence the achievable lap time for a given fuel and battery budget. Therefore, in this paper we propose a detailed low-level mathematical model of the Formula 1 powertrain suited for numerical optimization, and solve the time-optimal control problem in a computationally efficient way. First, we describe the powertrain dynamics by means of first principle modeling approaches and neural network techniques, with a strong focus on the low-level actuation of the internal combustion engine and its coupling with the energy recovery system. Next, we relax the integer decision variable related to the gearbox by applying outer convexification and solve the resulting optimization problem. Our results show that the energy consumption budgets not only influence the fuel mass flow and electric boosting operation, but also the gearshift strategy and the low-level engine operation, e.g., the intake manifold pressure evolution, the air-to-fuel ratio or the turbine waste-gate position.

1. Introduction

Formula 1 (F1) is the top category for single-seater circuit racing and features the fastest racing cars in the world. In the past, these were powered by V8, V10 or even V12 engines. Since 2014, however, the Fédération Internationale de l’Automobile (FIA) [1,2] has prescribed the parallel hybrid powertrain topology shown in Figure 1. In this so-called power unit, the main source of propulsive power is still the internal combustion engine, but in the form of a downsized and turbocharged V6 with a displacement volume of $1.6\mathrm{L}$. However, the engine is supplemented by an electric recovery system composed of two electric motor generator units, denoted by motor generator units (MGUs), and an electric energy storage in the form of a battery. The MGU-K (K for kinetic) is coupled to the engine’s crankshaft and is used to increase the propulsive power of the car or to recuperate energy during braking. To increase the powertrain’s efficiency, the FIA allows the installation of a second electric motor on the turbocharger shaft, referred to as MGU-H (H for heat). As a consequence, the turbine is oversized and in general produces more power than absorbed by the compressor, such that part of the energy contained in the hot exhaust gases can be recuperated by the MGU-H and fed to the battery or directly to the MGU-K. Additionally, the MGU-H can be used to speed up the compressor to avoid turbo lag. While the FIA restricts the MGU-K power to $\pm 120\mathrm{kW}$, no limits are imposed on the MGU-H. Finally, the power unit is connected to the wheels by an eight-speed sequential gearbox and a differential. In addition to the finite capacity of the battery of $4\mathrm{MJ}$, the FIA also limits the fuel consumption for an entire race (≈300 km) to $110\mathrm{kg}$ and constrains the maximum fuel mass flow as a function of the engine speed. Every F1 driver’s goal is to win the race, i.e., to complete the defined number of laps in the shortest amount of time. Therefore, controlling the powertrain in a time optimal way is of utmost importance. The complex interactions between the thermal and the electrical part, which are coupled not only through the drivetrain but also through the electrified turbocharger, together with the limited available energy imposed by the FIA rules, call for a rigorous model-based optimization of the powertrain operation. Since the behavior of several internal combustion engine components highly depends on the engine speed, any offline optimization must also take into account the gearshift strategy. Therefore, in this paper we present a framework that jointly computes the time-optimal control and gearshift strategy for a detailed low-level model of the current hybrid electric F1 powertrain. In the following paragraphs, we review the relevant literature, formulate the research statement and outline our contribution, and give an overview on the paper structure.

Given that the F1 power unit is hybrid electric, the methodology applied to optimize its energy management can be related to hybrid electric vehicle research: for road-going vehicles, such as passenger cars, the main objective is to minimize the fuel consumption, whilst for race cars the key performance indicator is the lap time. Therefore, we identify two streams of research related to our topic. The first one deals with the fuel-optimal control of hybrid electric vehicles, where in the majority of the cases the power request and the vehicle speed are known and given, i.e., the vehicle is assumed to follow a specific driving cycle. Non-causal optimal control strategies were investigated using dynamic programming [3,4,5], convex optimization [6,7,8] and Pontryagin’s minimum principle [9,10,11]. The sizing of the powertrain components [12], of the drivetrain [13,14] or the design of the powertrain topology [15] were also discussed in previous works. For feedback control, equivalent consumption minimization strategies have been applied for fuel consumption minimization [16,17] considering also pollutant emissions [18], the battery state of health [19] or turbocompounding [20,21]. Furthermore, the integer nature of gears and the engine on/off choice were tackled with iterative algorithms [7], Pontryagin’s minimum principle [22,23], dynamic programming [8,23], outer convexification [24,25,26,27,28] and shooting or bisection methods [29,30]. However, all these approaches mainly rely on quasi-static system models and ignore the dynamic behavior of the internal combustion engine, e.g., the intake manifold or the turbocharger dynamics. To properly capture the engine dynamics, mean value engine models are usually formulated using nonlinear differential equations [31,32,33,34] and optimized offline [35,36,37,38,39,40] or online in a model predictive control fashion [41,42,43,44]. Recently, the air path of a forced-intake engine was modeled in a quasi-static manner as a basis for minimum fuel control [45,46]. However, in these works several assumptions about the air-to-fuel ratio and the turbocharger operation were made to render the resulting model convex. In Ref. [47] instead, a piecewise affine low-level model of a high-performance hybrid electric power unit is presented. The proposed model includes the most relevant low-level engine dynamics, as well as the integer nature of the cylinder deactivation, and presents a methodology to fit the highly nonlinear turbocharger maps by means of neural networks techniques. While the resulting fuel-optimal mixed integer linear program has global optimality guarantees, it relies on measured engine speed and power request trajectories, i.e., the gearshift command and the vehicle speed are not optimized but are given. The second stream of research concerns the time-optimal control of race cars. In Ref. [48,49,50], the authors solved the general problem of evaluating the vehicle lateral and longitudinal control that yields to the minimum time maneuvering around a racetrack. Structural components that additionally influence the achievable lap time were included in the optimization framework: the suspensions [51] and the tyres [52]. The computation of the optimal driving trajectory was further extended with the inclusion of power flows [53], of the F1 specific power flows and regulations [54], or of the power flows and therewith associated thermal management of a Formula E race car [55]. Finally, the authors of [56] show that the time-optimal energy management marginally influences the race car’s velocity evolution in the corners, whilst it significantly impacts its speed on the straights. The most important assumption made in the previous works on F1 in our group, is to decouple the energy management from the driving path: we assume the driver to be responsible for following the time-optimal path on the racetrack, whilst we put the focus on the optimal usage of the available energy [57,58,59,60,61,62]. The time-optimal control strategies for race and qualifying scenarios were solved numerically in [57,60], while real-time feedforward and feedback control architectures were derived analytically and tested in [58,59]. A two-level model predictive control architecture was also proposed as an online controller in [62] and afterwards extended with adaptation algorithms in [61]. Finally, in [63] we combined convex optimization, dynamic programming and Pontryagin’s minimum principle in an iterative scheme to solve the energy management and gearshift control of a modern F1 powertrain. However, the majority of the models used for minimum time optimization are high-level models that mainly rely on quasi-static formulations, ensuring low computational times in the order of seconds to minutes. Thus, they can serve a rapid parametrization of the energy management but neglect important effects such as transient motor generator unit (MGU)-H operation during the gearshifts. Moreover, they do not consider the dynamics of the engine, e.g., the intake manifold and the turbocharger, they do not model the engine cylinder deactivation and do not include the air-to-fuel ratio and the spark advance influence on the achievable engine power.

To the best of the authors’ knowledge, there exists no framework so far that jointly computes the time-optimal control and gearshift strategy for a detailed low-level model of a F1 hybrid electric powertrain. The impact of both the gearshift strategy and the targets for fuel and battery usage on the engine power generation and on the low-level dynamics, e.g., on the air-to-fuel ratio or on the intake manifold pressure evolution, has not been studied yet for time-optimal control purposes. Therefore, in this paper we propose a detailed low-level mathematical formulation of the F1 powertrain suited for numerical optimization. Beside the specific focus on the engine speed dependency of each powertrain’s component, we propose a continuous nonlinear model to describe the engine cylinder deactivation and reformulate non-smooth inequality constraints such that they can be accommodated by off-the-shelf nonlinear solvers. We relax the resulting mixed integer nonlinear model using outer convexification techniques, thereby reducing computational times for solving the time-optimal control problem by two orders of magnitude. The presented framework can be used as a tool to benchmark the low-level actuation of the internal combustion engine components in the context of the time-optimal energy management. Specifically, the suboptimality of the control strategy adopted during a race week-end can be assessed, reference trajectories for online control can be generated and control heuristics can be inferred. Moreover, the framework allows to perform parametric studies aimed at the design of new powertrain components. For the purpose of this work, only a part of the FIA sporting and technical regulations [1,2] are considered and modeled. As an example, specific limits on the energy transfer between the battery and the electric motor connected to the crankshaft are not implemented. This is due to confidentiality reasons and also facilitates the results analysis given the complexity of the system.

The paper is structured as follows: in Section 2 we present a mixed integer nonlinear model of the F1 hybrid electric powertrain and the car’s longitudinal dynamics, especially focusing on the engine speed dependency, the dynamic behavior of the intake manifold and the turbocharger, as well as the engine cylinder deactivation. Thereafter, we formulate the time-optimal control problem and apply outer convexification. In Section 4 we present parameter studies where we compute and discuss the time-optimal control strategies for different energy budgets and operational limits, both over portions of lap and entire lap. Finally, conclusions are drawn in Section 5.

2. Modeling of the Formula 1 Powertrain

In this section we identify and validate a detailed low-level model of the F1 hybrid electric powertrain. A schematic of the system is shown in Figure 1 and its mathematical mean-value model formulation suited for numerical optimization is mostly based on or inspired by [31,47,57].

In addition to the energy storages onboard, i.e., the battery and the fuel tank, the power unit consists of an internal combustion engine equipped with an electrified turbocharger and two motor generator units: the MGU-K is mounted on the engine shaft and is used for electric boosting or regenerative braking, whilst the MGU-H is an electric motor that converts part of the heat contained in the hot exhaust gases into electrical energy or reduces the turbo lag. The engine system is composed of an intake manifold equipped with an intercooler and a throttle valve, a 1.6 L gasoline direct injection engine with cylinder deactivation capabilities, an exhaust manifold and a turbocharger. The turbocharger consists of a radial compressor, a radial turbine equipped with a waste-gate valve and the MGU-H mounted on the turbocharger shaft. Finally, the remaining part of the powertrain consists of a gearbox and a differential that connect the engine and the MGU-K to the wheels.

A schematic representation of how the state and the input variables of the considered system relate to each other, is displayed in the causality diagram shown in Figure 2. The system’s state variables are the car’s velocity v and its position s on the racetrack, the intake manifold pressure $p_{\mathrm{im}}$, the turbocharger rotational kinetic energy $E_{\mathrm{tc}}$, the battery energy $E_{\mathrm{b}}$ and fuel energy $E_{\mathrm{f}}$ consumed. In contrast to the intake manifold dynamic behavior, the exhaust manifold is modeled as a static subsystem. This modeling assumption will be motivated in Section 2.4. The system’s input variables are the power of the friction brakes $P_{\mathrm{brk}}$, the MGU-K power $P_{\mathrm{k}}$, the MGU-H power $P_{\mathrm{h}}$, the engaged gear g, the throttle valve position $u_{\mathrm{th}}$, the spark-advance position $u_{\mathrm{sa}}$, the cylinder fuel mass flow ${\dot{m}}_{\mathrm{f}},\mathrm{cyl}$, the cylinder deactivation variable ${\mathsf{\Psi }}_{\mathrm{e}}$ and the waste-gate valve position $u_{\mathrm{wg}}$. The powertrain model is partly taken or inspired by the works conducted at our institute or published by ourselves. Specifically, the energy recovery system model in Section 2.5 and the longitudinal dynamics in Section 2.6 are taken from [57]. The linear waste-gate model in Section 2.3 and the static approximation of the exhaust manifold in Section 2.4 are taken from our previous work [47], where the engine speed over the entire mission was given and we made use of integer variables and linear equations to model a high-performance power unit. For confidentiality reasons, the model parameters cannot be disclosed and the data shown on the plots are normalized and denoted by ${(·)}_0$. For instance, the normalized value of the turbine mass flow ${\dot{m}}_{\mathrm{t}}$ is denoted by ${\dot{m}}_{\mathrm{t}},0$. To quantify the modeling mismatch, we use the normalized mean absolute error defined as

$$e_{\mathrm{norm}}=\frac{1}{\max |x_{\mathrm{meas}}|·N}·\sum _{i=1}^N\parallel x_{\mathrm{meas}}\left[i\right]-x_{\mathrm{model}}\left[i\right]{\parallel }_1,$$

(1)

where N is the number of data points considered, $x_{\mathrm{model}}$ is the modeled value and $x_{\mathrm{meas}}$ is the measured one. The data displayed in the Figures 3, 5–10 are telemetry data acquired during a race week-end on a current F1 hybrid electric race car by our industrial partner. For confidentiality reasons, we are not allowed to disclose the exact experimental setup used to obtain them, i.e., sensor specification and placement. The majority of the data are acquired directly by means of sensors, e.g., intake manifold pressure, or are computed online by means of maps that were recorded on dedicated test-benches, e.g., compressor and turbine mass flow maps. Finally, all the validation plots shown in the following sections cover roughly a third of the racetrack and include three straights and two corners.

2.1. Intake Manifold

The intake manifold is modeled as an isothermal receiver [31], since the intercooler operation leads to temperature fluctuations over one lap of less than $1\%$, as Figure 3 shows. As a consequence, we assume a constant temperature ${\vartheta }_{\mathrm{im}}$ inside the intake manifold.

We then model the dynamics of the intake manifold pressure $p_{\mathrm{im}}$ as

$$\frac{\mathrm{d}}{\mathrm{d}t}{p}_{\mathrm{im}}\left(t\right)=\frac{R_{\mathrm{air}}·{\vartheta }_{\mathrm{im}}}{V_{\mathrm{im}}}·\left({\dot{m}}_{\mathrm{c}}\left(t\right)-{\dot{m}}_{\beta }\left(t\right)\right),$$

(2)

where $V_{\mathrm{im}}$ is the intake manifold volume, $R_{\mathrm{air}}$ is the specific gas constant of air, ${\dot{m}}_{\mathrm{c}}$ is the compressor mass flow and ${\dot{m}}_{\beta }$ is the engine air mass flow. The engine air mass flow is modeled as

$${\dot{m}}_{\beta }\left(t\right)=\frac{p_{\mathrm{im}}\left(t\right)}{R_{\mathrm{air}}·{\vartheta }_{\mathrm{im}}}·\frac{{\omega }_{\mathrm{e}}\left(t\right)}{4\pi }·V_{\mathrm{d}}·{\lambda }_{\mathrm{vol}}\left({\omega }_{\mathrm{e}}\left(t\right)\right)·u_{\mathrm{th}}\left(t\right),$$

(3)

where ${\omega }_{\mathrm{e}}$ is the engine speed, $V_{\mathrm{d}}$ is the engine displacement volume, ${\lambda }_{\mathrm{vol}}$ is the volumetric efficiency and $u_{\mathrm{th}}\in [0,1]$ is the throttle position. Since the throttle is placed after the plenum, the intake manifold pressure dependency [31] of the volumetric efficiency is neglected. The map ${\lambda }_{\mathrm{vol}}\left({\omega }_{\mathrm{e}}\right)$ that relates the volumetric efficiency to the engine speed is identified using the neural network technique presented in Appendix C. For confidentiality reasons, it cannot be displayed.

2.2. Engine Power and Cylinder Deactivation

The internal combustion engine power $P_{\mathrm{e}}$ is the sum of three contributions: the engine power stemming from the fuel combustion $P_{\mathrm{e},\mathrm{c}}$, the engine power losses due to the pistons’ friction $P_{\mathrm{e},\mathrm{fr}}$ and the engine pumping power $P_{\mathrm{e},\mathrm{p}}$ arising from the gas exchange in the intake and exhaust manifold. Therefore we model the engine power as

$$P_{\mathrm{e}}\left(t\right)=P_{\mathrm{e},\mathrm{c}}\left(t\right)+P_{\mathrm{e},\mathrm{fr}}\left(t\right)+P_{\mathrm{e},\mathrm{p}}\left(t\right).$$

(4)

Since the engine features cylinder deactivation, the engine combustion power $P_{\mathrm{e},\mathrm{c}}$ is the sum of all the cylinder-individual combustion power values:

$$P_{\mathrm{e}},\mathrm{c}\left(t\right)=\sum _{i=1}^{N_{\mathrm{cyl}}}{P}_{\mathrm{e},\mathrm{c},i}\left(t\right),$$

(5)

where $P_{\mathrm{e},\mathrm{c},i}$ denotes the engine combustion power delivered by cylinder $i\in \{1,\dots ,N_{\mathrm{cyl}}\}$. It is defined as

$$P_{\mathrm{e},\mathrm{c},i}\left(t\right)={\eta }_{{\omega }_{\mathrm{e}}}\left({\omega }_{\mathrm{e}}\left(t\right)\right)·{\eta }_{\mathrm{AF},i}\left(t\right)·{\eta }_{\mathrm{SA},i}\left(t\right)·{\dot{m}}_{\mathrm{f},i}\left(t\right)·H_{\mathrm{l}},$$

(6)

where $H_{\mathrm{l}}$ is the fuel lower heating value and ${\dot{m}}_{\mathrm{f},i}$ is the amount of fuel injected in cylinder i. The combustion efficiency ${\eta }_{{\omega }_{\mathrm{e}}}$ is engine speed dependent, whereas ${\eta }_{\mathrm{AF},i}$ and ${\eta }_{\mathrm{SA},i}$ are cylinder individual correction factors that capture the influence of the air-to-fuel ratio and spark-advance position, respectively. Since all cylinders share the same engine speed and, if active, we assume the same amount of fuel injected ${\dot{m}}_{\mathrm{f},\mathrm{cyl}}$, air-to-fuel ratio and spark-advance, (6) can be rewritten as

$$P_{\mathrm{e},\mathrm{c}}\left(t\right)={\eta }_{{\omega }_{\mathrm{e}}}\left({\omega }_{\mathrm{e}}\left(t\right)\right)·{\eta }_{\mathrm{AF},\mathrm{cyl}}\left(t\right)·{\eta }_{\mathrm{SA},\mathrm{cyl}}\left(t\right)·[{\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)·{\dot{m}}_{\mathrm{f},\mathrm{cyl}}\left(t\right)]·H_{\mathrm{l}},$$

(7)

where the cylinder deactivation variable ${\mathsf{\Psi }}_{\mathrm{e}}\in [0,N_{\mathrm{cyl}}]$ defines the number of active cylinders. This variable is assumed to be continuous to reduce the calculation complexity (more on that in Section 3.2). The total fuel mass flow injected ${\dot{m}}_{\mathrm{f}}$ is therefore

$${\dot{m}}_{\mathrm{f}}\left(t\right)={\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)·{\dot{m}}_{\mathrm{f},\mathrm{cyl}}\left(t\right),$$

(8)

while the change in fuel energy $E_{\mathrm{f}}$ consumption is modeled as

$$\frac{\mathrm{d}}{\mathrm{d}t}{E}_{\mathrm{f}}\left(t\right)={\dot{m}}_{\mathrm{f}}\left(t\right)·H_{\mathrm{l}}.$$

(9)

The maximum amount of fuel injected per cylinder is limited by the FIA regulations [1,2]:

$$0⩽{\dot{m}}_{\mathrm{f},\mathrm{cyl}}\left(t\right)⩽{\dot{m}}_{\mathrm{f}}^{\max }\left({\omega }_{\mathrm{e}}\left(t\right)\right)/N_{\mathrm{cyl}},$$

(10)

where the speed-dependent maximum allowed fuel injection in $[\mathrm{kg}/\mathrm{h}]$ is

$${\dot{m}}_{\mathrm{f}}^{\max }\left({\omega }_{\mathrm{e}}\left(t\right)\right)=\left\{\begin{array}{cc}{f}_{\mathrm{FIA},{\omega }_{\mathrm{e}}}\left({\omega }_{\mathrm{e}}\left(t\right)\right)=(0.009·{\omega }_{\mathrm{e}}\left(t\right)\left[\mathrm{rpm}\right]+5.5)\hfill & \mathrm{if}{\omega }_{\mathrm{e}}\left(t\right)⩽{\omega }_{\mathrm{e}}^{\mathrm{FIA}},\hfill \\ f_{\mathrm{FIA},\mathrm{const}}=100\hfill & \mathrm{if}{\omega }_{\mathrm{e}}\left(t\right)>{\omega }_{\mathrm{e}}^{\mathrm{FIA}},\hfill \end{array}\right.$$

(11)

where ${\omega }_{\mathrm{e}}^{\mathrm{FIA}}=10,500rpm$ is the engine speed threshold imposed by the FIA regulations. To implement this non-smooth constraint in the optimization framework, we rewrite and relax (11) and introduce the lifting variable ${\dot{\tilde{m}}}_{\mathrm{f}}^{\max }$ in the following manner:

$$\begin{array}{cc}\hfill {\dot{\tilde{m}}}_{\mathrm{f}}^{\max }\left(t\right)& ⩾0,\hfill \\ \hfill {\dot{\tilde{m}}}_{\mathrm{f}}^{\max }\left(t\right)& ⩽f_{\mathrm{FIA},\mathrm{const}},\hfill \\ \hfill {\dot{\tilde{m}}}_{\mathrm{f}}^{\max }\left(t\right)& ⩽f_{\mathrm{FIA},{\omega }_{\mathrm{e}}}\left({\omega }_{\mathrm{e}}\left(t\right)\right).\hfill \end{array}$$

(12)

Finally, using (12) we can rewrite (10) as

$$0⩽{\dot{m}}_{\mathrm{f},\mathrm{cyl}}\left(t\right)⩽{\dot{\tilde{m}}}_{\mathrm{f}}^{\max }\left(t\right)/N_{\mathrm{cyl}}.$$

(13)

The reformulation (12) is schematically drawn in Figure 4 and can be interpreted as follows: if the engine speed is below $10,500rpm$ the maximum amount of fuel that can be injected is constrained by the engine-speed dependent ramp $f_{\mathrm{FIA},{\omega }_{\mathrm{e}}}$; if the engine speed is above $10,500rpm$, the maximum amount of fuel is limited by the constant constraint $f_{\mathrm{FIA},\mathrm{const}}$.

The spark-advance efficiency depends on the spark advance position $u_{\mathrm{SA}}$ and the engine speed ${\omega }_{\mathrm{e}}$. The relationship is usually defined by a map ${\eta }_{\mathrm{SA}}=f_{\mathrm{SA}}(u_{\mathrm{SA}},{\omega }_{\mathrm{e}})$. We assume that for any engine speed we can realize a spark-advance efficiency ${\eta }_{\mathrm{SA}}\in [{\eta }_{\mathrm{SA}}^{\min },{\eta }_{\mathrm{SA}}^{\max }]$ and in post-processing compute the spark-advance position $u_{\mathrm{SA}}$ which realizes the desired efficiency. The engine speed efficiency is modeled as

$${\eta }_{{\omega }_{\mathrm{e}}}\left({\omega }_{\mathrm{e}}\right)=k_{1,\mathrm{e}}·{\omega }_{\mathrm{e}}^2\left(t\right)+k_{2,\mathrm{e}}·{\omega }_{\mathrm{e}}\left(t\right)+k_{3,\mathrm{e}},$$

(14)

where $k_{1,\mathrm{e}}$, $k_{2,\mathrm{e}}$ and $k_{3,\mathrm{e}}$ are coefficients subject to identification. Next, we introduce the definition of the air-to-fuel ratio in the cylinder [31]

$${\lambda }_{\mathrm{AF},\mathrm{cyl}}\left(t\right)=\frac{{\dot{m}}_{\beta }\left(t\right)/N_{\mathrm{cyl}}}{{\dot{m}}_{\mathrm{f},\mathrm{cyl}}\left(t\right)}·\frac{1}{{\sigma }_0},$$

(15)

where ${\sigma }_0$ is the stoichiometric constant. The main drawback of (15) is its singularity if no fuel is injected, i.e., for ${\dot{m}}_{\mathrm{f},\mathrm{cyl}}=0$. For this reason, we introduce the fuel-to-air ratio ${\varphi }_{\mathrm{AF}},\mathrm{cyl}$, which is defined as

$${\varphi }_{\mathrm{AF},\mathrm{cyl}}\left(t\right)=\frac{1}{{\lambda }_{\mathrm{AF},\mathrm{cyl}}\left(t\right)}=\frac{{\dot{m}}_{\mathrm{f},\mathrm{cyl}}\left(t\right)}{{\dot{m}}_{\beta }\left(t\right)/N_{\mathrm{cyl}}}·{\sigma }_0.$$

(16)

The fuel-to-air ratio ${\varphi }_{\mathrm{AF},\mathrm{cyl}}$ is always defined and bounded, since the engine air mass flow cannot be zero during operation. Furthermore, the fuel-to-air ratio mixture has an impact on the combustion efficiency, which is defined by the map ${\eta }_{\mathrm{AF},\mathrm{cyl}}=f_{\mathrm{AF}}\left({\varphi }_{\mathrm{AF},\mathrm{cyl}}\right)$. Similar to the engine speed efficiency maps, the fuel-to-air ratio efficiency maps usually present a peak efficiency at the optimal fuel-to-air ratio value and gradually decrease in an asymmetric way if the fuel-to-air ratio is larger or smaller than that [31]. To properly capture the asymmetric form of the map, the fuel-to-air ratio efficiency is fitted using the neural network techniques presented in Appendix C. If a cylinder is activated, i.e., ${\mathsf{\Psi }}_{\mathrm{e}}\ne 0$, the fuel-to-air ratio must lie in the allowed range [31], leading to the case distinction

$${\varphi }_{\mathrm{AF},\mathrm{cyl}}\left(t\right)\in \left\{\begin{array}{cc}[{\varphi }_{\mathrm{AF}}^{\min },{\varphi }_{\mathrm{AF}}^{\max }]\hfill & \mathrm{if}{\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)\ne 0,\hfill \\ \mathbb{R}\hfill & \mathrm{if}{\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)=0.\hfill \end{array}\right.$$

(17)

It is worth noticing that ${\varphi }_{\mathrm{AF},\mathrm{cyl}}$ can be any value if ${\mathsf{\Psi }}_{\mathrm{e}}=0$. Ideally, when the engine is completely shut off, i.e., ${\mathsf{\Psi }}_{\mathrm{e}}=0$, the fuel-to-air ratio should be zero because no fuel should be injected, i.e.,

$${\dot{m}}_{\mathrm{f},\mathrm{cyl}}\left(t\right)\in \left\{\begin{array}{cc}\left[0,{\dot{m}}_{\mathrm{f}}^{\max }\left({\omega }_{\mathrm{e}}\left(t\right)\right)/N_{\mathrm{cyl}}\right]\hfill & \mathrm{if}{\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)\ne 0,\hfill \\ 0\hfill & \mathrm{if}{\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)=0.\hfill \end{array}\right.$$

(18)

Since it is computationally expensive to implement (18), we neglect it for the following reason: If ${\mathsf{\Psi }}_{\mathrm{e}}\ne 0$, the fuel injected in the cylinders is constrained by both the FIA regulations (11) and the fuel-to-air ratio constraints (17); If ${\mathsf{\Psi }}_{\mathrm{e}}=0$, the fuel injected in the cylinders is only constrained by the FIA regulations, but it will not have an impact on any algebraic or dynamic equation of the model, e.g., on the combustion engine power in (7) or the fuel energy dynamics in (9), since ${\dot{m}}_{\mathrm{f},\mathrm{cyl}}$ is always multiplied with ${\mathsf{\Psi }}_{\mathrm{e}}$. Finally, we make use of so-called vanishing constraints [64] to rewrite (17) as

$$\begin{array}{cc}\hfill 0⩽& {\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)·\left({\varphi }_{\mathrm{AF},\mathrm{cyl}}\left(t\right)-{\varphi }_{\mathrm{AF}}^{\min }\right),\hfill \\ \hfill 0⩽& {\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)·\left({\varphi }_{\mathrm{AF}}^{\max }-{\varphi }_{\mathrm{AF},\mathrm{cyl}}\left(t\right)\right).\hfill \end{array}$$

(19)

The engine friction power $P_{\mathrm{e},\mathrm{fr}}$ is modeled as a quadratic function of the engine speed as

$$P_{\mathrm{e},\mathrm{fr}}\left(t\right)=k_{1,\mathrm{fr}}·{\omega }_{\mathrm{e}}^2\left(t\right)+k_{2,\mathrm{fr}}·{\omega }_{\mathrm{e}}\left(t\right)+k_{3,\mathrm{fr}},$$

(20)

where $k_{1,\mathrm{fr}}$, $k_{2,\mathrm{fr}}$ and $k_{3,\mathrm{fr}}$ are coefficients subject to identification. Finally, the engine pumping power $P_{\mathrm{e},\mathrm{p}}$ is modeled as

$$P_{\mathrm{e},\mathrm{p}}\left(t\right)=\left(p_{\mathrm{im}}\left(t\right)·u_{\mathrm{th}}\left(t\right)·k_{\mathrm{p}}-p_{\mathrm{em}}\left(t\right)\right)·{\omega }_{\mathrm{e}}\left(t\right)·\frac{V_{\mathrm{d}}}{4\pi },$$

(21)

where $p_{\mathrm{em}}$ is the exhaust manifold pressure and $k_{\mathrm{p}}$ is a coefficient subject to identification. Figure 5 shows the fitting of the engine combustion, friction and pumping power. Whilst all three power contributions match precisely the measurements data on the straights, the friction and pumping powers partially differ during the corners, i.e., where the engine is not operated at full load. Regarding the friction power, this might be attributable to the simple model which lacks additional (input or state) dependencies, whilst the difference in the pumping power might be due to physical phenomena that zero-dimensional models are not able to capture properly.

2.3. Turbocharger and Waste-Gate

The turbocharger’s state variable is its rotational kinetic energy $E_{\mathrm{tc}}$ defined as

$$E_{\mathrm{tc}}\left(t\right)=\frac{1}{2}·{\Theta }_{\mathrm{tc}}·{\omega }_{\mathrm{tc}}^2\left(t\right),$$

(22)

where ${\Theta }_{\mathrm{tc}}$ is the moment of inertia of the entire turbocharger and ${\omega }_{\mathrm{tc}}$ is the rotational speed of the shaft. The rotational kinetic energy evolves as a function of the turbine power $P_{\mathrm{t}}$, compressor power $P_{\mathrm{c}}$ and MGU-H power $P_{\mathrm{h}}$ acting on the shaft:

$$\frac{\mathrm{d}}{\mathrm{d}t}{E}_{\mathrm{tc}}\left(t\right)=P_{\mathrm{t}}\left(t\right)-P_{\mathrm{c}}\left(t\right)+P_{\mathrm{h}}\left(t\right).$$

(23)

In Ref. [47], it was shown that the compressor and turbine mass flows and powers can be accurately described as a function of the turbocharger rotational energy $E_{\mathrm{tc}}$, the intake manifold pressure $p_{\mathrm{IM}}$, the exhaust manifold pressure $p_{\mathrm{em}}$ and exhaust manifold temperature ${\vartheta }_{\mathrm{EM}}$. Therefore, these quantities can be represented with the following multi-dimensional maps $\mathcal{M}$:

$$\begin{array}{cc}\hfill {\dot{m}}_{\mathrm{c}}\left(t\right)& ={\mathcal{M}}_{{\dot{m}}_{\mathrm{c}}}\left(E_{\mathrm{tc}}\left(t\right),p_{\mathrm{im}}\left(t\right)\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill P_{\mathrm{c}}\left(t\right)& ={\mathcal{M}}_{P_{\mathrm{c}}}\left(E_{\mathrm{tc}}\left(t\right),p_{\mathrm{im}}\left(t\right)\right),\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\dot{m}}_{\mathrm{t}}\left(t\right)& ={\mathcal{M}}_{{\dot{m}}_{\mathrm{t}}}\left(E_{\mathrm{tc}}\left(t\right),p_{\mathrm{em}}\left(t\right),{\vartheta }_{\mathrm{em}}\left(t\right)\right),\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill P_{\mathrm{t}}\left(t\right)& ={\mathcal{M}}_{P_{\mathrm{t}}}\left(E_{\mathrm{tc}}\left(t\right),p_{\mathrm{em}}\left(t\right),{\vartheta }_{\mathrm{em}}\left(t\right)\right).\hfill \end{array}$$

(27)

This is a direct consequence of the assumption that the pressure and temperature at the compressor entrance, as well as the pressure after the turbine, are assumed to be constant and equal to the ambient conditions [47].

Both compressor and turbine are subject to operational limits. While the mechanical, blocking and choke limits of the compressor can be accounted for by introducing constant constraints on the turbocharger rotational speed and compressor mass flow, the compressor surge limit varies depending on the turbocharger rotational speed:

$$p_{\mathrm{surge}}\left(t\right)={\mathcal{M}}_{p_{\mathrm{surge}}}\left(E_{\mathrm{tc}}\left(t\right)\right).$$

(28)

In order to avoid compressor operation in the surge region, we introduce the constraint

$$p_{\mathrm{im}}\left(t\right)⩽p_{\mathrm{surge}}\left(t\right).$$

(29)

All the turbocharger maps are fitted using the neural network modeling technique presented in Appendix C. Since the turbine and the waste-gate mass flows are well approximated as orifices [31] and share the same exhaust manifold pressure and temperature, the waste-gate mass flow ${\dot{m}}_{\mathrm{wg}}$ is modeled as [47]

$${\dot{m}}_{\mathrm{wg}}\left(t\right)=u_{\mathrm{wg}}\left(t\right)·{\dot{m}}_{\mathrm{wg},\max },\max \left(t\right),$$

(30)

where $u_{\mathrm{wg}}\in [0,1]$ is the normalized waste-gate position and ${\dot{m}}_{\mathrm{wg},\max }$ is defined as

$${\dot{m}}_{\mathrm{wg},\max }\left(t\right)=k_{1,\mathrm{wg}}·{\dot{m}}_{\mathrm{t}}\left(t\right)+k_{2,\mathrm{wg}},$$

(31)

with $k_{1,\mathrm{wg}}$ and $k_{2,\mathrm{wg}}$ being parameters subject to identification. Figure 6 and Figure 7 show that the fitting precision of the turbocharger quantities and the waste-gate maximum mass flow is satisfactory, with very good agreement between the model and the measurements.

2.4. Exhaust Manifold

The exhaust manifold of the considered engine displays a small gas residence time and its volume is much smaller compared to the one of the intake manifold. As a consequence, we approximate the exhaust manifold to be quasi-static and adiabatic [31]. In Ref. [47], this simplification was shown to satisfactorily approximate the dynamic nature of a small plenum. Therefore, the exhaust manifold temperature ${\vartheta }_{\mathrm{em}}$ is modeled as an algebraic function

$${\vartheta }_{\mathrm{em}}\left(t\right)=\frac{1}{N_{\mathrm{cyl}}}·\sum _{i=1}^{N_{\mathrm{cyl}}}{\vartheta }_{\mathrm{e},\mathrm{out},i}\left(t\right),$$

(32)

where ${\vartheta }_{\mathrm{e},\mathrm{out},i}$ is the temperature of the gases exiting cylinder i. Since we assume all the active cylinders to share the same combustion conditions, and the inactive cylinders to share the same non-combustion conditions, (32) can be rewritten as

$${\vartheta }_{\mathrm{em}}\left(t\right)=\frac{1}{N_{\mathrm{cyl}}}·\left({\mathsf{\Psi }}_{\mathrm{e}}\left(t\right)·{\vartheta }_{\mathrm{comb}}^{\mathrm{ON}}\left(t\right)+(N_{\mathrm{cyl}}-{\mathsf{\Psi }}_{\mathrm{e}}\left(t\right))·{\vartheta }_{\mathrm{comb}}^{\mathrm{OFF}}\left(t\right)\right),$$

(33)

where the variables ${\vartheta }_{\mathrm{comb}}^{\mathrm{ON}}$ and ${\vartheta }_{\mathrm{comb}}^{\mathrm{OFF}}$ represent the engine-out temperature if the cylinder is active and inactive, respectively. The engine-out temperature of the active cylinders ${\vartheta }_{\mathrm{comb}}^{\mathrm{ON}}$ is modeled by superposing an enthalpic balance and a temperature increase induced by the spark-advance retardation [65,66,67]

$${\vartheta }_{\mathrm{comb}}^{\mathrm{ON}}\left(t\right)={\displaystyle \frac{{\tilde{H}}_{\mathrm{comb},\mathrm{enth}}^{\mathrm{ON}}\left(t\right)+{\tilde{H}}_{\mathrm{comb},\mathrm{SA}}^{\mathrm{ON}}\left(t\right)}{{\dot{m}}_{\mathrm{f}},\mathrm{cyl}\left(t\right)+{\dot{m}}_{\beta ,\mathrm{cyl}}\left(t\right)}},$$

(34)

with

$${\tilde{H}}_{\mathrm{comb},\mathrm{enth}}^{\mathrm{ON}}\left(t\right)={\dot{m}}_{\beta ,\mathrm{cyl}}\left(t\right)·{\vartheta }_{\mathrm{IM}}+{\dot{m}}_{\mathrm{f}},\mathrm{cyl}\left(t\right)·\left({\vartheta }_{\mathrm{fuel}}+{\eta }_{\mathrm{EM},\mathrm{enth}}·\frac{H_{\mathrm{l}}}{c_{\mathrm{p}}}\right),$$

(35)

where ${\vartheta }_{\mathrm{fuel}}$ is the temperature of the injected fuel, $c_{\mathrm{p}}$ is the heat capacity at constant pressure of the exhaust gases and ${\eta }_{\mathrm{EM},\mathrm{enth}}$ is a fitting coefficient that represents how much combustion enthalpy goes into the exhaust manifold in form of heat. The term accounting for the temperature increase due to the ignition delay is

$${\tilde{H}}_{\mathrm{comb},\mathrm{SA}}^{\mathrm{ON}}\left(t\right)={\eta }_{\mathrm{EM},\mathrm{SA}}·\frac{{\dot{m}}_{\mathrm{f}},\mathrm{cyl}\left(t\right)·H_{\mathrm{l}}·\left(1-{\eta }_{\mathrm{SA}}\left(t\right)\right)}{c_{\mathrm{p}}},$$

(36)

where ${\eta }_{\mathrm{EM},\mathrm{SA}}$ is a fitting coefficient that represents how much the exhaust manifold temperature is increased by the extra heat arising from the spark-advance retardation.

Finally, the engine-out temperature of an inactive cylinder is modeled as in [47] by

$${\vartheta }_{\mathrm{comb}}^{\mathrm{OFF}}\left(t\right)=k_1,{\mathrm{em}}^{\mathrm{OFF}}+k_2,{\mathrm{em}}^{\mathrm{OFF}}·{\dot{m}}_{\beta ,\mathrm{cyl}}\left(t\right),$$

(37)

where $k_1,{\mathrm{em}}^{\mathrm{OFF}}$ and $k_2,{\mathrm{em}}^{\mathrm{OFF}}<0$ are coefficients subject to identification. The model (37) captures that the engine-out temperature is colder for a higher engine air-mass flow. Figure 8 shows the comparison of the exhaust manifold temperature model to measurement data. The proposed model properly captures the behavior on the straights where the temperature is high, especially during the upshifts (recognizable as little peaks), whereas it overestimates the temperature in the last part of the corner. This might be attributable to the quasi-static modeling choice, which entails the absence of the receiver’s inertia.

Following the quasi-static modeling approach, the pressure in the exhaust manifold is implicitly defined through the algebraic mass balance

$${\dot{m}}_{\beta }\left(t\right)+{\dot{m}}_{\mathrm{f}}\left(t\right)={\dot{m}}_{\mathrm{t}}\left(t\right)+{\dot{m}}_{\mathrm{wg}}\left(t\right).$$

(38)

In Section 2.3, it was shown that the waste-gate mass flow can be approximated as a function of the turbine mass flow ${\dot{m}}_{\mathrm{t}}$ and the waste-gate position $u_{\mathrm{wg}}$, i.e.,

$${\dot{m}}_{\mathrm{wg}}\left(t\right)={\mathcal{M}}_{{\dot{m}}_{\mathrm{wg}}}\left({\dot{m}}_{\mathrm{t}}\left(t\right),u_{\mathrm{wg}}\left(t\right)\right),$$

(39)

and that the turbine mass flow ${\dot{m}}_{\mathrm{t}}$ is a function of the turbocharger rotational kinetic energy $E_{\mathrm{tc}}$, exhaust manifold pressure $p_{\mathrm{em}}$ and exhaust manifold temperature ${\vartheta }_{\mathrm{EM}}$, i.e.,

$${\dot{m}}_{\mathrm{t}}\left(t\right)={\mathcal{M}}_{{\dot{m}}_{\mathrm{t}}}\left(E_{\mathrm{tc}}\left(t\right),p_{\mathrm{em}}\left(t\right),{\vartheta }_{\mathrm{EM}}\left(t\right)\right).$$

(40)

If we combine (38) with (39), we obtain

$${\dot{m}}_{\mathrm{t}}\left(t\right)={\tilde{\mathcal{M}}}_{{\dot{m}}_{\mathrm{t}}}\left({\dot{m}}_{\beta }\left(t\right)+{\dot{m}}_{\mathrm{f}}\left(t\right),u_{\mathrm{wg}}\left(t\right)\right).$$

(41)

Finally, combining (41) with (40) the exhaust manifold pressure is found to be

$$p_{\mathrm{em}}\left(t\right)={\mathcal{M}}_{p_{\mathrm{em}}}\left({\dot{m}}_{\beta }\left(t\right)+{\dot{m}}_{\mathrm{f}}\left(t\right),u_{\mathrm{wg}}\left(t\right),E_{\mathrm{tc}}\left(t\right),{\vartheta }_{\mathrm{EM}}\left(t\right)\right).$$

(42)

The exhaust manifold pressure map is identified using the neural network technique presented in Appendix C and its fitting is shown in Figure 9. Despite its quasi-static modeling, the exhaust manifold pressure partially displays a dynamic behavior: this is due to its dependency on the engine air mass flow, which is a function of two dynamic variables, i.e., the intake manifold pressure and the car’s velocity.

2.5. Energy Recovery System

In this section, we model the energy recovery system using a high-level modeling approach [57]. The power conversion losses of the motor generator units (MGUs) from mechanical-to-electrical conversion and vice-versa are approximated using the quadratic functions

$$\begin{array}{cc}\hfill P_{\mathrm{k},\mathrm{dc}}\left(t\right)& ={\alpha }_{\mathrm{k}}·P_{\mathrm{k}}^2\left(t\right)+P_{\mathrm{k}}\left(t\right),\hfill \\ \hfill P_{\mathrm{k},\mathrm{dc}}\left(t\right)& ={\alpha }_{\mathrm{h}}·P_{\mathrm{h}}^2\left(t\right)+P_{\mathrm{h}}\left(t\right),\hfill \end{array}$$

(43)

where $P_{\mathrm{k},\mathrm{dc}}$ and $P_{\mathrm{h},\mathrm{dc}}$ are the electrical powers of the MGUs, $P_{\mathrm{k}}$ and $P_{\mathrm{h}}$ the mechanical power of the MGUs, and ${\alpha }_{\mathrm{k}}$ and ${\alpha }_{\mathrm{h}}$ coefficients subject to identification. The MGUs electrical actuator limits are considered, namely

$$\begin{array}{cc}\hfill P_{\mathrm{h},\mathrm{dc}}^{\min }& ⩽P_{\mathrm{h},\mathrm{dc}}\left(t\right)⩽P_{\mathrm{h},\mathrm{dc}}^{\min },\hfill \\ \hfill P_{\mathrm{k},\mathrm{dc}}^{\min }& ⩽P_{\mathrm{k},\mathrm{dc}}\left(t\right)⩽P_{\mathrm{h},\mathrm{dc}}^{\min }.\hfill \end{array}$$

(44)

The battery electrical losses can also be approximated using a quadratic model

$$P_{\mathrm{i}}\left(t\right)={\alpha }_{\mathrm{b}}·P_{\mathrm{b}}^2\left(t\right)+P_{\mathrm{b}}\left(t\right),$$

(45)

where $P_{\mathrm{i}}$ is the internal battery power, $P_{\mathrm{b}}=P_{\mathrm{h}},\mathrm{dc}+P_{\mathrm{k}},\mathrm{dc}+P_{\mathrm{aux}}$ is the terminal battery power, ${\alpha }_{\mathrm{b}}$ a coefficient subject to identification and $P_{\mathrm{aux}}$ a constant auxiliary power.

Figure 10 shows the satisfactory fitting precision of the MGU-H, MGU-K and battery power conversions. Finally, we assume a constant battery voltage [68,69] and model the battery energy dynamics as an open integrator [47,57], i.e.,

$$\frac{\mathrm{d}}{\mathrm{d}t}{E}_{\mathrm{b}}\left(t\right)=-P_{\mathrm{i}}\left(t\right).$$

(46)

2.6. Longitudinal Dynamics and Gearbox

In this section we present the longitudinal dynamics and the gearbox of the race car. The position s of the car on the circuit evolves according to its velocity v, hence

$$\frac{\mathrm{d}}{\mathrm{d}t}s\left(t\right)=v\left(t\right).$$

(47)

The lateral dynamics of the car are condensed [57] into a track-dependent maximum velocity profile $v_{\max }$ which constrains the velocity at each position of the racetrack

$$v\left(s\left(t\right)\right)⩽v^{\max }\left(s\left(t\right)\right).$$

(48)

The achievable velocity when the driver is not fully pressing the accelerator pedal, i.e., in the corners, is inferred from the telemetry data or using sophisticated simulations of the car’s lateral dynamics, while the maximum velocity constraint on the straight is set constant and equal to a value that is not achievable, e.g., $400\mathrm{km}/\mathrm{h}$. As an example, in Figure 11 we depict a measurement of the car’s velocity on the layout of the Bahrain International Circuit.

The longitudinal dynamics are modeled as [57]

$$m_{\mathrm{car}}·\frac{\mathrm{d}}{\mathrm{d}t}v\left(t\right)=F_{\mathrm{tr}}\left(t\right)-F_{\mathrm{a}}\left(t\right)-F_{\mathrm{r}}$$

(49)

with

$$\begin{array}{cc}\hfill F_{\mathrm{tr}}\left(t\right)& =\frac{P_{\mathrm{tr}}\left(t\right)}{v},\hfill \\ \hfill F_{\mathrm{a}}\left(t\right)& =\left(k_{\mathrm{a}}+k_{\mathrm{curve}}·\gamma \left(s\left(t\right)\right)-k_{\mathrm{drs}}·{\xi }_{\mathrm{drs}}\left(s\left(t\right)\right)\right)·v^2\left(t\right),\hfill \\ \hfill F_{\mathrm{r}}& =m_{\mathrm{car}}·G·k_{\mathrm{r}},\hfill \end{array}$$

(50)

where $m_{\mathrm{car}}$ is the car mass, $k_{\mathrm{a}}$, $k_{\mathrm{curve}}$, $k_{\mathrm{drs}}$ and $k_{\mathrm{r}}$ are the aerodynamic, curvature, drag reduction (DRS) and rolling coefficient subject to identification, $\gamma \left(s\left(t\right)\right)$ is a position-dependent curvature parameter that accounts for the aerodynamic change during cornering of open wheeled cars, ${\xi }_{\mathrm{drs}}\left(s\left(t\right)\right)$ is a position-dependent boolean (where ${\xi }_{\mathrm{drs}}=1$ in DRS zones, otherwise 0), G is the gravitational acceleration and $P_{\mathrm{tr}}$ is the traction power. The traction power is defined as

$$P_{\mathrm{tr}}\left(t\right)=P_{\mathrm{p}}\left(t\right)-P_{\mathrm{brk}}\left(t\right),$$

(51)

where $P_{\mathrm{brk}}\in \left[0,P_{\mathrm{brk}}^{\max }\right]$ is the power dissipated in the friction brakes and $P_{\mathrm{p}}$ is the propulsive power modeled as [57]

$$P_{\mathrm{p}}\left(t\right)=k_{1,\mathrm{sl}}·P_{\mathrm{u}}^2\left(t\right)+k_{2,\mathrm{sl}}·P_{\mathrm{u}}\left(t\right),$$

(52)

where $k_{1,\mathrm{sl}}$ and $k_{2,\mathrm{sl}}$ are the slip coefficients subject to identification and $P_{\mathrm{u}}$ is the total power unit power, i.e., $P_{\mathrm{u}}=P_{\mathrm{e}}+P_{\mathrm{k}}$. More details on the longitudinal model and its validation can be found in [57]. Finally, the power unit is equipped with an eight-speed gearbox, i.e., $N_{\mathrm{g}}=8$. The overall transmission of the gear g from the engine shaft to the car’s wheel is denoted by $\Gamma \in \{{\Gamma }_1,\dots ,{\Gamma }_{N_{\mathrm{g}}}\}$. Therefore the engine speed whit the gear $g\in \{1,\dots ,N_{\mathrm{g}}\}$ engaged is related to the velocity by

$${\omega }_{\mathrm{e}}\left(t\right)={\Gamma }_g·v\left(t\right)$$

(53)

and the engine speed is constrained for drivability and regulatory reasons [1,2] as

$${\omega }_{\mathrm{e}}^{\min }⩽{\omega }_{\mathrm{e}}\left(t\right)⩽{\omega }_{\mathrm{e}}^{\max }.$$

(54)

3. Optimal Control Problem

In this section, we state the optimal control problem. The state variables are the position s, the car’s velocity v, the intake manifold pressure $p_{\mathrm{im}}$, the turbocharger’s rotational kinetic energy $E_{\mathrm{tc}}$, the battery energy $E_{\mathrm{b}}$ and the fuel energy $E_{\mathrm{f}}$. The input variables are the throttle position $u_{\mathrm{th}}$, the waste-gate position $u_{\mathrm{wg}}$, the cylinder fuel mass flow ${\dot{m}}_{\mathrm{f}},\mathrm{cyl}$, the cylinder deactivation variable ${\mathsf{\Psi }}_{\mathrm{e}}$, the spark-advance efficiency ${\eta }_{\mathrm{SA}}$, the MGU-K power $P_{\mathrm{k}}$, the MGU-H power $P_{\mathrm{h}}$, the braking power $P_{\mathrm{brk}}$ and the engaged gear g.

3.1. Objective and Space Domain

The objective of the optimal control problem is to minimize the lap time

$$min{\int }_0^T\mathrm{d}t,$$

(55)

subject to the system dynamics and constraints presented in Section 2, as well as the energy consumption targets imposed on the fuel $\Delta E_{\mathrm{f}},\mathrm{target}$ and on the battery $\Delta E_{\mathrm{b},\mathrm{target}}$ as terminal constraints:

$$\begin{array}{cc}\hfill E_{\mathrm{b}}\left(T\right)& ⩾E_{\mathrm{b}}\left(0\right)+\Delta E_{\mathrm{b}},\mathrm{target},\hfill \\ \hfill E_{\mathrm{f}}\left(T\right)& ⩽E_{\mathrm{f}}\left(0\right)+\Delta E_{\mathrm{f}},\mathrm{target}.\hfill \end{array}$$

(56)

As presented in Section 2.6, the track-dependent parameters, e.g., the maximum velocity profile $v^{\max }\left(s\left(t\right)\right)$ or the track curvature $\gamma \left(s\left(t\right)\right)$, are given as a function of the position $s\left(t\right)$ on the track. Therefore, we reformulate the optimal control problem in space domain with the following change of variable based on (47)

$$v\left(t\right)=\frac{\mathrm{d}s\left(t\right)}{\mathrm{d}t}⟶\mathrm{d}t\left(s\right)=\frac{\mathrm{d}s}{v\left(s\right)},$$

(57)

where $\mathrm{d}s$ and $\mathrm{d}t$ are the differential variables in space and in time domain, respectively. While the algebraic equations and constraints do not change, the differential equations change as follows:

$$\frac{\mathrm{d}}{\mathrm{d}t}x\left(t\right)=F\left(x\left(t\right),u\left(t\right),\Gamma \left(t\right)\right)⟶\frac{\mathrm{d}}{\mathrm{d}s}x\left(s\right)=\tilde{F}\left(x\left(s\right),u\left(s\right),\Gamma \left(s\right)\right)=\frac{F\left(x\left(s\right),u\left(s\right),\Gamma \left(s\right)\right)}{v\left(s\right)},$$

(58)

where $F\left(x(.),u(.),\Gamma (.)\right)$ is the right-hand side of the differential equation in time domain and $\tilde{F}\left(x(.),u(.),\Gamma (.)\right)$ is its counterpart in space domain. Finally, the objective function changes to

$$min{\int }_0^T\mathrm{d}t⟶min{\int }_0^S\frac{\mathrm{d}s}{v\left(s\right)},$$

(59)

where T is the time necessary to cover the considered track length S.

3.2. Outer Convexification and Rounding

The model presented in Section 2 contains an integer variable, i.e., the engaged gear. Since mixed-integer optimization is computationally expensive [70], we convexify the model with respect to the integer variables using the outer convexification methodology [71,72,73]. Specifically, we introduce a control variable ${\delta }_g\left(s\right)\in \{0,1\}$ for each gear choice $\Gamma \in \{{\Gamma }_1,\dots ,{\Gamma }_{N_{\mathrm{g}}}\}$, with the bijection

$$\Gamma \left(s\right)={\Gamma }_g⟺{\delta }_g\left(s\right)=1.$$

(60)

Furthermore, we impose the variables ${\delta }_1,{\delta }_2,\dots ,{\delta }_{N_{\mathrm{g}}}$ to fulfill the special ordered set type 1 property (SOS-1), i.e.,

$$\sum _{g=1}^{N_{\mathrm{g}}}{\delta }_g\left(s\right)=1,{\delta }_g\left(s\right)\in \{0,1\},\forall g.$$

(61)

Next, we relax the problem by making the integer variable continuous ${\widehat{\delta }}_g$, i.e.,

$${\widehat{\delta }}_g\left(s\right)\in [0,1]\forall g.$$

(62)

The outer convexified differential equations of the system read as

$$\frac{d}{ds}x\left(s\right)=\sum _{g=1}^{N_{\mathrm{g}}}{\widehat{\delta }}_g\left(s\right)·\tilde{F}\left(x\left(s\right),u\left(s\right),{\Gamma }_g\right),$$

(63)

while the inequality constraints of the form

$$0⩽h\left(s,x\left(s\right),u\left(s\right),\Gamma \left(s\right)\right),$$

(64)

read as

$$0⩽\sum _{g=1}^{N_{\mathrm{g}}}{\widehat{\delta }}_g\left(s\right)·h\left(s,x\left(s\right),u\left(s\right),{\Gamma }_g\right).$$

(65)

The detailed formulation of the outer convexified optimal control problem is presented in Appendix A. Finally, to construct an integer ${\delta }_g$ sequence of the gear command from the relaxed optimal one ${\widehat{\delta }}_g$, a rounding strategy is applied. Due to the conditions imposed by the equality constraint in (61), the rounding strategy should not violate the SOS-1 property. Therefore, we use the rounding methodology proposed in [74] that reads as

$${\alpha }_{g,i}=\sum _{m=0}^i{\widehat{\delta }}_{g,m}\Delta s_m-\sum _{m=0}^{i-1}{\delta }_{g,m}\Delta s_m$$

(66)

$${\delta }_{g,i}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{\alpha }_{g,i}⩾{\alpha }_{k,i}\forall k\ne g\mathrm{and}g<k\forall k:{\alpha }_{g,i}={\alpha }_{k,i},\hfill \\ 0\hfill & \mathrm{else},\hfill \end{array}\right.$$

(67)

where ${\alpha }_{g,i}$ is an auxiliary variable of gear g at point i introduced to make (67) more readable, $\Delta s_m=s_{m+1}-s_m$ is the distance between two space points and i is the space point at which the rounding is sought. This recursive methodology simultaneously takes into account all $N_{\mathrm{g}}$ gears up to $i-1$ and favors the choice of the smallest gear g with the largest ${\alpha }_{g,i}$.

Regarding the engine cylinder deactivation, the variable ${\mathsf{\Psi }}_{\mathrm{e}}$ is treated as continuous, although it is integer by nature. This is done to reduce the calculation time and complexity, i.e., to avoid introducing $N_{\mathrm{g}}\times N_{\mathrm{cyl}}$ relaxed integers at each time step instead of $N_{\mathrm{g}}$ only. In addition to that, the solution to this problem can be used as a benchmark to compare the results obtained with an integer number of active cylinders and assess the resulting suboptimality. Finally, for the sake of completeness, we compare in Appendix B the relaxed and rounded nonlinear solution to the original mixed-integer nonlinear program solution. The results highlight that the trajectories of the two solutions are close to identical and that the lap time difference is negligible, thereby justifying the outer convexification approach.

3.3. Numerical Solution Method

We discretize the optimization problem using the multiple shooting method and the standard explicit Euler integration method, which approximates the continuous derivative of the variable $x\left(s\right)$ as

$$x[k+1]=x\left[k\right]+\Delta s·\tilde{F}\left(x\left(k\right),u\left(k\right)\right),$$

(68)

where $\Delta s$ is the space domain discretization length, $k=\{0,1,2,\dots ,N-1\}$ is the discrete space variable, $N=\mathrm{round}(S/\Delta s)$ is the number of discrete steps and S is the considered track length. The optimization problem is implemented in MATLAB [75] using the symbolic framework CasADi [76] and solved with IPOPT [77]. Furthermore, the mixed-integer nonlinear program (MINLP) solution presented in Appendix B is computed using BONMIN with the branch-and-bound algorithm [78]. In [47], it was shown and validated that sampling frequencies larger than $15\mathrm{Hz}$ were able to properly capture the intake manifold and turbocharger dynamics. Since the optimization problem is described in space domain, we pre-compute a variable space domain discretization length $\Delta s\left[k\right]$ for a specific (time domain) sampling frequency $f_{\mathrm{s}}$ using a nominal car’s velocity trajectory $v_{\mathrm{nom}}$ measured during a race week-end on a particular circuit:

$$\Delta t\left[k\right]=\frac{\Delta s\left[k\right]}{v_{\mathrm{nom}}\left[k\right]}=\frac{1}{f_{\mathrm{s}}}⟶\Delta s\left[k\right]=\frac{1}{f_{\mathrm{s}}}·v_{\mathrm{nom}}\left[k\right].$$

(69)

While it is most likely that the optimized velocity profile will not be the same as the nominal one in each point of the track, the sampling frequency does not deviate more than 0.9–1.2 Hz on average for entire lap optimizations. The computational time to solve both the relaxed and the rounded nonlinear program solutions on a desktop computer varies from minutes to hours and depends on the chosen boundary conditions and the space horizon length. As an example, if only a small portion of the lap is considered, e.g., a straight with two corners, the computational time is between 5 and 20 min. If the entire lap is considered, the computational time varies between 2 and 4 h. Finally, the MINLP displays computational times up to 100 times slower compared to the presented methodology.

4. Results

In this section, we show the results of the time-optimal control problem presented in Section 3 for racing scenarios on the Bahrain International Circuit, a lap of which takes roughly $90\mathrm{s}$. The racetrack under consideration influences the vehicle dynamics through its specific maximum velocity profile, path curvature and DRS zones. Our non-causal optimization framework is not based on a driving cycle, since the velocity profile and the propulsive power are subject to optimization. First, we analyze over a portion of the lap the optimal gearshift and power unit control from a low-level point of view. Second, we extend the analysis to entire lap results.

4.1. Detailed Analysis on a Portion of the Lap

In this subsection we present several case studies where we apply the proposed optimization framework to practically relevant strategy and design decisions. To facilitate a detailed analysis, we consider only a portion of the lap. We start by comparing the optimal control of the power unit for two different battery recharge targets, which is a strategy parameter that can be changed by the driver during a race. Thereafter, we add a constraint on the compressor outlet temperature to showcase how the presented framework can support the design choice of the intercooler dimension.

Figure 12 shows the most significant optimal state and input trajectories for two different battery recharge targets and the same fuel consumption target. The considered portion of lap is the section between the points A and B in Figure 11. It can be noticed that the gearshift strategy is different: not only is the 8th gear solely engaged for the battery discharge case ($\Delta E_{\mathrm{b},\mathrm{target}}=-0.3\mathrm{MJ}$), but also several upshifts occur earlier compared to the battery recharge case ($\Delta E_{\mathrm{b},\mathrm{target}}=0.5\mathrm{MJ}$). This trend can be attributed to the abundance of electrical energy at disposal: the gears are selected to maximize the overall engine power and therefore the propulsive power. For the battery recharge case, however, the trade-off between the engine power and the MGU-H recuperation is more pronounced. The fact of not shifting into the 8th gear keeps the engine speed high and therefore the pressure in the intake manifold lower (compared to the battery discharge case). As a consequence, the compressor power is lower and more power can be recuperated by the MGU-H. The power unit control around the upshifts occurs each time in a similar manner: the MGU-H recuperates less to accelerate the turbocharger shaft, whilst the throttle valve is operated to achieve the optimal fuel-to-air ratio value. The battery recharge target also influences the waste-gate operation. In conventional gasoline turbocharged engines, the waste-gate is operated to diminish the turbine power and thereby controls the intake manifold pressure [79]. By contrast, for the considered high-performance power unit the waste-gate is used to increase the engine power by diminishing the exhaust manifold pressure. The turbine power extraction is reduced and therefore the MGU-H power recuperation goes roughly to zero (when no upshift is taking place). The intake manifold pressure is not affected, the engine pumping power increases and so does the overall engine power. In addition to that, the MGU-K operation is also influenced by the battery availability. Not only the MGU-K cut occurs later on the straight if more battery energy can be deployed, but also the slope is different: the analytical explanation for this can be found in [58]. Finally, the optimal behavior in the corner is analyzed. The optimal cylinder deactivation is mainly the same, except at the corner exit. In this region, for the battery discharge case, slightly more fuel is injected. At the same time the spark-advance is slightly retarded and produces a lower spark-advance efficiency. This causes a small increase in the exhaust manifold temperature and pressure, allowing the MGU-H to recuperate more. The procedure of converting a small portion of fuel energy into electrical energy through the increase of the exhaust manifold temperature is possible in the cases where the battery recharge targets allow it. For the battery recharge case, this operation only occurs for a very short amount of time, since every drop of fuel counts and the gearshift strategy is not aimed at maximizing the engine efficiency, but rather finding the optimal trade-off between the internal combustion engine and the energy recovery system.

Next, we investigate the suboptimality entailed by a smaller intercooler choice in Figure 13. In such a case, to keep the same intake manifold temperature, the compressor outlet temperature needs to be lowered. As an example we include a constraint to decrease the maximum outlet temperature of the compressor by $3\%$ compared to the unconstrained case, i.e., ${\vartheta }_{\mathrm{c}}^{\lim }⩽0.97·max\left({\vartheta }_{\mathrm{c}}^{\mathrm{unlim}}\right)$. The compressor outlet temperature was modeled using the assumptions presented in Section 2.3 and fitted with the methodology of Appendix C, i.e., ${\vartheta }_{\mathrm{c}}\left(s\right)={\mathcal{M}}_{{\vartheta }_{\mathrm{c}}}\left(E_{\mathrm{tc}}\left(s\right),p_{\mathrm{im}}\left(s\right)\right)$. To better visualize the transient phenomena during the gearshifts, thin vertical lines have been inserted at the upshift locations. The battery recharge target and fuel consumption target are the same for both scenarios. It can be noticed at first glance that the intake manifold and turbocharger speed trajectories are different and almost everywhere lower compared to the unconstrained case. This is mainly due to their dependency on the compressor outlet temperature. In addition, the engine air and fuel mass flows, as well as the gearshift strategy, differ markedly to case studies previously analyzed. Since the intake manifold pressure is (indirectly) limited, the gearshifts occur at higher engine speeds to avoid large engine air mass flow deficits after the upshift. Despite the throttle valve being completely opened after the upshifts, the fuel mass flow is limited by the available engine air mass flow and the allowed fuel-to-air ratio range. Consequently, the injected fuel displays a transient behavior induced by the intake manifold inertia. The influence of the gradual increase of both the fuel injection and the air mass flow are visible in the exhaust manifold quantities: The pressure is lower because of a lower total mass flow through the turbine, whilst the temperature is higher due to the larger fuel-to-air ratio. Finally, the lap time loss is roughly $12\mathrm{ms}$ and arises from the engine power deficits occurring after the upshifts.

4.2. Trends in Entire Lap Solutions

In this subsection we look at entire lap solutions with a focus on the gearshift strategy, power unit actuation and lap time trends. First, we discuss the engine speed thresholds at which the upshifts occur. Second, we extend the investigations of Section 4.1 for entire lap results. Third, we analyze the power unit optimal control on the straights for a set of battery recharge targets.

The engine speed thresholds at which the upshifts occur for each gear are shown in Figure 14 for three different battery recharge targets. In the left plot the raw data are depicted, in the right plot the averaged values are displayed. Regarding the raw data, the points are spread and no clear pattern on the engine speed threshold can be inferred. Instead, for the averaged data, two regions can be distinguished. In the first region, the upshift thresholds from 2nd to 3rd and from 3rd to 4th gear are the same, irrespective of the battery budget. This could be attributed to the fact that these two specific upshifts always occur in the corners, where the power unit is (for the considered energy budgets) controlled in a similar way, as also discussed in Section 4.1. In the second region on the other hand, the upshift thresholds are related to the energy budget. This trend could be explained by the fact that upshifting at higher engine speeds favors the electrical recuperation via MGU-H at the expenses of lower engine efficiency, as seen and discussed in the case study of Figure 12.

In Section 4.1 it was analyzed in detail how the compressor outlet temperature influences the achievable lap time. For this reason, the lap time losses over an entire lap for the limited compressor outlet temperature for five different battery recharge targets, are shown in Figure 15. The trend shows that irrespective of the battery recharge target, limiting the compressor outlet temperature can yield a significant lap time loss. This result clearly quantifies how lap time sensitive the choice of the compressor material or the intercooler design can be.

The average internal combustion engine and energy recovery system powers for different battery recharge targets but the same fuel target are shown in Figure 16. While the internal combustion engine and the energy recovery system powers in the corners influence the energy at disposal to the power unit when exiting the corners, they do not have a direct impact on the achievable lap time. For this reason, each power entry is averaged over all straights. The trends displayed in Figure 16 are in line with the ones seen in Section 4.1. The battery recharge target has an influence on the time interval the waste-gate is kept open, and therefore directly influences the resulting engine pumping power. The latter is shown to have the largest impact on the overall engine power, since the sum of the engine combustion and friction power are, on average, constant for all the battery budgets (for the considered fuel consumption target). Opening the waste-gate comes at the expense of a lower MGU-H recuperation: in fact, the average MGU-H recuperation (on straight) is roughly zero if the waste-gate is kept always open. A similar trend is observed for the MGU-K operation, whose average value becomes larger if more battery energy can be depleted. Overall, the power drawn from the battery, i.e., the sum of the MGU-K and MGU-H powers, becomes very small in magnitude if an aggressive battery recharge target has to be met: The power recuperated from the extra enthalpy contained in the hot exhaust gases is directly fed from the MGU-H to the MGU-K. From a lap time point of view, it can be noticed how the lap time sensitivity changes: the lap time loss from $\Delta E_{\mathrm{b},\mathrm{target}}=0.0\mathrm{MJ}$ to $\Delta E_{\mathrm{b},\mathrm{target}}=2.0\mathrm{MJ}$ is more than twice compared to the one from $\Delta E_{\mathrm{b},\mathrm{target}}=-2.0\mathrm{MJ}$ to $\Delta E_{\mathrm{b},\mathrm{target}}=0.0\mathrm{MJ}$. Among other effects this can be explained by the nonlinearity of the lap time integration in the objective and the longitudinal dynamics.

Finally, the lap time loss between $\Delta E_{\mathrm{b},\mathrm{target}}=-2.0\mathrm{MJ}$ and $\Delta E_{\mathrm{b}},\mathrm{target}=2.0\mathrm{MJ}$ for several fuel consumption targets is shown in Figure 17. It can be noticed how the lap time loss increases nonlinearly if less fuel is used over the lap. This is a further proof of the complexity and inherent cross-couplings between the internal combustion engine and the energy recovery system of the considered hybrid electric powertrain.

5. Conclusions

In this paper we presented an optimization framework to compute the time-optimal low-level control and gearshift strategies for the Formula 1 hybrid electric powertrain. We put the focus on the engine speed dependency of the internal combustion engine components and proposed a continuous nonlinear model to describe the engine cylinder deactivation. Furthermore, we applied convex relaxations to model the non-smooth maximum fuel mass flow regulation and captured the highly nonlinear turbocharger maps with neural network techniques. Finally, we solved the time-optimal control problem efficiently by applying outer convexification to the integer decision variable introduced by the gear selection. To showcase the optimization framework, we performed several case studies for different energy budgets. The results underlined the significant existing coupling between the internal combustion engine and the energy recovery system: The gearshift strategy, the MGU-K and MGU-H operation, the air-to-fuel ratio operational range, the waste-gate and throttle valve actuation, as well as the intake manifold pressure evolution were found to be strongly dependent on the energy budgets. Finally, we demonstrated the effectiveness of our framework as a design tool. For example, if the compressor outlet temperature has to be reduced by $3\%$ due to a change in the intercooler design, the lap time can increase by up to $130ms$ for the considered conditions. The results presented in this work underline that the framework can be used as a tool to benchmark real-time control strategies, to generate reference trajectories and to develop control heuristics.