Multiphase mixed-integer nonlinear optimal control of hybrid electric vehicles ✩

This study considers the problem of computing a non-causal minimum-fuel energy management strategy for a hybrid electric vehicle on a given driving cycle. Specifically, we address the multiphase mixed-integer nonlinear optimal control problem that arises when the optimal gear choice, torque split and engine on/off controls are sought in off-line evaluations. We propose an efficient model by introducing vanishing constraints and a phase specific right-hand side function that accounts for the different powertrain operating modes. The gearbox and driveability requirements translate into combinatorial constraints. These constraints have not been included in previous research; however, they are part of the algorithmic framework for this investigation. We devise a tailored algorithm to solve this problem by extending the combinatorial integral approximation ( CIA ) technique that breaks down the original mixed-integer nonlinear program into a sequence of nonlinear programs and mixed-integer linear programs, followed by a discussion of its approximation error. Finally, numerical results illustrate the proposed algorithm in terms of solution quality and run time.


Introduction
Automotive manufacturers and research centers have been significantly investing resources and efforts into the development of alternative powertrain technologies to lower fuel consumption and pollutant emissions in passenger and commercial vehicles. Hybrid electric vehicles (HEVs) represent a concrete answer to address these problems, since they can reduce greenhouse gas emissions and fuel consumption while providing a high-quality ride. Notwithstanding this, the growing complexity and degrees of freedom of current hybrid powertrain architectures impose a tailored supervisory energy management strategy (EMS) to unleash the full potential of the HEV in terms of the fuel economy and driveability. Therefore, to gain insight on how to implement the control strategy in on-line applications, or to compare different powertrain architectures in a fair and unbiased way, the off-line EMS based on the solution of an optimal control problem ✩ The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Peng Shi under the direction of Editor Thomas Parisini. so-called vanishing constraints and combinatorial constraints for integer controls.
The above problem will be specified with all of its variables and constraints in more detail in Sections 2 and 3 .
Related work. Similar problems to Problem 1 were investigated with a fixed time horizon and an on-line EMS setting based on past research (Guzzella & Sciarretta, 2013;Rizzoni & Onori, 2015) that serves as a comprehensive overview.
This paper focuses on off-line applications. In this context, researchers have successfully solved nonlinear OCPs applied to the EMS of HEVs through optimization techniques ranging from dynamic programming to direct and indirect methods. Discrete dynamic programming has been extensively used for solving the nonlinear EMS problem due to its straightforward implementation and global optimality guarantees, with prominent examples in Elbert et al. (2015) and Wang et al. (2015). However, the curse of dimensionality prevents solving problems with a number of states that are generally greater than one or two. Indirect methods stem from the application of Pontryagin's maximum principle, which provides the necessary conditions for optimality. The solution of the OCP is down to the solution of a two-point boundary value problem or a multi-point boundary value problem in the case of active path constraints; however, this requires a priori knowledge of the constrained arcs. As a result, this restricts indirect methods to problems that are generally not influenced by the active path constraints. Examples can be found in Kim et al. (2011), Salazar et al. (2017) and Sciarretta et al. (2015). Recent advancements in the field of direct methods led to a widespread use of these techniques to solve OCPs. This is due to the possibility to treat a general problem class that is characterized by hundreds of differential and algebraic states. Convex optimization methods have been proven to be beneficial to solve the off-line EMS problem (Murgovski et al., 2012;Robuschi et al., 2019). However, their main drawback is the simplification of nonlinearities when it is applied as a linear or quadratic model. To cope with the nonlinear effects, in Limebeer et al. (2014) and Wei et al. (2017) the authors propose the direct transcription of the OCP into a nonlinear program (NLP) that can be solved by using standard NLP solvers. Both approaches deal with mixed-integer problems, which arise when both continuous and discrete variables are embedded into the OCP. This leads to NP-hard problems that are computationally intractable for standard solvers when considering long time horizons. Bengea and DeCarlo (2005) and Sager (2006Sager ( , 2009 have independently developed the embedding transformation technique that is also called combinatorial integral approximation (CIA) decomposition for solving MIOCPs. It consists of solving the NLP with a dropped integrality constraint before approximating the relaxed controls in the CIA problem, which is a mixed-integer linear program (MILP). Combinatorial constraints such as the maximum switching number conditions (Sager et al., 2011) and the dwell time constraints (Ali & Egerstedt, 2018;Zeile et al., 2020) have been discussed in various settings and they can be applied in the CIA context as additional constraints. There are several studies about control theory in the automotive field that builds upon the CIA decomposition (Böhme & Frank, 2017;Kirches et al., 2013;Meyer et al., 2016); however, most neglect the combinatorial constraints by using rounding schemes such as Sum-Up-Rounding (SUR). A rare example application of multiphase MIOCP can be found in Bonami et al. (2013).
Contributions. This study investigates Problem 1 under realworld requirements, specifically: • The powertrain operates in different modes depending on a given speed profile, which imposes the multiphase setting of the ODE. • The dual clutch gearbox allows only a specific switching structure that this study proposes for mode transition constraints.
• Switching between the electric and hybrid driving mode during arbitrarily short periods of time is impossible, which translates into minimum dwell time constraints.
We design a novel generalized CIA decomposition algorithm that uses several NLP and MILP steps with the idea to construct a feasible solution with a promising objective value for complex MIOCPs that entail multiphase, vanishing, state, and combinatorial constraints. To achieve this, the algorithm inherits the property of a vanishing optimality gap with a vanishing grid length from the original CIA technique.
Paper structure. Section 2 describes the powertrain model with its variables and constraints. Section 3 introduces the combinatorial constraints. An algorithmic framework to solve multiphase MIOCPs is presented in Section 4. Finally, the numerical case studies are discussed followed by the conclusions in Sections 5 and 6 , respectively.

