A Robust Design of the Model-Free-Adaptive-Control-Based Energy Management for Plug-in Hybrid Electric Vehicle

: This paper proposes a robust design approach based on the Design for Six Sigma (DFSS), to promote the robustness of our previous model-free-adaptive-control-based (MFAC-based) energy management strategy (EMS) for the plug-in hybrid electric vehicles (PHEVs) in real-time application. First, the multi-island genetic algorithm (MIGA) is employed for a deterministic design of the MFAC-based EMS, and the Monte Carlo simulation (MCS) is utilized to evaluate the sigma level of the strategy with the deterministic design results. Second, a DFSS framework is formulated to reinforce the robustness of the MFAC-based EMS, in which the velocity and the vehicle mass are considered external disturbances whilst the terminal state of charge (SOC) of the battery and the fuel consumption (FC) are conducted as responses. In addition, real-time SOC constraints are incorporated into Pontryagin’s minimum principle (PMP) to conﬁne the ﬂuctuation of battery SOC in MFAC-based EMS to make it closer to the solution of the dynamic programming (DP). Finally, the effectiveness of the robust design results is assessed by contrasting with other strategies for various combined driving cycles (including velocity, vehicle mass, and road slope). The comparisons demonstrate the remarkable promotion of the robust design in terms of the energy-saving potential and the performance against external disturbance. The average improvement of the FCs can reach up to a considerable 19.66% and 9.79% in contrast to the charge-depleting and charge-sustaining (CD-CS) strategy as well as the deterministic design of MFAC-based EMS. In particular, the energy-saving performance is comparable to DP, where there is only a gap of − 1.68%. designed, the MFAC-based EMS can have a satisfactory performance in stochastic future driving cycles. The results demonstrate


Introduction
Electric vehicles (EVs) have been accepted as a global consensus owing to the shortage of petroleum resources and environmental concerns [1,2]. Nonetheless, the cost and energy density of batteries have considerably limited the development of EVs. Therefor, plug-in hybrid electric vehicles (PHEVs) have become a promising solution due to their significantly lower capacity of battery requirement and no range anxiety compared to pure battery electric vehicles (BEVs) [3]. To provide a preferable energy-saving potential, a scientific energy management strategy (EMS) is necessary to online optimize the distribution of the demanded power between the internal combustion engine (ICE) and electric motors (EMs) of the PHEVs [4]. However, it is still a challenging issue owing to the high nonlinearity of the hybrid powertrain, and is attracting numerous investigations in the literature [5].
Various methodologies have been provided to address the energy management problems, which can be categorized into rule-based strategies, optimization-based strategies, and learning-based strategies [6]. The most representative rule-based EMSs can be deterministic strategy [7], fuzzy logic strategy [8], and feedback control strategy [9], as well as charge-depleting and charge-sustaining (CD-CS) strategy, especially for PHEVs [10]. All of them can be easily implemented in real time due to their lower computational load and the battery state-of-charge (SOC) on battery voltage and resistance is neglected [44]. The optimal co-state (or equivalent factor) can be successfully derived based on prior known driving data. However, it may not be applicable for an unknown driving cycle. Therefore, adaptive or predictive approaches are commonly employed for the real-time application of PMP (or ECMS) [14,21,22,45]. Yet, the adaptive or predictive controller is only based on the typical historical driving cycles without considering the randomness of the driving cycle and vehicle mass, thereby leading to an unreasonable solution for real-time execution.
Owing to the straightforward control logic, feedback control is considered to be one of the most effective methods for real-time control of the EMSs [5,7,15,19,21,35], in which battery SOC is usually deployed as the feedback parameter [3,18,46]. Hence, it is crucial to provide a reasonable reference SOC trajectory close to the optimal battery SOC trajectory. Linear reference SOC planning is the simplest method to design a reference SOC trajectory. Once the whole trip distance and the battery SOC limitation (i.e., initial SOC and terminal SOC) are specified, the reference SOC trajectory can be easily attained [18,21,23]. Nevertheless, it may be extremely different from the optimal SOC trajectory since the actual battery SOC trajectory is strongly nonlinear. Hence, numerous nonlinear reference SOC planning approaches are provided to forecast the reference SOC trajectory [17,24,28,29,46]. However, a desirable estimation result is difficult to obtain unless the complexity of the method is increased, for which a higher computational burden may be aroused. This would be a significant drawback for real-time applications. On the other hand, trajectory tracking is another critical issue once the reference SOC trajectory is successfully designed. The PID-based (or other iterative-searching-based) controllers are constantly employed for tracking problems, whilst acceptable tracking performance can be provided [16,26,34,47]. Nonetheless, when the output parameter of the controller is designed to be the co-state of the PMP, unwanted fluctuations may be evoked, which can also lead to undesired increases in the energy consumption. Thus, a model-free-adaptive-control-based (MFAC-based) methodology was recommended to smooth the co-state in our previous research, in which the results demonstrated a meaningful improvement in the fuel economy of PHEVs [48]. However, it still needs to be strengthened. First, the reference SOC in the previous MFACbased EMS is designed to be a linear function, which may be distant from the optimal SOC trajectory, thus limiting the energy-saving potential. Second, the impact on the robustness of EMS due to stochastic driving cycles and vehicle mass is not yet considered. This may lead to unpredictable increases in energy consumption at different actual driving cycles with external disturbances in real time. Inspired by Refs. [49,50], this paper employs the Design for Six Sigma (DFSS) method to optimize the MFAC-based EMS in order to strengthen the robustness of the strategy for stochastic driving conditions in real-time application.
The main innovation and contribution are summarized as follows. (1) A deterministic optimization method based on multi-island genetic algorithm (MIGA) is proposed to optimize the MFAC-based EMS, and its robustness is evaluated by Monte Carlo simulation (MCS). (2) A robust design framework is proposed on the basis of the DFSS to enhance the adaptive capability of the MFAC-based EMS, in which a real-time SOC constraint derived from the linear reference SOC is incorporated into the PMP to further reinforce the energy-saving performance of the PHEV whilst operating in various driving conditions. The remainder of the paper is organized as follows. Section 2 presents the studied PHEV powertrain models and vehicle dynamic models. In Section 3, MFAC-based EMS is described, and optimized by MIGA, followed by the MCS analysis of the deterministic optimization results. The robust design framework of the MFAC-based EMS based on DFSS is formulated in Section 4, whilst the SOC constraint is also designed according to linear reference SOC. Section 5 introduces the actual driving cycles together with DP and rule-based strategy to verify the robustness of the proposed strategy. Finally, conclusions are summarized in Section 6.

