Planning Speed Mode of All-Wheel Drive Autonomous Vehicles Considering Complete Constraint Set

Abstract

The study aims to improve the technique of motion planning for all-wheel drive (AWD) autonomous vehicles (AVs) by including torque vectoring (TV) models and extended physical constraints. Four schemes for realizing the TV drive were considered: with braking internal wheels, using a rear-axle sport differential (SD), with braking front internal wheel and rear-axle SD, and with SDs on both axles. The mathematical model combines 2.5D vehicle dynamics model and a simplified drivetrain model with the self-locking central differential. The inverse approach implies optimizing the distribution of kinematic parameters by imposing a set of constraints. The optimization procedure uses the sequential quadratic programming (SQP) technique for the nonlinear constrained minimization. The Gaussian N-point quadrature scheme provides numerical integration. The distribution of control parameters (torque, braking moments, SDs’ friction moment) is performed by evaluating linear and nonlinear algebraic equations inside of optimization. The technique proposed demonstrates an essential difference between forecasts built with a pure kinematic model and those considering the vehicle’s drive/control features. Therefore, this approach contributes to the predictive accuracy and widening model properties by increasing the number of references, including for actuators and mechanisms.

1. Introduction

Usually, planning and tracking stages are separated as hierarchical components; however, a complex approach is also possible [1]. The planning quality reduces the time cost needed for optimizing the control at the tracking stage, as well as minimizes the probability of AV’s unstable motion modes. The variant of decomposing path and speed planning is caused by the need to decrease the number of space-time variations during optimization when the number of constraints increases. Thus, in studies such as [2], first, a visible foresighted trajectory was generated, and the speed plan was then optimized with respect to space. Afterward, the optimal control problem was formulated in the space instead of the time domain.

Models. The problem of smoothed speed arose as a need for high-quality control preventing the appearance of sudden loads and unstable transients. Some studies as [3] attempted to compose speed plans based on cubic polynomials in such a way as to ensure the continuity of accelerations at the transition nodes under conditions of minimizing the jerk. However, the function inflections did not occur directly in the nodes, which led to some excesses of the preset limits followed by the need for complicating the model to compensate for this fluctuation effect. This approach allows obtaining continuous speed values at any point, considering the initial and final values.

The range of models used to represent the AV at the planning stage is quite wide and depends on the control objective. Kinematic models are the simplest and meet the requirements of maximum performance. The bicycle model [4] is traditional and has proven itself well when using model predictive control (MPC). The best quality corresponds to the state-space kinematic models with the control actions such as jerk and the steered wheels’ angles [5].

Dynamic models can be both complex and simple to control a single aspect. For example, in [6] a longitudinal dynamic model is identified to tune a low-level speed-tracking controller. In [7] a simplified constrained model uses a second-order integrator to match the feasible AV dynamic behavior. This also includes speed planning to automatically adjust the vehicle’s capabilities relative to the road conditions. The yaw, roll, and pitch of AV body’s DOFs, each wheel’s dynamics based on slip ratios, and Pacejka’s combined slip tire model are computed. The comparison of kinematic and dynamic models concerning experimental data [4] shows different statistics of forecast errors, including the effect of discretization error accumulation. In [8], a search-based method solves the problem with nonlinear and unstable dynamics.

Objective functions determine the planning and control strategies. Most studies proceed from the control strategy regarding minimizing the motion time while keeping the vehicle within roadway boundaries [9]. Thus, in [10], the time-optimal trajectory around a racetrack is obtained by solving a minimum lap time problem (MLTP) for a racing vehicle. In [11], the path planning strategy includes preventing rear-end collisions during overtaking while minimizing traveling time.

Controllers. Different types of controllers are used for both planning and tracking. The class of controllers based on MPC can be considered the most used that allows obtaining high-quality forecasts due to a combination of hard and soft constraints [5]. Another popular class is the LQR-type controllers. For increasing performance and overcoming physical limitations, a technique can be represented by the augmented Lagrangian framework [6], which refers to iterative LQR (ILQR) and Constrained Iterative LQR (CILQR), respectively. In [12], the sliding mode control (SMC) calculates the total driving force for longitudinal control. In [7], the reference path is followed by low-level tracking using PID controllers. In [13], a learning-based Interaction Point Model (IPM) describes the interaction between agents. In [9], the Nonlinear Model Predictive Control (NMPC) strategy was aimed at controlling a small-scale car model for autonomous racing competitions. Study [6] presented a hierarchical framework with neural physics-driven models to enable the online planning and tracking of minimum-time maneuvers. A lateral speed prediction model for high-level motion planning was considered with economic nonlinear model predictive control (E-NMPC). In [1], a linear parameter-varying (LPV) MPC was deployed for AV trajectory tracking and compared with the linear MPC. The outcomes satisfy the processing rate and high-precision criteria, including safely avoiding obstacles. In [14], the path-following control strategy was based on the linear quadratic regulator (LQR) to compare the performance of models. In [12], the MPC controller was used to calculate the steering wheel angle and the total yaw moment for lateral control, and the sliding mode control (SMC) was used to calculate the total driving force for longitudinal control.

ADAS. Note that the problem of speed planning is typical not only for AVs but also for ADAS systems. In [5], it was proposed to form a speed profile for functioning the adaptive cruise control (ACC) based on MPC. The study [15] considered a preview servo-loop speed control algorithm to achieve smooth, accurate, and computationally inexpensive speed tracking for connected automated vehicles (CAVs). A classic PID with an optimal controller was implemented in an automated vehicle platform for lowering speed-tracking errors, mitigating operations, and smoothing brake/throttle activations. The paper [11] proposed a system for speed planning using MPC to estimate the vehicle overtake safety while controlling the maneuver for a dynamic vehicle model.

Paper [16] presented an AV’s trajectory planning and speed control algorithm for static traffic agent avoidance in multi-vehicle urban environments. Decision-making and motion planning frameworks were presented considering both the preceding static and surrounding vehicles. Paper [17] considers hierarchical architecture for decision making, motion planning, and control for an autonomous racing vehicle. The supervisor determined whether the subject vehicle should stay behind the preceding vehicle or start overtaking. A high-level trajectory planner generated the desired path and velocity profile in a receding horizon mode. A low-fidelity kinematic vehicle model was utilized for planning, and a low-level trajectory tracker based on MPC and dynamic vehicle model generated the lateral and longitudinal control inputs.

Some studies implemented the use of four-wheel models with increased degrees of freedom, refined tire models, actuators, etc. Thus, in [14], the question of tire dynamics’ influence on vehicle dynamics’ control was raised. The influence of tire nonlinearity on path offset was revealed. In [12], the longitudinal and lateral motion control for the four-wheel independent drive intelligent vehicle was designed: the upper layer was the motion controller, and the lower layer was the control distributor. In [8], complex dynamics in a predictive manner were applied to achieve optimal robot behavior in dynamic scenarios. In paper [18], a new real-time motion planning technique for an autonomous mobile robot is proposed considering actuator capacities’ limits and tire-ground adhesion constraints.

Actuators. Study [9] includes a simple drivetrain model in the optimization problem to limit the lateral and longitudinal forces acting on a car. An Autonomous Emergency Braking Pedestrian (AEB-P) was introduced in [19] to prevent collisions between vehicles and pedestrians. The emergency brake planner generates vehicle deceleration followed by tracking the trajectory, where the PI-controller was adapted to provide the optimum braking force. In [15], the brake/throttle control laws were introduced in five parts: three feedback controls of system states and two feedforward items previewing road slope and target speed. Study [20] proposes a clothoid-curve-based trajectory tracking control method to solve the problem of tracking errors caused by the discontinuous curvature of the control curve calculated by simple algorithms. The parameters of clothoid reference curves are tied to the vehicle’s safe motion constraints.

Despite the intensive studies on motion planning, in our opinion, there still are some issues that can be conditionally divided into three parts: vehicle model, restrictions, and planning techniques.

Vehicle model. Note that most studies in the field of AV motion planning are focused on generating reference curves using relatively simple (often kinematic bicycle) models. However, the motion planner’s high performance does not yet mean high forecast quality. The most important indicators are the prognosis feasibility and ensuring good accuracy in reflecting features of a specific vehicle design. Considering the intensive use of control systems in modern vehicles (for example, traction distribution), it is often not enough to obtain the trajectory and kinematic parameters’ references. Moreover, the desirable parameters themselves may significantly depend on the technologies used in a real vehicle. For example, the distribution of drive torques between the same axle’s wheels can be carried out through a symmetrical differential, through a sport differential (torque vectoring), and by activating the inner wheels’ brakes. All three options will cause different axle drive power, different traction and lateral forces, the tire–road adhesion degree, and, as a result, affect differently on forming the yaw moment. In a classic uncontrolled AWD drivetrain, a wheel with the worst adhesion conditions limits the traction potential. And vice versa, in transmissions design with individual wheel control, the maximum realization of total adhesion potential is possible. In the case of using a standard vehicle model, there will be no difference between these two variants when speed is the planning object. If an extended model is used, including the drive type and redistribution of vertical reactions, then, as shown in our previous articles, the difference between predictions of kinematic parameters may be significant especially for cases with low tire adhesion. Thus, this study’s target consists of composing an inverse vehicle dynamic model that allows estimating physical factors based on the kinematic parameters.

A set of constraints is a critical factor in reflecting the realism of forecasting and the optimization procedure’s performance. On the one hand, increasing the number of restrictions complicates the mathematical optimization model; on the other hand, narrowing the boundaries of AV motion parameters contributes to diminishing solution iterations and time costs.

When using a conventional kinematic model, constraints are also basically kinematic (linear and angular velocities, linear and angular accelerations, jerk, etc.), which does not always correspond to natural physical limits. For example, many researchers use a constant upper acceleration limit without relating it to the throttle response properties (speed–acceleration) for a specific vehicle design. At the same time, the traction potential decreases with speeding even for the case of electric drive, which especially should be considered for a medium-speed range.

The physical limitations are the most important but require an extended vehicle model. First, these restrictions are responsible for preventing motion instability, critical modes of vehicle-road interaction, and are associated with the design solutions for distributing the power by wheels. Controlling the tires’ adhesion limits is the most important aspect, requiring predicting both a redistribution of normal reactions and torques on wheels. In turn, the torque depends on the drive design and control. In addition, pure kinematic models do not allow direct operation with longitudinal and lateral slip. Nevertheless, it is possible to indirectly estimate tires’ slip angles and lateral deformations followed by including them in the constraints.

Safety restrictions. Most often, this includes speed limits when approaching moving and static obstacles. Within this planning technique, we suppose that a trajectory considering obstacles and speed mode has already been built. Therefore, in the framework of this study, we assume that the AV uses a safe space for realizing a maneuver. As a safety measure, the body roll angle caused by lateral inertia forces is applied. Restricting the maximum roll angle, limitations on lateral acceleration are automatically imposed.

Thus, in advance, we proceed from the fact that growth in the number of restrictions positively affects the rapidity and quality of forecasting the AV motion.

Planning technique. Here, the most important requirements concern ensuring the smoothness and unambiguousness of predicted functions. The most common approach implies using polynomials that describe the change in speed concerning time along a trajectory section. However, this does not always lead to a qualitative consistency between the speed and curvature derivatives involved in the formation of kinematic parameters such as jerk and angular acceleration. That is, there may be disruptions in the piecewise representation of these functions. This, in turn, does not contribute to the possibility of using inverse dynamics to determine the force factors acting on the vehicle. In this regard, our method consists of using an inverse approach, when instead of speed its second derivative is directly modeled by smooth functions. In this case, the continuity and smoothness of the vehicle’s linear and angular accelerations will be guaranteed.

Figure 1 shows a general approach scheme reflecting the modeling stages and the parameters used. The SQP optimization ensures iterative process, within which kinematic, dynamic, and physical parameters are evaluated to form optimization criteria and constraints. To model vehicle kinematics, the inverse method is used, allowing the representation of the longitudinal speed function via finite elements. The external data are the preset boundary values of the kinematic parameters, initial conditions, and the precalculated trajectory model. The SQP-loop at each iteration generates a vector of nodal parameters for the longitudinal speed’s second derivative. A simple and smooth form of the speed function is obtained by double numerical integration. After evaluating accelerations, the distribution of normal reactions over the wheels can be found. Using the AWD transmission model, vehicle steerability and nonlinear tire models, and algorithms for controlling traction distributing actuators, the inverse dynamics problem is solved, followed by determining the needed longitudinal and lateral forces and estimating degrees of longitudinal and total tire adhesion use. Knowing the lateral forces and accelerations allows estimating the equivalent sideslip angles, tire lateral deformations, and the body roll angle, including them in a set of nonlinear constraints. As a result, the optimization output contains kinematic, dynamic, and physical parameters, as well as the necessary control for actuators.

2. Existing 2.5D Vehicle Dynamics Model

2.1. Model Scheme

Let us present the schemes to reflect the active forces, geometric, and kinematic parameters needed for parametrizing the AV model (Figure 2).

2.2. Trajectory-Based Vehicle Kinematics

Assuming the maneuver trajectory is preoptimized, we outline the main points of the ideal vehicle kinematics providing better controllability and motion stability. At the same time, the main task consists of linking the kinematic parameters with the trajectory’s geometric ones, ensuring the continuity of their derivatives (Figure 2a).

Mass center speed V is directed by the trajectory tangent and may be expressed as follows

$$\overrightarrow{V}=V_{\zeta }{\overrightarrow{u}}_{\zeta }+V_{\mu }{\overrightarrow{u}}_{\mu }=V_x{\overrightarrow{u}}_x+V_y{\overrightarrow{u}}_y$$

(1)

where ${\overrightarrow{u}}_{\zeta }$, ${\overrightarrow{u}}_{\mu }$, V_ζ, V_μ = unit vectors and speed components of the vehicle local coordinate system ζμ; ${\overrightarrow{u}}_x$, ${\overrightarrow{u}}_y$, V_x, V_y = unit vectors and speed components of the fixed (global) coordinate system xy.

The absolute speed V may be decomposed with projections V_x and V_ζ tied through the tangent angle α and central slip angle β. Then,

$$V=\frac{ds}{dt}=\frac{ds}{dx}\frac{dx}{dt}=\frac{V_x}{cos\left(\alpha \right)}=\frac{V_{\zeta }}{\mathit{cos}\left(\beta \right)},\hspace{1em}{V}_x=V_{\zeta }\frac{\mathit{cos}\left(\alpha \right)}{\mathit{cos}\left(\beta \right)}$$

(2)

The first and second derivatives of V_ζ may be found as follows

$$\frac{d{V}_{\zeta }}{dt}=\frac{d{V}_{\zeta }}{dx}\frac{dx}{dt}=\frac{d{V}_{\zeta }}{dx}{V}_x,\hspace{1em}\frac{d^2{V}_{\zeta }}{d{t}^2}=\frac{d^2{V}_{\zeta }}{d{x}^2}{V}_x^2+\frac{d{V}_{\zeta }}{dx}\frac{d{V}_x}{dx}{V}_x$$

(3)

The derivative of V_x concerning the x-coordinate is

$$\frac{d{V}_x}{dx}=\left(\frac{d{V}_{\zeta }}{dx}+V_{\zeta }\left(\frac{d\beta }{\mathit{dx}}\mathit{tan}\left(\beta \right)-\frac{d\alpha }{\mathit{dx}}\mathit{tan}\left(\alpha \right)\right)\right)\frac{\mathit{cos}\left(\alpha \right)}{\mathit{cos}\left(\beta \right)}$$

(4)

Lateral Speed V_μ is geometrically tied with the longitudinal speed V_ζ by the central slip angle β

$${V_{\mu }=V}_{\zeta }\mathit{tan}\left(\beta \right)$$

(5)

Its time derivative

$$\frac{d{V}_{\mu }}{dt}=\frac{d{V}_{\mu }}{dx}\frac{dx}{dt}=\frac{d{V}_{\mu }}{dx}{V}_x,\hspace{1em}\frac{d{V}_{\mu }}{dx}=\frac{d{V}_{\zeta }}{dx}\mathit{tan}\left(\beta \right)+V_{\zeta }\frac{d\beta }{\mathit{dx}}{\mathit{sec}}^2\left(\beta \right)$$

(6)

Yaw Rate is the yaw angle ϕ = α – β derivative in the current global coordinates. Thus,

