Development of Global Optimization Algorithm for Series-Parallel PHEV Energy Management Strategy Based on Radau Pseudospectral Knotting Method

Abstract

The powertrain model of the series-parallel plug-in hybrid electric vehicles (PHEVs) is more complicated, compared with series PHEVs and parallel PHEVs. Using the traditional dynamic programming (DP) algorithm or Pontryagin minimum principle (PMP) algorithm to solve the global-optimization-based energy management strategies of the series-parallel PHEVs is not ideal, as the solution time is too long or even impossible to solve. Chief engineers of hybrid system urgently require a handy tool to quickly solve global-optimization-based energy management strategies. Therefore, this paper proposed to use the Radau pseudospectral knotting method (RPKM) to solve the global-optimization-based energy management strategy of the series-parallel PHEVs to improve computational efficiency. Simulation results showed that compared with the DP algorithm, the global-optimization-based energy management strategy based on the RPKM improves the computational efficiency by 1806 times with a relative error of only 0.12%. On this basis, a bi-level nested component-sizing method combining the genetic algorithm and RPKM was developed. By applying the global-optimization-based energy management strategy based on RPKM to the actual development, the feasibility and superiority of RPKM applied to the global-optimization-based energy management strategy of the series-parallel PHEVs were further verified.

1. Introduction

Hybrid electric vehicles (HEVs) combine the advantages of traditional internal combustion engine vehicles and pure electric vehicles [1]. The advantages of HEVs, such as enabling long driving distances, being environment-friendly, and energy savings, make them one of the most practical ways to solve energy and pollution problems [2]. Plug-in hybrid electric vehicles (PHEVs) are improved versions of traditional HEVs with a long single-charge range. Among them, the series-parallel hybrid power system has the advantages of both series and parallel connections, making it capable of adapting to working conditions and reducing fuel consumption.

PHEVs have an inherent problem involving power distribution and mode switching because they have multiple power sources. Therefore, it is necessary to introduce a specific energy management strategy that can rationally distribute the energy flow of the power system. At present, energy management strategies of PHEVs can be divided into two categories: (1) rule-based strategy and (2) optimization-method-based strategy [3].

Rule-based energy management strategies mainly include deterministic rule-based strategies [4,5] and fuzzy logic strategies [6,7,8]. Optimization-method-based strategies can be divided into real-time optimization strategies and global optimization strategies [9,10]. Real-time optimization strategies mainly include the equivalent consumption minimum strategy (ECMS) [11] and model predictive control (MPC) strategies [12,13].

According to the optimal control theory, the sum of local optima is not necessarily the global optimum. For PHEV energy management strategies, it is usually expected to obtain the global optimal solution for fuel economy and emission under the whole driving cycle. The energy management strategy based on global optimization needs to predict all the information regarding driving conditions, and its calculated amount is large, which is difficult to perform in real time [14]. However, global optimization provides theoretical optimal control and is suitable for providing references and evaluation criteria for other energy management strategies [15,16]. The study of energy management strategies based on global optimization is conducive to understanding the shortcomings of rules and developing better real-time control rules [17,18]. At present, research on energy management strategies with global optimization is primarily focused on dynamic programming (DP) [19,20,21], stochastic dynamic programming (SDP) [22], evolutionary algorithms [23,24], and Pontryagin minimum principle (PMP) [25,26,27,28]. However, the PMP algorithm requires the Hamilton function to be continuously differentiable for state variables, which is not suitable for solving strongly nonlinear problems [29]. As the complexity of the model increases, the DP algorithm will face a serious dimension disaster problem, so that the calculation takes a very long time. Moreover, there is no theoretical guidance for choosing the weight coefficient of performance index in DP algorithm [30]. Many scholars have made relevant researches on how to reduce the computational burden of global-optimal-based energy management strategies to shorten the solution time [31,32,33]. The energy management strategies based on global optimization have become the focus of research institutions [34]

In the development of hybrid powertrains, it is necessary to consider both the global optimization of the energy management strategy and the global optimization of the powertrain parameters to ensure the best overall vehicle performance [35]. The energy management strategies based on global optimization have become the important supporting conditions for assembly matching of major automobile manufacturers [36,37]. Chief engineer of hybrid system need to observe the variation tendency between fuel consumption and the size of system’s component from more than 10,000 different components parameter schemes. In an increasingly competitive market, producers hope to shorten the information feedback cycle to shorten the development cycle. Therefore, the solution time of global-optimal-based energy management strategy should be shortened as much as possible. However, as mentioned above, the common global optimization algorithms, the DP algorithm, and the PMP algorithm were computationally expensive, which means these requirements are difficult to realize. Chief engineers of hybrid system are therefore in urgent need of a tool to quickly solve global-optimization-based energy management strategy.

To solve the problems of existing global optimization algorithms and meet the needs of chief engineers of hybrid system, a pseudospectral knotting-based method that solves the global optimization problem of PHEV energy management strategies to improve computational efficiency is proposed in this paper. A bi-level nested component-sizing method combining the genetic algorithm (GA) and Radau pseudospectral knotting method (RPKM) was used to determine the series-parallel PHEV power system scheme with the best fuel economy. It showed that the hybrid-powertrain optimization matching program combined with the RPKM can better meet the needs of the actual development process.

The remainder of this paper is organized as follows. Section 2 discusses modeling of the power demand and hybrid powertrains of PHEVs. Section 3 describes the optimal energy management problem in PHEV powertrains with the objective of minimizing fuel consumption. Section 4 describes the RPKM-based transcription from the optimal control problem to a non-linear programming (NLP) problem and the comparisons between DP and RPM. Section 5 presents the development of a bi-level nested component sizing strategy framework for optimal power management of series-parallel PHEVs running under the new European drive cycle (NEDC) to further demonstrate the high calculation efficiency of the RPKM algorithm. Conclusions are provided in Section 6.

2. Dynamic Model for the Series-Parallel PHEV

The PHEV considered in this study adopts a series-parallel hybrid power system configuration; its specific structure is shown in Figure 1. The configuration is mainly composed of the engine, integrated starter generator (ISG) motor, tracking motor (TM), and power battery. The power coupling mechanism is a simple gear transmission system.

By controlling the operating state of the tracking motor and the ISG motor and the release and combination of the clutch, the hybrid system shown in Figure 1 can achieve different working modes, including pure electric drive mode, series drive mode, pure engine drive mode, parallel drive mode, brake energy recovery mode, and parking charge mode. The series-parallel mode in which the clutch is combined and the ISG motor operates is not considered. Because the general parallel mode can meet all high load requirements, this mode cannot be used in actual operation. Moreover, in this mode, it is difficult for the ISG motor to work at a suitable working point; therefore, charging efficiency cannot be guaranteed.

The values of the key parameters of the hybrid system are given in

Table 1. Key parameters of series-parallel PHEV power system.

Parameter	Value
Complete quality of whole vehicle/kg	1645
Wheel radius/m	0.317
Gliding quality (sliding resistance parameter)/kg	1710
Sliding resistance curve/N	F = 0.0353 × u2 + 1.4114 × u + 114.7712
Efficiency of transmission system	0.965
Engine maximum output power/kW	65
Maximum power of ISG motor/kW	40
Maximum power of tracking motor/kW	110
Power cell capacity/(AŸH)	34
Speed ratio of engine to wheel end	3.7
Speed ratio of tracking motor to wheel end	7.5

.

2.1. Engine Model

In the global optimization study of PHEV energy management strategies, we only focused on the fuel consumption of the PHEV under specific target conditions, regardless of the complex physical and dynamic response processes inside the engine. Ignoring the dynamic characteristics of the engine and the effect of the heat exchange process on fuel consumption, the instantaneous fuel injection rate of the engine depends only on its output speed and output torque. According to the experimental data of a typical engine, an analytical model of the universal and external characteristics of the engine was established using the least squares method. To avoid the morbid phenomenon in the fitting process, the engine speed and torque were normalized before fitting. The instantaneous fuel injection rate after two-dimensional and five-polynomial fitting and normalization can be expressed as follows:

$$F_e(T_e,n_e)={\displaystyle \sum _{i=0}^5{\displaystyle \sum _{j=0}^i{A}_{1+j+{\displaystyle {\sum }_0^{i}n}}}}{(\frac{T_e}{T_{e\_\max }})}^{i-j}{(\frac{n_e}{n_{e\_\max }})}^j$$

(1)

where $F_e$ is the instantaneous fuel injection rate of the engine, g/s; $A_{1+\mathrm{j}+{{\displaystyle \sum }}_0^{\mathrm{i}}\mathrm{n}}$ is the fitting coefficient; $T_e$ is the output torque of the engine, $\mathrm{N}\cdot \mathrm{m}$; $T_{e\_\max }$ is the maximum output torque of the engine, $\mathrm{N}\cdot \mathrm{m}$; $n_e$ is the output speed of the engine, r/min; $n_{e\_\max }$ is the maximum output speed of the engine, r/min.

The R-square of the above polynomial fitting process is 0.9998, which is very close to 1, indicating that the fitting effect is very good. The normalized engine fuel consumption rate (g/kWh) is shown in Figure 2.

The external torque characteristic of the engine was fitted with a fifth degree polynomial, and the fitted polynomial is as follows:

$$T_{\max }(n_e)={\displaystyle \sum _{i=0}^5{B}_i(\frac{n_e}{n_{e\_\max }}}{)}^i$$

(2)

where $B_i$ is the fitting coefficient. The R-square of this process was 0.9991.

2.2. ISG Motor Model

As with engine modeling, the ISG motor was modeled using experimental modeling, focusing only on the input and output characteristics of the motor model, while ignoring the complex internal physical processes. The motor body, inverter, and motor controller were considered as a whole during the modeling process. Ignoring the dynamics and the effect of the heat exchange process on efficiency, the operating efficiency of an ISG motor depends only on its output speed and torque. After two-dimensional and five-polynomial fitting and normalization, the efficiency analysis formula for the ISG motors can be expressed as:

$${\eta }_{\mathrm{ISG}}(T_{\mathrm{ISG}},n_{\mathrm{ISG}})={\displaystyle \sum _{i=0}^5{\displaystyle \sum _{j=0}^i{C}_{1+j+{\displaystyle {\sum }_0^{i}n}}}}{(\frac{T_{\mathrm{ISG}}}{T_{\mathrm{ISG}\_\max }})}^{i-j}{(\frac{n_{\mathrm{ISG}}}{n_{\mathrm{ISG}\_\max }})}^j$$

(3)

where ${\eta }_{\mathrm{ISG}}$ is the working efficiency of the ISG motor; $C_{1+j+{{\displaystyle \sum }}_0^{i}n}$ is the fitting coefficient; $T_{\mathrm{ISG}}$ is the output torque of the ISG motor, $\mathrm{N}\cdot \mathrm{m}$; $T_{\mathrm{ISG}\_\max }$ is the maximum output torque of the ISG motor, $\mathrm{N}\cdot \mathrm{m}$; $n_{\mathrm{ISG}}$ is the output speed of the ISG motor, r/min; $n_{\mathrm{ISG}\_\max }$ is the maximum output speed of the ISG motor, r/min.

The R-square of the polynomial fitting process was 0.9932. The efficiency of the ISG motor after fitting and normalization is shown in Figure 3.

2.3. Tracking Motor Model

The tracking motor of the series-parallel PHEV can not only be used as an electric motor in the driving state, but also as a generator to recover the kinetic energy of the vehicle in the braking state. The modeling process was consistent with that of the ISG motor. After two-dimensional and five-polynomial fitting and normalization, the efficiency analysis formula for the tracking motor can be expressed as:

$${\eta }_{\mathrm{TM}}(T_{\mathrm{TM}},n_{\mathrm{TM}})={\displaystyle \sum _{i=0}^5{\displaystyle \sum _{j=0}^i{D}_{1+j+{\displaystyle {\sum }_0^{i}n}}}}{(\frac{T_{\mathrm{TM}}}{T_{\mathrm{TM}\_\max }})}^{i-j}{(\frac{n_{\mathrm{TM}}}{n_{\mathrm{TM}\_\max }})}^j$$

(4)

where ${\eta }_{\mathrm{TM}}$ is the working efficiency of the tracking motor; $D_{1+\mathrm{j}+{{\displaystyle \sum }}_0^{\mathrm{i}}\mathrm{n}}$ is the fitting coefficient; $T_{\mathrm{TM}}$ is the output torque for the tracking motor, $\mathrm{N}\cdot \mathrm{m}$; $T_{\mathrm{TM}\_\max }$ is the maximum output torque for the tracking motor, $\mathrm{N}\cdot \mathrm{m}$; $n_{\mathrm{TM}}$ is the output speed of the tracking motor, r/min; $n_{\mathrm{TM}\_\max }$ is the maximum output speed of the tracking motor, r/min.

The R-square of the polynomial fitting process was 0.9831. The efficiency of the tracking motor after fitting and normalization is shown in Figure 4.

2.4. Battery Model

The state of charge (SOC) is an important parameter of the battery pack. In this study, the influence of temperature was ignored; the internal resistance model based on charge-discharge test data was used to model the power battery [38,39]. The equivalent circuit of the power battery is shown in Figure 5.

In the figure above, $U_{oc}$ is the open circuit voltage (OCV) of the power battery, V; $U$ is the terminal voltage of the power battery; $R_{\mathrm{int}}$ is the internal resistance of the power battery, Ω; $I$ is the current of the power battery (discharge is positive, charging is negative), A.

According to the equivalent circuit diagram shown in Figure 5, the relationship between the open circuit voltage and the terminal voltage of the power battery can be expressed by the following equation:

$$U=U_{oc}(soc(t))-I(t)\times R_{\mathrm{int}}(soc(t))$$

(5)

The output power of power battery can be expressed as:

$$P_{bat}(t)=U\cdot I(t)$$

(6)

where $P_{bat}$ is the output power of the power battery (discharge is positive, charge is negative), kW.

The SOC of the power battery is calculated by the ampere-hour integration method [40]:

$$soc(t)=soc(t_0)-\frac{{\displaystyle \int I(t)dt}}{Q_{nom}}$$

(7)

where soc(t₀) represents the SOC of the power battery at the initial time; $Q_{nom}$ represents the nominal capacity of the power battery, As.

Combined with Equations (5)–(7), the first-order differential of the power battery SOC with time can be expressed as:

$$s\dot{o}c(t)=-\frac{U_{oc}(soc(t))-\sqrt{U_{oc}^2(soc(t))-4P_{bat}(t)R_{\mathrm{int}}(soc(t))}}{2Q_{nom}{R}_{\mathrm{int}}(soc(t))}$$

(8)

Based on the battery module, the relationships between the OCV and internal resistance and SOC are shown in Figure 6 and Figure 7, respectively.

2.5. Clutch Model

In the series-parallel hybrid system used in this study, the clutch was mainly used to switch between the series and parallel working modes. Therefore, in the backward simulation model, only the separated and combined states are considered; the dynamic characteristics of the clutch are ignored, and the separation and coupling process of the clutch is considered to be completed in an instant.

2.6. Longitudinal Dynamic Model of whole Vehicle

In the problem of optimal control of the energy management of the series-parallel PHEV, the fuel economy index of the whole vehicle is the main concern; therefore, the dynamic model of the whole vehicle only involves the longitudinal dynamic part. The running resistance includes rolling resistance, air resistance, slope resistance, and acceleration resistance, which can be expressed by the following formula:

$$F=mgf\cos \left(\arctan (i)\right)+mg\sin \left(\arctan (i)\right)+\frac{C_{D}A}{21.15}{u}_a^2+\delta m\frac{d{u}_a}{dt}$$

(9)

where F is the sum of all motion resistance forces, $\mathrm{N}$; m is the vehicle mass, kg; g is the gravity constant, $\mathrm{m}/{\mathrm{s}}^2$; f is the rolling resistance coefficient; i is the slope of the road; $C_D$ is the aerodynamic drag coefficient; A is the frontal area of the vehicle, ${\mathrm{m}}^2$; $u_a$ is the vehicle speed, $\mathrm{m}/\mathrm{s}$; $\delta$ is the rotational inertial coefficient. The torque demand of the driving wheel $T_{req}$ can be calculated as follows:

$$T_{req}=Fr$$

(10)

where r is the tire effective radius, m. The expression for torque balance of the whole vehicle is:

$$T_{req}/{\eta }_t=T_e\cdot i_e+T_{\mathrm{TM}}\cdot i_{\mathrm{TM}}+T_{brake}$$

(11)

where $T_{brake}$ is the output torque for mechanical brakes.

2.7. Summary

The backward simulation model of the series-parallel PHEV includes the above engine model, ISG motor, tracking motor, power battery, clutch, and vehicle longitudinal dynamics models. The establishment of the backward simulation model was completed by programming in Matlab. The simulation process is shown in Figure 8.

3. Optimal Control Problem of Series-Parallel PHEV Energy Management Strategy

Although the PHEV has an infinite number of combinations of control variables, there is always a method for power distribution for obtaining optimal fuel economy under certain target driving cycles. Therefore, the essence of the optimal control problem of PHEV energy management is finding the best combination of control sequences to minimize the performance index under the target driving cycle. The optimal control problem framework of series-parallel PHEV energy management strategy is shown in Figure 9.

3.1. Performance Index

When the series-parallel PHEV needs to maintain the charge-discharge balance of the power battery under certain working conditions, its performance index can be expressed as:

$$J={\displaystyle \sum _{k=1}^N[Q_{eng\_k}}+\alpha {(so{c}_k-so{c}_f)}^2]$$

(12)

where $J$ is the total fuel consumption under the target driving cycle, g; $Q_{eng\_k}$ represents the fuel consumption of the engine at stage k, g; $\alpha$ is the penalty factor; $so{c}_f$ is the expected SOC value at the terminal moment. By adjusting the penalty factor reasonably, the SOC of the power battery can be kept unchanged before and after the driving cycle.

3.2. State Equation

The series-parallel PHEV has multiple working modes, and the energy distribution modes of different working modes are also different. In the pure electric drive mode and pure engine drive mode, there is only one power source that provides energy. In the series drive mode, power is distributed between the auxiliary power unit (APU) and power battery. In the parallel drive mode, torque is distributed between the engine and tracking motor. In the braking energy recovery mode, torque distribution is performed between the tracking motor and mechanical brake. To clearly represent the subtle differences in the equations of state in different working modes, the distribution coefficient $\mu$ was chosen as the only control variable for the system. In addition, according to the structural characteristics of the series-parallel PHEV, the SOC of the power battery was selected as the only state variable of the system.

The state equation for a series-parallel PHEV can be written as:

$$so{c}_{k+1}=so{c}_k-\frac{U_{oc}(so{c}_k)-\sqrt{U_{oc}^2(so{c}_k)-4P_{bat\_k}{R}_{\mathrm{int}}(so{c}_k)}}{2Q_{nom}{R}_{\mathrm{int}}(so{c}_k)}$$

(13)

where $so{c}_k$ and $so{c}_{k+1}$ respectively represent the SOC of the power battery at moment $k$ and moment $k+1$; $U_{oc}(so{c}_k)$ represents the OCV of the power battery at moment $k$, V; $R_{\mathrm{int}}(so{c}_k)$ represents the internal resistance of the power battery at moment $k$, Ω; $P_{bat\_k}$ is the output power of the power battery at moment $k$, kW.

It should be pointed out that the difference in the state equations corresponding to different working modes is that the expression for the output power of the power battery is different.

In the pure electric drive mode and series drive mode (the clutch is released), the output power of the power battery can be expressed as:

