Open AccessArticle

Influence of the Road Model on the Optimal Maneuver of a Racing Motorcycle

Jan Biniewicz

and

Mariusz Pyrz

Faculty of Automotive and Construction Machinery Engineering, Institute of Vehicles and Construction Machinery Engineering, Warsaw University of Technology, 02-524 Warsaw, Poland

Author to whom correspondence should be addressed.

Appl. Sci. 2024, 14(10), 4006; https://doi.org/10.3390/app14104006

Submission received: 12 April 2024 / Revised: 2 May 2024 / Accepted: 5 May 2024 / Published: 8 May 2024

(This article belongs to the Section Transportation and Future Mobility)

Download

Browse Figures

Versions Notes

Abstract

Motorcycle motion is largely influenced by the road geometry, which alters the allowable accelerations in longitudinal and lateral directions and influences the vertical wheel loads. Recently, a method for three-dimensional road reconstruction and its incorporation into transient and quasi-steady-state (QSS) minimum lap time simulations (MLTSs) has been proposed. The main purpose of this work is to demonstrate how significantly different results from a minimum lap time optimal control problem can be obtained when using inappropriate elevation data sources in the track reconstruction problem, and how the road model reconstructed using poor input data can lead to misleading conclusions when analyzing real vehicle and driver performances. Two road models derived from high- and low-resolution digital elevation models (DEMs) are compared and their impact on the optimal maneuver of a racing motorcycle is examined. The essentials of track identification are presented, as well as vehicle positioning on the 3D road and the generalized QSS motorcycle model. Obtained 3D and 2D road models are analyzed in detail, on a case example of the Road Atlanta racetrack, and used in minimum lap time simulations, which are validated by the experimental data recorded on the Supersport motorcycle. The comparative analysis showed that great care should be taken when selecting the elevation dataset in the track reconstruction process, and that the 1 m resolution local DEMs seem to be sufficient to obtain MLTS results close to the measured ones. The example of using the 3D free-trajectory QSS minimum lap time problem to localize the track segments where real driver actions can be improved is also presented. The differences between simulation results on different road models of the same racetrack can be large and influence the interpretation of optimal maneuver.

Keywords:

optimal maneuver; road modeling; three-dimensional road; minimum lap time simulation; optimal control; motorcycles; experimental data; trajectory optimization; free-trajectory

1. Introduction

The influence of road geometry on vehicle motion is something we experience every day. When driving motorized or human-powered vehicles, we observe differences in acceleration capabilities when going uphill or downhill. We sometimes experience a sense of weightlessness when driving over the tops of hills, and we notice that corners are banked, even on public roads. Therefore, to accurately reproduce real vehicle motion on a racetrack, it is necessary to consider the road model, especially in minimum lap time simulations, where we expect to use results as reference data to improve the real vehicle and driver performance.

Historically, the minimum lap time simulations (MLTSs) were solved assuming a predefined path, e.g., [1,2,3,4]. Later on, the MTLS problem was formulated as an optimal control problem and the assumption of the free-trajectory has been introduced. In optimal maneuver problems, we can distinguish two vehicle models: transient (e.g., [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]) and quasi-steady-state [20,21,22,23,24]. Transient vehicle models provide information about the transient behavior of the vehicle. However, they are more computationally expensive than the free-trajectory quasi-steady-state (QSS) approach. First presented in [20], the QSS approach is characterized by a concise system dynamic. The vehicle performance is constrained by g-g diagrams describing the vehicle–driver performance. They are computed before solving the optimal control problem (OCP), which results in a lower computational complexity of the OCP itself. According to [20], the free-trajectory QSS approach is about ten times faster than when the dynamic vehicle model is used. The idea of g-g diagrams has been discussed, for example, in [25]. An experimentally-based method to simulate the transient behavior in QSS simulations has been presented in [24].

An alternative approach to solving the minimum lap time problems is a model predictive control (MPC) framework, in which the control is based on a finite preview horizon [26,27,28,29]. In [29], cascade optimization was used to mitigate the human driver’s ability to learn the track and to minimize the global maneuvering time.

There are also many publications related to road geometry and the trajectories of autonomous vehicles. These works include problems such as trajectory planning, tracking control, and obstacle avoidance, for example in [30,31,32,33,34]. The problem presented in this work concerns the specific and limited problem of developing an optimal trajectory for race motorcycles.

A detailed historical overview of MLTS problems, along with the evolution of vehicle models and a brief description of methods used to solve OCPs, has been provided in the review article [35]. Another comprehensive study on MLTSs is presented in [36].

The three-dimensional road model was introduced in [10]. The road reconstruction procedure from real track measurements was discussed in detail in [13]. The road was modeled as a ribbon and its geometry was determined using an additional OCP, which allowed filtering out of the noise present in the measurements. Then, in [14], the optimal maneuver of the F1 car on the Barcelona-Catalunya circuit was determined and compared with the results obtained when assuming a “flat” road model. The track model was built based on elevation data from the global digital elevation model.

The track reconstruction procedure was also described in [37], where the performances of direct and indirect methods for solving OCPs were directly compared.

A more sophisticated road model has been presented in [21], when lateral-position-dependent camber, typical for some racetracks with large lateral curvature, has been introduced. The Darlington Raceway racetrack was reconstructed using high-resolution LIDAR measurements. The free-trajectory QSS simulation of NASCAR car was presented and compared with a “flat” road model. The track reconstructed in [21] was also used in [19], where the single rigid body dynamic model of the NASCAR car was presented. The computed optimal maneuver was compared with calculation results from an industrial simulation tool and telemetry data recorded under race conditions.

The optimal maneuver on a 3D road was also investigated in [17,18], where the road model presented in [13] was adopted. The Whipple bicycle model has been used in the former, while the full dynamic multibody model of a sport motorcycle has been used in the latter. In [17], the optimal maneuver on the Mugello and Barcelona-Catalunya 3D racetracks is presented and compared with 2D simulation results and experimental data recorded by the on-board data acquisition system. In addition, the influence of road geometry on vehicle movement was discussed. In [18], the Mugello racetrack was considered. The source of the elevation data used in the track reconstruction procedure was not specified.

The three-dimensional free-trajectory QSS simulation has been presented in [22], together with the generalized vehicle model and the associated g-g-g diagram. The optimal maneuvers of a Formula 1 car and a MotoGP motorcycle on the Mugello and Barcelona-Catalunya racetracks have been analyzed and compared with simulations on a 2D road. The track geometry was taken from [17].

In [23], a similar approach as in [22] was presented. In addition, an optimal planning approach was presented, which finds a new time-optimal solution when the vehicle has to deviate from the original race line due to track conditions, for example, obstacle avoidance or overtaking maneuvers. The Mount Panorama Circuit in Bathurst and the Las Vegas Motor Speedway were considered. The source of the elevation data was not specified.

The three-dimensional fixed-trajectory approach has been presented in [38], with numerical examples provided for a car and a motorcycle.

In the presented work, the MLTS results are compared with experimental data from an on-board data acquisition system. A comparison of the simulation results with experimental data can be found in [6,8,15,17,19,24]. In [6], optimal gear ratios were searched. The calculated optimal maneuver for the Superbike motorcycle was compared with the telemetry data recorded on the Mugello and Adria racetracks. In [8], the optimal position of the vehicle’s center of gravity was studied on dry and wet roads. A 125cc motorcycle was used and the simulation results were compared with experimental data recorded at the Valencia circuit. The free-trajectory QSS simulation was validated with the measured data in [24], where the lightweight minibike was considered. The experimental data were also used in [15], where the qualifying lap in GP2 car at Barcelona-Catalunya circuit was compared with an optimal maneuver computed using the multibody car model.

A comprehensive review on road models for vehicular control is presented in [39], where the evolution of road models in MLTSs and comparative analysis between them are presented. However, none of the available literature addresses the importance of the elevation input data used in the 3D road identification, which can greatly influence the computed road model and the related optimal maneuver. This work focuses on the use of the elevation datasets that are available in the public domain, which makes them attractive for practical applications. The field of motorcycle racing is an example of such a need. Two high- and low-resolution digital elevation models are used to derive the elevation of the track edges of the real track and the ribbon approach is used to numerically reconstruct the road. The differences between the computed road models are highlighted and the simulation results are compared in detail with the data recorded under race conditions.

The main objective of this work is to demonstrate how significantly different optimal maneuvers can be obtained when using inappropriate elevation data sources, and how this can lead to misleading conclusions when analyzing real vehicle and driver performance. Section 2 describes vital information about road modeling, contains the formulation of the track reconstruction OCP, and discusses low-cost and no-cost elevation measurement methods. In Section 3, the essentials of vehicle model and vehicle positioning on a 3D road are presented. The minimum lap time optimal control problem is formulated. Section 4 presents the example of Road Atlanta racetrack reconstruction. Two three-dimensional road models, obtained using elevation data from low and high-resolution DEMs, are compared in detail. In Section 5, the MLTS results for the Supersport motorcycle are validated with experimental data recorded during second race of the 2021 MotoAmerica championship. An example of the use of 3D free-trajectory QSS simulation results to improve real driver actions is also presented. Conclusions are summarized in Section 6.

2. Road Modeling

The problem to be analyzed requires the formulation of two separate optimal control problems. In the first one, the three-dimensional road is identified from measured data, while in the second one, the optimal maneuver of the racing motorcycle is determined.

2.1. Track Identification OCP

