Research on Cooperative Braking Control Algorithm Based on Nonlinear Model Prediction

: To address the coordinated distribution of motor braking and friction braking for the regenerative braking system, a cooperative braking algorithm based on nonlinear model predictive control (NMPC) is proposed, with braking energy recovery power, tire slip rate, and motor torque variation as the optimization objectives, and online optimization of the coordinated distribution of motor braking and friction braking. Using the ofﬂine model built in Matlab/Simulink, the cooperative braking algorithm is tested for energy efﬁciency and braking safety. The results show that when based on World Light Vehicle Test Cycle (WLTC), the energy recovery rate can reach 30.4%, and with a single high braking intensity, the braking safety can still be ensured.


Introduction
As the automotive industry gradually enters the era of electrification and intelligentization, the fully decoupled braking system has become a key research object for major OEMs and parts manufacturers [1]. The increase in the degree of decoupling of the braking system has laid the hardware foundation for the implementation of advanced braking control algorithms [2,3]. Among them, the cooperative braking control algorithm has been a hot topic of research in the industry, which is mainly used to solve the problem of coordinated distribution of motor braking and friction braking for the regenerative braking system and has an important impact on braking safety and braking energy recovery effect [4][5][6].
The cooperative braking control strategy mainly includes two types of parallel and series, and the schematic diagram is shown in Figure 1. For the parallel control strategy, motor braking and friction braking are independent of each other [7][8][9]. Motor braking is additionally superimposed on the drive shaft while keeping the conventional front and rear axle friction braking ratio unchanged. It can be easily achieved in the conventional undecoupled hardware system [10][11][12]. The additional motor braking power should not be too large; otherwise, it will affect the braking feeling and safety, so the regenerative braking system with a parallel control strategy will only recover less braking energy. For the series control strategy, friction braking and motor braking can be coordinatively controlled so that the motor braking capacity can be fully utilized while ensuring braking feeling and safety [13,14]. The regenerative braking system with a series control strategy can recover more braking energy, but its algorithm is more complex and must be implemented in the decoupled hardware system [15].
In recent years, many advanced cooperative braking algorithms have been derived based on these two basic types. The paper [16] takes the Economic Commission of Europe (ECE) braking safety regulation as the constraint condition and develops a new regenerative braking torque distribution strategy. To improve energy recovery efficiency, an optimization algorithm of energy recovery efficiency is proposed for parallel hydraulic hybrid systems (PHHS) using dynamic programming (DP) [17]. The paper [18] designs an energy-recovery control strategy based on the Kalman Filter, which can optimize the model accuracy and improve the robustness of the system. Based on logic threshold control algorithm with the input of z and I, a distribution relationship between electric and mechanical braking is determined, which is proved to be effective on road conditions with different slopes [19]. Based on the ECE regulation curve and ideal braking force distribution (I curve), a braking force distribution strategy of the front and rear axles is designed to improve the recovered braking energy [20]. In recent years, many advanced cooperative braking algorithms have been derived based on these two basic types. The paper [16] takes the Economic Commission of Europe (ECE) braking safety regulation as the constraint condition and develops a new regenerative braking torque distribution strategy. To improve energy recovery efficiency, an optimization algorithm of energy recovery efficiency is proposed for parallel hydraulic hybrid systems (PHHS) using dynamic programming (DP) [17]. The paper [18] designs an energy-recovery control strategy based on the Kalman Filter, which can optimize the model accuracy and improve the robustness of the system. Based on logic threshold control algorithm with the input of z and I, a distribution relationship between electric and mechanical braking is determined, which is proved to be effective on road conditions with different slopes [19]. Based on the ECE regulation curve and ideal braking force distribution (I curve), a braking force distribution strategy of the front and rear axles is designed to improve the recovered braking energy [20].
In summary, despite the fact that various advanced braking control algorithms have been proposed, it is mostly true that they only optimize a single control objective, while control of other objectives is rather vague and thus cannot systematically address the problem of multiple objectives and constraints in regenerative braking systems. In addition, the control objectives cannot be optimized in real time according to the change of vehicle state, which means the regenerative braking system may not be able to be fully exploited. Therefore, based on the above analysis, a cooperative braking algorithm based on NMPC is proposed for the fully decoupled braking system of front-drive vehicles in this paper, with braking energy recovery power, tire slip rate, and motor torque variation as the optimization objectives, and online optimization of the coordinated distribution of motor braking and friction braking to achieve good energy efficiency while ensuring braking safety [19,21,22].
The remainder of this paper is organized as follows. Section 2 describes the cooperative braking control algorithm based on NMPC. Simulation verification and results analysis are discussed in Section 3. Lastly, conclusions are drawn in Section 4.

Overview of NMPC
The control process that is based on the system dynamic model to find the optimal solution of the optimization problem in a finite prediction time domain and utilizes the rolling optimization to complete the optimization of the whole control process, which is defined as model predictive control (MPC) [23,24].
Classical MPC usually uses linear models for prediction, which works better for linear systems. But for nonlinear systems, linear models cannot accurately describe their dy-  In summary, despite the fact that various advanced braking control algorithms have been proposed, it is mostly true that they only optimize a single control objective, while control of other objectives is rather vague and thus cannot systematically address the problem of multiple objectives and constraints in regenerative braking systems. In addition, the control objectives cannot be optimized in real time according to the change of vehicle state, which means the regenerative braking system may not be able to be fully exploited. Therefore, based on the above analysis, a cooperative braking algorithm based on NMPC is proposed for the fully decoupled braking system of front-drive vehicles in this paper, with braking energy recovery power, tire slip rate, and motor torque variation as the optimization objectives, and online optimization of the coordinated distribution of motor braking and friction braking to achieve good energy efficiency while ensuring braking safety [19,21,22].
The remainder of this paper is organized as follows. Section 2 describes the cooperative braking control algorithm based on NMPC. Simulation verification and results analysis are discussed in Section 3. Lastly, conclusions are drawn in Section 4.

Overview of NMPC
The control process that is based on the system dynamic model to find the optimal solution of the optimization problem in a finite prediction time domain and utilizes the rolling optimization to complete the optimization of the whole control process, which is defined as model predictive control (MPC) [23,24].
Classical MPC usually uses linear models for prediction, which works better for linear systems. But for nonlinear systems, linear models cannot accurately describe their dynamic characteristics, making it difficult to apply them directly to nonlinear systems.
MPC is also a closed-loop optimization method that can correct the error caused by model mismatch through the closed-loop process. In practice, MPC can be applied in nonlinear processes by linearizing the nonlinear model. This control process is known as NMPC [25,26].
As a branch of MPC, NMPC can achieve predictive control based on the nonlinear model and perform closed-loop optimization. The nonlinear model can provide a more accurate description of the nonlinear process, solving the problem of insufficient accuracy of the linear model [27,28]. NMPC can only find the optimal solution in each finite time domain but not the global optimal solution of the whole operation process. In engineering practice, the target state of the vehicle in the future cannot be predicted in advance, so the global optimum cannot be found. Therefore, solving the local optimal results can already satisfy the optimization requirements of the regenerative braking process [29].

The Implementation Method of NMPC
Different NMPC algorithms differ in terms of prediction models and control performance, but all have the same implementation. In each control cycle, with the current state as the starting point, the control algorithm can solve the optimization problem in the finite time domain by a prediction model to obtain the optimal control sequence, the first element of which is used as the control input of the current state. The implementation of the algorithm consists of three parts: solving the prediction model, solving the optimization problem in the finite time domain, and feedback correction.
The first step in implementing the NMPC algorithm is to build a discrete prediction model. The predictive model can respond to system dynamic characteristics and be used to predict system future states and outputs. By solving the prediction model, each set of system states and outputs in the prediction time domain can be obtained.
To solve the optimization problem in the finite time domain, it is first necessary to design the optimization objective function. The optimization objective function contains the outputs and state variables to be optimized and contains the weights for each of these optimization objectives. Solving an optimization problem in the finite time domain is actually to find the input matrix that minimizes the value of the optimization objective function in the prediction time domain. The time range of the input matrix is called the control time domain, and usually, the control time domain is less than or equal to the prediction time domain. To ensure that the states and inputs can be solved for each moment in the prediction time domain, the last set of inputs of the input matrix can be taken as the inputs outside the control time domain.
In the actual control process, there is a common problem of model mismatch. Using closed-loop feedback, the NMPC algorithm can correct control errors caused by the model mismatch. The NMPC algorithm can correct control errors caused by model mismatch through closed-loop feedback.
The feedback correction consists of two parts: first, when solving the optimization problem at each moment, the prediction model starts predicting based on the actual feedback of the system at the previous moment, which can reduce the global cumulative error caused by the accumulation of model mismatch over time; second, the control algorithm can compare the prediction data before the current moment with the actual measurement results to obtain the prediction error caused by model mismatch, and use the prediction error to correct the output of the future prediction model at each moment, which can further reduce the local prediction error caused by the model mismatch.
When in the regenerative braking state, the algorithm control flow chart is shown in Figure 2.

Build Predictive Model
The braking energy recovery process is influenced by the wheels, the hydraulic braking system, the motor, and the battery; to accurately calculate the motion of the wheels, the load transfer of the sprung mass also needs to be predicted. Therefore, the desired predictive model should include models for the sprung mass, tires, hydraulic braking system, motor, and battery. To ensure the timeliness of online rolling optimization, the predictive model should be simplified as much as possible.
The reference vehicle force analysis curve is shown in Figure 3, and the basic parameters of the vehicle are shown in Table 1.  Figure 2. The nonlinear model predictive control process.

Build Predictive Model
The braking energy recovery process is influenced by the wheels, the hydraulic braking system, the motor, and the battery; to accurately calculate the motion of the wheels, the load transfer of the sprung mass also needs to be predicted. Therefore, the desired predictive model should include models for the sprung mass, tires, hydraulic braking system, motor, and battery. To ensure the timeliness of online rolling optimization, the predictive model should be simplified as much as possible.
The reference vehicle force analysis curve is shown in Figure 3, and the basic parameters of the vehicle are shown in Table 1. Where is the sprung mass, is the front wheel ground braking force, is the rear wheel ground braking force, is the air resistance, is the wheelbase, is the horizontal distance from the center of the sprung mass to the front axle, is the horizontal distance from the center of the sprung mass to the rear axle, ℎ is the height of the centroid, is the front wheel rotational speed, is the front wheel rotational inertia,

Build Predictive Model
The braking energy recovery process is influenced by the wheels, the hydraulic braking system, the motor, and the battery; to accurately calculate the motion of the wheels, the load transfer of the sprung mass also needs to be predicted. Therefore, the desired predictive model should include models for the sprung mass, tires, hydraulic braking system, motor, and battery. To ensure the timeliness of online rolling optimization, the predictive model should be simplified as much as possible.
The reference vehicle force analysis curve is shown in Figure 3, and the basic parameters of the vehicle are shown in Table 1. Where is the sprung mass, is the front wheel ground braking force, is the rear wheel ground braking force, is the air resistance, is the wheelbase, is the horizontal distance from the center of the sprung mass to the front axle, is the horizontal distance from the center of the sprung mass to the rear axle, ℎ is the height of the centroid, is the front wheel rotational speed, is the front wheel rotational inertia,  Where m sm is the sprung mass, F b1 is the front wheel ground braking force, F b2 is the rear wheel ground braking force, F a is the air resistance, l is the wheelbase, a is the horizontal distance from the center of the sprung mass to the front axle, b is the horizontal distance from the center of the sprung mass to the rear axle, h g is the height of the centroid, ω 1 is the front wheel rotational speed, J w1 is the front wheel rotational inertia, J m is the motor rotational inertia, i t is the total rotation ratio, T m is the motor torque, R t is the tire rolling radius, T f 1 is the front wheel friction braking torque, ω 2 is the rear wheel rotational speed, J w2 is the rear wheel rotational inertia, and T f 2 is the rear wheel friction braking torque.
The sprung mass is modeled as a planar single degree of freedom model. The sprung mass moves along a straight line with velocity v. Its motion equation is: The equation of air resistance is: where C d is the air resistance coefficient, A F is the windward area, and ρ a is air density.
where g is gravity.
Changes in vertical load and slip rate affect the longitudinal force of the tire, and to more accurately reflect this one, the Pacejka model is used to model the tire. So, regarding the relationship between the tire longitudinal force and the slip rate, the vertical load can be expressed as: where K B , K C , K D , and K E , the coefficients of the Pacejka model, are constants, F zi is the vertical load of each wheel, and s i is the slip rate.
The slip rate is related to the rotational speed and longitudinal speed of the wheel. Assuming that the wheel is rigidly connected to the sprung mass, its longitudinal speed is the same as the sprung mass and the rotational speed can be calculated from the torque acting on the wheel and its rotational inertia. The front wheel is the driving wheel and is subjected to the torque formed by the motor torque, friction braking torque, and groundbraking force acting on the wheel, so the motion equation of the front wheel can be expressed as follows: The rear wheel is the driven wheel, so the applied torque does not include the motor torque. Assuming that the front and rear wheels have the same rolling radius, the motion equation of the rear wheel can be expressed as: Because the friction braking torque ratio of the front and rear wheel brakes is constant, the driven wheel motion equation can also be expressed as: The total braking demand torque T total is equal to the sum of the total friction braking torque and the motor braking torque; the formula is: The total braking demand torque is obtained through the braking intent recognition module, and it is a constant while solving for the optimization results. Therefore, the front wheel friction braking torque can also be expressed by T m ; the relationship equation is: By establishing the relationship between the friction braking torque and the motor braking torque, the input dimension of the model can be reduced and the optimization problem can be simplified.
The front wheel is rigidly connected to the drive motor through the transmission system, so the motor speed is known from the front wheel speed; the relationship equation is: When the speed and torque are given, the motor output current is: where Loss(T m , ω m ) is the power loss of the motor, which can be obtained by bench testing the motor at different speeds and torques; V acc is the voltage of the power battery, which is related to the battery current and state of charge (SOC), but for each control process, the battery voltage does not change much, so V acc is considered as a constant.
To apply the above model to the NMPC algorithm, it needs to be discretized. With the sampling period of ∆t, the discretized model is: The air resistance after discretization is: The front axle load and rear axle load after discretization are: World Electr. Veh. J. 2021, 12, 173 7 of 16 The tire longitudinal force after discretization is: The slip rate after discretization is: The motor output current after discretization is: The input in the above model is: At the moment k, the input will be optimized online by the NMPC algorithm and the first element of the optimal control sequence will be used as the input to the system at the current moment. The above model involves the state variable at moment k-1, but when k = 0, it does not exist. In practical engineering, the controller can always monitor the state variable online, so when the braking starts, the controller can input the state variable measured at the previous moment as the state variable at moment k-1 into the model, so that the model can work normally.

Design Optimization Objective Function
The optimization objective function reflects the objective variables and weights that need to be optimized for the optimization problem. The optimization objectives of the control strategy proposed in this paper are the braking energy recovery power, tire slip rate, and motor torque variation. To recover as much braking energy as possible, the braking energy recovery power should be increased as much as possible while ensuring that the tire slip rate is within the safe range. Throughout the braking process, the braking torque is provided by both the motor and the friction braking system. The motor braking torque can be changed quickly by adjusting the current, but the change rate of the friction braking torque is smaller. Therefore, the change of the motor torque should be as small as possible to ensure smooth deceleration during the braking process.
According to the above optimization objectives, the optimization objective function can be expressed as: where ∆u k is the motor torque variation at moment k, s 1 (k + i|k) is the slip rate at moment k + i that is predicted at moment k, and P m (k + i|k) is the braking energy recovery power at moment k + i that is predicted at moment k.
Since the motor braking torque is superimposed on the front wheels, which corresponds to a reduced slope of the β line, the front wheels are more likely to hold and the rear wheels are less likely to hold, so only the slip rate of the front wheels is optimized. w u , w s , and w P , the weights of each, are constants. P m is the braking energy recovery power that can be expressed as: In the model prediction process, the prediction time domain h p is usually longer than the control time domain h c , and inputs beyond the control time domain h c are still using the last input element of h c . Since the vehicle dynamics system is very sensitive to changes in inputs, inputs beyond h c are not helpful to the optimization of the system, so the prediction time domain in this paper is the same as the control time domain to ensure that there is an independent control input at each moment.

