Development of Driving Cycle Construction for Hybrid Electric Bus: A Case Study in Zhengzhou, China

A driving cycle is important to accomplish an accurate depiction of a vehicle’s driving characteristics as the traction motor’s flexible response to stop and start commands. In this paper, the driving cycle construction of an urban hybrid electric bus (HEB) in Zhengzhou, China is developed in which a measurement system integrating global positioning and inertial navigation function is used to acquire driving data. The collected data are then divided into acceleration, deceleration, uniform, and stop fragments. Meanwhile, the velocity fragments are classified into seven state clusters according to their average velocities. A transfer matrix applied to reveal the transfer relationship of velocity clusters can be obtained with statistical analysis. In the third stage, a three-part construction method of driving cycle is designed. Firstly, according to the theory of Markov chain, all the alternative parts that satisfy the construction’s precondition are selected based on the transfer matrix and Monte Carlo method. The Zhengzhou urban driving cycle (ZZUDC) could be determined by comparing the performance measure (PM) values subsequently. Eventually, the method and the cycle are validated by the high correlation coefficient (0.9972) with original data of ZZUDC than that of the other driving cycle (0.9746) constructed with traditional micro-trip and as well by comparing several statistical characteristics of ZZUDC and seven international cycles. Particularly, with around 20.5 L/100 km fuel and approximately 12.8 kwh/100 km electricity consumption, there is a narrow gap between the energy consumption of ZZUDC and WVUCITY, and their characteristics are similar.


Introduction
Recently, hybrid electric buses (HEB) with better energy consumption and environmental protection performance compared with traditional fuel engine buses have been adopted in many Chinese cities and become a hot research topic [1].
As an evaluation standard and a detection criterion of vehicle emissions and fuel consumption, driving cycle is an important design basis for vehicle power matching and fuel economy optimization [2,3]. Driving cycle is a velocity-time trace that describes driving characteristics of specific vehicles in specific traffic environments (such as urban road or expressway). Usually, it is a representative running condition which is extracted from a large amount of measured driving data using mathematical and statistical methods. Different regions have different driving characteristics; hence, developing a typical driving cycle for each place according to its situation is essential [4]. The current dominant methods of driving cycle construction include the micro-trip method, fixed step analysis method, velocity-acceleration analysis method, wavelet analysis method, and Markov analysis method [5][6][7][8]. Many scholars have investigated the development of urban driving cycles

The HEB Original Driving Data Acquisition
The HEB original driving data can reflect the urban traffic conditions and HEB driving characteristics, which is the source of data for driving cycle construction. The more details of the original driving data used in the construction process, the more objective and representative the driving cycle's characteristics would be. Two representative HEB lines of Zhengzhou City in China are selected as the experimental objects of the driving data acquisition test after thorough consideration based on the test route, choosing rule proposed below. Note that we assume that HEB can work in pure electric driving mode at low velocity.

Test Route Selection
The rules for choosing the test route are as follows: 1. Hybrid electric bus is the only choice to drive on the test route.
2. The test route should include both traditional and rapid bus routes.
Bus rapid transit (BRT) is a new type of public transport system running between the railway and traditional bus transportation. BRT uses the special bus route to achieve rail transit mode operation, and it has been applied widely in China. The driving states of buses in the BRT routes are totally different from those in the traditional bus routes. Therefore, the construction of an urban driving cycle which fuses traditional bus route and rapid bus route is consistent with the actual operation situation of Chinese urban buses. 3. The selected bus route should cover congested and fluent urban regions.

The HEB Original Driving Data Acquisition
The HEB original driving data can reflect the urban traffic conditions and HEB driving characteristics, which is the source of data for driving cycle construction. The more details of the original driving data used in the construction process, the more objective and representative the driving cycle's characteristics would be. Two representative HEB lines of Zhengzhou City in China are selected as the experimental objects of the driving data acquisition test after thorough consideration based on the test route, choosing rule proposed below. Note that we assume that HEB can work in pure electric driving mode at low velocity.

Test Route Selection
The rules for choosing the test route are as follows:

1.
Hybrid electric bus is the only choice to drive on the test route.

2.
The test route should include both traditional and rapid bus routes. Bus rapid transit (BRT) is a new type of public transport system running between the railway and traditional bus transportation. BRT uses the special bus route to achieve rail transit mode operation, and it has been applied widely in China. The driving states of buses in the BRT routes are totally different from those in the traditional bus routes. Therefore, the construction of an urban driving cycle which fuses traditional bus route and rapid bus route is consistent with the actual operation situation of Chinese urban buses. 3.
The selected bus route should cover congested and fluent urban regions.
According to the above three rules, the BRT route B17 and traditional bus route 906 are selected as the test routes, which are shown in Figure 2 and Table 1.  As can be seen in Figure 2, BRT route B17 runs from east to west, passing through crowded areas including railway stations and commercial centers, as well as residential areas, factories, schools, parks, and urban landmark scenic spots. Meanwhile, bus route 906 runs from south to north, passing through bus stations, universities, government offices, and other important urban facilities. The two lines cover all levels of roads in the urban area of Zhengzhou and account for roughly the same proportion for all levels of roads.