The three-dimensional road is typically modeled as a ribbon generated by a spine curve, and this approach has already been described in detail in [13]. The road can be described by three coordinates of the spine curve

{[x (s), y (s), z (s)]}^{T}

, width

t_{w} (s)

, and three curvatures: relative torsion

Ω_{x} (s)

, normal curvature

Ω_{y} (s)

, and geodesic curvature

Ω_{z} (s)

, which are expressed in the ribbon (road) frame and describe its angular velocity. The road frame, called Darboux frame, is a moving trihedron located on the track spine curve in such a way that the

x

y

plane represents the road plane and the

x

-axis points along the tangent to the spine curve. The above quantities are functions of the curvilinear abscissa

s

, which represents the arc length along the spine curve. The curvilinear abscissa

s

varies between

s = 0

at the start of the road and

s = L

at the road endpoint.

The road curvatures

Ω_{B} = {[Ω_{x}, Ω_{y}, Ω_{z}]}^{T}

can be expressed in terms of Euler angles employed in the

z y x

convention:

θ (s)

ϕ (s)

, and

μ (s)

, called yaw (track heading), roll (road camber), and pitch (road inclination), respectively. The road curvatures

Ω_{B}

can be obtained using the following relationship:

Ω_{B} = [\begin{matrix} Ω_{x} \\ Ω_{y} \\ Ω_{z} \end{matrix}] = [\begin{matrix} 1 & 0 & - \sin μ \\ 0 & \cos ϕ & \cos μ \sin ϕ \\ 0 & - \sin ϕ & \cos μ \cos ϕ \end{matrix}] [\begin{matrix} ϕ' \\ μ' \\ θ' \end{matrix}] .

(1)

The three successive rotations using the Euler angles

θ

ϕ

, and

μ

describe the orientation of the ribbon frame relative to the inertial frame by the following rotation matrix

R = R_{z} (θ) R_{y} (μ) R_{x} (ϕ) = [\begin{matrix} c_{θ} c_{μ} & c_{θ} s_{μ} s_{ϕ} - s_{θ} c_{ϕ} & c_{θ} s_{μ} c_{ϕ} + s_{θ} s_{ϕ} \\ s_{θ} c_{μ} & s_{θ} s_{μ} s_{ϕ} + c_{θ} c_{ϕ} & s_{θ} s_{μ} c_{ϕ} - c_{θ} s_{ϕ} \\ - s_{μ} & c_{μ} s_{ϕ} & c_{μ} c_{ϕ} \end{matrix}],

(2)

where “s” and “c” symbolize the sine and cosine, respectively, and

R_{x} (ϕ)

R_{y} (μ)

, and

R_{z} (θ)

are elementary rotations about body-fixed axes given by

R_{x} (ϕ) = [\begin{matrix} 1 & 0 & 0 \\ 0 & \cos ϕ & - \sin ϕ \\ 0 & \sin ϕ & \cos ϕ \end{matrix}], R_{y} (μ) = [\begin{matrix} \cos μ & 0 & \sin μ \\ 0 & 1 & 0 \\ - \sin μ & 0 & \cos μ \end{matrix}], R_{z} (θ) = [\begin{matrix} \cos θ & - \sin θ & 0 \\ \sin θ & \cos θ & 0 \\ 0 & 0 & 1 \end{matrix}] .

(3)

The road geometry can be identified from three-dimensional coordinates of points located on the road edges. The position of these points can be determined by various methods, each of which will be characterized by a certain level of noise. The noisy road model in the related optimal maneuver problem would negatively affect the accuracy of the solution and would entail the need to reduce the integration step lengths, thus increasing the computation time.

To achieve smooth road characteristics, the track identification problem can be formulated as an optimal control problem, allowing for noise filtering. The cost function in OCP is given by

\min_{u \in U} I = \int_{s_{0}}^{s_{f}} F_{1} (s, X (s)) + F_{2} (s, v (s)) d s,

(4)

where

s_{0}

and

s_{f}

are the track lengths.

The state vector

X

consists of ten variables: the spine curve coordinates (

x, y, z)

, the Euler angles and their first derivatives, and the track width

X = {[x y z θ μ ϕ θ^{'} μ^{'} ϕ^{'} t_{w}]}^{T} .

(5)

The dynamics of the system are described by the following state equations:

X^{'} = f (s, X (s), v (s)) = [\begin{matrix} c_{θ} c_{μ} \\ s_{θ} c_{μ} \\ \begin{matrix} - s_{μ} \\ θ' \\ \begin{matrix} μ' \\ ϕ' \\ \begin{matrix} θ^{″} \\ μ^{″} \\ \begin{matrix} ϕ^{″} \\ t_{w}^{'} \end{matrix} \end{matrix} \end{matrix} \end{matrix} \end{matrix}] .

(6)

The control vector

v

consists of the second derivatives of the Euler angles and the first derivative of the track width

v = {[θ^{″} μ^{″} ϕ^{″} t_{w}^{'}]}^{T} = {[u_{θ} u_{μ} u_{ϕ} u_{w}]}^{T} .

(7)

The component

F_{1} (s, X (s))

under the integral in Equation (4) is equal to

F_{1} (s, X (s)) = {(x_{r} - x_{r_{0}})}^{2} + {(y_{r} - y_{r_{0}})}^{2} + {(z_{r} - z_{r_{0}})}^{2} + {(x_{l} - x_{l_{0}})}^{2} + {(y_{l} - y_{l_{0}})}^{2} + {(z_{l} - z_{l_{0}})}^{2},

(8)

and is intended to minimize the error between the output road edges and the provided dataset. The variables

x_{r_{0}}, y_{r_{0}}, z_{r_{0}}

represent coordinates of measured points on the right road edge, while

x_{l_{0}}, y_{l_{0}}, z_{l_{0}}

refer to the points located on the left road boundary. The coordinates of the numerically reconstructed edges,

x_{r}, y_{r}, z_{r}

and

x_{l}, y_{l}, z_{r}

, are related to the spine curve coordinates

x, y, z

by the following relations

x_{r} = [\begin{matrix} x \\ y \\ z \end{matrix}] + R [\begin{matrix} 0 \\ \frac{t_{w}}{2} \\ 0 \end{matrix}],

(9)

x_{l} = [\begin{matrix} x \\ y \\ z \end{matrix}] - R [\begin{matrix} 0 \\ \frac{t_{w}}{2} \\ 0 \end{matrix}],

(10)

where

R

is given by Equation (2).

The component

F_{2} (s, v (s))

in Equation (4) is expressed as

F_{2} (s, v (s)) = {w_{θ} u_{θ}}^{2} + w_{μ} u_{μ}^{2} + w_{ϕ} u_{ϕ}^{2} + w_{w} u_{w}^{2} .

(11)

The weights

w_{θ}

w_{μ}

w_{ϕ}

w_{w}

are used to control the input data filtering.

The cyclic boundary conditions on the states

X (s_{0}) = X (s_{f})

are needed when modeling the closed racetrack.

2.2. Acquiring and Preparing the Input Dataset

Each point in the input dataset must be described by three coordinates expressed in the global, fixed inertial frame. The process of acquiring and preparing the input dataset can be divided into three steps.

In the first step, track boundaries must be discretized. The points on the road edges can be located by GPS measurements, or by using computer programs that provide tools for exploring satellite images.

The second step is to assign the elevation data to the previously specified points. Low-cost or no-cost methods of obtaining elevation data include ellipsoidal height measurements using GPS and datasets available on the Internet, known as digital elevation models. A DEM is a discrete representation of terrain elevations, along with an interpolation algorithm that allows the calculation of elevations at any point within the area for which the model was created. There are two types of DEMs based on the resolution of the data:

Global DEMs covering large areas of the world, typically characterized by a resolution of 30 m or less;
Local DEMs provided by individual countries or regions, often with a high resolution of 5 m, 1 m, or even higher.

There are many free global DEMs available, and their use is facilitated by online services that automatically assign elevations to user-provided sets of points with known geographic coordinates. The use of local DEMs requires more effort due to differences between datasets, in terms of the file formats and reference systems used. Another challenge is related to the selective coverage of the world map, or sometimes the need to purchase access to the database.

The final step in the input data preparation process is to transform the geographic coordinates of the measured points into a rectangular coordinate system using a selected cartographic projection. In this step, the number of points on both road edges is standardized, the points on the track centerline are determined, and the arc length of the spine curve is estimated.

3. Three-Dimensional Free-Trajectory Quasi-Steady-State MLTS

3.1. Vehicle Positioning on the Three-Dimensional Road and Formulation of the Minimum Lap Time Problem

Vehicle positioning on a three-dimensional road is illustrated in Figure 1 and Figure 2. The

A_{\hat{x} \hat{y} \hat{z}}

frame is associated with the vehicle, while the

O_{t n m}

is a road frame (Figure 1). The

\hat{x}

-axis points in the direction of the tangent to the vehicle trajectory. In the following relationships, variables expressed in the

A_{\hat{x} \hat{y} \hat{z}}

frame are marked by a hat above the symbols. The vehicle frame is rotated by an angle

\hat{χ}

relative to the ribbon frame. Assuming that

Ω_{x}

and

Ω_{y}

are small, the vehicle accelerations

{\hat{a}}_{x}

and

{\hat{a}}_{y}

expressed in the vehicle frame are given by the following equations, provided, for example, in [22]:

{\hat{a}}_{x} = \dot{V},

(12)

{\hat{a}}_{y} = V (\dot{\hat{χ}} + Ω_{z} \dot{s}),

(13)

where

\dot{s}

is the speed of travelling along the spine curve, given by

\dot{s} = \frac{V \cos \hat{χ}}{1 - n Ω_{z}} .

(14)

Vehicle motion on a 3D road is described by the following state equations [21,22]:

V^{'} = \frac{{\hat{a}}_{x}}{\dot{s}},

(15)

n^{'} = \frac{1}{\dot{s}} V \sin \hat{χ},

(16)

{\hat{χ}}^{'} = \frac{1}{\dot{s}} \frac{{\hat{a}}_{y}}{V} - Ω_{z},

(17)

where

V

is the absolute velocity along the trajectory,

n

is the vehicle position along the ribbon frame unit vector

n

, and

\hat{χ}

is the relative orientation of

V

, with respect to the tangent to the road spine curve (Figure 2). The Equations (15)–(17) are expressed in the space domain and have been rewritten from the time domain using the following derivation rule:

x^{'} = \frac{d x}{d s} \frac{d t}{d t} = \frac{d x}{d t} \frac{d t}{d s} = \frac{\dot{x}}{\dot{s}} .

(18)

To control the vehicle, acceleration derivatives called longitudinal jerk

{\hat{j}}_{x}

and lateral jerk

{\hat{j}}_{y}

were introduced by the authors, as in [24]. As a result, two additional state equations were obtained:

{\hat{a}}_{x}^{'} = {\hat{j}}_{x},

(19)

{\hat{a}}_{y}^{'} = {\hat{j}}_{y} .

(20)

As a result, the state vector

x

consists of five entries:

x = {[V, n, \hat{χ}, {\hat{a}}_{x}, {\hat{a}}_{y}]}^{T},

(21)

while the control vector

u

has two elements

u = {[{\hat{j}}_{x}, {\hat{j}}_{y}]}^{T} .

(22)

The goal of OCP is to minimize the maneuver time. After spatial reformulation the target

L

is defined as

L = \int_{t_{0}}^{t_{f}} 1 d t = \int_{s_{0}}^{s_{f}} \frac{1}{\dot{s}} d s,

(23)

where

t_{0}

and

t_{f}

are the start and end times of control, respectively, while

s_{0}

and

s_{f}

are the track lengths.

The vehicle must respect the road boundaries

- t_{w} / 2 \leq n \leq t_{w} / 2,

(24)

where

t_{w} = t_{w} (s)

is a variable track width.

According to [24], the maximum allowable values of longitudinal and lateral jerks can be limited, respectively, by the following inequalities:

{\hat{J}}_{x_{d}} \leq {\hat{j}}_{x}^{} \leq {\hat{J}}_{x_{a}},

(25)

{\hat{j}}_{y}^{2} \leq {\hat{J}}_{y}^{2},

(26)

where

{\hat{J}}_{x_{d}} = {\hat{J}}_{x_{d}} (V)

{\hat{J}}_{x_{a}} = {\hat{J}}_{x_{a}} (V)

, and

{\hat{J}}_{y} = {\hat{J}}_{y} (V)

are the speed-dependent hyperbolic functions given by general equation

\frac{d}{d s} ({\hat{a}}_{r}) = {\hat{J}}_{r} = \frac{β_{0}^{r}}{V} + β_{1}^{r} .

(27)

The

r

takes the form

x_{a}

x_{d}

, and

y

depending on selected constraint (25) and (26). The function

{\hat{J}}_{x_{d}}

applies to the negative gradients of the longitudinal acceleration

{\hat{a}}_{x}

, while

{\hat{J}}_{x_{a}}

is a constraint for positive gradients of

{\hat{a}}_{x}

. The coefficients

β_{0}^{r}

and

β_{1}^{r}

can be determined based on the jerks recorded during the actual vehicle’s motion. This method has been validated and described by the authors in detail in [24].

3.2. Vehicle Model and Related Vehicle Performance Constraint

Assuming quasi-steady-state conditions, vehicle performance is limited by the envelope of the g-g diagram, referred to as the g-g-g diagram in the case of spatial motion on a 3D road. The g-g or g-g-g diagram represents the limit of tire adhesion when lateral and longitudinal accelerations occur simultaneously, taking into account the vehicle’s ability to accelerate and decelerate. The shape of the g-g diagram envelope is speed-dependent, primarily due to aerodynamic drag and lift.

Figure 3 shows a generalized rigid-body model of a vehicle. The equations describing the model are the same as in [22]. In addition, the model includes rolling resistance and considers a speed-dependent driving force,

F_{x r}

, analogous to [24]. Yaw and side-slip accelerations are omitted, the suspension is rigid, the steering angle is zero, the tires are infinitely narrow discs, and the driver is rigidly attached to the chassis. The longitudinal friction coefficient

μ_{x}

and the lateral friction coefficient

μ_{y}

are the same for the front and rear tires. The vehicle parameters introduced in Figure 3 are explained in Table 1, where their values for the Supersport motorcycle used in the simulations are given. Vehicle parameters were estimated based on measurements made on a motorcycle in very similar specification to those used in the MotoAmerica championship. Aerodynamic properties and tire friction coefficients were estimated from recorded data.

The variables

{\tilde{a}}_{x}, {\tilde{a}}_{y},

and

\tilde{g}

shown in Figure 3, and the vehicle speed

V

, allow the parameterization of the g-g-g diagram. The relationship between

{\tilde{a}}_{x}, {\tilde{a}}_{y}, \tilde{g},

and the accelerations

{\hat{a}}_{x}

and

{\hat{a}}_{y}

expressed in the vehicle frame can be determined by analyzing the motion of the vehicle’s CoG. Assuming constant CoG position relative to the road surface equal to

{[0, 0, - h]}^{T}

, the following relationships can be obtained:

{\tilde{a}}_{x} = {\hat{a}}_{x} - g (\cos μ \sin ϕ \sin \hat{χ} - \sin μ \cos \hat{χ}),

(28)

{\tilde{a}}_{y} = {\hat{a}}_{y} - g (\sin μ \sin \hat{χ} + \cos μ \sin ϕ \cos \hat{χ}),

(29)

\tilde{g} = g \cos μ \cos ϕ + V {\hat{ω}}_{y},

(30)

where

{\hat{ω}}_{y}

is the angular velocity of the vehicle frame about the

\hat{y}

-axis. The angular velocity

\hat{ω}

of the vehicle frame can be determined from a matrix equation

\hat{ω} = [\begin{matrix} {\hat{ω}}_{x} \\ {\hat{ω}}_{y} \\ {\hat{ω}}_{z} \end{matrix}] = R {(e_{z}, \hat{χ})}^{T} [\begin{matrix} Ω_{x} \\ Ω_{y} \\ Ω_{z} \end{matrix}] \dot{s} + [\begin{matrix} 0 \\ 0 \\ \dot{\hat{χ}} \end{matrix}] = [\begin{matrix} (Ω_{x} \cos \hat{χ} + Ω_{y} \sin \hat{χ}) \dot{s} \\ (Ω_{y} \cos \hat{χ} - Ω_{x} \sin \hat{χ}) \dot{s} \\ Ω_{z} \dot{s} + \dot{\hat{χ}} \end{matrix}],

(31)

where

R {(e_{z}, \hat{χ})}^{T}

is the rotation matrix used to express the angular velocity

Ω_{B}

of the road frame in the vehicle frame. Equation (31) is the sum of the angular velocity

Ω_{B}

of the ribbon frame and the angular velocity of the vehicle frame relative to the ribbon frame equal to

{[0, 0, \dot{\hat{χ}}]}^{T}

The effect of road geometry on the g-g diagram is illustrated in Figure 4a. Road camber

ϕ

affects the maximum lateral acceleration by shifting the envelope of the g-g diagram along the abscissa axis. Road pitch

μ

affects the allowable longitudinal acceleration by shifting the g-g diagram along the ordinate axis. The term

{\hat{ω}}_{y}

represents variations in normal road curvature (which may be associated with driving through a valley bottom or over the crest of a hill) and scales the envelope of the graph.

The g-g-g diagrams generated for different vehicle speeds

V

and

\tilde{g}

can be represented using hypersurfaces shown in Figure 4b. The parameters

\tilde{ρ}

and

\tilde{α}

, called the adherence radius and its orientation, respectively, are related to the description of the g-g diagram envelope in polar coordinates and are given by

\tilde{ρ} = \sqrt{{(\frac{{\tilde{a}}_{x}}{g})}^{2} + {(\frac{{\tilde{a}}_{y}}{g})}^{2},}

(32)

\tilde{α} = \arctan 2 ({\tilde{a}}_{y}, {\tilde{a}}_{x}) .

(33)

The vehicle’s performance in the OCP is limited using the following relationship:

\tilde{ρ} \leq {\tilde{ρ}}_{m a x} (\tilde{α}, V, \tilde{g}),

(34)

where

{\tilde{ρ}}_{m a x}

represents boundary of the g-g-g diagram.

A more detailed description of the construction of g-g diagrams and related hyper-surfaces can be found in the examples of [20,22,40].

4. Results of Track Reconstruction

The influence of the road model on the optimal maneuver is discussed using the example of the Supersport motorcycle and the Road Atlanta circuit located in the United States. The simulation results are divided into two sections. The first section includes a comparison of the road model obtained using global and local DEMs. In the second section, a comparison of the optimal maneuver computed on adopted road models is presented. The 3D simulations are compared with experimental data and the simulation on a “flat” two-dimensional road.

The formulated OCPs were solved using the orthogonal collocation method implemented in the GPOPS-II package [41] linked to the MATLAB R2019a software. The IPOPT and the ADiGator [42] were used as the optimizer and automatic differentiation tools, respectively. Numerical differentiation was employed in the three-dimensional minimum time maneuver problem.

Comparison between Results Obtained Using Local and Global DEMs

The track edges were discretized using Google Earth Pro 7.3.6 software. A total of 1355 points were manually selected (652 on the right and 703 on the left side of the road). Due to discrepancies between satellite images taken in different years, the image closest to the orthophoto map available in [43] was selected based on analysis of the location of several checkpoints spread over different sections of the track.

Three road models have been built:

In the first model, later referred to as 3D-local, the elevation of the points was determined using a local DEM provided in TIFF file format by the United States Geological Survey [43]. This database has a resolution of 1 m and 0.1 m root mean square elevation error;
In the second model, hereafter referred to as 3D-global, the elevation data were assigned using a tool available at [44]. The assigned elevation data were derived from the NED1 global DEM, characterized by 1 arc-second horizontal resolution;
In the third case, the simple two-dimensional track model was built—later referred to as 2D. It was obtained by omitting variables $z, μ, ϕ, μ^{'}, ϕ'$ in the state vector and $μ^{″}, ϕ ″$ in the control vector in the track identification OCP formulated in Section 2.1.

The geographical coordinates of the points located on the track edges were converted to a rectangular coordinate system using the Universal Transverse Mercator projection [45].

The Road Atlanta track contains numerous hills and valleys, banked corners, and turns with large and small curvatures. The latter caused difficulties in choosing a universal value for

w_{θ}

, which is mainly responsible for smoothing the geodesic curvature

Ω_{z}

. A small value of

w_{θ}

allowed accurate representation of corners with small radii of curvature, but led to an increase in simulation time in the associated minimum lap time problem. As the road geometry became more detailed, the number of mesh points required increased. On the other hand, a large value of

w_{θ}

made it possible to obtain a smooth geodesic curvature

Ω_{z}

, but this was at the expense of accuracy in the representation of corners with small curvature radii. As a result, the error in the tight corners between the numerically identified edges and the provided data increased significantly.

It was therefore decided to “manually” combine the results obtained for a penalty factor

w_{θ} = 4 \times 10^{4}

and

w_{θ} = 4 \times 10^{6}

. The results for

w_{θ} = 4 \times 10^{4}

have been used within turns from 3a to 3c, 10a, 10b, and, in turn, 12. A better representation of the track boundaries was achieved without any negative effect on the computation time in the associated lap time minimization problem. Other penalty factors were chosen as:

w_{ϕ} = w_{μ} = 4 \times 10^{6}

and

w_{w} = 1 \times 10^{4}

. The discrepancies between provided input data and reconstructed track boundaries in tight corners, depending on the chosen value of

w_{θ}

, are shown in Figure 5.

Figure 6a shows a two-dimensional map of the numerically reconstructed Road Atlanta track, while Figure 6b shows a comparison of computed curvatures

Ω_{x}, Ω_{y}

, and

Ω_{z}

. The Road Atlanta circuit is a clockwise racetrack, consisting of ten right-hand turns and seven left-hand corners. The longest straight is the main straight, which is about 500.0 m. The Road Atlanta track is characterized by numerous uphill and downhill sections. The most famous part of the track is the downhill section between turns 3 and 5 called ‘The Esses’.

The elevation zero point in both 3D road models is located at the left end of the finish line. The elevation of the left boundary is shown in Figure 7a, while the right boundary altitude is presented in Figure 7b. The height difference between track edges is shown in Figure 8. The difference between the lowest and the highest point of the track is 37.7 m in the 3D-local case and 40.8 m in the 3D-global case. The lowest point of the track is in turn 4c, while the highest point is in turn 11.

The differences between the computed road models are clearly visible in Figure 9, where the road camber

ϕ

and road inclination

μ

are shown. The maximum value of

ϕ

in the 3D-local case is

5.3 °

(turn 6), while the minimum is

- 1.6 °

at turn 5. In the 3D-global case, the maximum and minimum camber angles

ϕ

are

10.5 °

at turn 11 and

- 10.3 °

after turn 7, respectively. The highest absolute difference between the road models is

10.6 °

, which occurs at about

s = 2280

m. The road camber in the 3D-global case has the opposite sign to the 3D-local case for the large part of the track.

The road inclination

μ

shows greater agreement between the models, although in the 3D-global case there are clearly visible fluctuations. The maximum and minimum value of

μ

are, respectively, equal to

6.0 °

(after turn 1) and

- 8.8 °

(after turn 10b) in the 3D-local case. The extremes of

μ

in the 3D-global model are in similar locations and are equal to

6.5 °

and

- 10.5 °

. The highest absolute difference between road models is

3.2 °

and occurs before turn 10a.

Magnified 3D views of selected areas of the track with the largest discrepancies between the models are shown in Figure 10.

Differences in the road camber

ϕ

between the road models result in variations in the calculated track width

t_{w}

, as shown in Figure 11, which also includes the two-dimensional case. Track widths range from 11.5 to 12.7 m, with the minimum at curve 4a. The maximum difference in track width between the 3D-local and 3D-global cases is equal to 0.06 m and occurs at turn 4d. The width of the 2D road is always smaller because the camber angle is omitted. There are also differences in the estimated length of the track spine curve, which is equal to 4144.5 m in the 3D-local case, 4145.5 m in the 3D-global case, and 4138.9 m in the 2D case. The longer spine curve in the 3D cases results from the non-zero inclination angle

μ

5. Optimal Maneuver Simulation Results

The simulation results are compared with experimental data from the Supersport motorcycle used in the MotoAmerica championship. These data are available on the website [46] of one of the riders participating in the series. The shared database includes GPS data, and measurements from the engine control unit and numerous analogue and digital sensors, including suspension potentiometers and wheel speed sensors.

The model of the Supersport motorcycle was created based on measurements taken on a motorcycle similar to the one used by the American rider. The vehicle parameters used in the simulations are shown in the Table 1. The longitudinal friction coefficient

μ_{x}

and the lateral friction coefficient

μ_{y}

of the tires were determined by trial and error, based on a comparative analysis between simulation results and experimental data.

Figure 12a shows the torque and power curves of the engine obtained on a chassis dynamometer in a transient condition. The driving force for gears 2–6 is shown in Figure 12b. Gear number 1 was omitted because it was not used during the actual vehicle motion. The stock gear ratios were assumed, while the unknown final drive ratio was estimated based on the relationship between engine speed and vehicle speed, as measured by GPS.

The coefficients

β_{0}^{r}

and

β_{1}^{r}

included in the hyperbolic constraints (26) and (27) on the allowable values of

{\hat{j}}_{x}

and

{\hat{j}}_{y}

were determined by the method described in [24] and are presented in Table 2. They have been calculated using the entire shared database, which includes measurements from several different racetracks. Furthermore, a semi-empirical g-g diagram envelope described in [24] was used, where the maximum negative longitudinal acceleration

{\hat{a}}_{x}

is obtained when driving straight forward.

Three simulations were performed for the 3D-local, 3D-global, and 2D road models, then compared with the best lap from the provided database (later referred to as AQ).

In the graphs presented later, the elapsed distance

s_{d i s t}

was used as the independent variable instead of the curvilinear abscissa

s

, to match the experimental data domain. The experimental longitudinal and lateral accelerations presented on the following graphs were determined using GPS data post-processing.

5.1. Comparison of Simulation Results on Adopted Road Models

The calculated optimal lap times are 89.658 s in the 3D-local case, 90.844 s in the 3D-global case, and 91.211 s on the 2D road model. The best recorded lap time was 90.485 s (AQ), which is 0.829 s slower than the lap time in the 3D-local simulation.

Figure 13 shows the measured and computed vehicle speed

V

. The 3D-local case is closest to the experiment. The highest difference occurs in turn 10a and is equal to 3.0 m/s. The deviations in other corners do not exceed 1.6 m/s. For the other two simulations, there are significant differences in numerous corners. In the case of the 3D-global model, the largest deviations are observed in corners 5, 6, chicane 10a–10b, and turn 12. They are equal to 3.2 m/s, −3.6 m/s, −3.0 m/s, and 3.6 m/s, respectively. For the two-dimensional road model, the largest deviations are observed in corners 1, 2, and 6, where they differ from the experiment by −1.9 m/s, 3.6 m/s, and −2.8 m/s, respectively. The three-dimensional road models allowed for a better representation of the vehicle speed on sections with significant road inclination

μ

(between turns 2 and 3, and turns 4 and 5).

The maximum speed is achieved between turns 9 and 10 and is equal to 74.2 m/s (AQ), 73.2 m/s (3D-local), and 72.9 m/s in the 3D-global and 2D cases. The minimum speed in the slowest corner number 7 is 20.5 m/s (AQ), 20.2 m/s (3D-local), 19.2 m/s (3D-global), and 20.5 m/s (2D). The measured maximum speed is 2 m/s higher than the average maximum speed of the other laps and was reached in the aerodynamic shadow behind another rider. The thin, black line in Figure 13 represents a speed trace from the lap with a maximum speed closest to the calculated average value. Due to the aerodynamic shadow, the real vehicle gains about 0.210 s. The vehicle in the 3D-local case gains the most advantage over the AQ during braking before turns 3a, 6, 10a, 12, and on the straight section between turns 5 and 6, which is analyzed in detail in Section 5.3. In the 3D-global case, the vehicle has an advantage (over AQ) of 0.606 s at

s_{d i s t} = 1850

m, but then loses 1.505 s until

s_{d i s t} = 3607

m. A significantly higher speed in the last corner reduces the time loss at the end of the lap to 0.359 s. The vehicle from the 2D simulation is slower than the AQ by no more than 0.350 s until turn 6. In corner 6 and 7, the vehicle loses approximately 0.400 s. After that, the time difference remains approximately constant and is equal to 0.726 s at the end of the lap.

Figure 14 shows the vehicle’s longitudinal acceleration

{\hat{a}}_{x}

and lateral acceleration

{\hat{a}}_{y}

. The good agreement with the experimental data, in terms of the rate of change of the accelerations in the transient phases, is a consequence of the hyperbolic constraints introduced on the controls. The three-dimensional road model allowed a good prediction of the longitudinal acceleration

{\hat{a}}_{x}

during the thrust phase. This can be attributed to the inclusion of 3D effects (road inclination

μ

in this case) and the speed-dependent driving force in the vehicle model. In the 2D road model case, higher values of longitudinal acceleration

{\hat{a}}_{x}

are observed in the initial acceleration phase. For both simulations on the three-dimensional road, greater deceleration occurs in sections where the road flattens (1808–1882 m, 3309–3424 m, 3780–3864 m). The “simulated optimal rider” takes advantage of the higher vertical tire load resulting from changes in the road normal curvature

Ω_{y}

. This will be explained in Section 5.2.

The largest deviation in lateral acceleration between 3D-local and AQ occurs at turn 2 and is equal to −

0.33

g. In general, the 3D-local road model was able to faithfully represent the maximum values of longitudinal acceleration in the corners.

The differences between the 3D-global and AQ lateral acceleration plots are directly translated into the obtained corner speed. In the corners where the largest difference in vehicle speed was observed, the lateral acceleration was underestimated by 0.39 g (turn 6) and overestimated by 0.58 g (turn 12).

The vehicle performance in the performed simulations can be conveniently compared by superimposing the simulation results on the experimental g-g diagram, as shown in Figure 15. Again, the best agreement was obtained for the 3D-local case (Figure 15a), especially in the left half of the g-g diagram. The extreme of the lateral acceleration

{\hat{a}}_{y}

in the left corners is equal to −1.25 g (turn 5), while in the right corners it is 1.47 g (turn 6, where the camber angle

ϕ

is maximal). The obtained results are close or equal to the experimental ones, which are −1.25 g and 1.46 g, respectively.

In the 3D-global case, the maximum lateral acceleration

{\hat{a}}_{y}

in left corners is equal to −1.64 g and 1.82 g in the right-hand turns. The peak of 1.82 g at

s_{d i s t} = 3887

m seen in both Figure 14 and Figure 15 corresponds to the peak of normal curvature

Ω_{y}

s = 3947

m in Figure 6b. These discrepancies result from the low resolution of the global DEM, which prevents a satisfactory reconstruction of the road camber.

In the case of the two-dimensional road, the g-g diagram has an elliptical shape, typical for the motion on a “flat” road (Figure 15c). The maximum lateral acceleration

{\hat{a}}_{y}

the vehicle is able to achieve when cornering is constant (when the lift force is neglected) and depends on the tire lateral friction coefficient

μ_{y}

Figure 16 shows the vehicle lateral position on the track and the radius of curvature modulus, which can be obtained using the following equation:

|Ρ| = \frac{V^{2}}{|{\hat{a}}_{y}|} .

(35)

The largest difference in vehicle lateral position between the simulations occurs on the finish straight, between turns 8 and 9, and in turn 11 when the vehicles are going over the hill (Figure 17a). The 2D case is closer to the 3D-local case than the simulation in which the 3D-global road model was used. According to the Figure 16b, which shows the radius of curvature modulus, the vehicle in the 3D-global case takes slightly different lines at the exit of turn 1, in turn 4, and in the chicane 10a–10b, which is shown in the enlarged view in Figure 17. The radius of curvature is significantly smaller in the first corner of the chicane (36.6 m vs. 62.5 m in the 3D-local case), but larger in the second corner (56.0 m vs. 48.3 m). Corner 10b in the 3D-global road model has a road camber of

- 0.1 °

, while in the 3D-local case the road banking is

1.7 °

. Therefore, a wider trajectory is required to achieve a similar speed at the corner apex. The vehicle in the 3D-global case prioritizes the exit speed in turn 10b, which is before the straight section of the track, over the corner speed in the first corner of the chicane.

5.2. Analysis of Tire Vertical Load

Road geometry affects the vertical tire load according to the relationship

m \tilde{g} = N_{r} + N_{f},

(36)

which belongs to the equations of the vehicle model. The variable

\tilde{g}

is called apparent gravity and is given by the Equation (30). Based on the vehicle model, the following equations can be derived for the vertical load on the front tire

N_{f}

and the vertical load on the rear tire

N_{r}

N_{r} = \frac{(m {\tilde{a}}_{x} h + F_{d} h_{p}) \cos \tilde{φ} + (w - b) m \tilde{g}}{w},

(37)

N_{f} = m \tilde{g} - N_{r},

(38)

where

\tilde{φ}

is a vehicle roll angle relative to the road surface given by

\tilde{φ} = \arctan (\frac{{\tilde{a}}_{y}}{\tilde{g}}) .

(39)

Variations in road geometry can significantly affect normal tire load, increasing the likelihood of front or rear wheel lift during hard acceleration or braking, as well as altering allowable tire longitudinal and lateral forces.

The calculated vertical load on the front wheel, shown in the upper graph of Figure 18, makes it possible to predict (in 3D-local case) that the front wheel lift can occur in two track segments during the ride over the top of the hills, after turns 5 and 11 (between the vertical dotted lines marked in Figure 18). The simulation results are confirmed by the measurements recorded by the data acquisition system. The front wheel lift is indicated by local dips in the front wheel speed and the near zero value of the front suspension deflection (middle graph on Figure 18). At the top of the hill after turn 5, the front suspension was in the position of maximum extension for more than 40 m, while the linear velocity of the front wheel speed decreased from 37.5 m/s to 32.8 m/s. In addition, the actual rider had to stop the pitching motion toward the rear wheel by partially closing the throttle at

s_{d i s t} = 1488

m (bottom plot on Figure 18).

In the second identified track segment, the same features mentioned above are visible in the experimental data. There are numerous dips in the front wheel speed, and the suspension deflection is close to zero. The rider and the “simulated driver” try to counteract the front wheel lift by driving over the hill with a non-zero lateral acceleration

{\hat{a}}_{y}

, as shown in Figure 14b. According to the g-g diagram shown in Figure 4a, leaning the motorcycle to a certain extent allows a higher longitudinal acceleration than in straight motion.

When comparing the OCP results for the adopted road models, it can be seen again that the 3D-local case is closest to the experiment. In the 3D-global case, the vertical front wheel load in the identified road segments also takes values close to zero, but occurs in the shorter distance. A “flat” road model gives significantly different results. The vertical load on the front wheel in the aforementioned road segments is significantly higher. By comparing the results for 3D and 2D roads, it is possible to determine whether the low vertical load on the front tire is only due to high longitudinal acceleration, or whether it is additionally enhanced by varying road geometry.

Another interesting conclusion can be drawn from Figure 19, which shows the total tire vertical load. Intuitively, the total tire vertical load on a two-dimensional road is constant and equal to

m g

. For the 3D road models, the vertical load on the tires increases when driving through a valley bottom, or decreases when driving over a hill, as in the road sections shown above. This is an effect of centripetal force in a vertical direction, resulting from road curvature

Ω_{y}

variations. The maximum total vertical force is equal to 3546 N (3D-local) and 4527 N (3D-global), while the minimum value is 1349 N (3D-local) and 523 N (3D-global). These large differences between extremes are due to significant variations in normal curvature

Ω_{y}

in the 3D-global road model. These variations are also visible on the total vertical force plot, as significantly larger fluctuations than in the 3D-local case.

5.3. Detailed Analysis of Selected Track Segment

The last example illustrates the possibility of using the results of the three-dimensional quasi-steady-state minimum lap time simulation to improve the actual rider’s actions. A detailed analysis of the track section between turns 5 and 6 is given. The experimental results are compared with the 3D-local simulation. The analyzed section is shown on the left side of the Figure 20, while the right side includes graphs of vehicle speed

V

, longitudinal acceleration

{\hat{a}}_{x},

lateral acceleration

{\hat{a}}_{y}

, curvature radius

Ρ

modulus, and selected measurements from the data acquisition system: engine speed

n_{e}

, throttle position, and selected gear.

The simulated bike is 1.0 m/s slower in turn 5, but reaches a more than 2.9 m/s higher top speed at the end of the straight between turns 5 and 6. The time difference between the vertical dotted lines shown in top-right plot is 0.166 s, in favor of the simulation.

Despite a fast and smooth throttle opening by the actual driver, and almost the same start of the thrust phase (indicated by asterisks on the left plot in Figure 20), the maximum recorded longitudinal acceleration is lower than in the optimal maneuver. The graph of the driving force (Figure 12b) shows that, for a given vehicle speed, the second gear should be used instead of the third gear by the actual driver. The smaller radius of curvature of the vehicle trajectory in the 3D-local case, together with the lower speed at the apex of turn 5, allowed the “simulated driver” to straighten the motorcycle earlier and use more of the available driving force.

Another factor in time lost is related to a premature shift from third to fourth gear after turn 5. The rider changed gear at the top of the hill when the front wheel lost contact with the road. The gear change was made at too low an engine speed, resulting in underutilization of available driving force after the top of the hill, when longitudinal acceleration was no longer limited by front wheel lift. In addition, the vehicle in 3D-local simulation travels slightly less distance, due to the straightening of the trajectory between

s_{d i s t} = 1500

m and

s_{d i s t} = 1700

Incorrect gear before turn 5 and a too early gear change after turn 5 resulted in a speed difference of almost 3.0 m/s at the end of the following straight section of the track.

6. Conclusions

The influence of the road model on the optimal maneuver of the racing motorcycle has been studied. The essentials of the road reconstruction problem have been presented, as well as vehicle positioning on the 3D road and generalized QSS vehicle model. Three road models were constructed: the “flat” road model and two three-dimensional road models, reconstructed using elevation data from 1 m and 1 arc-second resolution digital elevation models. A 3D free-trajectory quasi-steady-state minimum lap time problem with experimentally bounded controls (jerks) was employed and compared in numerical examples to experimental data from the Supersport motorcycle used in the MotoAmerica championship.

The use of digital elevation models is convenient, because it does not require any additional work to be done on the reconstructed track, nor does it (usually) incur additional costs. In general, DEMs can be divided according to the resolution of the data provided. The use of low-resolution global DEMs is encouraged by numerous free tools that automate the process of assigning elevation data to datasets containing geographic coordinates. Local DEMs require more user effort due to differences in file formats and reference systems used.

Significant differences were found between road models of the Road Atlanta racetrack reconstructed using DEMs of different resolutions. Differences were particularly noticeable in the road camber angle and relative torsional curvature, while better agreement was observed in road inclination and normal curvature. The road model derived from the global DEM exhibited notable fluctuations in road geometry parameters, except for lane width and geodesic curvature, which were less dependent on elevation data.

The comparative analysis between the computed optimal maneuvers and the experimental data showed that the results obtained from the track, reconstructed using the local DEM, closely matched the experimental data compared to two other road models. This agreement is particularly true for predicted corner speed, lateral acceleration extremes, and g-g diagram shape. The largest discrepancies in corner speed occurred in the simulation using the global DEM-based road model, exceeding those observed in the 2D case. The 1 arc-second resolution of the global DEM proved insufficient for accurate road camber identification. On tracks with significant variations in road geometry, such as the Road Atlanta racetrack analyzed, the “flat” road model provides limited information on vehicle motion. The maximum lateral acceleration is constant and depends only on the lateral friction coefficient.

Both three-dimensional road models allowed better prediction of vehicle speed on inclined track segments, with superior results observed in the simulation using the local DEM road model.

The total wheel load is greatly influenced by the road geometry as presented in the article. In the 2D simulation, the total vertical force acting on the wheels is also constant and equal to the product of the standard gravity and the sum of the vehicle and rider masses, while in the 3D simulations it is influenced by variations in the road’s normal curvature. The total vertical load is higher when the vehicle travels through a valley bottom and lower when the vehicle travels over the crest of a hill. Both three-dimensional road models enabled the prediction of track segments with a high probability of front wheel lift resulting from normal curvature variations, which was consistent with suspension position and wheel speed measurements taken during actual vehicle motion. Again, the local DEM road model provided results closer to the measured ones.

The last numerical example showed that the three-dimensional quasi-steady-state simulation on a detailed 3D road, in the presence of the experimentally limited controls (jerks), can be used to localize the track segments where the real rider actions can be improved. Based on the comparative analysis of the OCP results and the recorded data, the reasons for the real rider’s time loss on the selected track section were identified. By incorporating the detailed road model and the jerk’s constraints into the minimum lap time optimal control problem, it was possible to locate the premature gear change and find the possibility of modifying the vehicle trajectory to fully open the throttle earlier when exiting the corner onto the following straight.

As presented in this work, great care should be taken when using the global DEMs in the track identification, as this may later lead to significant differences in the results of the related optimal maneuver problem. Therefore, the OCP results obtained using a poor track model may lead to erroneous conclusions in the comparative analysis with the experimental data. The presented examples showed that the 1 m resolution local DEMs seem to be sufficient to obtain MLTS results close to the measured ones, and these results have a high potential to be used as reference data in the process of vehicle and driver performance improvement.

Author Contributions

Conceptualization, J.B. and M.P.; methodology, J.B. and M.P.; software, J.B.; validation, J.B.; formal analysis, J.B. and M.P.; investigation, J.B.; resources, J.B.; data curation, J.B.; writing—original draft preparation, J.B.; writing—review and editing, J.B. and M.P.; visualization, J.B.; supervision, M.P.; project administration, J.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The digital elevation models used, and the experimental data presented are available on the Internet (accessed on 6 April 2024), and the corresponding URLs are included in the references [43,46].

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

CoG	center of gravity
CoP	center of pressure
DEM	digital elevation model
MLTS	minimum lap time simulation
MPC	model predictive control
OCP	optimal control problem
QSS	quasi-steady-state
List of symbols:
${\hat{a}}_{x}$	longitudinal acceleration
${\tilde{a}}_{x}$	longitudinal acceleration in generalized vehicle model
${\hat{a}}_{y}$	lateral acceleration
${\tilde{a}}_{y}$	lateral acceleration in generalized vehicle model
$b$	longitudinal position of CoG
$C_{d} A$	drag area coefficient
$F_{d}$	drag force
$F_{x f}, F_{x r}$	longitudinal force acting on the front and rear wheel, respectively
$F_{y f}, F_{y r}$	lateral force acting on the front and rear wheel, respectively
$g$	standard gravity
$\tilde{g}$	apparent gravity
$h$	height of CoG
$h_{p}$	height of CoP
$i_{g}$	gear ratio
$i_{p}$	primary transmission ratio
$i_{s}$	secondary transmission ratio
${\hat{J}}_{x_{a}} (V)$	longitudinal jerk constraint for positive gradients of the longitudinal vehicle acceleration ${\hat{a}}_{x}$
${\hat{J}}_{x_{d}} (V)$	longitudinal jerk constraint for negative gradients of the longitudinal vehicle acceleration ${\hat{a}}_{x}$
${\hat{J}}_{y} (V)$	lateral jerk constraint
${\hat{j}}_{x}$	longitudinal jerk
${\hat{j}}_{y}$	lateral jerk
$I, L$	cost function
$L$	track length
$m$	total mass (vehicle + rider)
$N_{f}, N_{r}$	vertical force acting on the front and rear wheel, respectively
$n$	lateral vehicle position with respect to the track spine curve
$n_{e}$	engine speed
$R$	rotation matrix
$r_{t}$	tire rolling radius
$s$	curvilinear abscissa
$s_{0}, s_{f}$	start and end track length
$s_{d i s t}$	elapsed distance
$\dot{s}$	speed of travelling along the spine curve
$t_{0}, t_{f}$	start and end time of control
$t_{w}$	track width
$u, v$	control vectors
$V$	absolute velocity along the trajectory
$w$	wheelbase
$w_{θ}, w_{μ}, w_{ϕ}, w_{w}$	weights in the cost function in the track reconstruction problem
$x, X$	state vectors
$x, y, z$	spine curve cartesian coordinates
$x_{l}, y_{l}, z_{l}$	left road border cartesian coordinates
$x_{r}, y_{r}, z_{r}$	right road border cartesian coordinates
$\tilde{α}$	angular coordinate (orientation)
$β_{0}, β_{1}$	coefficients used in jerks’ constraints
$θ$	yaw (track heading)
$μ$	pitch (road inclination)
$μ_{x}$	longitudinal friction coefficient
$μ_{y}$	lateral friction coefficient
$Ρ$	radius of curvature
$ρ_{a}$	air density
$\tilde{ρ}$	adherence radius
${\tilde{ρ}}_{m a x}$	maximum value of adherence radius
$ϕ$	roll (road camber)
$\tilde{φ}$	vehicle roll angle relative to the road surface
$\hat{χ}$	relative orientation of vehicle speed $V$ with respect to the tangent to the track spine curve
$Ω_{B}$	ribbon frame angular velocity
$Ω_{x}$	relative torsion
$Ω_{y}$	normal curvature
$Ω_{z}$	geodesic curvature
$ω$	vehicle frame angular velocity

References

Siegler, B.; Deakin, A. Lap Time Simulation: Comparison of Steady State, Quasi-Static and Transient Racing Car Cornering Strategies. In Proceedings of the 2000 SAE Motorsports Engineering Conference and Exposition, Dearborn, MI, USA, 13 November 2000. [Google Scholar] [CrossRef]
Brayshaw, D.L.; Harrison, M.L. Use of numerical optimization to determine the effect of the roll stiffness distribution on race car performance. Proc. Inst. Mech. Eng. Part D J. Automob. Eng. 2005, 219, 1141–1151. [Google Scholar] [CrossRef]
Brayshaw, D.L.; Harrison, M.L. A quasi steady state approach to race car lap simulation in order to understand the effects of racing line and centre of gravity location. Proc. Inst. Mech. Eng. Part D J. Automob. Eng. 2005, 219, 725–739. [Google Scholar] [CrossRef]
Völkl, T.; Muehlmeier, M.; Winner, H. Extended Steady State Lap Time Simulation for Analyzing Transient Vehicle Behavior. SAE Int. J. Passeng. Cars–Mech. Syst. 2013, 6, 283–292. [Google Scholar] [CrossRef]
Cossalter, V.; Da Lio, M.; Lot, R.; Fabbri, L. A General Method for the Evaluation of Vehicle Manoeuvrability with Special Emphasis on Motorcycles. Veh. Syst. Dyn. 1999, 31, 113–135. [Google Scholar] [CrossRef]
Bobbo, S.; Cossalter, V.; Massaro, M.; Peretto, M. Application of the “Optimal Maneuver Method” for Enhancing Racing Motorcycle Performance. SAE Int. J. Passeng. Cars–Mech. Syst. 2009, 1, 1311–1318. [Google Scholar] [CrossRef]
Casanova, D.; Sharp, R.S.; Symonds, P. Minimum Time Manoeuvring: The Significance of Yaw Inertia. Veh. Syst. Dyn. 2000, 34, 77–115. [Google Scholar] [CrossRef]
Cossalter, V.; Lot, R.; Tavernini, D. Optimization of the centre of mass position of a racing motorcycle in dry and wet track by means of the “Optimal Maneuver Method”. In Proceedings of the 2013 IEEE International Conference on Mechatronics (ICM), Vicenza, Italy, 27 March 2013; pp. 412–417. [Google Scholar] [CrossRef]
Lot, R.; Evangelou, S.A. Lap Time Optimization of a Sports Series Hybrid Electric Vehicle. In Proceedings of the World Congress on Engineering, London, UK, 3–5 July 2013. [Google Scholar] [CrossRef]
Lot, R.; Biral, F. A Curvilinear Abscissa Approach for the Lap Time Optimization of Racing Vehicles. IFAC Proc. Vol. 2014, 47, 7559–7565. [Google Scholar] [CrossRef]
Limebeer, D.J.N.; Perantoni, G.; Rao, A.V. Optimal control of Formula One car energy recovery systems. Int. J. Control 2014, 87, 2065–2080. [Google Scholar] [CrossRef]
Perantoni, G.; Limebeer, D.J.N. Optimal control for a Formula One car with variable parameters. Veh. Syst. Dyn. 2014, 52, 653–678. [Google Scholar] [CrossRef]
Perantoni, G.; Limebeer, D.J.N. Optimal Control of a Formula One Car on a Three-Dimensional Track—Part 1: Track Modeling and Identification. J. Dyn. Syst. Meas. Control 2015, 137, 051018. [Google Scholar] [CrossRef]
Perantoni, G.; Limebeer, D.J.N. Optimal Control of a Formula One Car on a Three-Dimensional Track—Part 2: Optimal Control. J. Dyn. Syst. Meas. Control 2015, 137, 051019. [Google Scholar] [CrossRef]
Dal Bianco, N.; Lot, R.; Gadola, M. Minimum time optimal control simulation of a GP2 race car. J. Automob. Eng. 2018, 232, 1180–1195. [Google Scholar] [CrossRef]
Massaro, M.; Marconi, E. The effect of engine spin direction on the dynamics of powered two wheelers. Veh. Syst. Dyn. 2018, 56, 604–620. [Google Scholar] [CrossRef]
Leonelli, L.; Limebeer, D.J.N. Optimal control of a road racing motorcycle on a three-dimensional closed track. Veh. Syst. Dyn. 2020, 58, 1285–1309. [Google Scholar] [CrossRef]
Marconi, E.; Massaro, M. The Effect of Suspensions and Racetrack Three-Dimensionality on the Minimum Lap Time of Motorcycles. In Advances in Dynamics of Vehicles on Roads and Tracks; Klomp, M., Bruzelius, F., Nielsen, J., Hillemyr, A., Eds.; Springer International Publishing: Cham, Switzerland, 2020; pp. 1368–1377. [Google Scholar] [CrossRef]
Limebeer, D.J.N.; Bastin, M.; Warren, E.; Fensham, H.G. Optimal control of a NASCAR—Specification race car. Veh. Syst. Dyn. 2022, 61, 1210–1235. [Google Scholar] [CrossRef]
Veneri, M.; Massaro, M. A free-trajectory quasi-steady-state optimal-control method for minimum lap-time of race vehicles. Veh. Syst. Dyn. 2020, 58, 933–954. [Google Scholar] [CrossRef]
Lovato, S.; Massaro, M.; Limebeer, D.J.N. Curved-ribbon-based track modelling for minimum lap-time optimisation. Meccanica 2021, 56, 2139–2152. [Google Scholar] [CrossRef]
Lovato, S.; Massaro, M. A three-dimensional free-trajectory quasi-steady-state optimal-control method for minimum-lap-time of race vehicles. Veh. Syst. Dyn. 2022, 60, 1512–1530. [Google Scholar] [CrossRef]
Rowold, M.; Ögretmen, L.; Kasalowsky, U.; Lohmann, B. Online Time-Optimal Trajectory Planning on Three-Dimensional Race Tracks. In Proceedings of the 2023 IEEE Intelligent Vehicles Symposium (IV), Anchorage, AK, USA, 4–7 June 2023; pp. 1–8. [Google Scholar] [CrossRef]
Biniewicz, J.; Pyrz, M. A quasi-steady-state minimum lap time simulation of race motorcycles using experimental data. Veh. Syst. Dyn. 2024, 62, 372–394. [Google Scholar] [CrossRef]
Biral, F.; Lot, R. An interpretative model of g-g diagrams of racing motorcycle. In Proceedings of the 3rd ICMEM International Conference on Mechanical Engineering and Mechanics, Beijing, China, 21–23 October 2009. [Google Scholar]
Timings, J.P.; Cole, D.J. Vehicle trajectory linearization to enable efficient optimisation of the constant speed racing line. Veh. Syst. Dyn. 2012, 50, 883–901. [Google Scholar] [CrossRef]
Timings, J.P.; Cole, D.J. Minimum maneuver time calculation using convex optimization. J. Dyn. Syst. Meas. Control 2013, 135, 031015-1–031015-9. [Google Scholar] [CrossRef]
Timings, J.P.; Cole, D.J. Robust lap-time simulation. Proc. Inst. Mech. Eng. Part D J. Automob. Eng. 2014, 228, 1200–1216. [Google Scholar] [CrossRef]
Anderson, J.; Ayalew, B. Modelling minimum-time manoeuvering with global optimisation of local receding horizon control. Veh. Syst. Dyn. 2018, 56, 1508–1531. [Google Scholar] [CrossRef]
Kanchwala, H. Path Planning and Tracking of an Autonomous Car With High Fidelity Vehicle Dynamics Model and Human Driver Trajectories. In Proceedings of the 2019 IEEE 10th International Conference on Mechanical and Aerospace Engineering (ICMAE 2019), Brussels, Belgium, 22–25 July 2019; pp. 306–313. [Google Scholar] [CrossRef]
Yang, B.; Zhang, H.; Jiang, Z. Path Planning and Tracking Control for Automatic Driving Obstacle Avoidance. In Proceedings of the 2019 International Conference on Robotics, Intelligent Control and Artificial Intelligence (RICAI 2019), Shanghai, China, 20–22 September 2019; pp. 300–304. [Google Scholar] [CrossRef]
Tan, Z.; Wei, J.; Dai, N. Real-time Dynamic Trajectory Planning for Intelligent Vehicles Based on Quintic Polynomial. In Proceedings of the 2022 IEEE 21st International Conference on Ubiquitous Computing and Communications (IUCC/CIT/DSCI/SmartCNS), Chongqing, China, 19–21 December 2022; pp. 51–56. [Google Scholar] [CrossRef]
Pan, R.; Jie, L.; Zhao, X.; Wang, H.; Yang, J.; Song, J. Active Obstacle Avoidance Trajectory Planning for Vehicles Based on Obstacle Potential Field and MPC in V2P Scenario. Sensors 2023, 23, 3248. [Google Scholar] [CrossRef] [PubMed]
Zhang, H.; Yao, J.; Tian, S. Model Predictive Control for Trajectory Planning Considering Constraints on Vertical Load Variation. Electronics 2024, 13, 1488. [Google Scholar] [CrossRef]
Massaro, M.; Limebeer, D.J.N. Minimum-lap-time optimisation and simulation. Veh. Syst. Dyn. 2021, 59, 1069–1113. [Google Scholar] [CrossRef]
Limebeer, D.J.N.; Massaro, M. Dynamics and Optimal Control of Road Vehicles; Oxford University Press: Oxford, UK, 2018. [Google Scholar] [CrossRef]
Dal Bianco, N.; Bertolazzi, E.; Biral, F.; Massaro, M. Comparison of direct and indirect methods for minimum lap time optimal control problems. Veh. Syst. Dyn. 2019, 57, 665–696. [Google Scholar] [CrossRef]
Lovato, S.; Massaro, M. Three-dimensional fixed-trajectory approaches to the minimum-lap time of road vehicles. Veh. Syst. Dyn. 2022, 60, 3650–3667. [Google Scholar] [CrossRef]
Limebeer, D.J.N.; Warren, E. A review of road models for vehicular control. Veh. Syst. Dyn. 2023, 61, 1449–1475. [Google Scholar] [CrossRef]
Massaro, M.; Lovato, S.; Veneri, M. An optimal control approach to the computation of g-g diagrams. Veh. Syst. Dyn. 2024, 62, 448–462. [Google Scholar] [CrossRef]
Patterson, M.A.; Rao, A.V. GPOPS-II: A MATLAB Software for Solving Multiple-Phase Optimal Control Problems Using hp-Adaptive Gaussian Quadrature Collocation Methods and Sparse Nonlinear Programming. ACM Trans. Math. Soft. (TOMS) 2014, 41, 1–37. [Google Scholar] [CrossRef]
Weinstein, M.J.; Rao, A.V. Algorithm 984: ADiGator, a Toolbox for the Algorithmic Differentiation of Mathematical Functions in MATLAB Using Source Transformation via Operator Overloading. ACM Trans. Math. Soft. 2017, 44, 1–25. [Google Scholar] [CrossRef]
United States Geological Survey. Available online: www.usgs.gov (accessed on 6 April 2024).
GPS Visualizer. Available online: www.gpsvisualizer.com (accessed on 6 April 2024).
Snyder, J.P. Map Projections: A Working Manual; U.S. Government Printing Office: Washington, DC, USA, 1987.
Nolan “Indy” Lamkin. Available online: www.nolanlamkinracing.com (accessed on 6 April 2024).

Figure 1. The ribbon frame

O_{t n m}

and the vehicle frame

A_{\hat{x} \hat{y} \hat{z}}

Figure 1. The ribbon frame

O_{t n m}

and the vehicle frame

A_{\hat{x} \hat{y} \hat{z}}

Figure 2. Vehicle positioning on the three-dimensional road.

Figure 3. Generalized rigid-body motorcycle model.

F_{x},

F_{y},

and

N

are the longitudinal, lateral, and vertical forces acting on the front (subscript

f

) and the rear (subscript

r

) tire.

F_{d}

is the drag force. The center of gravity and center of pressure are indicated by the acronyms CoG and CoP, respectively.

Figure 3. Generalized rigid-body motorcycle model.

F_{x},

F_{y},

and

N

are the longitudinal, lateral, and vertical forces acting on the front (subscript

f

) and the rear (subscript

r

) tire.

F_{d}

is the drag force. The center of gravity and center of pressure are indicated by the acronyms CoG and CoP, respectively.

Figure 4. Vehicle performance envelope: (a) exemplary g-g diagram of motorcycle at selected speed

V

modified by road camber

ϕ

, road inclination

μ

and road normal curvature variations represented by the term

{\hat{ω}}_{y} V

; (b) hypersurfaces of the adhesion radius generated for various values of apparent gravity

\tilde{g}

Figure 4. Vehicle performance envelope: (a) exemplary g-g diagram of motorcycle at selected speed

V

modified by road camber

ϕ

, road inclination

μ

and road normal curvature variations represented by the term

{\hat{ω}}_{y} V

; (b) hypersurfaces of the adhesion radius generated for various values of apparent gravity

\tilde{g}

Figure 5. Enlarged view of turn 12. Measured points on track edges (circles), with the reconstructed track boundaries using penalty factor

w_{θ}

equal to

4 \times 10^{4}

(thick line in light blue) and

4 \times 10^{6}

(thin line in purple). The dashed lines symbolize the calculated spine curves.

Figure 5. Enlarged view of turn 12. Measured points on track edges (circles), with the reconstructed track boundaries using penalty factor

w_{θ}

equal to

4 \times 10^{4}

(thick line in light blue) and

4 \times 10^{6}

(thin line in purple). The dashed lines symbolize the calculated spine curves.

Figure 6. Track reconstruction results: (a) plain view of Road Atlanta circuit, (b) track curvatures.

Figure 7. Track boundaries altitude: (a) left edge, (b) right edge.

Figure 8. Height difference between track boundaries.

Figure 9. Comparison of the Euler angles: (a) road camber

ϕ

(b) road inclination

μ

Figure 9. Comparison of the Euler angles: (a) road camber

ϕ

(b) road inclination

μ

Figure 10. Discrepancies in the road geometry between 3D-local and 3D-global cases: (a) track sector between turns 5 and 8, (b) turn 11 and chicane composed of corners 10a and 10b. The z-axis is exaggerated by a factor of three to emphasize the differences between the models.

Figure 11. Track width in particular road models.

Figure 12. Plots related to the vehicle powertrain: (a) engine torque and power measured at the rear wheel, (b) driving force for gears 2–6 and total resistance force.

Figure 13. Velocity profiles (top) and time difference (bottom) as a function of elapsed distance: (a) comparison between experimental data and 3D-local case, (b) comparison between GPS speed and vehicle speed computed in 3D-global and 2D cases.

Figure 14. Longitudinal acceleration (top) and lateral acceleration (bottom). Experimental data compared with (a) 3D-local simulation, (b) 3D-global and 2D cases.

Figure 15. g-g diagrams. Experimentally measured accelerations compared with (a) 3D-local simulation, (b) 3D-global simulation, (c) 2D road case.

Figure 16. Comparison of (a) vehicle lateral position, (b) radius of curvature modulus.

Figure 17. Comparison of the optimal trajectories; enlarged view of (a) turn 11, (b) chicane 10a–10b.

Figure 18. From top to bottom: vertical force acting on the front wheel, front suspension deflection and (in blue) wheel speed, throttle position.

Figure 19. Total vertical force depending on the adopted road model.

Figure 20. Enlarged view of the track section including turn 5 and the straight between turns 5 and 6 (left). The right side of the figure shows (from top to bottom) vehicle speed, longitudinal acceleration, lateral acceleration, radius of curvature modulus, throttle position, front brake pressure, engine speed, and gear selected.

Table 1. Supersport motorcycle parameters.

Symbol	Description	Value
$g$	Standard gravity	9.81 m/s²
$m$	Total mass (vehicle + rider)	255 kg
$w$	Wheelbase	1.40 m
$b$	Longitudinal position of CoG	0.69 m
$h$	Height of CoG	0.66 m
$h_{p}$	Height of CoP	0.66 m
$f_{w}$	Rolling resistance coefficient	0.02
$C_{d} A$	Drag area coefficient	0.28 m²
$ρ_{a}$	Air density	1.20 kg/m³
$μ_{x}$	Longitudinal friction coefficient	1.18
$μ_{y}$	Lateral friction coefficient	1.13
$i_{p}$	Primary transmission ratio	2.07
$i_{g_{2}}$	Second gear ratio	2.00
$i_{g_{3}}$	Third gear ratio	1.67
$i_{g_{4}}$	Fourth gear ratio	1.44
$i_{g_{5}}$	Fifth gear ratio	1.29
$i_{g_{6}}$	Sixth gear ratio	1.15
$i_{s}$	Secondary transmission ratio	2.88
$r_{t}$	Tire rolling radius	0.330 m

Table 2. Values of the coefficients

β_{0}^{r}

and

β_{1}^{r}

used in simulations.

Table 2. Values of the coefficients

β_{0}^{r}

and

β_{1}^{r}

used in simulations.

Symbol	Value
$β_{0}^{j_{x_{a}}}$	15.440
$β_{1}^{j_{x_{a}}}$	−0.148
$β_{0}^{j_{x_{d}}}$	−28.440
$β_{1}^{j_{x_{d}}}$	0.148
$β_{0}^{j_{y}}$	32.559
$β_{1}^{j_{y}}$	−0.378

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Biniewicz, J.; Pyrz, M. Influence of the Road Model on the Optimal Maneuver of a Racing Motorcycle. Appl. Sci. 2024, 14, 4006. https://doi.org/10.3390/app14104006

AMA Style

Biniewicz J, Pyrz M. Influence of the Road Model on the Optimal Maneuver of a Racing Motorcycle. Applied Sciences. 2024; 14(10):4006. https://doi.org/10.3390/app14104006

Chicago/Turabian Style

Biniewicz, Jan, and Mariusz Pyrz. 2024. "Influence of the Road Model on the Optimal Maneuver of a Racing Motorcycle" Applied Sciences 14, no. 10: 4006. https://doi.org/10.3390/app14104006

APA Style

Biniewicz, J., & Pyrz, M. (2024). Influence of the Road Model on the Optimal Maneuver of a Racing Motorcycle. Applied Sciences, 14(10), 4006. https://doi.org/10.3390/app14104006

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Influence of the Road Model on the Optimal Maneuver of a Racing Motorcycle

Abstract

1. Introduction

2. Road Modeling

2.1. Track Identification OCP

2.2. Acquiring and Preparing the Input Dataset

3. Three-Dimensional Free-Trajectory Quasi-Steady-State MLTS

3.1. Vehicle Positioning on the Three-Dimensional Road and Formulation of the Minimum Lap Time Problem

3.2. Vehicle Model and Related Vehicle Performance Constraint

4. Results of Track Reconstruction

Comparison between Results Obtained Using Local and Global DEMs

5. Optimal Maneuver Simulation Results

5.1. Comparison of Simulation Results on Adopted Road Models

5.2. Analysis of Tire Vertical Load

5.3. Detailed Analysis of Selected Track Segment

6. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

Abbreviations

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI