Parallel Hybrid Electric Vehicle Modelling and Model Predictive Control

: This paper presents the modelling and calculations for a hybrid electric vehicle (HEV) in parallel conﬁguration, including a main electrical driving motor (EM), an internal combustion engine (ICE), and a starter/generator motor. The modelling equations of the HEV include vehicle acceleration and jerk, so that simulations can investigate the vehicle drivability and comfortability with different control parameters. A model predictive control (MPC) scheme with softened constraints for this HEV is developed. The new MPC with softened constraints shows its superiority over the MPC with hard constraints as it provides a faster setpoint tracking and smoother clutch engagement. The conversion of some hard constraints into softened constraints can improve the MPC stability and robustness. The MPC with softened constraints can maintain the system stability, while the MPC with hard constraints becomes unstable if some input constraints lead to the violation of output constraints.


Introduction
Controllers for HEVs powertrains and speeds can be included model-free or modelbased. Model-free controllers are mostly used with heuristic, fuzzy, neuro, AI, or human virtual and augmented reality. The use of model-free methods will be presented in the next part of this study. Model-based controllers can be used with a conventional adaptive PID, H 2 , H ∞ , or sliding mode. However, all conventional control methods cannot include the real-time dynamic constraints of the vehicle physical limits, the surrounding obstacles, and the environment (road and weather) conditions. Therefore, a MPC with horizon state and open loop control prediction subject to dynamic constraints are mainly used to control as real-time the HEV speeds and torques. Due to the limit size of this paper, we have reviewed some of the most recent research of MPC applications for HEVs. In this paper, vehicle dynamic formulas and calculations are referred to the reference [1].
A recent modelling and control of the dual clutch transmission for HEVs are presented in [2], in which a new controller was designed for synchronizing the dual clutch transmission (DCT) with higher performances and lower fuel consumption. Another MPC for an autonomous driving vehicle is developed in [3], in which the MPC was used to drive the HEV to exactly track the given feasible trajectories. Moreover, a controller for a hybrid dual clutch transmission powertrain for a HEV is introduced in [4], in which the ICE and EM were driven by a DCT powertrain. A MPC for a HEV with linear parameters and a varying model is presented in [5], in which the MPC controller was designed to improve the fuel economy of the power-split HEV.
A MPC for a HEV is necessary not only to control the torque and speed, but also to control the gas emission and improve the fuel economy. The authors of [6] developed a

