A Mixed Logical Dynamical-model Predictive Control (mld-mpc) Energy Management Control Strategy for Plug-in Hybrid Electric Vehicles (phevs)

Plug-in hybrid electric vehicles (PHEVs) can be considered as a hybrid system (HS) which includes the continuous state variable, discrete event, and operation constraint. Thus, a model predictive control (MPC) strategy for PHEVs based on the mixed logical dynamical (MLD) model and short-term vehicle speed prediction is proposed in this paper. Firstly, the mathematical model of the controlled PHEV is setup to evaluate the energy consumption using the linearized models of core power components. Then, based on the recognition of driving intention and the past vehicle speed data, a nonlinear auto-regressive (NAR) neural network structure is designed to predict the vehicle speed for known driving profiles of city buses and the predicted vehicle speed is used to calculate the total required torque. Next, a MLD model is established with appropriate constraints for six possible driving modes. By solving the objective function with the Mixed Integer Linear Programming (MILP) algorithm, the optimal motor torque and the corresponding driving mode sequence within the speed prediction horizon can be obtained. Finally, the proposed energy control strategy shows substantial improvement in fuel economy in the simulation results.


Introduction
Focusing on plug-in hybrid electric vehicles (PHEVs), which are a typical hybrid system (HS) [1], an accurate system model is mandatory to study their control strategy.Commonly used hybrid system modeling methods include hierarchical model [2], Petri net model based on the proposed network diagram [3], the "unified model" of hybrid dynamic systems by Branicky et al. [4] and the Piecewise Affine (PWA) switching model [5].All these models are relatively complex and inconvenient compared to the conventional control technology.In [6], Bemporad and Morari proposed the mixed logical dynamical (MLD) model that includes interdependent physical laws, logic rules and operational constraints, to apply the conventional control methods (e.g., optimization and prediction) to the hybrid system.The MLD model also uses differential equations to avoid the complex conditions in hybrid systems and could be transferred into mixed integer linear programming (MILP) [7] for a solution.
The control strategy for PHEVs aims to find an appropriate power and torque distribution between the combustion engine and the electric motor to improve the fuel economy and emission performance [8].Among various control strategies, the rule-based control method [9] is the most widely used, but it relies on experience to a greater extent, which makes it difficult to adapt to dynamic changes in the driving cycles and to optimize the fuel economy of the vehicle.The optimization-based control strategy [10] includes instantaneous optimization and global optimization.The instantaneous optimization control strategy [11] for an equivalent minimal fuel consumption does not give sufficient consideration to the vehicle state in the future; while for the global optimization [12], complete statistical information on driving cycles must be provided, including road conditions, which makes it hard to compute fast and it lacks adaptability, and eventually cannot be applied directly to real-time control of the vehicle.If the vehicle state information is predicted accurately beforehand, global optimization within this time could be applied to overcome the disadvantages of inadaptability.This could also form the basis of the Model Predictive Control (MPC) method for a hybrid electric vehicle.In [13], a generalized predictive optimal control framework is proposed, and MILP methodology is utilized, which enables the control strategy to get better economy than instantaneous strategies.However, it assumes the traffic pattern and terrain are known in advance without any prediction.
In this paper, the required vehicle torque is calculated based on the predicted vehicle speed.Commonly used vehicle speed prediction methods are the model-based method [14], Kalman filter [15], hidden Markov models [16] and neural network models [17], etc.The non-linear auto-regressive (NAR) model neural network [18,19] has advantages of solving time variant and nonlinear problems, making it suitable for vehicle speed prediction.For training the neural network, most researchers have only used global positioning system (GPS) data from or statistical analysis data from the vehicle driving cycle [20][21][22], but ignored the driver's operation such as accelerator pedal and braking pedal operations which can reflect the future speed variation trends.Therefore, a vehicle speed prediction method using a NAR neutral network based on the combination of past speed data and driving intention data recognized by fuzzy inference is proposed.
Based on the analysis above, this paper focuses on a parallel plug-in hybrid electric city bus, and the powertrain model is set-up firstly.Then, the driving intention is recognized by using fuzzy inference to analyze the pedal operation.A vehicle speed prediction method using a NAR neutral network based on the combination of driving intention and the past speed data is proposed to predict the future short-term vehicle running state and calculate the required vehicle torque.Afterwards, the MLD model is established for the PHEV powertrain.Combined with MILP algorithm and aiming at the minimal equivalent fuel consumption, a mixed logic dynamical-model predictive control (MLD-MPC) strategy based on vehicle speed prediction is proposed.Finally, the simulation results are presented to validate the proposed energy control strategy.