Introduction of PHEV Powertrain
This study focused on a single-shaft parallel hybrid powertrain equipped with a multi-speed transmission. As shown in Figure 1, the hybrid powertrain of the PHEV is composed of an internal combustion engine (ICE), a diaphragm spring clutch, an electric motor integrated with an automated mechanical transmission (AMT), a final drive, and a lithium iron phosphate battery pack that can be charged from the power grid.

Engine
Clutch Electric Motor

AMT Final drive
Battery Pack Plug-in Since the electric motor can not only be utilized as a driving motor but also as a generator, the vehicle can be operated in five different modes to regulate the working state of the ICE and electric motor through the separation and combination of the clutch. The flexible working modes will be of great significance for the energy management of the PHEV to strengthen its energy-saving potential. In this work, the hybrid powertrain is applied to a bus that is 12 m long with a curb weight of 12,500 kg, and the key parameters of the PHEV are provided in Table 1. To design the EMS for PHEV, the control-oriented models including engine, motor, AMT, battery, and longitudinal dynamics are established for simulation.

Powertrain Models
Both the engine and motor are considered as the steady-state numerical model in the exploitation of an EMS. The displacement of the engine is 6.75 L, and its peak power is 162 kW. The brake special fuel consumption (BSFC) map of the engine is plotted in Figure 2a, according to the steady-state experiment data. The fuel consumption (FC) rate of the engine can be derived from the BSFC by the following equation [41,51].
whereṁ f represents the engine FC rate. b e can be obtained by looking up the BSFC map in accordance with the engine rotational speed n e and the engine torque T e . ρ f denotes the diesel density. The electric motor is a permanent magnet synchronous motor with a peak power of 130 kW and a peak rational speed of 4500 r/min. The efficiency map of the electric motor is indicated in Figure 2b. Since the electric motor can be served as both a driving motor and a generator, its efficiency for the driving mode and regenerative braking mode is considered to be the same whilst the torque is, respectively, defined as positive and negative. Accordingly, the efficiency of the electric motor can be expressed as the parametric function of the rotational speed and torque.
where n m and T m denote the rational speed and torque of the electric motor, respectively, and η m represents the mapping relationship of the motor efficiency. Then, the power of the electric motor P m can be expressed as follows [41,51].
where sgn(·) is a sign function, which can be defined by According to the characteristics of the hybrid powertrain, the rotation speed of the engine and motor should be consistent when it is operated in the hybrid mode. The relationship of the rotational speed between driving wheels and that of the two power sources is expressed as follows [48,51].
where n wheel is the rational speed of the driving wheels, and i g and i 0 denote the gear ratio of AMT and final drive, respectively. The AMT can not only assist in adjusting the working points of the engine and motor, but also downsize the output torque of the engine and motor. The AMT model can be described as follows.
where T out and T in denote the torque output and input from the AMT, respectively. η AMT represents the efficiency of the AMT for different gear. The lithium iron phosphate battery pack with a rated voltage of 354 V is utilized for the studied PHEV, and its nominal capacity is 50 A·h. As shown in Figure 3a, the internal-resistance equivalent model is utilized to simplify the dynamic characteristics of the battery pack due to that it has been widely adopted to describe the complicated dynamics of the charge and discharge process [10,12,17,21,48]. In this model, the open circuit voltage and the internal resistance are generally considered as the function of the battery SOC. The nonlinear characteristics of the battery dynamics are indicated in Figure  3b, according to the experimental data.  The equivalent circuit model of the battery pack is described as follows.
where P batt is the battery power consisting of the electrical power input or output of the battery. U oc , I, and R 0 represent the open circuit voltage, the battery current, and the internal resistance, respectively. By solving Equation (7) for the current I and according to the definition for the variation of the battery SOC, the dynamics of the battery pack system can be expressed as follows.
where SȮC represents the variation of the battery SOC, and Q batt denotes the nominal capacity of the battery.