Modelling of the Parallel HEV
In parallel hybrid electric vehicles, both combustion engine and main electric motor are installed in parallel and work in independent configurations. The vehicle can run and switch in four driving modes: pure main electric motor (EV1) at a low speed and/or low load; pure combustion engine (EV2) at a high speed and/or high load; both main electric motor and combustion engine (EV3) at a very high and/or very high load; and a combination of the main electric motor, generator motor, and combustion engine (EV4) at extreme high speed and/or load. In 2021, Hyundai introduced a new version of Sonata Hybrid, the middle size family passenger HEV that combines updated technologies for a typical parallel electric vehicle, as shown in Figure 1.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 16 switch in four driving modes: pure main electric motor (EV1) at a low speed and/or low load; pure combustion engine (EV2) at a high speed and/or high load; both main electric motor and combustion engine (EV3) at a very high and/or very high load; and a combination of the main electric motor, generator motor, and combustion engine (EV4) at extreme high speed and/or load. In 2021, Hyundai introduced a new version of Sonata Hybrid, the middle size family passenger HEV that combines updated technologies for a typical parallel electric vehicle, as shown in Figure 1. This parallel HEV consists of one internal combustion engine (ICE) with four cylinders with multiple point injection, a volume of 2.4 litters, a maximum power of 156 kW at 6000 rpm, and a peak torque of 265 Nm; one electric motor starter (EM2) with a maximum power of 8 kW and a maximum torque of 43 Nm; the main electric motor (EM1) with a maximum power of 35 kW and a maximum torque of 205 Nm; the battery HEV Li-ion with a capacity of 6.1 Ah; a transmission gearbox with fully automated 6 speeds; and a friction clutch engagement. The vehicle curb weight is 1569 kg. This Sonata Hybrid vehicle is used to simulate our system modelling and test the new MPC scheme with softened constraints.
The schematic architecture of the 2021 Hyundai Sonata Hybrid in Figure 1 can be modelled with a simple drivetrain and is shown in Figure 2. The first part of this mechanical structure consists of an internal combustion engine (ICE) and the electric starter/generator motor (EM2) can be grouped into one inertia , including the left clutch disk, the sharp 1, EM2, and ICE. is modelled as one rigid inertia. and are the torques on the ICE and EM2, respectively. and are the angular position and the velocity of the sharp 1, respectively. Similarly, 2 J is modelled as the lumped rigid inertia of the main electric motor EM1 and the right clutch disk. and are the angular position and velocity of the sharp 2, respectively. The third powertrain part connecting the gearbox and the vehicle's driven wheels can be modelled by a gear ratio via a torsional spring and damper with , , and as the stiffness, damping, and acceleration coefficient, respectively, of which the acceleration has not been studied before. The third part, the lumped inertia , consists of the rest of the vehicle, including the gearbox, differential gear, shaft This parallel HEV consists of one internal combustion engine (ICE) with four cylinders with multiple point injection, a volume of 2.4 litters, a maximum power of 156 kW at 6000 rpm, and a peak torque of 265 Nm; one electric motor starter (EM2) with a maximum power of 8 kW and a maximum torque of 43 Nm; the main electric motor (EM1) with a maximum power of 35 kW and a maximum torque of 205 Nm; the battery HEV Li-ion with a capacity of 6.1 Ah; a transmission gearbox with fully automated 6 speeds; and a friction clutch engagement. The vehicle curb weight is 1569 kg. This Sonata Hybrid vehicle is used to simulate our system modelling and test the new MPC scheme with softened constraints.
The schematic architecture of the 2021 Hyundai Sonata Hybrid in Figure 1 can be modelled with a simple drivetrain and is shown in Figure 2. The first part of this mechanical structure consists of an internal combustion engine (ICE) and the electric starter/generator motor (EM2) can be grouped into one inertia J 1 , including the left clutch disk, the sharp 1, EM2, and ICE. J 1 is modelled as one rigid inertia. M ICE and M EV1 are the torques on the ICE and EM2, respectively. θ 1 and ω 1 are the angular position and the velocity of the sharp 1, respectively. Similarly, J 2 is modelled as the lumped rigid inertia of the main electric motor EM1 and the right clutch disk. θ 2 and ω 2 are the angular position and velocity of the sharp 2, respectively. The third powertrain part connecting the gearbox and the vehicle's driven wheels can be modelled by a gear ratio i via a torsional spring and damper with k θ , k β , and k α as the stiffness, damping, and acceleration coefficient, respectively, of which the acceleration has not been studied before. The third part, the lumped inertia J 3 , consists of the rest of the vehicle, including the gearbox, differential gear, shaft 3, and the driven wheels. θ 3 and ω 3 are the angular position and velocity of shaft 3, respectively. rr is the rolling radius of the vehicle's wheels.
3, and the driven wheels. and are the angular position and velocity of shaft 3, respectively. is the rolling radius of the vehicle's wheels. In this paper, all vehicle dynamic formulas and constraints were taken from the technical book in [1]. The vehicle resistance torque is the approximation of the air density , air drag coefficient , the vehicle crossing area , the wheel rolling radius , vehicle friction resistant coefficient , natural gravity , vehicle mass , and the polynomial coefficients of , and The vehicle rolling resistance torque can be calculated as: In Equation (1), the additional road conditions, such as the road dynamics, the road increase, and other environment conditions, can be added as disturbances that lead to some reduction of or increase in the vehicle rolling resistance torque. Changes of vehicle velocity, depending on the road conditions as well as the vehicle dynamic constraints between the vehicle speed and vehicle steering wheel, are referred to in [1].
At a low speed of less than 40 km/h, the clutch is open, and only the main electric motor EM1 propels the HEV. The contribution of some other exponential coefficients is small and can be ignored. The vehicle rolling resistance torque at a low speed can be simplified as: where is the initial resistance constant of air drag and rolling friction. is a linear coefficient that depends on the gear ratio.
On the first part, the torque applied is: This torque can be calculated as: where is the torque from ICE; is the torque from motor ME2; and C M is the torque from the clutch.
When the clutch is locked, the clutch torque is the maximum static clutch friction: where is the clutch radius; is the normal force; and µS is the clutch friction coefficient. When the clutch is in the transitional period, < , the clutch torque is: where is the clutch slipping coefficient. In this paper, all vehicle dynamic formulas and constraints were taken from the technical book in [1]. The vehicle resistance torque is the approximation of the air density ρ, air drag coefficient c w , the vehicle crossing area A, the wheel rolling radius rr, vehicle friction resistant coefficient f r , natural gravity g, vehicle mass m, and the polynomial coefficients of a 0 , a 1 and a 2 The vehicle rolling resistance torque M v can be calculated as: In Equation (1), the additional road conditions, such as the road dynamics, the road increase, and other environment conditions, can be added as disturbances that lead to some reduction of or increase in the vehicle rolling resistance torque. Changes of vehicle velocity, depending on the road conditions as well as the vehicle dynamic constraints between the vehicle speed and vehicle steering wheel, are referred to in [1].
At a low speed of less than 40 km/h, the clutch is open, and only the main electric motor EM1 propels the HEV. The contribution of some other exponential coefficients is small and can be ignored. The vehicle rolling resistance torque at a low speed can be simplified as: where M v0 is the initial resistance constant of air drag and rolling friction. k v is a linear coefficient that depends on the gear ratio. On the first part, the torque applied is: This torque can be calculated as: where M ICE is the torque from ICE; M M2 is the torque from motor ME2; and M C is the torque from the clutch. When the clutch is locked, the clutch torque M C is the maximum static clutch friction: where r C is the clutch radius; F NC is the normal force; and µ S is the clutch friction coefficient. When the clutch is in the transitional period, M C < M Static f max , the clutch torque is: where µ K is the clutch slipping coefficient. On the second part, the torque applied on the main motor ME1 is: The sum of inertias is calculated as: The torque changing is calculated as: .
The balance of torque M 2o is calculated as: where η is the transmission efficiency of the gearbox and the differential gear. The above torque equations can be transformed to the following dynamic equations: .
The angular acceleration of the shaft 1 is calculated as: .
where k β1 is the shaft 1 friction coefficient. The angular acceleration of the shaft 2 is calculated as: .
where k β2 is the shaft 2 friction coefficient. Finally, the angular acceleration of the shaft 3 is calculated as: .
where k β3 is the shaft 3 friction coefficient. The jerk on the drivetrain is calculated as: . .
The torque generated on the main motor is calculated as: where M DC_MOTOR is the main motor torque; k T is the motor constant, k T = M Torque I Current (Nm/A); k E is the electromotive force (EMF) constant, k E = k T ; R I is the resistance; V is the voltage supply; and ω is the angular velocity. Now we proceed and transform all the above equations into a first order linear system as: . . . . . . .
If we put the space vector as for the torque on the combustion engine (ICE), the input voltage for motor EM1 and EM2, torque on clutch, and the initial air-drag load, a linear space state of the vehicle dynamics system can be expressed as: The linear first order state space model in Equation (24) can be used to form the MPC algorithms in the next part. The system in Equation (24)   When the HEV runs in a low speed of less than 40km/h, only the main motor EM1 is activated. Considering the inputs of M ICE = 0, V 1 = 0, and M C = 0, and the state variables of θ 1 = 0 and ω 1 = 0, then the above linear system can be simplified as: where the states The output T Torque3 is the unmeasured torque at shaft 3. In this When the HEV runs in a high speed greater than 40 km/h, the starter motor EM2 activates the combustion engine ICE while the friction clutch is still open, the state equations of the first part can be written as: . .
where ς is the additional coefficient for starting motor EM2 as a compensation load for the starting period. The linear state space system in the first part is: .
x e = A e x e + B e u e y e = C e x e + D e u e , with where x e = [ θ 1 ω 1 ] , u e = V 1 M ICE , and y e = ω 1 T Torque1 . The output T Torque1 is the unmeasured torque at shaft 1.

