Real-Time Predictive Control of Path Following to Stabilize Autonomous Electric Vehicles Under Extreme Drive Conditions

A novel real-time predictive control strategy is proposed for path following (PF) and vehicle stability of autonomous electric vehicles under extreme drive conditions. The investigated vehicle configuration is a distributed drive electric vehicle, which allows to independently control the torques of each in-wheel motor (IWM) for superior stability, but bringing control complexities. The control-oriented model is established by the Magic Formula tire function and the single-track vehicle model. For PF and direct yaw moment control, the nonlinear model predictive control (NMPC) strategy is developed to minimize PF tracking error and stabilize vehicle, outputting front tires’ lateral force and external yaw moment. To mitigate the calculation burdens, the continuation/general minimal residual algorithm is proposed for real-time optimization in NMPC. The relaxation function method is adopted to handle the inequality constraints. To prevent vehicle instability and improve steering capacity, the lateral velocity differential of the vehicle is considered in phase plane analysis, and the novel stable bounds of lateral forces are developed and online applied in the proposed NMPC controller. Additionally, the Lyapunov-based constraint is proposed to guarantee the closed-loop stability for the PF issue, and sufficient conditions regarding recursive feasibility and closed-loop stability are provided analytically. The target lateral force is transformed as front steering angle command by the inversive tire model, and the external yaw moment and total traction torque are distributed as the torque commands of IWMs by optimization. The validations prove the effectiveness of the proposed strategy in improved steering capacity, desirable PF effects, vehicle stabilization, and real-time applicability.


Introduction
The rapid development of vehicle-to-everything technologies provides more opportunities to promote vehicle intelligence, enabling multiple-source information exchange to Academic Editor: Weichao Zhuang. improve vehicle safety, traffic mobility, and fuel economy [1]. To support these applications, the path following (PF) of autonomous vehicles is a crucial and non-negligible implementer, i.e., bridging the higher-level decision-making and lower-level vehicle dynamics. The prime target of PF is to effectively track the reference path as close as possible to the premise of vehicle stability. This is not a facile task since the control effects are sensitive to the vehicle dynamics, especially under extreme driving conditions [2]. The force features of tires under extreme cycles exhibit high nonlinearity, which potentially induces the risks of vehicle instability. This also leads to design difficulties in controller developments.
To prevent vehicle instability, the first question is how to evaluate the vehicle stability status. To this end, the phase plane portrait approach is widely adopted, which can intuitively yield the variants of vehicle motion phase trajectories by giving different initial states and driving conditions. In general, two types of phase planes are mostly studied in the last decade, i.e., sideslip angle-sideslip angle rate ( -̇ ), and sideslip angle-yaw rate ( -). In Refs. [3,4], the -̇ phase plane was analyzed, and the diamond vehicle stability region is defined by two saddle points and two calibrated interactions of zero sideslip angle. In Refs. [5,6], the two-lines method was developed to determine the stability region in -̇ phase planes. To avoid potential tire skids, the yaw rate limits of the maximum lateral acceleration are considered in -̇ phase planes, and a modified stability region based on the two-lines method was proposed in Ref. [7]. The above approaches are verified to be effective in assessing vehicle stability, but entailing enormous work burdens for vertices calibrations or polynomial fittings. Unlike -̇ phase plane portraits, the -one is more accurate and easier to operate since the stable bounds can be directly expressed by the limits of yaw rate and tires' sideslip angle. In Refs. [8,9], the parallelogram envelope is proposed to define the vehicle stability region, which is derived from the maximum values of yaw rate and sideslip angle of rear tires. Based on the parallelogram envelope, Refs. [10,11] concluded that the bounds regarding limited yaw rate, and front and rear tires' sideslip angles, will meet at one saddle point. Hence, the steering angle envelope can be deduced for vehicle stability assessments. Further, Ref. [10] points out that the above parallelogram envelop will cause small overshoots of yaw rate when converging to the equivalent point during steadystate steering. To solve that, the limited front and rear tires' sideslip angle, as well as yaw acceleration nullclines, are combined to develop the stable region by -phase plane analysis. Although the above research is successful, there still exists a deficiency for PF issues under extreme conditions. That is their limited yaw rate bounds, derived from maximum lateral acceleration, are set up under the assumptions of zero lateral velocity differential. The zero-sideslip angle is generally pursued in stability control, making zero lateral velocity differential, but it is not always helpful for PF control under extreme conditions. For instance, during vehicle turning in high velocity and low tire-road friction, the reasonably increase/decrease of lateral velocity is beneficial to minimize PF errors, inducing the nonzero lateral velocity differentials. If the yaw rate is restricted by the bounds with zero lateral velocity differential, the vehicle steering capacity may be confined, entailing the deteriorations in PF effects.
After defining the vehicle stability regions, the outstanding PF controller should be developed. There are a series of PF control approaches proposed in the last decade [12,13], e.g., fuzzy logic control [14], sliding mode control (SMC) [15], robust H ∞ control [16], and nonlinear Lyapunov control [17]. To avoid the disturbances on path curvature, a vision-based lane keep controller is developed by two proportional-integral nested loops [18]. In Ref. [19], an exponential-like-sliding-mode fuzzy type-2 neural network strategy was designed, which is verified to yield the desirable convergence and robustness. In Ref. [15], considering the system uncertainties, a sliding mode control (SMC)-based composite nonlinear feedback controller was proposed to reject the lumped disturbance. In Ref. [20], the Takagi-Sugeno fuzzy technique was adopted to consider the time-varying vehicle velocity in modeling, and output feedback gainscheduled PF strategy was proposed that can improve the closed-loop transient performance by D-stability concept. To handle the uncertainties of impaired driver behaviors, a linear-parameter-varying based robust H ∞ strategy was designed for PF of driver assistance control system [21]. For active steering and direct yaw moment control (DYC) issues, an integrated controller was developed by Lyapunov nonlinear method, which actualized the function-decoupling control and guarantees the closed-loop system stability [17].
However, the aforementioned research mainly focuses on transient performance or robustness but neglects constraint handling performance. There are various system constraints in real-world PF applications, such as limits of vehicle stability, actuator, and driving area. Often, they are ignored in the controller design phase and then introduced through post-process operations (like imposing saturation limits regarding inputs and/or control barrier function [22]), while this may deteriorate the control optimality. Therefore, it is significant to think about the constraints during controller development.
Alternatively, model predictive control (MPC) is an effective and appropriate choice. Unlike other methods, in MPC, the constraints can be easily considered by explicit expressions, and the system states in the predictive horizon can be optimized by algorithms, which guarantees control optimality [23,24]. In Ref. [11], an MPC control strategy based on a linear single-track model was proposed for PF and vehicle stabilization. In Ref. [25], the linear-time-varying (LTV) MPC was developed for active steering control, while with poorer performance than nominal nonlinear MPC (NMPC). Indeed, considering the high nonlinearity of tires, the NMPC is preferable to linear MPC but suffers from huge calculation burdens. To this end, the explicit MPC was proposed in Ref. [26]. It is verified to show a similar control performance with nominal NMPC in the premise of a smaller grid, but leading to high memory requirements that limit its wider application. Benefited from the parallel process capacity of field-programmable gate array (FPGA) chip, the NMPC by particle swarm optimization algorithm is validated to be real-time applicable under hardware-in-the-loop test [27]. That said, it is hard to have a large-scale application due to the high expense of FPGA chips.
Except for the real-time capability, the closed-loop stability is also a crucial property of the control system. Because the control commands in MPC are determined by optimization, they are implicit in the system state, bringing huge difficulties in analyzing the closed-loop stability, in particular for nonlinear systems. The popular methods are to first linearize the system and then introduce the terminal constraints and the terminal cost [28]. However, as mentioned previously, the linearized MPC is hard to take the nonlinearities of vehicle dynamics into account, entailing poorer performance.
To fill the above gaps, this paper proposes a real-time NMPC strategy for desirable effects of PF and vehicle stability. The distributed drive electric vehicle (DDEV) is selected as the study object in this paper, which brings more control degrees of freedom (DOF) for improved effects but with greater complexities [29]. The main contributions are summarized: (1) For the expected PF effects and application requirements, an accurate but not too complicated model for controller design is indispensable, which is generally named a "control-oriented model". In this paper, an accurately nonlinear control-oriented model combined with the Magic Formula tire model and the singletrack vehicle model is established to feature the lateral dynamics of vehicles. Unlike lots of previous research [30] that apply front steering angle, the lateral force of front tires is adopted as the control variable, to reduce control complexity but reserve the front tires' nonlinearity. (2) The lateral velocity differential is considered in the maximum lateral acceleration, and the phase plane analysis of lateral velocity and yaw rate are carried out. It is found that for constant steering angle and tire-road friction, the maximum lateral acceleration is varied with vehicle states and markedly different from the yaw rate bounds without considering lateral velocity differential. Moreover, there exists a steering angle under critical stability cases, and meanwhile, the yaw rate of the maximum lat-eral acceleration will converge to that without considering the lateral velocity differential. Given the phase plane analysis, the stable bounds of the steering angle are proposed and transformed as front tires' lateral forces according to the vehicle model. The proposed bounds can effectively guarantee vehicle stability. Owing to the consideration of lateral velocity differential in the maximum lateral acceleration, it can also improve the vehicle steering capacity, which is very helpful for PF control under extreme conditions. (3) To guarantee the closed-loop stability in PF, the Lyapunov-based constraint is developed for the NMPC strategy. Its sufficient conditions for recursive feasibility and closed-loop stability are explicitly derived. (4) The C/GMRES algorithm is proposed for real-time optimization in the NMPC strategy. To handle the inequality constraints, the relaxation function method is adopted to transform the inequality constraints as the equivalent cost in the optimization objective. The proposed C/GMRES algorithm-based NMPC strategy is verified to implement superior control effects computationally.
The rest of this paper is organized as follows. Section 2 describes the control-oriented model. Section 3 introduces the control framework of the proposed strategy. Section 4 illustrates the C/GMRES algorithm-based NMPC controller. The inversive tire model's application and torque allocation method are shown in Sect. 5. Section 6 gives the validations and analysis, followed by the conclusions in Sect. 7.

