Stochastic model predictive braking control for heavy-duty commercial vehicles during uncertain brake pressure and road profile conditions

When heavy-duty commercial vehicles (HDCVs) must engage in emergency braking, uncertain conditions such as the brake pressure and road profile variations will inevitably affect braking control. To minimize these uncertainties, we propose a combined longitudinal and lateral controller method based on stochastic model predictive control (SMPC) that is achieved via Chebyshev–Cantelli inequality. In our method, SMPC calculates braking control inputs based on a finite time prediction that is achieved by solving stochastic programming elements, including chance constraints. To accomplish this, SMPC explicitly describes the probabilistic uncertainties to be used when designing a robust control strategy. The main contribution of this paper is the proposal of a braking control formulation that is robust against probabilistic friction circle uncertainty effects. More specifically, the use of Chebyshev–Cantelli inequality suppresses road profile influences, which have characteristics that are different from the Gaussian distribution, thereby improving both braking robustness and control performance against statistical disturbances. Additionally, since the Kalman filtering (KF) algorithm is used to obtain the expectation and covariance used for calculating deterministic transformed chance constraints, the SMPC is reformulated as a KF embedded deterministic MPC. Herein, the effectiveness of our proposed method is verified via a MATLAB/Simulink and TruckSim co-simulation.


Introduction
Advanced vehicle control features, including traction control, anti-lock braking systems, electronic stability control, adaptive cruise control, lane-keeping assistance, etc., have been topics of intense study over many years, and a variety of such functions have already been implemented to assist drivers. These and even more advanced driving automation functions are currently being readied for application to passenger cars, buses, and large vehicles whose handling B Kazuma Sekiguchi ksekiguc@tcu.ac.jp Ryota Nakahara 2081037@tcu.ac.jp 1 Graduate School of Integrative Science and Engineering, Tokyo City University, Tamazutsumi, Setagaya-ku, Tokyo 1588557, Japan 2 Nabtesco Automotive Co., Ltd., Hirakawacho, Chiyoda-ku, Tokyo 1020093, Japan requires superior driving skills. In particular, in emergency situations, advanced braking control is required to regulate the vertical vehicle control of vehicles of buses or trucks with large inertia and bring them to safe stops. Today, most buses and trucks are equipped with pneumatic braking systems that use compressed air as their energy medium. However, pneumatic braking systems have characteristics that make control design difficult [1][2][3][4]. First, large time delays occur due to the compressibility of air, and the dynamics of pneumatic braking systems inevitably mean that the relationship between the air medium and pressure flow is non-linear. Hence, when combined with the vehicle's longitudinal dynamics, pneumatic braking systems are prone to various uncertainties. These include supply pressure variations due to brake application and release, temperature increases resulting from frequent braking, brake system component wear, vehicle load fluctuations, and changes in road surface conditions due to rain or snow.
When considering such factors, model predictive control (MPC) [5], which considers dynamic coupling [6,7], topo-graphical information [8], and horizontal vehicle dynamics, provides an optimal control method for safe vehicle control during emergency stop conditions. For that reason, MPC is incorporated into a wide variety of automatic control systems such as train pneumatic braking systems [9], bus braking energy efficiency systems [10], and the steering systems of all-terrain cranes with slow dynamic response levels [11]. The high utility of MPC is its ability to solve the constrained optimization problem of a finite-time interval for each sampling step based on the model to be controlled. Then the controller applies the first component of the calculated input series in each sampling step. MPC makes it possible to calculate control inputs that explicitly consider functional, physical, and safety restrictions by applying constraints for input and state. However, MPC is vulnerable to the influence of actual motion factors that are not considered in models. Therefore, although preset constraints may be satisfied in MPC predictions, those constraints may be exceeded in actual conditions. MPC-based control has also been applied to heavy-duty commercial vehicle (HDCV) wheel slip control [12,13], fuel consumption systems [14][15][16][17], automatic hydraulic transmission systems [18], actuator time delay controllers [19], pneumatic actuator controllers [20]. The challenge of the controller design of HDCVs are robustness against model uncertainty due to actuator or loading. The methods proposed in [21] is not considered against the uncertainty of the plant model. Also, the Kalman filer (KF) embedded MPC brake system [22] is available to use state covariance for the controller. To solve that, robust antislip control systems for different road surface environments [23] and systems that compensate for vehicle specification changes caused by load variations [24,25] is reported. How-ever, the deterministic constraints result in aggressive inputs when the model has stochastic uncertainty. This may result in the model not complying with the constraints as theory suggests.
In contrast, stochastic model predictive control (SMPC) is a robust control method that considers such uncertainties. In the SMPC approach, which is shown in Fig. 1, we can see that SMPC inputs a stochastic constraint called a "chance constraint" based on the probability distribution of the prediction error and designates it as the upper and lower limits of the state. This allows SMPC to explicitly describe the probabilistic uncertainties in designing the robust control strategy. For example, let the deterministic constraint be x ≤ x max , with x as the states and x max as the upper limit states. If x is uncertain, the probability will not reach fifty percent. To avoid creating too large a value, the deterministic constraint needs to converted into a probabilistic (also called chance) constraint, in which Pr(x ≤ x max ) ≥ 1−α, with α ∈ (0, 0.5] is the allowable probability. In many applications, system uncertainties are often considered to be stochastic processes, which allows the incompleteness and uncertainties of set optimization problems to be taken into consideration. This suppresses erroneous control that can result from discrepancies between theory and an actual environment. SMPC can also improve both robustness and control performance using statistical information related to disturbances. For that reason, SMPC is also applied to uncertainties in automobile control [26][27][28]. On the other hand, external disturbances due to road environments, such as gradients and road surface friction coefficients, affect dynamic vehicle behaviors during emergency stops [22]. In such cases, wheel slip suppression Fig. 1 Overview of SMPC approach. SMPC explicitly describes the probabilistic uncertainties that need consideration when designing a robust control strategy using a probabilistic (also called chance) constraint. Here, Pr(x ≤ x max ) ≥ 1 − α, with α ∈ (0, 0.5] is the allowable probability improvements can be anticipated by considering the vehicle uncertainties, such as the brake pressure and road profile.
To facilitate this, this paper proposes a braking control method that considers the road profile outlined in Reference [29] as a stochastic uncertainty, and via which SMPC is applied to the control of HDCVs during emergency braking. The results of this study show that our proposed method improves braking system robustness against fluctuations in temperature, air brake system measurement values, and road surface features, all of which are mathematical programming problems involving random variables that are generally considered challenging to solve. In our method, we reformulate an SMPC into a deterministic equivalent model based on the expectations and covariances of random variables, which we achieve using the KF algorithm. This allows HDCV uncertainties to be controlled via the SMPC by considering disturbance-related statistical information. Hence, the primary contribution of this paper is the proposal of a braking control formulation that is robust against the probabilistic friction circle uncertainties affect. The remainder of this paper is organized as follows. The vehicle model, including an HDCV airbrake system, is described in Sect. 2, while road profiles are discussed in Sect. 3. Next, the KF-based state estimation is described in Sect. 4, while the SMPCbased controlled design is discussed in Sect. 5 and numerical simulation results are shown in Sect. 6. Finally, Sect. 7 concludes the paper.
Notations This paper adopts the following notations. The subscripts fl, fr, rl, and rr refer to front left, front right, rear left, and rear right, respectively. A := B refer to define A as B. N ∼ N (a, b) refer to the random variable N follows a normal distribution N (a, b) with standard deviation a and variance b. x P := √ x T P x represents the weighted norm. R represents set of real numbers. Additionally, A O, A O denotes A is a positive definite matrix and semi-positive definite matrix, respectively. Pr(·) represents probability, R n represents an n-dimensional real vector, R n×m represents an n × m real matrix, 0 n is n-dimensional zero vector, and O n×m represents n × m zero matrix. The subscript k + i | k refer to the predicted/estimated value of the variable for the time k + i, which is given at time k.

Vehicle model
This section describes the equation of motion of HDCV's kinematic and dynamic maneuvers in the horizontal X -Y plane, as shown in Fig. 2. To facilitate its adaptation as a control model, some assumptions are simplified. The vehicle kinematics model is expressed as follows: where γ is the yaw rate, and v x and v y are, respectively, the longitudinal and lateral velocities in the vehicle frame. Let β ref = 0, reference vehicle kinematics is expressed as follows: where ρ ref is road curvature. To incorporate lane-keeping, it is necessary to transform the lateral vehicle dynamics into path tracking error dynamics. In such cases, the relative errors between the ego and reference vehicles can be expressed as follows: Next, assuming the vehicle side-slip angle β and the front-wheel steering angle δ are small, the equations for the longitudinal, lateral, and yaw dynamics in the vehicle frame are expressed, respectively, as follows: where M is total vehicle mass, l f and l r are the distances from the center of gravity (CoG) to the front and rear axles, I z is the vehicle body moment of inertia about the z-axis, and w f and w r are the widths of the front and rear track, respectively. In addition, F xi j and F yi j are the longitudinal and lateral tire forces, s d is the vehicle driving distance, and e y and e ψ denote the heading error and the lateral deviation values, respectively.

Longitudinal dynamics
When the slip ratio of each wheel is small (|λ| 1), brake torque is approximated as follows: where r w is the effective wheel rolling radius. Since this paper focuses on emergency vehicle stops, it deals solely with braking and does not consider driving-related aspects. We also assume that the pressure-braking force characteristics are linear, as outlined below: where k b is the proportional relationship constant of the brake pressure to the braking force characteristics. Since good mechanical lubrication is assumed, we ignore the effects of mechanical and viscous actuator friction. Hence, from (4), (5) and (6), the longitudinal dynamics can be expressed as follows: with x lon = s d V T , u lon = P fl P fr P rl P rr T ,

Lateral dynamics
Tire side-slip angles α i are expressed as follows [30,31]: where l f and l r are the CoG distances to the front and rear axle, respectively, and w f and w r are the front and rear track widths, respectively. Assuming |δ|, |β|, |α f |, and |α r | 1, α i is approximated as follows: Next, assuming that the characteristics of the left and right tires are equal and that lateral tire forces F yi are approximated as linear functions of the tire side-slip angles α i , the lateral tire forces F yi are approximated as follows: where K yf and K yr are the cornering stiffness values. From (4) and (10), the lateral dynamics in the vehicle-fixed frame are expressed as follows:

Air-brake model
The air-brake system structure is shown in Fig. 3. Assuming good mechanical lubrication, we ignore the effects of mechanical and viscous friction of actuators. This system consists of an air tank, a brake chamber, a solenoid, and a valve. The brake chamber is supplied with compressed air through the relay valve. We assume that the supply source is always filled with compressed air. The air pressure is adjusted by opening and closing the valve by the command voltage applied to the solenoid.
In this paper, we assume that a time delay occurs due to air propagation and that the time delay is modeled as a first-order delay system as follows: x air = A air x air + B air u air , with x air = u lon , u air = u fl u fr u rl u rr T , where τ i and K Pi are the time constant and gain of the airbrake system, respectively, which parameters are same as [21]. u i is the command voltage, and P i j is the brake pressure.

Vehicle load transfer
The formula used for the vehicle load transfer estimate is as follows [32]: with , where L = l f + l r is the wheelbase, m s is the vehicle sprung mass, m uf and m ur are the vehicle front and rear spring mass, respectively, g is the acceleration due to gravity, h g is the CoG height, and r w is the effective tire rolling radius.