Model Predictive Control for the HEV
A MPC is an open loop, infinite horizon prediction, and optimization subject to dynamic constraints. The continuous first order linear space state equation in Equation (24) can be discretized into time intervals with a discrete k and k + 1 = k + ∆t (∆t is the computer scanning speed or the time sampling interval). Now, the continuous time form in Equation (24) can be discretized into: which is subject to the states, inputs, outputs, and the input increase constraints A MPC calculates the open loop input and output prediction horizon. For the calculation simplicity, we assumed the input prediction length is always equal to the output prediction length, or N u = N y . The objective function of the MPC for the HEV is: By substituting subject to the linear matrices inequality (LMI), GU ≤ W + Ex(t), where the column vector U u k , . . . , u k+N−1 ∈ R s , s mN u is the optimization vector, H = H > 0, and H, F, Y, G, W, E are obtained from Q, R and in (9) as only the optimizer vector U is needed, the term involving Y is usually removed from (10). The optimization problem (10) is a quadratic program (QP). The MPC optimizer will calculate the optimal input vector U ∆u k , . . . , ∆u k+N u−1 subject to the hard constraints of the inputs, u k ∈ U , and u k+i ∈ [umax min ]; of the outputs y k ∈ Y, and y k+i|k ∈ [ymax min ]; and of the input increments ∆u k+i ∈ [∆umax min ]. But only the first input increment, ∆u k , is taken into the implementation. Then, the optimizer will update the outputs and states variables with the new update input and repeat the calculation for the next time interval. Therefore, the MPC is also called as the receding time horizon control. A diagram control system for this NMPC is shown in Figure 3. . But only the first input increment, Δ , is taken into the implementation. Then, the optimizer will update the outputs and states variables with the new update input and repeat the calculation for the next time interval. Therefore, the MPC is also called as the receding time horizon control. A diagram control system for this NMPC is shown in Figure 3. The MPC scheme for the HEV in Figure 3 calculates the real-time optimal control action, Δ , and feeds into the vehicle dynamic equations and updates the current states, inputs, and outputs. The updated states, inputs, and outputs will feedback and compare to the reference desired trajectory data for generating the next optimal control action, Δ , in the next interval.
When the system is non-linear and has a general derivative nonlinear form, it is calculated as: where x is the state variables and u is the inputs. The non-linear equation in (33) can be approximated in a Taylor series at referenced positions of ( , ) for = ( , ), so that: where . and . are the Jacobian function calculating approximation of and , respectively, moving around the referenced positions ( , ).
The linearized system in Equation (35) can be used as the linear system in Equation (24) for the MPC calculation. However, the MPC real-time optimal control action Δ | Figure 3. Diagram of the MPC system.
The MPC scheme for the HEV in Figure 3 calculates the real-time optimal control action, ∆u k , and feeds into the vehicle dynamic equations and updates the current states, inputs, and outputs. The updated states, inputs, and outputs will feedback and compare to the reference desired trajectory data for generating the next optimal control action, ∆u k , in the next interval.
When the system is non-linear and has a general derivative nonlinear form, it is calculated as: .
where x is the state variables and u is the inputs. The non-linear equation in (33) can be approximated in a Taylor series at referenced positions of (x r , u r ) for .
X r = f (x r , u r ), so that: .
where f x.r and f r.x are the Jacobian function calculating approximation of x and u, respectively, moving around the referenced positions (x r , u r ). Substituting Equation (34) for .
X r = f (x r , u r ), we can obtain an approximation linear form in continuous time (t): .
The linearized system in Equation (35) can be used as the linear system in Equation (24) for the MPC calculation. However, the MPC real-time optimal control action ∆u k+i|k must Appl. Sci. 2021, 11, 10668 9 of 18 be fed into the original non-linear system in Equation (33) for the updated states, outputs, and inputs.

The MPC with Softened Constraints for the HEV
The conventional MPC objective function in Equation (31) subject to the constraints in Equation (30) regarding states, outputs, inputs, and input increase may deal with so many hard constraints. The MPC optimizer may not find out a solution that satisfies all constraints. Thus, we considered to widen the MPC feasibility by converting some possible hard constraints from Equation (30) into softened constraints to increase the possibility of finding a solution. The new MPC scheme subject to the softened constraints has the following form: where µ is assigned as big values as a weighting factor (µ > 0), and ε i is the constraints penalty terms (ε i ≥ 0) added into the MPC objective function. X and z i are the corresponding matrix of the hard constraints.
The new items in Equation (37) are softened constraints selected from hard constraints in u k ∈ U , and u k ∈ U , ∆u k+i ∈ [∆umax min ], for i = 0, 1, . . . , N u − 1, y k ∈ Y, and y k+i|k ∈ [ymax min ], for i = 0, 1, . . . , N y − 1, ∆u k = u k − u k−1 ∈ ∆U , and ∆u k+i = 0, for i ≥ N u , x k|k = x(k), x k+i+1|k = A(k)x k+i|k + B(k)u k+i , u k+i|k = u k+i−1|k + ∆u k+i|k , y k+i|k = C(k)x k+i|k , where, ε i (k) = ε y ; ε u , y k+i|k y max min , and u k+i|k u max min ; and Λ = Λ ≥ 0 is the additional penalty matrix (generally Λ > 0 and assign to small values). In this new MPC scheme, the penalty term of the softened constraints subject to G S U S ≤ W S + E S x(k), where U S u k , u k+1 , · · · , u k+N p −1 , ε k , ε k+1 , · · · , ε k+N p is the new optimization input vector; H S = H 0 0 M ; F S = F µ ; and the matrices for inequality constraints H, F, G, W, and E are obtained from Equation (38), To illustrate the ability of this controller, we test the two MPC schemes in Equations (31) and (36) by the following simple example as considering the non-linear system shown below: .
It is assumed that the system in Equation (39) Figure 4 shows the performance of two NMPC schemes: this initial state position, x 0 , does not lead to any violation of states and input (x min = −1 −1 and −2 ≤ u ≤ 2). In this x 0 , the solutions of the two control schemes are always available. We can see that the NMPC with a softened state approaches the asymptotic point faster than the hard constraints. It means that, if we somehow loosen some of the constraints, the optimizer can generate easier optimal inputs and the system will be more stable.
hard constraints. It means that, if we somehow loosen some of the constraints, the optimizer can generate easier optimal inputs and the system will be more stable.
It is interesting to see in Figure    If we select x 0 = −0.9 −0.8 , this initial condition will lead to the violations of the state and the input constraints as x 1 min = −1.0441 and u max = 2.2303. These violations will make the RMPC with hard constraints infeasible. Meanwhile, the RMPC scheme with softened constraints is still running well and it is still easy to find optimal input solutions, as shown in Figure 5. After a very short transitional period, the fully constrained solution is returned, or there is no more constrained violation.
If we select = −0.9 −0.8 , this initial condition will lead to the violations of the state and the input constraints as = −1.0441 and = 2.2303. These violations will make the RMPC with hard constraints infeasible. Meanwhile, the RMPC scheme with softened constraints is still running well and it is still easy to find optimal input solutions, as shown in Figure 5. After a very short transitional period, the fully constrained solution is returned, or there is no more constrained violation. The new MPC scheme with softened constraints for the HEV will be further analyzed and simulated in the next section.