$$\omega =\frac{d\varphi }{dt}=\frac{d\varphi }{dx}\frac{dx}{dt}=\frac{d\varphi }{dx}{V}_x,\hspace{1em}\frac{d\omega }{dx}=\frac{d}{dx}\left(\frac{d\varphi }{dx}{V}_x\right)=\frac{d^2\varphi }{d{x}^2}{V}_x+\frac{d{V}_x}{dx}\frac{d\varphi }{dx}$$

(7)

Angular acceleration ε is derived from the yaw rate ω concerning time

$$\epsilon =\frac{d\omega }{dt}=\frac{d\omega }{dx}\frac{dx}{dt}=\frac{d\omega }{dx}{V}_x=\frac{d^2\varphi }{d{x}^2}{V}_x^2+\frac{d{V}_x}{dx}\omega$$

(8)

Accelerations in the local vehicle coordinate system ζμ are

$$\overrightarrow{a}=\frac{d\overrightarrow{V}}{dt}={\left(\begin{array}{c}{a}_{\zeta }\\ a_{\mu }\end{array}\right)}^T\left(\begin{array}{c}{\overrightarrow{u}}_{\zeta }\\ {\overrightarrow{u}}_{\mu }\end{array}\right),\hspace{1em}\left(\begin{array}{c}{a}_{\zeta }\\ a_{\mu }\end{array}\right)=\left(\begin{array}{c}d{V}_{\zeta }/dt\\ d{V}_{\mu }/dt\end{array}\right)+\omega \left(\begin{array}{c}-V_{\mu }\\ V_{\zeta }\end{array}\right)$$

(9)

Jerks in the vehicle coordinate system are

$$\overrightarrow{j}=\frac{d\overrightarrow{a}}{dt}={\left(\begin{array}{c}{j}_{\zeta }\\ j_{\mu }\end{array}\right)}^T\left(\begin{array}{c}{\overrightarrow{u}}_{\zeta }\\ {\overrightarrow{u}}_{\mu }\end{array}\right)={\left(\begin{array}{c}\frac{d^2{V}_{\zeta }}{d{t}^2}-\left(2\frac{d{V}_{\mu }}{dt}+V_{\zeta }\omega \right)\omega -\epsilon V_{\mu }\\ \frac{d^2{V}_{\mu }}{d{t}^2}+\left(2\frac{d{V}_{\zeta }}{dt}-V_{\mu }\omega \right)\omega +\epsilon V_{\zeta }\end{array}\right)}^T\left(\begin{array}{c}{\overrightarrow{u}}_{\zeta }\\ {\overrightarrow{u}}_{\mu }\end{array}\right)$$

(10)

Velocities and accelerations at the wheels’ centers are needed for assessing the forces acting in the tire–road contacts. Let us introduce designations

$$\begin{gathered} U=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right),\hspace{1em}UU=-\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)=-E_2,\hspace{1em}{r}_j=\left(\begin{array}{c}{r}_{j\zeta }\\ r_{j\mu }\end{array}\right),\hspace{1em}v=\left(\begin{array}{c}{V}_{\zeta }\\ V_{\mu }\end{array}\right),\hspace{1em}{v}_j=\left(\begin{array}{c}{v}_{\zeta j}\\ v_{\mu j}\end{array}\right), \\ H\left(·\right)=\left(\begin{array}{cc}\mathit{cos}\left(·\right)& \mathit{sin}\left(·\right)\\ -\mathit{sin}\left(·\right)& \mathit{cos}\left(·\right)\end{array}\right),\hspace{1em}\left(\begin{array}{c}{\overrightarrow{u}}_{\zeta j}\\ {\overrightarrow{u}}_{\mu j}\end{array}\right)=H\left({\theta }_j\right)\left(\begin{array}{c}{\overrightarrow{u}}_{\zeta }\\ {\overrightarrow{u}}_{\mu }\end{array}\right),\hspace{1em}\left(\begin{array}{c}{\overrightarrow{u}}_{\zeta }\\ {\overrightarrow{u}}_{\mu }\end{array}\right)=H^T\left({\theta }_j\right)\left(\begin{array}{c}{\overrightarrow{u}}_{\zeta j}\\ {\overrightarrow{u}}_{\mu j}\end{array}\right) \end{gathered}$$

(11)

where E₂ = identity matrix of dimension 2 × 2, r_j = vector of j-th wheel center’s coordinates in the AV local coordinate system, $v_j$ = vector of j-th wheel center’s velocities in the wheel local coordinate system, H = rotational matrix, θ_j = angle of j-th wheel turn, and ${\overrightarrow{u}}_{\zeta j}$, ${\overrightarrow{u}}_{\mu j}$ = unit vectors of the j-th wheel local coordinate system.

Speed and Acceleration vectors in the center of j-th wheel’s local coordinates

$$v_j=\left(\begin{array}{c}{v}_{\zeta j}\\ v_{\mu j}\end{array}\right)=H\left({\theta }_j\right)\left(v+\omega U^T{r}_j\right),\hspace{1em}{a}_j=\left(\begin{array}{c}{a}_{\zeta j}\\ a_{\mu j}\end{array}\right)=H\left({\theta }_j\right)\left(a+\left(\epsilon U^T-{\omega }^2{E}_2\right)r_j\right)$$

(12)

Rotation angles of steered wheels can be determined as Ackerman’s angles subject to that the trajectory curvature K is much less than the maximum possible. Thus,

$${\theta }_1=arccot\left(\frac{\mathit{cos}\left(\beta \right)}{K{B}_{12}}-\frac{B_{13}}{2B_{12}}\right),\hspace{1em}{\theta }_3=arccot\left(\frac{\mathit{cos}\left(\beta \right)}{K{B}_{12}}+\frac{B_{13}}{2B_{12}}\right)$$

(13)

2.3. Vehicle Dynamics

2.3.1. Body Pseudo-Roll

To enhance the influence of transversal forces on vehicle safety and redistributing vertical reactions, the body roll phenomenon may also be used. It can be estimated using the simplest scheme, as depicted in Figure 2b. Assuming that the vehicle body mass center can deviate on the angle ψ relative to a roll center, the dynamics equation for the balance of moments is derived as follows

$${\overrightarrow{M}}_{\psi }={\overrightarrow{r}}_{\psi }\times \left({\overrightarrow{P}}_{\mu }+m_b\overrightarrow{g}\right)={\overrightarrow{r}}_{\psi }\times \left({\overrightarrow{a}}_{\mu }+\overrightarrow{g}\right)m_b$$

(14)

where m_b = mass of the vehicle body, r_ψ = roll lever, g = gravity acceleration constant.

Performing a cross product, obtain

$$M_{\psi }{\overrightarrow{u}}_{\psi }=r_{\psi }{\overrightarrow{u}}_r\times \left({\overrightarrow{u}}_{\mu }{a}_{\mu }-{\overrightarrow{u}}_{z}g\right)m_b=r_{\psi }\left(a_{\mu }\mathit{cos}\left(\psi \right)+g\frac{B_{12}}{B_r}\mathit{sin}\left(\psi \right)\right)m_b{\overrightarrow{u}}_{\psi }$$

(15)

where ${\overrightarrow{u}}_{\psi }$, ${\overrightarrow{u}}_r$ = unitary vectors along the roll axis and roll lever, respectively,

$$B_r=\sqrt{B_{12}^2+{\left({h_{24}-h}_{13}\right)}^2}$$

At the same time, the moment modulus is proportional to the angle ψ and the angular stiffness coefficient k_ψ

$$M_{\psi }=k_{\psi }\psi$$

(16)

Equating the modules of Equations (15) and (16), we obtain for the angle ψ

$$\psi =\frac{{m_{b}r}_{\psi }}{k_{\psi }}\left(a_{\mu }\mathit{cos}\left(\psi \right)+g\frac{B_{12}}{B_r}\mathit{sin}\left(\psi \right)\right)$$

(17)

Supposing the angle ψ to be small, the expression of Equation (17) may be reduced to the linear form. If the constant part is denoted as p_ψ, then the roll angle will be a linear function of the lateral acceleration

$$\psi =p_{\psi }{a}_{\mu },\hspace{1em}\hspace{1em}{p}_{\psi }=1/\left(\frac{k_{\psi }}{{m_{b}r}_{\psi }}-g\frac{B_{12}}{B_r}\right)$$

(18)

The roll lever r_ψ (Figure 2c) can be estimated based on dependencies [21]

$$r_{\psi }=h_g\frac{B_{12}}{B_r}-\left(h_{24}c+h_{13}b\right)/B_r\approx h_g-\left(b{h}_{13}+c{h}_{24}\right)/B_{12}$$

(19)

where h₁₃, h₂₄ = ordinates of the roll centers of the front and rear suspensions, respectively (double wishbone suspensions are usually slightly below the wheel-center axis).

The angular stiffness of suspensions can be defined as follows

$$k_{\psi }=k_{s13}+k_{s24},\hspace{1em}{k}_{s13}=\frac{1}{2}{B}_{13}^2{k}_{sf},\hspace{1em}{k}_{s24}=\frac{1}{2}{B}_{24}^2{k}_{sr}$$

(20)

k_sf, k_sr = reduced stiffness of independent front and rear suspensions, respectively.

The suspension stiffness can be estimated from the lowest natural frequency of free oscillations. Then,

$$k_{sf,r}={4{\pi }^{2}m}_{sf,r}{f}_{sf,r}^2$$

(21)

m_sf,r = sprung masses distributed on every single front and rear suspensions, respectively, f_sf,r = natural oscillation frequencies of sprung masses (for modern cars f_s = 0.6…1.3 Hz).

2.3.2. Redistribution of Vertical Reactions

The distribution of vertical reactions is the basic factor determining the wheel’s potential in transmitting longitudinal and transversal forces. The nature of the longitudinal and transverse accelerations generated during optimization should be such as to provide an approximate constancy of adhesive properties without violating the imposed restrictions. Considering (Figure 2a), the balance of moments in the vehicle’s longitudinal plane can be derived as follows

$$\left({\mathit{l}}_x\mathit{cos}\left(\gamma \right)+h_g\mathit{e}\left(\mathit{sin}\left(\gamma \right)+\frac{a_{\zeta }}{g}\right)\right)mg+\mathit{e}{P}_{ad}{h}_a+\left(\begin{array}{cc}\mathit{e}& \mathit{e}\end{array}\right){\mathit{M}}_{\mu }-{\mathit{R}}_{za}{B}_{12}={\mathbf{0}}_2$$

(22)

where γ = angle of road slope, and

$${\mathit{R}}_{za}=\left(\begin{array}{c}{R}_{zf}\\ R_{zr}\end{array}\right),\hspace{1em}{\mathit{M}}_{\mu }=\left(\begin{array}{c}{M}_{\mu f}\\ M_{\mu r}\end{array}\right),\hspace{1em}{\mathit{l}}_x=\left(\begin{array}{c}b\\ c\end{array}\right),\hspace{1em}\mathit{e}=\left(\begin{array}{c}-1\\ 1\end{array}\right),\hspace{1em}{B}_{12}=b+c$$

(23)

Provided that for most of the moving time γ ≈ 0, Equation (22) may be reduced to

$$\left({\mathit{l}}_x+h_g\mathit{e}\frac{a_{\zeta }}{g}\right)mg+\mathit{e}{P}_{ad}{h}_a+\left(\begin{array}{cc}\mathit{e}& \mathit{e}\end{array}\right){\mathit{M}}_{\mu }-{\mathit{R}}_{za}{B}_{12}={\mathbf{0}}_2$$

(24)

Wheel rolling resistance moments M_r can be partially estimated by omitting kinematic losses (no sliding) and seeing only power losses. However, this requires a normal load value on each wheel. To simplify, we may consider the moments from conditionally paired front and rear wheels. Then,

$$M_{rf\left(r\right)}=R_{zf\left(r\right)}{r}_{wf\left(r\right)}{f}_{rf\left(r\right)}$$

(25)

where r_wf(r) = dynamic radius of front (rear) wheel, f_rf(r) = coefficient of front (rear) rolling resistance.

Using the steered wheels’ average turning angle θ_f(r), we obtain

$$M_{\mu f\left(r\right)}=M_{rf\left(r\right)}\mathit{cos}\left({\theta }_{f\left(r\right)}\right)$$

(26)

Then, the rolling resistance moments can be derived in terms of vertical reactions

$${\mathit{M}}_{\mu }=\left(\begin{array}{c}{R}_{zf}{r}_{wf}{f}_{rf}\mathit{cos}\left({\theta }_f\right)\\ R_{zr}{r}_{wr}{f}_{rr}\mathit{cos}\left({\theta }_r\right)\end{array}\right)=\left(\begin{array}{cc}{r}_{wf}{f}_{rf}\mathit{cos}\left({\theta }_f\right)& 0\\ 0& r_{wr}{f}_{rr}\mathit{cos}\left({\theta }_r\right)\end{array}\right)\left(\begin{array}{c}{R}_{zf}\\ R_{zr}\end{array}\right)={\mathit{m}}_{\mu }{\mathit{R}}_{za}$$

(27)

Thus, Equation (23) may be rewritten in the form

$${\mathit{R}}_{za}={\left(B_{12}{E}_2-\left(\begin{array}{cc}\mathit{e}& \mathit{e}\end{array}\right){\mathit{m}}_{\mu }\right)}^{-1}\left(\left({\mathit{l}}_x+h_g\mathit{e}\frac{a_{\zeta }}{g}\right)mg+\mathit{e}{P}_{ad}{h}_a\right)$$

(28)

Assuming that the increments of normal reactions on the right and left sides are equal due to the vehicle symmetry relative to the longitudinal plane, we obtain the vector of dynamic reactions caused by the longitudinal acceleration

$${\mathit{R}}_{z\zeta }=\left(\begin{array}{c}{\mathsf{\Delta }R}_{z\zeta 1}\\ {\mathsf{\Delta }R}_{z\zeta 2}\\ \begin{array}{c}{\mathsf{\Delta }R}_{z\zeta 3}\\ {\mathsf{\Delta }R}_{z\zeta 4}\end{array}\end{array}\right)=\frac{1}{2}\left(\begin{array}{cc}1& 0\\ 0& 1\\ 1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}{R}_{zf}\\ R_{zr}\end{array}\right),\hspace{1em}{\mathit{R}}_{z\zeta }=\frac{1}{2}\left(\begin{array}{c}{E}_2\\ E_2\end{array}\right){\mathit{R}}_{za}$$

(29)

The distribution of vertical reactions along the sides, considering Figure 2d, may be derived from the equations of moments’ equilibrium relative to each wheel’s contact center. Designate

$${\mu }_w=\frac{1}{2}\left(\begin{array}{c}{B}_{13}\\ B_{24}\\ \begin{array}{c}{-B}_{13}\\ {-B}_{24}\end{array}\end{array}\right),e=\left(\begin{array}{c}1\\ 1\\ \begin{array}{c}1\\ 1\end{array}\end{array}\right),{\mathit{R}}_{z\mu }=\left(\begin{array}{c}{\mathsf{\Delta }R}_{z\mu 1}\\ {\mathsf{\Delta }R}_{z\mu 2}\\ \begin{array}{c}{\mathsf{\Delta }R}_{z\mu 3}\\ {\mathsf{\Delta }R}_{z\mu 4}\end{array}\end{array}\right),\mathit{\mu }=e{\mu }_w^T,\mathit{Y}=\left(\mathit{\mu }-{\mathit{\mu }}^T\right),{\mathit{M}}_{\zeta }=e\left(M_{\zeta r}+M_{\zeta l}\right)$$

(30)

where M_ζr, M_ζl = overturning moment from the right and left vehicle sides.

Since the tire deformations in the lateral direction are omitted at this stage, only the rolling resistance of moment components is taken as the wheel overturning moment. Considering them along the axes, we obtain

$$M_{\zeta f\left(r\right)}=M_{rf\left(r\right)}\mathit{sin}\left({\theta }_{f\left(r\right)}\right),\hspace{1em}{\mathit{M}}_{\zeta }=e\left(M_{rf}\mathit{sin}\left({\theta }_f\right)+M_{rr}\mathit{sin}\left({\theta }_r\right)\right)=e{M}_{rf}\mathit{sin}\left({\theta }_f\right)$$

(31)