Control-Oriented Model Formulation
The studied vehicle is an autonomous DDEV with frontsteering and four-wheel independently-driven functionalities. Before modeling, some assumptions are arranged: (1) the vehicle velocity v x , the lateral velocity v y , the yaw rate , and the tire-road friction factor , can be accurately observed or estimated by advanced sensors/algorithms; (2) assuming that the rotational inertia of wheels is small, there is F x ≈ T∕r w , where F x and T are the longitudinal force and the torque for one tire and related IWM, respectively, and r w is the tire effective radius.
The schematic of the PF issue is depicted in Fig. 1. Defining the lateral position error as y e and heading error as e [31]: where and r are the vehicle heading angle and its reference value, respectively. To determine the appropriate (1) ̇y e = v y e + v ẋ e =̇−̇r = −̇r waypoints, a preview method in the previous research [32] is employed, which can violate the scene "take a shortcut" and contribute to the expected effects.
The single-track vehicle model is applied to feature vehicle motions, as the schematic in Fig. 2: where m is the vehicle mass. l a and l b are the distances from front axle and rear axle to body center of gravity (CG), respectively. I z is the body yaw moment of inertia. F yf and F yr are the total lateral tire force of front and rear wheels, respectively, which can be represented by F yf = F yf l + F yfr F yr = F yrl + F yrr .
In this paper, the subscripts of "fl", "fr", "rl", and "rr", indicate that the corresponding variables are related to left front, right front, left rear, and right rear wheels, respectively. Unlike traditional powertrains, IWMs in DDEV can independently provide the traction torque at each wheel, so that the external yaw moment ΔM z is produced surrounding the vertical axle of vehicle CG to govern yaw motion, as shown in the second equation of Eq. (2). The sideslip angles of front and rear tires, i.e., f and r , are expressed: where is the steering angle of the front wheels. To characterize the tire lateral forces, the Magic Formula model with similarity theory under pure-slip conditions [33] is employed: where F z• is the vertical force of one wheel. B y is the stiffness factor, and D y is the peak factor. Here D y = To sum up, the control-oriented model can be sorted as where x = y e , e , v y , T is t he st ate var iable, u = F yf , ΔM z T is the control variable, and d =̇r is the external disturbance. y is the system output. C is the output matrix and equals to diag{[1, 1, 0, 0]} . In this paper, v x is fed back from the vehicle system. Remark 1: Rather than as the control variable for general models, F yf is chosen, and command can be acquired by the inversive tire model (as introduced in Sect. 5). In this manner, the high nonlinearity from couples between F yf and is wiped out, which lowers the model complexity and reduces the calculation loads in the controller.
Remark 2: Due to the study emphasis in this paper, the longitudinal motion is arranged as a constant velocity tracking issue. That is the velocity variant of each tire will not be greatly changed, resulting in a small tire slip rate. Hence, the Magic Formula model under pure-slip conditions can represent the tire lateral force with acceptable accuracy.