Vehicle Dynamic Model
The longitudinal dynamic model is normally deployed to solve the energy consumption problem while the impact of the vehicle's lateral dynamics and vertical dynamics are neglected. The simplified vehicle dynamic model is described as follows [51].
where M represents the vehicle gross mass, consisting of the curb mass m r and vehicle load mass m p . T d , r wh , C d , A, ρ d , and θ represent the wheel torque, wheel radius, air drag coefficient, front area, air density, and road slope, respectively. v is the vehicle velocity. g, f r , and η T , respectively, represent the gravity acceleration, rolling resistance coefficient, and the efficiency of the transmission. T b is the mechanical braking torque on wheels provided by conventional friction brakes, while regenerative braking is insufficient to ensure the desired braking torque.

Description of the MFAC-Based EMS
The MFAC-based EMS is deployed as a hierarchical control architecture consisting of two layers. The lower layer is the PMP-based EMS, which is used to distribute the desired torque between engine and motor to optimize energy consumption. The upper layer based on MFAC is employed to online adjust the co-state of the PMP to achieve satisfactory energy-saving performance. It primarily consists of three modules, i.e., reference SOC planning module, PMP module, and MFAC module, as shown in Figure 4.

Reference SOC Planning
The reference SOC planning is simplified as a linear function of the traveled distance in the MFAC-based EMS, considering that the battery power is normally depleted at the end of the driving trip. Then, the reference SOC of the vehicle can be estimated according to the initial SOC and terminal SOC of the battery, which is expressed as follows.
where SOC re f , SOC i , and SOC t represent the linear reference SOC, initial SOC, and terminal SOC of the battery, respectively. D dist and D total represent the traveled distance and the whole trip distance, respectively.

