Computationally Efficient Nonlinear One- and Two-Track Models for Multitrailer Road Vehicles

This paper presents nonlinear mathematical models of one- and two-track multitrailer vehicles. We derive nonlinear equations of motion in the form of a system of implicit ordinary differential equations (ODEs) by using Lagrangian mechanics. The system of ODEs has the minimum number of states and equations that enables efficient computations yet maintains the most important nonlinear vehicle dynamic behavior and allows actuator coordination and energy consumption evaluation. As examples, we build different models of a 4-unit long combination vehicle, i.e., two-track 11-axle and single-track 6-axle nonlinear models as well as a linear single-track 6-axle model. We compare the performance of these models to experimental data of different driving maneuvers. The nonlinear single-track model demonstrates close dynamic behavior to the experiment, which makes it an efficient alternative to the two-track model. The vehicle equations can be generated automatically by using the code provided in this paper and subsequently used for conducting frequency analysis, evaluating energy consumption, deriving performance measures from simulations, and facilitating optimal control applications that involve combined steering, braking and propulsion control.


I. INTRODUCTION
Approximately 1.35 million road traffic fatalities worldwide were reported in 2016 [1], and the trend continues to increase. In Europe, heavy vehicles were associated with approximately 10% of road accident fatalities in 2016 [2], [3]. In recent years, advances in vehicle safety, including antilock braking systems (ABS), electronic braking systems (EBS), automatic traction control systems (TCS), and electronic stability control (ESC), have been used in commercial vehicles and have contributed to a reduction in fatalities related to these vehicles of approximately 40% between 2007 and 2016 in Europe [3]. Moreover, driving automation systems [4] show potential to further increase road safety [5]- [8], e.g., by providing adaptive cruise control, collision avoidance, automatic lane centering, and/or by performing a complete dynamic driving task. However, compared to those The associate editor coordinating the review of this manuscript and approving it for publication was Bohui Wang . regarding passenger cars, there are few publications related to active safety and driving automation of heavy vehicles, despite differences in their dynamic behavior due to differences in dimension, height of center of gravity (COG), and articulation [9].
In addition, growing demands of road freight transport [10] entail deployment of longer and heavier vehicles, i.e., long combination vehicles (LCVs) on roads [11]. Compared to tractors and semitrailers, LCVs occupy less space on roads for the same amount of transported freight and exhibit reduced fuel consumption and total cost of ownership on average of approximately 17% and 30%, respectively [12], [13]. However, there are only a few countries in the world that allow LCVs on public roads, as the dynamic behavior of LCVs must be better controlled and investigated to ensure safe performance on roads [14]- [16].
One reason why LCVs and articulated vehicles have drawn less attention in the literature than passenger cars is the difficulty in the derivation of dynamic equations of motion.
The motion of a vehicle in a specific maneuver, driving scenario, and road condition can be described by a mathematical model of vehicle dynamics, as an alternative to performing real-world experiments [17]- [23]. The fidelity and complexity of mathematical models vary depending on the vehicle model application and the trade-offs between model accuracy and computational efficiency. For example, with the objective of driving automation, including combined steering, braking and propulsion control [24]- [26], certain assumptions about system degrees of freedom (DOFs) and linearization are needed to address the most important vehicle behavior within the available computation time. Vehicle models are usually built from scratch, requiring much time and effort, and depending on the model complexity, in the case of articulated vehicles, the use of advanced simulation tools might be inevitable. Advanced simulation tools, however, require expert knowledge and detailed parameterization and incur a high computational burden. Moreover, these tools do not provide mathematical models in the form of differential algebraic equations (DAE) or ordinary differential equations (ODE) that are transferable between different solvers, e.g., with the purpose of active safety and closed-loop vehicle motion control, unless by requiring a commercial license. An example is the functional mock-up interface (FMI) toolbox [27].
Depending on the required dynamic responses of interest and with the purpose of state prediction, researchers derive their own mathematical models with the necessary level of complexity. For control design and optimization, it is often important to make the models as simple as possible, often in terms of linearity [28]- [39]. The problem then is that models in vehicle dynamics are often linearized too early, already in symbolic form, by assuming a linear tire model, as well as small articulation, steering and slip angles, and by decoupling longitudinal and lateral equations. In that case, nonlinear terms that describe essential dynamic behavior might be neglected, and the boundaries of model validity become too narrow, e.g., enforcing overly short prediction horizons in model predictive control applications. Examples of dynamic behaviors that can be lost during such a linearization are the relation between longitudinal and lateral dynamics, the phenomena from non-Ackermann geometry (between the left and right wheels or between axles), the accuracy in the lateral acceleration and articulation angle estimation during driving on a curved road and the accuracy in open-loop simulations given the input trajectories.
For vehicle on-road applications, where the vehicle body motion in the yaw plane is of interest, e.g., vehicle motion optimal control [24], [26], [40], [40], [41], and evaluations of the performance measures [42], [43], we provide nonlinear DAE and/or ODE of the dynamical system in the form of symbolic equations that can be directly used in any solver/integrator and in any programming language, e.g., in ACADO [44] and the CasADi [45] toolkit. The system dynamic equations can, therefore, be linearized either around the operating point [24]- [26], [40] with the purpose of closed-loop control or around a reference trajectory [46] with the purpose of model predictive control in long prediction horizons and real-time iteration.
An example of a real-world application, where our approach of vehicle modeling is essential, is an LCV with an electric dolly converter [47], where the optimal control of the electric propulsion by using a long prediction horizon is necessary for energy minimization purposes. In that problem, the coupling between lateral and longitudinal dynamics is necessary, together with the knowledge of actuation limits of the electric propulsion in the whole prediction horizon; i.e., electric actuator action must be planed optimally over the horizon, considering the battery state of charge and the tire and lateral stability limits. Our model has enough fidelity to evaluate the electric dolly for different operating cycles and includes dynamical limits of the dolly propulsion. Moreover, the tire forces generated by our model can act as inputs to a powertrain model for evaluating the energy consumption.
In the literature, the widely used linearized models for articulated vehicles are single-track linear (STL) models. Usually, the derived models are not compared with high-fidelity models or experimental data. Recent publications [48], [49] compared a couple of different STL models with high-fidelity models, and in [50], a nonlinear vehicle model of constant speed was compared with test data. In this paper, however, we compared the derived nonlinear equations of two-track and single-track models with experimental data including constant and varying speed maneuvers, whereas longitudinal and lateral dynamics were coupled. Moreover, according to [49], the STL model derived in [31] was one of the models that demonstrated the least error compared to a high-fidelity vehicle model. Therefore, we also compared our model with a similar STL model. The experimental data consisted of different maneuvers, including multiple single-lanechange, multiple double-lane-change, steady-state cornering, brake-in-curve, random steering, and step-steering.
Vehicle models in this paper include single-and two-track multiunit vehicles, e.g., a passenger car, a multiaxle bus or a truck, and vehicles with several multiaxle units, for example, a tractor-semitrailer and an A-double, i.e., a tractorsemitrailer-dolly-semitrailer. We did not model roll, suspension dynamics or body frame flexibility in this paper, thereby reducing the system dynamics to the two-dimensional yaw plane, which is mostly suitable for studying the motion of articulated vehicles in the normal driving condition and low lateral acceleration and/or for vehicles with a low COG height. However, we considered the effect of load transfer of vertical loads between vehicle sides and between fore and aft due to lateral and longitudinal accelerations.
We used Lagrangian mechanics to describe the system equations of motion. The benefit of using the Lagrange equations rather than Newton equations is that the constraint equations are not needed and that the system ODEs can be directly obtained in the minimum order, which can then be integrated efficiently. This means that the number of obtained equations equals the number of DOFs of the system for the fewest number of chosen independent generalized coordinates. In this case, inclusion of a vehicle unit adds only a single DOF to the system. A nonlinear multiunit vehicle model based on Newton's equations that has four DOFs for each vehicle unit, i.e., three DOFs in the yaw plane and one roll angle, can be found in [51], [52]. Generalized coordinates are any set of variables that describes the system position in the space that is consistent with constraint relations. This property of the Lagrange equations is particularly useful for modeling multibody systems comprising many bodies and force elements. The drawbacks are, however, that finding the best set of generalized coordinates is not straightforward and that the reaction forces cannot be directly calculated in the case of a minimum-order system and the absence of the Lagrange multipliers. Therefore, their calculation must be performed separately by postprocessing the ODE solution, if needed, for example, to constrain their value in an optimal control problem, which we provided in this paper as an Appendix.
In short, our main contributions comprise the following: providing system of nonlinear DAEs and/or ODEs of the dynamical system of one-and two-track multitrailer vehicles in the form of symbolic equations that can be directly used in any solver/integrator and in any programming language, whereas the derived nonlinear system of ODEs has the minimum number of states and equations that enables efficient computations yet maintains the most important nonlinear vehicle dynamic behavior; comparison of the derived nonlinear vehicle models with experimental data including constant and varying speed maneuvers, whereas longitudinal and lateral dynamics were coupled; performance comparison of an STL model with the nonlinear model and experimental data and providing a generic code that constructs the equations based on the user preferences.
The provided vehicle models can be used in various applications related to vehicles with human drivers, driving automation systems, and automated guided vehicles including evaluation of the energy consumption, efficient evaluation of performance-based characteristics, e.g., off-tracking and rearward amplification in driving cycles with different topographies and curvatures as well as motion planning and control involving combined steering, braking and propulsion control, including optimal path following at low and high speeds, trajectory planning, control allocations for maintaining stability, reverse-motion planning, steering by braking, off-tracking minimization in LCVs, and frequency analysis. We appended a MATLAB code to build the nonlinear equations of this paper as supplementary material; however, simulation examples can be found in [53].
This paper is organized as follows: Section II presents a summary of the vehicle model and its definition. Section III provides the steps of vehicle modeling. Section IV outlines the tire models and forces acting on the vehicle body. Section V describes the experimental setup, driving maneuvers, and parameter tuning. Section VI presents the approach for model linearization. Section VII provides the simulation results and compares the performance of different models with the experimental data, followed by the discussion and conclusion. Finally, the Appendix provides models for longitudinal and lateral load transfer, together with the test-vehicle data.

