A Novel Simulation-Based Optimization Method for Autonomous Vehicle Path Tracking with Urban Driving Application

: Autonomous driving technology heavily depends on accurate and smooth path tracking. Facing complex urban driving scenarios, developing a suite of high-performance and robust parameters for controllers becomes imperative. This paper proposes a stochastic simulation-based optimization model for optimizing the Proportional–Integral–Diﬀerential (PID) controller parameters, with tracking accuracy and smoothness as bi-objectives, and solves it using a domination-measure-based eﬃcient global optimization (DMEGO) algorithm. In this model, the tracking accuracy and smoothness are indexed by the normalized dynamic time warping (NDTW) and the mean absolute lateral acceleration (MALA), respectively. In addition, we execute the PID controller in a realistic simulation environment using a CARLA simulator, which consider various city scenes, diverse routes, diﬀerent vehicle types, road slopes, etc., to provide a comprehensive and reliable evaluation for the designed PID controller. In the DMEGO method, each solution undergoes evaluation using a ﬁxed number of costly simulations. Then, utilizing the solutions and their estimated bi-objective values, two surrogate models for the bi-objectives are constructed using the Gaussian process (GP) model. The preliminary nondominated solutions can be obtained by optimizing the two surrogate models. Finally, a novel performance metric known as the domination measure is employed to evaluate the quality of each solution. This metric is then integrated with the crowding distance to selectively retain a candidate solution exhibiting superior performance and good diversity for the next iteration. In our numerical experiments, we ﬁrst test the DMEGO algorithm against three other counterparts using a stochastic FON benchmark. The proposed approach is then employed to optimize the PID parameters considering the complexity and uncertainty of urban traﬃc. The numerical results demonstrate that the nondominated solutions obtained by DMEGO exhibit excellent performance in terms of tracking accuracy and smoothness under limited simulation budgets. Overall, the proposed approach may be a viable tool for solving multi-objective simulation-based optimization problem under uncertainties.


Introduction
In recent years, the self-driving car industry has witnessed a notable proliferation of commercial ventures, generating significant interest in autonomous driving [1][2][3][4].The existing advanced automatic driving systems operate in a modular framework encompassing various components.Apart from the hardware aspect, these systems comprise a global navigation module [5], a perception module [6,7], a local planning module [8], and a control module [9][10][11].This study concentrates on the control module, an integral component that ensures safety, comfort, and energy efficiency in autonomous vehicles via a meticulously crafted controller that precisely adheres to a predefined trajectory [12][13][14].
The design of a controller for autonomous driving must meet two essential requirements.Firstly, it should be capable of real-time operation, necessitating a simple structure that enables rapid response to feedback [15].Secondly, the controller must exhibit high accuracy [16], ensuring that autonomous vehicles can diligently adhere to the intended trajectory while keeping their tracking errors within a permissible range.Nowadays, PID controllers have emerged as a widely adopted control mechanism in autonomous driving due to their simplicity and real-time nature.With the proportional, integral, and differential components, PID controllers enable vehicles to respond rapidly to environmental changes, providing excellent lateral and longitudinal control to ensure a smooth and safe driving experience.Despite the above advantages of PID controllers, their effectiveness in complex urban scenarios within the real world heavily relies on appropriate parameter selection.It is crucial to note that suboptimal parameter choices can result in inadequate control performance or even system failure.Figure 1 illustrates the utilization of manually set PID parameters for tracking the planning path depicted in Figure 1a.However, it becomes evident that during the driving process, the vehicle experiences undesirable shaking, resulting in a substantial deviation between the tracked trajectory and the planning path.Consequently, such erratic behavior compromises the safety of the driving operation.Furthermore, the vehicle's oscillations contribute to a rapid escalation in absolute lateral acceleration, significantly impacting the ride quality.Overall, the complexity and uncertainties inherent in real-world traffic scenarios make determining the PID parameters for self-driving vehicles still a challenging task.Some existing works have a empted to achieve a well-designed tracking controller for autonomous vehicles by determining the optimal parameters of the PID controller.Manual PID tuning involves adjusting the controller gains based on the observed system response and the tuner's experience and knowledge.Notably, this manual process might be time-consuming and require a deep understanding of the autonomous vehicle system.With the advancement of technology, more and more researchers try to determine the parameters of PID controllers automatically using intelligent optimization algorithms.These advanced techniques exhibit superior precision and alleviate the labor-intensive nature inherent in manual tuning processes.Consequently, they are highly preferred when a mathematical model of the system is accessible or can be derived, thereby facilitating their practical implementation in academic and industrial se ings.Zhao et al. [17] introduced a genetic algorithm (GA)-based PID controller for trajectory tracking in self-driving vehicles.Utilizing a Fractional Order PID (FOPID) controller, Al-Mayyahi et al. [18] focused on controlling a non-holonomic autonomous ground vehicle's movement to follow a specific reference path, with its parameters refined using particle swarm optimization (PSO).A soft computing evolutionary method, specifically the GA, was applied by Abajo et al. [19] to find the optimal PID regulator se ings for varying routes.Gao et al. [20] proposed a segmented fuzzy PID controller enhanced by a hybrid PCAG algorithm that combines adaptive PSO and GA.Hajjami et al. [21] utilized the bu erfly optimization algorithm (BOA) to establish an optimal PID controller for improving the lateral dynamics of AVs.Loucif et al. [22] determined the optimal parameters for a robot manipulator controller by employing an innovative evolutionary algorithm, namely the whale optimization algorithm (WOA).Ma'ani and Nazaruddin [23] designed a longitudinal controller applying the flower pollination algorithm (FPA) for controller parameter optimization.The proposed longitudinal controller's efficacy in maintaining the desired speed while navigating non-straight paths was evaluated via a simulation study using the CARLA simulator.
Although the studies mentioned above have demonstrated promising outcomes within controlled laboratory se ings, they still suffer from some notable limitations that hinder their practical applicability in path tracking for autonomous vehicles.Firstly, most existing work focuses on optimizing the PID parameters for a specific pre-planned path.While this approach may ensure satisfactory performance on the designated trajectory, it may fail to guarantee comparable performance when generalized to different routes or scenarios.Secondly, most existing studies on PID parameter optimization rely on MATLAB or mathematical simulations, where researchers often employ simplified vehicle dynamics and overlook external environmental factors.As a result, a significant disparity exists between their findings and real-world scenarios.Nevertheless, autonomous vehicles are highly intricate dynamic systems wherein numerous external environmental factors significantly influence path-tracking behavior.Consequently, it is complex and challenging to optimize PID parameters that well fit the environmental perturbations caused by the complexity of vehicle dynamics, road slope, weather, and vehicle type.Lastly, current works mainly a empt to minimize the tracking error in path tracking using a single-objective optimization technique.However, achieving a balance between smoothness and tracking error control is challenging for the path-tracking problem.In summary, Table 1 presents a comprehensive comparison between the relative literature and our work, offering a detailed overview of the key distinctions and contributions.To our knowledge, there is a research gap in the literature regarding the design of a robust PID controller that can effectively mitigate environmental perturbations while ensuring both comfort and accuracy in control.
For the above problems, this study developed a realistic and comprehensive simulation platform to design a robust PID controller for AV path tracking.However, the high simulation resolution and comprehensive environmental modeling also result in increased computational demand.When dealing with computationally intensive and timeconsuming simulations, the heuristic evolutionary algorithms are not applicable because of their low efficiency in exploring the search space with a large number of function evaluations.The simulation-based optimization (SO) method provides a promising tool for addressing the complexities of global optimization problems that involve intractable functions, requiring expensive evaluations.SO techniques are wide-ranging, including random search, stochastic approximation, and especially Bayesian Optimization (BO) [24].Among these, BO is widely recognized as the dominant algorithm, with applications spanning drug molecular design [25,26], material discovery [27], and automatic hyperparameter tuning [28,29].BO originated with [30] and gained further prominence with the introduction of Efficient Global Optimization (EGO) by [31].BO stands out by reducing the number of function evaluations by utilizing a surrogate model and sampling methods guided by elaborate selection criteria.These selection criteria can be categorized into two types: sequential selection, which suggests a single optimal sample for evaluation in each iteration, and batch selection, which selects multiple samples to be evaluated in parallel during each iteration.More specifically, sequential selection has garnered significant scholarly interest, as evidenced by its exploration in the works of [32][33][34].While this method is associated with improved performance due to fewer evaluations, it involves a trade-off of increased iterations and reduced convergence speed.In contrast, batch selection sacrifices some efficiency to accommodate parallel evaluations, aligning with current trends toward time-efficient experimental practices [35,36].For the scope of this paper, our emphasis is placed on the sequential selection scheme.
The optimization of computationally expensive problems using the BO method has become commonplace.However, extending BO to address multi-objective optimization problems is not straightforward.ParEGO [37] stands out as a pioneering and widely used multi-objective Bayesian Optimization (MOBO) algorithm.It transforms the multi-objective problem into a single-objective problem using random scalarization and selects a sample based on the maximum expected improvement (EI).EHI [38] extends ParEGO by introducing the expected improvement of the hypervolume and provides a direct computation procedure for the integral expression.Similarly, SUR [39] introduces an enhanced MOBO algorithm that leverages stepwise uncertainty reduction principles.In this approach, the optimization process is characterized by progressively minimizing the volume of excursion sets below the currently known optimal solutions.In order to improve the efficiency of algorithms, complementary research aims to enhance the reliability of surrogate models.He et al. [40] introduce a bi-objective efficient global optimization (BOEGO) algorithm that utilizes regressive kriging models to predict bi-objectives under stochastic disturbances.However, existing studies have not considered maintaining diversity in the performance space within MOBO despite its clear significance in practical applications.
The remainder of this paper is organized as follows.In Section 2, we develop a realistic simulation environment for AV path tracking, which includes a vehicle dynamic model, longitudinal and lateral PID controllers, external environmental perturbations, and two evaluation indices describing accuracy and smoothness, respectively.Section 3 presents the DMEGO method.In Section 4, we test a synthetic function and study a PID parameter optimization problem.Finally, some insightful conclusions are reached.