Road profile
Because there is such a wide variety of road surface environments, including asphalt, cobblestones, and unpaved ground, the International Organization for Standardization (ISO) has classified road surface conditions into eight levels (from A to H) in the ISO 8608 [33] standard. This standard also classifies road surface height z g (s d ), which is expressed as follows: with where A n is the complex Fourier coefficient, f 0 = 1/(2π), and f n is the spatial frequency, S( f n ) is the power spectral density (PSD) function, and θ ∼ U(0, 2π) is the phase with a uniform distribution. We consider the finite partial sum up to the N term for the infinite series of (16) [34]. First, from the Fourier spectrum, Ψ ( f n ) of S( f n ) is defined as follows: where Ψ * is the conjugate complex number of Ψ , and the Fourier spectrum Ψ is the effective value of the complex Fourier coefficient A n , i.e., Ψ = A n / √ 2. Accordingly, if we let L road be the overall road length and assume the bandwidth is f n ≈ n f = 2π/L road , S( f n ) can be approximated as follows: Then, we substitute (18) into (16), the finite sum up to the N term can be expressed as follows:   Table 1. The A and B classes correspond approximately to asphalt, while the D-E classes correspond to irregular roads. The road profile (with respect to the driving distance represented by (16)) is shown in Fig. 4, where L road = 250 m, and N = 2500. Increased k values indicate increased vertical height. Next, we will examine the distribution of road profile characteristics. The D-E class road profile histogram is shown in Fig. 6. To compare similarities between normal and sample distributions, the quantile-quantile (Q-Q) plot of the D-E class road profile is shown in Fig. 5 where the horizontal axis is the quantiles of the standard normal distribution, and the vertical axis is the sample quantiles. As the quantiles increase in similarity, the tendency of the blue dots to follow the straight line indicated by the red alternating long-short dashed line increases. The green line shows the interquartile. As can be seen in Fig. 5, the road profile characteristics differ from the normal distribution.

Kalman filter
In this paper, the output is expressed as follows: and the unobtainable state is estimated by KF. First, the estimation model discretized (15) by the zero-order hold at the sampling time t is expressed as follows: where A d , B d , and E d collectively define the coefficient matrix of a discrete-time system, w ∼ N (0, w ) and v ∼ N (0, v ) are the process and observation noises, respectively, and w O and v O are the covariance of the process and observation noises, respectively. Then, the prediction steps are calculated aŝ wherex k|k−1 and P k|k−1 are the a priori state estimate and error covariance, respectively. The filtering steps are calculated as where M k is the Kalman gain, andx k|k and P k|k are the a posteriori state estimate and error covariance, respectively.

Stochastic MPC-based vehicle control
This section describes the formulation of the optimization problem designed to achieve vehicle stability. First, we describe the cost function required to balance the input and state priority necessary to achieve the desired performance, after which we describe the constraints. In particular, we examine how uncertainty is considered via the chance constraint in the optimization problem.

Cost function
In this paper, the evaluation function of the following equation is defined as follows: where u k is the current time input,

Constraints
By designing a brake pressure P i constraint that satisfies the friction circle F 2 xi + F 2 yi ≤ μ i F zi , vehicle slip occurrence is suppressed. Here, μ i and F zi are the road surface friction coefficient and vertical force, respectively. From the relational expression of the circle of forces and (6), the limit value of the brake pressure P i is expressed by the following inequality: In addition, the speed constraint for not taking into account backward movement, the upper and lower limits for the command voltage of the air brake system, and the lateral deviation e y constraint for preventing lane departure are set as follows: Next, we will explain how to transform the chance constraint into a deterministic constraint to solve the stochastic optimization problem as an equivalent deterministic problem [35]. First, let g j be a jth basis vector, and let the jth component of the state be given by x j = g T j x k . Then, the chance constraint of x j is where η j := 1−α j and α j is allowable probability. Since we focus solely on design brake pressure chance constant here, α j is expressed as follows: Note that α = 1/2 is equivalent to the deterministic constraint. As can be seen in Fig. 6, the road profile is an asymmetrical and multimodal distribution, so an approximation as a Gaussian distribution is not reasonable. Therefore, from the Chebyshev-Cantelli inequality [36], the expected value E(x k ) =x k|k , and the variance σ 2 j = g T j P k|k g j , the chance constraint is transformed as follows: where x margin j Even though it is a conservative constraint, it is possible to use the Chebyshev-Cantelli inequality to deal with general distributions.