II. SUMMARY OF THE VEHICLE MODELING
Based on Lagrangian dynamics and the user-defined vehicle types and parameters, we create a mathematical symbolic model in the form of a nonlinear system of DAE (or implicit ODE) in the time domain: or, in space domain where x, u, and y denote the states, inputs, and variables used in algebraic equations, respectively; t denotes the time; and s denotes the space, or the distance traveled. Algebraic equations appear if, for example, longitudinal load transfer is included in the model. In this paper, similar to [31], [39], [54], the set of generalized coordinates q that describes a single-unit vehicle position in the yaw plane includes the global coordinates of the unit COG, X 1 and Y 1 , and its yaw angle φ 1 in an inertia global frame; i.e., The inclusion of the second unit, e.g., a trailer, adds a single new generalized coordinate sufficient for describing the system position in the space, that is, the articulation angle θ 1 , as shown in Fig. 1. For a vehicle consisting of four units, the set of generalized coordinates becomes where θ i denotes the angle between units i and i + 1 relative to unit i + 1, i.e., θ i = φ i − φ i+1 , according to ISO-8855 [55], where φ i is the yaw angle of the i th unit in the inertia global frame.
The system states comprise the generalized coordinates and their first derivatives; e.g., in the time domain x = [q,q], the terms up to the second derivative of the generalized coordinates appear in the equations. The input variables are the longitudinal wheel forces, i.e., braking and propulsion forces, and the wheel steering angles. A single axle or all axles can be selected to be steerable, where, for a two-track vehicle model, each axle has two steering angles for left and right wheels. An axle-end can have more than one tire that helps the vehicle carry heavier loads. Dual tires, however, receive only one longitudinal force as an input, and if they are selected to be steerable, which is not common, they receive a single steering angle input. For example, for the single-track A-double shown in Fig. 2, the input vector u is where F xwij denotes the longitudinal force on the j th lumped axle of the i th unit, caused by braking or propulsion in the  One-track representation of an A-double together with dimensions and forces acting on the vehicle (coupling forces and moments are excluded). The forces acting on the vehicle include axle forces F xw caused by propulsion/braking actuation, axle rolling resistance F wRR , air resistance F air , and body forces acting on COG of vehicle units caused by gravity and road grades. wheel hub coordinate system, and δ ij denotes the steering angle of the j th lumped axle of the i th unit.
The tire model can be selected to be either linear or nonlinear, either with or without a friction ellipse combined slip model [56], [57]. The nonlinear tire model used in this paper is based on the Pacejka tire model [57]. We tuned both the linear and nonlinear tire models to fit the real-world experimental data by solving optimization problems. The simulation results of the vehicle model by using the tuned tire models are in a very good agreement with the experimental results we evaluated in terms of the root mean squared error (RMSE) and correlation coefficient [58] of yaw rates and lateral accelerations of different units in different maneuvers. The test vehicle was an A-double heavy vehicle combination consisting of four units shown in Fig. 1. VOLUME 8, 2020 To reduce the computational demand for evaluating the performance measures or for an optimal motion control application, e.g., generating a reference for control allocation, a single-track model can be used instead of a two-track model, maintaining most of the dynamic behavior of the original model. The model and code provided in this paper build the equation of motion of a single-track articulated vehicle. An example of a single-track A-double heavy vehicle combination is depicted in Fig. 2, where an axle group is replaced with a single lumped axle. The single-track vehicle model also exhibits very similar dynamic behavior to the tested real-world vehicle.
Force elements that act on the vehicle units comprise the wheel local longitudinal force caused by actuation, the wheel local lateral force, the air resistance opposite to the local x-direction of the first unit, the rolling resistance on each wheel in the local wheel longitudinal direction, and the body forces on the COG of each unit as the result of road pitch and banking. The force elements are shown in Figs. 1 and 2.
The vehicle operational design domain can be any navigable hilly path with curvature. Even though the generated nonlinear symbolic equations are long, the simulation time is approximately 100-500 times faster than real-time driving, i.e., less than 1 s on a 5 km hilly curved road for a four-unit vehicle simulated on a laptop with an Intel Core i7 2.3 GHz processor in MATLAB. The simulation time varies depending on the driving cycle and maneuvers.