Control Framework
The control framework of the proposed strategy is illustrated in Fig. 3. It is assumed that there is a higher-level path planner that can generate the desirable path based on the road environment, and the path information is applicable. At each control instant, the path information, stable bounds of front tires' lateral force, Lyapunov-based stability constraint, and feedback states, are provided to the NMPC controller. In the controller, the relaxation function is applied to transform the original problem into an unconstrained one for the inequality constraint handling. The optimization problem is solved by the C/GMRES algorithm in a computationally efficient way, which outputs the expected front tires' lateral force and external yaw moment. With the target total front lateral force, the steering angle command is acquired by the inversive tire model and inputs into the vehicle plant. A torque allocation method based on Karush-Kuhn-Tucker (KKT) conditions is also employed to gain the expected IWMs' torque commands, to minimize tire workload and satisfy virtual torques (i.e., total traction torque and external yaw moment). Here a proportional-integral controller is applied to output the total traction torque for longitudinal velocity tracking, which is widely used and verified to be effective [7].

Nonlinear Model Predictive Control Problem
The optimization objective is defined as (14), (16), and (17), which are elaborated in Sect. 4.2. N p is the predictive horizon length.
are the weight matrixes of outputs and controls, respectively. For ease of weight parameters' tuning, the quadratic cost item of u in Eq. (6) is normalized via dividing r 1 and r 2 by 2 F zf 2 and ΔM z 2 , respectively. While for y e and e , the controller is considered to effectively track the target path, and the normalization is not imposed.

Vehicle Stability Constraint of Front Tires' Lateral Force
For vehicle stability, one popular parallelogram envelope is proposed by Stanford University's laboratory [8,34,35]: where Eqs. (7), (8) and (9) are derived from the maximum lateral acceleration and maximum rear wheel sideslip angle, respectively. Equation (9) is to limit f within a stable range since F yf is positively related to f when f ≤ f ≤ f . The envelope of Eqs. (7) to (9) can also be explained in the v yphase plane portraits. The v y -phase plane profiles and the bounds of Eqs. (7) to (9) are shown in the upper subfigures of Fig. 4 (a) to (g). Here, the symbols of solid dark "○" and "◇"are stable and saddle points, respectively; the green solid and dotted lines are the maximum and minimum bounds regarding f and f , respectively; the blue solid and dotted lines are the maximum and minimum bounds regarding r and r , respectively; and the red solid and dotted lines are the maximum and minimum bounds for = ± g∕v x , respectively. With the increasing , the stable points (i.e., the coordinate position v y * , * ) move towards the minus v y direction. When reaches 12 deg or -12 deg, the stable point disappears, and only a saddle point on one side of -axis exists. To determine the stability zones, the unstable portions are removed [8,34,35], as the regions outside blue lines by Eq. (8), red lines by Eq. (7), and green lines by Eq. (9). However, it should note that the maximum bounds in Eq.
Hereinafter, we name the bounds of static = − static = g∕v x (9) F yf ≤ F yf ≤ F yf dynamic bounds are varied with vehicle states and markedly different from static ones. Hence it is difficult to directly introduce the dynamic bounds into a stable region of vehicle by phase planes. Herein, the stability bounds of front tires' lateral force are proposed considering dynamic bounds, and the derivations are: (1) Figure 4 illustrates that with increasing/decreasing , the stable point moves towards the negative/positive saddle point, respectively, and the stable zones disappear after the overlapping of stable point and one saddle point. Under the critical situation to reserve selfstable zones, the bounds derived from f , r , and static _ , interact at the saddle point of positive v y , and those derived from f , r , and static , are crossed over at the saddle point of minus v y . Meanwhile, the lines of dynamic and dynamic are convergent to points of v y * , static and v y * , static , respectively. This is because when v y = v y * , ̇v y = 0 , and the maximum acceleration Therefore, the and for vehicle, stability can be gained by f , r , dynamic , and f , r , dynamic . Considering the angle relations for front and rear axles, there are By introducing ̇v y in Eqs. (10) and (11), and can be adaptively adjusted by the vehicle states to improve the vehicle steering property for PF under extreme conditions. Moreover, Eqs. (10) and (11) are effective to stabilize vehicles since they are calculated for the limited steering of vehicle stability (see Fig. 4 (b), and (f)). If we adopt static and static to calculate and ,̇v y = 0 , and it is a special case of v y = v y * .
(2) To reserve the benefits from F yf control (see Remark 1), Eqs. (10) and (11) are changed to the bounds of f : as static bounds, and the bounds of dynamic = g −̇v y ∕v x and dynamic = − g +̇v y ∕v x dynamic bounds. Although the zero-sideslip angle is generally pursued in vehicle stability control, making ̇v y ≈ 0 , it is not always beneficial for PF issues under extreme conditions. From Eq. (1), sometimes, the reasonably increasing/decreasing v y is good to simultaneously reduce y e and e , but bringing varying ̇v y . This case is similar to that under drift cycles by professional race drivers when a relatively great v y is helpful for extreme turning. More intuitively, the bottom subfigures of Fig. 4a-g illustrate the phase plane trajectories of static and dynamic bounds, where the symbols of hollow "○" represent the convergence points for dynamic bounds. For constant , the  14), the vehicle stability and steering capacity can be guaranteed for extreme drive cases. Its effectiveness and benefits will be demonstrated in the following section of validations.
Remark 3: By the proposed bounds in Eq. (14), the constraint number is reduced from 3 to 1 to reduce the control complexity. Moreover, Eq. (14) is related to control input (i.e., the lateral force of front tires), rather than state constraints, which contributes to two advantages as follows. First, it can reduce computing burdens for some controllers, such as MPC. Second, the proposed bounds are applicable to stabilize vehicles for general state feedback controllers (such as proportional integral and SMC) with ease by adding a saturation process of control inputs. One should note that it is difficult to merge the state constraints into general state feedback controllers with superior performance. For example, the control barrier function is delegated for state constraints handling in Ref. [22], but entails small smoothness in controls.