The MPC for the HEV in Pure Electrical Drive
The main motor ME1 was used to run the HEV at a low speed. In this mode, the clutch is open. ICE and ME2 are off. We run the MPC in this mode with the discrete time interval of 0.05 s. ME1 has a maximum power of 35 kW, a maximum torque of 205 Nm, rigidity torque = 1158, inertia = 1, constants = = 10, inertia = 2, gear ratio = 2.34, damping = 0.5 and = 12, and resistance = 5.
Some softened constraints were converted, as input constraints for the DC voltage applied for the vehicle is | | ≤ 300 V, ( ) = +/−5 V. The output softened constraints are also set on the shaft with the shear strength for carbon steel of = 25 MPa or N/mm 2 . The output torque on the shaft 2 is constrained as | | = , with the diameter = 0.05 m. Then, the torque softened constraint on shaft 2 is | | = 455 Nm.
The MPC parameters were set up with the predictive horizon of = = = 5; the weighting matrices are set at = 1 0 0 1 and = [1]. The MPC performance with softened constraints is shown in Figure 6. The new MPC scheme with softened constraints for the HEV will be further analyzed and simulated in the next section.

The MPC for the HEV in Pure Electrical Drive
The main motor ME1 was used to run the HEV at a low speed. In this mode, the clutch is open. ICE and ME2 are off. We run the MPC in this mode with the discrete time interval of 0.05 s. ME1 has a maximum power of 35 kW, a maximum torque of 205 Nm, rigidity torque k θ = 1158, inertia J 2 = 1, constants k E2 = k T2 = 10, inertia J 3 = 2, gear ratio i = 2.34, damping k β2 = 0.5 and k β3 = 12, and resistance R I2 = 5.
Some softened constraints were converted, as input constraints for the DC voltage applied for the vehicle is |V 2 | ≤ 300 V, u(t) = +/ − 5 V. The output softened constraints are also set on the shaft with the shear strength for carbon steel of τ = 25 MPa or N/mm 2 . The output torque on the shaft 2 is constrained as |T| = τπ d 3 16 , with the diameter d = 0.05 m. Then, the torque softened constraint on shaft 2 is |T 2 | = 455 Nm.
The MPC parameters were set up with the predictive horizon of N u = N y = N p = 5; the weighting matrices are set at Q = 1 0 0 1 and R = [1]. The MPC performance with softened constraints is shown in Figure 6. It is worth noting that the weighting matrix for output Q and input R can be varied according to the desired variation on outputs or inputs. If we want to limit the errors or keep the output variation in a small value, we have to pay for more input energy or increase the input variation. With this aim in mind, we increased Q and reduced R. It means that any small variation in output will lead to a big penalty amount adding to the MPC objective function. Figure 7 shows the MPC for the HEV performance with Q = 100 and R = 1. Appl. Sci. 2021, 11, x FOR PEER REVIEW 12 of 17 It is worth noting that the weighting matrix for output Q and input R can be varied according to the desired variation on outputs or inputs. If we want to limit the errors or keep the output variation in a small value, we have to pay for more input energy or increase the input variation. With this aim in mind, we increased Q and reduced R. It means that any small variation in output will lead to a big penalty amount adding to the MPC objective function. Figure 7 shows the MPC for the HEV performance with Q = 100 and R = 1.  As shown in Figures 6 and 7, we set up softened constraints on the input voltage of | | ≤ 300 V. The MPC allows some input voltage violation at the starting time to ensure the controller stability and feasibility. Then, after a very short transitional period, the solution is returned without constraint violation. In these cases, the MPC with hard constraints becomes infeasible and unstable.