A. MODEL LIMITATIONS
The limitations of the vehicle model comprise the following: -The model is limited to yaw-plane motion since roll, wheel lift, and chassis frame flexibility are not modeled, which makes the model invalid for roll-over accident analysis. -Articulation joints are assumed to be moment-free in the Z-direction and to have no friction. -The powertrain is not modeled, and the inputs are wheel or axle forces on the ground plane rather than engine or electric motor torques. -Kinetic energy in wheel rotation is not considered; i.e., wheel rotational inertia is assumed to be zero. -It is assumed that there exists a longitudinal wheel slip controller that provides the requested longitudinal force by choosing a proper longitudinal slip, up to the grip limit of the road, i.e., up to the maximum available longitudinal force. Hence, the longitudinal slip does not appear as a variable in the model. -Air resistance force acts only on the first unit at ground height in the local x-direction. No side or lift air forces were considered. -Actuator dynamics, e.g., actuator delays, are not modeled. -Vehicle longitudinal speed cannot be zero.
To overcome the limitations above, the corresponding differential and algebraic equations must be added to the system, which in turn increase the computational demand for solving the system of equations.
Moreover, the experiment used for validating the mathematical model in this paper included normal driving maneuvers with a lateral acceleration mostly less than 2.5 m/s 2 . The lateral load transfer, therefore, was not considerable, and it was not included in the vehicle model that was compared to the experimental data. However, we implemented a steady-state lateral load transfer model to capture the possible effects in scenarios different from those provided by the experimental data. The steady-state lateral load transfer model with the assumption of roll-free couplings and a longitudinal load transfer model can be found in Appendix A. The code gives users an option to include the lateral load transfer. The longitudinal load transfer, however, is not coded.

B. VEHICLE DEFINITION AND PARAMETERS 1) VEHICLE PARAMETERS
First, the number of units and axles must be defined, which is accomplished by defining a binary matrix ua. For the A-double shown in Fig. 1, Each row of the matrix represents a unit, and each column represents the existence of an axle. If ua ij = 0, then the j th axle in the i th unit does not exist. The number of columns of matrix ua corresponds to the maximum number of axles in a unit. Similarly, matrix sa defines the steerable axles. For the A-double shown in Fig. 1, If sa ij = 1, then the j th axle in the i th unit is steerable. Other vehicle parameters related to the i th unit are the unit mass m i , the unit yaw inertia around COG in the vertical Z-direction J i , the unit height of the COG h i , and the position of the coupling points relative to the COG of the unit x ci1 and x ci2 for the front and rear coupling, respectively. The parameters related to the j th axle of the i th unit are the axle position relative to the unit's COG x aij , the number of tires per axle n taij , and the axle track width t wij . Furthermore, if the lateral load transfer is included in the model, the related parameters are the axle roll center height h rij and the axle roll stiffness c rij .

2) OPERATIONAL PARAMETERS
The other set of parameters relates to the environment and driving cycle. These parameters might vary in time or distance traveled and thus can be different for different axles. They include the road pitch angle λ p (s), the road banking or roll angle λ b (s), and the road friction coefficient µ(s). Other 203858 VOLUME 8, 2020 physical parameters are g, A f , c d , ρ a , and f r , which denote the gravitational acceleration, the vehicle front area, the air drag coefficient, the air density, and the rolling resistance coefficient, respectively. The parameters related to the tire model are explained in Section IV. Table 1 provides some of the important notations used in this paper.
In the examples considered in this paper, a simple driver model for steering, similar to pure-pursuit [41], [59], is used to follow the road by looking ahead and steering toward the aim point. For such a driver to work, the global road coordinates X (s)-Y (s) should also be defined as a function of the distance traveled. We also used a simple proportional driver model for controlling the longitudinal speed. More details are provided in Section VII.

III. VEHICLE BODY MOTION MODELING
According to [60], the Lagrange equation is given by where T is the system kinetic energy, V denotes the system potential energy, n q = 3 + n u − 1 is the total number of generalized coordinates of the system comprising n u units, and a dot (˙) above a variable represents the time derivative. The generalized coordinates q are selected as where X 1 and Y 1 denote the position of the COG of the first unit in the global inertia coordinate system, φ 1 denotes the global yaw angle of the first unit, and θ i denotes the articulation angle between the i th and the (i + 1) th units. Q l denotes the generalized force given by where n f denotes the total number of all force elements acting on the system; F Xk and F Yk denote the X and Y components of the k th force element expressed in the global inertia coordinate system, respectively; and P Xk and P Yk are the X and Y positions expressed in the global inertia coordinate system where the k th force element acts, respectively. The force elements include the longitudinal and lateral tire forces, the rolling resistance forces, and the body and air resistance forces, as shown in Figs. 1 and 2. The potential energy is zero in our vehicle on-road system. Instead, we consider the gravitational force as a force element acting on the COG of each unit. The kinetic energy is given by where i is the vehicle unit index, φ i is the vehicle unit angle in global coordinates, and v Xi and v Yi are the velocity components of the COG of the i th unit. It does not matter which coordinate system the velocity is expressed in since its magnitude remains the same irrespective of the coordinate system.
To be able to calculate (10) and (11), all force elements, the positions at which they act, and the positions of the COGs VOLUME 8, 2020 and their velocities must be determined in terms of generalized coordinates, input variables and system parameters. By using the system kinematics consistent with the coupling positions, the position of the COGs of the vehicle units can be found as The vehicle unit velocities in the global system of coordinates can be calculated by taking the time derivative of the COG positions according to The global positions of the wheels, i.e., the positions where the wheel force elements act, is the main difference between two-track and single-track vehicle models. First, let M be the transformation matrix from vehicle local to global coordinates defined by The global positions of the wheels for a single-track vehicle can be calculated as (16) and in a two-track vehicle, the global position of the right wheels can be calculated as where n a is the number of axles of the unit with the most axles in the vehicle. Similarly, the left-wheel global positions are Index k denotes the wheel index of each unit ordered starting from the right wheels, from front to rear, and then the left wheels, from front to rear, in the case of a two-track vehicle.
In the same manner, the position of the air force in the both single-track and two-track vehicles can be calculated as Having calculated all forces, their locations and the velocity of the COGs, the equation of motion (8) can be written in the final form where state vector x is It should be noted that generalized coordinates X 1 and Y 1 do not appear in the final form of the equations of motion, and to obtain them, v X1 and v Y1 must be integrated. Moreover, equation (20) is an implicit ODE or a DAE; i.e., it is not solved for the state derivatives. The reason is that the explicit equations in the formẋ = f (x, u) include a substantially greater number of terms than the implicit formulation, which makes the function evaluation more expensive. The equations have algebraic parts if the longitudinal load transfer is included.
The input vector u includes the wheel forces in the x-direction in the wheel frames F xwik and the steering angle of each wheel of the steerable axles δ ik .
The symbolic differentiation can be performed by using the MATLAB symbolic toolbox, as performed in the code provided with this paper, or any other toolbox, for example, in Mathematica, Maple, or SymPy.

A. LAGRANGE EQUATIONS IN VEHICLE LOCAL COORDINATE SYSTEM
In vehicle dynamic applications, it is more interesting to express variables and equations in a coordinate system that rotates with the first vehicle unit, i.e., in the vehicle local coordinate system. The velocity expressed in such a system is the same velocity that a driver observes. Moreover, such a system makes derivation of tire forces easier, especially in vehicles with a single unit, as tire slips depend on lateral and longitudinal velocities expressed in the vehicle (or wheel) local coordinate system.
The conversion from the global coordinate system to the vehicle local coordinate system, known as quasi-coordinates, can be accomplished by the change in variables from the 203860 VOLUME 8, 2020 global velocity of the first vehicle unit to the local velocity according to where M 1 is given by (15), and v x1 and v y1 are the longitudinal and lateral velocities in the first unit's local system of coordinates, respectively. It should be noted that the change in the variable shown in (22) makes the kinetic energy a function of the local velocities. In contrast, the local velocities are functions ofẊ 1 , Therefore, chain rule differentiation must be performed to differentiate the kinetic energy in (8). Hence, whereḡ(φ 1 ) is a part of the kinetic energy that is a function of φ 1 but not a function of v x1 and v y1 . However, always ∂ḡ(φ 1 ) The new set of states in the first vehicle unit local coordinate system is

B. CONVERSION TO THE SPACE DOMAIN
In the case of performing the simulation of the vehicle model on actual driving cycles with varying grades and curvatures and/or for energy minimization purposes, it is more convenient to use the distance traveled s as an independent variable rather than the time since the position on the road and its grade and curvature are better described as functions of the distance traveled. The change in the independent variable is performed according to and where a dot (·) represents any variable. We assumed ds dt ≈ v x1 , as the resulting error in the calculation of the distance traveled by performing this assumption is negligible.
The resulting equation of motion in the space domain is in the form where the state vector x is

IV. FORCE ELEMENTS AND TIRE MODELS
There are different kinds of forces that act on the vehicle units: the tire longitudinal force F xw caused by propulsion/braking actuation in the ground plane, the tire lateral force F yw as a result of the lateral slip, the tire rolling resistance F wRR , the air resistance F air , and the gravitational body forces acting on the COGs of the vehicle units caused by gravity and road grades. The longitudinal tire forces F xw , caused by propulsion or braking, are parts of the inputs u to the model. The tire lateral forces, i.e., F yw , are, however, caused by the tire lateral slip and must be calculated by using a tire model.

A. SLIP
Tire longitudinal and lateral slips, i.e., s x and s y , are defined as where R denotes the tire radius, ω denotes the tire rotational speed, v wyik and v wxik denote the wheel hub velocity components in the wheel local coordinate system, and i and k denote the unit and wheel indices (or axle, if single-track), respectively. In addition, there exist other definitions of the slips, but those shown in (32) are derived from the brush model, which is a physically motivated derivation of a tire model. Generally, both s x and s y can be nonzero at the same time, which is an operating condition called a ''combined slip''. However, if s x = 0 and s y = 0, we call the operating condition a ''pure lateral slip'' given by The wheel velocity local components can be calculated by differentiating their global position with respect to time by using (16)- (18) and then by transforming to the wheel local coordinate system, as follows: Two different tire models, with and without the combined slip model, have been considered in this paper and in the attached code. More details are provided in the following sections. The linear tire model, i.e., the linear relation between the lateral force and the lateral slip, is given by where F ywik | (s x =0) is the tire lateral force in the wheel coordinate system at zero longitudinal slip s x , and C yij is the tire cornering stiffness. We can assume that C yij is proportional to the tire vertical force with the proportionality coefficient C cyij . We refer to (35) as a linear tire model regardless of considering the combined slip, which makes the actual tire model nonlinear.

2) NONLINEAR TIRE MODEL
The nonlinear tire model, i.e., the nonlinear relation between the lateral force and the lateral slip, is inspired by the Pacejka magic tire model [57] that also considers the vertical force acting on the tire as a nonlinear function. The tire model is developed in the project ''Performance-based standards II'' [61], for the lateral motion simulation of heavy vehicle combinations, initially reported by [62]. A modified version of that model is presented here, where the tire lateral force F yw is given by where F zik is the tire normal force, C cyik is the cornering coefficient, and C is defined as where u 2 = 0.8 is the ratio between the road friction coefficient at a large slip and the road friction coefficient at zero slip; where u yg is the maximum lateral force gradient, which is a tuning variable that ranges between −0.3 and −0.1, and F z0ik is the nominal tire normal force. The cornering coefficient C cy is defined as where C cy0 is the cornering coefficient at the nominal tire normal force. The cornering coefficient C cy0 must be estimated based on the experimental data. In this paper, its value was tuned in the range between 1 and 20 for two different types of tires, i.e., the tires of the steerable axles and the tires of the unsteerable axles.

C. COMBINED SLIP
It is important to consider the tire combined slip model in driving cycles with curvature and combined steering together with braking and propulsion. For example, if both the lateral and the longitudinal slip are approximately 0.1, then neglecting the combined slip in the tire force calculations results in overestimation of the total force by approximately 20% on a high-friction road. This result was obtained by analyzing the combined slip model of an anisotropic tire brush model with a parabolic pressure and stick and slip friction, which is one of the tire models that can be explained by a physics-based model [56].

1) BRUSH MODEL
We provided here the combined slip brush tire model as a reference for validating a simpler combined slip model, i.e., the friction ellipse. The anisotropic tire brush model with a parabolic pressure and stick and slip friction can be derived as in [56] and reads where We removed the unit and axle indices for the sake of brevity. In the above equations, C x , C y , µ slip , µ stick , L, and F z denote the longitudinal tire stiffness, the lateral tire stiffness, the road friction at nonzero slip, the road friction at zero slip, the road-tire contact length, and the tire vertical load, respectively. Equations (40) and (41) are valid if the tire is rotating and ωv x ≥ 0. Fig. 3 shows the tire forces resulting from the combined slip based on the brush model. It should be noted that the left-hand side of (40), i.e., F xw , accounts for the total tire longitudinal force, whereas we neglected the rolling resistance.