Test Bus Selection
The HEB that drives on bus route 906 is the same as that on BRT route B17. Thus, there is no difference between buses in the process of the fusing of the aforementioned two bus routes for urban driving cycle construction. The specific parameters of the test HEB are shown in Table 2.

The Driving Data Acquisition System
A global positioning and inertial navigation system (GPS/INS integrated system) named OXTS inertial+ is used for the HEB driving data acquisition. The measurement parameters include velocity, transient acceleration, and road slope. The location precision of measurement is 2 cm 1σ, and velocity  As can be seen in Figure 2, BRT route B17 runs from east to west, passing through crowded areas including railway stations and commercial centers, as well as residential areas, factories, schools, parks, and urban landmark scenic spots. Meanwhile, bus route 906 runs from south to north, passing through bus stations, universities, government offices, and other important urban facilities. The two lines cover all levels of roads in the urban area of Zhengzhou and account for roughly the same proportion for all levels of roads.

Test Bus Selection
The HEB that drives on bus route 906 is the same as that on BRT route B17. Thus, there is no difference between buses in the process of the fusing of the aforementioned two bus routes for urban driving cycle construction. The specific parameters of the test HEB are shown in Table 2.

The Driving Data Acquisition System
A global positioning and inertial navigation system (GPS/INS integrated system) named OXTS inertial+ is used for the HEB driving data acquisition. The measurement parameters include velocity, transient acceleration, and road slope. The location precision of measurement is 2 cm 1σ, and velocity precision is 0.05 km/hRMS. As for the acceleration precision, the error is 10 mm/s 2 1σ, linearity is 0.01%, and the total range is 100 m/s 2 .

Velocity and Acceleration Measurement
Compared with the traditional method, in which the driving data are read by the anti-lock braking system (ABS) speed signal, the GPS receiver can not only measure bus velocity with higher precision but also acquire the bus altitude and location for road slope estimation. The bus acceleration information is measured by the inertial navigation system, and the measurement error will not increase with the enhancement of sampling frequency. In order to retain more details of velocity fragments, the test sampling frequency is set to 5 Hz.

Road Slope Estimation
As indicated, other than velocity, road slope is also a particularly important characteristic for HEBs in control and design. This paper estimates road slope by calculating the ratio for the altitude difference between two sampling times to the horizontal distance traveled by the HEB during the sampling period and then taking the arctangent, which follows the equation where ∆h is the altitude change between the previous and current sampling time, v x is the bus horizontal velocity, T is the sampling period. The road slope can be calculated with Equation (1); however, the computed preliminary results have high magnitude and some of the points obviously violate the longitude slope design requirements of the China code for design of urban road engineering [24]. Hence, the road slope preliminary values need to be filtered, the steps of which are as follows:

1.
Primarily, eliminate slope calculated values with low velocity to avoid the denominator of Equation (1) being close to 0.

2.
Carry out filtering. If the calculated slope values are higher than 9% [25], then replace them with those of the previous sampling time. 3.
Perform average filtering against the processed slope values to smooth out the sawtooth signals.
The road slope after filtering is partially shown in Figure 3.
Sustainability 2020, 12, x FOR PEER REVIEW 5 of 18 precision is 0.05 km/hRMS. As for the acceleration precision, the error is 10 mm/s² 1σ, linearity is 0.01%, and the total range is 100 m/s².

Velocity and Acceleration Measurement
Compared with the traditional method, in which the driving data are read by the anti-lock braking system (ABS) speed signal, the GPS receiver can not only measure bus velocity with higher precision but also acquire the bus altitude and location for road slope estimation. The bus acceleration information is measured by the inertial navigation system, and the measurement error will not increase with the enhancement of sampling frequency. In order to retain more details of velocity fragments, the test sampling frequency is set to 5 Hz.

Road Slope Estimation
As indicated, other than velocity, road slope is also a particularly important characteristic for HEBs in control and design. This paper estimates road slope by calculating the ratio for the altitude difference between two sampling times to the horizontal distance traveled by the HEB during the sampling period and then taking the arctangent, which follows the equation where ∆ is the altitude change between the previous and current sampling time, is the bus horizontal velocity, is the sampling period. The road slope can be calculated with Equation (1); however, the computed preliminary results have high magnitude and some of the points obviously violate the longitude slope design requirements of the China code for design of urban road engineering [24]. Hence, the road slope preliminary values need to be filtered, the steps of which are as follows: 1. Primarily, eliminate slope calculated values with low velocity to avoid the denominator of Equation (1) being close to 0. 2. Carry out filtering. If the calculated slope values are higher than 9% [25], then replace them with those of the previous sampling time. 3. Perform average filtering against the processed slope values to smooth out the sawtooth signals.
The road slope after filtering is partially shown in Figure 3.

Original Data Preprocessing
The original data preprocessing consists of four steps: velocity fragment division, fragment state cluster distribution analysis, the Markov transfer matrix calculation, and characteristic parameter analysis. Among these four steps, the transfer matrix calculation is the most essential step of data preprocessing and driving cycle construction. The velocity fragment division and its state cluster distribution analysis have a direct influence on the accuracy of the Markov transfer matrix, and analysis of characteristic parameters is a criterion used to judge whether the characteristics of the constructed driving cycle agree with that of the original driving data.