Modeling of Powertrain System
In this part, the vehicle powertrain is modeled.It is used to calculate the total vehicle torque demand, and to evaluate the engine and motor energy consumption quickly for the control strategy, similar torque-speed modelling method could also be found in literature [12,23,24].Parallel powertrain and transmission system of PHEV [25] adopted in this paper are shown in Figure 1.The mathematical models of each part in powertrain are described as follows.The speed and torque models of wheels are described in Equation (1): In Equation ( 1), v is the vehicle velocity in m/s, ρ a is the air density (1.2258 kg/m 3 ), w r is the radius of wheel in m, ω w is the wheel rotating speed in rad/s, W T is the wheel torque in m N ⋅ , vehicle slope angle in rad, d d v / t is the acceleration in m/s 2 .Model of transmission torque and speed is given in Equation (2), where ω in is the required input rotating speed of transmission in rad/s; 0 i is the total transmission ratio; in T is the required input torque of transmission in N m ⋅ ; η is the efficiency from transmission input to tire [26].The model of engine fuel consumption rate is given in Equation (3), where f m  is the engine fuel consumption rate, which is a function of its rotation speed ω e and torque e T : According to the engine fuel consumption data, the relationship between fuel consumption rate f m  and torque under different rotating speed is curve-fitted using the least squares method and is plotted in Figure 2. In Equation ( 4), 0 a and 1 a are the fitting coefficients in linear expression (see Appendix A).In Equation ( 1), v is the vehicle velocity in m/s, ρ a is the air density (1.2258 kg/m 3 ), r w is the radius of wheel in m, ω w is the wheel rotating speed in rad/s, T W is the wheel torque in N • m, A f is the frontal area in m 2 , C d is the drag coefficient, f r represents tire rolling resistance coefficient, m is the vehicle mass in kg, J represents total vehicle rotational inertia in kg • m 2 , α is vehicle slope angle in rad, dv/dt is the acceleration in m/s 2 .Model of transmission torque and speed is given in Equation ( 2), where ω in is the required input rotating speed of transmission in rad/s; i 0 is the total transmission ratio; T in is the required input torque of transmission in N • m; η is the efficiency from transmission input to tire [26].The model of engine fuel consumption rate is given in Equation ( 3), where .m f is the engine fuel consumption rate, which is a function of its rotation speed ω e and torque T e : .
According to the engine fuel consumption data, the relationship between fuel consumption rate .m f and torque under different rotating speed is curve-fitted using the least squares method and is plotted in Figure 2. In Equation (4), a 0 and a 1 are the fitting coefficients in linear expression (see Appendix A).In Equation ( 1), v is the vehicle velocity in m/s, ρ a is the air density (1.2258 kg/m 3 ), w r is the radius of wheel in m, ω w is the wheel rotating speed in rad/s, W T is the wheel torque in  2), where ω in is the required input rotating speed of transmission in rad/s; 0 i is the total transmission ratio; in T is the required input torque of transmission in N m ⋅ ; η is the efficiency from transmission input to tire [26].The model of engine fuel consumption rate is given in Equation ( 3), where f m  is the engine fuel consumption rate, which is a function of its rotation speed ω e and torque e T : According to the engine fuel consumption data, the relationship between fuel consumption rate f m  and torque under different rotating speed is curve-fitted using the least squares method and is plotted in Figure 2. In Equation ( 4), 0 a and 1 a are the fitting coefficients in linear expression (see Appendix A).  . .
Equation ( 5) is the model of equivalent fuel consumption rate of electric motor, where .m m is the rate of fuel consumption which is equivalent to the electric power the motor used, R is the mass calorific value of diesel (43,000 kJ/kg), ω m is the motor speed in rad/s, and T m is the motor torque in N•m, U 0 is the open-circuit voltage in V, I is the battery output current in A. λ is the equivalent fuel factor of the battery power, it is decided by the ratio of electricity and fuel price as can be seen in Equation ( 6), P o is the fuel price (5.6 RMB/L), P e is the electricity price (1 RMB/kWh), ρ is the fuel density (835 kg/L), η c is the efficiency of the charging pile.
According to the battery output current [27] given by Equation (7), the relationship between motor torque and the rate of change of state of charge (SOC) of the battery under a certain motor speed is plotted in Figure 3, and the approximate linear function obtained using least square method is given by Equation (8), where .
x is the rate of change of SOC, b 1 and b 0 are the fitting coefficients (see Appendix A).R i is the equivalent resistance of the battery in Ω, and Qmax is the capacity of the battery in C: .
Energies 2017, 10, 74 4 of 18 λ ( / 3600 / 1000 / η ) / ( / / ρ) According to the battery output current [27] given by Equation (7), the relationship between motor torque and the rate of change of state of charge (SOC) of the battery under a certain motor speed is plotted in Figure 3, and the approximate linear function obtained using least square method is given by Equation (8), where x  is the rate of change of SOC, 1 b and 0 b are the fitting coefficients (see Appendix A). i R is the equivalent resistance of the battery in Ω , and max Q is the capacity of the battery in C:

Driving Intention Identification
Driving intention depends on the driving environment, the vehicle running state and the driver's driving habits, etc.It is a typical empirical model with characteristics of fuzziness.Given the obvious advantages of fuzzy theory in dealing with the empirical models, a fuzzy logic inference system was developed to identify the driving intention.