The MPC for the HEV in High Speed with ICE
When the HEV runs at a high speed, the starter/generator ME2 starts the ICE. Depending on the required output torque, the ICE alone, or the ICE and ME1, or all ICE, As shown in Figures 6 and 7, we set up softened constraints on the input voltage of |V 2 | ≤ 300 V. The MPC allows some input voltage violation at the starting time to ensure the controller stability and feasibility. Then, after a very short transitional period, the solution is returned without constraint violation. In these cases, the MPC with hard constraints becomes infeasible and unstable.

The MPC for the HEV in High Speed with ICE
When the HEV runs at a high speed, the starter/generator ME2 starts the ICE. Depending on the required output torque, the ICE alone, or the ICE and ME1, or all ICE, ME1 and ME2 will be running and together providing torque.
In this mode, we assumed that the vehicle runs at ω 3 = 2000 rpm, and the torque of the air drag resistance at this speed of M v0 = 30 Nm. The parameters of the starter motor EM2 are as constants k E2 = k T2 = 15, inertia J 1 = 1, damping coefficient k β1 = 0.5, resistance R I1 = 7, compensation ς = 0.5, and a discrete time of 0.05 s.
The softened constraints were imposed on the input voltage constraints for the starter of |V 1 | ≤ 48 V, u(t) = +/ − 5 V/interval, and the output constrained torque on shaft 1 of |T 2 | ≤ 455 Nm.
For the MPC parameters, we selected the predictive horizon length of N u = N y = N p = 5 and the weighting matrices Q = 10 0 0 10 and R = 1 0 0 1 . The MPC performance with starting EM2 is shown in Figure 8. the controller stability and feasibility. Then, after a very short transitional period, the solution is returned without constraint violation. In these cases, the MPC with hard constraints becomes infeasible and unstable.