$$\begin{array}{l}{P}_{bat\_k}=\frac{{\mu }_k{T}_{\mathrm{TM}\_k}{n}_{\mathrm{TM}\_k}}{9550{\eta }_{\mathrm{TM}\_k}(T_{\mathrm{TM}\_k},n_{\mathrm{TM}\_k})}\\ T_{\mathrm{TM}\_k}=\frac{T_{d\_k}}{i_{\mathrm{TM}}{\eta }_t}\end{array}$$

(14)

where ${\eta }_{\mathrm{TM}\_k}(T_{\mathrm{TM}\_k},n_{\mathrm{TM}\_k})$ represents the efficiency of the tracking motor at moment $k$; ${\mu }_k$ is the distribution coefficient at moment $k$, and its value range is ${\mu }_k\le 1$; $T_{\mathrm{TM}\_k}$ represents the output torque of the tracking motor at moment $k$.

When ${\mu }_k=1$, the required power of the vehicle is completely provided by the power battery, and the PHEV works in the pure electric drive mode. When ${\mu }_k<1$, PHEV works in the series drive mode, and the output power of the APU can be expressed as:

$$P_{\mathrm{APU}\_k}=\frac{(1-{\mu }_k)T_{\mathrm{TM}\_k}{n}_{\mathrm{TM}\_k}}{9550{\eta }_{\mathrm{TM}\_k}(T_{\mathrm{TM}\_k},n_{\mathrm{TM}\_k})}$$

(15)

In the series drive mode, when $0<{\mu }_k<1$, the PHEV works in the combined drive mode with the APU and power battery; when ${\mu }_k=0$, the PHEV works in the APU separate drive mode; when ${\mu }_k<0$, the PHEV works in the APU drive charging mode, where the output power of the power battery is negative.

In the pure engine drive mode and the parallel drive mode (the clutch is engaged), the output power of the power battery can be expressed as:

$$\begin{array}{l}{P}_{bat\_k}=\frac{T_{\mathrm{TM}\_k}{n}_{\mathrm{TM}\_k}}{9550{\eta }_{\mathrm{TM}\_k}(T_{\mathrm{TM}\_k},n_{\mathrm{TM}\_k})}\\ T_{\mathrm{TM}\_k}=\frac{{\mu }_k{T}_{d\_k}}{i_{\mathrm{TM}}{\eta }_t}\end{array}$$

(16)

The value range of the distribution coefficient is $0\le {\mu }_k<1$. When ${\mu }_k=0$, the PHEV works in the pure engine drive mode; when $0<{\mu }_k<1$, the PHEV works in parallel drive mode. In these two working modes, the output torque of the engine $T_{e\_k}$ at moment $k$ can be expressed as:

$$T_{e\_k}=\frac{(1-{\mu }_k)T_{d\_k}}{i_e{\eta }_t}$$

(17)

In the braking energy recovery mode (independent of the clutch state), the output power of the power battery is negative and can be expressed as:

$$\begin{array}{l}{P}_{bat\_k}=\frac{T_{\mathrm{TM}\_k}{n}_{\mathrm{TM}\_k}}{9550}\cdot {\eta }_{\mathrm{TM}\_k}(T_{\mathrm{TM}\_k},n_{\mathrm{TM}\_k})\\ T_{\mathrm{TM}\_k}=\frac{{\mu }_k{T}_{d\_k}}{i_{\mathrm{TM}}{\eta }_t}\end{array}$$

(18)

The value range of the distribution coefficient is $0\le {\mu }_k<0.3$. The reason why the distribution coefficient does not exceed 0.3 is that the highest efficiency of the PHEV in the braking energy recovery process is generally not more than 30%.

Table 2. The output power of the power battery at the different working modes. APU: auxiliary power unit.

Clutch States	Expression for the Output Power of the Power Battery	Mode	$\mathbf{Value}\mathbf{Range}\mathbf{of}{\mathit{\mu }}_{\mathit{k}}$
clutch is released	$\begin{array}{l}{P}_{bat\_k}=\frac{{\mu }_k{T}_{\mathrm{TM}\_k}{n}_{\mathrm{TM}\_k}}{9550{\eta }_{\mathrm{TM}\_k}(T_{\mathrm{TM}\_k},n_{\mathrm{TM}\_k})}\\ T_{\mathrm{TM}\_k}=\frac{T_{d\_k}}{i_{\mathrm{TM}}{\eta }_t}\end{array}$	pure electric drive mode	${\mu }_k=1$
		series drive mode	$0<{\mu }_k<1$
		APU separate drive mode	${\mu }_k=0$
		APU drive charging mode	${\mu }_k<0$
clutch is engaged	$\begin{array}{l}{P}_{bat\_k}=\frac{T_{\mathrm{TM}\_k}{n}_{\mathrm{TM}\_k}}{9550{\eta }_{\mathrm{TM}\_k}(T_{\mathrm{TM}\_k},n_{\mathrm{TM}\_k})}\\ T_{\mathrm{TM}\_k}=\frac{{\mu }_k{T}_{d\_k}}{i_{\mathrm{TM}}{\eta }_t}\end{array}$	parallel drive mode	$0<{\mu }_k<1$
		pure engine drive mode	${\mu }_k=0$
independent of clutch state	$\begin{array}{l}{P}_{bat\_k}=\frac{T_{\mathrm{TM}\_k}{n}_{\mathrm{TM}\_k}}{9550}\cdot {\eta }_{\mathrm{TM}\_k}(T_{\mathrm{TM}\_k},n_{\mathrm{TM}\_k})\\ T_{\mathrm{TM}\_k}=\frac{{\mu }_k{T}_{d\_k}}{i_{\mathrm{TM}}{\eta }_t}\end{array}$	braking energy recovery mode	$0\le {\mu }_k<0.3$

showed the expression for the output power of the power battery at the different working modes.

3.3. Constraints

In the solution process, limited by actual mechanical properties, each power component of the PHEV needs to meet the following constraints:

$$\{\begin{array}{l}{n}_{e\_\min }\le n_{e\_k}\le n_{\max }\\ T_{e\_\min }\le T_{e\_k}\le T_{e\_\max }\\ n_{\mathrm{ISG}\_\min }\le n_{\mathrm{ISG}\_k}\le n_{\mathrm{ISG}\_\max }\\ T_{\mathrm{ISG}\_\min }\le T_{\mathrm{ISG}\_k}\le T_{\mathrm{ISG}\_\max }\\ n_{\mathrm{TM}\_\min }\le n_{\mathrm{TM}\_k}\le n_{\mathrm{TM}\_\max }\\ T_{\mathrm{TM}\_\min }\le T_{\mathrm{TM}\_k}\le T_{\mathrm{TM}\_\max }\\ so{c}_{\min }\le so{c}_k\le so{c}_{\max }\\ P_{\mathrm{APU}\_\min }\le P_{\mathrm{APU}\_k}\le P_{\mathrm{APU}\_\max }\\ P_{b\_\min }\le P_{b\_k}\le P_{b\_\max }\end{array}$$

(19)

where $n_{e\_\min }$ and $n_{e\_max}$ represent the minimum and maximum speeds of the engine, respectively, r/min; $T_{e\_\min }$ and $T_{e\_max}$ represent the minimum and maximum torques of the engine, respectively, N·m; $n_{\mathrm{ISG}\_\min }$ and $n_{\mathrm{ISG}\_max}$ represent the minimum and maximum speeds of the ISG motor, respectively, r/min; $T_{\mathrm{ISG}\_\min }$ and $T_{\mathrm{ISG}\_max}$ represent the minimum and maximum torques of the ISG motor, respectively, N·m; $n_{\mathrm{TM}\_\min }$ and $n_{\mathrm{TM}\_max}$ represent the minimum and maximum speeds of the tracking motor, respectively, r/min; $T_{\mathrm{TM}\_\min }$ and $T_{\mathrm{TM}\_max}$ represent the minimum and maximum torques of the tracking motor, respectively, N·m; $so{c}_{\min }$ and $so{c}_{max}$ represent the minimum and maximum allowable SOC values of the power battery, respectively; $P_{\mathrm{APU}\_\min }$ and $P_{\mathrm{APU}\_max}$ represent the minimum and maximum output power of the APU, respectively, kW; $P_{b\_\min }$ and $P_{b\_max}$ represent the minimum and maximum output power of the power battery, respectively kW.

4. PHEV Global Optimization Algorithm Based on RPKM

The series-parallel PHEV has multiple working modes, and the state equations of different working modes are also different. Therefore, the energy management optimal control problem is essentially a nonlinear switching-type optimal control problem. For such optimal control problems, to ensure the accuracy and efficiency of the solution, staged processing is required. This section introduces the stage division principle of the Radau pseudospectral knotting method; then, on the basis of the optimal control problem discussed in Section 3, it describes the specific steps for solving the problem by using the Radau pseudospectral knotting method.

4.1. Stage Division Principle

According to the acceleration, the driving conditions of the car can be divided into acceleration, uniform speed, and deceleration (braking) conditions. In theory, the pure electric drive mode, pure engine drive mode, series drive mode, and parallel drive mode of the series-parallel PHEV can work in both accelerated and constant speed conditions. Therefore, its working modes in accelerated and uniform working conditions have infinite combinations. However, in practical engineering applications, to balance the requirements of comfort and life of mode-switching components such as clutches, the working mode of the PHEV does not normally switch frequently (just as conventional internal combustion engines do not frequently shift gears). In addition, the switching time between different working modes in the PHEV energy management optimal control problem is often unknown, and too many switching time points to be optimized will reduce the efficiency of the Radau pseudospectral knotting method. Therefore, before introducing the RPKM, the typical working condition is shown in Figure 10 as an example to determine the stage division principle.

The mode division rules were set such that only when the speed was above 30 km/h can the engine start to work. In addition, only the pure electric drive mode and parallel drive mode (including pure engine driving mode) were considered under accelerated working conditions; only the pure electric drive mode and series drive mode were considered under uniform speed conditions. Under deceleration (braking) conditions, only the braking energy recovery mode was considered.

4.2. Numerical Solution of RPKM

The solution principle of RPKM is as follows: firstly, converting the time domain of Bolza-type optimal control problem to Legendre orthogonal multi-interval. Secondly, the Lagrange global interpolation polynomials are used to approximate control variables and state variables at Legendre orthogonal multi-intervals. Then, the Gauss-Radau numerical integral and the differential matrix are used to transform the integral type performance function and the differential type state equation, respectively, to approximate the derivative of the state variable to time in the dynamic equation. Based on this, calculus operations can be transformed into algebraic operations of interpolation polynomials. Finally, the original Bolza-type optimal control problem is transformed into the NLP problem with the discrete points of the control variables and state variables as the parameters to be optimized. The NLP problem can be solved using mature large-scale nonlinear optimization software (such as Sparse Nonlinear OPTimizer (SNOPT)). The optimal solution for the global optimization problem can be obtained.

According to the solution principle of RPKM, assuming that the number of stages of the PHEV in a driving cycle is Q. The initial time, the time at the $Q-1$ segmentation point and the terminal time were recorded as $T_0,T_1,\cdot \cdot \cdot ,T_Q$, and the time point satisfied $T_0<T_1<\cdot \cdot \cdot <T_Q$. The solution process of RPKM includes the following steps.

(1) Time domain conversion

Time intervals for each stage $[T_{q-1},T_q]$($q\in [1,Q]\cap Z^{+}$) are converted to the defined interval of the Legendre orthogonal polynomial:

$$\tau =\frac{2t-(T_q+T_{q-1})}{T_q-T_{q-1}},\tau \in [-1,1]$$

(20)

(2) Collocation and discretization

The Legendre–Gauss–Radau (LGR) point used by RPKM can be valued in the half-open interval $\tau \in (-1,1]$ or $\tau \in [-1,1)$. The two collocation methods are equivalent in solving the problem. In this study, we take the LGR collocation of interval $\tau \in (-1,1]$, which is the root of polynomial $P_N(\tau )-P_{N-1}(\tau )$, where $P_N(\tau )$ is the N-order Legendre orthogonal polynomial:

$$P_N(\tau )=\frac{1}{2^{N}N!}\frac{d^N}{d{\tau }^N}[{({\tau }^2-1)}^N],N=0,1,2,\cdot \cdot \cdot$$

(21)

The time history curves of state variables and control variables in each stage can be approximated using Lagrange interpolation polynomials:

$$so{c}_q(\tau )\approx {\displaystyle \sum _{i=0}^{N_q}{L}_{q,i}(\tau )\cdot so{c}_{q,i}}$$

(22)

$${\mu }_q(\tau )\approx {\displaystyle \sum _{i=1}^{N_q}{\tilde{L}}_{q,i}(\tau )\cdot {\mu }_{q,i}}$$

(23)

where $L_{q,i}(\tau )$ and ${\tilde{L}}_{q,i}\left(\tau \right)$ are Lagrangian interpolation basis functions, defined as follows:

$$L_{q,i}(\tau )={\displaystyle \prod _{j=0,j\ne i}^{N_q}\frac{\tau -{\tau }_{q,j}}{{\tau }_{q,i}-{\tau }_{q,j}}}$$

(24)

$${\tilde{L}}_{q,i}(\tau )={\displaystyle \prod _{j=1,j\ne i}^{N_q}\frac{\tau -{\tau }_{q,j}}{{\tau }_{q,i}-{\tau }_{q,j}}}$$

(25)

(3) Dynamic equation constraint transformation into algebraic constraints

After collocation and discretization, the state variable was approximated using the global interpolation polynomial; thus, its derivative could be expressed using the derivative of the interpolation polynomial (22):

$$s\dot{o}{c}_q({\tau }_{q,k})={\displaystyle \sum _{i=0}^{N_q}{\dot{L}}_{q,i}}({\tau }_{q,k})\cdot so{c}_{q,i}={\displaystyle \sum _{i=0}^{N_q}{D}_{ki}^q}\cdot so{c}_{q,i}$$

(26)

where ${\tau }_{q,k}$ represents the collocation point of stage q, with $k=1,2,\cdots ,N_q$; $D^q$ is the differential matrix of order $N_q\times (N_q+1)$, representing the differential values of the Lagrangian interpolation basis function at each LGR collocation point in stage q, which can be expressed by the following formula:

$$D_{ki}^q={\dot{L}}_i({\tau }_k)=\{\begin{array}{ll}\frac{g({\tau }_{q,k})}{({\tau }_{q,k}-{\tau }_{q,i})g({\tau }_{q,i})},& k\ne i\\ \frac{g({\tau }_{q,i})}{2g({\tau }_{q,i})},& k=i\end{array}$$

(27)

where $g({\tau }_{q,i})=(1+{\tau }_{q,i})[P_{N_q}({\tau }_{q,i})-P_{N_q-1}({\tau }_{q,i})]$.

Combined with Equations (8), (26), and (27), the dynamic equation constraints of the serial-parallel PHEV energy management optimal control problem in stage q can be transformed into algebraic constraints at collocation point ${\tau }_{q,k}$:

$${\displaystyle \sum _{i=0}^{N_q}{D}_{ki}^q}\cdot so{c}_{q,i}=\frac{T_q-T_{q-1}}{2}\frac{\sqrt{U_{oc}^2(so{c}_{q,k})-4P_{bat}({\tau }_{q,k})R_{\mathrm{int}}(so{c}_{q,k})}-U_{oc}(so{c}_{q,k})}{2Q_{nom}{R}_{\mathrm{int}}(so{c}_{q,k})}$$

(28)

where $P_{bat}({\tau }_{q,k})$ represents the output power value of the power battery at collocation point ${\tau }_{q,k}$.

The expressions for the output power of the power battery corresponding to each working mode can be obtained from Equations (14)–(18).

(4) Performance index transformation

In the optimal control problem of series-parallel PHEV energy management strategy, the performance index includes two parts, namely engine fuel consumption and power battery power consumption, as shown in Equation (29), which corresponds to the Bolza-type performance index:

$$J=\beta (so{c}_0-so{c}_f)+{\displaystyle {\int }_0^{t_f}{F}_e(T_e},n_e)dt$$

(29)

where $\beta$ is the oil-to-electricity conversion coefficient, which converts electricity consumed by the power battery into the equivalent engine fuel consumption; $so{c}_0$ and $so{c}_f$ represent the initial and terminal SOC values of the power battery, respectively; $F_e(T_e,n_e)$ is the instantaneous fuel injection rate of the engine, g/s.

For stage q, the form of the performance index depends on the working mode of the PHEV. For the pure electric drive mode and braking energy recovery mode, as the engine does not participate in the work, the performance index at this time only has a nonintegral term (Mayer-type performance index), expressed as follows:

$$J_q=\beta (so{c}_{q,0}-so{c}_{q,N_q})$$

(30)

For the pure engine drive mode and parallel drive mode, as the engine always produces power, its performance index contains an integral term (Lagrange performance index), which can be approximated using the Gauss-Radau integral method:

$$J_q=\beta (so{c}_{q,0}-so{c}_{q,N_q})+\frac{T_q-T_{q-1}}{2}{\displaystyle \sum _{k=1}^{N_q}{\omega }_{q,k}}{F}_e(T_e({\tau }_{q,k}),n_e({\tau }_{q,k}))$$

(31)

where ${\omega }_{q,k}$ is the integral weight, defined as:

$${\omega }_{q,k}=\frac{1-{\tau }_{q,k}}{N_q^2\cdot {[P_{N_q-1}({\tau }_{q,k})]}^2}$$

(32)

For the series drive mode, the APU can always operate at the optimum efficiency point owing to decoupling of the APU from the driving conditions. Considering the working efficiency of the engine and ISG motor, the optimal fuel consumption rate of the APU can be expressed as a function of the output power [41]:

$$F_{\mathrm{APU}}(P_{\mathrm{APU}})={\displaystyle \sum _{i=0}^2{E}_i}{P}_{\mathrm{APU}}^i$$

(33)

where $F_{APU}$ is the optimal fuel consumption rate of the APU under different powers, g/s; $E_i$ is the fitting coefficient; $P_{APU}$ is the output power of the APU, kW. The performance index in the series drive mode can be expressed as:

$$J_q=\beta (so{c}_{q,0}-so{c}_{q,N_q})+\frac{T_q-T_{q-1}}{2}{\displaystyle \sum _{k=1}^{N_q}{\omega }_{q,k}}{F}_{\mathrm{APU}}(P_{\mathrm{APU}}({\tau }_{q,k}))$$

(34)

(5) NLP problem

To facilitate the unification of state equations and performance index of different working modes, $m1$, $m2$, $m3$, $m4$, and $m5$ were used to represent the pure electric drive mode, pure engine drive mode, series drive mode, parallel drive mode, and braking energy recovery mode, respectively. After performing the above steps, the optimal control problem of serial-parallel PHEV energy management was transformed into the NLP problem as shown below.

$$\begin{array}{l}J=\beta (so{c}_{1,0}-so{c}_{Q,N_Q})+{\displaystyle \sum _{q=1}^Q\frac{T_q-T_{q-1}}{2}}{\displaystyle \sum _{k=1}^{N_q}{\omega }_{q,k}}F\\ F=\{\begin{array}{ll}0,& q=m1,m5\\ F_e(T_e({\tau }_{q,k}),n_e({\tau }_{q,k})),& q=m2,m4\\ F_{\text{APU}}(P_{\text{APU}}({\tau }_{q,k})),& q=m3\end{array}\\ s.t.\\ {\displaystyle \sum _{i=0}^{N_q}{D}_{ki}^q}\cdot so{c}_{q,i}=\frac{T_q-T_{q-1}}{2}\frac{\sqrt{U_{oc}^2(so{c}_{q,k})-4P_{bat}({\tau }_{q,k})R_{\mathrm{int}}(so{c}_{q,k})}-U_{oc}(so{c}_{q,k})}{2Q_{nom}{R}_{\mathrm{int}}(so{c}_{q,k})}\\ m1,m3:\{\begin{array}{l}{P}_{bat}({\tau }_{q,k})=\frac{{\mu }_{q,k}{T}_d({\tau }_{q,k})n_{\text{TM}}({\tau }_{q,k})}{9550i_{\text{TM}}{\eta }_t\cdot {\eta }_{\text{TM}}(T_d({\tau }_{q,k})/(i_{\text{TM}}{\eta }_t),n_{\text{TM}}({\tau }_{q,k}))}\\ P_{\text{APU}}({\tau }_{q,k})=\frac{(1-{\mu }_{q,k})T_d({\tau }_{q,k})n_{\text{TM}}({\tau }_{q,k})}{9550i_{\text{TM}}{\eta }_t\cdot {\eta }_{\text{TM}}(T_d({\tau }_{q,k})/(i_{\text{TM}}{\eta }_t),n_{\text{TM}}({\tau }_{q,k}))}\end{array}\\ m2,m4:\{\begin{array}{l}{P}_{bat}({\tau }_{q,k})=\frac{{\mu }_{q,k}{T}_d({\tau }_{q,k})n_{\text{TM}}({\tau }_{q,k})}{9550i_{\text{TM}}{\eta }_t\cdot {\eta }_{\text{TM}}({\mu }_{q,k}{T}_d({\tau }_{q,k})/(i_{\text{TM}}{\eta }_t),n_{\text{TM}}({\tau }_{q,k}))}\\ T_e({\tau }_{q,k})=\frac{(1-{\mu }_{q,k})T_d({\tau }_{q,k})}{i_e{\eta }_t}\end{array}\\ m5:P_{bat}({\tau }_{q,k})=\frac{{\mu }_{q,k}{T}_d({\tau }_{q,k})n_{\text{TM}}({\tau }_{q,k})}{9550i_{\text{TM}}{\eta }_t}\cdot {\eta }_{\text{TM}}({\mu }_{q,k}{T}_d({\tau }_{q,k})/(i_{\text{TM}}{\eta }_t),n_{\text{TM}}({\tau }_{q,k}))\\ T_d({\tau }_{q,k})=[mgf+mg\cdot i({\tau }_{q,k})+\frac{C_{D}A}{21.15}{u}^2({\tau }_{q,k})+\delta m{a}^2({\tau }_{q,k})]\cdot r\\ n_{\text{TM}}({\tau }_{q,k})=\frac{i_{\text{TM}}\cdot u({\tau }_{q,k})}{0.377r}\\ n_e({\tau }_{q,k})=\frac{i_e\cdot u({\tau }_{q,k})}{0.377r}\\ C(so{c}_{q,k},{\mu }_{q,k})\le 0\end{array}$$

(35)

In Equation (35), $q=1,2,\cdots \dots ,Q$; $k=1,2,\dots ,N_q$; F is the instantaneous fuel injection rate, g/s. To simplify the expression in Equation (35), $C(so{c}_{q,k},{\mu }_{q,k})\le 0$ was used to represent the path constraint of the NLP problem. The detailed expression is given in Equation (19). The optimized variables of the NLP problem are the SOC value at the LGR points, allocation coefficient, and switching time points between different modes, which can be solved using the mature high-dimensional sparse NLP solver SNOPT. Its core algorithm is the sequence quadratic programming method with local superlinear convergence based on line search.

4.3. Result Analysis

To verify the feasibility of applying the Radau pseudospectral method to the global optimization of the series-parallel PHEV energy management strategy, we used the DP algorithm and RPKM to solve the optimal control problem of energy management under the typical working conditions shown in Figure 11. The working conditions comprise acceleration, uniform speed, and deceleration. The total time was 146 s and the maximum speed was 90 km/h. Figure 12 showed the power of each component of the dynamic system of serial-parallel PHEV under this working condition.

It can be seen from

Table 3. Comparison of optimization results between DP and RPKM.

Comparison Items		DP	RPKM
Before post-processing	Fuel consumption/g	107.99	109.91
	100 km fuel consumption/(L/100 km)	5.97	6.08
	Relative error of two algorithms for 100 km fuel consumption/(%)	1.81	—
After post-processing	Fuel consumption/g	109.32	109.44
	100 km fuel consumption/(L/100 km)	6.047	6.054
	Relative error of two algorithms for 100 km fuel consumption/(%)	0.12	—
Solving time/s		32477.19	17.98

that the fuel consumptions of the two global optimization algorithms are close, with a relative error of only 1.81%. It is evident that the two algorithms have very close global optimization capabilities; therefore, the applicability of the RPKM was verified to optimize the energy management strategy of the series-parallel PHEV.

It is also clear from Table 3 that when solving the same optimal control problem, the DP algorithm takes 32,477 s (9.06 h), while the RPKM only needs 17.98 s. The calculation efficiency of the RPKM is approximately 1806 times higher than that of the DP algorithm, which shows that the RPKM is considerably fast in solving the optimal control problem in series-parallel PHEV energy management.

RPKM treats the problem as a continuous process in each sub-stage and searches for the optimal solution through optimal gradient descent. Compared with DP algorithm, which discretizes all the state variables and control variables of the global optimization problem, and then searches for the optimal solution through a method similar to exhaustive method, RPKM avoids the “dimension disaster” problem; thus, its solving efficiency is higher.

4.4. Error Analysis

The solution errors of RPKM mainly include the approximate errors of state variables and control variables, the differential errors caused by the derivation of Lagrange global interpolation polynomials, and the integral errors in performance indexes. Since the truncation error of Lagrangian global interpolation polynomial is n + 1 times and the algebraic accuracy of Gauss numerical integral is 2n + 1 times, the pseudo-spectral method has a high computational accuracy.

The solution accuracy of DP is greatly affected by factors such as time step (number of stages), the discrete step size of state variables and control variables, and the penalty factors. When the DP algorithm was used to solve the instantaneous fuel consumption corresponding to different state variables at each stage, the instantaneous oil injection rate under the current state is multiplied by the time step to approximate the integral of instantaneous oil injection rate with time, so the integral error of this process was greatly affected by the time step. Theoretically, when the discrete step size of the state variables and the control variables are fixed, the smaller the time step of DP algorithm is, the better the solution accuracy is. However, as the time step decreases, the solving time of DP algorithm will increase exponentially.

Therefore, the following methods are proposed to reduce the integral error in the DP algorithm: the optimal control variable sequence obtained by the DP algorithm is re-substituted into the original optimal control problem. Then the new fuel consumption value was obtained by the forward solution in a time step of 10⁻⁶ s. This post-processing method can not only avoid too long solving time, but also ensures the solution accuracy of DP algorithm.

As shown in Table 2, after the same treatment of RPKM, the relative error decreases from 1.81% to 0.12%, indicating that the two algorithms have very close global optimization ability.

5. Bi-Level Nested Component-sizing Method Based on GA and RPKM

In the process of parameter sizing of the PHEV dynamic system, the influence of power system parameters and energy management strategies on vehicle performance is coupled. The parameter optimization of the dynamic system and global optimization of the energy management strategy cannot be considered separately. Therefore, a bi-level nested parameter optimization method based on the GA and RPKM was developed to optimize the final dynamic system of the PHEV.

5.1. Description of the Component-Sizing Problem

The component-sizing problem of the series-parallel PHEV power system is to find the best power system scheme with the best fuel economy while satisfying the power performance index of the vehicle. This component-sizing problem can be expressed using the following mathematical model:

$$\begin{array}{l}\min f(X)\\ s.t.g_j(X)\ge 0,j=1,2,\cdots ,m\\ X=[x_1,x_2,\cdots ,x_n]\\ x_i^{low}\le x_i\le x_i^{up},i=1,2,\cdots ,n\end{array}$$

(36)

where $f(X)$ is the optimization objective; $X$ is an n-dimensional vector, representing a set of feasible solutions to the above optimization problem; $g_j(X)\ge 0$ is the constraint representing the dynamic performance of the vehicle; $m$ is the number of constraints; $n$ is the number of parameters to be optimized; $x_i^{low}$ and $x_i^{up}$, respectively, represent the lower and upper limits of the i-th parameter to be optimized, $x_i$.

The goal of power system parameter optimization is to achieve the best energy consumption economy. Energy consumption in the optimization process includes the power consumption of the battery and fuel consumption of the engine. The power consumption of the battery and fuel consumption of the engine were converted into corresponding operating costs (where we assumed 1 kWh as equivalent to 300 g of gasoline). The operating cost of the driving cycle was taken as the optimization objective $f(X)$ of component sizing of the power system.

The dynamic performance and fuel economy of series-parallel PHEVs are influenced by the engine, ISG motor, tracking motor, power battery, and transmission system. To obtain the best vehicle performance, all relevant parameters need to be optimized. Therefore, the maximum output power of the engine, peak power of the ISG motor, peak power of the tracking motor, total capacity of the power battery, speed ratio of the tracking motor to the wheel, and speed ratio of the engine to the wheel were selected as the optimization variables in this study.

The test for determining the energy consumption of the PHEV was carried out under the NEDC, as shown in Figure 13.

Here, the CD–CS energy management strategy was adopted. The CD strategy was adopted in the urban working condition, in which the tracking motor was only provided with the demand power of the vehicle. The CS strategy was adopted in the suburban working condition, in which the SOC of the battery remained the same before and after applying the working condition. In this condition, the mode division was also based on the rules given in Section 4.1.

5.2. Solution for Bi-Level Nested Component-Sizing Method

The GA is a global optimization random search algorithm simulating biological inheritance and evolution processes. It has proved to be an effective strategy to solve complex engineering optimization problems characterized by nonlinear, multimodal, and nonconvex objective functions.

The specific optimization steps of this bi-level nested method for PHEV component sizing are shown in Figure 14.

• The optimization variables of the PHEV parameter optimization problem were converted into individuals in the genetic space, and then a group of initial populations was randomly generated in the feasible domain of the problem.
• The power components and transmission parameters represented by each individual in the population were replaced in the series-parallel PHEV backward simulation model, and the consumption of electricity in the urban driving cycle of NEDC was calculated. The optimal fuel consumption under the suburban driving cycle of NEDC was determined using the RPKM. The electricity consumption and the optimal fuel consumption were converted to the corresponding transportation costs to determine the suitability of each individual (the total operating cost of NEDC).
• The next generation of the population was selected according to the size of the fitness. The smaller the adaptability, the larger the probability of being selected; the individuals with large adaptability were eliminated in the process of evolution.
• Crossover and mutation operations of individuals in the new species group produced the next generation of individuals.
• Steps 2 to 4 were repeated until the termination condition was satisfied. At this time, the algorithm converged to the best chromosome, obtaining the optimal solution to the problem.

In this method, the backward simulation model of the PHEV was embedded into the GA program as a module. The inner layer calculated its adaptability according to the component parameters of the individual; the outer layer used the fitness of each individual according to the inner layer to continuously search and optimize parameters in the multidimensional parameter space, guaranteeing the global optimality of component-sizing results.

5.3. Component-Sizing Results

Figure 15 was developed based on the relationship between optimal fitness and iteration number after 50 iterations based on the GA and RPKM. When the fitness is small and remains unchanged, the optimal value of GA is considered to be the global optimal solution. It can be seen from Figure 15 that, when the iteration reached 20 generations, the fitness tended to be stable. At this point, we can assume that the parameters obtained by GA are the global optimal solution, that is, the optimal solution of component-sizing. This process took 16.7 h; however, if it were based on the GA and DP algorithm, it would require approximately hundreds of hours.

The result of parameter optimization presented in

Table 4. Parameter optimization results.

Optimization Items	Values
$P_{e\max }$/kW	69.99
$P_{\mathrm{ISG}\max }$/kW	65.14
$P_{\mathrm{TM}\max }$/kW	120.38
$Q_{nom}$/A·h	39.46
$i_m$	7.65
$i_e$	3.49

can provide a reference basis for the selection of power system schemes for trial vehicles.

6. Conclusions

To solve the problem that the traditional DP algorithm and PMP algorithm have a long solution time or even cannot solve due to the complexity of the series-parallel PHEVs powertrain model, this paper proposes to use the Radau pseudospectral knotting method (RPKM) to solve the global-optimization-based energy management strategy of the series-parallel PHEVs to improve computational efficiency. A backward simulation model of the series-parallel PHEVs for optimal control was established, and the optimal control problem of energy management strategy was established. The optimal control problem was transformed into an NLP problem by using the RPKM. Simulation results showed that compared with DP algorithm, the global-optimization-based energy management strategy based on the RPKM improves the computational efficiency by 1806 times with a relative error of only 0.12%. Based on the high calculation efficiency of RPKM, a bi-level nested component-sizing method combining the GA and RPKM was developed. The results indicated that RPKM has lower computational cost and can provide optimization results for chief engineers of hybrid system more efficiently. Through the above research, the feasibility and superiority of RPKM applied to the global-optimization-based energy management strategy of the series-parallel PHEVs were verified.