Design Constraints
During the braking process, the motor torque is also influenced by the motor external characteristics and the charging current limit of the power battery. In engineering, motor drives usually limit the motor torque based on the motor external characteristic curve but do not consider the effect of battery charging current limit. Therefore, when the control strategy is designed, if the above limitation is still not considered, it will cause damage to the power battery or cause a serious mismatch in the model. To avoid the above situation, the above limitation should be considered when designing the regenerative braking control strategy. So, the constraints on the state variables and inputs of the system can be expressed as: where T mmin (ω m ) is the maximum braking force of the motor on the external characteristic curve of the fourth quadrant; s max , the limit of the wheel slip rate, is a constant; and I max is the charge current limit of the battery, which varies with SOC and temperature. Since there is an order of magnitude difference between the time duration of a control moment and the time duration of the various time of I max , it can be considered as a constant when solving the optimization problem.

Build Simulation Model
The whole vehicle simulation model is established based on the Simscape toolbox in Matlab/Simulink to test the effectiveness of the above NMPC-based cooperative braking algorithm. The simulation model includes the driving cycle, longitudinal driver, environment, controllers, passenger car, and visualization modules.
The original controller module included only the powertrain controller. To verify the above NMPC-based cooperative braking algorithm, a brake controller with brake energy recovery was built independently, as shown in Figure 4.  The brake controller model includes a sub-model based on the NMPC cooperative braking algorithm, as shown in Figure 5. The input signal is the total braking demand torque, and the optimal motor braking torque is output through the model prediction and optimal solution procedure of NMPC algorithm, which can achieve a good energy recovery effect while ensuring braking safety, and at the same time, the vehicle speed and front The brake controller model includes a sub-model based on the NMPC cooperative braking algorithm, as shown in Figure 5. The input signal is the total braking demand torque, and the optimal motor braking torque is output through the model prediction and optimal solution procedure of NMPC algorithm, which can achieve a good energy recovery effect while ensuring braking safety, and at the same time, the vehicle speed and front and rear wheel speeds are fed back to NMPC algorithm for rolling optimization of the prediction model. The optimization objectives are the braking energy recovery power, tire slip rate, and motor torque variation, and the constraints are the motor torque limited by the motor external characteristic curve, power battery charging current, and wheel slip rate limits.

Simulation and Result Analysis Based on WLTC
The energy recovery effect of the proposed control algorithm will be tested based on WLTC, and the energy recovery effect will be evaluated using the braking energy recovery rate.

Select Driving Cycle
Currently, there are a variety of driving cycles used for vehicle energy efficiency testing in various countries around the world, including NEDC, FTP75, UDDS, and JC08. For China, the evaluation of the regenerative braking system started late and did not develop its test cycle conditions, usually using NEDC as the standard. NEDC is divided into two parts: urban conditions and suburban conditions. Urban conditions consist of four UDC conditions, suburban conditions for one EUDC condition. The target speed curve for NEDC is shown in Figure 6a.

Simulation and Result Analysis Based on WLTC
The energy recovery effect of the proposed control algorithm will be tested based on WLTC, and the energy recovery effect will be evaluated using the braking energy recovery rate.

Select Driving Cycle
Currently, there are a variety of driving cycles used for vehicle energy efficiency testing in various countries around the world, including NEDC, FTP75, UDDS, and JC08. For China, the evaluation of the regenerative braking system started late and did not develop its test cycle conditions, usually using NEDC as the standard. NEDC is divided into two parts: urban conditions and suburban conditions. Urban conditions consist of four UDC conditions, suburban conditions for one EUDC condition. The target speed curve for NEDC is shown in Figure 6a. The overall NEDC is relatively stable, and most of the time the vehicle is in a constant speed state. However, NEDC does not take into account the frequent stopping of vehicles in urban traffic jams, there is a large deviation from the actual road traffic in China. Therefore, in recent years, the conditions-testing standard in China is transitioning from NEDC to WLTC. WLTC is a standard jointly developed by Japan, the United States, and the European Union, which is obtained by processing the real driving conditions data collected worldwide. Compared with NEDC, it fully takes into account factors such as vehicle rolling resistance, gearing, and vehicle weight, and its speed fluctuates greatly, with fewer idling conditions and covering a wider speed range, without cyclical acceleration and deceleration, better reflecting the frequent changes in speed under different road congestion states. Therefore, the energy-saving evaluation effect obtained by using WLTC is more practical.
WLTC is divided into four parts: low speed, medium speed, high speed, and superhigh speed. Compared with NEDC, its test cycle increases from the 1180 s to 1800 s; the duration of each part is 589 s, 433 s, 455 s, and 323 s, respectively, and the average speed increases from 34 km/h to 46 km/h. WLTC has a longer test period and higher average speed, which is more in line with the actual driving conditions. The broader speed range for the overall performance of the vehicle also puts forward higher requirements. The target speed curve for WLTC is shown in Figure 6b.
Based on the above analysis, WLTC was selected for the simulation test in this paper.