Driving Intention Identification
Driving intention depends on the driving environment, the vehicle running state and the driver's driving habits, etc.It is a typical empirical model with characteristics of fuzziness.Given the obvious advantages of fuzzy theory in dealing with the empirical models, a fuzzy logic inference system was developed to identify the driving intention.
Driving intention, in general, can be divided into acceleration intention and braking intention.Based on the urgency of accelerating, accelerating intention is divided into slow acceleration, relatively slow acceleration, normal acceleration, relatively rapid acceleration, and rapid acceleration.When no operation on the accelerator pedal it can be considered as braking intention and is divided into general braking and coasting based on whether the braking pedal is operated or not.Based on the urgency of general braking, general braking intention is divided into slow braking, relatively slow braking, normal braking, relatively rapid braking, and rapid braking.
Drivers operate the accelerator pedal and braking pedal directly according to the driving environment and the running state of the vehicle, so as to realize their driving intention.Therefore, the acceleration pedal travel and the braking pedal travel are used as the two main identification parameters.However, the pedal travel cannot completely reflect the urgency of driving intention.Hence the rate of change of acceleration pedal travel and braking pedal travel are also considered to identify the driving intention.
The acceleration intention is identified by fuzzy inference, and its identification parameters are acceleration pedal travel and rate of change of acceleration pedal travel, membership functions of which are shown in Figure 4a,b respectively.The membership functions of the acceleration intention are shown in Figure 4c.The fuzzy rules for the acceleration intention identification are shown in Table 1.
Driving intention, in general, can be divided into acceleration intention and braking intention.Based on the urgency of accelerating, accelerating intention is divided into slow acceleration, relatively slow acceleration, normal acceleration, relatively rapid acceleration, and rapid acceleration.When no operation on the accelerator pedal it can be considered as braking intention and is divided into general braking and coasting based on whether the braking pedal is operated or not.Based on the urgency of general braking, general braking intention is divided into slow braking, relatively slow braking, normal braking, relatively rapid braking, and rapid braking.
Drivers operate the accelerator pedal and braking pedal directly according to the driving environment and the running state of the vehicle, so as to realize their driving intention.Therefore, the acceleration pedal travel and the braking pedal travel are used as the two main identification parameters.However, the pedal travel cannot completely reflect the urgency of driving intention.Hence the rate of change of acceleration pedal travel and braking pedal travel are also considered to identify the driving intention.
The acceleration intention is identified by fuzzy inference, and its identification parameters are acceleration pedal travel and rate of change of acceleration pedal travel, membership functions of which are shown in Figure 4a,b respectively.The membership functions of the acceleration intention are shown in Figure 4c.The fuzzy rules for the acceleration intention identification are shown in Table 1.The braking intention is identified by fuzzy inference, and its identification parameters are braking pedal travel and rate of change of braking pedal travel, membership functions of which are showed in Figure 5a,b respectively.The membership functions of the braking intention are shown in Figure 5c.The fuzzy rules for the braking intention identification are shown in Table 2.The braking intention is identified by fuzzy inference, and its identification parameters are braking pedal travel and rate of change of braking pedal travel, membership functions of which are showed in Figure 5a,b respectively.The membership functions of the braking intention are shown in Figure 5c.The fuzzy rules for the braking intention identification are shown in Table 2.In this paper, the United Kingdom Bus (UKBUS) driving cycle is considered as an example, and the driving intention identification is carried out.Figure 6 shows the speed segment of UKBUS driving cycle.The acceleration intention identification results are in the range of 0-1, the closer that this value is to 1, the more urgent the acceleration intention is.The braking intention identification results are in the range of −1-0, the closer that this value is to −1, the more urgent the braking intention is.Besides, the driving intention is identified as coasting when there are no operations on both accelerator pedal and braking pedal.The resultant values are a continuous time sequence ranging from −1 to 1, based on which the driving intention can be distinguished between acceleration or braking at each and every moment, as well as the urgency of the driving intention.Then the obtained driving intention time sequence, as well as the vehicle speed time sequence, are given as inputs to the NAR neural network.The acceleration calculated through the real velocity is normalized in −1-1 and it is used as the real driving intention intensity.The comparison of the driving intention identified by fuzzy inference and the real driving intention intensity is shown in Figure 7.The mean absolute error (MAE) of the driving intention recognition result is 0.16 when the data ranges from −1 to 1.This fuzzy recognition result with sufficient accuracy could be used as input for the speed prediction neural network to improve the prediction precision.In this paper, the United Kingdom Bus (UKBUS) driving cycle is considered as an example, and the driving intention identification is carried out.Figure 6 shows the speed segment of UKBUS driving cycle.The acceleration intention identification results are in the range of 0-1, the closer that this value is to 1, the more urgent the acceleration intention is.The braking intention identification results are in the range of −1-0, the closer that this value is to −1, the more urgent the braking intention is.Besides, the driving intention is identified as coasting when there are no operations on both accelerator pedal and braking pedal.The resultant values are a continuous time sequence ranging from −1 to 1, based on which the driving intention can be distinguished between acceleration or braking at each and every moment, as well as the urgency of the driving intention.Then the obtained driving intention time sequence, as well as the vehicle speed time sequence, are given as inputs to the NAR neural network.The acceleration calculated through the real velocity is normalized in −1-1 and it is used as the real driving intention intensity.The comparison of the driving intention identified by fuzzy inference and the real driving intention intensity is shown in Figure 7.The mean absolute error (MAE) of the driving intention recognition result is 0.16 when the data ranges from −1 to 1.This fuzzy recognition result with sufficient accuracy could be used as input for the speed prediction neural network to improve the prediction precision.

Short-Term Vehicle Speed Prediction Using Nonlinear Auto-Regressive Neural Network
The NAR neural network is a kind of dynamic neural network, which could add the output memory to the input automatically to calculate the output at the next instance.As the driver's operation can reflect the future speed variation trend, the obtained time series of driving intention are provided as input to the network together with the past vehicle speed data.
Definition of the prediction model is given in Equation ( 9): ( 1) ( ( ), ( 1), ) where is the output at a time t, f is the nonlinear mapping, y d is the output memory order.
The output at the time (t + 1), i.e., ( 1) t + Y depends on the former y d steps of output ( ) t Y , as the predicted values could be fed back to the input, the long-term predictions can be achieved.
The NAR neural network structure, consisting a hidden layer, an output layer, and an output feedback layer, is given in Figure 8.As shown in Figure 8, the hidden and the output layer include 15 and 1 neurons respectively.The order of output memory is 5. 1 w is the connection weight between input and output layer neuron, 1 b is the threshold of the layer, 1 f is the transfer function of the hidden layer, 2 w is the connection weight between hidden layer and output layer, 2 b is the threshold of output layer and 2 f is the transfer function of output layer.
This paper uses a Bayesian algorithm to train the network, choosing Tan-sigmoid transfer function as hidden layer function, and purelin linear function for the output layer.In this paper, a plug-in hybrid electric city bus is studied which is running on pre-determined routes, thus the speed prediction method is for known driving profiles.Taking the UKBUS driving cycle as an example, the driving intention time series corresponding to the conditions are obtained by the aforementioned driving intention identification method.A total of 3288 samples of vehicle speed-intention data, choosing 2788 samples as training samples among which 75% is considered as training data, and the remaining 25% is used as validation data and test data.The remaining 500 samples which didn't participate in the training are used to test the network's predictive ability.

Short-Term Vehicle Speed Prediction Using Nonlinear Auto-Regressive Neural Network
The NAR neural network is a kind of dynamic neural network, which could add the output memory to the input automatically to calculate the output at the next instance.As the driver's operation can reflect the future speed variation trend, the obtained time series of driving intention are provided as input to the network together with the past vehicle speed data.
Definition of the prediction model is given in Equation ( 9): ( 1) ( ( ), ( 1), ) where is the output at a time t, f is the nonlinear mapping, y d is the output memory order.
The output at the time (t + 1), i.e., ( 1) t + Y depends on the former y d steps of output ( ) t Y , as the predicted values could be fed back to the input, the long-term predictions can be achieved.
The NAR neural network structure, consisting a hidden layer, an output layer, and an output feedback layer, is given in Figure 8.As shown in Figure 8, the hidden and the output layer include 15 and 1 neurons respectively.The order of output memory is 5. 1 w is the connection weight between input and output layer neuron, 1 b is the threshold of the layer, 1 f is the transfer function of the hidden layer, 2 w is the connection weight between hidden layer and output layer, 2 b is the threshold of output layer and 2 f is the transfer function of output layer.
This paper uses a Bayesian algorithm to train the network, choosing Tan-sigmoid transfer function as hidden layer function, and purelin linear function for the output layer.In this paper, a plug-in hybrid electric city bus is studied which is running on pre-determined routes, thus the speed prediction method is for known driving profiles.Taking the UKBUS driving cycle as an example, the driving intention time series corresponding to the conditions are obtained by the aforementioned driving intention identification method.A total of 3288 samples of vehicle speed-intention data, choosing 2788 samples as training samples among which 75% is considered as training data, and the remaining 25% is used as validation data and test data.The remaining 500 samples which didn't participate in the training are used to test the network's predictive ability.