Model description
This section presents the powertrain shown in Fig. 1. It consists of an internal combustion engine (ICE), an electric motor (EM) that provides boosting and regenerative braking, and a second electric motor (EM2) connected to the ICE through a belt. This can be used to recharge the battery when the vehicle stands still. The engine is connected to a 7-speed dual clutch gearbox, while the EM is coupled to the output shaft of the gearbox with an additional gear set. The final drive (FD) and the differential transmit the propulsive power to the wheels. The fuel tank and the battery are used for on-board energy storage.
In order to correctly evaluate the fuel consumption and the battery's state-of-charge while retaining a simple and fast estimation, we use a backward quasi-static modeling approach (Gao et al., 2007;Guzzella & Sciarretta, 2013) to describe the noncausal relationships between the powertrain subsystems. By making this choice, the number of states needed to describe the powertrain were reduced. We consider the speed profile v(t) of a given driving cycle as an exogenous variable and drop the driver model that would have otherwise been necessary to follow a reference speed profile; thus, reducing the complexity of the HEV model. The efficiencies and parameters of the main subsystems were introduced by means of look-up tables; hence, making it possible to implement a model with nonlinear data. Furthermore, we cast the model into a multiphase problem, in which a different set of RHS functions applies for each phase.

Dividing the time horizon into phases
Given a speed profile v(t) ∈ R ≥0 for t ∈ [0, T ], we assume a sampling time of one second and discretize the driving cycle accordingly with N intervals and the grid set . In each period ∆t j , we consider the first-order differential equation describing the vehicle's longitudinal dynamics (Guzzella & Sciarretta, 2013): where m eq is the equivalent mass of the vehicle, F t (t) is the traction force, F a (t) is the aerodynamic drag force, and F r (t) is the rolling resistance force. We discretize Eq. (1) and approximate dv dt (t) with the explicit Euler scheme that is applied for v(t). We note that this integration scheme is a simple approximation; however, it appears to be appropriate when the grid length is small. By rearranging (1) in terms of F t (t), we identify three possible operating modes for each interval ∆t j : (1) F t > 0, traction; (2) F t < 0, braking; (3) F t = 0, stand-still.
Thereafter, we collect all of the intervals that were subject to the same operating mode in the n p = 3 model phases (see Fig. 2), assuming disjoint subintervals [t j−1 , t j ), j ∈ [N], of the time horizon. Let the function map each time point t to its associated phase.