Model
This section first discusses the vehicle dynamic model and introduces the longitudinal and lateral PID controller for AV path tracking.Subsequently, a realistic simulation environment based on the CARLA simulator is developed to execute path-tracking tasks while accounting for external environment perturbations, such as weather conditions, road slopes, and various vehicle types.Finally, we present two evaluation indicators of tracking accuracy and smoothness and formulate the PID controller parameter optimization model under uncertainties.

Vehicle Dynamic Model
Vehicle dynamics models are sophisticated tools that accurately simulate various vehicular actions, such as acceleration, braking, steering, and suspension responses.These models integrate critical parameters like vehicle mass, inertia, tire properties, and aerodynamic forces to predict the vehicle behavior under various conditions.Within the established vehicle dynamics models, the bicycle model is often adopted due to its simplified yet effective representation of AVs.This model conceptualizes the AV as a singular mass point, linking the front and rear wheel assemblies with a hypothetical axle.It provides a kinematic approximation that is sufficiently accurate for the synthesis of advanced control algorithms, which are critical for precise path tracking.
This paper employs the two-degrees-of-freedom bicycle vehicle model, as depicted in Figure 2. The control inputs for the entire model can be simplified as ( , ), where represents the vehicle's acceleration.Positive acceleration is a ained by engaging the accelerator, while negative acceleration is achieved by applying the brake pedal.Furthermore, represents the angle formed between the wheels and the body, serving as a reliable indicator of the steering wheel's turning angle.The chassis is aligned with the vehicle's center of gravity (CG).We define the state vector [ , , , , ] at the CG, with and as the CG's global coordinates, as the roll angle, and and as the immediate coordinates of the vehicle's body.O represents the instantaneous center of rotation (ICR), which is a notional point or axis around which a vehicle rotates at any given moment.Herein, we present the dynamic differential equations governing the vehicle's behavior.

̈= ̇× ̇+
(1) where ̇ and ̇ represent the longitudinal and lateral velocities of the vehicle at the CG, respectively.Furthermore, denotes the vehicle mass, stands for the vehicle rotational inertia, represents the vehicle lateral offset angle, and and indicate the distances from the CG to the respective front and rear axles.The tire forces acting laterally on both the front and rear tires are expressed as , , which can be determined using the linear tire model described by Formula (7).In this model, is referred to as the tire cornering stiffness.
represents the deflection angle of the tire, which refers to the deviation between the direction of the tire and its velocity vector.For a more comprehensive understanding of vehicle dynamics, please refer to [41].

PID Controller
A PID controller based on traditional feedback control eliminates the need for a derived model and allows for quick and straightforward deployment using computed tracking errors.To achieve precise path tracking, both the longitudinal and lateral PID controllers are employed in this study.The primary role of the longitudinal PID controller is to regulate the accelerator and brake pedals, thereby maintaining a constant speed at the target speed.Simultaneously, the lateral controller governs the steering wheel angle, ensuring the vehicle can traverse a predetermined path.
As illustrated in Figure 3, the longitudinal controller exclusively considers the vehicle's current speed.The longitudinal tracking error (i.e., ) is defined as the difference between the target speed and the current speed.Moreover, the lateral controller considers both the direction vector of the vehicle speed and the direction vector of the waypoint.The angle between these vectors serves as the longitudinal tracking error (i.e., ).Based on the above principles, the longitudinal and lateral PID controllers can be designed as follows: where , , , , and , are the parameters, respectively, corresponding to the proportion, integration, and differentiation components of the longitudinal PID controller.Likewise, in the lateral PID controller, , , , , and , play analogous roles, signifying the proportion, integration, and differentiation terms.and denote the outputs of the longitudinal and lateral controllers, i.e., the opening of the gas or brake pedal and the steering wheel angle, respectively, and both are normalized to [−1, 1].Due to the feedback control mechanism, the PID controller can promptly respond to the tracking error and exhibit stable driving behavior if all the parameters have been well tuned.
Nevertheless, determining an optimal set of parameters poses a significant challenge, and the remaining sections of this paper are dedicated to addressing this issue.

Simulation Platform
This study develops a comprehensive simulation environment for AV path tracking, considering various environmental perturbations.Specifically, this study employs CARLA [42], an open-source high-fidelity AV simulation platform that enables the design, testing, and validation of AV algorithms in realistic urban driving environments.The proposed simulation environment is shown in Figure 4a, which is designed to utilize a client/server (CS) architectural framework and incorporates external vehicular control hardware.In this simulation environment, CARLA, operating as the server, provides 3D static city scenes, realistic vehicle dynamic models, and various onboard sensors.The main interface for the client is Python scripts.By using a Python API, it is possible to select maps, change the weather, spawn various vehicles, and design vehicular control algorithms.Additionally, we have introduced Logitech G29's virtual driving hardware [43] to observe specific vehicle operations visually.In this paper, CARLA version 0.9.13 is utilized, and Figure 4b illustrates the AV pathtracking flow chart in our simulation environment.First of all, our simulation environment initializes the 3D urban map, spawns the autonomous vehicle, deploys onboard sensors, and defines the origin and destination (OD).Following these initial setup steps, the global path planning module leverages the A-star algorithm [44] to generate the global planning route that connects the specified origin and destination point.The A-star algorithm, also known as A*, is a popular pathfinding algorithm used in various applications such as AI for games, robotics, and autonomous driving systems.It is used to efficiently find the shortest path between two points in a graph-or grid-based environment, taking into account the obstacles and costs associated with different paths.In CARLA, the A-star algorithm can be employed to calculate optimal paths for autonomous vehicles to navigate from their current location to a desired destination.Within the control loop, the Global Navigation Satellite System (GNSS) is employed to accurately determine the vehicle's position, and the Inertial Measurement Unit (IMU) supplies essential data regarding the vehicle's speed and a itude.Simultaneously, by utilizing the vehicle's state information, including position, speed, and heading angle, the PID controller generates control signals for the vehicle, i.e., the thro le and steering wheel angle commands, which dictate the desired actions to be executed.It is noted that the proposed environment encompasses a variety of perturbations, including city scenes, planning routes, weather conditions, road slopes, and various vehicle types.A further depiction of these perturbations is provided in Table 2.As illustrated in Table 2, various environmental perturbations are introduced to rigorously test the capabilities of the AV path-tracking system.Our study selected three distinct urban cities in the CARLA simulator.These urban locales encompass Town 1, a European-style town composed of solely one-line roads and T junctions.Moving on, Town 3 emerges as a more expansive urban landscape, mirroring the a ributes commonly found in bustling downtown areas.Lastly, Town 10 hosts many junction arrangements with various lane markings, crossings, and signal types to challenge the autonomous vehicles in negotiating with other traffic.We curated 150 intricately designed routes in these three

Environmental Perturbations Examples
Town Routes Road slope

Weather
Vehicle type cities, encompassing various navigational complexities such as straight paths, turning maneuvers, and lane-changing scenarios.In addition to accounting for route diversity, our study also incorporates three distinct road types: flat terrain, uphill gradients, and downhill slopes.Recognizing weather's impact on road surface adhesion coefficients and vehicle wind resistance, we have incorporated four distinct weather categories into our experimental framework: sunny, light rain, moderate rain, and heavy rain.Lastly, acknowledging that the vehicle type significantly influences its maneuvering capabilities, we provide a selection comprising an SUV, a compact sedan, and a mini-car for experimentation purposes.To the best of our knowledge, this paper is the first to synthesize multiple environmental perturbations for an in-depth examination of AV path tracking.The procedure of AV path tracking is shown in Algorithm 1.
Algorithm 1: PID control for AV path tracking.Input: , , , , , , , , , , , #parameters of longitudinal and lateral PID controller , , , , # configuration for PID control algorithm Output: the global planning route , the historical track point set , and the historical lateral acceleration set 1.: Initialize the environment, select the town, vehicle type, OD pair, and weather, and set = { }, = {} 2.: Utilize A-star algorithm to generate global planning route = , , , , , , , … , , , , where each point is a reference point, and the distance of adjacent reference points is set as obtain the vehicle position , , , using GNSS, then ← ∪ , , , obtain the vehicle velocity, acceleration, and heading angle using IMU 7.
execute the controller output and update the simulation 11.