2) FRICTION ELLIPSE MODEL
In this paper, the tire longitudinal forces serve as inputs, and it is assumed that there is a longitudinal slip controller on-board that regulates the longitudinal slip to deliver the requested longitudinal force. Therefore, we did not calculate the longitudinal slip. A model that accounts for the combined slip effect in the lateral force by knowing the longitudinal force instead of the longitudinal slip is the friction circle (or ellipse) model [56], [57], based on the Pythagorean theorem, i.e., where we neglected the rolling resistance, and e is a scaling factor that defines the so-called conservativity or the shape of the combined slip model; i.e., if e = 1, then the model is a friction circle, and if e > 1, then the model is an ellipse that represents the anisotropy, with more force available longitudinally than laterally. The lateral force at zero slip, i.e., F ywik | (s x =0) , is given by either (35) or (36). The friction circle combined slip model is a conservative model, for small values of s x , with respect to the reduction in the lateral grip; i.e., the reduction in the lateral force is higher for an increase in s x than a more advanced combined slip model, such as the combined slip model of the anisotropic tire brush model with a parabolic pressure and stick and slip friction. The comparisons of the two models are shown in Figs. 4 and 5. Fig. 4 compares the two models for different values of e and illustrates the region where the friction ellipse model is not conservative. Here, we calculated F yw based on the friction ellipse model by using (42), whereas F xw is given by (40). Fig. 5 shows the magnitude of the total force caused by the longitudinal and lateral slips for the two combined slip models.
Finally, the tire forces in the global coordinate system are

D. GRAVITATIONAL FORCE, AIR AND ROLLING RESISTANCE FORCES
If the differences for the road pitch and banking angles between the two successive axles are negligible, which is the case in the code provided with this paper, then the effect of the gravitational force can be simplified considering two forces that act on the COG of the vehicle units  given by where the road pitch angle λ pi is positive downhill in front of the vehicle unit, and the road roll (or banking) angle λ bi is positive downhill at the left side of the vehicle unit. If the differences in the road angles between the two successive axles are not negligible, then the effect of the gravitational force should be considered through the vertical forces acting on the axles rather than the COGs of the vehicle units according to where λ pij and λ bij represent the road pitch and banking angles in the global positions of the axles, respectively, and F Zik is the wheel vertical load in the global system of coordinates. The calculation of F Zik is explained in Appendix A-B. The forces calculated above should be expressed in the global system of coordinates, for example, for the forces acting on the COGs, according to where φ ri denotes the road yaw angle in the global system of coordinates and at the global position of the COG of the i th unit. If it is assumed that the vehicle units remain parallel to the road, then φ ri = φ i , which is the assumption made in the code provided with this paper. The air resistance force acting at the front of the first vehicle unit in the vehicle local system of coordinates is given by and in the global coordinates is where v x1 is either a state in the case of quasi-coordinates or given by (22) if the quasi-coordinates are not used. The rolling resistance forces F wRR are defined as and in the global coordinates as Physically, the rolling resistance of a tire on a hard ground is a torque on each wheel, counter-directed to the wheels' rotational velocity. It appears very similar to a bearing toque loss. Hence, it adds up to the wheel torque, which is the sum of the propulsion torque, brake torque and rolling resistance. In this paper, however, we modeled the propulsion and brake torques as input forces. Therefore, the rolling resistance was also modeled as a force.

V. EXPERIMENTAL DATA AND TUNING THE TIRE MODELS A. EXPERIMENTAL DATA
To validate the physics-based modeling and applicability of the made assumptions, we compared the simulation result to the real-world experiment for certain driving maneuvers. The tested vehicles included an A-double as shown in Fig. 6. The experiment was conducted as a part of the projects Optimal Distributed Propulsion [63] and Performance-Based Standards II [61] in collaboration with Chalmers University of Technology, REVERE lab, and Volvo Group Truck Technology. The experimental tractor was a Volvo FH16, 6x4T. Both semitrailers were equipped with the two inertial measurement units at the front and rear of the unit. The localization was performed by using an OxTS RT3000 GNSS inertial system [64] positioned in the tractor cabin. The data were collected at 100 Hz frequency on dry asphalt. The vehicle parameters and load distribution can be found in Appendix B

B. TUNING THE TIRE PARAMETERS
The vehicle parameters in the vehicle model are the same as the measured/estimated parameters of the experimental vehicle. The tire model parameters, however, needed to be tuned. The advantage of the nonlinear tire model presented in section IV-B is that the number of parameters to tune is low, i.e., four parameters in the case of a vehicle model with the lateral load transfer, including C cy0 and u yg for two types of tires, i.e., the steering tires and the other tires, and two parameters in the case of a vehicle model without considering the lateral load transfer, i.e., C cy0 for the steering and other tires. The parameter e can also be tuned for a brakein-curve maneuver as the combined slip occurs during such a maneuver. However, in this paper, e = 1. We performed the tuning of the tire parameters by solving optimization problems in the form find C s cy0 , u s yg , C ns cy0 , u ns yg to minimize i∈ [2,4] 1 t e −t 0   yaw rate of the i th unit measured in the experiment; and t 0 and t e represent the start and end tuning time, respectively. The cost function is the sum of the root mean squared of the difference between the measured yaw rate and the modeled value in the power of three. The power of three is used to give a higher weight to larger values of the yaw rate. The tuning was performed for the yaw rate of the fourth unit (i = 4) and the second unit (i = 2) of the A-double vehicle.
The reference experimental data used for tuning the tire parameters included the multiple single-lane-change maneuver of the vehicle shown in Fig. 6, at a speed of 80 km/h. Figure 7 shows the reference steering angle inputs for the first axle of the tractor, i.e., the steering wheel angle divided by the steering wheel ratio (18.6), used during the experiment. We assumed that the left and right road wheel steering angles of the tested vehicle were equal, i.e., δ 11ref = δ 14ref . The measured road wheel steering angle served as an input to the open-loop simulation, whereas we used a simple proportional closed-loop controller to control the vehicle speed, where the brake and propulsion forces served as controlled inputs.
The four tire parameters were tuned for the two-track A-double vehicle model shown in Fig. 1 with a nonlinear tire model excluding the load transfer. We assumed that the road friction coefficient is one, i.e., µ = 1, since the experiment  was performed on dry asphalt. The tuning parameters are shown in Table 2.
We used a similar optimization problem as (53) to tune the linear tire parameters, i.e., the six-axle-group cornering stiffness, for both the vehicles shown in Figs. 1 and 2. The optimization design variables were then changed to C yg , where g represented the axle-group index, g = 1, . . . , 6. The wheel cornering stiffness C yij could then be calculated by dividing C yg by the number of wheels per axle group. The tuning was performed for the 11-axle two-track model. Table 3 shows the tuned parameters.