Short-Term Vehicle Speed Prediction Using Nonlinear Auto-Regressive Neural Network
The NAR neural network is a kind of dynamic neural network, which could add the output memory to the input automatically to calculate the output at the next instance.As the driver's operation can reflect the future speed variation trend, the obtained time series of driving intention are provided as input to the network together with the past vehicle speed data.
Definition of the prediction model is given in Equation ( 9): where Y(t) is the output at a time t, f is the nonlinear mapping, dy is the output memory order.The output at the time (t + 1), i.e., Y(t + 1) depends on the former dy steps of output Y(t), as the predicted values could be fed back to the input, the long-term predictions can be achieved.
The NAR neural network structure, consisting a hidden layer, an output layer, and an output feedback layer, is given in Figure 8.As shown in Figure 8, the hidden and the output layer include 15 and 1 neurons respectively.The order of output memory is 5. w 1 is the connection weight between input and output layer neuron, b 1 is the threshold of the layer, f 1 is the transfer function of the hidden layer, w 2 is the connection weight between hidden layer and output layer, b 2 is the threshold of output layer and f 2 is the transfer function of output layer.
This paper uses a Bayesian algorithm to train the network, choosing Tan-sigmoid transfer function as hidden layer function, and purelin linear function for the output layer.In this paper, a plug-in hybrid electric city bus is studied which is running on pre-determined routes, thus the speed prediction method is for known driving profiles.Taking the UKBUS driving cycle as an example, the driving intention time series corresponding to the conditions are obtained by the aforementioned driving intention identification method.A total of 3288 samples of vehicle speed-intention data, choosing 2788 samples as training samples among which 75% is considered as training data, and the remaining 25% is used as validation data and test data.The remaining 500 samples which didn't participate in the training are used to test the network's predictive ability.The predicted speed segmentation with a receding horizon of 5 s is shown in Figure 9.As can be observed, the NAR-based speed predictor predicts the micro-trip behaviors effectively.Root mean square error (RMSE) is used in this paper to estimate the error between the prediction value and the actual value.Smaller the RMSE means higher prediction accuracy.In Figure 10, the comparison between the real value and the single-step prediction value of vehicle velocity segment is given.It can be seen in the figure that the error falls into (−2.5, 2.5) km/h and has good prediction results, based on which the length of the step is increased further.The RMSE results within different prediction horizon are shown in Table 3.Along with the increase in estimation time of vehicle velocity, the prediction error is in the increasing trend.And in order to verify the effect of introducing driving intention data as NAR input, the RMSE results of the NAR prediction model with the same parameters which only use the past speed data are also shown in Table 3.As can be seen, the RMSE results of prediction model without using driving intention data is higher, which reflects that the introduction of driving intention recognition results as NAR input can improve the prediction accuracy.The predicted speed segmentation with a receding horizon of 5 s is shown in Figure 9.As can be observed, the NAR-based speed predictor predicts the micro-trip behaviors effectively.The predicted speed segmentation with a receding horizon of 5 s is shown in Figure 9.As can be observed, the NAR-based speed predictor predicts the micro-trip behaviors effectively.Root mean square error (RMSE) is used in this paper to estimate the error between the prediction value and the actual value.Smaller the RMSE means higher prediction accuracy.In Figure 10, the comparison between the real value and the single-step prediction value of vehicle velocity segment is given.It can be seen in the figure that the error falls into (−2.5, 2.5) km/h and has good prediction results, based on which the length of the step is increased further.The RMSE results within different prediction horizon are shown in Table 3.Along with the increase in estimation time of vehicle velocity, the prediction error is in the increasing trend.And in order to verify the effect of introducing driving intention data as NAR input, the RMSE results of the NAR prediction model with the same parameters which only use the past speed data are also shown in Table 3.As can be seen, the RMSE results of prediction model without using driving intention data is higher, which reflects that the introduction of driving intention recognition results as NAR input can improve the prediction accuracy.Root mean square error (RMSE) is used in this paper to estimate the error between the prediction value and the actual value.Smaller the RMSE means higher prediction accuracy.In Figure 10, the comparison between the real value and the single-step prediction value of vehicle velocity segment is given.It can be seen in the figure that the error falls into (−2.5, 2.5) km/h and has good prediction results, based on which the length of the step is increased further.The RMSE results within different prediction horizon are shown in Table 3.Along with the increase in estimation time of vehicle velocity, the prediction error is in the increasing trend.And in order to verify the effect of introducing driving intention data as NAR input, the RMSE results of the NAR prediction model with the same parameters which only use the past speed data are also shown in Table 3.As can be seen, the RMSE results of prediction model without using driving intention data is higher, which reflects that the introduction of driving intention recognition results as NAR input can improve the prediction accuracy.The predicted speed segmentation with a receding horizon of 5 s is shown in Figure 9.As can be observed, the NAR-based speed predictor predicts the micro-trip behaviors effectively.Root mean square error (RMSE) is used in this paper to estimate the error between the prediction value and the actual value.Smaller the RMSE means higher prediction accuracy.In Figure 10, the comparison between the real value and the single-step prediction value of vehicle velocity segment is given.It can be seen in the figure that the error falls into (−2.5, 2.5) km/h and has good prediction results, based on which the length of the step is increased further.The RMSE results within different prediction horizon are shown in Table 3.Along with the increase in estimation time of vehicle velocity, the prediction error is in the increasing trend.And in order to verify the effect of introducing driving intention data as NAR input, the RMSE results of the NAR prediction model with the same parameters which only use the past speed data are also shown in Table 3.As can be seen, the RMSE results of prediction model without using driving intention data is higher, which reflects that the introduction of driving intention recognition results as NAR input can improve the prediction accuracy.

Modeling of Mixed Logic Dynamical Predictive Control Strategy
The short-term future vehicle speed in each receding horizon can be obtained through the proposed speed prediction method in Section 3. Then the predicted vehicle speed is used to calculate the vehicle required torque at each moment within the receding horizon.The proposed MLD-MPC control strategy aims to distribute the vehicle required torque between motor and engine reasonably at each moment within the prediction horizon to achieve optimal fuel economic.Thus MLD model is established to tackle this optimization problem with the objective function of equivalent fuel consumption, the control variable of motor torque u(k) and the state variable of SOC x(k).The optimal control variable sequence are obtained by minimizing the total equivalent fuel consumption in the receding horizon through MLD model and only the first control variable need to be applied in the controlled PHEV model at the present moment.The MLD-MPC receding horizon control procedure is described as follows: (1) at time k, predict the short-term future speed profile for the current control horizon (k-k + N, where N is the receding horizon) and calculate the corresponding required torque through Equations ( 1) and ( 2); (2) the MLD model calculates the optimal control policy (u(k)-u(k + N)) for the current prediction horizon (k-k + N); (3) apply the first time-step of the optimal control policy u(k) in the controlled PHEV model; (4) update the state variable and system constraints, repeat the computation procedure 1-3 at the next time instant (time k + 1).
MLD model is suitable for solving problems combining binary variables and continuous variables.In [28], a constrained optimal problem is formulated and solved based on MPC using the simplified MLD model for DC-DC boost converter.Similarly, in our PHEV energy management control strategy, based on the simplified powertrain model in Section 2, the state transition equation and evaluation equation of MLD model can be described in Equation (10): (10) where, where x(k) is the SOC of battery at time k, y(k) is the equivalent rate of fuel consumption at time k, δ(k) is the operation mode matrix at time k, which is a 6 × 1 matrix whose entries are either 0 or 1 (logic variables), and only one of the six elements of the matrix is 1 and others are all 0. z(k) is the auxiliary variable, z(k) = δ(k) • u(k), in which u(k) is the control variable which represents the torque of electric motor, and its inequality constraints are described in Equation (11).i t is the transmission ratio of the torque coupler which connects the engine to the motor.In the equation, the variables are obtained based on the linearized engine and battery model in Section 2. B 1 , B 2 is used to calculate the SOC change corresponding to the specific motor torque, b 0 (k) and b 1 (k) are the fitting coefficients for the specific motor speed at time k, and they are obtained from Equation (8).Similarly, D 1 and D 2 is used to calculate the sum of engine fuel consumption and motor equivalent fuel consumption, a 0 (k) and a 1 (k) are from Equation (4): The inequality constraint of the six operation modes for PHEV presented in Table 4 includes electric power mode, fuel consumption mode, mixed driven mode, drive charging mode, regenerate braking mode and stopping mode.The constraint describes the relationship of motor torque u and the vehicle required torque T in in each operation mode.In Table 4, M m and M in are the maximum motor torque and vehicle required torque respectively; T e_max and T m_max are the maximum engine torque and motor torque respectively; T dec_max is the maximum regenerative braking torque; T dec_max = max(T in /i t , −T m_max ); ε is the machine accuracy, ε = 0.0001; T in is the required torque, which can be can be calculated from Equations ( 1) and ( 2); δ i is the logic variable representing the operation mode of the vehicle (if δ i = 1, it means that the vehicle is running under the operation mode i); Consider the inequality constraints above into time sequence k and transfer them into the unified matrix inequality as follows, where E 1 , E 2 , E 3 and E 4 are all the coefficients in the form of Matrix converted from Table 4.

Electric Power Mode
Drive Charging Mode Applying receding horizon optimization based predictive control theory to the MLD model, and the optimal model is set-up as given in Equations ( 13) and ( 14), where the optimization metric is the sum of equivalent fuel consumption in the receding horizon at time k, the state transition and evaluation equality constraint are from Equation (10), the inequality constraints are taken from Equation ( 12): min within it, B 1 , B 2 , D 1 , D 2 is the same as Equation (10), and E 1 , E 2 , E 3 , E 4 , is defined as follows: In the equation, N is the prediction horizon, k is the specific second in the prediction horizon, x min (k) and x max (k) are the upper and lower limit of SOC accessible domains respectively at time k.Considering the values of matrices B 1 , B 2 , D 1 , D 2 , E 1 , E 2 , E 3 and E 4 at every sampling time k and according to the objective function of minimal equivalent fuel consumption, Equations ( 13) and ( 14) can be converted into a MILP problem.

Solution of Mixed Integer Linear Programming
The MILP problem can be modeled using YALMIP toolbox in the Matlab platform and solved using the Gurobi optimizer.The optimal control variable sequence (u(k)-u(k + N)) in the prediction horizon (k-k + N) can be calculated at time k, and the first value of the optimal control variable sequence u(k) is applied to the controlled object.Then utilizing the receding horizon optimization method, the former solving steps are repeated to get a new control value at time k + 1 to realize the PHEV MPC.In the process of updating the model state variable, SOC, the original nonlinear model as shown in Equation ( 6) is used instead of the simplified linear equality constraint model in MILP to obtain a better estimate of the true SOC value at the end of the control time domain.The obtained value is considered as the initial state value to the next receding horizon optimization, so as to realize the feedback correction in predictive control.
Heuristic rules for some special driving conditions are set to reduce the range of feasible solution of MILP and improve the real-time performance of the control strategy.The introduced heuristic rules are presented in Table 5.While the condition in the first column is satisfied, the corresponding operation mode δ(k) and motor torque can be decided respectively.Total required torque T in can be used as a measure to judge whether the car is braking or stopping.Based on the value of total required torque, whether it is less than or equal to zero, the vehicle can choose regenerate braking mode or stop mode directly, and the corresponding control variable value can be determined directly.When the vehicle is in urgent acceleration or at high speed, and the total inquired torque is larger than the maximum engine torque at that speed, the vehicle can choose mixed driven mode directly.Otherwise, the control strategy can switch among four modes, i.e., electric power mode, fuel consumption mode, mixed driven mode, drive charging mode.The control variable value in electric power mode and fuel consumption mode can be directly obtained.The control variable in the mixed driven mode and drive charging mode can be selected through the optimization of the objective function.These heuristic rules can significantly reduce the search space of MILP and thus improve the efficiency of the algorithm.

Driving Condition Operation Mode Motor Torque Decision
When SOC approaches to the lower limit, once the required torque is larger than the maximum torque of the engine and the objective function in MLD-MPC is not modified, the motor needs to output the power which would make the SOC of the battery to exceed the lower limit.This violates the constraints defined in MILP problem and cannot find a feasible solution.Thus, when SOC is around the lower limit, the heuristic rule is applied to solve the control problem.First, estimate the highest SOC at each step in the prediction horizon under the assumption that the engine outputs its maximum torque.If the estimated highest SOC is higher than the lower limit of SOC, the current SOC constraints are imposed in MILP, otherwise the corresponding SOC constraints will be removed.When the current SOC is lower than the lower limit of SOC, the algorithm will stop solving the MILP problem and change to the heuristic solver mode.In this mode, the maximum charging torque (or the minimum discharging torque) of the motor is calculated directly while the engine is set to output its maximum torque.After the SOC rises to a certain value higher than the lower limit of SOC, the MILP solver will be restarted.