End while
In Algorithm 1, the PID control algorithm is implemented on the CARLA simulation platform for autonomous vehicle trajectory tracking.The process involved three main steps: environment initialization, global path planning, and PID control.Initially, the CARLA environment is set up, encompassing the selection of towns, weather conditions, and vehicle types.Subsequently, we define the origin and destination (OD) for the selfdriving trip and utilize the A-star algorithm to generate a global planning route = , , , , , , , … , , , , wherein the distance between successive waypoints in is denoted as .In this study, we set to 2 m.Finally, we employ longitudinal and lateral PID controllers to track the global path .The parameters of the two controllers are denoted as , , , , , for the longitudinal controller, and , , , , , for the lateral controller.Trajectory tracking is achieved by sequentially traversing the reference points in .If the vehicle's distance from the target point (denoted as ) is below the predefined threshold (set to 1.5 m in this study), the vehicle's target point is updated to the subsequent point in .This iterative process continues until either the distance between the vehicle and the destination becomes less than (set to 0.2 m in this study) or the number of iterations exceeds the maximum value .At that point, the PID control is terminated.The outputs of Algorithm 1 consist of the global planned route, the vehicle's historical track point set, and the vehicle's historical lateral acceleration set.In the subsequent section, we will utilize these results to evaluate the accuracy and smoothness for AV trajectory tracking.

PID Parameter Optimization Model
This paper a empts to achieve a trade-off between tracking accuracy and smoothness in AV path tracking.To this end, two corresponding evaluation indicators are introduced below.

Normalized Dynamic Time Warping
In AV path tracking, the tracking accuracy can be measured by the similarity between the globally planned path and the actual driving trajectory controlled by a PID algorithm.Trajectory similarity measures have garnered significant a ention in various domains, leading to the development of numerous measurement techniques in the literature [45,46].One of the most straightforward methods to measure the similarity between two trajectories is by calculating the spatial distance between corresponding points, such as the initial two points of each trajectory, the subsequent two points, and so forth.Subsequently, more sophisticated measures a empt to compare trajectories in more intricate ways.Dynamic time warping (DTW) stands as a classical dynamic programming algorithm that has been successfully utilized for measuring the similarity between trajectories [47,48].In this paper, DTW is computed by aligning elements of a reference trajectory (i.e., the globally planned path ) and a query trajectory (i.e., the actual driving trajectory ) while preserving their order.Specifically, DTW finds the optimal ordered alignment between the two trajectories by minimizing the cumulative cost: ( , ) = min ∈ , ( , )∈ (10) where is some distance function mapping pairs of elements from the two trajectories to real non-negative numbers and The principle of NDTW can be expressed in Figure 5, where the black lines in Figure 5a,b represent the point where the two tracks match and align.NDTW has many desirable properties for evaluating actual driving trajectories.It measures the similarity between the entirety of two trajectories, softly penalizing deviations., , , , , , , , , , , multiple independent simulation evaluations will yield varying and values owing to the perturbations in the simulation environment (e.g., randomly selected towns, routes, and weather conditions, etc.).Accordingly, the two objectives are represented as expectations relative to and , and the PID parameter optimization model is formulated as It is evident that model (13) seeks to simultaneously minimize the expectations related to or .In the next section, the DMEGO algorithm will be utilized to solve this model.

Methodology
We first provide a detailed introduction to the EGO algorithm.Readers already acquainted with the Gaussian process and expected improvement may proceed directly to Section 3.2.In Section 3.2, we propose a domination-measure-based efficient global optimization algorithm (DMEGO) for a multi-objective simulation optimization problem.

Efficient Global Optimization
As a prelude, a powerful framework namely EGO was proposed for single-objective global optimization problems with implicit functions that are expensive to evaluate.With the help of a surrogate model and acquisition function, the framework consists of three pivotal steps: (1) Utilizing a cost-effective surrogate model to estimate the implicit function based on solutions evaluated using true objective functions or simulations; (2) Designing an acquisition function that guides the optimization process toward promising solutions with improved performance and reduced uncertainty; (3) Selecting a candidate solution via the optimization of the acquisition function, evaluating it by the true objective function, and updating the surrogate model.

Gaussian Process
A Gaussian process (GP) [51] is a statistical approach to model probability distributions over functions.It is concisely defined by a prior mean function denoted as (•) and a covariance function represented as (•,•).As a widely favored surrogate model, the GP provides an estimation of the real objective function using a set of previously evaluated solutions.Let : = [ , … , ] be a finite collection of evaluated solutions, where ∈ ℝ , = 1, … , , and their corresponding objective function values are : , where = ( ).Thus, the GP prior is expressed with a multivariate normal distribution: : ~ : , ( : , : ) Consider : to be the vector containing the mean function values, with each corresponding to ( ).Let ( : , : ) be the covariance matrix, and the element at position ( , ) is denoted by , , quantifying the covariance between and .When dealing with situations where the objective function's form is entirely unknown, a typical approach involves utilizing a zero or constant as the prior mean.For the covariance function, the power exponential kernel is often preferred, defined by the following expression: In Equation ( 15), − represents the Euclidean distance between and , while and ℓ are parameters of the power exponential kernel.The analytical tractability of the GP is one of the factors that contributes to its popularity.As indicated in Equation (14), the true objective function values (i.e., : ) follow a joint Gaussian distribution based on the GP prior.Similarly, for a new sample , the estimate of the objective function ( ) is also expected to follow a jointly Gaussian distribution with the GP prior, which can be mathematically expressed as follows: where We can consider the solutions : and their true objective function values : as historical observations, while and ( ) represent a newly sampled solution and its estimate for an objective function, respectively.By conditioning on them, we can restrict the distribution of possible functions predicted by the GP, which is referred to as the GP posterior.When there are noisy observations (e.g., Gaussian noise with variance ), the predicted distribution-based observations follow a Gaussian distribution characterized by a posterior mean value ( ) and a posterior variance value ( ): ( )| : , , … , ~ ( ), ( ) ( ) = ( ) + ( , : )[ ( : , : ) + ] : − ( ) ( ) = ( , ) − ( , : )[ ( : , : ) + ] ( : , ) where is the identity matrix.The inverse operation on the matrix [ ( : , : ) + ] demands computational effort scaling with ( ), which can become computationally intensive when the observation count is substantial (for instance, upwards of 10 observations).Nevertheless, it is essential to highlight that in the context of SO, the number of observations typically remains comparatively limited, thereby mitigating concerns related to computational scalability.