External Yaw Moment Torques Constraint
The external yaw moment is limited considering the physical constraints of IWMs: . T max is the IWMs' torque limit, and d s is the track width. Remark 4: Benefited by ΔM z , the system has more control DOFs to adjust vehicle motion in PF, which also influences the F yf bounds in Eq. (14). For example, under an extreme drive case, F yf and ΔM z contribute to change vehicle motion for expected PF, while moving away from the equivalent point of vehicle stability. The range of bounds in Eq. (14) will shrink according to the feedback states. Therefore, although the vehicle stability is considered in the control of F yf , F yf and ΔM z can synergistically complete the PF and vehicle stabilization.

Closed-loop Stability Constraint
To make the closed-loop stability of PF control, the Lyapunovbased constraint is added: where u is the auxiliary control law, and V is the corresponding Lyapunov function. The introduction of Eq. (17) allows the NMPC controller to inherit the stability properties of u (0) , so as to analytically characterize a guaranteed region of attraction (ROA).
The process to gain V and u are shown: Define one Lyapunov function as V 1 = 1 2 y e 2 + 1 2 e 2 , and make ̇y e = −k 1 y e , ̇e = −k 2 e , where k 1 > 0 and k 2 > 0 are two specific gain constants. By Eq. (1), if the virtual controls are it yields: However, the real values for v y and are not the same with the virtual controls v yr and r , respectively. Thus, the second Lyapunov function is defined: , (22), and (23), it has The recursive feasibility and closed-loop stability for the established NMPC problem are analyzed: Assumption 1: The desired path is smooth, and the corresponding target heading angle together with its derivatives is smooth and bounded, holding that

Assumption 2:
The command F yf is close to the real value, and their difference is small and can be ignored [8].
Lemma 1: Given that the vehicle is controlled by auxiliary control law u , lateral velocity v y , yaw rate , lateral position error rate ̇y e , and heading angle error rate ̇e , are the NMPC problem admits recursive feasibility.

Constraint Handling Method and Problem Reformulation
The original C/GMRES algorithm cannot directly handle the inequality constraints, thus the relaxation function is applied [36]: and the NMPC problem can be changed to an unconstrained one where j is the jth constraints, j = 1, 2, 3 , corresponding to Eqs. (14), (16), and (17). z denotes the limited variable. z max and z min are the maximum and minimum for z , respectively. is a penalty coefficient, and is a positive variable that determines the varying rate of function value in Eq. (33). In Eq. (33), the shape of function value is influenced by and . That is for different z , z max and z min , the orders of magnitudes for (z) may differ greatly, which will worsen the optimization effects. One method for handling this issue is to specially design two adjustment factors and for each constraint, making similar orders of magnitudes among each (z) , but bringing the huge parameter calibration burdens. To this end, Eq. (33) is modified by normalization: where ẑ = 2z−(z max +z min ) z max −z min , ẑ max = 1,ẑ min = −1 . Figure 5 shows an example of Eq. (35) for the allowable range of z ∈[-4,10] and =10. When −1 < � z < 1 , (z) closes to zero. Instead, if ẑ approaches to the bounds, (z) smoothly enlarges and reaches . Therefore, in the NMPC controller, when the optimized variable closes to bounds, the objective cost will increase, and the algorithm must find the best solution to minimize the objective cost, thus satisfying the constraints.
To discuss the reasonable range, Fig. 6 illustrates the (z) profiles for =10 and different . Generally, it is expected to avoid interferences for optimization until the variables are close to bounds. Therefore, ≥ 10 is considered a good choice.