Control variables
We introduce the integer control variable γ (t) ∈ [n b ] 0 with n b = 7 and the continuous control variable u(t) ∈ [u min , 1], where u min < 0, models the powertrain's mixed-integer nature. The variable γ (t) can help determine whether to operate the HEV in the electric mode (EM is the only power source, the ICE is turned off, and the clutch is disengaged) or in hybrid mode (EM and ICE are simultaneously used to power the vehicle). γ (t) receives a value of 0 whenever the vehicle is required to operate in the electric mode. It can also take on values in the set [n b ] when it operates in the hybrid mode with selected gears Γ ranging from the 1st to the 7th, respectively. The control variable u(t) allows to regulate the torque split between the ICE and the EM in each hybrid operating mode. Specifically, by varying the control u, we identify three different hybrid configurations: (1) if u(t) ∈ [u min , 0): load point shift (LPS). The operating point of the ICE is shifted toward higher loads and the exceeding power recharges the battery; (2) if u(t) = 0: ICE mode. ICE is the only power source and it propels the vehicle; (3) if u(t) ∈ (0, 1]: boost. ICE and EM can cooperate to fulfill the power requirements at the wheels. Fig. 3 illustrates these scenarios during the traction/braking phase. The control u is dropped in the stand-still phase since only two different scenarios are applicable: either the ICE is turned off (γ (t) = 0), or the ICE is turned on (γ (t) ∈ [n b ]) so that the battery can be recharged by means of EM2 that is operated with a fixed value for the speed and torque provided by the ICE.

Differential states
We model the powertrain's dynamics with the fuel mass flow rateṁ f (t), the battery state-of-charge ξ (t), and the ICE cooling water temperature θ(t) as differential states. The latter is needed to account for the higher fuel consumption during the ICE heating-up transient (van Berkel et al., 2014;Lescot et al., 2010).
We express the dependencies of the ODE for t ∈ [0, T ] as: For a detailed description of the smooth functions f(·) we refer to Robuschi (2019). We group the differential states into vectors x(·) and their RHS functions into f(·) as proposed in Problem 1. (3)

Outer convexification
The idea of outer convexification applied to MIOCPs (Sager, 2009) is to reformulate the ODE with affine entering binary controls b(t) ∈ {0, 1} n b +1 in order to get rid of the integer control γ (t). This reformulation opens the door for a canonical relaxation b(t) ∈ [0, 1] n b +1 that leads to an NLP after discretization. The solution of the latter can be used as part of a rounding problem for constructing a solution to the MIOCP. Therefore, we introduce holds. In this way, we reformulate Eq. (3) and obtain the outer convexified dynamics for t ∈ [0, T ]: where (4b) is needed because the definition of the integer controls γ implies mutually exclusive operation modes.

Path and vanishing constraints
The state-of-charge has to fulfill path constraints in order to preserve durability and reliability of the electric buffer. The choice of these limits ξ min , ξ max ∈ [0, 1] is arbitrary and is generally expressed as: The operating points for the ICE and EM torque and the internal current for the battery have to be within a realistic range. This restricts the choices of the continuous and binary controls. We model these restrictions by mode specific lower and upper bounds u lb,i , u ub,i ∈ [u min , 1], i ∈ [n b ] 0 , for u(·) and obtain so-called vanishing constraints: To avoid numerical issues, we relax (5) by replacing zero with the parameter ε = −1e −4 . We chose the above indicator formulation due to its tight relaxation compared with other formulations such as the Big M method, please see Jung et al. (2018) for further details. The next section summarizes the above path and vanishing constraints together with combinatorial constraints in the constraint function g from Problem 1.

Combinatorial constraints
Technical requirements in realistic scenarios imply combinatorial constraints. We already introduced the combinatorial constraint Eq. (4b), which ensures that exactly one mode is active for all time points. This section discusses the further restrictions that includes prefixing, the dwelling time, and mode transition constraints. Before we define these conditions, let us consider the details for the problem discretization because the combinatorial constraints appear more intuitive for the discrete setting.