Acquisition Function
The acquisition function utilizes the posterior mean and variance derived from the GP to determine the next solution for evaluation.From a minimization perspective, choosing a solution with a low posterior mean function value indicates exploitation, whereas opting for a solution with a high posterior variance signifies exploration.As such, the acquisition function encompasses both (•) and (•).To determine the next solution for evaluation, it is necessary to maximize the acquisition function.
Research on the BO algorithm has introduced several acquisition functions formulated using analytical expressions, including but not limited to the probability of improvement [52], expected improvement [53], and the upper confidence bound [54].In this paper, we employ the expected improvement (EI) acquisition function for the minimization problem, which is defined as follows: Where = ( ) where represents the smallest estimated objective function value among all evaluated solutions.The function (•) and (•) in Equation ( 23

Domination-Measure-Based Efficient Global Optimization
There exist several approaches to expand the scope of the EGO algorithm for addressing multi-objective simulation-based optimization problems.The most straightforward approach, referred to as ParEGO, entails transforming multi-objective values into a single objective value by employing a parameterized scalarizing weight vector.By selecting a different weight vector at each search iteration, it becomes feasible to progressively generate an approximation of the entire true Pareto front.The BOEGO algorithm extends the EGO to address multi-objective scenarios by integrating two regressive kriging models and employing an EI infill strategy.Nevertheless, based on our current knowledge, none of the prior methods consider diversity in the performance space, which can be critical for numerous multi-objective problems encountered in real-world scenarios.
Without losing generality, we assume a SO problem that seeks to minimize objective functions.Firstly, we use GP models to approximate true objective functions, respectively.Then, the key insight of our acquisition function is that instead of finding an -dimensional Pareto front corresponding to objectives, it finds a 2 -dimensional Pareto front where dimensions correspond to the performance objectives (each given by ( ) , = 1, … , In Equation ( 26), finding a solution within the feasible region that simultaneously minimizes all the objectives is extremely rare.Therefore, a compromised target is to identify all non-dominated solutions in based on their objective values.More specifically, the domination rule, Pareto optimal solutions, and Pareto front for problem (26)  set, and its representation in the objective space is referred to as the Pareto front.
Inspired by [55], this study introduces a unary and parameter-free performance metric called domination measure, to measure the quality of a solution under multiple evaluation criteria with no distinction between objectives.Specifically, this measure can be understood as the likelihood that a randomly chosen solution from the feasible region dominates that solution according to a predefined probability measure.So, ∀ ∈ , the domination measure for is defined as follows: where represents the collection of solutions that dominate , while (•) corresponds to the underlying Lebesgue measure, and (•) denotes the probability distribution on induced by (•).The indicator function { } equals 1 when event E is true and 0 otherwise.
[•] represents the expectation with respect to (•) .Several properties of the domination measure in can be straightforwardly established: (1) For any ∈ , it holds that 0 ≤ D( ) ≤ 1; (2) If ≺ , then D( ) ≤ D( ); (3) For any Pareto optimal solution * ∈ * , where * is the Pareto optimal solution set in , then D( * ) = 0. We present a numerical example demonstrating the advantages of using the domination measure to evaluate the solution quality in a multi-objective scenario.[56]  With no further elaboration here, the Pareto optimal set for Example 1 is ∩ ( [5,25] ∪ [60,85]), and it is highlighted with green and purple bold lines in Figure 7. Directly utilizing the two objective functions to locate the Pareto optimal set for Example 1 is difficult and unintuitive, so we take two methods to rigorously quantify the quality of a solution, respectively, i.e., the domination measure and the linear weighted sum method.Firstly, the domination measure ( ) and the linear combination of the biobjectives (i.e., w ( ) + w ( )) are calculated for each solution in the finite solution space, and the quality of a solution is easily obtained from Figure 8.Note that the la er method is parameter-based: without a loss of generality, we set five groups of weighting coefficients to represent the different preference degrees for the bi-objective functions.As shown in Figure 8, although both methods can strictly evaluate the solution quality using a deterministic scalar, the domination measure exhibits significant advantages in the following two aspects: (1) We can readily and directly identify the Pareto optimal set in Figure 8a by considering the fundamental principle that a smaller domination measure indicates a superior solution quality.In other words, solutions with a domination measure of zero precisely constitute the Pareto optimal set for Example 1.In Figure 8b, even a solution dominates another solution ; the linear weighted sum of is not always lower than that of .This inconsistency arises from the potential loss of information regarding the dominance relationship due to applying a weighting scheme to the bi-objective functions.Thus, the linear weighted sum method cannot effectively characterize the dominance relationship between solutions, which results in relatively good solutions rather than the Pareto optimal set.( 2) As a parameter-free performance metric, the domination measure does not depend on any fixed parameters.However, the linear weighted sum method requires first customizing the weighting coefficients, which would considerably impact evaluating the quality of a solution.Specifically, the larger the w in Figure 8b, the more the linear weighted sum method prefers the solutions in ∩ [60,85] and vice versa for solutions in ∩ [5,25].Consequently, the quality evaluation of solutions using the linear weighted sum method varies considerably under different weighting coefficients, making the selection of suitable coefficients challenging.In contrast to the linear weighted sum method, the domination measure does not necessitate parameter tuning and offers a rigorous and precise quantification of solution quality in multi-objective problems.

Example 1. Consider a bi-objective optimization problem
In this paper, we focus on a sequential selection strategy, which suggests a single optimal sample for evaluation in each iteration.However, the maximization problem for the acquisition function ( 26) aims to find a set of solutions rather than a single solution.Therefore, we propose a domination-measure-based efficient global optimization algorithm to solve multi-objective SO problems in a sequential selection manner.Specifically, the DMEGO algorithm iteratively carries out the following four main steps: (1) independently construct separate GP models for each objective function using evaluated observations; (2) employ the widely used NSGA-II algorithm to solve the maximization problem for the acquisition function, and obtain a set of promising candidate solutions; (3) utilize the domination measure and the crowding distance to selectively retain a candidate solution exhibiting superior performance and good diversity for the next iteration; (4) evaluate the ultimately selected candidate solution using true objective functions.In summary, the procedure of the DMEGO algorithm is shown in Algorithm 2.  , where * is the smallest integer that is greater than or equal to * .

9.
Reconstruct the set of candidate Pareto optimal solutions by = { | ( ) ≤ , ∈ }.Then, find a solution * in that has the largest crowding distance for .10.
Call the expensive evaluation function (•) with the solution * .Then, ← ∪ * , ← ∪ ( * ).11.By utilizing the fast nondominated sorting method with , the Pareto optimal solution set and its corresponding multi-objective function values can be obtained.

End while
The DMEGO algorithm is proposed to address the multi-objective simulation-based optimization problem.Algorithm 2 provides a general overview of the entire flow of DMEGO.However, when applying the DMEGO algorithm to solve the PID parameter optimization problem (referred to as problem ( 13)), it becomes necessary to clarify specific details.Initially, the decision variable dimension, , is set to 6, indicating three parameters each for the longitudinal and lateral controllers.Subsequently, there are only two objective functions: NDTW and MALA.Lastly, the true evaluation function (•) is composed of Algorithm 1 and Equations ( 11) and ( 12), representing the complex mapping relationship between the NTWR, MLSA, and PID control parameters.In Section 4.2, the DMEGO algorithm will be employed to solve the PID control parameter optimization problem.

Numerical Experiments
This section begins by conducting a comprehensive evaluation of the proposed DMEGO algorithm.Subsequently, the DMEGO algorithm is applied to tackle the PID controller parameter optimization problem.To test the effectiveness of the DMEGO algorithm and perform a comparative analysis with other competing algorithms, we employ the inverted generational distance (IGD) metric [57] and the diversity metric (Υ) [58] to quantify the convergence and diversity of DMEGO.More specifically, the IGD metric, defined in the objective space, is expressed as: denotes the average of these distances.In essence, Υ quantifies the evenness of the distribution of solutions in the solution space.Clearly, a smaller Υ value shows be er diversity.

A Benchmark Test Problem
This study employs a challenging benchmark test problem (30), essentially a stochastic adaptation of the FON function proposed by [59].The test problem (30) introduces stochastic perturbations into each objective function, thereby complicating the relationship between potential solutions and the expected outcomes of the bi-objective functions.Moreover, these stochastic variations also make the dominance relationships between any pair of solutions less apparent.In light of these challenges, we compare the performance of DMEGO against three other state-of-the-art algorithms, namely NSGA-II [60], ParEGO, and BOEGO, in addressing the problem (30).To facilitate a fair comparison among these four algorithms, we set a fixed budget of 200 function evaluations for this numerical case.Obviously, Figure 10 demonstrates the superiority of the nondominated solutions obtained using the DMEGO algorithm in terms of diversity and closeness to the TPF compared to the results obtained by the other three algorithms.Specifically, the nondominated solutions obtained using NSGA-II exhibit a considerable distance from the TPF, indicating the limited applicability of NSGA-II in addressing the problem (30) under restricted computational resources.As for ParEGO, although most nondominated solutions lie on the TPF, they exhibit sparsity and insufficient spread.BOEGO generates more nondominated solutions lying on the TPF.However, there are still some solutions that deviate from the TPF.Our proposed DMEGO method yields more nondominated solutions than the other algorithms, and all solutions exhibit excellent spread and precise alignment with the TPF.More detailed comparison results are presented in Table 3, which reveal that, despite the limited computational costs, DMEGO achieves a perfect Pareto front with superior convergence and diversity (evidenced by a smaller IGD value and a smaller Υ value).In conclusion, due to limited evaluation resources, the other three algorithms inefficiently utilize a large portion of samples by focusing on expanding a single region and encountering difficulties in escaping local minima.In contrast, our proposed DMEGO algorithm efficiently explores and identifies all regions comprising the TPF.Additionally, we visualize the iterative process in the DMEGO algorithm to validate its superiority and effectiveness.Figure 11a,b demonstrate that, after about 125 function evaluations, the IGD and HV using DMEGO remain relatively stable and notably outperform NSGA-II, ParEGO, and BOEGO.Finally, we employ K-fold cross-validation (KFCV) to evaluate the estimation errors of the two GP surrogate models for the bi-objectives during the iterations of DMEGO.The errors are quantified using two widely recognized indices, namely the mean absolute percent error (MAPE) and root mean square error (RMSE).Figure 11c illustrates that estimation errors are initially high with a small number of samples but gradually decrease as the sequential selection scheme generates more data points.After about 50 function evaluations, the RMSE and MAPE decrease to as low as approximately 0.09 and 5%, respectively.This demonstrates that the Gaussian process model consistently improves in accuracy and reliability throughout the iterative process, effectively guiding the DMEGO algorithm toward the appropriate exploration region.

PID Parameter Optimization Problem Using the DMEGO Method
In this section, we utilize the DMEGO method to solve the PID parameter optimization problem (13) for AV path tracking.The PID controller parameters have six dimensions, with three dedicated to the longitudinal controller and the remaining three for the lateral controller.NDTW and MALA serve as two objective functions for evaluating the tracking error and smoothness, and their true function evaluations require expensive simulations to collect tracking results (refer to Algorithm 1).Our experiment is tested using Intel As illustrated in Figure 12, the DMEGO algorithm successfully obtains three Pareto nondominated solutions (represented by red, green, and blue circles) after conducting 600 costly simulation evaluations (indicated by gray pentagrams).Specifically, in Figure 12a, it is evident that the optimization process gradually improves the NDTW metric from an initial value of approximately 40 to within 2. Additionally, the MALA metric is optimized from its initial value of approximately 3.5 to within 1. Figure 12b,c present a comprehensive depiction of the three nondominant solutions, referred to as Plan 1, Plan 2, and Plan 3. Notably, an inverse relationship is evident: as the NDTW value ascends, the MALA value descends.This observation underscores the trade-off between the tracking accuracy and smoothness, highlighting their mutually exclusive nature in AV path tracking.More detailed results for these three nondominated solutions are presented in Table 4.In addition to presenting the three optimized plans obtained by the DMEGO method, we also explore the influence of the PID parameters on the path-tracking performance of autonomous vehicles.After 600 simulation evaluations using the DMEGO algorithm, we select the last 50 solutions.Their distribution in the solution space is depicted in Figure 13.Each line in Figure 13 represents a set combination of six parameters of the lateral and longitudinal PID controller, and the line's color reflects the tracking performance (the more purple the color, the be er the performance).It can be seen from the figure that the vast majority of runs have bi-objective values below 2, which reflects the excellent stability of the DMEGO algorithm.Furthermore, as shown in Figure 13a, it is evident that higher values of , and , lead to a substantial improvement in tracking accuracy.This improvement can be a ributed to the fact that increasing the , parameter enables the vehicle to achieve faster acceleration or deceleration, thus effectively tracking the desired speed.The , parameter is utilized to mitigate the accumulation of speed deviations and ensure the vehicle maintains a stable speed.Thus, increasing the , parameter could enhance the speed stability.Likewise, as depicted in Figure 13b, it is evident that reducing the values of , and , considerably enhances the tracking smoothness.Nevertheless, raising the , parameter amplifies the control system's sensitivity to lateral errors, enabling the vehicle to rectify deviations promptly.It is important to note that excessively high values of , may result in overshooting and oscillation.Thus, the , value should not be set too high.On the other hand, it is crucial to avoid excessively high , values as they may introduce sensitivity into the noise, thereby diminishing the robustness of the PID controller.To strike a balance between the tracking accuracy and smoothness, we choose Plan 2 as the final scheme.Finally, we conduct comprehensive tests on the PID controller configured with Plan 2 in our proposed simulation environment, taking into account the presence of stochastic noise.This environment encompasses different city scenes, diverse routes, and various vehicle types, among other factors.Figure 14 presents a 3D visualization of our experimental scenarios.The green points depict the global route planned using the A-star algorithm for the autonomous vehicle, while the red points represent the target points that have been traversed.We select three towns, namely Town 1, Town 3, and Town 10, to create a diverse testing environment encompassing various urban environments and road characteristics.Specifically, Town 1 represents a relatively simple urban environment, making it suitable for initial testing and performance evaluation.Town 3 represents a moderately complex urban environment with straight and curved roads, varying road widths, and complex intersections, which may involve traffic lights, roundabouts, or multi-lane junctions.Town 10 represents a highly challenging and complex urban environment characterized by straight roads, sharp turns, complex junctions, and narrower streets.Furthermore, we have considered three distinct weather conditions: light rain, moderate rain, and heavy rain.The varying levels of ground wetness associated with these weather conditions result in variations in the adhesion coefficient of the road surface, posing a significant challenge to the path tracking of automated driving systems.In this paper, we have examined three different vehicle types in this study: van, car, and mini car.These vehicles exhibit varying body lengths and maneuvering stability levels, which can impact the performance of the PID controller.For a more intuitive demonstration of the experimental effects of Plan 2 discussed in this paper, a video version of the AV path tracking can be found in Supplementary Materials Video S1.As we can see from the video, our optimized PID control parameters can be adapted to various challenging scenarios, such as sharp turns, narrow roads, and even roundabouts.Notably, the controller performs well even in heavy rain and can be applied to a wide range of vehicles, highlighting the robustness of our optimized PID parameters against environmental perturbations in real-world conditions.Finally, several representative yet challenging routes (as shown in Figure 15) are selected to evaluate the path-tracking performance.Obviously, Plan 2 obtained using DMEGO achieves excellent results both in terms of tracking accuracy and smoothness.Specifically, across a variety of routes encompassing straight segments, sharp turns, and lane changes, Plan 2 consistently achieves a trajectory (i.e., the blue do ed line) that precisely aligns with the global planning path (i.e., the gray solid line), which ensures the safety of the autonomous vehicle.In addition, even when encountering sharp turns, Plan 2 effectively maintains the lateral acceleration below 30 m/s , thereby ensuring smooth operation and enhancing passenger comfort.In summary, the method proposed in this paper can obtain a set of robust and excellent PID parameters, which can be adapted to complex urban road environments.

Conclusions
In this paper, a novel simulation-based optimization method for AV path tracking was proposed.This method fully considers the complexity and uncertainty of the urban road environment and provides a robust solution for AV path tracking.This study develops a comprehensive simulation environment for AV path tracking, considering various environmental perturbations.Then, we present two evaluation indicators of tracking accuracy and smoothness and formulate the PID controller parameter optimization model under uncertainties.To solve this model, a simulation-based optimization algorithm named DMEGO is proposed and mainly includes three steps: (1) independently construct separate GP models for each objective function using evaluated observations; (2) employ the widely used NSGA-II algorithm to solve the maximization problem for the acquisition function and obtain a set of promising candidate solutions; (3) utilize the domination measure and the crowding distance to selectively retain a candidate solution exhibiting superior performance and good diversity for the next iteration.In the numerical experiments, the DMEGO algorithm is first tested using a stochastic FON function, and the numerical results validate the be er performance of DMEGO compared with NSGA-II, ParEGO, and BOEGO.Finally, the proposed approach is applied to optimize the PID control parameters considering the complexity and uncertainty of urban roads.The numerical results show that the nondominated solutions from DMEGO perform well in terms of the tracking accuracy and smoothness under limited simulation budgets.
Future enhancements and directions will focus on the following aspects: (1) incorporating a parallel computing mechanism to evaluate multiple solutions concurrently, thereby reducing the computational time even further; (2) adopting a more advanced machine learning or deep learning model as the surrogate model and implementing batch selection strategies to enhance the algorithm performance; (3) exploring the application of reinforcement learning techniques to achieve online dynamic tuning of PID parameters; and (4) although all experiments in this paper were conducted in a simulator, it is crucial to evaluate the practical performance in real-world scenarios in the future.

Figure 1 .
Figure 1.The result of the PID controller parameters by hand.Note: (a) represents the vehicle's tracking trajectory, and (b) represents the lateral acceleration of the vehicle.

Figure 3 .
Figure 3.The illustration of the longitudinal and lateral PID controllers.Note: (a) denotes the control of the vehicle's speed using the longitudinal controller, while (b) signifies the control of the vehicle's direction using the lateral controller.

Figure 4 .
Figure 4.The proposed simulation environment for AV path tracking.(a) Sequence diagram for information synchronization among server, client, and hardware.(b) The detailed flow chart for AV path tracking in CARLA.
= , …| | is a warping with = ( , ) ∈ [1: | |] × [1: | |] , respecting the step-size ( − ∈ (1,1), (1,0), (0,1) ), boundaries ( = (1,1), and | | = (| |, | |)) conditions.In AV path tracking, the cost function is the Euclidean distance between the coordinates of any two points in space.The ideal metric for evaluating the tracking error must be sensitive to the scale and density of the nodes along the actual driving trajectory.In order to address this, we propose a normalized DTW metric, which is defined as the sum of at least | | distance terms, by normalizing it by a factor of | | .This normalization process allows for a fair comparison and consistent evaluation across different scenarios.Directly derived from DTW, the NDTW metric is shown in Equation(11).

Figure 5 .
Figure 5. Illustration of two pairs of reference trajectory and query trajectory, and the optimal warping between them (black lines in (a,b)) obtained by computing NDTW.

( 13 )
s. t.Θ( , ζ) ⟼ { ( , ζ), ( , ζ)}In Equation(13), the feasible solution set is defined by = { ∶ ≤ ≤ }.The exogenous parameters ζ encompass both fixed (e.g., vehicle dynamic model and tire model) and stochastic (e.g., uncertainties in towns, routes, weather, and vehicle types) factors within the simulation environment.The simulation platform, denoted by Θ( , ζ) , processes inputs and ζ with the mapping operator ⟼, yielding the simulation outcomes.( , ) and ( , ) denote the resultant values for and , obtained from the simulation given the inputs and ζ.The expectation [•] considers the stochastic parameters in ζ and is calculated with respect to or .

Figure 6 Figure 6 .
Figure 6 illustrates a 1D optimization problem solved using the GP posterior and acquisition function.The green dots represent the observed data points used for fi ing the GP, resulting in the posterior mean function represented by the red solid line and the posterior variance depicted by the shaded region.The true objective function is illustrated by the black solid line, while the EI acquisition function is represented by the blue solid line.From the figure, it can be observed that the EI acquisition function a ains its peak values in regions where a low posterior mean function value coincides with a large posterior variance.Using iterative steps involving the construction of the GP model, the ) and other dimensions correspond to the uncertainty of those objectives ( ( ), = 1, … , ).In other words, our acquisition function (•) simultaneously minimizes the predictive mean functions and maximizes the posterior variance based on GP models: * = arg min ∈ ( ) = { ( ), − ( )}, ∈ [1, ] with a nonconvex Pareto optimal set: min ( ) = 0.001 ( − 10)( − 60)( − 100) + 1000 min ( ) = 0.001 ( − 70)( − 100)( − 200) + 6000 s. t. ϵ ∩ [0,100]

Figure 7 .
Figure 7. Two objective functions and the Pareto optimal set for Example 1.

Figure 8 .
Figure 8.The rigorous quantification of the quality of a solution in Example 1 using (a) domination measure and (b) linear weighted sum method.

Let be a designated
reference set containing | | evenly distributed points from the true Pareto front (TPF).The function vector (•) = (•), (•), … , (•) comprises the collection of true objective functions.Meanwhile, Z denotes the approximate Pareto front (APF) produced by the evaluated algorithm.The IGD metric quantifies the average proximity between the members of the reference set and the approximate Pareto front, essentially measuring how well the APF mirrors the TPF.Thus, a lower IGD metric indicates superior convergence and diversity.Before discussing the diversity metric Υ in the solution space, it is essential to sequence the | | approximate Pareto optimal solutions , , … , | | produced by the evaluated algorithm, in ascending order based on the (•) criterion, such that ( ) ≤ ( ) ≤ ⋯ ≤ | | .Furthermore, ( ) and (| |) denote the extremities of the APF, termed the leftmost and rightmost boundary points, correspondingly.Similarly, and represent the extremities of the TPF.The diversity metric Υ is then described below: |) ) is the distance between the leftmost (rightmost) boundary point of the APF and the leftmost (rightmost) boundary point of the TPF.= ( ) − ( ) refers to the distance between an approximate Pareto solution and its closest neighbor, whereas ̅ = 1/(| | − 1) × ∑ | |

Figure 10 .
Figure 10.Comparison of the nondominated solutions obtained using various algorithms.

Figure 11 .
Figure 11.Numerical results in term of (a) the iterative curve for IGD, (b) the iterative curve for Υ, and (c) the estimation errors of two GP models for bi-objective functions.

Figure 12 .
Figure 12.Optimization results for PID parameter optimization problem, where (a) represents all explored solutions in the objective space, (b) denotes the Pareto front, and (c) denotes the bi-

Figure 13 .
Figure 13.Parallel coordinate plots of last 50 solutions for NDTW (a) and MALA (b).

Figure 15 .
Figure 15.The results of path tracking achieved by Plan 2. Note: (a,c,e,g) represent the vehicle's tracking trajectory, and (b,d,f,h) represent the lateral acceleration of the vehicle.

Table 1 .
Summary of the PID parameter optimization work.

Table 2 .
Environmental perturbations for path tracking.

Table 3 .
Comparison results for four algorithms.

Table 4 .
Three optimized plans obtained using the DMEGO method.