Since the model assumes a pseudo-roll, which causes some displacement of the mass center (Figure 2b), the gravity effect can be clarified with the moment m_b·g·r_ψ·sin(ψ) ≈ m_b·g·r_ψ·ψ, where the pseudo-roll angle can be preliminarily estimated from the lateral acceleration. The matrix equation of the moment balance relative to the centers of the wheels’ contacts is represented as follows

$$\mathit{Y}{\mathit{R}}_{z\mu }+{ma}_{\mu }{h}_g+{\mathit{M}}_{\zeta }+m_{b}g{r}_{\psi }\psi =0$$

(32)

The matrix Y is singular because of the vehicle symmetry. Therefore, the pseudo inverse matrix Y⁺ is used for redistributing vertical reactions caused by the lateral acceleration. Then,

$${\mathit{R}}_{z\mu }=-{\mathit{Y}}^{+}\left({ma}_{\mu }{h}_g+m_{b}g{r}_{\psi }\psi +{\mathit{M}}_{\zeta }\right)$$

(33)

Thus, the vector of vertical reactions is a combination of both distributions

$${\mathit{R}}_z={\mathit{R}}_{z\zeta }+{\mathit{R}}_{z\mu }$$

(34)

2.3.3. Pseudo-Tires’ Slip Angles

Since a kinematic model is used for speed planning, which does not imply the tire’s side elasticity and sideslip, the motion trajectory may acquire a certain excessively ideal character. For partial leveling of this drawback, it is possible to estimate in the first approximation the values of tires’ lateral sideslips and deformations stipulated by the speed V_ζ and acceleration a_μ. Since the lateral accelerations and speeds at the wheels’ centers are known by Equation (12), it is possible to evaluate the lateral forces required for stable motion. To do this, briefly consider the main factors affecting the lateral sideslip.

When the slip angel δ changes from zero to some value b (segment 0b, Figure 3a), which is different for various tires, normal loads, and tire–road friction coefficients, the dependence F_μ = f(δ) is almost linear,

$$F_{\mu }=k_{\delta }\delta$$

(35)

This segment corresponds to F_μ, at which the slippage is small.

The segment bc corresponds to the F_μ values, at which the slip occurs on a significant part of the tire’s contact area, and the more intensively the larger δ is. At point c, the force F_μ reaches its maximum value according to the adhesion condition, and in segment cd, it is determined by the equality F_μmax = R_zφ_μmax, where φ_μmax is the coefficient of transverse adhesion. Conditionally, in segment 0c, the wheel lateral movement caused by the force F_μ is called the sideslip, and in segment cd it is called the pure slip. The angle value δ_c, at which sliding begins, depends on the tire structure, normal load, coefficient φ_μmax, and other factors. Usually, for a dry road surface δ_c = 12–20°. In terms of kinematics, it does not matter for what reason the velocity vector deviation from the rotation plane occurs; therefore, the angle δ is formed by the components of the wheel’s longitudinal and transverse speeds over the entire range 0d.

The vehicle steerability is largely defined by the dependence k_δ = f(R_z) (Figure 3b). For passenger vehicles, k_δ has a maximum value k_δmax when the force R_zopt is close to that corresponding to the vehicle’s gross weight R_z0. The coefficient k_δ value depends on the coefficient φ_μmax. In segment 0b, the coefficient k_δ is practically independent of the coefficient φ_μmax. However, the smaller φ_μmax, the smaller the angle δ value corresponding to the point d. On a dry road surface with a force R_z corresponding to standards for a given tire, it can be taken equal to δ_b ≈ 3–4°. Formula (35) is also valid for the section bc, at which k_δ = f(δ) or k_δ = f(F_μ) (Figure 3c). In this segment, the smaller φ_μmax, the smaller the values of k_δ.

During the vehicle motion, the change in value and direction of the wheel’s longitudinal force R_ζ affects the dependence F_μ = f(δ). With an increase in R_ζ/R_z, to obtain the same angles δ in the traction mode, a smaller force F_μ is needed. In the braking mode, at small R_ζ/R_z, its increase leads to a slight growth in the force F_μ, and at larger R_ζ/R_z—the F_μ decreases. If R_ζ/R_z is small, then the effect of longitudinal forces F_μ = f(δ) is insignificant. The longitudinal forces exert the greatest influence at values of the R_ζ forces that are close in value to the maximum adhesion.

Slippage starts when R_max = R_zφ_max. In this case, the adhesion coefficient φ_max depends on the slip direction of the contact patch relative to the road surface, which coincides with the direction of the reaction R. When sliding in the rotation plane, the adhesion coefficient is φ_ζ, and in the transverse direction φ_μ.

With the onset of lateral slip, slippage in the wheel’s rotation plane also occurs. The ratio between the speeds of the lateral and longitudinal slips is the same as between the lateral and longitudinal reactions.

Antonov D. A. [21] proposed considering the influence of various factors on the coefficient k_δ by multiplying k_δmax with several correction factors. With the wheel’s straight-line steady rolling on an even, nondeformable road surface

$$k_{\delta }=k_{\delta max}{q}_{\phi }{q}_z{q}_{\zeta }{q}_p$$

(36)

where k_δmax = the driven wheel’s slip resistance coefficient on the dependence F_μ = f(δ) linear section at the maximum values of the dependences k_δ = f(R_z) or k_δ = f(p_a). q_φ = correction factor considering the slip resistance dependence on the angle δ while moving on roads with different φ (non-linear dependence F_μ = f(δ)), q_z = correction factor considering the effect of normal load deviation from the optimum one, q_ζ = correction factor considering the influence of longitudinal forces on k_δ, and q_p = correction factor considering the effect of the tire’s air pressure p_a deviation from the optimal pressure (supposedly neglected and equals 1).

The coefficients can be calculated from the following dependences

$$q_{\phi }=\frac{arctan\left({\lambda }_{\phi }\left(\delta -{\delta }_0\right)\right)}{{\lambda }_{\phi }\left(\delta -{\delta }_0\right)},\hspace{1em}{\lambda }_{\phi }=\frac{\pi k_{\delta 0}}{2{\phi }_{\mu }{R}_z}$$

(37)

where k_δ0 = the slip resistance factor at a given value of R_z in a linear section F_μ = f(δ), and δ₀ = the slip angle corresponding to the transition from a linear section to a nonlinear section (regarding b). Usually, δ₀ = 0.025–0.035 rad.

$$q_z=0.4{\lambda }_z^3-1.8{\lambda }_z^2+2.4{\lambda }_z,\hspace{1em}{\lambda }_z=\frac{R_z}{R_{zopt}}\approx \frac{R_z}{R_{z0}}$$

(38)

$$q_{\zeta }=\frac{\sqrt{1-{\left(R_{\zeta }/R_{\zeta max}\right)}^2}}{1+0.375{\phi }_{\zeta }},\hspace{1em}{\phi }_{\zeta }=\frac{R_{\zeta }}{R_z},\hspace{1em}{R}_{\zeta max}=R_z{\phi }_{\zeta max}$$

(39)

The dependence F_μ = f(δ) corresponds to the steady-state wheel rolling when F_μ = const and the trajectory of the wheel’s center is a straight line. In most cases, this dependence can also be used to study controlled motion at F_μ = f(t).

As mentioned, the linear segment of the tire’s lateral force characteristic continues to the slip angle of about δ_b ≈ 3–4° under a static load and dry surface at the wheel’s driven mode, providing the maximum potential of φ_μ. Since a decrease in φ_μ leads to a reduction in k_δ0, the angle δ_b also drops. Assuming a linear dependence, the initial angle δ₀ of the linear zone can be calculated as

$${\delta }_0\approx 0.4R_z\frac{{\phi }_{\mu }}{k_{\delta 0}}\le {\delta }_b$$

(40)

Thus, the lateral force, considering the load, adhesion, and traction use, can be expressed as a dependence on the slip angle as

$$F_{\mu }=f\left(\delta \right)=-\left\{\begin{array}{cc}{k}_{\delta 0}\delta ,\hfill & \hfill 0\le \left|\delta \right|\le {\delta }_0\\ k_{\delta 0}{q}_{\phi }\delta ,\hfill & \hfill \left|\delta \right|>{\delta }_0\end{array}\right.,\hspace{1em}{k}_{\mu 0}=k_{\delta max}{q}_z{q}_{\zeta }{q}_p$$

(41)

This nonlinear characteristic will be used for determining the wheels’ side reactions.

2.3.4. Pseudo-Tires’ Lateral Deformations

The tire lateral deformation also directly affects the wheel slip angle and the trajectory offset relative to the rigid wheel option. At the same time, to ensure motion stability, it is expedient to provide the following condition:

$$\frac{m_1}{k_{\mu 01}}=\frac{m_2}{k_{\mu 02}}$$

(42)

where k_μ01,2 = the lateral stiffness coefficient at static deformation [22] (generally depends on deformation and tire pressure) and m₁, m₂ = gross masses distributed on every single front and rear suspension, correspondingly.

Then, for each wheel, the lateral deformation can be defined as

$${\mathsf{\Delta }}_{\mu }=-{\mathsf{\Delta }F}_{\mu }/k_{\mu }$$

(43)

In the process of tire vertical deformation, the height of its profile changes and, accordingly, the lateral stiffness change too. Usually, the lateral stiffness factor is specified for quasi-static deformation at nominal wheel load. If a shear model is used to estimate the lateral deformation, then the lateral stiffness in the first approximation can be represented as follows

$$k_{\mu }=\frac{G_S}{h_z},\hspace{1em}{k}_{\mu 0}=\frac{G_S}{h_{z0}},\hspace{1em}{h}_z=r_w-\frac{D_r}{2},\hspace{1em}{r}_w=r_{w0}-\mathsf{\Delta }r$$

(44)

where G_S = the so-called rigidity of the cross-sectional area (supposed to be constant), r_w0 = tire’s free radius, Δr = tire’s radial deformation, h_z0 = tire’s profile height under static load, h_z = tire’s profile dynamic height, and D_r = the diameter of wheel’s rim.

Then

$$k_{\mu }=k_{\mu 0}\frac{h_{z0}}{h_z}=k_{\mu 0}\frac{r_{w0}-\mathsf{\Delta }{r}_0-D_r/2}{r_{w0}-\mathsf{\Delta }r-D_r/2}$$

(45)

Δr₀ = tire’s radial static deformation.

Obviously, at zero dynamic deformation, the lateral stiffness is equal to the initial k_μ0 one. As the deformation Δr decreases, it leads to a decrease in lateral stiffness, and vice versa.

2.3.5. Evaluation of Lateral Reactions

Let us define the wheel’s traction force. Generally, the balance equation of rotational dynamics for the j-th wheel is [21]

$$R_{\zeta j}=T_j/r_{wj}-R_{zj}{f}_{rj}-I_{wj}{a}_{\zeta j}/r_{ej}/r_{wj}$$

(46)

where for the j-th wheel: T_j = drive torque, r_wj = dynamic radius, r_ej = effective radius, f_rj = rolling resistance, and I_wj = inertia of rotating masses tied and reduced to a wheel.

The rolling resistance at each wheel can be calculated from the dependence [23]

$$f_{rj}=q_{sy1}+q_{sy3}\left|v_{\zeta j}\right|+q_{sy4}{\left(v_{\zeta j}/v_m\right)}^4$$

(47)

$v_{\zeta j}$ = longitudinal speed in j-th wheel local coordinate system, $v_m$ = speed at which the empirical measurements were made, and q_sy1, q_sy3, q_sy4 = coefficients [23].

The dynamic radius is evaluated by expressions [23]

$$r_{wj}=r_{w0}-\mathsf{\Delta }{r}_j,\hspace{1em}\mathsf{\Delta }{r}_j=R_{zj}/C_{fzj},\hspace{1em}{C}_{fzj}=\frac{F_{zj0}}{r_{w0}}\sqrt{q_{fz1}^2+4q_{fz2}^2}$$

(48)

where Δr_j = radial deformation of the j-th wheel, C_fzj = the j-th tire’s radial stiffness, F_zj0 = static load on the j-th wheel, and q_fz1, q_fz2 = coefficients [23].

Centrifugal increment of the j-th wheel radius [23]

$$d{r}_j=q_{v1}{r}_{w0}{\left(\frac{{\omega }_{wj}{r}_{w0}}{v_m}\right)}^2$$

(49)

The effective radius of the j-th wheel [22]

$$r_{ej}=r_{w0}+d{r}_j-\frac{F_{zj0}}{C_{fz}}\left(D_{ref}arctan\left(B_{ref}\mathsf{\Delta }{r}_j\frac{C_{fz}}{F_{zj0}}\right)+F_{ref}d{r}_j\frac{C_{fz}}{F_{zj0}}\right)$$

(50)

where B_ref, D_ref, F_ref = coefficients [23].

Denote vectors

$${\mathit{T}}_w=\left(\begin{array}{c}{T}_1\\ \begin{array}{c}{T}_2\\ T_3\\ T_4\end{array}\end{array}\right),{\mathit{I}}_w=\left(\begin{array}{c}{I}_{w1}\\ \begin{array}{c}{I}_{w2}\\ I_{w3}\\ I_{w4}\end{array}\end{array}\right),{\mathit{f}}_r=\left(\begin{array}{c}{f}_{r1}\\ \begin{array}{c}{f}_{r2}\\ f_{r3}\\ f_{r4}\end{array}\end{array}\right),{\mathit{a}}_{w\zeta }=\left(\begin{array}{c}{a}_{\zeta 1}\\ \begin{array}{c}{a}_{\zeta 2}\\ a_{\zeta 3}\\ a_{\zeta 4}\end{array}\end{array}\right),{\mathit{r}}_e=\left(\begin{array}{c}{r}_{e1}\\ \begin{array}{c}{r}_{e2}\\ r_{e3}\\ r_{e4}\end{array}\end{array}\right),{\mathit{r}}_w=\left(\begin{array}{c}{r}_{w1}\\ \begin{array}{c}{r}_{w2}\\ r_{w3}\\ r_{w4}\end{array}\end{array}\right)$$

(51)

Then, for all wheels based on Equation (46), we can write

$${\mathit{R}}_{\zeta }={\mathit{i}}_{rw}{\mathit{T}}_w-{\mathit{i}}_{fr}{\mathit{R}}_z-{{\mathit{i}}_{rw}\mathit{i}}_{re}{{\mathit{D}}_{Iw}\mathit{a}}_{w\zeta }={\mathit{i}}_{rw}{\mathit{T}}_w-{\mathit{i}}_{rw}\mathit{N}$$

(52)

where

$$\begin{gathered} {\mathit{i}}_{rw}={diag\left({\mathit{r}}_w\right)}^{-1},\hspace{1em}{\mathit{i}}_{fr}=diag\left({\mathit{f}}_r\right),\hspace{1em}{\mathit{i}}_{re}={diag\left({\mathit{r}}_e\right)}^{-1},\hspace{1em}{\mathit{D}}_{Iw}=diag\left({\mathit{I}}_w\right), \\ \mathit{N}=\left({{\mathit{i}}_{rw}^{-1}\mathit{i}}_{fr}{\mathit{R}}_z+{\mathit{i}}_{re}{{\mathit{D}}_{Iw}\mathit{a}}_{w\zeta }\right) \end{gathered}$$

Now, let us consider the balance of dynamics. Since the transmission torque T determines the traction forces R_ζj on the wheels, and the lateral reactions R_μ are unknown, we need five equations. The first three equations express the balances of forces and moments relative to the vehicle axes. Denote

$$\begin{gathered} {\mathit{k}}_{\zeta }=\left(\begin{array}{ccc}\begin{array}{c}{c}_1\\ s_1\\ {\zeta }_1{s}_1-{\mu }_1{c}_1\end{array}& \begin{array}{c}1\\ 0\\ -{\mu }_2\end{array}& \begin{array}{cc}\begin{array}{c}{c}_3\\ s_3\\ {\zeta }_3{s}_3-{\mu }_3{c}_3\end{array}& \begin{array}{c}1\\ 0\\ -{\mu }_4\end{array}\end{array}\end{array}\right),\hspace{1em}{\mathit{k}}_a=diag\left(\begin{array}{c}m\\ m\\ I\end{array}\right), \\ {\mathit{k}}_{\mu }=\left(\begin{array}{ccc}\begin{array}{c}{-s}_1\\ c_1\\ {\zeta }_1{c}_1+{\mu }_1{s}_1\end{array}& \begin{array}{c}0\\ 1\\ {\zeta }_2\end{array}& \begin{array}{cc}\begin{array}{c}{-s}_3\\ c_3\\ {\zeta }_3{c}_3+{\mu }_3{s}_3\end{array}& \begin{array}{c}0\\ 1\\ {\zeta }_4\end{array}\end{array}\end{array}\right),\hspace{1em}\mathit{a}=\left(\begin{array}{c}{a}_{\zeta }\\ a_{\mu }\\ \epsilon \end{array}\right),\hspace{1em}{\mathit{k}}_{ad}=-\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right) \end{gathered}$$