VI. MODEL LINEARIZATION
In vehicle dynamics, it is common to simplify the models by removing the insignificant terms, using engineering judgments regarding the application of the vehicle model, and with the purpose of model linearization. For example, model linearization can be performed either by using a linear tire model, assuming small angles and other simplifying physical assumptions [31], [48], or by removing terms that include the products of inputs and states and their derivatives [39]. The validity of such a linearization is limited to, for example, applications where the lateral dynamics are exclusively studied neglecting the brake and propulsion forces and assuming a small lateral acceleration, constant longitudinal high speed, low curvature, etc. However, if the problem studied allows such limitations, the models that are linearized by using the physical assumptions can be solved efficiently.
The inclusion of the longitudinal forces in the equations makes the model nonlinear since the products of input forces and states appear in the equations. Moreover, the longitudinal speed should not be a state or an input but a known constant parameter (or a time-varying variable), since the products of VOLUME 8, 2020 the longitudinal speed and states appear in many terms that cannot be assumed negligible.
In this paper, the single-track model shown in Fig. 2 was linearized by making the simplifying assumptions as explained above. The equation concerning the longitudinal dynamics was neglected, and instead, the longitudinal speed was considered as a parameter. Moreover, all force elements except the lateral tire forces, i.e., the gravitational, rolling and air resistance forces, were set to zero. A straightforward approach for obtaining the linear model was to use the first-order Taylor expansion around zero states and inputs, except for v x , which was assumed to be a known parameter. The obtained lateral linear vehicle equation of motion was similar to the STL model presented in [31], and it readṡ where the states are x = [φ 1 , θ 1 , θ 2 , θ 3 , v y1 ,φ 1 ,θ 1 ,θ 2 ,θ 3 ], and the inputs are u = [δ 11 , δ 31 ].

VII. RESULTS AND COMPARISONS
We tested the correct implementation of the Lagrange equations against a similar A-double vehicle model derived based on the Newtonian mechanics in Modelica. For the same given parameters, the simulation of both models resulted in the same solutions. Fig. 8 illustrates the trajectory and the yaw rate of the 4 th unit of a nonlinear 6-axle group and a single-track A-double vehicle obtained by solving the Newtonian and Lagrangian mechanics. It should be noted that we did not present the Newton equations to limit the scope of this paper. Moreover, we compared the three basic vehicle models with each other and with the experimental data for different maneuvers. The three vehicle models were • a nonlinear 11-axle and a two-track A-double, shown in Fig. 1, with the nonlinear tire model including the combined slip, without lateral load transfer; • a nonlinear 6-axle group and a single-track A-double, shown in Fig. 2, with the nonlinear tire model including the combined slip, without lateral load transfer; and • a linear 6-axle group and a single-track A-double with the linear tire model without the combined slip and lateral load transfer.  We used the same nonlinear tire parameters as given in Table 2 in both nonlinear vehicle models. The tire parameters, i.e., the axle-group cornering stiffness of the linear vehicle model, were according to Table 3.
The yaw rate and lateral acceleration for different vehicle models in the multiple single-lane-change maneuver are illustrated in Figs. 9 and 10. We used the same maneuver for tuning the tire parameters.
The RMSE of the yaw rate and lateral acceleration obtained by simulating the mathematical models compared with the experimental data are shown in Table 4 for different maneuvers. We chose unit 4 to compare the signals because it shows the highest yaw rate and lateral acceleration compared  to those of the other units. In addition, we used the correlation coefficient r (−1 ≤ r ≤ 1) as a dimensionless index to measure the correlation of the two signals [58]. The correlation coefficient is given by where y n denotes the sampled experimental data,ȳ denotes the mean of all the samples of the experimental data, f n denotes the model output, andf denotes the mean of all samples of the model output. The correlation coefficients of the yaw rate and lateral acceleration are shown in Table 5. Furthermore, the steering and speed inputs together with the yaw rate and the lateral acceleration of the brake-in-curve FIGURE 11. The road wheel steering input in both the experiment and the mathematical model. The longitudinal speed was followed by a closed-loop proportional controller. The maneuver is the brake-in-curve. maneuver are shown in Figs. 11-13 to demonstrate the performance of the nonlinear vehicle models in a combined steering, braking, and propulsion condition. In this maneuver, for the linear model, we treated the forward velocity as a known time-varying variable.
We also compared the vehicle models and the experimental data in the frequency domain by using the random steering maneuver. The input-output gain and phase difference are illustrated in Fig. 14, where the input is the road wheel steering angle, and the output is the yaw rate. The resonance frequency of unit 4, according to the experimental data, is approximately 0.5 Hz, excluding the noisy high-frequency response. However, the simulation shows a resonance frequency of approximately 0.42 Hz. The difference between the simulated and measured resonance frequency is due to the delay in the actual vehicle yaw rate response to the steering input that was not modeled. The delay could be due to the tire relaxation [48], [56]. The tire relaxation refers to a transient process between two steady-state conditions of the tire. Tire relaxation can be modeled by describing the rate of change of the lateral slip s y using a first-order nonlinear differential VOLUME 8, 2020  The input-output gain and phase difference comparison between the different vehicle models and measured experimental data. The input is the road wheel steering angle, and the output is the yaw rate of unit 4.

equation according tȯ
where l r is the relaxation length, which is typically 0.25-0.5 of the tire circumference. Note that ifṡ y = 0, then s y is the same as in (33). Consideration of the tire relaxation improved the yaw rate response in the frequency domain, as shown in Fig. 15, for l r = 0.3. The effect of considering the tire relaxation on the RMSE in the time domain is shown in Table 6.

A. SIMULATION RESULTS OF A DRIVING CYCLE WITH VARYING ROAD GRADE AND CURVATURE
All experimental maneuvers were performed on a flat road. However, it is possible to use the nonlinear vehicle models to  predict the vehicle motion and behavior in long, curved and hilly driving cycles. To be able to run closed-loop simulations on curved and hilly roads, a driver model is required to control the steering, propulsion and braking. For the propulsion and braking, we used a simple proportional controller that generates positive and negative force inputs to follow the reference speed. The positive traction force was distributed between the driven axles, whereas the negative total brake force was distributed between all axles. For the steering control, a simple driver model involves steering the wheels toward an aim point on the road, i.e., similar to the pure-pursuit method [41], [59]. The steering angle of the first unit front axle can be calculated as where X h (t) and X h (t) are the coordinates of the aim point on the road. The absolute value of the steering angle cannot be 203868 VOLUME 8, 2020   larger than 0.8 rad. The first term in (57) must be modified to give a correct angle in the case of a circular road with more than 180 degrees of rotation. The motion of the vehicle shown in Fig. 2 is simulated on a hilly circular road shown in Fig. 16, given the inputs shown in Fig. 17, and the simulation results are shown in Fig. 18.