Problem discretization
We follow a first-discretize-then-optimize approach in the sense that we discretize the ODE with the direct collocation method and Lagrange interpolation polynomials; see Betts (2010) and Biegler (2010) for more details. For the control discretization, we assume piecewise constant controls on the time grid G N .
Hence, we introduce the discrete control variables u ∈ [u min , 1] N , b ∈ {0, 1} (n b +1)×N×np , where b incorporates variables for the phases and the intervals because of our algorithmic idea; see Section 4. Discretized controls can change values only on the grid We apply an equidistant grid G N with a grid length ∆t j = 1s so that we neglect ∆t j in the formulation of combinatorial constraints.

Prefixing constraints
We restrict the set of feasible gear choices to satisfy the minimum and maximum ICE speed. Since the velocity profile is known a priori, it is possible to pre-calculate the allowed gears for each interval j ∈ [N]. Therefore, we exclude some options for all p ∈ [n p ]:

Minimum dwell time constraints
Minimum dwell or min-up time constraints describe the requirement to stay in a certain mode for a chosen time duration after activation. Problem 1 requires min-up times for the electric and hybrid mode that can overlap different phases. Therefore, we introduce sets of the dwell time coupled controls S e the electric and hybrid specific min-up times U S e c , U S h c ∈ N. The constraints are defined for S c ∈ {S e c , S h c }, j = 2, . . . , N − 1, l = j + 1, . . . , min{N, j + U Sc } as: Fig. 4 illustrates an example of the min-up time and mode transition constraints.

Mode transition constraints
By mode transition restrictions we refer to the situation in which the activation of one control b i 1 ,j,p excludes some control indices i 2 from activation in the time step j + l, l ≥ 1: These restrictions are motivated by the dual-clutch gearbox at hand. This can switch one gear up or down per second, which is independent from the active phase. In practice, the driver can use the time during the electric mode to change the gear setting; however, it is limited by one gear shift per second. For all index tuples, (i a , p a ), i a ∈ [n b ], and p a ∈ [n p ] representing the active mode and phase, we introduce the allowed control indices ''neighborhood'' for l = 1, . . . , 5 as: Then, we define the constraint for j = 1 + l, . . . , N as:

Solving combinatorial constrained multiphase mixedinteger control problems
Our algorithmic idea is to apply the theory of CIA decompositions (Sager, 2009;Zeile et al., 2018) for MIOCPs. We provide an intuition of this decomposition before providing a problem specific CIA algorithm.

Intuition of CIA decomposition theory
The difference of Problem 1 compared to the usual MIOCP (Sager, 2009) is the dependency of f i on the current (given) phase p(t). If we set f i (p, ·) to zero outside the corresponding phase p for all i ∈ [n b ] 0 , we obtain an equivalent MIOCP with n p · (n b + 1) modes. In this way, each phase induces n b + 1 additional modes.
Thus, we can equivalently transform a multiphase MIOCP into an MIOCP by requiring for the binary controls: We note that the above constraint corresponds to a simple variable fixing in Problem 3 because p(t j−1 ) is known beforehand from Eq. (2). Discretized MIOCPs result in mixed-integer nonlinear programs (MINLPs), which are generally NP-hard (Belotti et al., 2013), and can only be solved with methods such as branch-and-bound (BnB) under great computational effort. The CIA decomposition suggests to approximate the resulting MINLP by a sequence of subproblems in which each one is less hard to solve than the original problem. In particular, we consider the canonical relaxation of the MINLP by dropping its integrality constraint, which results in the following NLP.
After solving (NLP), we obtain relaxed values a, which we want to round in an appropriate way in order to construct a feasible solution for Problem 1 that exhibits a promising objective value. Minimizing the maximum accumulated control deviation, i.e., results in such a rounding, as proven in Sager et al. (2011). This problem neglects the system's dynamics and can be reformulated into an MILP that takes into account combinatorial constraints. This MILP constitutes the combinatorial integral approximation (CIA) problem, which we now introduce formally in the multiphase setting.
Problem 3 (Multiphase CIA). Let a ∈ [0, 1] np×(n b +1)×N be given together with min-up times U Sc for each coupled control set S c .
The first constraint in (MCIA) ensures that exactly one mode and phase is active for all intervals, whereas the second constraint is a reformulation of Eq. (10). The standard CIA decomposition consists of three steps: (1) Solve (NLP) → x, u, a; (2) Solve (MCIA) → b; (3) Solve (NLP) with fixed b → x, u.
The objective value of the second (NLP) serves as an approximation of the exact discretized MIOCP. Nevertheless, the approach is favorable for application here due to its enormous reduction of the run time (Zeile et al., 2018) and since Theorem 6.3 from Sager (2009) established a convergence result of the approximative to the optimal solution by refining the discretization grid. In the next subsection, we generalize the CIA decomposition and prove that this convergence result still holds.

A tailored CIA decomposition algorithm
The large number of binary variables in ''all-at-once'' rounding in (MCIA) can lead to an infeasible NLP related to terminal, path, or vanishing constraints that cannot be met. Therefore, our idea is to apply more than one rounding step in order to achieve more freedom with respect to achieving a feasible solution. To this end, we propose to solve a sequence of alternating (NLP) and (MCIA) problems where the number of fixed binary variables is gradually increased. The next definition formalizes this idea.
Definition 4 (Binary Subset CIA-NLP Sequence). Let S 1 := {(i, p) | i ∈ [n b ] 0 , p ∈ [n p ]} be the index set of all binary control variables. We denote by S j , j = 2, . . . , n dec , with n dec ≤ n p ·(n b +1) a chosen sequence of binary control index subsets S n dec ⊂ . . . ⊂ S 2 ⊂ S 1 . We define CIA(S j ), j = 1, . . . , n dec −1 as (MCIA), where the binary variables with indices out of S j are optimized with the according input values a. Analogously, NLP(S j ) refers to (NLP), where we relax all b i,p , (i, p) ∈ S j and all b i,p , (i, p) ∈ S 1 \S j are considered to be fixed with values from CIA(S j−1 ). Furthermore, let S n dec := ∅.
We present in Algorithm 1 a tailored version of the CIA decomposition that consists of solving n dec NLPs and n dec − 1 CIA problems with gradually decreased number of free binary controls as given in Definition 4. We illustrate this algorithm in Fig. 5, where we also give reference to the used solver software.
Assumption 5. The values a obtained from solving NLP(S i ), i = 1, . . . , n dec − 1 are almost combinatorial in the sense that each CIA(S i ), i = 1, . . . , n dec − 1 has an objective ζ i that is bounded by: [N] ∆t j , with constants C i > 0 so that ζ i vanishes with ∆t → 0.  5. The concept of Algorithm 1 is to solve MIOCPs and thereby Problem 1. After outer convexification and discretization, we obtain an MINLP, which we relax by removing the integrality constraint. The algorithm involves solving a sequence of alternating CIA and NLP problems. The number of free variables in these problems are represented by S j according to the binary subset CIA-NLP sequence in Definition 4 and is gradually reduced until all variables a are fixed in (NLP(S n dec )). The objective value of the latter problem serves as an approximation to (MINLP); however, it can be calculated much faster by this algorithm. We achieve local minima of the NLPs by using IPOPT (Wächter & Biegler, 2006), while we apply the BnB solver of pycombina (Bürger et al., 2020) to obtain the global optimal solutions of the CIA problems.
Algorithm 1 CIA decomposition for solving Problem 1 Input: Discretized Problem 1 with time grid G N .
Output: (Local) optimal variables x * , u * , b * , with objective value Remark 6. Assumption 5 is not restrictive for the case without dwell time or mode transition constraints, see Sager et al. (2011). Since we include these constraints in (MCIA), this assumption is critical. Therefore, we argue that it is advantageous in technical practice to adhere to mode transition and dwell time rules and the proposed model transfers these technical relationships to the structure of a.
We adopt Theorem 6.3 from Sager (2009) and provide the theoretical justification for Algorithm 1. When we apply Algorithm 1 with a vanishing grid length and under Assumption 5, the rounding gap also vanishes, i.e.: proof. Note first that the functions f i are Lipschitz continuous since they are assumed to be smooth. We denote by C L > 0 their maximal Lipschitz constant. From Theorem 4.1 in Sager (2009), the differential state approximation error for a single rounding step is bounded in the sense of: where ζ is the objective value of the CIA problem, x rel is the differential state trajectory received by using the relaxed controls a and x bin the one using binary controls b. We transform this result to our setting in which x S i denotes the trajectory based on solving NLP(S i ), i ∈ [n dec ]. Thus, we obtain: In the penultimate approximation step, we took advantage of Assumption 5. Finally, the claim follows by the continuity of the objective function J. ■ Theorem 7 states that Problem 1 can be approximated arbitrarily well by Algorithm 1 using grid refinement. Compared to the original CIA decomposition, the rounding error is slightly larger due to (possibly) several rounding steps. In return, we obtain the desired effect in relation to constructing a feasible solution.