Formulation of PMP
Since PMP can convert global optimization into an instantaneous one, it has been widely adopted to improve EMSs [14,15,17]. The optimization objective of the PMP for a PHEV can be represented only by the engine FC due to the electric power being commonly used up at the end of a driving mission [18,21,41]. The cost function is expressed as follows.
J= min where x(t) and u(t) are the state and control action, represents the instantaneous fuel consumption rate. Note that the battery SOC and battery power are critical parameters for EMS; they are taken as the state variable and control variable, respectively. Accordingly, the Hamiltonian function is presented as follows.
where λ is the co-state variable. The PMP can simplify the global optimization problem in capturing the optimal control sequence u * (t) at each instant of time to minimize the Hamiltonian function, which can be described as follows.
while the state variable and the co-state evolve in the following expression.
Furthermore, the constrained boundaries should be respected for minimizing the Hamiltonian, in which the initial and the terminal of the battery SOC, namely, SOC(t 0 ) and SOC(t f ), are prescriptive, respectively. Considering the physical limitations of the powertrain components, some boundaries are also designed to guarantee the desired power at the wheels, which can be described as follows: where n e (t) and n m (t) denote the rotational speed of the engine and the motor, whose lower and higher boundaries are limited by n e min (t), n e max (t), n m min (t), and n m max (t), respectively. P e (t) and P m (t) denote the output power of the engine and the motor, where the lower and higher boundaries are, respectively, confined by P e min (t), P e max (t), P m min (t), and P m max (t).

Formulation of MFAC
For PMP-based EMSs, the co-state is a significant parameter for achieving an outstanding control performance. Inspired by Ref. [52], a model-free-adaptive-control (MFAC) method is employed to online adjust the co-state of PMP based on a feedback control architecture, where the battery SOC is considered as the state variable. The MFAC consists of the control module and estimating module, and both of them can be converted into a discrete-time state. Thus, it can be utilized for online adjustment of the co-state for the PMP-based EMS in real time, and the control system is described as follows.
where SOC(k) and λ(k) denote the system state and output variable at time k, representing the current SOC and co-state, respectively. ϕ(k) is the bounded pseudo partial derivative (PPD) of the system. Then, regarding the one-step-ahead prediction error function that can address the excessive control effort of the MFAC, the control law of the system is expressed as follows.
where ρ k and λ k represent the step factor and weighting factor to ensure the universality of the method and confine the variation of the control variable, respectively. SOC * (k + 1) represents the reference SOC at k + 1 time.φ(k) is the estimation of the ϕ(k), which is determined by the following formula.
where η k and µ k represent the step factor and weighting factor, respectively. ε is a small positive constant.φ(1) is the initial value ofφ(k).

Deterministic Optimization of the MFAC-Based EMS
Once the parameters of the MFAC are determined, the MFAC-based EMS can be utilized for a real-time application. A deterministic method may be appropriate to calibrate the characteristic parameters of the MFAC to acquire a desirable co-state for PMP. Thereby, a multi-island genetic algorithm (MIGA) and a combined driving cycle (shown in Figure 5) containing velocity, road slope, and the random variation of the passenger numbers are employed to optimize the parameters of the MFAC. Note that the stochastic vehicle mass is composed of the curb mass and the passenger mass, which can be calculated by the numbers of the passenger, where the passenger quality is defined as 70 kg [53]. The combined driving cycle is acquired from an actual bus route. The main parameters and the constraints of MFAC are listed in Table 2. Table 2. Main parameters of the MFAC.

Descriptions Parameters Constraints
Step factor of the control law ρ k (0, 1] Weighting factor of the control law λ k (0, 1] Step factor of the PPD estimation η k (0, 1] Weighting factor of the PPD estimation µ k (0, 1] In the optimization process, the four parameters of the MFAC are considered as the design variables, while the fuel consumption of the engine is taken as the objective. Furthermore, the terminal SOC at the end of the trip is designed as the constraint, which is defined as a soft constraint of [0.25, 0.32] to avoid the overdischarge of the battery. Thus, the deterministic optimization model is governed by the following equation.
where FC dt represents the fuel consumption of the engine for deterministic optimization. The detailed description of the optimization process to obtain the optimal parameters of MFAC based on MIGA is illustrated in Table 3. Table 3. The deterministic optimization process.