The MPC for the HEV in High Speed with ICE
When the HEV runs at a high speed, the starter/generator ME2 starts the ICE. Depending on the required output torque, the ICE alone, or the ICE and ME1, or all ICE, ME1 and ME2 will be running and together providing torque.
In this mode, we assumed that the vehicle runs at   In the next simulation, we run the EM2 and the ICE to track the speed desired setpoints and ignite the clutch engagement. It was assumed that the main motor EM1 runs at 1500 rpm and the starter EM2 starts the ICE and is engaged into the system. The clutch engagement must take place at ≥ or = 1.05 * for driving comfortability and a low jerk. The ICE and ME2 must track on the EM1 speed at +5% offset. The full engagement is completed in 2.3 s, as shown in Figure 9.  In the next simulation, we run the EM2 and the ICE to track the speed desired setpoints and ignite the clutch engagement. It was assumed that the main motor EM1 runs at 1500 rpm and the starter EM2 starts the ICE and is engaged into the system. The clutch engagement must take place at ω 1 ≥ ω 2 or ω 1 = 1.05 * ω 2 for driving comfortability and a low jerk. The ICE and ME2 must track on the EM1 speed at +5% offset. The full engagement is completed in 2.3 s, as shown in Figure 9. In Figure 9, we see the ICE and ME2 tracking ME1 on the desired setpoints in 1.9 s. In the normal mode, at speed higher than 40 km/h, the starter ME2 ignites the ICE and is turned as a generator charging to battery. The main motor EM1 now is also turned off. Only the ICE propels the HEV.
In Figure 10, it can be seen that the ME2 is turned off and becomes the generator after igniting the ICE. The main motor EM1 is also turned off, and the ICE alone propels the HEV. The HEV reaches and tracks the desired speed setpoints after 3.5 s. Finally, we compared the performances of the MPC with hard constraints and he MPC with softened constraints. We run the MPC with hard constraints in Equation (31) and the MPC with softened constraints in Equation (36) to track the desired speed setpoints, as shown in Figure 11. Figure 11 shows that the MPC with hard constraints generates smaller inputs and, hence, it needs a longer time to track into the speed setpoint. The MPC with hard constraints reaches the speed setpoint after 4.5 s, while the MPC with softened constraints needs only 3.5 s to fully track into the speed setpoint. In Figure 9, we see the ICE and ME2 tracking ME1 on the desired setpoints in 1.9 s. In the normal mode, at speed higher than 40 km/h, the starter ME2 ignites the ICE and is turned as a generator charging to battery. The main motor EM1 now is also turned off. Only the ICE propels the HEV.
In Figure 10, it can be seen that the ME2 is turned off and becomes the generator after igniting the ICE. The main motor EM1 is also turned off, and the ICE alone propels the HEV. The HEV reaches and tracks the desired speed setpoints after 3.5 s.  In Figure 9, we see the ICE and ME2 tracking ME1 on the desired setpoints in 1.9 s. In the normal mode, at speed higher than 40 km/h, the starter ME2 ignites the ICE and is turned as a generator charging to battery. The main motor EM1 now is also turned off. Only the ICE propels the HEV.
In Figure 10, it can be seen that the ME2 is turned off and becomes the generator after igniting the ICE. The main motor EM1 is also turned off, and the ICE alone propels the HEV. The HEV reaches and tracks the desired speed setpoints after 3.5 s. Finally, we compared the performances of the MPC with hard constraints and he MPC with softened constraints. We run the MPC with hard constraints in Equation (31) and the MPC with softened constraints in Equation (36) to track the desired speed setpoints, as shown in Figure 11. Figure 11 shows that the MPC with hard constraints generates smaller inputs and, hence, it needs a longer time to track into the speed setpoint. The MPC with hard constraints reaches the speed setpoint after 4.5 s, while the MPC with softened constraints needs only 3.5 s to fully track into the speed setpoint. Finally, we compared the performances of the MPC with hard constraints and he MPC with softened constraints. We run the MPC with hard constraints in Equation (31) and the MPC with softened constraints in Equation (36) to track the desired speed setpoints, as shown in Figure 11.