(53)

where m = vehicle gross mass, I = vehicle inertia relative to the mass center, and ζ_j, μ_j = coordinates of the j-th wheel center in the vehicle coordinate system.

Then, the plane dynamics equations can be represented in the vector-matrix form

$${\mathit{k}}_{\zeta }{\mathit{R}}_{\zeta }+{\mathit{k}}_{\mu }{\mathit{R}}_{\mu }+{\mathit{k}}_{ad}{P}_{ad}={\mathit{k}}_a\mathit{a}$$

(54)

where P_ad = aerodynamic drag force.

Due to the lack of equations, it is necessary to consider the redistribution of the lateral reactions between the same axle’s wheels. Considering Equation (41), we can introduce the ratios

$$\frac{R_{\mu 1}}{R_{\mu 3}}=\frac{f\left({\delta }_1\right)}{f\left({\delta }_3\right)}=k_{\mu 13}\left({\delta }_1,{\delta }_3\right),\hspace{1em}\frac{R_{\mu 2}}{R_{\mu 4}}=\frac{f\left({\delta }_2\right)}{f\left({\delta }_4\right)}=k_{\mu 24}\left({\delta }_2,{\delta }_4\right)$$

(55)

Either

$$R_{\mu 1}-R_{\mu 3}{k}_{\mu 13}\left({\delta }_1,{\delta }_3\right)=0,\hspace{1em}{R}_{\mu 2}-R_{\mu 2}{k}_{\mu 24}\left({\delta }_2,{\delta }_4\right)=0$$

(56)

Denote

$${\mathit{R}}_{\mu }=\left(\begin{array}{c}{R}_{\mu 1}\\ \begin{array}{c}{R}_{\mu 2}\\ R_{\mu 3}\\ R_{\mu 4}\end{array}\end{array}\right),\hspace{1em}\mathit{\delta }=\left(\begin{array}{c}{\delta }_1\\ \begin{array}{c}{\delta }_2\\ {\delta }_3\\ {\delta }_4\end{array}\end{array}\right),\hspace{1em}{\mathit{k}}_{R\mu }=\left(\begin{array}{ccc}\begin{array}{c}1\\ 0\end{array}& \begin{array}{c}0\\ 1\end{array}& \begin{array}{cc}\begin{array}{c}-k_{\mu 13}\\ 0\end{array}& \begin{array}{c}0\\ {-k}_{\mu 24}\end{array}\end{array}\end{array}\right),\hspace{1em}\mathrm{and}\hspace{1em}{\mathit{k}}_{R\mu }{\mathit{R}}_{\mu }={\mathbf{0}}_2$$

(57)

Note that the coefficients of lateral forces’ ratios depend nonlinearly on the distribution of vertical loads, traction forces, and slip angles.

2.4. Variants for Organizing Distributed Traction

2.4.1. Variants of Wheel Torque Control

As a part of the objective of distributing the torque in the transmission of an all-wheel drive autonomous vehicle, four options can be implemented (Figure 4). Let us assume that the vehicle is equipped with a self-locking center differential (Audi Quattro), which does not require electronic control and distributes torques in a certain ratio depending on the output shafts’ resistances. Let us assume that at some motion moment, the longitudinal and lateral accelerations are positive, which causes a redistribution of vertical reactions to the rear axle and right side. We also make an important suggestion that the maximum friction coefficient φ_max is the same for all wheels, that is, an unstable movement mode of any wheel occurs when the geometric sum of the longitudinal and lateral adhesion coefficients approaches the limit. Within the motion planning, we do not consider here the cases of control when the maximum possible adhesions on all wheels are different.

Let us consider the structural elements of transmissions in Figure 4 in more detail. The scheme in Figure 4a shows a variant of all-wheel drive with two symmetrical inter-wheel differentials, where the carriers’ drive torques are conditionally divided in half between the wheels’ semi-axles. In this scheme, to create an additional yaw moment stabilizing the vehicle steerability, it is necessary to slightly brake the wheels being in the worst adhesion conditions, that is, internal in this case. Due to the braking torques B_fl, B_rl, it is possible to increase the drive moments of the axles and realize greater traction on the outer wheels.

The scheme in Figure 4b shows the case of using a sport differential (Tork Vectoring). Here, the front wheels are driven through a symmetrical differential, and the yaw moment is provided by redistributing the rear axle’s torques with a controlled (sport) differential. Since, for the case described, the outer rear wheel has more traction potential, by activating the hydraulic clutch larger torque is transferred to the wheel moving with the higher angular speed. At the same time, the torque on the inner wheel’s shaft is reduced. The difference in traction forces on the rear axle wheels produces a needed yaw moment.

The option in Figure 4c combines the rear axle’s sport differential and the front inner wheel’s brake activation. The variant in Figure 4d implements the full torque vectoring by using two inter-axle sport differentials. Structurally, this option is the most complex and expensive; however, obviously, it allows obtaining the best efficiency in torque redistribution.

2.4.2. Distribution Devices

Let us consider the transmission mechanisms ensuring the torque redistribution along with their dynamic and kinematic models.

Inter-axle (center) differential. Figure 5a shows the Audi Quattro self-locking differential belonging to the limited-sip (LSD) class. The different contact radii R_f, R_r of the output shafts’ crown gears (Figure 5b) give the initial asymmetry to the differential mechanism (DM) design. This leads to distributing the gearing torques with a ratio equal to 40/60 [24]. In normal conditions (both output shafts rotate with the same angular speed), the rear axle’s driving shaft transmits about 60% of the total torque, and the front axle −40%, correspondingly. Each crown gear has a friction clutch connecting it with the differential’s carrier. The relative sliding of friction elements generates moments that are redistributed according to the mechanism’s kinematic state. Thus, the output shaft’s torque may be reduced or increased. The clutch packs may be installed with precompression, stipulating the static friction torque by the pressing washers. With a decrease of one shaft’s resistance, its angular speed becomes greater than the carrier’s angular speed. The excessive power flow from an outrunning shaft returns to the carrier by the friction moment and increases the lagging shaft’s torque. A feature of this DM is the dynamic redistribution of friction moments due to the axial component of the gearing reaction. The greater the satellite’s force (Figure 5c), the greater the compression, and in turn, the greater the moment that can be passed by friction clutches.

The general kinematic relations for the asymmetric DM depicted in Figure 5, using the well-known Willis formula, may be reflected as follows

$${\omega }_D=\frac{{\omega }_f+g_{fr}{\omega }_r}{1+g_{fr}}$$

(58)

where ω_D, ω_f, and ω_r = angular velocities of the differential carrier, front, and rear output shafts, respectively; and g_fr = R_r/R_f = 60⁄40 = 3⁄2—algebraic ratio from the front drive axle to the rear drive axle.

The front and rear axles’ driving shafts’ angular speeds ω_f, ω_r can be estimated based on the symmetrical DM’s kinematic properties. Assuming that all wheels are rotating in approximately the same mode relative to the tires’ longitudinal slip, we can calculate with sufficient accuracy

$${\omega }_f=i_f\frac{{\omega }_{fl}+{\omega }_{fr}}{2}=\frac{i_f}{2}\left(\frac{v_{\zeta 1}}{r_{e1}}+\frac{v_{\zeta 3}}{r_{e3}}\right),\hspace{1em}{\omega }_r=i_f\frac{{\omega }_{rl}+{\omega }_{rr}}{2}=\frac{i_f}{2}\left(\frac{v_{\zeta 2}}{r_{e2}}+\frac{v_{\zeta 4}}{r_{e4}}\right)$$

(59)

where ω_fl, ω_fr = angular speeds of the front left and right wheels, respectively, ω_rl, ω_rr = angular speeds of the rear left and right wheels, respectively, i_f = final gear ratio, $v_{\zeta j}$ = longitudinal speed of the j-th wheel center in the wheel’s local coordinate system, and r_ej = wheel’s effective radius, j = [1, …, 4].

Thus, if during curvilinear motion, each wheel follows its trajectory without critical longitudinal slip, then, consequently, a certain ratio of output moments is established in the inter-axle (center) differential. This ratio is determined by the difference in angular speeds ω_D, ω_f, ω_r. Note that the friction discs operate in an oily environment, which contributes to a smooth redistribution of output torques without the effects of sticking, and jamming, which are typical for cases close to dry friction. That is, to describe the frictional interaction, we may use an expression such as

$${\mu }_{Fk}={\mu }_k\mathit{tanh}\left(c_k\left({\omega }_D-{\omega }_k\right)\right),$$

(60)

where μ_Fk = actual friction coefficient between friction pairs, μ_k = modulus of maximum friction coefficient, c_k = intensity factor, k = f (front), and r (rear).

For each k = f, r shaft’s clutch, the dynamic friction torque F_kd is modeled as follows

$$F_{dk}=n_s{P}_{Nk}{n}_{Fk}{R}_{Fk}{\mu }_{Fk}=n_s{P}_{Nk}{R}_k\frac{R_{Fk}}{R_k}{n}_{Fk}{\mu }_{Fk}$$

(61)

where n_s = number of satellites, P_Nk = axial gearing force, n_Fk = number of friction couples, R_k = average radius of side gear, and R_Fk = average friction radius.

Since the axial component P_N of the gearing contact reaction, P depends neither on the direction of the shaft rotation nor on the direction of the tangential component P_T, the friction moment’s dynamic part is written as follows

$$F_{dk}=n_s{R}_k\left|P_{Tk}\right|\mathit{tan}\left(\alpha \right)\frac{R_{Fk}}{R_k}{n}_{Fk}{\mu }_{Fk}=\left|M_k\right|{\eta }_k,\hspace{1em}k=f,r$$

(62)

where P_Tk = tangential gearing force, α = gearing angle, M_k = side gear torque, and η_k = friction factor.

$$M_k=n_s{R}_k{P}_{Tk},\hspace{1em}{\eta }_k=\mathit{tan}\left(\alpha \right)\frac{R_{Fk}}{R_k}{n}_{Fk}{\mu }_{Fk}$$

(63)

This DM also allows precompressing the friction clutches, which provides a static friction moment

$$F_{sk}=P_{N0k}{R}_{Fk}{n}_{Fk}{\mu }_{Fk}$$

(64)

where P_N0k = preliminary established axial compression force.

The preset friction moment can also be expressed in terms of the equivalent moment M_k0 on the side gears.

$$F_{sk}=\frac{P_{N0k}{R}_k}{\mathit{tan}\left(\alpha \right)}\mathit{tan}\left(\alpha \right)\frac{R_{Fk}}{R_k}{n}_{Fk}{\mu }_{Fk}=\frac{P_{N0k}{R}_k}{\mathit{tan}\left(\alpha \right)}{\eta }_k={M_{0k}\eta }_k,\hspace{1em}k=f,r$$

(65)

Consider the distribution of moments without inertial forces for simplification. Let us introduce vector and matrix denotations

$$\begin{gathered} \mathit{T}=\left(\begin{array}{c}{T}_f\\ T_r\end{array}\right),\hspace{1em}\mathit{M}=\left(\begin{array}{c}{M}_f\\ M_r\end{array}\right),\hspace{1em}{\mathit{F}}_d=\left(\begin{array}{c}{F}_{df}\\ F_{dr}\end{array}\right),\hspace{1em}{\mathit{F}}_s=\left(\begin{array}{c}{F}_{sf}\\ F_{sr}\end{array}\right),\hspace{1em}{\mathit{M}}_0=\left(\begin{array}{c}{M}_{0f}\\ M_{0f}\end{array}\right),\hspace{1em}{\mathit{t}}_D=\left(\begin{array}{c}0.4\\ 0.6\end{array}\right), \\ {\mathit{E}}_2=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\hspace{1em}\mathit{\eta }=\left(\begin{array}{cc}{\eta }_f& 0\\ 0& {\eta }_r\end{array}\right),\hspace{1em}{\mathit{t}}_M=\left(\begin{array}{cc}{\mathit{t}}_D& {\mathit{t}}_D\end{array}\right) \end{gathered}$$

(66)

where vector t_D is associated with initial torque redistribution.

The torques on the output shafts

$$\mathit{T}=\mathit{M}+{\mathit{F}}_d+{\mathit{F}}_s=\left({\mathit{E}}_2+\mathit{\eta }\right)\mathit{M}+{\mathit{\eta }\mathit{M}}_0$$

(67)

On the other hand, the drive torques on the side gears can be expressed in terms of the moment on the differential carrier T_D

$$\mathit{M}={\mathit{t}}_D{T}_D-{\mathit{t}}_M{\mathit{F}}_d-{\mathit{t}}_M{\mathit{F}}_s={\mathit{t}}_D{T}_D-{\mathit{t}}_M\mathit{\eta }\mathit{M}-{\mathit{t}}_M{\mathit{\eta }\mathit{M}}_0$$

(68)

Whence

$$\left({\mathit{E}}_2+{\mathit{t}}_M\mathit{\eta }\right)\mathit{M}={\mathit{t}}_D{T}_D-{\mathit{t}}_M{\mathit{\eta }\mathit{M}}_0,\hspace{1em}\mathit{M}={\left({\mathit{E}}_2+{\mathit{t}}_M\mathit{\eta }\right)}^{-1}\left({\mathit{t}}_D{T}_D-{\mathit{t}}_M{\mathit{\eta }\mathit{M}}_0\right)$$

(69)

Substituting M in Equation (68), obtain

$$\mathit{T}=\left({\mathit{E}}_2+\mathit{\eta }\right){\left({\mathit{E}}_2+{\mathit{t}}_M\mathit{\eta }\right)}^{-1}\left({\mathit{t}}_D{T}_D-{\mathit{t}}_M{\mathit{\eta }\mathit{M}}_0\right)+{\mathit{\eta }\mathit{M}}_0$$

(70)

Denote

$${\mathit{t}}_{Ef}=\left({\mathit{E}}_2+\mathit{\eta }\right){\left({\mathit{E}}_2+{\mathit{t}}_M\mathit{\eta }\right)}^{-1}$$

(71)

Thus, if the kinematic distribution of the wheels’ angular speeds is given, the relations between the output shafts’ torques and the differential carrier’s torque

$$\mathit{T}={\mathit{t}}_{Ef}{\mathit{t}}_D{T}_D+\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{\mathit{M}}_0$$

(72)

Note that the sensitivity of the moment distribution depends on the coefficient c_k.

Inter-wheel Sport DM must ensure transmitting a larger torque to an outrunning axle, which causes additional cornering (yaw moment). Such a solution can be obtained using the design scheme and functioning of Audi Sport Differential shown in Figure 6.

The drive is carried out over the final gear pinion 1. The ring gear is rigidly connected to the differential carrier 2. Satellites 3, rotating around the axis fixed in the differential carrier, interact with side gears that drive the output axles 4 and 5 by the slots. The differential corps (carrier) has gear rims on the end sides for driving by internal gearing the coupling halves 6 and 7, in which rotational axes are, respectively, shifted relative to the carrier rotational axis. Using the frictional packs, the half-clutches 6 and 7 can drive the half-clutches 8 and 9, respectively, which are connected by internal gearing to the output shafts 4 and 5. Toroidal hydraulic cylinders 10 and 11 are installed from each side to act on the clutch packs using the pressure p₁₀ and p₁₁. Thus, by activating the required hydraulic cylinder, part of the carrier torque may be passed to the needed semi-axle using the frictional adhesion between the half-couplings over the two-step internal gearing (i₂₈ = 0.8752, i₂₆ = 1.2258, i₈₄ = 0.714).

Since the Sport DM is formed based on a symmetrical one, it is possible to relate the moments on the output semi-axes T_l, T_r with the moment T_d on the DM carrier and the moment ΔT_lr of the friction clutch.

$$\left(\begin{array}{c}{T}_l\\ T_r\end{array}\right)={\mathit{t}}_d\left(\begin{array}{c}{T}_d\\ \mathsf{\Delta }{T}_{lr}\end{array}\right),\hspace{1em}{\mathit{t}}_d=\frac{1}{2}\left(\begin{array}{cc}1& -1\\ 1& 1\end{array}\right),\hspace{1em}\left(\begin{array}{c}{T}_d\\ \mathsf{\Delta }{T}_{lr}\end{array}\right)={\mathit{t}}_s\left(\begin{array}{c}{T}_l\\ T_r\end{array}\right),\hspace{1em}{\mathit{t}}_s=\left(\begin{array}{cc}1& 1\\ -1& 1\end{array}\right)$$

(73)

Thus, the required friction clutch’s moment is equal to the difference of driving torques on the wheels’ semi-axles. Knowing the difference in the output shafts’ angular speeds, it is possible to calculate the required effect of the clutch drive.

2.4.3. Drivetrain Dynamics

Now, let us define the distribution and realization of traction forces on the wheels, considering the transmission type (Figure 4).

Since the drive of all wheels is realized by the same transmission, the drive torques can be linked to each other depending on the drive type. Let us consider options for using drives with torque vectoring (TV) between axles and wheels (Figure 4). Obviously, maximum vehicle controllability and stability are ensured when all the driving wheels work at the approximately same adhesion conditions. If we assume that each wheel realizes the same potential of longitudinal adhesion, then, consequently, all wheels will have a close potential of adhesion in the lateral direction. Let us search for distributing the transmission torques based on this ground.

Rewrite Equation (52) in the form regarding moments.

$$T_j=R_{zj}{\phi }_{\zeta j}{r}_{wj}+N_j,\hspace{1em}{N}_j=R_{zj}{f}_{rj}{r}_{wj}+I_{wj}{a}_{\zeta j}/r_{ej}$$

(74)

where j = [1, 3] for the front axle, and j = [2, 4] for the rear axle.

Then, the ratios k_f, k_r of the moments on the front and rear axles’ wheels

$$\frac{R_{z1}{\phi }_{\zeta 1}{r}_{w1}}{R_{z3}{\phi }_{\zeta 3}{r}_{w3}}=\frac{T_1-N_1}{T_3-N_3}=k_f,\hspace{1em}\frac{R_{z2}{\phi }_{\zeta 2}{r}_{w2}}{R_{z4}{\phi }_{\zeta 4}{r}_{w4}}=\frac{T_2-N_2}{T_4-N_4}=k_r$$

(75)

As mentioned, we are trying to reach such a distribution of moments at which the adhesion potentials are the same and do not depend on the wheels’ loads. Then, optimally, φ_ζ1 = φ_ζ3 = φ_ζf, which, after reduction, leads to equations

$$\frac{R_{z1}{r}_{w1}}{R_{z3}{r}_{w3}}=k_f,\hspace{1em}{T}_1-k_f{T}_3=N_1-k_f{N}_3$$

(76)

Similarly, for the rear axle, based on the same principle that φ_ζ2 = φ_ζ4 = φ_ζr,

$$\frac{R_{z2}{r}_{w2}}{R_{z4}{r}_{w4}}=k_r,\hspace{1em}{T}_2-k_r{T}_4=N_2-k_r{N}_4$$

(77)

Note that with known distributions of speed and vertical reactions, the components N_j, as well as the coefficients k_f, k_r, can be obtained in advance.

Note that the wheels’ drive torques T₁, T₃ are generally not equivalent to the drive torques of the output shafts of the symmetrical differential T_fl = T_fr, nor are the torques T₂, T₄ with respect to the torques T_rl = T_rr.

Now, let us form a vector of driving moments for the cases using different DMs.

Case A. Front and rear symmetrical differentials share the torques T_f, T_r equally between the wheels’ semi-axles, i.e., T_fl = T_fr, T_rl = T_rr. The correction is fulfilled by activating the brake mechanisms, which stipulates the braking moments’ vector B. Then, the wheels’ torques will be related to semi-axle torques as follows

$${\mathit{T}}_w={\mathit{t}}_p{\mathit{T}}_a-\mathit{B}={\mathit{t}}_p{\mathit{t}}_T{i}_f\mathit{T}-\mathit{B}$$

(78)

where t_T is responsible for distributing the torque between the drive semi-axles, and t_p is the matrix for permuting the T_a vector’s elements

$${\mathit{T}}_w=\left(\begin{array}{c}{T}_1\\ \begin{array}{c}{T}_2\\ T_3\\ T_4\end{array}\end{array}\right),\hspace{1em}{\mathit{T}}_a=\left(\begin{array}{c}{T}_{fl}\\ \begin{array}{c}{T}_{fr}\\ T_{rl}\\ T_{rr}\end{array}\end{array}\right),\hspace{1em}{\mathit{t}}_p=\left(\begin{array}{ccc}\begin{array}{c}1\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}1\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}0\\ 1\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\end{array}\right),$$

(79)

$${\mathit{t}}_T=\frac{1}{2}\left(\begin{array}{cc}1& 0\\ 1& 0\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}1\\ 1\end{array}\end{array}\right),\hspace{1em}\mathit{B}=\left(\begin{array}{c}{B}_1\\ \begin{array}{c}{B}_2\\ B_3\\ B_4\end{array}\end{array}\right),\hspace{1em}\mathit{T}=\left(\begin{array}{c}{T}_f\\ T_r\end{array}\right)$$

At the same time, such a control vector B must be provided that satisfies the requirements of redistributing the moments from Equations (76) and (77) and minimizing the control in general.

$$\left\{\begin{array}{c}{T}_1-k_f{T}_3=\left(T_{fl}-B_1\right)-k_f\left(T_{fr}-B_3\right)\\ T_2-k_r{T}_4=\left(T_{rl}-B_2\right)-k_r\left(T_{rr}-B_4\right)\\ B_1{B}_3=0,B_1\ge 0,B_3\ge 0\\ B_2{B}_4=0,B_2\ge 0,B_4\ge 0\end{array}\right.$$

(80)

Wheel torques in Equation (80) may be replaced according to Equations (76) and (77).

$$\left\{\begin{array}{c}{N}_1-k_f{N}_3=\left(T_{fl}-B_1\right)-k_f\left(T_{fr}-B_3\right)\\ N_2-k_r{N}_4=\left(T_{rl}-B_2\right)-k_r\left(T_{rr}-B_4\right)\end{array}\right.$$

(81)

Then, the combination of additional linear and nonlinear conditions for braking torques will have the form

$${\mathit{t}}_{AB}\mathit{B}={\mathit{t}}_A\left({\mathit{t}}_p{\mathit{T}}_a-\mathit{N}\right)={\mathit{t}}_A\left({\mathit{t}}_p{\mathit{t}}_T{i}_f\mathit{T}-\mathit{N}\right)$$

(82)

where

$${\mathit{t}}_A=\left(\begin{array}{ccc}\begin{array}{c}-1\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ -1\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}{k}_f\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ k_r\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}\right),\hspace{1em}{\mathit{t}}_B=diag\left(\begin{array}{c}0\\ \begin{array}{c}0\\ B_1\\ B_2\end{array}\end{array}\right),\hspace{1em}{\mathit{t}}_{AB}={\mathit{t}}_A+{\mathit{t}}_B$$

(83)

Note that we are only interested in positive values of the braking torques.

Substituting the expression from Equation (72) to Equation (78), we obtain

$${\mathit{T}}_w={\mathit{t}}_p{\mathit{t}}_T{i}_f\left({\mathit{t}}_{Ef}{\mathit{t}}_D{T}_D+\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{\mathit{M}}_0\right)-\mathit{B}$$

(84)

Denote