Numerical results
We perform a case study of Problem 1 and Algorithm 1 applied on the Worldwide Harmonized Light Vehicles Test Procedure (WLTP) driving cycle, which represents a real and challenging optimization problem due to a long time horizon and the frequent activation of vanishing constraints. Before that we show exemplarily how the relaxed and binary controls behave as part of the CIA problem in connection with the combinatorial constraints.

Used hardware and software
All computations were conducted on a Dell XPS15 desktop PC with an Intel Core i7-6700HQ CPU and 16 GB RAM running Ubuntu 16.04. To parse the NLPs we used CasADi 3.4.5 (Andersson et al., 2019) within the Python 2.7 environment, while the solution is provided by the sparse NLP solver IPOPT 3.12.3 (Wächter & Biegler, 2006), running the linear solver MA97 from HSL (2018). We applied pycombina 2 (Buerger et al., 2019) to solve the CIA problems. pycombina is an open-source software package designed for such MILP problems since its BnB algorithm exploits the specific min-max structure of Eq. (10) and has been proven to be up to three orders of magnitude faster than standard MILP solvers (Bürger et al., 2020).

Exemplary CIA rounding step
We illustrate the functionality of the CIA problem with the driving cycle from Fig. 2. After solving NLP(S 1 ), we obtain the relaxed binary control values a, which we depict in Fig. 6 with In the next step, we solve CIA(S 1 ) so that we obtain binary values b that fulfill all combinatorial constraints, as depicted with the gray lines in Fig. 6. The binary values approximate the relaxed ones quite well and the few differences are mainly due to the min-up time. For instance, this can occur during the second activation of the electric mode.