C/GMRES Algorithm Application
According to Eqs. (5) and (34), the Hamiltonian function is written as H = L + Tẋ , where denotes the co-state vector. Based on the minimum principle [37], the necessary conditions to gain the optimal control sequences are expressed: where H x and H u represent the partial derivative of H regarding x and u , respectively. Here the optimized variables are defined as U = u(0), ..., u(k), ..., u N p − 1 . By the Hamiltonian function and Eq. (36), the optimization problem can be rewritten as The above root-finding problems can be solved by the iterative Newton method, but they suffer from the expensive computation in evaluating Jacobians, Hessian, and inverses [38]. To avoid these, the C/GMRES algorithm is proposed. Equation (37) can be approximated as a dynamic system Ḟ (U, x, t) = A s F(U, x, t) , where A s is a stable matrix to stabilize F at the origin. Then, the forward difference approximation is used, i.e., F UU = A s F − F xẋ − F t , where F U , F x , and F t represent the partial derivative of F regarding U , x , and time t , respectively. Given a small positive value , there is.
Finally, the residuals can be defined as Minimizing the residuals between both sides of Eq. (37) is equivalent to minimizing the approximated one in Eq. (38). To do so, the forward approximation GMRES (FDGMRES) algorithm is applied for fast solving, i.e., where U 0 and U 0 are the initial guesses with respect to U and U , respectively. e tol is the allowable maximum error toleration, and k max is the maximum iterative number. Considering the limited space, the descriptions of FDGMRES are not illustrated here. Interested readers can find more details in Refs. [29,[38][39][40]. For the first sample time in NMPC, the initial guesses U 0 and U 0 are given, and the optimal control sequence increment U * is gained by the FDGMRES algorithm. For the given sample period Δt , the optimal control sequence is calculated by U * = U 0 +U * Δt . Herein, the "warm start-up" mechanism is applied for the expected convergence of optimization [41].

Inversive Tire Model Application and Torque Allocation
With the F yf optimized by NMPC controller, the control command can be determined according to Eq. (4): where (⋅) is the affine function of F yf , , and F zf , and there  Figure 7 yields the illustration of (⋅) for cases of F zf = 8000 N and F zf = 2400 N. Given the where u = T f l , T fr , T rl , T rr T is the torque vector.
is the weight matrix for tire workloads minimization. P 2 = diag p T tot , p ΔM z is the weight vector for satisfying the virtual torques, where where T tot is the total traction torque. The constraints are arranged as, which is imposed to meet the limits of motor torque and tire friction circle, simultaneously [42]. In this paper, a KKT conditions-based optimization method is applied for fast torque allocation, which is simple, globally optimal, and insensitive to initial solutions. Due to the limited space and study emphasis, please find more details in the previous research [43].

Validations and Analysis
The =10 -3 , and k 1 =k 2 =k 3 =k 4 =10 -3 . Moreover, two methods are adopted for comparisons: (1) Linear quadratic regulator (LQR) method: LQR, derived from the optimal control theory, is a state feedback control method. Since nominal LQR cannot handle the inequality constraints, the saturation operations of control constraints Eqs. (9) and (16) are imposed. The state update function of Eq. (5) is applied, and the objective function is set as y Qy + u Ru . Moreover, given that the nominal LQR can only apply to linear models, the cornering stiffness of rear tires is linearized. To adapt the influences of vertical tire forces and tire-road friction factor, the lateral force of rear tires is calculated by F yr = 2C yr , where C yr (⋅) is the mapping relationship about cornering stiffness of rear tires when r = 0 , as shown in Fig. 8. At each sampling moment, the LQR control gains are calculated online by solving Riccati equation according to the feedback states and the calculated F yr . Then, the controls of F yf and ΔM z can be determined.

(2) NMPC by SQP Algorithm (named as "SQP-NMPC"):
The cost function is the same as the proposed strategy, i.e., Eq. (6), and the sequential quadratic programming (SQP) algorithm is used for rolling optimization. Unlike the proposed one, this method adopts the parallelogram envelope and limits of F yf and ΔM z as constraints, namely Eqs. (7) to (9) and (16). Since the SQP algorithm can handle the inequality constraints, the relaxation function is not imposed. For fairness, in LQR and SQP-NMPC controllers, all the parameters are set the same as that of the proposed strategy, and the KKT conditions-based optimization method is employed for torque allocation.

Path Following Effects
The double lane change (DLC) and sinusoid cycles are applied, as shown in Fig. 9, and two extreme drive conditions are given: (1) target v x =100 km/h and =0.8, and (2) target v x =80 km/h and =0.4. Figures 10 and 11 illustrate the PF results under the DLC cycle of target v x =100 km/h and =0.8, and under a sinusoid cycle of target v x =80 km/h and =0.4, respectively. The PF errors by the LQR method are greater than the other two methods, because LQR is essentially a linear controller that cannot optimize the system's nonlinear dynamics in a rolling time horizon. From Figs. 10a and 11a, the proposed strategy can achieve a smaller magnitude of lateral error, compared with SQP-NMPC. Moreover, from the zoomed-in subfigures in Figs. 10 and 11, the proposed strategy yields a much faster and smoother convergence than the SQP-NMPC.
It is noteworthy that theoretically, although the approximations are carried out in the proposed C/GMRES algorithm, the control effects by the two MPC methods should be similar [29], which is inconsistent with the results in Figs. 10 and 11. Figure 12 illustrates the Lyapunov stability constraint profile of the proposed strategy, and all the values are below zero. That means the Lyapunov stability constraint is satisfied, and the ROA for the designed control law u is inherited in the proposed strategy. Given this, superior tracking performance is yielded by the proposed strategy. For example, combining Fig. 10a with Fig. 12a, for the position of global X 50, 80, 110, 140, and 170 m, the stability constraints in the proposed strategy contribute to restricting the Lyapunov function for improved PF convergence, leading to the differences between two MPCs' control effects. Moreover, the relaxation function is verified to be effective in constraint handling in Fig. 12.
More intuitively, Fig. 13 showcases the quantized results of PF errors for three methods. The mean square errors (MSE) regarding lateral errors and heading errors by the proposed strategy are much smaller than that of LQR, making up reductions of 68.22% and 90.31%,  respectively. Roughly, compared with the SQP-NMPC method, the proposed strategy conducts similar MSEs in lateral errors and heading errors. In particular, under DLC cycles, the MSE of lateral errors by the proposed method is significantly smaller than SQP-NMPC's, which reaches more than 30% performance improvements. This is because, with the higher cornering requirements of the DLC cycle than the sinusoids cycle, the restrictions of Lyapunov stability constraint in the proposed strategy are more intense, bringing more desirable tracking convergence to reduce PF errors.

Vehicle Stability Effects
To verify the effectiveness of the proposed stable bounds, a test cycle, that is more extreme than above, is arranged here, i.e., the DLC cycle of target v x =120 km/h and =0.8. Figure 14 illustrates the PF results. The proposed strategy can effectively track the target path with smooth transient performance, while the SQP-NMPC approach exhibits greater errors in lateral position and heading angle. Figures 15 and 16 demonstrate F yf and vehicle state results by the proposed and SQP-NMPC strategies, respectively. For the proposed strategy in Fig. 15, F yf is effectively restricted within the bounds, so , f , r , and a y stay within the allowable range, indicating the vehicle stability. Moreover, since the lateral velocity differential ̇v y is introduced into bounds in the proposed strategy, the greater allowable range of is provided compared with the traditional one (i.e., ± g∕v x =±0.2354 ) during steering, as the results of around 2, 2.7, 4, 4.5, and 5.5 s in Fig. 15 (b). Therefore, the vehicle steering capacity is improved, leading to the expected PF performance in Fig. 14. Additionally, from Fig. 15a, it is effective to apply the proposed relaxation function for constraint handling.
With the constraints Eqs. (7) to (9), the vehicle steady drives by SQP-NMPC strategy, but F yf , , and f reach the boundaries several times. The vehicle steering capacity is limited, resulting in greater PF errors in Fig. 14. For example, at a moment of around 3 s, greater is expected to realize the more satisfactory PF effects. However, to guarantee vehicle stability, F yf and are limited within their boundaries. Meanwhile, there exists no more control degree of freedom to reduce PF errors since ΔM z can only adjust vehicle yaw rate from Eq. (5). Hence, the PF errors cannot be eliminated well, as shown in Fig. 14. This explanation is also applicable for cases of other moments, e.g., at around 3, 5, and 5.7 s, as yielded in Fig. 16a, b. Finally, the PF errors by the SQP-NMPC method are accumulated and rise greater than that of the proposed strategy.  Figure 17 depicts the external yaw moment results. In Fig. 17a, the smaller magnitude of ΔM z is implemented by the proposed method. The SQP-NMPC controller repeatedly adjusts the external yaw moment with greater magnitude for expected PF. To illustrate, at around 5.7 s, F yf reaches the bound, and ΔM z increase to stay within the allowable range, as shown in Fig. 16a, b. Even so, due to its adopted static yaw rate constraints limiting steering capacity, the PF errors are hard to be eliminated, as shown in Fig. 14a. Figure 18 yields the phase plane portrait under a parallelogram envelope. The phase trajectories by the proposed strategy are outside the parallelogram envelope; however, the vehicle is stable according to the results in Fig. 15. This is because the parallelogram envelope is developed on the basis of v y = 0 and only considers the static yaw rate bounds g∕v x . The maximum allowable for vehicle stability is derived from  Figure 19 shows the results of calculation time and iterations by the proposed and SQP-NMPC strategies. Under three drive cycles, the calculation time by the proposed strategy is limited to below 0.01 s and much smaller than the sample time of 0.02 s, testifying the real-time applicability. In contrast, the maximum calculation time by the SQP algorithm is distinctly higher than that by the C/GMRES algorithm. Under the DLC cycle of v x =120 km/h and =0.8, the SQP-NMPC strategy repeatedly adjusts the controls to satisfy the static yaw rate limits (see Fig. 16b), leading to the increased iterations and calculation time to meet the error toleration value of 10 -3 . This is why its calculation time per sample reaches up to around 0.3 s in Fig. 19c.

Conclusions
A real-time NMPC-based strategy is proposed for path following (PF) and vehicle stability under extreme steering conditions. By v y -phase plane analysis, the novel lateral force bounds of front tires are developed and adopted in the NMPC controller, which can define the stable zone of vehicles and will not limit the steering capacity under extreme conditions. Additionally, the Lyapunov-based constraint is developed to guarantee the closed-loop system stability, and its recursive feasibility is proved. The NMPC problem is established with the objective of the minimum PF errors, and the constraints are arranged by the proposed stable bounds of F yf , ΔM z bounds, and the Lyapunov stability constraint. For fast optimization, the C/GMRES algorithm is proposed, and the relaxation function is adopted to handle the constraints. After NMPC optimization, the target lateral force of front tires and external yaw moment is acquired. The target lateral force of front tires is transformed to be the steering angle command by the inverse tire model. Moreover, according to the total traction torque and external yaw moment, the KKT conditions-based allocation method is employed to gain the IWMs' torque commands. The proposed strategy is validated and compared with LQR and SQP-NMPC methods, which implement superior performance and real-time applicability in PF and vehicle stability. Moreover, the developed stable bounds are verified to be effective in guaranteeing stability and improving steering capacity under extreme drive conditions. The conclusions are summed up: (1) By the proposed strategy, the PF performance is greatly superior to that of the LQR controller. Owing to the developed Lyapunov stability constraint, the proposed strategy demonstrates much faster and smoother convergence than the SQP-NMPC in PF. (2) Compared with the parallelogram envelop [8], the proposed stable bounds of lateral forces can realize vehicle stability while improving the steering capacity. Owing to the consideration of lateral velocity differential, the proposed stable bounds can effectively relieve the conflicts between the PF and vehicle sta-bility for extreme drive conditions, benefiting from superior PF effects. (3) The computing time of below 0.01 s per sample is conducted by the proposed C/GMRES algorithm, much smaller than the set sample time of 0.02 s. Instead, by the SQP algorithm, the optimization per sample is completed by up to around 0.3 s, which is hard to be applied in realtime.
The future works are the experimental validations of real-world vehicles and the improvements of the strategy's transient performance as well as robustness.

Appendix 1
With the definition of (t) = y e , e , v y − v yr , − r T in Lemma 1, the Lyapunov function V can be formulated as 1

Appendix 2
It is evident that u is always feasible for NMPC problem if u follows the constraints of F yf and ΔM z . That is u can be seen as a feasible solution in NMPC optimization to analyze the recursive feasibility, as follows. F yf and F yr are coupled, and at limited driving cases, ̇v y = 0 , ̇= 0 , and l a F yf = l b F yr [9]. If the maximum F yr reaches, F yf equals to its maximum, which indicates that the vehicle state locates at one saddle point in phase plane. By Eq. (25) and taking infinity norm on both sides, can be obtained ‖F yf ‖ ∞ = F . Since the constraint of F yf is not symmetrical regarding origin, the bounds of F yf should be expressed as can be obtained Eq. (31). Similarly, taking infinity norm on