$${\mathit{k}}_{TD}={\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_T{\mathit{t}}_{Ef}{\mathit{t}}_D{i}_f,\hspace{1em}{\mathit{k}}_{M0}={\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_T\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{i}_f$$

(85)

Substituting coefficients in Equation (52) yields

$${\mathit{R}}_{\zeta }={\mathit{i}}_{rw}{\mathit{T}}_w-{\mathit{i}}_{rw}\mathit{N}={\mathit{k}}_{TD}{T}_D+{\mathit{k}}_{M0}{\mathit{M}}_0-{\mathit{i}}_{rw}\mathit{N}-{\mathit{i}}_{rw}\mathit{B}$$

(86)

Let us transform additional conditions

$${\mathit{t}}_{AB}\mathit{B}={\mathit{t}}_A{\mathit{t}}_p{\mathit{t}}_T{i}_f\left({\mathit{t}}_{Ef}{\mathit{t}}_D{T}_D+\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{\mathit{M}}_0\right)-{\mathit{t}}_A\mathit{N}$$

(87)

Denote

$${\mathit{k}}_{BT}={\mathit{t}}_A{\mathit{t}}_p{\mathit{t}}_T{\mathit{t}}_{Ef}{\mathit{t}}_D{i}_f,\hspace{1em}{\mathit{k}}_{BM}={\mathit{t}}_A{\mathit{t}}_p{\mathit{t}}_T\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{i}_f$$

(88)

Substituting the expression from Equation (88) to Equation (87), we obtain additional conditions compactly

$${\mathit{t}}_{AB}\mathit{B}={\mathit{k}}_{BT}{T}_D+{\mathit{k}}_{BM}{\mathit{M}}_0-{\mathit{t}}_A\mathit{N}$$

(89)

Thus, five unknowns are formed, such as torque T_D and braking torques B.

Case B. The front symmetrical differential divides the torque T_f equally between the semi-axles without further correction, i.e., T_fl = T_fr = T₁ = T₃. The rear sport differential redistributes the drive torque T_r according to Equation (73). Then, the drive moments on the wheels will be related as follows

$${\mathit{T}}_f=\left(\begin{array}{c}{T}_{fl}\\ T_{fr}\end{array}\right)={\mathit{E}}_{2\times 1}\frac{i_f}{2}{T}_f,\hspace{1em}{\mathit{E}}_{2\times 1}=\left(\begin{array}{c}1\\ 1\end{array}\right),\hspace{1em}{\mathit{T}}_r=\left(\begin{array}{c}{T}_{rl}\\ T_{rr}\end{array}\right)={\mathit{t}}_d\left(\begin{array}{c}{T}_r{i}_f\\ \mathsf{\Delta }{T}_r\end{array}\right)$$

(90)

Then, the vector of moments on the wheels’ drive shafts

$${\mathit{T}}_a=\left(\begin{array}{c}{\mathit{T}}_f\\ {\mathit{T}}_r\end{array}\right)=\underset{{\mathit{t}}_{fr}}{\underbrace{\left(\begin{array}{cc}{\frac{1}{2}\mathit{E}}_{2\times 1}& {\mathbf{0}}_{2\times 2}\\ {\mathbf{0}}_{2\times 1}& {\mathit{t}}_d\end{array}\right)}}\left(\begin{array}{c}{i}_f{T}_f\\ \begin{array}{c}{i_{f}T}_r\\ \mathsf{\Delta }{T}_r\end{array}\end{array}\right)={\mathit{t}}_{fr}{\mathit{T}}_t,$$

(91)

where the matrix t_fr = responsible for proportions of the distribution of moments, and

$${\mathit{T}}_t=\left(\begin{array}{c}{i}_f{T}_f\\ \begin{array}{c}{i_{f}T}_r\\ \mathsf{\Delta }{T}_r\end{array}\end{array}\right)=\left(\begin{array}{c}{i}_f\mathit{T}\\ T_{rr}-T_{rl}\end{array}\right)=\left(\begin{array}{c}{i}_f\mathit{T}\\ \left(1-k_r\right)T_{rr}{-N}_2+k_r{N}_4\end{array}\right)$$

(92)

Using Equation (92) and expanding T_t by components, we obtain

$${\mathit{T}}_a={\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\mathit{T}+{\mathit{t}}_a{\mathit{T}}_a+{\mathit{t}}_N\mathit{N}\right)$$

(93)

where, considering E₂ to be a 2 × 2 identity matrix,

$${\mathit{t}}_T=\left(\begin{array}{c}{\mathit{E}}_2\\ {\mathbf{0}}_{1\times 2}\end{array}\right),\hspace{1em}{\mathit{t}}_a=\left(\begin{array}{ccc}\begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 1-k_r\end{array}\end{array}\end{array}\right),\hspace{1em}{\mathit{t}}_N=\left(\begin{array}{ccc}\begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ -1\end{array}& \begin{array}{cc}\begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ k_r\end{array}\end{array}\end{array}\right),\hspace{1em}\mathit{N}=\left(\begin{array}{c}{N}_1\\ \begin{array}{c}{N}_2\\ N_3\\ N_4\end{array}\end{array}\right)$$

Thus, the vector of semi-axles’ moments is expressed as

$${\mathit{T}}_a={\left({\mathit{E}}_4-{\mathit{t}}_{fr}{\mathit{t}}_a\right)}^{-1}{\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\mathit{T}+{\mathit{t}}_N\mathit{N}\right)={\mathit{t}}_e{\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\mathit{T}+{\mathit{t}}_N\mathit{N}\right)$$

(94)

Denote

$${\mathit{t}}_e={\left({\mathit{E}}_4-{\mathit{t}}_{fr}{\mathit{t}}_a\right)}^{-1}$$

(95)

where vector t_p ensures permuting the wheels’ numbering order, E₄ = 4 × 4 identity matrix.

Then, the vector of wheels’ moments

$${\mathit{T}}_w={\mathit{t}}_p{\mathit{T}}_a={\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\mathit{T}+{\mathit{t}}_N\mathit{N}\right)$$

(96)

Substituting instead of T its expression from Equation (72) for the inter-axle differential drive, we obtain

$$\begin{array}{c}{\mathit{T}}_w={\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\left({\mathit{t}}_{Ef}{\mathit{t}}_D{T}_D+\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right){\mathit{\eta }\mathit{M}}_0\right)+{\mathit{t}}_N\mathit{N}\right)\hfill \\ ={{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\mathit{t}}_{Ef}{\mathit{t}}_D{i}_f{T}_D+{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{i}_f{\mathit{M}}_0+{{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}\mathit{t}}_N\mathit{N}\hfill \end{array}$$

(97)

Considering Equations (52) and (97), traction forces on the wheels can be calculated as follows

$$\begin{gathered} {\mathit{R}}_{\zeta }={\mathit{i}}_{rw}{\mathit{T}}_w-{\mathit{i}}_{rw}\mathit{N}={\mathit{i}}_{rw}{{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\mathit{t}}_{Ef}{\mathit{t}}_D{i}_f{T}_D \\ +{\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{i}_f{\mathit{M}}_0+{\mathit{i}}_{rw}\left({{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}\mathit{t}}_N-{\mathit{E}}_4\right)\mathit{N} \end{gathered}$$

(98)

Denote

$$\begin{gathered} {\mathit{k}}_{TD}={\mathit{i}}_{rw}{{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\mathit{t}}_{Ef}{\mathit{t}}_D{i}_f,\hspace{1em}{\mathit{k}}_{M0}={\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{i}_f, \\ {\mathit{k}}_N={\mathit{i}}_{rw}\left({{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}\mathit{t}}_N-{\mathit{E}}_4\right) \end{gathered}$$

(99)

Then, in the compact form

$${\mathit{R}}_{\zeta }={\mathit{k}}_{TD}{T}_D+{\mathit{k}}_{M0}{\mathit{M}}_0+{\mathit{k}}_N\mathit{N}$$

(100)

Case C. The front symmetrical differential divides the moment T_f equally between the semi-axles, i.e., T_fl = T_fr. The correction is made by activating the front brake mechanisms with braking torques B₁, B₃. The rear sport differential redistributes the drive torque T_r according to Equation (73). Then, the drive moments across the semi-axles will be related as follows

$${\mathit{T}}_w={\mathit{t}}_p{\mathit{T}}_a-{\mathit{t}}_b\mathit{B}={\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\mathit{T}+{\mathit{t}}_N\mathit{N}\right)-{\mathit{t}}_b\mathit{B}$$

(101)

$${\mathit{T}}_a={\mathit{t}}_e{\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\mathit{T}+{\mathit{t}}_N\mathit{N}\right)$$

(102)

By analogy with Equation (80)

$$\left\{\begin{array}{c}{T}_1-k_f{T}_3=\left(T_{fl}-B_1\right)-k_f\left(T_{fr}-B_3\right)\\ B_1{B}_3=0,B_1\ge 0,B_3\ge 0\\ B_2=B_4=0\end{array}\right.,\hspace{1em}\mathit{B}=\left(\begin{array}{c}{B}_1\\ B_3\end{array}\right),\hspace{1em}{\mathit{t}}_b=\left(\begin{array}{cc}1& 0\\ 0& 0\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}1\\ 0\end{array}\end{array}\right)$$

(103)

Wheel moments in Equation (103) can be replaced according to Equation (76)

$$N_1-k_f{N}_3=\left(T_{fl}-B_1\right)-k_f\left(T_{fr}-B_3\right)$$

(104)

Then, the combination of additional linear and nonlinear conditions for braking moments will have the form

$${\mathit{t}}_{AB}{\mathit{t}}_b\mathit{B}={\mathit{t}}_A{\mathit{t}}_p{\mathit{T}}_a-{\mathit{t}}_A\mathit{N}={\mathit{t}}_A{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\mathit{T}+{\mathit{t}}_N\mathit{N}\right)-{\mathit{t}}_A\mathit{N}$$

(105)

where

$${\mathit{t}}_A=\left(\begin{array}{ccc}\begin{array}{c}-1\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}{k}_f\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\right),\hspace{1em}{\mathit{t}}_B=\left(\begin{array}{ccc}\begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}0\\ B_1\end{array}& \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\right),\hspace{1em}{\mathit{t}}_{AB}={\mathit{t}}_A+{\mathit{t}}_B$$

(106)

Substituting the expression from Equation (72) to Equation (101), we obtain

$${\mathit{T}}_w={\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\left({\mathit{t}}_{Ef}{\mathit{t}}_D{T}_D+\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{\mathit{M}}_0\right)+{\mathit{t}}_N\mathit{N}\right)-{\mathit{t}}_b\mathit{B}$$

(107)

Denote

$$\begin{gathered} {{\mathit{i}}_{rw}\mathit{T}}_w={\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T{\mathit{t}}_{Ef}{\mathit{t}}_D{i}_f{T}_D+{\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{i}_f{\mathit{M}}_0 \\ +{\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_N\mathit{N}-{\mathit{i}}_{rw}{\mathit{t}}_b\mathit{B}, \\ {\mathit{k}}_{TD}={\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T{\mathit{t}}_{Ef}{\mathit{t}}_D{i}_f,{\mathit{k}}_{M0}={\mathit{i}}_{rw}{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{i}_f, \\ {\mathit{k}}_N={\mathit{i}}_{rw}\left({\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_N\mathit{N}-{\mathit{E}}_4\right),\hspace{1em}{{\mathit{k}}_B=\mathit{i}}_{rw}{\mathit{t}}_b \end{gathered}$$

(108)

Substituting the coefficients in Equation (52), we have

$${\mathit{R}}_{\zeta }={\mathit{i}}_{rw}{\mathit{T}}_w-{\mathit{i}}_{rw}\mathit{N}={\mathit{k}}_{TD}{T}_D+{\mathit{k}}_{M0}{\mathit{M}}_0+{\mathit{k}}_N\mathit{N}-{\mathit{k}}_B\mathit{B}$$

(109)

Denote

$$\begin{gathered} {\mathit{k}}_{BT}={\mathit{t}}_A{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T{\mathit{t}}_{Ef}{\mathit{t}}_D{i}_f,\hspace{1em}{\mathit{k}}_{BM}={\mathit{t}}_A{\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_T\left({\mathit{E}}_2-{\mathit{t}}_{Ef}{\mathit{t}}_M\right)\mathit{\eta }{i}_f, \\ {\mathit{k}}_{BN}={\mathit{t}}_A\left({\mathit{t}}_p{\mathit{t}}_e{\mathit{t}}_{fr}{\mathit{t}}_N-{\mathit{E}}_4\right) \end{gathered}$$

(110)

Then, the additional conditions

$${\mathit{t}}_{AB}{\mathit{t}}_b\mathit{B}={\mathit{k}}_{BT}{T}_D+{\mathit{k}}_{BM}{\mathit{M}}_0+{\mathit{k}}_{BN}\mathit{N}$$

(111)

Case (d). In this scheme, the front and rear sport differentials redistribute the drive torques T_f, T_r following Equation (73). Then, the drive moments along all the semi-axles will be related as follows

$${\mathit{T}}_a=\left(\begin{array}{c}{\mathit{T}}_f\\ {\mathit{T}}_r\end{array}\right)={\mathit{t}}_{fr}{\mathit{T}}_t,\hspace{1em}{\mathit{t}}_{fr}=\left(\begin{array}{cc}{\mathit{t}}_d& {\mathbf{0}}_{2\times 2}\\ {\mathbf{0}}_{2\times 2}& {\mathit{t}}_d\end{array}\right),\hspace{1em}{\mathit{T}}_t=\left(\begin{array}{c}{i_{f}T}_f\\ \begin{array}{c}\mathsf{\Delta }{T}_f\\ {i_{f}T}_r\\ \mathsf{\Delta }{T}_r\end{array}\end{array}\right)$$

(112)

where

$${\mathit{T}}_t=\left(\begin{array}{c}{i_{f}T}_f\\ \begin{array}{c}{T}_{fr}-T_{fl}\\ {i_{f}T}_r\\ T_{rr}-T_{rl}\end{array}\end{array}\right)=\left(\begin{array}{c}{i_{f}T}_f\\ \begin{array}{c}\left(1-k_f\right)T_{fr}-N_1+k_f{N}_3\\ {i_{f}T}_r\\ \left(1-k_r\right)T_{rr}{-N}_2+k_r{N}_4\end{array}\end{array}\right)$$

(113)

Using Equation (113) and expanding T_t by components, we obtain

$${\mathit{T}}_a={\mathit{t}}_{fr}\left({\mathit{t}}_T{i}_f\mathit{T}+{\mathit{t}}_a{\mathit{T}}_a+{\mathit{t}}_N\mathit{N}\right)$$

(114)

where

$${\mathit{t}}_T=\left(\begin{array}{cc}1& 0\\ 0& 0\\ 0& 1\\ 0& 0\end{array}\right),\hspace{1em}{\mathit{t}}_a=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 1-k_f& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1-k_r\end{array}\right),\hspace{1em}{\mathit{t}}_N=\left(\begin{array}{cccc}0& 0& 0& 0\\ -1& 0& k_f& 0\\ 0& 0& 0& 0\\ 0& -1& 0& k_r\end{array}\right)$$

The subsequent technique corresponds to Equations (94)–(100).

Now, for drive variants in Figure 4, let us form systems of dynamics equations with additional conditions based on Section 2.2. Thus, each variant will have its own set of equations describing the traction dynamics and braking effect. Nevertheless, some equations describing the lateral forces will be common for all the schemes.

Case A. The system is nonlinear; therefore, we reformulate the equations of dynamics from Equation (54), substituting the vector R_ζ from Equation (86).

$${\mathit{E}\mathit{q}}_1={\mathit{k}}_{\mu }{\mathit{R}}_{\mu }-{\mathit{k}}_{\zeta }{\mathit{i}}_{rw}\mathit{B}+{\mathit{k}}_{\zeta }{\mathit{k}}_{TD}{T}_D+{\mathit{k}}_{\zeta }{\mathit{k}}_{M0}{\mathit{M}}_0-{{\mathit{k}}_{\zeta }\mathit{i}}_{rw}\mathit{N}-{\mathit{k}}_a\mathit{a}+{\mathit{k}}_{ad}{P}_{ad}={\mathbf{0}}_3$$

(115)

The additional conditions for braking torques in the nonlinear form

$${\mathit{E}\mathit{q}}_2={\mathit{t}}_{AB}\mathit{B}-{\mathit{k}}_{BT}{T}_D-{\mathit{k}}_{BM}{\mathit{M}}_0+{\mathit{t}}_A\mathit{N}={\mathbf{0}}_4$$

(116)

Conditions for nonlinear tire’s lateral reactions based on Equation (41)

$${\mathit{E}\mathit{q}}_3={\mathit{k}}_{R\mu }{\mathit{R}}_{\mu }={\mathbf{0}}_2,\hspace{1em}{\mathit{E}\mathit{q}}_4={\mathit{R}}_{\mu }-f\left(\mathit{\delta }\right)={\mathbf{0}}_4$$

(117)

Then, the system of nonlinear equations and vector of variables are

$$\mathit{E}\mathit{q}=\left(\begin{array}{c}{\mathit{E}\mathit{q}}_1\\ {\mathit{E}\mathit{q}}_2\\ {\mathit{E}\mathit{q}}_3\\ {\mathit{E}\mathit{q}}_4\end{array}\right)=\left(\begin{array}{c}{\mathbf{0}}_3\\ {\mathbf{0}}_4\\ {\mathbf{0}}_2\\ {\mathbf{0}}_4\end{array}\right)={\mathbf{0}}_{13},\hspace{1em}\mathit{x}=\left(\begin{array}{c}{\mathit{R}}_{\mu }\\ \mathit{B}\\ T_D\\ \mathit{\delta }\end{array}\right)$$

(118)

Cases B and D. These variants are distinguished by the generality of expressions and the absence of braking moments as unknowns. The equation of dynamics in the nonlinear form, given Equations (54) and (100).

$${\mathit{E}\mathit{q}}_1={\mathit{k}}_{\zeta }{\mathit{k}}_{TD}{T}_D+{\mathit{k}}_{\zeta }{\mathit{k}}_{M0}{\mathit{M}}_0+{\mathit{k}}_{\zeta }{\mathit{k}}_N\mathit{N}+{\mathit{k}}_{\mu }{\mathit{R}}_{\mu }+{\mathit{k}}_{ad}{P}_{ad}-{\mathit{k}}_a\mathit{a}={\mathbf{0}}_3$$

(119)

In addition, we use Equation (117). Then, the system of nonlinear equations and vector of variables are

$$\mathit{E}\mathit{q}=\left(\begin{array}{c}{\mathit{E}\mathit{q}}_1\\ {\mathit{E}\mathit{q}}_3\\ {\mathit{E}\mathit{q}}_4\end{array}\right)=\left(\begin{array}{c}{\mathbf{0}}_3\\ {\mathbf{0}}_2\\ {\mathbf{0}}_4\end{array}\right)={\mathbf{0}}_9,\hspace{1em}\mathit{x}=\left(\begin{array}{c}{\mathit{R}}_{\mu }\\ T_D\\ \mathit{\delta }\end{array}\right)$$

(120)

Case C. The system is also nonlinear, therefore, we reformulate the equations of dynamic Equation (109), substituting the vector R_ζ from Equation (54)

$${\mathit{k}}_{\mu }{\mathit{R}}_{\mu }-{\mathit{k}}_{\zeta }{\mathit{k}}_B\mathit{B}+{{\mathit{k}}_{\zeta }\mathit{k}}_{TD}{T}_D+{\mathit{k}}_{\zeta }{\mathit{k}}_{M0}{\mathit{M}}_0+{{\mathit{k}}_{\zeta }\mathit{k}}_N\mathit{N}-{\mathit{k}}_a\mathit{a}+{\mathit{k}}_{ad}{P}_{ad}={\mathbf{0}}_3$$

(121)

Also, rewrite Equation (111) in the form

$${\mathit{t}}_{AB}{\mathit{t}}_b\mathit{B}-{\mathit{k}}_{BT}{T}_D-{\mathit{k}}_{BM}{\mathit{M}}_0-{\mathit{k}}_{BN}\mathit{N}={\mathbf{0}}_{\mathbf{2}}$$

(122)

Then, the system of equations for unknown lateral reactions, front wheels’ braking moments and drive torque have a form similar to Equation (118). Note that the solution of nonlinear systems can be performed using the MATLAB lsqnonlin function.

3. Optimization

3.1. Generalized Approach to Speed Model

We proceed from the fact that some trajectory for the AV along a distance D has already been built. This question was extensively reflected in our previous studies. Using the general approach, let us compose models for obtaining the speed in finite elements (FE). The distance D may be divided on n FEs, each i-th of which is represented by the length L_i and parameter ξ ∈ [0, 1]. Thus, the current linear space is x = ξL_i. The basis functions F_ξ are the form functions corresponding to the nodal degree of freedoms (DOFs) Q_v. The number of nodal DOFs depends on the degree p of Lagrangian polynomial: d = (p + 1)/2. Then, any function y(x) within i-th FE may be expressed as follows:

$${\mathit{l}}_i=\left(\begin{array}{ccc}{L}_i^0& \cdots & L_i^d\end{array}\right),\hspace{1em}{\mathit{L}}_i=diag\left(\begin{array}{cc}{\mathit{l}}_i& {\mathit{l}}_i\end{array}\right),\hspace{1em}{y}_i\left(x\right)=y_i\left({\mathit{Q}}_{vi},L_i,\xi \right)={\mathit{Q}}_{vi}^T{\mathit{L}}_i{\mathit{F}}_{\xi }$$

(123)

We will use the second speed derivative as the basic model to reduce the polynomial extent p, the quantity of nodal unknowns, and speed up computations. Then, it can be written

$$\frac{d^2{V}_{\zeta i}}{d{x}^2}={\mathit{Q}}_{vi}^T{\mathit{L}}_i{\mathit{F}}_{\xi },\hspace{1em}\frac{d^3{V}_{\zeta i}}{d{x}^3}={\mathit{Q}}_{vi}^T\left(\frac{{\mathit{L}}_i}{L_i}\right)\frac{d{\mathit{F}}_{\xi }}{d\xi }$$

(124)

where nodal and FE parameters are

$${\mathit{q}}_{v\left(i\right)}={\left(\begin{array}{cc}{\left.\frac{d^2{V}_{\zeta }}{d{x}^2}\right|}_i& {\left.\frac{d^3{V}_{\zeta }}{d{x}^3}\right|}_i\end{array}\right)}^T,\hspace{1em}{\mathit{Q}}_{vi}=\left(\begin{array}{c}{\mathit{q}}_{v\left(i\right)}\\ {\mathit{q}}_{v\left(i+1\right)}\end{array}\right)$$

(125)

Thus, d = 2 and we need the cubic Lagrangian polynomial to provide the smooth conjugating functions between all the segments. The first derivative of the longitudinal speed within i-th segment can be derived from the first integral

$$\frac{{dV}_{\zeta i}}{dx}=\int \frac{d^2{V}_{\zeta i}}{d{x}^2}dx+{\left.\frac{{dV}_{\zeta i}}{dx}\right|}_0$$

(126)

where dV_ζi/dx(0) = integration constant defined from initial conditions.

The antiderivative is determined as

$$\int \frac{d^2{V}_{\zeta i}}{d{x}^2}dx=\int {\mathit{Q}}_{vi}^T{\mathit{L}}_i{\mathit{F}}_{\xi }dx={\mathit{Q}}_{vi}^T{\mathit{L}}_i\int {\mathit{F}}_{\xi }{L}_{i}d\xi ={\mathit{Q}}_{vi}^T{\mathit{L}}_i{L}_i\int {\mathit{F}}_{\xi }d\xi$$

(127)

The longitudinal speed V_ζi in the i-th segment is found by repeating integration of Equation (126)

$$V_{\zeta i}=\iint \frac{d^2{V}_{\zeta i}}{d{x}^2}d{x}^2+{\left.\frac{{dV}_{\zeta i}}{dx}\right|}_{0}x+V_{\zeta i0}$$

(128)

where V_ζi0 = integration constant corresponding to the initial speed for i-th segment.

Components of Equation (128) are found as follows

$$\iint \frac{d^2{V}_{\zeta i}}{d{x}^2}d{x}^2={\mathit{Q}}_{vi}^T{\mathit{L}}_i{L}_i^2\iint {\mathit{F}}_{\xi }d{\xi }^2,\hspace{1em}{\left.\frac{{dV}_{\zeta i}}{dx}\right|}_{0}x={L_i\left.\frac{{dV}_{\zeta i}}{dx}\right|}_0\xi$$

(129)

3.2. Optimization Technique

Minimization of an objective function S has the quadratic form and, considering nonlinear constraints, uses the SQP method. It can be written as

$$\underset{\mathit{q}}{\mathit{min}}S\left(\mathit{q}\right)\mathrm{subject}\mathrm{to}\left\{\begin{array}{c}{\mathit{c}}_{eq}\left(\mathit{q}\right)=0\\ {\mathit{A}}_{eq}\mathit{q}={\mathit{b}}_{eq}\\ {\mathit{q}}_L\le \mathit{q}\le {\mathit{q}}_U\end{array}\right.,\hspace{1em}\mathit{q}=\left(\begin{array}{c}{\mathit{q}}_1\\ ⋮\\ {\mathit{q}}_{n+1}\end{array}\right),\hspace{1em}{\mathit{q}}_i=\left(\begin{array}{c}{q}_{d\left(i-1\right)+1}\\ ⋮\\ q_{d\left(i-1\right)+d}\end{array}\right)$$

(130)

where q = vector of nodal parameters; c_eq(q) = vector function of nonlinear equality constraints; A_eq, b_eq = matrix and vector of linear equality constraints, respectively; q_L, q_U = lower and upper limits; and i $\in$ [1, n] = segment number.

The objective function may be a sum of integrals across all FEs. The integrals can be solved numerically using an N-point Gaussian quadrature scheme. Then, any integrand z_i(x), replacing x = ξL, within an interval [x_i−1, x_i] can be evaluated as follows

$${\int }_{x_{i-1}}^{x_i}{z}_i\left(x\right)dx=L_i{\int }_0^1{z}_i\left({\mathit{Q}}_{vi},L_i,\xi \right)d\xi \approx L_i{\sum }_{k=1}^N{w}_k{z}_i\left(\xi \left({\vartheta }_k\right)\right)det\left(J\left({\vartheta }_k\right)\right)$$

(131)

where w_k = integration weight in the k-th point; ϑ_k = k-th point in the master–element coordinate system; J = Jacobian; k $\in$ [1, N]; and N = number of integration points.

For one-dimensional FE

$$\xi \left({\vartheta }_k\right)=\frac{\left({\xi }_1-{\xi }_0\right)}{2}\left({\vartheta }_k+1\right)+{\xi }_0=\frac{1}{2}\left({\vartheta }_k+1\right),\hspace{1em}det\left(J\right)=\frac{{\xi }_1-{\xi }_0}{2}=\frac{1}{2}$$

(132)

Using vector designations,

$$\mathit{\vartheta }=\left(\begin{array}{c}{\vartheta }_1\\ ⋮\\ {\vartheta }_N\end{array}\right),\hspace{1em}\mathit{w}=\left(\begin{array}{c}{w}_1\\ ⋮\\ w_N\end{array}\right),\hspace{1em}\mathit{\xi }={\left(\begin{array}{c}{\xi }_1\\ ⋮\\ {\xi }_N\end{array}\right)}^T,\hspace{1em}{\mathit{z}}_i={\left(\begin{array}{c}{z}_{i,1}\\ ⋮\\ z_{i,N}\end{array}\right)}^T$$

(133)

we obtain a short expression for calculating the integral of Equation (71) along all n segments

$${\sum }_{i=1}^n{L}_i{\int }_0^1{z}_i\left(\xi \right)d\xi \approx \frac{1}{2}{\sum }_{i=1}^n{L}_i{\mathit{z}}_i\mathit{w}=\frac{1}{2}{\mathit{L}}_s\mathit{z}\mathit{w}$$

(134)

where L_s = vector of segment lengths; z = matrix of integrand values of n × N size; and

$${\mathit{L}}_s={\left(\begin{array}{c}{L}_1\\ ⋮\\ L_n\end{array}\right)}^T,\hspace{1em}\mathit{z}=\left(\begin{array}{c}{\mathit{z}}_1\\ ⋮\\ {\mathit{z}}_n\end{array}\right)=\left(\begin{array}{ccc}{z}_{1,1}& \dots & z_{1,N}\\ ⋮& \ddots & ⋮\\ z_{n,1}& \dots & z_{n,N}\end{array}\right)$$

(135)

Here, we use N = 6-point scheme enough for quality and performance.

Now, let us choose optimization criteria.

Longitudinal Speed Deviation relative to a preset upper-level V_ζU value along the path s.

$$I_v={\int }_0^{\mathit{s}}{\left(V_{\zeta U}-V_{\zeta }\right)}^{2}ds={\int }_0^{\mathit{D}}{\left(V_{\zeta U}-V_{\zeta }\right)}^2{s}_x^{\prime }dx={\sum }_{i=1}^n{L}_i{\int }_0^1{z}_{vi}\left(\xi \right)d\xi$$

(136)

Denoting the integrand z_vi(ξ), using a set of FE speed parameters Q_vi from the model of Equation (124), a preset of the trajectory (curvature’s derivative) parameters Q_ti, and the approach of Equation (134), we have

$$\begin{gathered} z_{vi}\left(\xi \right)={\left(V_{\zeta U}-V_{\zeta }\left({\mathit{Q}}_{vi},L_i,\xi \right)\right)}^2{s}_x^{\prime }\left({\mathit{Q}}_{ti},L_i,\xi \right), \\ {{\mathit{z}}_{vi}=z}_{vi}\left({\mathit{Q}}_{ti},{{\mathit{Q}}_{vi},L}_i,\xi \left({\mathit{\vartheta }}^T\right)\right),\hspace{1em}{I}_v\approx \frac{1}{2}{\mathit{L}}_s{\mathit{z}}_v\mathit{w} \end{gathered}$$

(137)

where the trajectory scale factor is

$$s_x^{\prime }=\sqrt{1+{\left(dy/dx\right)}^2}$$

(138)

Third Derivative of Longitudinal Speed

$$I_{d3v}={\int }_0^{\mathit{s}}{\left(\frac{d^3{V}_{\zeta }}{d{x}^3}\right)}^{2}ds={\int }_0^{\mathit{D}}{\left(\frac{d^3{V}_{\zeta }}{d{x}^3}\right)}^2{s}_x^{\prime }dx={\sum }_{i=1}^n{L}_i{\int }_0^1{z}_{d3vi}\left(\xi \right)d\xi$$

(139)

Denoting the integrand z_d3vi(ξ) and considering the approach above, we obtain

$$\begin{gathered} z_{d3vi}\left(\xi \right)={\left(\frac{d^3{V}_{\zeta }}{d{x}^3}\left({\mathit{Q}}_{vi},L_i,\xi \right)\right)}^2{s}_x^{\prime }\left({\mathit{Q}}_{ti},L_i,\xi \right), \\ {{\mathit{z}}_{d3vi}=z}_{d3vi}\left({\mathit{Q}}_{ti},{{\mathit{Q}}_{vi},L}_i,\xi \left({\mathit{\vartheta }}^T\right)\right),\hspace{1em}{I}_{d3v}\approx \frac{1}{2}{\mathit{L}}_s{\mathit{z}}_{d3v}\mathit{w} \end{gathered}$$

(140)

Fourth Derivative of Longitudinal Speed

$$I_{d4v}={\int }_0^{\mathit{s}}{\left(\frac{d^4{V}_{\zeta }}{d{x}^4}\right)}^{2}ds={\int }_0^{\mathit{D}}{\left(\frac{d^4{V}_{\zeta }}{d{x}^4}\right)}^2{s}_x^{\prime }dx={\sum }_{i=1}^n{L}_i{\int }_0^1{z}_{d4vi}\left(\xi \right)d\xi$$

(141)

Denoting the integrand z_d4vi(ξ) and by analogy with the previous, we obtain

$$\begin{gathered} z_{d4vi}\left(\xi \right)={\left(\frac{d^4{V}_{\zeta }}{d{x}^4}\left({\mathit{Q}}_{vi},L_i,\xi \right)\right)}^2{s}_x^{\prime }\left({\mathit{Q}}_{ti},L_i,\xi \right), \\ {{\mathit{z}}_{d4vi}=z}_{d4vi}\left({\mathit{Q}}_{ti},{{\mathit{Q}}_{vi},L}_i,\xi \left({\mathit{\vartheta }}^T\right)\right),\hspace{1em}{I}_{d4v}\approx \frac{1}{2}{\mathit{L}}_s{\mathit{z}}_{d4v}\mathit{w} \end{gathered}$$

(142)

The speed’s objective function S_v is derived as the sum of the weighted criteria. Then, the following must be satisfied

$$S_v=S_v\left({\mathit{q}}_v\right)=W_v^T{\mathit{I}}_v\left({\mathit{q}}_v\right)\to min$$

(143)

where q_v = vector of speed’s second derivative’s nodal parameters (DOFs); I_v = vector of objective criteria integrals; and W_v = vector of weight factors.

$${\mathit{W}}_v=\left(\begin{array}{c}{W}_v\\ W_{d3v}\\ W_{d4v}\end{array}\right),\hspace{1em}{\mathit{I}}_v={\mathit{I}}_v\left({\mathit{q}}_v\right)=\left(\begin{array}{c}{I}_v\left({\mathit{q}}_t,{\mathit{q}}_v\right)\\ I_{d3v}\left({\mathit{q}}_t,{\mathit{q}}_v\right)\\ I_{d4v}\left({\mathit{q}}_t,{\mathit{q}}_v\right)\end{array}\right)$$

(144)

where W_v, W_d3v, W_d4v = weight coefficients for quadratic velocity deviations and its third and fourth derivatives, respectively; and q_t = vector of curvature’s (trajectory’s) derivative nodal parameters (DOFs).

4. Constraints

General Integral Approach to Nonlinear Equality Constraints. Since the kinematic, dynamic, and physical vehicle motion parameters have been formed, let us consider the integral technique of composing equality constraints. Suppose that smooth piecewise polynomial functions describe all the parameters based on the nodal DOFs of the speed and curvature derivatives. Since the numerical integration based on the Gaussian scheme is applied for optimization, the same scheme will be used to form restrictions. Assume that some parameter Ψ changes along the path s so that it does not exceed the upper Ψ_U and lower Ψ_L boundaries. Then, the sum of the areas between the upper limit Ψ_U and the function Ψ and between the lower limit Ψ_L and the function Ψ must be strictly equal to the area within boundaries. That is, for i-th FE and along the trajectory

$${\int }_{x_{i-1}}^{x_i}\left({\Psi }_U-{\Psi }_L\right)s_{xi}^{\prime }dx={\int }_{x_{i-1}}^{x_i}\left|{\Psi }_U-\Psi \right|s_{xi}^{\prime }dx+{\int }_{x_{i-1}}^{x_i}\left|\Psi -{\Psi }_L\right|s_{xi}^{\prime }dx$$

(145)

Integrals can be calculated numerically, for which we introduce the following denotations.

For integrals between bounds (UL), upper bound and function (UF), and function and lower bound (FL), we may denote:

$$z_{\Psi i}^{\left(UL\right)}\left(\xi \right)=\left({\Psi }_U-{\Psi }_L\right)s_{xi}^{\prime },\hspace{1em}{z}_{\Psi i}^{\left(UF\right)}\left(\xi \right)=\left|{\Psi }_U-\Psi \right|s_{xi}^{\prime },\hspace{1em}{z}_{\Psi i}^{\left(FL\right)}\left(\xi \right)=\left|\Psi -{\Psi }_L\right|s_{xi}^{\prime }$$

(146)

Then, for the i-th path segment, using Equation (131), vectors of Equation (146) functions are represented as follows

$${\mathit{z}}_{\Psi i}^{\left(k\right)}=z_{\Psi i}^{\left(k\right)}\left(\xi \left({\mathit{\vartheta }}^T\right)\right),\hspace{1em}k=UL,UF,FL$$

(147)

Combining all the segments, the integrals can be calculated numerically as

$${\mathit{I}}_{\Psi }^{\left(k\right)}\approx \frac{1}{2}{\mathit{L}}_s{\mathit{z}}_{\Psi }^{\left(k\right)}\mathit{w},\hspace{1em}k=UL,UF,FL$$

(148)

Thus, the requirement of nonlinear equality constraint along all the segments for a parameter Ψ, according to Equation (145) is expressed as follows

$$c_{\Psi }={\mathit{I}}_{\Psi }^{\left(UL\right)}-{\mathit{I}}_{\Psi }^{\left(UF\right)}-{\mathit{I}}_{\Psi }^{\left(FL\right)}=0$$

(149)

Kinematic Parameters’ Constraints. Using this scheme above, we form a vector Ψ_k of nonlinear integral constraints for a series of kinematic parameters such as longitudinal speed, yaw rate, angular acceleration, and longitudinal jerk, where each j-th element corresponds to Ψ_kj.

$${\mathit{\Psi }}_k=\left(\begin{array}{c}{V}_{\zeta }\\ \omega \\ \begin{array}{c}\epsilon \\ j_{\zeta }\end{array}\end{array}\right),\hspace{1em}{\mathit{\Psi }}_{kU}=\left(\begin{array}{c}{V}_{\zeta U}\\ {\omega }_U\\ \begin{array}{c}{\epsilon }_U\\ j_{\zeta U}\end{array}\end{array}\right),\hspace{1em}{\mathit{\Psi }}_{kL}=\left(\begin{array}{c}{V}_{\zeta L}\\ {\omega }_L\\ \begin{array}{c}{\epsilon }_L\\ j_{\zeta L}\end{array}\end{array}\right),\hspace{1em}{\mathit{c}}_k=\left(\begin{array}{c}{c}_{V\zeta }\\ c_{\omega }\\ \begin{array}{c}{c}_{\epsilon }\\ c_{j\zeta }\end{array}\end{array}\right)={\mathbf{0}}_4$$

(150)

where Ψ_kU, Ψ_kL = upper and lower limits of kinematic parameters, c_k = vector of constraints corresponding to kinematic parameters.

Dynamic Parameters’ Constraints. The vector Ψ_d of dynamic parameters includes slip angles δ, tire lateral deformations Δ_μ, roll angle ψ, and longitudinal acceleration a_ζ as a function of speed. The maximum allowable roll angle at 0.4 g should not exceed 7°. For passenger cars, the recommended limit value is ψ_U = (10.8 − 4.3 B₂₄/h_g/2)°. We may accept 7.5°. Thus, introducing the limits ψ_L = −ψ_U. Let us use symmetrical limits.

A particular case is the limitation of the vehicle’s traction potential. The vehicle’s maximum acceleration depends on the speed and is due to design features. Thus, if the vehicle speed–acceleration characteristic is known, the condition must be met

$${\int }_{V_{\zeta \left(i\right)min}}^{V_{\zeta \left(i\right)max}}\left(a_{\zeta U}-a_{\zeta L}\right)d{V}_{\zeta }={\int }_{V_{\zeta \left(i\right)min}}^{V_{\zeta \left(i\right)max}}\left|a_{\zeta U}-a_{\zeta i}\right|d{V}_{\zeta }+{\int }_{V_{\zeta \left(i\right)min}}^{V_{\zeta \left(i\right)max}}\left|a_{\zeta i}-a_{\zeta L}\right|d{V}_{\zeta }$$

(151)

where a_ζU, a_ζL = upper and lower limit values of acceleration potentially implemented by the vehicle.

The tire’s lateral deformation limits are Δ_μL = −Δ_μU. Combining the parameters, we obtain the vectors

$${\mathit{\Psi }}_d=\left(\begin{array}{c}\mathit{\delta }\\ {\mathsf{\Delta }}_{\mu }\\ \begin{array}{c}\psi \\ a_{\zeta }\end{array}\end{array}\right),\hspace{1em}{\mathit{\Psi }}_{dU}=\left(\begin{array}{c}{\mathit{\delta }}_U\\ {\mathsf{\Delta }}_{\mu U}\\ \begin{array}{c}{\psi }_U\\ a_{\zeta U}\end{array}\end{array}\right),\hspace{1em}{\mathit{\Psi }}_{dL}=\left(\begin{array}{c}{\mathit{\delta }}_L\\ {\mathsf{\Delta }}_{\mu L}\\ \begin{array}{c}{\psi }_L\\ a_{\zeta L}\end{array}\end{array}\right),\hspace{1em}{\mathit{c}}_d=\left(\begin{array}{c}{\mathit{c}}_{\delta }\\ {\mathit{c}}_{\mathsf{\Delta }}\\ \begin{array}{c}{c}_{\psi }\\ c_{a\zeta }\end{array}\end{array}\right)={\mathbf{0}}_{10}$$

(152)

Adhesion Constraints. Having calculated the necessary traction forces on the wheels according to Equation (54), including Equations (55) and (56), we determine the degree of using the longitudinal and lateral adhesions on each j-th wheel.

$${\phi }_{\zeta j}=R_{\zeta j}/R_{zj}<{\phi }_{max},\hspace{1em}{\phi }_{\mu j}=R_{\mu j}/R_{zj}$$

(153)

Then, the basic condition regarding the j-th wheel’s adhesion potential can be formed as follows

$${\phi }_j=\sqrt{{\phi }_{\zeta j}^2+{\phi }_{\mu j}^2}\le {\phi }_{max}$$

(154)

Using Equation (154), it is possible to impose physical restrictions on the conditions of all tire–road adhesions. Then

$${\mathit{\Psi }}_a=\left(\begin{array}{c}{\phi }_1\\ {\phi }_2\\ \begin{array}{c}{\phi }_3\\ {\phi }_4\end{array}\end{array}\right),\hspace{1em}{\mathit{\Psi }}_{aU}={\phi }_{max}\left(\begin{array}{c}1\\ 1\\ \begin{array}{c}1\\ 1\end{array}\end{array}\right),\hspace{1em}{\mathit{\Psi }}_{aL}=\left(\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right),\hspace{1em}{\mathit{c}}_a=\left(\begin{array}{c}{c}_{\phi 1}\\ c_{\phi 2}\\ \begin{array}{c}{c}_{\phi 3}\\ c_{\phi 4}\end{array}\end{array}\right)={\mathbf{0}}_4$$

(155)

Boundary Parameters’ Constraints. Another type of constraint determines the boundary conditions of kinematic parameters. Thus, one can require, for example, that the initial (0) and final (f) values of the predicted acceleration and jerk must correspond to preset constant values A_ζ0(f) and J_ζ0(f). That is,

$${\mathit{\Psi }}_b=\left(\begin{array}{c}{a}_{\zeta 0}\\ j_{\zeta 0}\\ a_{\zeta f}\\ j_{\zeta f}\end{array}\right),\hspace{1em}{\mathit{\Psi }}_{b0}=\left(\begin{array}{c}{A}_{\zeta 0}\\ J_{\zeta 0}\\ 0\\ 0\end{array}\right),\hspace{1em}{\mathit{\Psi }}_{bf}=\left(\begin{array}{c}0\\ 0\\ A_{\zeta f}\\ J_{\zeta f}\end{array}\right),\hspace{1em}{\mathit{c}}_b={\mathit{\Psi }}_b-{\mathit{\Psi }}_{b0}-{\mathit{\Psi }}_{bf}={\mathbf{0}}_4$$

(156)

The desired final speed may also be included.

Total restricting conditions. Thus, the complete sets of parameters, bounds, and nonlinear constraints are

$$\mathit{\Psi }=\left(\begin{array}{c}{\mathit{\Psi }}_k\\ {\mathit{\Psi }}_d\\ {\mathit{\Psi }}_a\\ {\mathit{\Psi }}_b\end{array}\right),\hspace{1em}{\mathit{\Psi }}_U=\left(\begin{array}{c}{\mathit{\Psi }}_{kU}\\ {\mathit{\Psi }}_{dU}\\ {\mathit{\Psi }}_{aU}\\ {\mathit{\Psi }}_{b0}\end{array}\right),\hspace{1em}{\mathit{\psi }}_L=\left(\begin{array}{c}{\mathit{\Psi }}_{kL}\\ {\mathit{\Psi }}_{dL}\\ {\mathit{\Psi }}_{aL}\\ {\mathit{\Psi }}_{bf}\end{array}\right),\hspace{1em}{\mathit{c}}_{eq}=\left(\begin{array}{c}{\mathit{c}}_k\\ {\mathit{c}}_d\\ {\mathit{c}}_a\\ {\mathit{c}}_b\end{array}\right)={\mathbf{0}}_{22}$$

(157)

5. Simulation

Trajectory. As an example, consider planning a speed mode on a curvilinear section shown in Figure 7. Let us assume that the maneuver trajectory has already been determined by the inverse method described in our previous studies.

Before optimizing the speed mode, let us set the values of weight coefficients W_v in Equation (144) and limit values.

$${\mathit{W}}_v={\left(\begin{array}{ccc}1& 3& 1\end{array}\right)}^T$$

(158)

Speed derivative model limits

$${\left.\frac{d{V}_{\zeta }}{dx}\right|}_U=5.7\left[\frac{1}{s}\right],\hspace{1em}{\left.\frac{d^2{V}_{\zeta }}{d{x}^2}\right|}_U=1\left[\frac{1}{ms}\right],\hspace{1em}{\left.\frac{d^3{V}_{\zeta }}{d{x}^3}\right|}_U=0.1\left[\frac{1}{m^{2}s}\right]$$

(159)

Initial and limiting values of kinematic parameters

$$\begin{gathered} V_{\zeta 0}=60\left[\frac{km}{h}\right],\hspace{1em}{V}_{\zeta U}=100\left[\frac{km}{h}\right],\hspace{1em}{V}_{\zeta L}=50\left[\frac{km}{h}\right], \\ {a}_{\zeta U}={\phi }_{max}g\left[\frac{m}{s^2}\right],\hspace{1em}{a}_{\zeta L}=-0.5\left[\frac{m}{s^2}\right],\hspace{1em}{a}_{\zeta 0}=a_{\zeta f}=0\left[\frac{m}{s^2}\right],\hspace{1em}{\omega }_U=0.5\left[\frac{rad}{s}\right], \\ {\epsilon }_U=3\left[\frac{rad}{s^2}\right],\hspace{1em}{j}_{\zeta U}=5\left[\frac{m}{s^3}\right],\hspace{1em}{j}_{\zeta L}=-2.5\left[\frac{m}{s^3}\right],\hspace{1em}{j}_{\zeta 0}=j_{\zeta f}=0\left[\frac{m}{s^3}\right] \\ {\delta }_U=12\left[°\right],\hspace{1em}{\mathsf{\Delta }}_{\mu U}=2.5\left[cm\right],\hspace{1em}{\psi }_U=7.5\left[°\right] \end{gathered}$$

(160)

We use the data of Audi A4 3.2 FSI [24] to represent the AV. All the calculations are accomplished by using MATLAB tools [25].

Comparison of drive cases’ outputs. Figure 8 shows the results of output parameters for the various torque vectoring schemes corresponding to the variants (A, B, C, D) in Figure 4. Note that the same initial conditions and vehicle data are applied for all variants. Let us consider a case of road conditions that provide the maximum adhesion coefficient φ_max = 0.5.

Longitudinal speed and acceleration. As seen in Figure 8a,b, all variants provide a smooth and similar shape increase in speed. The lowest peak values of speed and acceleration belong to variant (B) with single rear sport DM. Variants (A) and (D) show close results, which are explained by the same strategy of using the longitudinal adhesion potential for the same axle wheels. In this case, it does not matter in which way the redistribution of drive torques on the same axle’s wheels is achieved—due to internal or external moments relative to the semi-axles of the inter-wheel differential. Finally, the greatest speed and acceleration are achieved with variant (C), combining the rear sport DM and front brakes. This is because at the moment of vehicle acceleration, the inner front wheel is braked, which leads to an increase in resistance on its semi-axle with a subsequent torque raising on the opposite outer semi-axle. This results in growing the torque on the inter-axle differential’s output shafts. That is, this variant provides not only the torque distribution over the same axle wheels but partially over the axles too.

Jerk and roll angle. As seen in Figure 8c, longitudinal jerks differ in amplitude but do not exceed the specified limits, while simultaneously meeting the requirements of zero initial and final values. Obviously, with an increase in acceleration amplitude, the jerk also grows; nevertheless, the curves’ unambiguity and smoothness are ensured for all drive variants. This indicates a certain positive effect of the torque vectoring on the smoothness and stability of the curvilinear motion.

Figure 8d shows pseudo-roll angles reflecting the potential vehicle’s body roll due to predominantly centrifugal forces. Since the trajectory curvature determines the lateral acceleration, which is the same for all variants, the roll angles differ little and do not exceed the critical values.

Yaw rate and angular acceleration. As seen in Figure 8e,f, the curves of yaw rates and angular accelerations are almost the same for different drive control options. This is because the trajectory for all variants is the same and only the motion speed character stipulates the differences. Owing to the unprecedented smoothness and continuity of the angular acceleration functions, it is possible to successfully apply the inverse modeling approach, i.e., to evaluate force and dynamic parameters using the predefined kinematic ones. Note that the smoothness of all reference functions is one of the most important factors for the AV’s stable control.

Inter-axle (central) DM torques. Figure 8g shows the torques on the central differential’s carrier. As expected, the highest power costs correspond to variant (A), since the drive torques on the inner wheels are limited by the braking moments. That is, despite the high-quality regulation, this option requires excessive energy costs. The lowest energy consumption is provided by variant (B) with a rear sport DM, which requires almost 250 Nm less than in the case of variant (A). Variant (C) also ensures less torque than (A) and, nevertheless, more speed. Variant (D) is characterized by the lowest specific power costs regarding the speed mode since the torque redistribution is provided by internal and not external moments relative to the wheels’ driving semi-axles.

Figure 8h shows the distribution of torques between the vehicle axles by the variants. Note that the main difference in the values of the axles’ drive torques is observed in the first phase of acceleration since the effect of redistributing the vertical reactions is especially sensitive for the rear axle. As the steering angles and acceleration begin to diminish when the curvature decreases (about 35 m of the path), the rear tires’ kinematic radii instantly become smaller than the same radii of the front wheels. As a result, this leads to the temporal effect of equalizing the axles’ drive torques. Subsequently, owing to decreasing the influence of redistributing the vertical reactions and approaching the trajectory to a straight line, the ratio of torques becomes close to the mechanically preset in the inter-axle differential (i.e., 40/60).

Case A. Figure 9 summarizes the results of optimizing the parameters for the variant by controlling the rear and front brakes. As seen (Figure 9a,b), the sideslip angles and the tires’ lateral deformations change smoothly and do not exceed the preset limits. The peak sideslip angles are reached by the front inner and rear outer wheels. In this case, the front inner wheel moves along the path with the greater curvature, and the rear outer wheel has the maximum vertical load (Figure 9c). The distribution of lateral forces (Figure 9d) is mostly proportional to the wheel loads. Due to individual adjustment, all longitudinal forces (Figure 9e) on the wheels are different. The braking moments are of greatest interest (Figure 9f).

As seen, the braking moments of the inner wheels change smoothly, and the front inner wheel’s braking moment is greater than that of the rear one. This is due to the need for withdrawing a larger moment from the front inner wheel, which is in worse adhesion conditions. At the same time, the braking moments of the outer wheels remain zero. Owing to the proposed strategy for distributing torques between the same axle wheels, the longitudinal adhesion potentials are pairwise the same for the wheels of the front and rear axles, respectively (Figure 9g). The lateral adhesion factor (Figure 9h) and total adhesion (Figure 9i) are represented by relatively narrow bunches of the common tendency curves. At the same time, none of the wheels violates the adhesion limit requirements. Figure 9j depicts the relation between the acceleration’s longitudinal and transversal components, showing continuity, smoothness, and good agreement with the trajectory curvature and the longitudinal speed’s derivatives as well.

Case B. Figure 10 shows the results of optimizing parameters for the variant with a controlled rear sport DM. As seen (Figure 10a,b), the sideslip angles and the tires’ lateral deformations have the smallest values among all the variants since the speed is also the lowest. In Figure 10e reflecting the longitudinal reactions, as noted, the traction forces on the front wheels are almost coincident, which is explained by the equality of the drive torques on the semi-axles of the symmetric (open) front DM. Figure 10f depicts the internal moment of the rear sport DM, which is positive for the outer wheel’s clutch and zero for the inner wheel’s clutch. In Figure 10g, the rear wheels’ adhesion potentials are the same, unlike the front wheels. Thus, the adjustment possibilities are limited by the adhesion of the inner front wheel. Given the drive design at the specified road and initial conditions, this variant partially underuses the adhesion limit but provides good smoothness of the vehicle dynamics.

Case C. Figure 11 reflects the results of optimizing the parameters for the combined control variant with the rear sport DM and front brakes. As seen (Figure 11a,b), the sideslip angles and tires’ lateral deformations reach the largest values among all options since the speed is also the highest. In Figure 11e the wheels’ traction forces differ significantly in the maneuver’s first phase and converge in the second phase after decreasing the curvature. Figure 11f depicts the moments required for the rear sport DM’s clutches and front brakes. These moments are positive for the rear outer wheel’s clutch and front inner wheel’s brake, and they are zero for the rear inner wheel’s clutch and front outer wheel’s brake. In Figure 11g, the adhesion potentials of the front and rear wheels are pairwise coincident, which emphasizes the strategy’s effectiveness. Thus, the control possibilities are expanded and make it possible to bring the wheels to modes close to the limiting adhesion conditions if needed.

Case D. Figure 12 shows the results of optimizing the parameters for the variant with two controlled inter-wheel sport DMs. As seen (Figure 12a,b), the sideslip angles and tires’ lateral deformations are characterized by smoothness and their values do not exceed the corresponding limits. In Figure 12e the front wheels’ traction forces for this variant play the greatest role in comparison with other schemes, which has a positive effect on the vehicle steerability. Figure 12f depicts the internal moments of the sport DMs’ clutches, which are positive for the outer wheels’ clutches and zero for the inner wheels’ clutches. In Figure 12g, the front and rear wheels’ adhesion potentials are pairwise coincident. Thus, the possibilities of regulating motion stability and vehicle controllability are maximized.

Note that this drive variant is characterized by the highest cost despite the high efficiency and economy since there is no external resistance applied to the wheels’ semi-axle.

6. Conclusions

The main goal of this study is to improve the inverse modeling techniques for tasks on planning AV’s speed modes, considering the design features of the drivetrain structure. Based on the results, the following may be highlighted:

• The proposed technique is highly efficient. The inverse approach provides stable, unambiguous functions for all kinematic, dynamic, and physical parameters, which are characterized by continuity and smoothness. In addition, the method’s performance is distinguished by its relative simplicity and versatility, confirmed by the successfully applied technique for predicting both trajectory and speed.
• The vehicle’s motion dynamics significantly influences the formation of motion planning reference curves. This is facilitated by various dynamic and physical constraints built on the kinematic parameters to be optimized. The number of restrictions positively, in general, affects forecasts’ realism, considering the design features of vehicle’s drive, as well as the speed mode’s safety within the limits of vehicle capabilities.
• Pseudo-parameters, such as sideslip, lateral tire deformation, and body roll, were used as extended criteria for motion safety, which allowed these physical parameters to be included in the nonlinear constraints of the AV kinematic model. Even though these parameters express equivalent values, their use increases the accuracy of forecasts and reduces the probability of critical regimes.
• Parameters such as sideslip and tire lateral elasticity are the most sensitive to the curvilinear motion mode. The graphs show these parameters do not exceed the preset limits. However, the tangent angles of these curves at points close to the limit levels may be critical for the subsequent curves’ inflection at the next planning cycle. In this regard, the use of restricted derivatives for slip angles and tire lateral deformation may also be expedient.
• The considered schemes of transmissions and control for distributing torques provide the effective capability of planning the actuators’ impact based on the AV kinematic model. At the same time, in accordance with the strategy of equal longitudinal adhesion potential of the same axle wheels, the curves’ smoothness of braking moments and friction clutches’ torques of sport DMs is ensured, which has a positive effect on maintaining the stability of both control and motion. The predetermined continuous kinematic parameters, such as yaw rate and angular acceleration make it possible to calculate with sufficient accuracy and smoothness a set of force parameters required for the vehicle model’s dynamic equilibrium.
• Using a nonlinear tire model in the proposed approach had advantages and disadvantages. On the one hand, this leads to more accurate predictions of lateral forces, and on the other hand, due to the inevitable use of iterations in solving nonlinear equations, the optimization performance is worsened. Thus, regarding further improving the approach, a compromise between the modeling quality and the forecasting rapidity is needed, which may require the use of linearization followed by transiting from direct parameters to their increments with subsequent accumulation.