Simulation and Result Analysis
The data from WLTC were input into the simulation model of regenerative braking system with NMPC-based cooperative braking algorithm for testing. To better verify the effectiveness of the proposed algorithm, the comparison test with the parallel control strategy was also performed. The comparison curve of target speed and actual speed is shown in Figures 7 and 8. Setting the initial SOC to 0.8, the battery SOC curve can be obtained, as shown in Figure 9. The motor braking torque curve and slip rate curve are shown in Figure 10 and Figure 11, respectively. The overall NEDC is relatively stable, and most of the time the vehicle is in a constant speed state. However, NEDC does not take into account the frequent stopping of vehicles in urban traffic jams, there is a large deviation from the actual road traffic in China. Therefore, in recent years, the conditions-testing standard in China is transitioning from NEDC to WLTC. WLTC is a standard jointly developed by Japan, the United States, and the European Union, which is obtained by processing the real driving conditions data collected worldwide. Compared with NEDC, it fully takes into account factors such as vehicle rolling resistance, gearing, and vehicle weight, and its speed fluctuates greatly, with fewer idling conditions and covering a wider speed range, without cyclical acceleration and deceleration, better reflecting the frequent changes in speed under different road congestion states. Therefore, the energy-saving evaluation effect obtained by using WLTC is more practical.
WLTC is divided into four parts: low speed, medium speed, high speed, and superhigh speed. Compared with NEDC, its test cycle increases from the 1180 s to 1800 s; the duration of each part is 589 s, 433 s, 455 s, and 323 s, respectively, and the average speed increases from 34 km/h to 46 km/h. WLTC has a longer test period and higher average speed, which is more in line with the actual driving conditions. The broader speed range for the overall performance of the vehicle also puts forward higher requirements. The target speed curve for WLTC is shown in Figure 6b.
Based on the above analysis, WLTC was selected for the simulation test in this paper.

Simulation and Result Analysis
The data from WLTC were input into the simulation model of regenerative braking system with NMPC-based cooperative braking algorithm for testing. To better verify the effectiveness of the proposed algorithm, the comparison test with the parallel control strategy was also performed. The comparison curve of target speed and actual speed is shown in Figures 7 and 8. Setting the initial SOC to 0.8, the battery SOC curve can be obtained, as shown in Figure 9. The motor braking torque curve and slip rate curve are shown in Figures 10 and 11          Then, the above-mentioned sets of comparison graphs are analyzed as follows. According to the speed curve, the NMPC algorithm can follow the target velocity better, so its control accuracy and response effect is higher than that of the parallel algorithm. The variation of SOC curve indicated that after finishing an WLTC cycle, the SOC dropped from 80% to 73.2% under the proposed NMPC algorithm. When compared to 68.5% under the parallel algorithm, the difference was 4.7%, so the NMPC algorithm consumes less power. The motor torque curve shows that compared to the parallel strategy, the NMPCbased motor braking torque is higher, up to 90 Nm, and the drive torque spikes are smaller, contributing to smoother motor operation and better comfort. Furthermore, according to the slip rate curve, the NMPC algorithm can better keep the slip rate within the safety interval.
The energy efficiency of the algorithm is also evaluated using the energy recovery rate , with the following equation:   Then, the above-mentioned sets of comparison graphs are analyzed as follows. According to the speed curve, the NMPC algorithm can follow the target velocity better, so its control accuracy and response effect is higher than that of the parallel algorithm. The variation of SOC curve indicated that after finishing an WLTC cycle, the SOC dropped from 80% to 73.2% under the proposed NMPC algorithm. When compared to 68.5% under the parallel algorithm, the difference was 4.7%, so the NMPC algorithm consumes less power. The motor torque curve shows that compared to the parallel strategy, the NMPCbased motor braking torque is higher, up to 90 Nm, and the drive torque spikes are smaller, contributing to smoother motor operation and better comfort. Furthermore, according to the slip rate curve, the NMPC algorithm can better keep the slip rate within the safety interval.
The energy efficiency of the algorithm is also evaluated using the energy recovery rate , with the following equation: Figure 11. The slip rate curve.
Then, the above-mentioned sets of comparison graphs are analyzed as follows. According to the speed curve, the NMPC algorithm can follow the target velocity better, so its control accuracy and response effect is higher than that of the parallel algorithm. The variation of SOC curve indicated that after finishing an WLTC cycle, the SOC dropped from 80% to 73.2% under the proposed NMPC algorithm. When compared to 68.5% under the parallel algorithm, the difference was 4.7%, so the NMPC algorithm consumes less power. The motor torque curve shows that compared to the parallel strategy, the NMPC-based motor braking torque is higher, up to 90 Nm, and the drive torque spikes are smaller, contributing to smoother motor operation and better comfort. Furthermore, according to the slip rate curve, the NMPC algorithm can better keep the slip rate within the safety interval.
The energy efficiency of the algorithm is also evaluated using the energy recovery rate η reg , with the following equation: where E reg is the recovered motor braking energy during the whole cycle, which can be obtained by integrating the charging power P c , and E brk is the total energy consumption during the whole cycle, which is the sum of the front wheels braking energy consumption E µ f , the rear wheels braking energy consumption E µr , and the motor braking energy consumption E µm . The expressions are as follows: where U c is the charge voltage, I c is the charge current, F µ f is the front wheels braking force, F µr is the rear wheels braking force, F µm is the motor braking force, and u a is the vehicle speed. Based on Equations (29)-(31), the calculation module is inserted in Simulink and the results are presented in Table 2. According to these simulation results, the proposed cooperative braking algorithm can more effectively recover the braking energy than the conventional parallel control strategy, and under the whole WLTC, the energy recovery rate reaches 30.4%, which is higher than 21.1% of the parallel algorithm.
Based on the above analysis, the NMPC algorithm is higher than the parallel algorithm in terms of control accuracy, comfort, safety, and energy saving, which proves the effectiveness of the proposed algorithm.

Simulation and Result Analysis Based on Single Braking Condition
To further verify the safety of the algorithm, a single high-intensity braking test was conducted, in which the vehicle was gradually accelerated to 100 km/h, followed by braking with a maximum braking intensity of 0.8, and the simulation result curves       According to the above curves, the slip rate is always kept within 15-25% during braking. When the braking intensity and the slip rate are high, the motor braking torque will be reduced rapidly, which can prepare for the possible wheel locking. Therefore, it can still ensure good braking safety with high braking intensity. According to the above curves, the slip rate is always kept within 15-25% during braking. When the braking intensity and the slip rate are high, the motor braking torque will be reduced rapidly, which can prepare for the possible wheel locking. Therefore, it can still ensure good braking safety with high braking intensity.

Conclusions
The NMPC-based cooperative braking control algorithm is proposed in this paper, with braking energy recovery power, tire slip rate, and motor torque variation as the optimization objectives, and online optimization of the coordinated distribution of motor braking and friction braking to achieve good energy efficiency while ensuring braking safety.
Using the whole-vehicle dynamics model and the braking energy recovery system model, which were built in Matlab/Simulink, when based on WLTC, the NMPC algorithm proved to be better than the conventional parallel algorithm in terms of control accuracy, comfort, safety, and energy saving. Meanwhile, the energy recovery rate was 30.4%, and the maximum motor braking torque was 90 Nm. Under single high-intensity braking, the slip rate was always within 15~25% and the motor quickly exited the braking process to prepare for the possible wheel locking, which further proves that the algorithm can achieve good braking safety.
Considering that the discrete model in the paper was obtained on the basis of the Eulerian method, which is the first-order integrator characterized by low accuracy, and that its stability can only be ensured by sufficiently small step-sizes, a comparative analysis with other higher-order numerical integration methods will be performed in future studies [30]. In addition, the study has only progressed to the model and simulation stage, but in a real application, there is also a need to consider the sensors and instrumentation required to measure the relevant data, including the delays imposed by each element, as well as an estimation of the total response time of the whole system, including the microprocessor. So, this will also need to be addressed in future studies.