MIGA Method to Optimize MFAC-Based EMS
1 Initialize the variables ρ k , λ k , η k , µ k and setting the constraints; 2 Output Λi through MFAC according to the prescribed ρ k , λ k , η k , µ k by MIGA at each iteration i; 3 Calculate the Hamilton function according to Equation (12) at each time j for each iteration i, and obtain the optimal control that minimizes H at each time j; 4 Update the state variable SOC based on the above optimized control variables for the time j + 1; 5 Repeat step 2 to 4 until j meets the time requirement whilst the state variable SOC also meets the predesigned limitation requirement; 6 Calculate the engine fuel consumption for iteration i; 7 Repeat step 2 to 6 until i meets the iteration requirement, and take the control parameters corresponding to the iteration with the minimum fuel consumption as the optimization results.
The process of the optimization can be seen in Figure 6, and the optimal parameters of MFAC are obtained corresponding to the minimum objective value where the optimization lasts for 668 iterations. The optimization results of the parameters are listed in Table 4.  Here, the optimal constant co-state (i.e., constant co-state) of the PMP with its corresponding results of the EMS are also calculated for comparison with the MFAC-based strategy. The performance of the MFAC-based EMS for the deterministic optimal solution can be seen in Figure 7.
The results indicate that the co-state can converge nicely to the optimal constant costate after an adaptive adjustment. This will be of significant advantage to decrease the fuel consumption of the PHEV. Similar to DP and optimal constant co-state control, the battery SOC of the MFAC-based EMS descends quickly at the initial stage since the vehicle is operated in the EV mode. With the adaptive adjustment of co-state, the vehicle switches to the hybrid mode, so that the descending speed of the battery SOC tends to be steady until the end of the trip. The terminal SOC of the MFAC-based strategy is 0.2526, and is slightly lower compared to the other two strategies, yet it still satisfies the limitation of the terminal SOC. As shown in Table 5, the FC of the MFAC-based EMS with a deterministic optimal solution is very close to the other two offline optimization methods. It means that the deterministic solution has considerable effective for the MFAC-based EMS. It is worth mentioning that the PMP with an optimal co-state can also achieve a perfect performance in energy saving in consistencywith DP, yet a little difference exists for their SOC trajectory and FC trajectory.