Simulation Experiments
Parameters of the PHEV are given in Table 6.The diesel engine displacement is 7.3 L with a rated power of 177 kW, and maximum engine speed is 2300 r/min.The type of motor is AC induction motor, with rated power of 124 kW, and maximum motor speed of 10,000 r/min.The type of battery used is a lithium battery with a rated capacity of 90 A • h.M m is 272 N•m, and M in is 1170 N•m.Fitting coefficients of engine and motor are listed in the Appendix A. ADVISOR (Advanced Vehicle Simulator, 2002, National Renewable Energy Laboratory, Golden, CO, USA) is used to simulate the UKBUS driving cycle as shown in Figure 11, with an initial SOC value of 0.7 (0.3 ≤ SOC ≤ 0.8), prediction horizon N of 15 s and sampling period of 1 s (the effect of predict time domain is discussed in Section 5.2).The simulation results are presented in Figure 12. Figure 12a is the comparison of the SOC curves between the MLD-MPC control strategy proposed in this paper and Charge Depleting-Charge Sustaining (CD-CS) control strategy.As can be observed, the CD-CS strategy consumes the battery energy within 8000 s, and then sustains SOC around 0.3 for the remainder of the trip.In the MLD-MPC control strategy, the battery energy is completely consumed by the end of the trip.Compared to CD-CS control strategy, the proposed MLD-MPC strategy can make more reasonable use of battery energy.Figure 12b is the comparison of the fuel consumption between the MLD-MPC control strategy presented in this paper and CD-CS control strategy.As seen in the picture, the CD-CS control strategy makes full use of battery energy and does not consume fuel within the first 8000 s, then the fuel consumption grow rapidly during the Charge-Sustaining mode, while the MLD-MPC control strategy maintains a relatively slow fuel consumption rate.Up to the end of the whole driving cycle, the fuel consumption is 12.93 L, and the end value of SOC is 0.3043 for the MLD-MPC control strategy, while the fuel consumption is 15.39 L, and the end value of SOC is 0.3125 for the CD-CS control strategy.It can be concluded from this results that under the same driving cycle conditions, the MLD-MPC control strategy reduces the fuel consumption of a PHEV. Figure 13 shows the comparison of engine operation points of these two control strategies.As can be observed, the engine operation points of the CD-CS control strategy spread out on a much larger area, and a lot of operation points in the low efficiency area.While the engine operation points of MLD-MPC control strategy is relatively centralized and most of them are in the high efficiency area which will improve fuel economic.
Besides the simulation of UKBUS driving cycle, driving cycles of Orange County Bus Cycle (OCC) and New European Drive Cycle (NEDC) are simulated as well.The results can be seen in Table 7, which shows that the fuel economy of PHEV under these driving cycles is improved to a noticeable extent.), prediction horizon N of 15 s and sampling period of 1 s (the effect of predict time domain is discussed in Section 5.2).The simulation results are presented in Figure 12. Figure 12a is the comparison of the SOC curves between the MLD-MPC control strategy proposed in this paper and Charge Depleting-Charge Sustaining (CD-CS) control strategy.As can be observed, the CD-CS strategy consumes the battery energy within 8000 s, and then sustains SOC around 0.3 for the remainder of the trip.In the MLD-MPC control strategy, the battery energy is completely consumed by the end of the trip.Compared to CD-CS control strategy, the proposed MLD-MPC strategy can make more reasonable use of battery energy.Figure 12b is the comparison of the fuel consumption between the MLD-MPC control strategy presented in this paper and CD-CS control strategy.As seen in the picture, the CD-CS control strategy makes full use of battery energy and does not consume fuel within the first 8000 s, then the fuel consumption grow rapidly during the Charge-Sustaining mode, while the MLD-MPC control strategy maintains a relatively slow fuel consumption rate.Up to the end of the whole driving cycle, the fuel consumption is 12.93 L, and the end value of SOC is 0.3043 for the MLD-MPC control strategy, while the fuel consumption is 15.39 L, and the end value of SOC is 0.3125 for the CD-CS control strategy.It can be concluded from this results that under the same driving cycle conditions, the MLD-MPC control strategy reduces the fuel consumption of a PHEV. Figure 13 shows the comparison of engine operation points of these two control strategies.As can be observed, the engine operation points of the CD-CS control strategy spread out on a much larger area, and a lot of operation points in the low efficiency area.While the engine operation points of MLD-MPC control strategy is relatively centralized and most of them are in the high efficiency area which will improve fuel economic.
Besides the simulation of UKBUS driving cycle, driving cycles of Orange County Bus Cycle (OCC) and New European Drive Cycle (NEDC) are simulated as well.The results can be seen in Table 7, which shows that the fuel economy of PHEV under these driving cycles is improved to a noticeable extent.

Influence of Prediction Horizon for Mixed Logic Dynamical-Model Predictive Control (MLD-MPC) Strategy
Two key steps are included in the MLD-MPC control strategy proposed in this paper, which are NAR neural network based vehicle speed prediction and minimal fuel consumption objective

Influence of Prediction Horizon for Mixed Logic Dynamical-Model Predictive Control (MLD-MPC) Strategy
Two key steps are included in the MLD-MPC control strategy proposed in this paper, which are NAR neural network based vehicle speed prediction and minimal fuel consumption objective