Optimization problem
Using a vehicle model, a cost function, and the relevant constraints, a deterministic equivalent constrained stochastic programming is expressed as follows: Deterministic equivalent SMPC: To avoid infeasibility, the state and input constraints are relaxed by the slack variable [37]. The relaxation degree is adjusted by Equal Concern for Relaxation (ECR) v max •,ECR . Note that v max •,ECR = 0 refers to a hard constraint that is relaxed as the value of the ECR increases.

Simulation result and discussion
The control system is shown in Fig. 7. Here, it can be seen that KF estimates the state acquired from the controlled object. The SMPC control input is calculated using the target value, a measurable disturbance, and the estimated value. The control target is a four-axle eight-wheel TruckSim truck model provided by Virtual Mechanics. In this study, we used the default Magic Formula (MF) model in TrcuckSim as our tire model. We also considered it difficult to identify the MF parameters from time to time due to the various changes in road conditions. Therefore, in our control model for MPC, we assume that the tire force is linear for the longitudinal and lateral slip of the tire as mentioned in Sects. 2.1 and 2.2. In SMPC, this  is considered to be a two-axle four-wheel vehicle model with l f = (l 1 + l 2 )/2, l r = (l 3 + l 4 )/2, where l i is the distance from the front to the i th axle and the CoG. The truck model inputs are the brake torque T i and the handle angle δ SW . The brake torque is calculated by (5) and distributed to each axle as T 1,2 = T f /2, T 3,4 = T r /2. The steering wheel angle δ SW is calculated by δ SW = n gear δ, where n gear is the steering gear ratio. The vehicle specifications used in the equation of state are listed in Table 2, while the SMPC control parameters are listed in Table 3 and the estimated KF parameters are as shown in Table 4. Chance constraints are designed for each of the four-wheel brake pressures, each of which had the same allowable probability α. In this study, the value of the speed weights is set to a large value to stop the vehicle quickly. In this way, the weights are designed to emphasize tracking to the target speed. In the driving scenario shown in Fig. 8, the road surface is a split-μ circular route with μ on the left and right sides in the direction of travel set at 0.6 and 0.9, respectively, and be known in this study. This road surface information is assumed to be known. The turning radius of the circular path is r = 152.4 m. Assuming emergency braking conditions, the target speed V ref is changed from 70 km/h to 0 km/h at time 2. Additionally, the target state is Hence, the left and right road surfaces are different, and the longitudinal and lateral force act at the same time. The control performance validity in this scenario is verified through emergency braking. Note that the brake pressure noise follows a normal distribution of N (0, σ 2 f,r ). In this study, braking performance verification is conducted on roads with uneven surfaces, as shown in the D-E class of Fig. 4 and be unknown in this study. The MPC results are also shown for comparison purposes. The SMPC/MPC iteration max is set to 50.
In this figure, the black-filled area shows the time evolution after the vehicle has stopped. A comparison between SMPC and MPC model results during braking on a split-μ road with a vertical road profile displacement is shown in Fig. 9. The red and blue solid lines show the cases of SMPC and MPC, respectively. The green dashed line shows the ref-     Fig. 9a, while Fig. 9b shows the time evolution of the vehicle longitudinal speed. The lateral deviation time evolution is shown in Fig. 9c. The chain line shows constraints. The time evolution of QP iterations is shown in Fig. 9d. Here, the yellow line shows that using the QP is infeasible, while the chain line shows max iterations. The time evolution of traveling distance is shown in Fig. 9e. The Key Performance Indicators (KPIs) is shown in Table 5 and Fig. 9f, where index factor is max deceleration, max corrective steering angle, max traveling distance, mean front wheels slip ratio, mean The yellow line shows that using QP is infeasible, while the chain line shows the max iterations rear wheels slip ratio, max yaw rate [13]. The time evolution of execution time is shown in Fig. 10, while Fig. 10a, c and Fig. 10b, d show the SMPC and MPC cases, respectively. The time evolution of brake pressure in the SMPC case is shown in Fig. 11, while Fig. 11a and b show the SMPC and MPC cases, respectively. Additionally, the slip ratio of each wheel in the SMPC case are shown in Fig. 12, while Fig. 12a shows the slip ratio of first and second axle wheels. Next, Fig. 12b shows the slip ratio of third and fourth axle wheels, while Fig. 13 shows the slip ratio of each wheel in the MPC case. The vehicle braked while maintaining lane tracking in each case, as shown in Fig. 9a. As shown in Fig. 9c, the lateral deviation constraint was violated in both the SMPC and MPC cases because it was necessary to allow the slack variable constraint to be violated to make the chance constraint of the friction circle feasible. In addition, the lateral deviation in the SMPC case was smaller than in the MPC case because the solution became infeasible in the MPC case when its optimality was lost at the 2-s mark, as shown in Fig. 9d. Hence, when considering the statistical information related to disturbances, the results of this study show that lane followability and running stability during braking on a split-μ circular route are improved by the use of SMPC. In Fig. 9b, we can see that the SMPC simulation stopping time was slightly longer than the MPC stopping time. However, As shown in Fig. 9e, stopping distance of SMPC is 71.8 m and MPC is 72.2 m. Also, as shown Table 5 and in Fig. 9f, SMPC KPIs are improved compared with MPC KPIs. The slip ratio during braking in the SMPC case, as shown in Fig. 12, was smaller than the MPC case, as shown in Fig. 13. This is because SMPC considers the confidence interval of the brake pressure to be a chance constraint and is more conservative than the original constraint, as shown in (31) and Fig. 11. As a result, slippage due to uncertainty is suppressed, as shown in Fig. 12. In the MPC case, the brake pressure uncertainty is not explicitly considered, so the ratio of time exceeding the tire friction circle constraint is higher than that of SMPC, which means the vehicle slipped more, as shown in Fig. 13. The slip ratio values in the black-filled area oscillated after the vehicle stopped as shown in Figs. 12 and 13. This is because the formula for calculating the vehicle skid angle has a term that includes velocity in the denominator, and the calculation becomes unstable due to division by zero as the vehicle is braked. From Fig.10, the longer the prediction horizon H p is, the more dimensionality of the matrices handled in the  The black-filled area shows the time evolution after the vehicle has stopped. a First and second left/right wheels, b third and fourth left/right wheels optimization calculation. Therefore, the average and peak execution times are longer. In all cases of H p = 10, 7, 5, the computation time peaked at 2 s, which is the start of stopping, but it was below the sampling time of MPC t = 100 ms. This is because the vehicle switched from constant driving, which does not reach the friction circle constraint, to braking, which does reach the friction circle constraint. As a result, the active constraints in the effective constraint method were updated, and the Hesse matrix of the Lagrangian function was re-evaluated. This resulted in a temporary increase in the run time by 2 s.

Conclusions
This paper proposed a method of SMPC that considers the uncertainty of brake pressure and road profiles in relation to the emergency braking of heavy-duty commercial vehicles (HDCVs). Herein, the proposed method was verified through a simulation involving turning braking on a split-μ road. The running stability of an HDCV ensured stable conditions in areas where the controlling braking force and the lateral force interacted.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.
Hiroki Hasebe received the B.Sc. and M.Sc. degrees from the Tokushima University, Tokushima, Japan, in 1999 and 2001, respectively. He is currently work at Nabtesco Automotive Co., Ltd., Kanagawa, Japan.
Kenichi Matsubara received the B.Sc. degree from the Muroran Institute of Technology, Hokkaido, Japan, in 1991. He is currently work at Nabtesco Automotive Co., Ltd., Kanagawa, Japan.