Monte Carlo Simulation for Deterministic Method
The MFAC-based EMS with a determined optimization has the acceptable performance for real-time energy management of the PHEV. However, the driving cycles cannot be known prior; they randomly change. Additionally, the randomness in the number of passengers may also have a considerable impact on the power requirements of the vehicle, which will lead to an inadaptability of the deterministic optimal solution. Consequently, the Monte Carlo simulation is deployed to estimate the reliability of the deterministic solution. As shown in Figure 8, both the driving cycles (composed of velocity and road slope) and vehicle mass (regarding the changing of the passengers' mass) are designed as the disturbance, where they are assumed as the normal distribution. After applying the disturbance to the MFAC-based EMS, the response variables are evaluated based on the designed feasible region, then the sigma level of the deterministic design can be derived.
In this paper, the terminal SOC is designed as the response, so that the reliability analysis can be transformed into determiningwhether the terminal SOC locates within the feasible region when the determined optimal solution is applied. Since the descriptive sampling method can not only reduce sampling times but also maintain the statistical quality of response analysis compared with the conventional random sampling method, it is selected to simulate the random disturbance in MCS. Then, the reliability of the determined design can be described as follows.
where R is the reliability, P f is the probability of failure, and N f and N t represent the number of points in the infeasible region and the total number of samples, respectively.  Moreover, some specific statistical variables, including the mean and standard deviation (i.e., Std. Dev.), of the terminal SOC are also evaluated by the following equations.

MFAC-based
where n is the total number of samples, and i represents the i th sample. The MSC results are shown in the following figures. From Figure 9a, it can be seen that many sampling points of the terminal SOC are located within the infeasible region, which implies that it may lead to the overdischarge or incomplete discharge of the battery. Both of them will have a disadvantageous impact on the vehicle operating. A similar conclusion can be obtained from Figure 9b. Although the mean value of the terminal SOC is 0.3, whilst the Std. Dev. is just 0.04, the reliability of the deterministic optimization results is still unacceptable owing to that the sigma level of the terminal SOC at the lower and upper boundary are only −1.244σ and 0.487σ. This also demonstrates that the deterministic optimal method has higher rates of failure probability, and it may be unsuitable for the characteristic parameters design of the MFAC-based EMS.

SOC Constraint Based on Linear Reference SOC
In an MFAC-based EMS, reference SOC is designed as a linear function. In fact, the actual battery SOC has strong fluctuant characteristics, which may give rise to promote the energy-saving potential. To investigate the battery SOC characteristics, different driving cycles are deployed (Figure 10), and the results for MFAC-based EMS and DP are shown in Figure 11.  The results reveal that the SOC fluctuation caused by different driving cycles and load changes is within a certain range, whether MFAC-based EMS or DP. Particularly, the interval can be thought of as parallel to the linear reference SOC, where the interval of the DP is much narrower in contrast to the MFAC-based EMS. Inspired by the conclusion, a battery SOC constraint based on linear reference SOC is designed to impose restrictions on the battery SOC caused by MFAC-based EMS. It means that the upper and lower boundary of the battery SOC for MFAC-based EMS can be defined according to the linear reference SOC. Therefore, the boundaries are expressed as follows.
where SOC high denotes the high boundary of the battery SOC, and α is an adjustment factor, herein α = 0.1.
where SOC low denotes the lower boundary of the battery SOC, and γ is utilized to adjust the lower boundary, which is defined as 0.15. Moreover, the minimum lower boundary is limited at 0.25, according to the terminal SOC of PMP. Then, a penalty factor is introduced in the Hamiltonian of PMP to confine the irregularfluctuation of battery SOC, for which an additive penalty function µ(SOC) is incorporated into Equation (12), which is rewritten as follows: where the penalty function is the form of an SOC piecewise function, which is described as where κ is a penalty factor to confirm the battery SOC trajectory within the constraints, for which it will be more close to the results of DP. Note that the value of µ(SOC) is zero when the battery SOC limits are satisfied, for which the original formulation (as shown in Equation (12)) will not be changed. Regardless of whether the battery SOC tends to go above or go below the boundary limit, the penalty factor κ will promote or impede the usage of electric energy to confirm the actual SOC within the predesigned interval.

Robust Design Based on DFSS
The diagram of the DFSS method for parameter optimization is shown in Figure 12. It is deployed as a two-layer framework. The outer layer is used for control parameters (i.e., ρ k , λ k , η k , µ k ) optimization of MFAC, for which the MIGA algorithm is employed. The inner layer is utilized for reliability analysis of the optimized results, which is constituted by MCS and MFAC-based EMS. It is necessary to note that the outer layer is the same as the deterministic optimization, which is described in detail in Section 3.2. The inner layer analyzes the robustness of the deterministic optimization results based on MCS, which can randomly generate external disturbance (i.e., stochastic driving cycles and random vehicle mass).  The optimization process is detailed as follows.
Step 1: The control parameters of MFAC are provided by the outer layer for each iteration while the disturbance is generated by the inner layer. A total of 200 experiments of the random vehicle mass, together with random velocity profile, are deployed by descriptive sampling method for each iteration.
Step 2: The fuel consumption and terminal SOCs are determined through MFACbased EMS based on the disturbance for each iteration.
Step 3: The sigma level of the terminal SOC together with the mean and the standard deviation of the FC are assessed by the MCS method. Then, they are considered as responses to be transmitted into the outer layer.
Step 4: Step 1 to Step 3 are continued until the number of the iterations is satisfied to the predesigned N. Then, the lowest F(X) among all of the solutions is considered as the optimal solution, whilst it should satisfy the constraints of the sigma level.
In this paper, the mean and the standard deviation of the FC are designed as the objective, while the sigma level of the terminal SOC is designed as a constraint. Hence, the optimization function can be expressed as follows.
where F(X) is the cost function, X is a vector composed of the control parameters of the MFAC and the external disturbance, µ FC and σ FC denote the mean and standard deviation of the fuel consumption, respectively, ω 1 and ω 2 denote the weighting factor of the µ FC and σ FC , and n is the sigma level, herein it is 6. It needs to be further explained that the control parameters of the MFAC and the constraints were defined in Equation (19). The external disturbance is mainly constituted by random velocity and vehicle mass. The random velocity is acquired from a collection of the driving cycles, as shown in Figure 13. Here, 35 groups of the velocity profile from an actual bus route are collected, for which 27 groups are used for DFSS design and the last 8 groups are employed to estimate the performance of the DFSS design results, subsequently. Accordingly, the variation of the driving cycle can be defined as follows.
where i should be a positive integer, and cycle i represents the i th driving cycle. The vehicle mass is randomly changing according to the variation of the bus stop interval. It is randomly generated by the optimal Latin hypercube design (Opt, LHD) method [48], for which the constraints are defined within a scope of [12,500,16,500] kg.

Robust Results
As shown in Figure 14, the terminal SOCs and the FCs corresponding to 200 experiments are returned at each iteration for MCS analysis. The variation of the terminal SOC has noteworthy distinctive characteristics at the iteration. The lower boundary can be significantly satisfied for each experiment, while the upper boundary is surpassed in many cases. It demonstrates that the SOC constraint has a remarkable promotion of the reliability for the minimum terminal SOC. Additionally, the terminal SOC also has a noticeable impact on the FC, for which the higher the terminal SOC, the larger the FC. Consequently, a robust design process lasts for 500 iterations to obtain the optimal solution. The changing of the terminal SOC for each iteration can be seen in Figure 15a. It is worth noting that the terminal SOC of each iteration is derived based on the mean of the external disturbance. That means that the mean value of the driving cycles and the vehicle mass based on the 200 experiments are employed to calculate the required results at each iteration. The robustness of the design is evaluated by the sigma level of the terminal SOC, as shown in Figure 15b. The lower boundary can be better guaranteed, while the sigma level of the upper boundary is unstablewith respect to the disturbance. Thus, the optimal results need to ensure that the cost function is minimal whilst the sigma level of the upper and lower boundaries are simultaneously satisfied. The optimal solution of the robust design appears at 412 iterations while the sigma level of the response is 8 for both the upper and lower boundary of the terminal SOC. Meanwhile, the minimum cost function can also be acquired at the same iteration (Figure 15c).
Accordingly, the optimal control parameters of the MFAC can be obtained at 412 iterations, whilst the optimal parameters are, respectively, 0.5476, 0.9401, 0.1490, and 0.7607, as shown in Figure 16.

Discussion
Eight groups of the combined driving cycles are deployed to evaluate the performance of the robust results whilst the road slope is also considered, according to Figure 5. As shown in Figure 17, the variation of the velocity and vehicle mass is considered as the main disturbance to the MFAC-based EMS.
To verify the effectiveness of the robust results, the DP results, as well as the solution based on the CD-CS, are also calculated and compared with the robust design.    Figure 17. The combined driving cycles for evaluation.

SOC Constraints
In this paper, the SOC constraints based on the linear reference SOC are designed to narrow the changing space of the battery SOC for a robust design. Figure 18 exhibits the performance of the SOC constraints concerning different driving cycles. The upper and the lower constraints of the battery SOC can be satisfactorily generated according to the linear reference SOC. Furthermore, they can successfully inhibit the battery SOC within the narrower boundary, thereby it will be more approximate to the optimal battery SOC trajectory. It should be pointed out that the constraints are integrated into the PMP whilst a penalty is added to determine whether the battery is working. Although the lower boundary is evidently less than 0.25 in some cases, it does not mean that the battery SOC will be below 0.25. It is worth noting that the lower boundary is designed to maintain at 0.25 near the end of the trip; this will be a benefit to ensure the terminal SOC at the expected value. On the other hand, the upper boundary is parallel to the linear reference SOC in most cases. Yet, it may have fluctuations in some cases (e.g., in NO. 6), thereby leading to an unsuccessful restriction to battery SOC. Nevertheless, the introduction of the SOC constraints can remarkably facilitate the robust design of the MFAC-based EMS.

Co-State
The comparisons of the co-state between deterministic design and robust design are shown in Figure 19. The results indicate that the adaptive ability of the robust design is considerably stronger than that of the deterministic design for which the convergence of the robust results is obviously better than deterministic solution under different driving cycles. It needs to be further clarified that the co-state of the MFAC-based EMS can gradually converge to the optimal constant co-state of the PMP which was confirmed in Ref. [48].
In other words, the performance of the MFAC-based EMS is significantly influenced by the convergence of the co-state with respect to different driving cycles. The better the convergence, the stronger the robustness. NO Figure 19. Comparisons of the co-state.

Battery SOC
The battery SOC can reflect the charge and discharge of the battery, which has a considerable impact on FC. Figure 20 provides the SOC trajectories of the robust design for various driving cycles. SOC trajectories for DP and CD-CS, as well as the deterministic design, are also deployed for comparison. DP has been accepted as the global optimization method, and the engine and motor will work more in blended mode during the whole trip. Therefore, the battery SOC shows a slow downward trend for all driving cycles, unlike the CD-CS. Regardless of whether for the robust design or the deterministic design,the overall trend of the battery SOC also descends until the end of the trip. However, there is more fluctuation for the deterministic design, which therefore may lead to a higher FC. Yet, the robust design can remarkably reduce the unexpected fluctuation, bringing it closer to the global optimum.
As shown in Table 6, the terminal SOC of the different strategies can satisfy the lower and upper limitation of 0.25 to 0.32. Nevertheless, the standard deviation of the terminal SOC for deterministic design and CD-CS is much larger than that of the robust design and DP. This reveals that the robust design has stronger stability to resist external disturbance, and its performance is closer to DP. Furthermore, the mean of the terminal SOC for the robust design is smaller than others. This also demonstrates that the electrical energy can be consumed more thoroughly, which will decrease fuel consumption.

Fuel Consumption
To examine the energy-saving performance of the robust design, the engine operating points of the robust design and deterministic design are plotted in Figure 21, as well as the DP and CD-CD for comparison. Four strategies are executed under different driving cycles, and the FCs are cumulatively computed based on the BSFC. It can be seen that the working points of the strategies are comparatively concentrated in the high-efficiency zones, except that of the CD-CS. It seems that the work points of the robust design are more concentrated compared to DP; however, this is mainly because the number of the working points is different. The MFAC-based EMS can lower the operating frequency of the engine in the low-efficiency region, and both the deterministic and robust designs of the MFAC-based EMS have fewer engine working points compared to DP. Thereby, the percentage of the engine working points in the high-efficiency zone is lower than DP.  As shown in Table 7, the percentages of engine operating points in the high-efficiency area with be ≤ 205g/kWh are calculated. The statistics are focused on the ratio of the working points within the efficient area to all operating points during the whole trip, so as to more intuitively reflect the FC of the engine. The results indicate that DP has the most energy-saving potential due to its globally optimal characteristic. The robust design of the MFAC-based strategy is slightly suboptimal compared to DP. Although the deterministic design may have an obvious advantage in contrast to the robust design in some cases (e.g., in NO.1), its stability for various driving cycles is weaker than the robust design. The robust design can be considered as the perfect approach for a real-time application compared with the deterministic design and CD-CS. The fuel consumption of the four strategies can be seen in Figure 22. The FCs of the DP are deployed as a benchmark to evaluate the performance of the other strategies. Note that the robust design can also perform successfully to ensure the energy-saving potential for different driving cycles, for which the FCs are more approximate to the DP solution.
Although the deterministic design can also have a better performance than CD-CS, its robustness still requires to be promoted to adapt to the variation of the driving cycles. The FCs of the robust design are, respectively, compared with other strategies, and the results are listed in Table 8. Overall, the FCs of the robust design for MFAC-based EMS are somewhat higher than DP whilst considerably lower compared to the deterministic design and CD-CS. The mean value of the promotion in FCs can reach up to 9.79% compared with the deterministic design, and a remarkable 19.66% compared withthe CD-CS. The most important thing is that it is only −1.68% less than DP for mean value of FCs for different driving cycles. This will be a significant advantage for a real-time application.

Conclusions
In this paper, a robust design framework based on DFSS is presented to strengthen the robustness of the MFAC-based EMS for real-time application, in which an SOC constraints method based on linear reference SOC is also incorporated to further plan the battery SOC. The robust control parameters of the MFAC-based EMS can be designed on the basis of the historical driving data, and the velocity and vehicle mass are considered as the external disturbance. Once the parameters are robustly designed, the MFAC-based EMS can have a satisfactory performance in stochastic future driving cycles. The results demonstrate that the SOC constraints have considerable advantages to limit the battery SOC, for which the actual battery SOC of the robust design is much closer to the DP. Additionally, the co-state of the PMP for robust design is more reasonable to decrease the fuel consumption in contrast to the deterministic design. Finally, the comparisons of the FCs also indicate that the robust design method has a remarkable FC decrease compared to the deterministic design, especially compared to the CD-CS strategy. The improvement can reach up to 9.79% and 19.66%, respectively. More importantly, the energy-saving potential is the same as that of DP, only −1.68% less than DP.
In the future, a test bench needs to be established to further access our recommended method, so as to be reliable for real-time application for PHEV. Moreover, an advanced electro-thermal-aging coupled battery model and the information based on V2I and V2V can also be integrated into our research to further promote the practicability of the MFACbased EMS.  Data Availability Statement: The data in this paper are from a real bus route and involve a confidentiality agreement. The dataset in this paper is not publicly available.

Conflicts of Interest:
The authors declare no conflicts of interest.