Influence of Prediction Horizon for Mixed Logic Dynamical-Model Predictive Control (MLD-MPC) Strategy
Two key steps are included in the MLD-MPC control strategy proposed in this paper, which are NAR neural network based vehicle speed prediction and minimal fuel consumption objective function-based MLD model prediction control.With an increase in prediction horizon, the accuracy of prediction results will decrease compared to the actual speed and thus affect the accuracy of the MPC solution and increase the fuel consumption.The MPC uses the online receding horizon optimization within the limited horizon, and when the optimal horizon increases, the optimization performance index also increases which leads to better results.Therefore, it is important to find out the proper prediction horizon for the optimal fuel consumption simulation.
Here, NEDC is used to study the effect of prediction horizon length on target equivalent fuel consumption.And the same NAR neural network training method as for UKBUS is utilized for NEDC.With a total of 1190 samples of vehicle speed-intention data, the first 900 samples are chosen as training samples, and the remaining 290 samples which didn't participate in the training are used to test the network's predictive ability.When the value of SOC does not drop to the lower limit value, the constraints of state variable SOC are not activated.The optimal control variable can be solved independently, and the function of prediction horizon cannot be reflected, as the state association of each second is not restricted by constraints within the prediction horizon.The function of prediction horizon is realized only when the value of SOC is going to cross the lower limit.Therefore, initial SOC value of 0.7 is selected and the lower limit of SOC is set to 0.67 to simulate the condition that SOC is going to reach the lower limit.In this phase, heuristic solver mode will be activated to recharge the batteries up to 0.68 when SOC falls below 0.67 and then MILP solver mode will be selected.
Figure 14 shows the battery SOC curve under different prediction horizons, and the corresponding time to reach the lower limit of SOC value is shown in Table 8.It can be seen that the time to reach the lower limit of SOC is extended as the prediction horizon increases.That is to say, in the process of optimization with MILP, it can predict in advance that SOC is about to reach the lower limit so the strategy can choose to save electricity in advance and improve the fuel efficiency.Figure 15 shows the equivalent fuel consumption per hundred kilometers corresponding to the different prediction horizons.
Energies 2017, 10, 74 15 of 18 function-based MLD model prediction control.With an increase in prediction horizon, the accuracy of prediction results will decrease compared to the actual speed and thus affect the accuracy of the MPC solution and increase the fuel consumption.The MPC uses the online receding horizon optimization within the limited horizon, and when the optimal horizon increases, the optimization performance index also increases which leads to better results.Therefore, it is important to find out the proper prediction horizon for the optimal fuel consumption simulation.
Here, NEDC is used to study the effect of prediction horizon length on target equivalent fuel consumption.And the same NAR neural network training method as for UKBUS is utilized for NEDC.With a total of 1190 samples of vehicle speed-intention data, the first 900 samples are chosen as training samples, and the remaining 290 samples which didn't participate in the training are used to test the network's predictive ability.When the value of SOC does not drop to the lower limit value, the constraints of state variable SOC are not activated.The optimal control variable can be solved independently, and the function of prediction horizon cannot be reflected, as the state association of each second is not restricted by constraints within the prediction horizon.The function of prediction horizon is realized only when the value of SOC is going to cross the lower limit.Therefore, initial SOC value of 0.7 is selected and the lower limit of SOC is set to 0.67 to simulate the condition that SOC is going to reach the lower limit.In this phase, heuristic solver mode will be activated to recharge the batteries up to 0.68 when SOC falls below 0.67 and then MILP solver mode will be selected.
Figure 14 shows the battery SOC curve under different prediction horizons, and the corresponding time to reach the lower limit of SOC value is shown in Table 8.It can be seen that the time to reach the lower limit of SOC is extended as the prediction horizon increases.That is to say, in the process of optimization with MILP, it can predict in advance that SOC is about to reach the lower limit so the strategy can choose to save electricity in advance and improve the fuel efficiency.Figure 15 shows the equivalent fuel consumption per hundred kilometers corresponding to the different prediction horizons.

Bottoming Time Prediction Horizon
First Bottoming Second Bottoming The change of RMSE of the predicted speed with different prediction horizons is shown in Figure 16.Compared to Figure 15, with an increase in prediction horizon, the equivalent fuel deceases continuously until the prediction horizon reaches 15 s, while the equivalent fuel consumption deceased only a little when the prediction horizon is increased from 15 s to 30 s.However, RMSE of the predicted speed is significant increased when the prediction horizon is increased from 15 s to 30 s.Therefore, a prediction horizon of 15 s is chosen in this paper.

Conclusions
A MPC strategy for PHEV based on a MLD model and vehicle speed prediction was proposed in this paper.Firstly, to convert the PHEV control problem into a linear programming problem, the energy consumption model of engine and motor were both linearized.After that, focusing on a city bus which usually has determined routes and known velocity profiles, a vehicle speed prediction method using a NAR neural network based on the combination of driving intention and historical speed data was proposed, in which the driving intention is recognized by using fuzzy inference to analyze the driver's pedal operation.Then, a vehicle MLD model is established with constraints in six possible driving modes and an objective function of minimum equivalent fuel consumption.With the predicted vehicle speed, a MPC method was utilized and the energy management problem was converted to a MILP problem, and the optimal motor torque and the corresponding driving mode sequence within the speed prediction horizon can be solved.The simulation was carried out under three driving cycles, and the effect of prediction horizon on the result of the simulation is analyzed, the simulation results verified the feasibility and effectiveness of the vehicle speed prediction based MLD-MPC strategy.DUT15LK13).The change of RMSE of the predicted speed with different prediction horizons is shown in Figure 16.Compared to Figure 15, with an increase in prediction horizon, the equivalent fuel deceases continuously until the prediction horizon reaches 15 s, while the equivalent fuel consumption deceased only a little when the prediction horizon is increased from 15 s to 30 s.However, RMSE of the predicted speed is significant increased when the prediction horizon is increased from 15 s to 30 s.Therefore, a prediction horizon of 15 s is chosen in this paper.The change of RMSE of the predicted speed with different prediction horizons is shown in Figure 16.Compared to Figure 15, with an increase in prediction horizon, the equivalent fuel deceases continuously until the prediction horizon reaches 15 s, while the equivalent fuel consumption deceased only a little when the prediction horizon is increased from 15 s to 30 s.However, RMSE of the predicted speed is significant increased when the prediction horizon is increased from 15 s to 30 s.Therefore, a prediction horizon of 15 s is chosen in this paper.