VIII. DISCUSSION
In the majority of maneuvers, the two-track 11-axle vehicle performs better than the other vehicle models. However, it does not perform as well as the other two models in the steady-state cornering. The reasons might be non-Ackermann steering geometry, disproportional steering ratio between the steering wheel and the road wheel, side force steering, scrubbing wheels of distributed axles, etc. The nonlinear single-track vehicle model, however, performs well in all maneuvers, despite having fewer tuning parameters compared to that of the linear vehicle.
The linear vehicle would behave poorly if we were to reduce the number of tuning parameters, as it would have fewer DOFs to match the experimental data. It should also be noted that we tuned the tire model by using a single driving maneuver excluding the lateral load transfer. Potentially, the inclusion of more maneuvers, a greater number of tire model parameters to tune, and the inclusion of the lateral load transfer in the tuning procedure could result in better vehicle model performance.
Moreover, the correlation coefficients of all models and the experimental data are low in the steady-state maneuver. The reason is that the signals are ''steady'' with an offset from each other.
The linear vehicle model performs better in estimating the yaw rate than the lateral acceleration. Therefore, with the purpose of evaluating the performance measures, such as rearward amplification, we recommend using the yaw rate instead of the lateral acceleration in the case of a linear model. Moreover, there exists a high variation in the lateral acceleration data that might be due to the sensor error and/or the suspension compliance, as the IMU sensors were mounted at a height of approximately 1.4 m.
The lateral and longitudinal load transfers used in this paper are steady-state models. The reason for choosing a steady-state load transfer was to avoid the additional states in the equation of motion that have a minor influence on the vehicle behavior under normal driving conditions but demand a high computational cost. For the same reason, the longitudinal load transfer was neglected in the code, as it adds algebraic equations to the system, even in the case of the steady-state model. However, for the sake of completeness, we derived and included the longitudinal load transfer equations in Appendix A. In addition, the steady-state lateral load transfer can be disabled in the code because it has a minor effect on the vehicle motion in the case of a low lateral acceleration.
For some maneuvers, the steering input had a high rate of change that caused the DAE solver to stop before reaching the final simulation time. This issue can be avoided by changing the solver's settings, e.g., by changing the absolute and relative tolerances. However, a large tolerance reduces the accuracy of the solution and thus increases the RMSE. An alternative approach is to smooth the steering input data.
The behavior of the nonlinear single-track model with the grouped and lumped axles is similar to that of the two-track 11-axle model, even though no additional parameter tuning was performed for the single-track model. Moreover, it has considerably fewer terms in the equations, and it is approximately 3-5 times faster to simulate than the two-track model. Therefore, the nonlinear single-track models can be used instead of the two-track models in applications involving combined steering, propulsion and braking, except in the applications where there exists a different propulsion and/or braking on the left and right wheels.
The linear model also performs well in maneuvers with low lateral acceleration and small steering angle, and it is computationally more efficient, approximately 5-8 times, than the nonlinear single-track model. Moreover, a linear vehicle model can be used to simulate the lateral vehicle behavior in combined steering, braking, and propulsion conditions, as we did for the brake-in-curve maneuver. However, the longitudinal and lateral dynamics are decoupled, and the forward velocity is treated as a known parameter in the lateral model. In this case, the effects of the forward acceleration/deceleration, the nonlinear tire behavior and the combined slip are lost in the lateral motion. These effects are important in maneuvers with high acceleration and/or low road friction.
The tire parameters were obtained under a specific friction coefficient, i.e., for µ = 1, since the experiment was performed exclusively on dry asphalt. More experiments on roads with different friction coefficients are needed to verify if the tuned tire parameters are also suitable for those conditions and, consequently, to improve them if needed. However, since the road friction coefficient serves as a parameter in the nonlinear tire model, we expect that the tire behavior can still be predicted to an acceptable level using the already tuned parameters.
The consideration of tire relaxation in the model improves the yaw rate response to the steering input in the frequency domain. In the time domain, however, it does not have a considerable influence on the yaw rate and/or lateral acceleration responses unless there are high-transient changes of the steering input, such as the case of the random steering according to Table 6 compared to Table 4. In these tables, it can be seen that the RMSE of the random steering maneuver drops considerably after including the tire relaxation, whereas there is a slight change in RMSE of the other maneuvers. Moreover, the tire relaxation adds a new single state for each of the wheels or each of the axle groups. Therefore, it adds to the computation time, and it can be neglected in applications where fewer states in the ODE are preferable.
The nonlinear vehicle models showed the ability to predict the vehicle motion behavior in all tested maneuvers and on a hilly curved road. In addition, the models can predict vehicle motion at low speed, large steering and articulation angles, and high propulsion and braking, up to the tire limits, in the forward and rearward motions. However, since the roll dynamics were not modeled, the validity of the vehicle models is limited to the low lateral acceleration below 2 m/s 2 .
By ''computationally efficient models'' we mean the following: 1) efficiency regarding the physics-based modeling, compared to a high-fidelity model, whereas the proposed models could be efficiently integrated, i.e., 100-500 times faster than real-time driving, and yet they maintained the most important vehicle behavior close to the experiment; 2) efficiency regarding the mathematical modeling, since the proposed model derivation using the Lagrangian mechanics ''directly'' yielded the minimum number of equations and states, compared to that of the Newtonian mechanics, where the minimum number of equations and states could still be obtained but ''indirectly'' by performing additional manipulation of the equations and solving/removing the algebraic equations.
Finally, the code provided with this paper generates the nonlinear equations of motion for vehicles with an arbitrary number of units based on the user-defined parameters. The single-track and two-track vehicle models can be generated. More examples together with the simulation code can be found in [53].

IX. CONCLUSION
Based on Lagrangian mechanics, this paper presented computationally efficient nonlinear one-and two-track models for multitrailer road vehicles. The generated models are systems of implicit ODEs comprising the fewest states necessary to predict the vehicle motion under different normal driving conditions. We neglected the roll dynamics and made certain assumptions about the lateral and longitudinal load transfer to minimize the number of states and equations yet maintain nonlinear vehicle behavior under combined steering, braking and propulsion operating conditions. We compared the vehicle model results to the experimental data for different maneuvers. The single-track nonlinear vehicle model demonstrated similar dynamic behavior to the experiment, with low RSMEs. The nonlinear equations can be used for vehicle motion planning and control within the model predictive control framework, for evaluating the performance measures, and for the combined optimal control of the energy usage and lateral dynamics. Furthermore, we provided the MATLAB code for an automatic generation of the multitrailer vehicle models.
The future work includes using the presented vehicle models in evaluating the performance measures, designing linear and nonlinear model predictive controllers for vehicle motion planning and trajectory tracking and predictive energy-usage optimization.