Conclusions
In this study, we have presented the modelling of the HEV and the MPC algorithms for controlling the HEV. In the HEV modelling, we have included the system acceleration and jerk into the equations to investigate and compare the vehicle driving comfortability with different control parameters. The MPC scheme with softened constraints has proved its superiority over the MPC with only hard constraints. The control system now becomes more flexible, stable, and robust against model uncertainties, time variant, and constraint violations. The new MPC scheme can control the HEV with faster clutch engagement and lower jerk reduction. The MPC with softened constraint still stable and robust, while the MPC with only hard constraints becomes unstable and infeasible because of the constraint violations. In a following study, we will investigate the control of the HEV friction clutch for a smooth and fast engagement with high comfortability and low jerk, and apply these algorithms in the real HEV.   Figure 11 shows that the MPC with hard constraints generates smaller inputs and, hence, it needs a longer time to track into the speed setpoint. The MPC with hard constraints reaches the speed setpoint after 4.5 s, while the MPC with softened constraints needs only 3.5 s to fully track into the speed setpoint.

Conclusions
In this study, we have presented the modelling of the HEV and the MPC algorithms for controlling the HEV. In the HEV modelling, we have included the system acceleration and jerk into the equations to investigate and compare the vehicle driving comfortability with different control parameters. The MPC scheme with softened constraints has proved its superiority over the MPC with only hard constraints. The control system now becomes more flexible, stable, and robust against model uncertainties, time variant, and constraint violations. The new MPC scheme can control the HEV with faster clutch engagement and lower jerk reduction. The MPC with softened constraint still stable and robust, while the MPC with only hard constraints becomes unstable and infeasible because of the constraint violations. In a following study, we will investigate the control of the HEV friction clutch for a smooth and fast engagement with high comfortability and low jerk, and apply these algorithms in the real HEV.