Conclusions
A MPC strategy for PHEV based on a MLD model and vehicle speed prediction was proposed in this paper.Firstly, to convert the PHEV control problem into a linear programming problem, the energy consumption model of engine and motor were both linearized.After that, focusing on a city bus which usually has determined routes and known velocity profiles, a vehicle speed prediction method using a NAR neural network based on the combination of driving intention and historical speed data was proposed, in which the driving intention is recognized by using fuzzy inference to analyze the driver's pedal operation.Then, a vehicle MLD model is established with constraints in six possible driving modes and an objective function of minimum equivalent fuel consumption.With the predicted vehicle speed, a MPC method was utilized and the energy management problem was converted to a MILP problem, and the optimal motor torque and the corresponding driving mode sequence within the speed prediction horizon can be solved.The simulation was carried out under three driving cycles, and the effect of prediction horizon on the result of the simulation is analyzed, the simulation results verified the feasibility and effectiveness of the vehicle speed prediction based MLD-MPC strategy.

Conclusions
A MPC strategy for PHEV based on a MLD model and vehicle speed prediction was proposed in this paper.Firstly, to convert the PHEV control problem into a linear programming problem, the energy consumption model of engine and motor were both linearized.After that, focusing on a city bus which usually has determined routes and known velocity profiles, a vehicle speed prediction method using a NAR neural network based on the combination of driving intention and historical speed data was proposed, in which the driving intention is recognized by using fuzzy inference to analyze the driver's pedal operation.Then, a vehicle MLD model is established with constraints in six possible driving modes and an objective function of minimum equivalent fuel consumption.With the predicted vehicle speed, a MPC method was utilized and the energy management problem was converted to a MILP problem, and the optimal motor torque and the corresponding driving mode sequence within the speed prediction horizon can be solved.The simulation was carried out under three driving cycles, and the effect of prediction horizon on the result of the simulation is analyzed, the simulation results verified the feasibility and effectiveness of the vehicle speed prediction based MLD-MPC strategy.

fA
is the frontal area in 2 m , d C is the drag coefficient, r f represents tire rolling resistance coefficient, m is the vehicle mass in kg, J represents total vehicle rotational inertia in 2 kg m ⋅ , α is

Figure 2 .Figure 1 .
Figure 2. Linear fitting curve of engine fuel consumption rate and torque.

2 m
, d C is the drag coefficient, r f represents tire rolling resistance coefficient, m is the vehicle mass in kg, J represents total vehicle rotational inertia in 2 kg m ⋅ , α is vehicle slope angle in rad, d d v / t is the acceleration in m/s 2 .Model of transmission torque and speed is given in Equation (

Figure 2 .
Figure 2. Linear fitting curve of engine fuel consumption rate and torque.

Figure 2 .
Figure 2. Linear fitting curve of engine fuel consumption rate and torque.

Figure 3 .
Figure 3. Fitted curves of motor torque and state of charge (SOC) change.

Figure 3 .
Figure 3. Fitted curves of motor torque and state of charge (SOC) change.

Figure 4 .
Figure 4. (a) Membership functions of acceleration pedal travel; (b) membership functions of rate of change of acceleration pedal travel; and (c) membership functions of acceleration intention.

Figure 4 .
Figure 4. (a) Membership functions of acceleration pedal travel; (b) membership functions of rate of change of acceleration pedal travel; and (c) membership functions of acceleration intention.

Figure 5 .
Figure 5. (a) Membership functions of braking pedal travel; (b) membership functions of rate of change of braking pedal travel; and (c) membership functions of braking intention.

Figure 5 .
Figure 5. (a) Membership functions of braking pedal travel; (b) membership functions of rate of change of braking pedal travel; and (c) membership functions of braking intention.

Figure 6 .
Figure 6.The speed segment of United Kingdom Bus (UKBUS) driving cycle.

Figure 7 .
Figure 7.The driving intention identification result comparison.

Figure 6 .
Figure 6.The speed segment of United Kingdom Bus (UKBUS) driving cycle.

Figure 6 .
Figure 6.The speed segment of United Kingdom Bus (UKBUS) driving cycle.

Figure 7 .
Figure 7.The driving intention identification result comparison.

Figure 7 .
Figure 7.The driving intention identification result comparison.

Figure 9 .
Figure 9. Short-term predicted speed compared with the real speed.

Figure 9 .
Figure 9. Short-term predicted speed compared with the real speed.

Figure 9 .
Figure 9. Short-term predicted speed compared with the real speed.

Figure 12 .Figure 13 .
Figure 12.(a) comparison of SOC curves of two strategies; and (b) comparison of the fuel consumption of two strategies.

Figure 12 .Figure 12 .Figure 13 .
Figure 12.(a) comparison of SOC curves of two strategies; and (b) comparison of the fuel consumption of two strategies.

Figure 13 .
Figure 13.(a) Engine operation points of charge depleting-charge sustaining (CD-CS) strategy; and (b) engine operation points of mixed logic dynamical-model predictive control (MLD-MPC) Strategy.

Figure 14 .
Figure 14.Battery SOC curve under different prediction horizons.

Figure 15 .
Figure 15.Equivalent fuel consumption under different prediction horizons.

Table 1 .
The fuzzy rules for identification of the acceleration intention.

Table 1 .
The fuzzy rules for identification of the acceleration intention.

Table 2 .
The fuzzy rules for identification of the braking intention.

Table 2 .
The fuzzy rules for identification of the braking intention.

Table 3 .
Prediction results comparison of NAR neural network at different prediction horizon.

Table 4 .
Constraints during various operation mode for PHEV.

Table 6 .
Parameters of the PHEV.

Table 6 .
Parameters of the PHEV.The type of motor is AC induction motor, with rated power of 124 kW, and maximum motor speed of 10,000 r/min.The type of battery used is a lithium battery with a rated capacity of 90 A h ⋅ .Mm is 272 N•m, and Min is 1170 N•m.Fitting coefficients of engine and motor are listed in the Appendix A. ADVISOR (Advanced Vehicle Simulator, 2002, National Renewable Energy Laboratory, Golden, CO, USA) is used to simulate the UKBUS driving cycle as shown in Figure11, with an initial SOC

Table 7 .
Fuel consumption under two control strategies.NEDC: New European Drive Cycle; OCC: Orange County Bus Cycle.

Table 7 .
Fuel consumption under two control strategies.NEDC: New European Drive Cycle; OCC: Orange County Bus Cycle.

Table 7 .
Fuel consumption under two control strategies.NEDC: New European Drive Cycle; OCC: Orange County Bus Cycle.

Table 8 .
SOC bottoming time under different prediction horizons.

Table 8 .
SOC bottoming time under different prediction horizons.
Figure 15.Equivalent fuel consumption under different prediction horizons.