APPENDIX A CALCULATION OF THE VERTICAL FORCE AND LOAD TRANSFER A. LONGITUDINAL LOAD TRANSFER
In the provided code, the load transfers in the longitudinal and lateral directions are neglected, and the static axle vertical forces must be provided as input parameters by the user. However, by using the Newtonian equations, we provided a procedure for inclusion of the load transfer for specific articulated vehicles with an arbitrary number of vehicle units, where the axles of the articulated units have equal vertical forces. Moreover, the solution includes the reaction coupling forces.
A more general vehicle model includes driving cycles where the road grades have high variations along the road. The high variation of the road grades is particularly important for modeling of LCVs, where each axle experiences a different road pitch and banking angle. Fig. 19 illustrates the free-body diagram of a vehicle on such a road.
The Newtonian force and moment dynamic equations are used for calculating the global vertical forces F Zij and the coupling forces F cXi1 , F cXi2 , F cZi1 , and F cZi2 . We made the following assumptions to model the longitudinal load transfer of the articulated vehicles: -The articulation angles are small, i.e., cos(θ i ) ≈ 1.
-The coupling moment at the articulation points is very small, i.e., M ci1 ≈ M ci2 ≈ 0. -The pitch angular acceleration of each of the vehicle units is very small. -The global vertical forces of all axles in an articulated vehicle unit are equal, e.g., for the second unit shown in Fig. 19, F Z21 = F Z22 = F Z23 . -A single-track vehicle is used to calculate the longitudinal load transfer. The distribution of the global vertical force between the left and right wheels is considered in the next section. Following the above simplifications, by using the Newtonian equations for the force equilibrium in the global X and Z directions, including the action and reaction forces at the coupling points, and considering the moment equilibrium around the COGs, it can be shown that the vertical and coupling forces can be found by solving an algebraic system of equations in the form where the unknowns are and F Zi represents the vertical force of the group axle belonging to the i th unit. Let α i be the pitch angle of the i th unit, h ci1 be the height of the front coupling point of unit i, h ci2 be the height of the rear coupling point of unit i, and n ai be the number of the axles of unit i. Matrix A can be derived by using the Newtonian equations. For the sake of brevity, the derivation steps are omitted here. Let B, C i , and D denote the submatrices defined (60) and (61), as shown at the bottom of the next page. and, If the vehicle consists of a single unit and if there exist two vehicle units and in the case of three units in the vehicle The same pattern can be seen if the vehicle has four units.
The right-hand side of (58) is defined as

B. LATERAL LOAD TRANSFER
The model of the lateral load transfer used in this paper is very similar to the model presented in [56] but is adopted for the articulated vehicles. The intended applications of this paper, however, are mainly normal driving scenarios with low lateral acceleration (below 2 m/s 2 ), so the lateral load transfer can be neglected. The model assumptions include the following: -The roll angular velocity and acceleration are assumed to be small, so the load transfer becomes steady state. -The coupling roll moments at articulation points are very small, so the lateral load transfer does not transfer between units. However, in truck-full trailer combinations, this assumption is realistic [31]. -The axles of all articulated units are grouped to a single axle. Let l r = − 1 n ai −1 n ai j=2 x a1j be the distance from the COG to the rear axle group of the first unit, L = x a11 + l r the wheelbase of the first vehicle unit, h rg1 = 1 n a1 −1 n a1 j=2 h r1j the roll center height of the rear axle group of the first unit, c rg1 = n a1 j=2 c r1j the rear axle-group roll stiffness, and F lfr , F lfl , F lfl , and F lfl the lateral load transfers of the front-right, front-left, rear-right, and rear-left wheels, respectively. Let us also define dh = h 1 − (l r h r11 − x a11 h rg1 )/L. The steady-state lateral load transfer of the first unit [56] can then be defined as wherev y1 is the unit COG's local lateral acceleration, i.e., in the unit's local system of coordinates. A similar approach should be applied for the other articulated units if they are full trailers. Otherwise, all axles of the articulated units should be grouped to a single axle. Let F lri and F lli be the lateral load transfer of the right and left sides of the axle group of unit i, respectively. The lateral load transfer is then defined as where h rgi = 1 n ai n ai j=1 h rij is the axle-group roll center height, and c rgi = n ai j=1 c rij is the axle-group roll stiffness. In the above equations, the track width of the grouped axles is assumed to be equal.
In addition to the load transfer caused by the lateral acceleration, the lateral load transfer due to the road banking should be calculated; see Fig. 20. It is assumed that the road banking angle is the same for the left and right wheels of the same axle. Therefore, by using (69) and (70) and the solution of (58) for the total axle-group forces, i.e., F Z11 and F Z12 for the first vehicle unit and F Zi for the other units, the left-and right-side global group-axle forces can be defined as F Zir = F Zi w tijl w tij + F lri , F Zil = F Zi w tijr w tij + F lli , i = 2, . . . , n u , j = 1, where the first term accounts for the load transfer due to the road banking, and w tijr = w tij cos(λ bij ) 2 The calculated global vertical forces on right and left sides of the group axles should be ungrouped, and their compo-  nents in the road local coordinate system should be found, at the given road position. The z-direction of the road local coordinate system is perpendicular to the road surface, and its y-direction is perpendicular to the road longitudinal direction with an angle given by the banking angle in the global system of coordinates. Following the wheel indexing used, e.g., in (18), the ungrouping of the group-axle global vertical forces to the tire forces F Zik is performed as , j = 2, . . . , n a , n aij = 0, k = j (78) , j = 2, . . . , n a , n aij = 0, k = j + n a , and for the other vehicle units where n taij denotes the number of tires per axle. The forces F Zik are then used in (46) and (47) to find their projections on the road surface in the road local x-and ydirections.

APPENDIX B VEHICLE PARAMETERS
The parameters of the tested vehicle are provided in Tables 7 and 8.

APPENDIX C SUMMARY OF THE PROVIDED CODE
We provided a MATLAB file as supplementary material with this paper. The file ''generateEqs.m'' generates the equations of motion of the user-defined vehicle in the form of implicit ODEs and stores them in ''Equations.m,'' which is created automatically.
The equations can be integrated by using a DAE solver in MATLAB. Alternatively, the ODEs can be used in a different environment or programming language with the purpose of simulation or optimal control. The users have options to define a vehicle comprising an arbitrary number of units and axles. The vehicle can be either single-track or two-track. The tire model can be either linear or nonlinear. Moreover, the users have options to include the lateral load transfer and combined slip models. See [53] for simulation examples and vehicle definitions with more user-defined options, as well as the model linearization.