Case study with the WLTP driving cycle
The velocity profile of the WLTP driving cycle is given in Fig. 7. We solve Problem 1 applied to this driving cycle with Algorithm 1 and n dec = 2, 3 decomposition steps. Moreover, we set u min = −1 as a lower bound for the torque split control. The case n dec = 2 refers to the standard CIA decomposition, where we used a predefined gearshift profile obtained by applying a heuristic algorithm 3 and we therefore optimize only the binary choice between the electric and hybrid mode. For the algorithmic case n dec = 3 we first optimize the gear and electric mode choices, i.e., S 1 = {(i, p) | i ∈ [n b ] 0 , p ∈ [n p ]}. Afterwards, we fix Table 1 Comparison of the normalized objective function J and the run time (CPU) for the solution of Problem 1 obtained by dynamic programming and the NLPs from Algorithm 1, with either optimized or predefined gear choices and varying dwell time constraints for the WLTP driving cycle. i.e., we set S 2 = {(i, p) | p ∈ [n p ], i = 0} to achieve optimization between the electric and hybrid mode and, complementary, we fix the gearshift pattern found in the previous step and use it as an exogenous variable in NLP(S 2 ). To compare the proposed algorithm with a method constructing a global optimal solution, we solved Problem 1 also with a backward dynamic programming approach (Wang et al., 2015); however, we skip detailed dwell-time scenarios since this is beyond the scope of this research. We collect in Table 1 the values of the normalized total fuel consumption of the three approaches with varied minimum dwell times from one to five seconds. The objective value increases progressively from the first to the last NLP required to solve the problem, as expected from Theorem 7. In addition, the fuel consumption increases with increasing dwell times, which involves a substantially decreased number of switches (from 54 to 41); thus, providing a better driveability. When comparing the predefined and optimized gearshift scenarios, we obtain for the latter savings of 13.45% and 12.05% for the fuel consumption when the dwell time constraint (DTC) is set to 5 s and 1 s, respectively. Note that this comes at the expense of an increased run time, since a total of 4325 s instead of 1385 s is required to optimize the gearshift when the DTC is set to 1 s. We left out the run times of the CIAs, because the tailored BnB feature of pycombina runs only for a few seconds. Dynamic programming is meant to provide globally optimal solutions. However, in a practical implementation, the is usually different from the values in the state space tabulation. This effect, possibly increased by using different integration schemes, is also the reason why in our implementation the objective function value of the dynamic programming solution has a higher objective function value than the one found by our direct optimization approach. Nevertheless, we see the similarity of the found solutions as an indication for the quality of our new approach. Fig. 7 presents the evolution of the state and control trajectories for both the predefined and optimized gearshift scenarios for Algorithm 1. It is worth noting how the variation of the state profiles between NLP(S 1 ) and NLP(S 2 ) (6th plot in Fig. 7) is more pronounced by considering the optimization of the gearshift. This is mainly due to the enforcement of combinatorial constraints after CIA(S 1 ). The difference from NLP(S 2 ) to NLP(S 3 ) is marginal since a negligible rearrangement of the switching is needed in CIA(S 2 ). Finally, it can be observed in Fig. 8 how the optimal Fig. 7. Top: velocity profile for the WLTP driving cycle. Plots 2-9: profiles of the state-of-charge for the battery ξ , the ICE cooling temperature θ , the torque split u and the gear choice Γ of the solutions obtained with the predefined and optimized gearshift approaches of Algorithm 1 and with a minimum dwell time equal to five seconds. We imposed the constraints ξ (t f ) = 0.6 and 0.4 ≤ ξ (t) ≤ 0.7. control solutions translate into different operating points of the electric drive (ED, electric motor plus inverter) and ICE. The optimized gearshift entails a higher gear selection in comparison to the predefined gears; thus, allowing the ICE to operate with a greater torque and a lower speed within a higher efficiency region. On the other hand, the difference among the operating points of the ED is negligible, since the EM is directly coupled through a constant transmission ratio to the final drive shaft.

Conclusions
This study presents a tailored CIA decomposition to address the solution of multiphase MIOCPs applied to the EMS for an HEV. Our algorithm is able to cope with a general problem class including multiphase, vanishing, state and combinatorial constraints. We proved that the algorithm constructs under certain assumptions a near-optimal solution since the optimal solution can be approximated arbitrarily well by refining the discretization grid. This study showcased the effectiveness of our approach for the realistic WLTP driving cycle. It also demonstrated accurate solutions with reasonable run times, which makes it possible to benchmark causal controllers or to objectively compare different powertrain architectures. The findings from this research can be beneficial to researchers and professionals that work in the field of hybrid electric vehicles. Future work may address the online application of the proposed algorithm, more sophisticated approaches to include terminal or path constraints and other gearbox settings. Dr. Robuschi's research activity was funded and carried out in cooperation with the GT R&D Department of Ferrari S.p.A., focusing on the optimal energy management of hybrid electric vehicles. His research interests include optimization algorithms for hybrid electric vehicles and high-fidelity simulation methods for non-smooth dynamical systems. His Ph.D. thesis was awarded the Cum Laude mention. Clemens Zeile studied mathematics at the Universities of Hamburg, Lund, Berkeley and Göttingen, where he also obtained his M.Sc. degree in 2015. He is currently a Ph.D. degree student at the Otto-von-Guericke Universität Magdeburg with research focus on the theory and applications of mixed-integer optimal control. Sebastian Sager studied mathematics in Heidelberg, where he also obtained his Ph.D. in 2006 and his habilitation in 2012. Since 2012 he is full professor for algorithmic optimization at the Otto-von-Guericke Universität Magdeburg. The focus of his work is on the development of optimization algorithms for problems that combine the properties integrality, nonlinearity, time-dependence, and uncertainty. He uses them in the estimation, control, experimental design, and analysis of dynamic systems and for machine learning. In recent years he has been addressing particular challenges that arise in clinical decision support. He received an ERC Consolidator Grant in 2015 and is spokesperson of the DFG 2297 Research Training Group ''Mathematical Complexity Reduction''.