Original Data Preprocessing
The original data preprocessing consists of four steps: velocity fragment division, fragment state cluster distribution analysis, the Markov transfer matrix calculation, and characteristic parameter analysis. Among these four steps, the transfer matrix calculation is the most essential step of data preprocessing and driving cycle construction. The velocity fragment division and its state cluster distribution analysis have a direct influence on the accuracy of the Markov transfer matrix, and analysis of characteristic parameters is a criterion used to judge whether the characteristics of the constructed driving cycle agree with that of the original driving data.

Velocity Fragment Division
A velocity fragment is the minimum unit of the driving cycle construction, which is abstracted from the original data according to specific rules.

Velocity Threshold Determination
The HEB generally drives in pure electric mode at a low velocity without an idling state. Therefore, it is necessary to consider the distinction between the shutdown state and driving state at a low velocity. The HEB operates in a shutdown state when the current velocity is smaller than a threshold value. However, when the velocity threshold is set to be extremely small, as shown in Figure 4, the fluctuant velocity in a small range will be cut into a great many shutdown fragments, which does not match the actual situation and will affect the analysis of characteristic parameters.

Velocity Fragment Division
A velocity fragment is the minimum unit of the driving cycle construction, which is abstracted from the original data according to specific rules.

Velocity Threshold Determination
The HEB generally drives in pure electric mode at a low velocity without an idling state. Therefore, it is necessary to consider the distinction between the shutdown state and driving state at a low velocity. The HEB operates in a shutdown state when the current velocity is smaller than a threshold value. However, when the velocity threshold is set to be extremely small, as shown in Figure 4, the fluctuant velocity in a small range will be cut into a great many shutdown fragments, which does not match the actual situation and will affect the analysis of characteristic parameters.
Based on the above consideration, the velocity sampling points below 4 km/h are separated from the original data and the distribution of different velocity intervals is analyzed. The corresponding statistical data are shown in Table 3.
It can be seen that the sampling points below 0.7 km/h can represent 67.33% of the population. Meanwhile, the velocity threshold value should be constricted as an appropriately small value. Thus, 0.7 km/h is selected as the velocity threshold of stop fragments after comprehensive consideration.

Acceleration Threshold Determination
The HEB is considered to travel at a uniform velocity when the bus's acceleration fluctuates within the threshold range. As in the previous discussion, when the velocity is below 0.7 km/h, the sampling points will be divided into shutdown fragments, and the acceleration of these sampling points should also theoretically fluctuate within a small range. The absolute values of acceleration distribution of shutdown fragments are analyzed, and the statistical data are shown in Table 4.
As shown in Table 4, around 92.01% of acceleration sampling points are in the interval of (0, 0.1 m/s 2 ), which indicates that the HEB drives at a uniform velocity when the acceleration is in the interval of (−0.1, 0.1). Considering that the road bump may magnify the measurement error of acceleration, the selection value of final acceleration threshold can approximately rise to 0.15 m/s 2 .  Based on the above consideration, the velocity sampling points below 4 km/h are separated from the original data and the distribution of different velocity intervals is analyzed. The corresponding statistical data are shown in Table 3. It can be seen that the sampling points below 0.7 km/h can represent 67.33% of the population. Meanwhile, the velocity threshold value should be constricted as an appropriately small value. Thus, 0.7 km/h is selected as the velocity threshold of stop fragments after comprehensive consideration.

Acceleration Threshold Determination
The HEB is considered to travel at a uniform velocity when the bus's acceleration fluctuates within the threshold range. As in the previous discussion, when the velocity is below 0.7 km/h, the sampling points will be divided into shutdown fragments, and the acceleration of these sampling points should also theoretically fluctuate within a small range. The absolute values of acceleration distribution of shutdown fragments are analyzed, and the statistical data are shown in Table 4. As shown in Table 4, around 92.01% of acceleration sampling points are in the interval of (0, 0.1 m/s 2 ), which indicates that the HEB drives at a uniform velocity when the acceleration is in the interval of (−0.1, 0.1). Considering that the road bump may magnify the measurement error of acceleration, the selection value of final acceleration threshold can approximately rise to 0.15 m/s 2 .

Basis of Velocity Fragment Division
In general, the values of both sampling velocity and acceleration should account for fragment division. The original data are divided into four types of velocity fragment: shutdown (stop) fragment, acceleration fragment, deceleration fragment, and uniform fragment. The basis of fragment division is shown in Figure 5.

Calculation of Required Power and Nominal Velocity
Every sampling time in the original driving data corresponds to the information of the velocity, acceleration, road slope, etc. To comprehensively depict the driving feature with fragment cluster distribution, the required power is also an important factor, which can be computed with Equation where P0(a > 0.15) denotes the required power of acceleration fragment when bus acceleration a > 0.15 m/s 2 , P0(a ≤ 0.15) represents the required power of non-acceleration fragment when a ≤ 0.15 m/s 2 . ua is the bus velocity, i is the road slope, G is the bus gravity, f is the rolling resistance coefficient, CD is the drag coefficient, A is the bus windward area, and ηT means the bus powertrain mechanical efficiency. Required power is an unintuitive parameter of fragment division and state cluster distribution; meanwhile, as position information, the road slope cannot be directly added into the velocity-time trace. Therefore, the nominal velocity with the same required power and acceleration is a suitable

Calculation of Required Power and Nominal Velocity
Every sampling time in the original driving data corresponds to the information of the velocity, acceleration, road slope, etc. To comprehensively depict the driving feature with fragment cluster distribution, the required power is also an important factor, which can be computed with Equations (2) and (3).
where P 0 (a > 0.15) denotes the required power of acceleration fragment when bus acceleration a > 0.15 m/s 2 , P 0 (a ≤ 0.15) represents the required power of non-acceleration fragment when a ≤ 0.15 m/s 2 . u a is the bus velocity, i is the road slope, G is the bus gravity, f is the rolling resistance coefficient, C D is the drag coefficient, A is the bus windward area, and η T means the bus powertrain mechanical efficiency. Required power is an unintuitive parameter of fragment division and state cluster distribution; meanwhile, as position information, the road slope cannot be directly added into the velocity-time trace. Therefore, the nominal velocity with the same required power and acceleration is a suitable method to involve the road slope information.
The specific operations are shown in Equations (4) and (5). When a > 0.15 m/s 2 , When a ≤ 0.15 m/s 2 , where u denotes the nominal velocity.

State Cluster Distribution
Calculate the average nominal velocity of each fragment, then distribute all velocity fragments into 7 state clusters according to the calculated average values, which are (−∞, 5], (5,15], (15,25), (25,35), (45, 55), and (55, +∞). Note that the nominal velocity cannot be directly applied to the driving cycle construction, and it only works in the fragment state cluster distribution where the road slope information is added into the traditional velocity-time trace and helps to make the state cluster distribution more reasonable.

Calculation of the Markov Transfer Matrix
The bus driving process can be regarded as a discrete Markov chain. The alteration of bus velocity can be regarded as a transfer between the current and next state clusters of velocity fragment, and the transfer process satisfies the probability distribution. Define X(n); n ≥ 0 as the Markov chain and S as state space, namely P X(n + 1) = j X(n) = i Sustainability 2020, 12, x FOR PEER REVIEW 8 of 18 where u denotes the nominal velocity.

State Cluster Distribution
Calculate the average nominal velocity of each fragment, then distribute all velocity fragments into 7 state clusters according to the calculated average values, which are (−∞, 5], (5,15], (15,25), (25,35), (45, 55), and (55, +∞). Note that the nominal velocity cannot be directly applied to the driving cycle construction, and it only works in the fragment state cluster distribution where the road slope information is added into the traditional velocity-time trace and helps to make the state cluster distribution more reasonable.

Calculation of the Markov Transfer Matrix
The bus driving process can be regarded as a discrete Markov chain. The alteration of bus velocity can be regarded as a transfer between the current and next state clusters of velocity fragment, and the transfer process satisfies the probability distribution. Define { ( ); 0} X n n ≥ as the Markov chain and S as state space, namely where pij(n) denotes the one step transfer probability at moment n.
The meaning of Equation (6) is the probability of model event to transfer from time n with state i to time n + 1 with state j, and the transfer matrix which represents the fragment state migration can be constructed via the transfer probabilities with different fragment states.
When state clusters are distributed, every velocity fragment has a state number. A vehicle's driving condition transfers from one fragment to another together with the corresponding state cluster migration. The transfer matrix is calculated by including the migration of state clusters in the original data. The elements in the transfer matrix, namely the transfer probability, can be calculated as follows: where Nij is the time of velocity fragment transferring from state i at time t − 1 to state j at time t. Take the first row of the transfer matrix, for example. As shown in Table 5, it indicates the moments of velocity fragment transfer from state 1 to other states, and the corresponding transfer probabilities have also been counted but are omitted due to limited space. The transfer probability of other states can be calculated in the same way, and the 7 × 7 state transfer matrix T could be constituted as follows: where p ij (n) denotes the one step transfer probability at moment n. The meaning of Equation (6) is the probability of model event to transfer from time n with state i to time n + 1 with state j, and the transfer matrix which represents the fragment state migration can be constructed via the transfer probabilities with different fragment states.
When state clusters are distributed, every velocity fragment has a state number. A vehicle's driving condition transfers from one fragment to another together with the corresponding state cluster migration. The transfer matrix is calculated by including the migration of state clusters in the original data. The elements in the transfer matrix, namely the transfer probability, can be calculated as follows: where N ij is the time of velocity fragment transferring from state i at time t − 1 to state j at time t.
Take the first row of the transfer matrix, for example. As shown in Table 5, it indicates the moments of velocity fragment transfer from state 1 to other states, and the corresponding transfer probabilities have also been counted but are omitted due to limited space. The transfer probability of other states can be calculated in the same way, and the 7 × 7 state transfer matrix T could be constituted as follows:

Characteristic Parameter Matrix of Original Data
A 1 × 16 original characteristic matrix M 0 is then constructed to represent the original data. The PM (performance measure) value is the differences in characteristic parameters between the driving cycle and original data, which is normalized and used to represent the differentiating degree.

Calculation of PM Value
The calculation of PM value consists of four steps:

1.
Count the characteristic parameters of driving cycle candidates n and construct the n × 16 matrix M.

2.
Obtain the absolute differentiating value matrix pm through each row of matrix M minus the original characteristic M 0 : 3. Normalize pm. Normalization is necessary for totally different parameter ranges, which can be deduced as where pm imax is the maximum value of the matrix pm's ith column, while pm imin is the minimum one.

4.
Sum the elements in the kth row in the matrix PM (k,i) to obtain the PM value of the driving cycle.

Sufficient Proof of the Original Data
The original data are the basis of driving cycle construction and the constructed result will be more credible, with more collected original data. The volume of original data can have a strong influence on the driving cycle's representation. In this paper, the data acquisition experiment lasted for two weeks, and more than 600,000 sampling points were acquired, which provided sufficient data for subsequent research.
As we have constructed the driving cycle based on Markov chain theory, the credibility of the Markov transition matrix can be enhanced with a larger original database. Therefore, the original data are firstly fused day by day; then, the F norm of the differentiating matrix, which can characterize the variation in Markov transition matrixes with the increase in data amount, is calculated by evaluating the difference in the two adjacent matrixes. The computing results are shown in Figure 6.

Driving Cycle Construction
As indicated, due to the Markov chain, the velocity of the sampling point at the next moment is completely dependent on the velocity of the sampling point at the current moment and the transition probability corresponding to this state. However, at the beginning and end of the driving cycle, HEB's driving condition has been determined since HEB must gradually accelerate or decelerate towards a standstill state from other driving conditions. Adopting the Markov transfer matrix to construct these two parts is not advisable. Therefore, the start, middle, and end parts of the driving cycle are constructed, respectively, and then they are spliced according to the time sequence to accomplish the final driving cycle construction.

Construction of Start Part
The construction of the start part contains 2 steps. First, select all start part candidates from the original data with proposed rules. Then, the PM value of every part is calculated. To this end, we select 16 parameters to characterize the original data, as shown in Table 6.
The time length of the start part selected in this paper is 75 s. The stop duration before the bus standing start is mainly determined by the traffic condition. In Table 6, the proportion of stop time accounts for 24.09% of the original data; accordingly, the stop time in the 75 s duration of start part is no more than 18 s. Moreover, durations of stop fragments of more than 5 s make up 54.21% of all stop fragment durations, while durations longer than 10 s make up 32.85%. The stop fragments that can be used for comparison decrease with the increase in duration, which develops in an opposite trend to that for the start part selection. Therefore, the stop duration of the stop part candidate is set to 5 s.
Generally, a long-time parking situation does not appear in the start part, so the time proportion of 4 kinds of velocity fragments, the maximum, average, and minimum velocity, will not be considered during the comparison of characteristic parameters in the PM value calculation process. Only the difference between the candidate start part and the original data, which represents the acceleration, road slope, and road power, should be considered.
Then, all the start part candidates are sorted according to the PM value; the one which has the lowest PM value will be the optimal selection.  The F norm decreases gradually, and it fluctuates around a stable value after the 4th day. This phenomenon demonstrates that the Markov transfer matrix is stable and the original driving data are sufficient.

Driving Cycle Construction
As indicated, due to the Markov chain, the velocity of the sampling point at the next moment is completely dependent on the velocity of the sampling point at the current moment and the transition probability corresponding to this state. However, at the beginning and end of the driving cycle, HEB's driving condition has been determined since HEB must gradually accelerate or decelerate towards a standstill state from other driving conditions. Adopting the Markov transfer matrix to construct these two parts is not advisable. Therefore, the start, middle, and end parts of the driving cycle are constructed, respectively, and then they are spliced according to the time sequence to accomplish the final driving cycle construction.

Construction of Start Part
The construction of the start part contains 2 steps. First, select all start part candidates from the original data with proposed rules. Then, the PM value of every part is calculated. To this end, we select 16 parameters to characterize the original data, as shown in Table 6.
The time length of the start part selected in this paper is 75 s. The stop duration before the bus standing start is mainly determined by the traffic condition. In Table 6, the proportion of stop time accounts for 24.09% of the original data; accordingly, the stop time in the 75 s duration of start part is no more than 18 s. Moreover, durations of stop fragments of more than 5 s make up 54.21% of all stop fragment durations, while durations longer than 10 s make up 32.85%. The stop fragments that can be used for comparison decrease with the increase in duration, which develops in an opposite trend to that for the start part selection. Therefore, the stop duration of the stop part candidate is set to 5 s.
Generally, a long-time parking situation does not appear in the start part, so the time proportion of 4 kinds of velocity fragments, the maximum, average, and minimum velocity, will not be considered during the comparison of characteristic parameters in the PM value calculation process. Only the difference between the candidate start part and the original data, which represents the acceleration, road slope, and road power, should be considered. Then, all the start part candidates are sorted according to the PM value; the one which has the lowest PM value will be the optimal selection.

Construction of Middle Part
The construction of the middle part is solved mainly with the Monte Carlo method. A random number is generated to combine the probability distribution of the Markov transfer matrix, and then the next velocity fragment's state cluster can be determined.
The time length of the middle part is designed to be more than 1000 s. If the start part has been constructed, the velocity and state cluster of the end point, which is also the initial point of the middle part, can be acquired.
The specific selection rules for the middle part are shown below: 1.
Confirm the current state cluster and the transfer probability, which are determined according to the Markov transfer matrix T. P X(n) = j X(n − 1) = i = P ij (11)

2.
Generate a random number r in the interval of (0, 1), if r meets the in-equation as follows: Then, the next state cluster of velocity fragment is q.

3.
Determine the velocity fragments whose velocity differences in current and next state clusters are below 0.5 km/h. The one which has the minimum velocity difference will be the next determined velocity fragment. 4.
The sampling points which are used in the start part construction will not be used in the middle part construction, and 10 middle part candidates are constructed for optimal selection.
Combined with the constructed start part for PM value calculation, in this section, all 16 characteristic parameters are used for calculation. We chose the middle part candidate with the lowest PM value to construct the middle part of the driving cycle.

Construction of End Part
The method of end part construction is similar to that of start part; it is noteworthy that the end part candidates should be combined with the constructed start and middle parts for the PM value calculation, and similarly, the one which has the lowest PM value will be chosen as the final end part of the driving cycle considering all 16 characteristic parameters.
The selection method of the end part is the same as that of the start part. The specific selection rules for the end part are shown below: 1.
The velocity of the sampling points that are 10 s prior to the checkpoint is above velocity threshold 0.7 km/h.

2.
The velocity sampling points which are 5 s later than the checkpoint are below the velocity threshold 0.7 km/h. 3.
If both the above conditions are met, set the checkpoint at the time t, then the sampling points in the time interval of (t − 114, t + 5) will be chosen as the end part candidates.

4.
Check all sampling points in the first 60 s of every end part candidate, determine 10 whose velocities in sampling points are closest to the terminal velocity of the constructed middle part, and then set these points as the start times of the 10 end part candidates.

Construction of Driving Cycle
The PM values of the 3 parts of the driving cycle are sorted in size, as shown in Table 7. The selection details of the 3 final parts of the driving cycle can be apparently represented with the PM values. As shown in the first column of the table, the value of candidate start parts is 913, and those of the other two are both 10. The first row of the table is the lowest PM value for every part's optimal selection. According to the time sequence, the start, middle, and end parts of the driving cycle, which have the lowest PM values in each group, are spliced to form a complete time-velocity curve, as shown in Figure 7, which indicates that the driving cycle of the hybrid electric bus based on Markov chain has been constructed. According to the time sequence, the start, middle, and end parts of the driving cycle, which have the lowest PM values in each group, are spliced to form a complete time-velocity curve, as shown in Figure 7, which indicates that the driving cycle of the hybrid electric bus based on Markov chain has been constructed.  The mileage of the constructed driving cycle can be acquired through the integration of the velocity in time domain. Through mapping the road slope values of corresponding sampling points, the distance-slope curve shown in Figure 8 is acquired. The figure indicates that the driving distance of one driving cycle is 5 km, and the road slope is less than 4 degrees, which meets the requirements for the designing of urban road engineering.
Moreover, the velocity-acceleration probability distribution (VAPD) is a decisive criterion for the representative evaluation of the constructed driving cycle [26]. The VAPDs of the original data and the constructed driving cycle are shown in Figure 9. It can be seen that the VAPD of the constructed driving cycle is similar to that of the original data, and both of them have a long stop The mileage of the constructed driving cycle can be acquired through the integration of the velocity in time domain. Through mapping the road slope values of corresponding sampling points, the distance-slope curve shown in Figure 8 is acquired. The figure indicates that the driving distance of one driving cycle is 5 km, and the road slope is less than 4 degrees, which meets the requirements for the designing of urban road engineering. velocity in time domain. Through mapping the road slope values of corresponding sampling points, the distance-slope curve shown in Figure 8 is acquired. The figure indicates that the driving distance of one driving cycle is 5 km, and the road slope is less than 4 degrees, which meets the requirements for the designing of urban road engineering.
Moreover, the velocity-acceleration probability distribution (VAPD) is a decisive criterion for the representative evaluation of the constructed driving cycle [26]. The VAPDs of the original data and the constructed driving cycle are shown in Figure 9. It can be seen that the VAPD of the constructed driving cycle is similar to that of the original data, and both of them have a long stop period with relatively gentle acceleration and deceleration.  Moreover, the velocity-acceleration probability distribution (VAPD) is a decisive criterion for the representative evaluation of the constructed driving cycle [26]. The VAPDs of the original data and the constructed driving cycle are shown in Figure 9. It can be seen that the VAPD of the constructed driving cycle is similar to that of the original data, and both of them have a long stop period with relatively gentle acceleration and deceleration.

Comparison between the Driving Cycles Based on Markov and Traditional Micro-Trip
The micro-trip method is a traditional way to construct a driving cycle [27]. This paper adopts the same original date to construct a 1200 s driving cycle based on the micro-trip method. It is noteworthy that the velocity and acceleration thresholds for micro-trip separation are the same as that of the Markov chain-based method. The results of comparing the driving cycles constructed with the two methods are shown in Figure 10 and Table 8. In Table 8, most of the characteristic parameters of the driving cycle constructed with the Markov method are closer to those of the original data than the micro-trip method.
Furthermore, the similarity degree of the constructed driving cycle and the original data is represented by the correlation coefficient of characteristic parameters in this paper. The calculating equation is expressed as where X represents the characteristic parameters of the original data, Y represents the characteristic parameters of the constructed driving cycle, Cov (X, Y) is the covariance of X and Y, and Var(X), Var(Y) are the variance of X and Y, respectively. According to Table 8 and Equation (13), the correlation coefficient of original data and driving cycle based on the Markov chain is 0.9972, which is higher than that of the traditional micro-trip, with

Comparison between the Driving Cycles Based on Markov and Traditional Micro-Trip
The micro-trip method is a traditional way to construct a driving cycle [27]. This paper adopts the same original date to construct a 1200 s driving cycle based on the micro-trip method. It is noteworthy that the velocity and acceleration thresholds for micro-trip separation are the same as that of the Markov chain-based method. The results of comparing the driving cycles constructed with the two methods are shown in Figure 10 and Table 8. In Table 8, most of the characteristic parameters of the driving cycle constructed with the Markov method are closer to those of the original data than the micro-trip method.

Dynamic Programming
The experiment vehicle is a plug-in hybrid bus [28], which weight 1,3100 kg and has driven on the constructed driving cycle ZZUDC. As for the parameters of every component-for example, engine, battery, motor, etc.-these are omitted, but a related paper [28] is a decent reference.
The numbers of every driving cycle are evaluated with the same constraints. The target function of dynamic programming (DP) is St. The constraints of this problem come from both state variables and control singles. Except for the restrictions in Equation (15), the maximum charging/discharging rate can also affect the slope of the alteration of SOC. Specific constraints of the control singles derive from the related components' limiting speed and torque [29].

Comparison between the Statistical Characters of Various Driving Cycles
In order to assess the performance of ZZUDC more comprehensively, we selected seven standard international driving cycles as comparisons, which are CTUB (China typical urban bus driving cycle), LA92, IM240, JN1015 (Japan 1015), UDDS (urban dynamometer driving schedule), MANHATTAN, and WVUCITY. As can be seen in Table 9, most of the selected international driving cycles are urban route types. The total cycle durations significantly vary between cycles, ranging from  Furthermore, the similarity degree of the constructed driving cycle and the original data is represented by the correlation coefficient of characteristic parameters in this paper. The calculating equation is expressed as where X represents the characteristic parameters of the original data, Y represents the characteristic parameters of the constructed driving cycle, Cov(X,Y) is the covariance of X and Y, and Var(X), Var(Y) are the variance of X and Y, respectively. According to Table 8 and Equation (13), the correlation coefficient of original data and driving cycle based on the Markov chain is 0.9972, which is higher than that of the traditional micro-trip, with 0.9746. The results indicate that the Markov chain-based driving cycle is more objective and effective.

Dynamic Programming
The experiment vehicle is a plug-in hybrid bus [28], which weight 13,100 kg and has driven on the constructed driving cycle ZZUDC. As for the parameters of every component-for example, engine, battery, motor, etc.-these are omitted, but a related paper [28] is a decent reference.
The numbers of every driving cycle are evaluated with the same constraints. The target function of dynamic programming (DP) is St.
where SOC (state of charge) is a single state variable, function f illustrates the time-variant model of the vehicle, SOC 0 and SOC N are the initial and final figures, respectively, and SOC k is the state value at kth discretized point. The control single is u, including five elements which are u{1}, u{2}, u{3}, u{4}, and u{5}. The meanings of control singles are torque of engine, torque of ISG, speed of engine, state of clutch, and torque of mechanical braking, respectively. All control singles are continuous types, except for the state of clutch (u{4}), which is a 0-1 discrete type. The constraints of this problem come from both state variables and control singles. Except for the restrictions in Equation (15), the maximum charging/discharging rate can also affect the slope of the alteration of SOC. Specific constraints of the control singles derive from the related components' limiting speed and torque [29].

Comparison between the Statistical Characters of Various Driving Cycles
In order to assess the performance of ZZUDC more comprehensively, we selected seven standard international driving cycles as comparisons, which are CTUB (China typical urban bus driving cycle), LA92, IM240, JN1015 (Japan 1015), UDDS (urban dynamometer driving schedule), MANHATTAN, and WVUCITY. As can be seen in Table 9, most of the selected international driving cycles are urban route types. The total cycle durations significantly vary between cycles, ranging from 241 s (IM240) to 1408 s (WVUCITY), and mainly distribute between 1000 s and 1500 s. The duration of ZZUDC is 1184 s, which is an average level. The major difference is the highest velocity, which ranges from around 40 km/h to above 108.15 km/h. This is because different cities can have different traffic conditions, and the driving cycles should reflect the distinctions. In addition, some of the driving cycles, such as LA92, are composite driving cycles and contain both urban route types and suburban route types. The suburban part would certainly amplify the parameter of highest velocity.
However, values of some assessment parameters of the MANHATTAN and WVUCITY are close to the constructed driving cycle ZZUDC, including duration, mileage, average velocity, maximum velocity, average acceleration, and deceleration. The traveling distances of WVUCITY and ZZUDC are both around 5 km, and the gaps of average acceleration and deceleration between these two cycles are within 10%. The maximum acceleration, maximum deceleration, and analogous root mean square of acceleration (RMS of acc) of ZZUDC are similar to those of the MANHATTAN driving cycle. Moreover, the energy consumption of these two cycles is at the same level, with around 21.3 L/100 km for fuel consumption and approximately 13.3 kwh/100 km for electricity consumption for MANHATTAN and around 20.5 L/100 km fuel consumption and approximately 12.8 kwh/100 km electricity consumption for ZZUDC. Furthermore, in order to visually depict the relative distributions of all the cycles in every statistical item, we can normalize the figures with every line according to (16) where x i is the specific figure in a vector. The results of all the driving cycles are displayed in Figure 11. The normalization range is between 0 and 1. If most of the figures converge in one certain pole (0 or 1) in any statistical indicator, this indicates that there is a stand-out. However, in Figure 11, the plots of eight cycles, ZZUDC included, are substantially decentralized in every item. Thus, ZZUDC is mediocre, rather than extreme, which verifies the rationality of the cycle. Besides this, as we mentioned before, the characteristics of ZZUDC are similar to those of MANHATTAN and WUVCITY. The trends of their broken lines prove this conclusion. As for the distinction in fuel consumption of ZZUDC and WUVCITY, the dominant reason is that the required power of the bus when it drives on WUVCITY is at a lower level than on ZZUDC, which can be reflected, to some extent, by the summit power. The normalization range is between 0 and 1. If most of the figures converge in one certain pole (0 or 1) in any statistical indicator, this indicates that there is a stand-out. However, in Figure 11, the plots of eight cycles, ZZUDC included, are substantially decentralized in every item. Thus, ZZUDC is mediocre, rather than extreme, which verifies the rationality of the cycle. Besides this, as we mentioned before, the characteristics of ZZUDC are similar to those of MANHATTAN and WUVCITY. The trends of their broken lines prove this conclusion. As for the distinction in fuel consumption of ZZUDC and WUVCITY, the dominant reason is that the required power of the bus when it drives on WUVCITY is at a lower level than on ZZUDC, which can be reflected, to some extent, by the summit power.

Discussion
The validation process of ZZUDC mainly conducted from the perspective of the effectiveness of the construction method and the driving cycle's statistical items. Firstly, the correlation coefficients of characteristic parameters were computed and VPAD figures were plotted, which indicated that the constructed driving cycle (ZZUDC) largely maintains the original data. Moreover, it is evident that, compared with the traditional micro-trip construction method, ZZUDC reflects most of the original data's characteristics better. Besides this, among ZZUDC and the other seven international driving cycles, ZZUDC is distributed in the middle of them and similar to MANHATTAN and WVUCITY with respect to the majority of normalized statistical items. Therefore, all results show that the constructed ZZUDC is reasonable and credible.

Conclusions
In this study, the driving cycle data of two HEB routes in Zhengzhou city are collected with a measurement system integrating global positioning and inertial navigation function. New methods for velocity fragment division and state cluster distribution, improved PM calculation method, and dynamic programming are applied to develop the Zhengzhou urban driving cycle for HEB. We construct the driving cycle parts separately, namely start part, middle part, and end part. During the construction process, characteristic matrices for each part of the driving cycle are established, using 16 characteristic parameters, including road slope, in which mathematical statistic theories related to Markov chain are used. Then, a special method of PM calculation, including the normalization process, is presented to determine the final driving cycle, which can eliminate the difference that is caused by the non-uniform unit of different characteristic parameters and guarantee the reasonability of velocity fragment selection.
By comparing the correlation coefficients of the driving cycles that are constructed based on Markov and traditional micro-trip, the driving cycle constructed in this paper is confirmed to be more objective and effective, with a higher coefficient with the original data. Furthermore, through analyzing the statistical characteristics of ZZUDC and seven international cycles, the target cycle is proven to be reasonable. As for the energy consumption, the values of ZZUDC and MANHATTAN are similar, combining dynamic programming. To be specific, fuel consumption is around 20.5 L/100 km and electricity consumption is approximately 12.8 kwh/100 km.
The established HEB driving cycle construction approach can provide a reference for developing a driving cycle for HEV that drives in pure electric mode at a low velocity without an idling state. Furthermore, this paper lays the foundation for using dynamic programming to compute fuel consumption and electricity consumption for HEBs. In future research, we will investigate the possibility of developing an energy consumption management strategy to improve HEB emission.