Machine Learning Approach to Develop a Novel Multi-Objective Optimization Method for Pavement Material Proportion

: Asphalt mixture proportion design is one of the most important steps in asphalt pavement design and application. This study proposes a novel multi-objective particle swarm optimization (MOPSO) algorithm employing the Gaussian process regression (GPR)-based machine learning (ML) method for multi-variable, multi-level optimization problems with multiple constraints. First, the GPR-based ML method is proposed to model the objective and constraint functions without the explicit relationships between variables and objectives. In the optimization step, the metaheuristic algorithm based on adaptive weight multi-objective particle swarm optimization (AWMOPSO) is used to achieve the global optimal solution, which is very efﬁcient for the objectives and constraints without mathematical relationships. The results showed that the optimal GPR model could describe the relationship between variables and objectives well in terms of root-mean-square error (RMSE) and R 2 . After the optimization by the proposed GPR-AWMOPSO algorithm, the comprehensive pavement performances were enhanced in terms of the permanent deformation resistance at high temperature, crack resistance at low temperature as well as moisture stability. Therefore, the proposed GPR-AWMOPSO algorithm is the best option and efﬁcient for maximizing the performances of composite modiﬁed asphalt mixture. The GPR-AWMOPSO algorithm has advantages of less computational time and fewer samples, higher accuracy, etc. over traditional laboratory-based experimental methods, which can serve as guidance for the proportion optimization design of asphalt pavement.


Introduction
Nowadays, it has been realized that ordinary asphalt pavements are usually unable to service the development of traffic, due to various forms of destruction including cracking, rutting etc. [1][2][3][4][5][6][7][8].Many researchers around the world tried to use different methods to enhance the service performances of asphalt pavement, and many additives like polymers, diatomite, fibres, etc. were selected to mix with asphalt materials [9][10][11][12][13].Studies demonstrated that diatomite has been effective in improving the anti-rutting performance and water damage resistance [14].In addition, fibres were proved to be a kind of effective material for the crack resistance of asphalt mixture [15].Thus, asphalt mixture proportion design is a key step in asphalt pavement design and application.
Conventionally, asphalt mixture proportions are optimized by laboratory tests including asphalt content, amount of aggregates, additive parameters, etc. to achieve a designed performance that meets requirements [16,17].When there is only one objective to be optimized, the traditional experimental-based method is generally feasible [18].However, in most cases, proportion optimization of asphalt mixture is a complex problem, which needs to optimize multiple objectives simultaneously, such as performances, cost and environmental pollution, etc.Moreover, the proportion optimization problem with multiple objectives for asphalt mixture is usually investigated based on traditional prescriptive approaches and design of experiments (DOE), such as orthogonality design method, factorial experiments and response surface methodology (RSM), etc. [19].This would result in an exponential increase in the number of laboratory trial specimens.These experimental optimization process would be both time-and resource-intensive.On the other hand, asphalt mixture proportions obtained by traditional prescriptive approaches and experiment-based methods are considered as practical solutions, not optimal solutions [20].
To solve the boundedness of traditional prescriptive approaches and experiment-based methods, many research studies have attempted computational design optimization, which is a mathematical approach for mixture proportions using advanced mathematical techniques and high-performance computing [21].In this computational design optimization process, there are three steps including the problem statement, modelling and optimization algorithm [22].The problem statement involves variables, objectives and constraints.The modelling step focuses on the mathematical relationships between variables and objectives.Nowadays, computational optimization methods based on machine learning (ML) technology are popularly used in a lot of fields of engineering for optimization and prediction [23].In most cases, the mixture proportion optimization problem needs to be forecasted without knowing the explicit relationships between variables and objectives [24].Therefore, the objective functions can be developed based on ML techniques.There are some ML models used for the optimization of mixture proportion, e.g., clustering, artificial neural networks, instance-based learning, regression-based methods, decision trees and support vector machines, etc. [25][26][27].After that, the optimization algorithm was employed to mathematically address the problem and obtain its optimal mixture proportion [28].Computational optimization methods based on metaheuristic algorithms are used for mixture optimization problem with the help of simplicity and efficiency.The metaheuristic algorithms known as "guided search" aim to find the global optimal solution for optimization problems, which is very useful for objectives and constraints without mathematical relationships.Many metaheuristic optimization algorithms were developed for dealing with optimization problems.In the field of mixture proportion design, there are multiple competing objectives, creating a rapidly growing area of optimization study [29].For multi-objective problems of mixture proportion design, metaheuristic optimization algorithms based on evolution, genetic and search algorithms are commonly used optimization tools [30].For instance, the optimization of mixture design needs to maximize physical performances as well as minimize both cost and environmental pollution, in which multiple objectives need to be optimized at once.The metaheuristic algorithms e.g., particle swarm optimization (PSO), evolutionary and genetic, etc. are iterative algorithms, which generate and search the optimal solutions based on the fitness function and iterations.
However, in the real design of asphalt mixture, the traditional experimental-based method is still widely employed, implying that it would require a large number of samples and spend more time on a feasible solution.The above ML models were used to deal with the uncertain relationships between variables and objectives of experiment-based methods.Compared to other ML models (such as support vector machine, neural networks, etc.), the Gaussian process regression (GPR) as a non-linear modelling method, is easy to implement and predict with a probability distribution.In addition, the metaheuristic algorithm based on multi-objective particle swarm optimization (MOPSO) is efficient and tractable, which could provide a sufficiently good solution to an optimization problem.
The main objective of this study is to investigate a metaheuristic optimization algorithm for material proportion design of asphalt mixture employing ML models.After the gradation of aggregates for asphalt mixtures is determined, the mixture proportion including asphalt content and additive parameters usually need to be optimized to balance the competing multiple objectives for the best comprehensive physical performances subjected to several constraints like volumetric properties, pavement performances, etc.Since these models between variables and objectives are not explicit expressions, the ML model based on GPR is chosen to train the system model of these objectives and constraints for further multi-objective problems.The next step is to use the metaheuristic algorithm based on MOPSO for searching the global best solutions to solve the multi-objective optimization (MOO) problem of asphalt mixture proportion with objectives and constraints modeled by ML.Therefore, this study applies both ML methods and metaheuristic algorithms to optimize the multi-objective problem of asphalt mixture proportions.Figure 1 shows the structure of this study, in which the metaheuristic optimization process for asphalt mixture proportion based on an ML model is described in Part I and two cases are presented in Part II.
cient and tractable, which could provide a sufficiently good solution to an optimization problem.
The main objective of this study is to investigate a metaheuristic optimization algorithm for material proportion design of asphalt mixture employing ML models.After the gradation of aggregates for asphalt mixtures is determined, the mixture proportion including asphalt content and additive parameters usually need to be optimized to balance the competing multiple objectives for the best comprehensive physical performances subjected to several constraints like volumetric properties, pavement performances, etc.Since these models between variables and objectives are not explicit expressions, the ML model based on GPR is chosen to train the system model of these objectives and constraints for further multi-objective problems.The next step is to use the metaheuristic algorithm based on MOPSO for searching the global best solutions to solve the multi-objective optimization (MOO) problem of asphalt mixture proportion with objectives and constraints modeled by ML.Therefore, this study applies both ML methods and metaheuristic algorithms to optimize the multi-objective problem of asphalt mixture proportions.Figure 1 shows the structure of this study, in which the metaheuristic optimization process for asphalt mixture proportion based on an ML model is described in Part I and two cases are presented in Part II.

Definition for Multi-Objective Optimization (MOO) Problem
The MOO problem is generally regarded as a mathematical optimization problem with multiple objective functions to be optimized simultaneously [31].Mixture proportion optimization involves determining specifications and contents of raw materials.Generally, some objectives are used to evaluate the mechanical and pavement performances, such as strength, stability, and so on.With the diversification of additives and modifiers, competing objectives, such as additive contents and performances are introduced when designing a proportion for road materials.At present, the mixture proportion optimization methods involving only a single mechanical performance have gradually lost recognition.Many researchers tried to construct multi-objective mathematical models to obtain better mixture proportion schemes.However, as far as we know, there

Definition for Multi-Objective Optimization (MOO) Problem
The MOO problem is generally regarded as a mathematical optimization problem with multiple objective functions to be optimized simultaneously [31].Mixture proportion optimization involves determining specifications and contents of raw materials.Generally, some objectives are used to evaluate the mechanical and pavement performances, such as strength, stability, and so on.With the diversification of additives and modifiers, competing objectives, such as additive contents and performances are introduced when designing a proportion for road materials.At present, the mixture proportion optimization methods involving only a single mechanical performance have gradually lost recognition.Many researchers tried to construct multi-objective mathematical models to obtain better mixture proportion schemes.However, as far as we know, there are few papers reported to investigate a general optimization model of asphalt mixture proportion.Thus, the MOO problem for mixture proportion is a complex combinatorial optimization problem with more challenges.
MOO problem could be defined to search the optimal set of input variables for minimizing and/or maximizing multiple objectives simultaneously.A typical MOO mathematical formulation can be expressed as: where x is the vector of design variables, x ∈ [x 1 , x 2 , . . ., T denotes the objective function vector containing m objective functions; g j (x) and p represent the constraint function and the number of constraints, respectively.
In general, a solution to the MOO problem in Equation ( 1) is known as a Pareto optimal solution without degrading some objective values, which is considered equally efficient.Thus, from different viewpoints, a set of Pareto optimal solutions could be found according to the goals, that could satisfy the subjective preferences.Due to multiple Pareto optimal solutions for MOO problems, these usually require a lot of evaluations.Therefore, some methods are used to convert the original MOO problem into a single-objective optimization problem by a weighted-average method as follows [32]: where ω q is the q-th weight coefficient; f q * is the normalized factor to the q-th objective function.However, the above weighted-average method has some disadvantages.The optimal solution by the method relies on a suitable weight coefficient [33].Moreover, the scaling coefficients are often chosen based on the ideal optimal solutions of each single optimization problem which can be rather time-consuming with increase in the number of objectives.It cannot find certain Pareto-optimal solutions in the case of a nonconvex objective space.For the above reasons, this study aims to investigate a novel MOO method by combining the GPR based on ML methods with PSO based on metaheuristic algorithms.The GPR method is utilized as an optimizer for the objective functions, and the PSO algorithm is adopted to search the optimal solution for MOO problems.

System Modeling Based on Machine Learning (ML) Method-Gaussian Process Regression (GPR)
A Gaussian process (GP) is a stochastic process, which consists of random variables with normal distributions.Thus, the GP distribution is a joint distribution of any number of random variables (X = [x 1 , x 2 , . . ., x n ] T ) [34].From the function-space view, a GP is denoted by a mean function µ(x) as well as covariance function k(x, x'), and as such, the GP could be expressed in Equation (3) as follows: where , in which the mean function is usually zero for notational simplicity.The covariance function is considered as a crucial ingredient that defines nearness or similarity.
In GPR, the squared exponential (SE) covariance function is probably the most commonly-used kernel, that corresponds to a Bayesian linear regression model.The SE covariance function is expressed as follows: where the amplitude σ f and characteristic length-scale l are two hyperparameters to model a smooth rising trend.
For the standard linear regression model with Gaussian noise, the Bayesian analysis is written as follows: where x is the input vector, y is the observed value considering noise pollution, f is the objective function, f (x) is the function value.It can be shown that y differs from f (x) by ε.
Here, suppose that ε follows an independent, uniformly Gaussian distribution (µ = 0, σ 2 n ), it has the following distribution: Subsequently, for multivariate input variables, the covariance function for the standard linear regression model is expressed in Equation ( 7): where n is the dimension of vector x (x s , x s ', in which s = 1, 2, . . ., n), l s is the characteristic length-scale of the s-th input variable.Assume a training sample set S = {(x i , y i ), i = 1, 2, . . ., t)}, t is the number of the training sample.According to the GP described in Equation ( 3), the function values in the form of a vector is F = [ f (x 1 ), f (x 2 ), . . . ,f (x t )], which abides by a joint Gaussian distribution.The prior distribution of y can be described in Equation (8): where I t is the identity vector of size t, K(x, x) donates the t × t order symmetric positive definite covariance vector.Subsequently, given a new test input x * , the test function value f * is generated from the joint posterior distribution by the corresponding mean and covariance vector.The joint distribution between the observed values and predictive values can be written as follows: where K(x, x * ), K(x * , x) and K(x * , x * ) are covariance matrix calculated through Equation (7), and here the expressions: K(x, x * ) = K(x * , x) T .Thus, by deriving the conditional distribution corresponding to Equation (9), the posterior distribution of the test function value f * for GPR can be expressed as [35]: where, the mean value: the variance of f * : the mean value f * is considered as a point estimate, and the variance cov(f * ) is used to estimate the uncertainty.Figure 2 presents the graphical model for GPR.
the mean value  f is considered as a point estimate, and the variance cov( )  f is used to estimate the uncertainty.Figure 2 presents the graphical model for GPR.

Search Information Updating Based on Metaheuristic Algorithm Adaptive Weight Multi-Objective Particle Swarm Optimization (AWMOPSO)
The PSO algorithm proposed by Kennedy and Eberhart is a metaheuristic algorithm in view of swarm intelligence inspired by the movement behaviour in a bird flock or fish cluster [36].PSO algorithms make few assumptions for the optimization problem and can search a large candidate solution space.It has been widely applied in optimization problems because of the simplicity and efficiency for gradient-free optimization algorithms.In PSO algorithm, each particle (bird or fish) searches for food in a cooperative way and learns from the experience of the swarm.Then, global collective behaviour emerges in the swarm to locate the optimal position of food by sharing information among the swarm.PSO algorithms try to improve candidate solutions in the corresponding search space for optimizing a problem by iteration.The PSO process starts with a stochastically initialized particle swarm as possible optimal solutions for a problem and every particle is represented with locations as well as velocity.All position vectors would be assessed by the fitness function, which could evaluate and judge a position [37].
Considering a swarm with p particles in an n-dimensional space, the position and velocity of the i-th particle at t iteration are represented 12 ( , ,..., ) ( , ,..., ) , respectively.These points of velocity and position of the i-th particle at t iteration can be updated and adjusted according to Equations ( 13) and ( 14).There are three components for each particle's movement at each iteration in a PSO algorithm, i.e., velocity update, individual cognition, and social learning. where x  are the velocity and position of the i-th particle at (t + 1) iteration for the n-th dimension.
For the first term of velocity update in Equation ( 13), w denotes the inertia weight parameter determining how much the particles' previous motion is preserved.In the individual cognition term, an individual-cognition parameter c1 is a positive constant; r1 is a random value within [0,1], avoiding premature convergences and increasing the most likely global optima; pbesti represents the i-th particle's own best position.As for the third one (social learning) meaning that all particles in the swarm share the information of the globally best point achieved (gbest), a social learning parameter c2 is also a positive constant; r2 plays exactly the same role as r1.In general, set c1 = c2 = 2 for acceleration, w = wmax − iter  (wmax − wmin)/itermax, the value of w decreases with the increase of iteration times.Besides, the update for a particle's velocity and position at t iteration is presented in Figure 3.

Search Information Updating Based on Metaheuristic Algorithm Adaptive Weight Multi-Objective Particle Swarm Optimization (AWMOPSO)
The PSO algorithm proposed by Kennedy and Eberhart is a metaheuristic algorithm in view of swarm intelligence inspired by the movement behaviour in a bird flock or fish cluster [36].PSO algorithms make few assumptions for the optimization problem and can search a large candidate solution space.It has been widely applied in optimization problems because of the simplicity and efficiency for gradient-free optimization algorithms.In PSO algorithm, each particle (bird or fish) searches for food in a cooperative way and learns from the experience of the swarm.Then, global collective behaviour emerges in the swarm to locate the optimal position of food by sharing information among the swarm.PSO algorithms try to improve candidate solutions in the corresponding search space for optimizing a problem by iteration.The PSO process starts with a stochastically initialized particle swarm as possible optimal solutions for a problem and every particle is represented with locations as well as velocity.All position vectors would be assessed by the fitness function, which could evaluate and judge a position [37].
Considering a swarm with p particles in an n-dimensional space, the position and velocity of the i-th particle at t iteration are represented x t i = (x i1 , x i2 , . . ., x in ) and v t i = (v i1 , v i2 , . . ., v in ), respectively.These points of velocity and position of the i-th particle at t iteration can be updated and adjusted according to Equations ( 13) and ( 14).There are three components for each particle's movement at each iteration in a PSO algorithm, i.e., velocity update, individual cognition, and social learning.
where v t+1 i and x t+1 i are the velocity and position of the i-th particle at (t + 1) iteration for the n-th dimension.
For the first term of velocity update in Equation ( 13), w denotes the inertia weight parameter determining how much the particles' previous motion is preserved.In the individual cognition term, an individual-cognition parameter c 1 is a positive constant; r 1 is a random value within [0,1], avoiding premature convergences and increasing the most likely global optima; pbest i represents the i-th particle's own best position.As for the third one (social learning) meaning that all particles in the swarm share the information of the globally best point achieved (gbest), a social learning parameter c 2 is also a positive constant; r 2 plays exactly the same role as r 1 .In general, set c 1 = c 2 = 2 for acceleration, w = w max − iter × (w max − w min )/iter max , the value of w decreases with the increase of iteration times.Besides, the update for a particle's velocity and position at t iteration is presented in Figure 3.

Figure 3. The update of velocity and position at t iteration regarding a bi-dimensional pro
As mentioned above, PSO could effectively solve single optimization p Based on PSO, MOPSO algorithm was developed for MOO problems.The perf of PSO largely depends on its control parameters, in which the inertia weight pa (w) is an important parameter.A larger value of w is beneficial to improve the alg global search capability, while a smaller value of w will enhance the algorithm search capability.The MOPSO algorithm is a continuation of the PSO algorithm tively solving single-objective problems.It adds non-dominated ideas to the sea ciple of the optimal solution.The improved adaptive weight MOPSO (AWMO gorithm considers the fitness of the overall particle, and establishes a new self-renewal strategy [38].The improved inertia weight is as follows: , , where wmax and wmin are the maximum and minimum inertia weights, respectivel current fitness value of the particle; Javg and Jmin are the average and minimum value, respectively.Similarly, the AWMOPSO seeks the global optimal position ing to the positions of an individual particle and particle swarm.Several globa exist that constitute the Pareto front and the AWMOPSO collects all the non-do particles into a repository.Finally, each particle selects its leader in the reposito on different decisions.The flowchart of the proposed method is illustrated in F which is the integration AWMOPSO with GPR (i.e., GPR-AWMOPSO).As mentioned above, PSO could effectively solve single optimization problems.Based on PSO, MOPSO algorithm was developed for MOO problems.The performance of PSO largely depends on its control parameters, in which the inertia weight parameter (w) is an important parameter.A larger value of w is beneficial to improve the algorithm's global search capability, while a smaller value of w will enhance the algorithm's local search capability.The MOPSO algorithm is a continuation of the PSO algorithm in effectively solving single-objective problems.It adds non-dominated ideas to the search principle of the optimal solution.The improved adaptive weight MOPSO (AWMOPSO) algorithm considers the fitness of the overall particle, and establishes a new weight self-renewal strategy [38].The improved inertia weight is as follows: where w max and w min are the maximum and minimum inertia weights, respectively; J is the current fitness value of the particle; J avg and J min are the average and minimum fitness value, respectively.Similarly, the AWMOPSO seeks the global optimal position according to the positions of an individual particle and particle swarm.Several global optima exist that constitute the Pareto front and the AWMOPSO collects all the non-dominated particles into a repository.Finally, each particle selects its leader in the repository based on different decisions.The flowchart of the proposed method is illustrated in Figure 4, which is the integration AWMOPSO with GPR (i.e., GPR-AWMOPSO).

Choice Criteria for MOO Problems
The technique for order of preference by similarity to ideal solution (TOPSIS) is a multi-criteria decision analysis method, which was originally developed by Ching-Lai Hwang and Yoon [39,40].The TOPSIS method compares and selects a set of alternatives by normalising scores for each criterion and calculating the geometric distance between each alternative and the ideal alternative, which is the best score in each criterion.Calculate the separation measures, using the k-dimensional Euclidean distance.The separations ( i d  and i d  ) of each solution from the positive ideal solution and negative ideal solution are given as Equations ( 16) and (17).The relative closeness coefficient (Ri) can be calculated by Equation (18).

Choice Criteria for MOO Problems
The technique for order of preference by similarity to ideal solution (TOPSIS) is a multi-criteria decision analysis method, which was originally developed by Ching-Lai Hwang and Yoon [39,40].The TOPSIS method compares and selects a set of alternatives by normalising scores for each criterion and calculating the geometric distance between each alternative and the ideal alternative, which is the best score in each criterion.Calculate the separation measures, using the k-dimensional Euclidean distance.The separations (d + i and d − i ) of each solution from the positive ideal solution and negative ideal solution are given as Equations ( 16) and (17).The relative closeness coefficient (R i ) can be calculated by Equation (18).
Appl.Sci.2021, 11, 835 where F ij is the normalized value; F + j and F − j are the ideal solution and nadir solution, respectively.By ranking the preference order, the solution with the highest value of closeness coefficient (R i ) is the best solution.

Computational Examples and Analysis
To verify the proposed GPR-AWMOPSO optimization method, two tri-objective optimization cases for modified asphalt mixtures with different gradations are discussed.The optimization results based on the traditional experimental design methods (i.e., RSM and orthogonal experimental design) and the proposed GPR-AWMOPSO optimization algorithm are compared.

Case Study I: Bi-Objective Optimization
A lot of research has shown that basalt fibre has a definite improvement influence on the crack resistance of asphalt materials [41].For basalt fibre-reinforced asphalt mixture, the detailed preparation parameters are discussed and analyzed including fibre length and its content as well as asphalt-aggregate ratio.At the same time, three preparation parameters are optimized based on the volumetric and strength properties.

Data Description
In the optimization process of SMA-13 containing basalt fibre designed by RSM, raw materials included SBS-modified asphalt, coarse and fine aggregates, mineral filler and basalt fibre, whose basic physical properties were summarized in the previous study [16].The experiment with three independent factors at three experimental levels was devised according to the face-centered central composite RSM.The preparation parameters are regarded as fibre content (X 1 ) and its length (X 2 ) as well as asphalt-aggregate ratio (X 3 ); and Marshall stability (MS) (Y 1 ) and flow (FV) (Y 2 ) are performance parameters, air voids (VA) (Y 3 ), voids in mineral aggregates (VMA) (Y 4 ) and voids filled with asphalt (VFA) (Y 5 ) represent volumetric properties.The experimental design and outputs with 19 groups of basalt fibre-reinforced asphalt mixtures are listed in detail in Table 1.In this work, a typical proportion optimization for asphalt mixture is discussed.A metaheuristic optimization model for asphalt mixture proportion based on an ML model is presented to maximize physical performances and satisfy the requirements of volumetric properties.Hereto, the mixture proportion problem is usually taken to be a complex optimization problem.
(1) Objective Function Following the Chinese specification JTG F40-2004, the targets of asphalt mixture proportion optimization should not only meet the volumetric properties but also maximize the strength performance and minimize cost.According to the modelled GPR of basalt fibre-reinforced asphalt mixture, the objective function is calculated as follows: where Y 1 represents the Marshall stability (MS); F 2 (X) is the cost, the prices per kilogram of basalt fibre and asphalt are 12 and 2.2, respectively.
(2) Constraints Following the Chinese specification JTG F40-2004, the range constraints and volumetric properties constraints for SMA-13 specimens reinforced with basalt fibre are summarized in Tables 2 and 3, respectively.The corresponding constraint expressions are written as: (1) System Models Based on GPR As mentioned above, GPR is adopted to model objective functions based on ML methods, which would be trained using the collected data.There are three inputs and five outputs, the whole dataset includes several groups selected randomly from 19 orthogonal groups and is also separated into two datasets for training along with testing.Then we can predict the volumetric and physical performances of asphalt mixture containing basalt fibre from 25% of the selected dataset by using the training dataset from 75% of the selected dataset randomly.Besides, the root-mean-square error (RMSE) in Equation ( 21) is selected to assess the difference between observed values and predictive values and the squared correlation coefficient (R 2 ) in Equation ( 22) is selected to evaluate the correlation between these two values [42].
where y * i and y i denote the predictive and observed values, y * and y are the mean values of the observed values and predictive values.Figure 5 shows the reliable results of GPR model versus size of training dataset.It can be observed that the reliability of the GPR model is better with the size of training dataset increasing.Therefore, in the following study, the training dataset consisting of 8 groups would be selected randomly from the experimental groups for the GPR model, which possesses a good prediction effect.
thogonal groups and is also separated into two datasets for training along with testing.Then we can predict the volumetric and physical performances of asphalt mixture containing basalt fibre from 25% of the selected dataset by using the training dataset from 75% of the selected dataset randomly.Besides, the root-mean-square error (RMSE) in Equation ( 21) is selected to assess the difference between observed values and predictive values and the squared correlation coefficient (R 2 ) in Equation ( 22) is selected to evaluate the correlation between these two values [42].The prediction results for the volumetric and physical performances of asphalt mixture containing basalt fibre are listed in Table 4. Figure 6 shows the predictive values by the GPR model versus the experimental values for training and testing sets.The regression analysis indicates that the optimum GPR model was successful in learning the relationship between the input and output variables.There are lower RMSE values and higher R 2 values for both training and testing datasets.It can be indicated that the GPR model can accurately model the relationship between input and output variables.Thus, it could be confirmed that the optimal GPR model has been well trained based on the comparison between the volumetric and physical properties of training and testing datasets.The prediction results for the volumetric and physical performances of asphalt mixture containing basalt fibre are listed in Table 4. Figure 6 shows the predictive values by the GPR model versus the experimental values for training and testing sets.The regression analysis indicates that the optimum GPR model was successful in learning the relationship between the input and output variables.There are lower RMSE values and higher R 2 values for both training and testing datasets.It can be indicated that the GPR model can accurately model the relationship between input and output variables.Thus, it could be confirmed that the optimal GPR model has been well trained based on the comparison between the volumetric and physical properties of training and testing datasets.(2) Proportion Optimization Results for Basalt Fibre-Modified Asphalt Mixture using GPR-AWMOPSO According to the above AWMOPSO algorithm, the optimal solution can be obtained for the proportion of asphalt mixture containing basalt fibre.The bi-objective modified asphalt mixture optimization results are shown in Figure 7a.GPR is used as the objective function for modelling proportion design of asphalt mixture.The Pareto front based on mechanical performance (Marshall stability) and cost is obtained.These points are distributed within a reasonable range, showing a good effectiveness of GPR-AWMOPSO.Meanwhile, it is observed that the asphalt mixture cost would increase with the increasing of mechanical performance.In total, several non-dominated solutions (the optimal mixture proportions) are selected.In Figure 7a, the solution with higher mechanical performance on the Pareto front is more expensive.In general, the optimal solution depends on the mixture design consideration.According to the closeness coefficient of TOPSIS on the Pareto front, the point A has the highest TOPSIS score.The solution B on the Pareto front is more expensive but has a higher mechanical performance, whereas the solution C has a lower mechanical performance with a correspondingly lower cost.It also shows that the solution A by TOPSIS has an intermediate cost and Marshall stability compared with the solutions B and C. With the mechanical performance as the primary goal, the point B marked by a red five-pointed star will be the best option to achieve the maximum mechanical performance.While minimizing cost is the main objective, the point C marked by a green five-pointed star will be the best option.
To verify the proposed GPR-AWMOPSO algorithm, the optimized proportion of asphalt mixture is compared with the results of RSM, traditional MOPSO [43][44][45] and meshadaptive direct search (MADS) [46][47][48].When the objective is to achieve the maximum Marshall stability, the three inputs and five outputs variables are summarized in Table 5.
From Table 5, it can be found that the optimal results yield higher mechanical properties without a noticeable increase in the mixture proportions.Meanwhile, these results are well satisfied with the constraints shown in Table 3 following the Chinese specification JTG F40-2004.[16], traditional multi-objective particle swarm optimization (MOPSO), mesh-adaptive direct search (MADS) and GPR-adaptive weight multi-objective particle swarm optimization (AWMOPSO) algorithm.Moreover, Figure 7b shows the swarm maximum value of the objective function (at the swarm best position) versus iteration for traditional MOPSO, MADS and GPR-AWMOPSO.After several numbers of iteration, it can be found that a significant increase in the maximum objective function value occurs, indicating that both traditional MOPSO, MADS and GPR-AWMOPSO algorithms are efficient in the proportion optimization of composite modified asphalt mixture.The traditional MOPSO, MADS and GPR-AWMOPSO algorithms also possess good precision.As for calculation speed, the proposed GPR-AWMOPSO in this paper starts to converge rapidly after about 6 iterations, while the convergence rate of the traditional MOPSO is slightly slower.In contrast, the MADS method also has better results.The average iteration number of the MADS method is slightly larger than that of GPR-AWMOPSO, while the convergence rate of GPR-AWMOPSO in this paper is stable and fast.The optimal solution generally relies on the considerations of mixture design.Compared with the RSM design and traditional MOPSO, the proposed GPR-AWMOPSO algorithm has higher mechanical performances.When the main objective is maximizing mechanical performances, the proposed GPR-AWMOPSO algorithm can be considered as the best option.le 5. Comparison of optimization results among RSM [16], traditional multi-objective particle swarm optimization PSO), mesh-adaptive direct search (MADS) and GPR-adaptive weight multi-objective particle swarm optimization MOPSO) algorithm.(

3) Performance Verification
To assess the comprehensive pavement performances of the asphalt mixture containing basalt fibre designed using RSM and the proposed GPR-AWMOPSO algorithm, the permanent deformation resistance, crack resistance as well as moisture stability were conducted at high and low temperatures.For the two groups of designed mixture proportions, three replicate specimens were prepared for the rutting experiment at 60 °C, (3) Performance Verification To assess the comprehensive pavement performances of the asphalt mixture containing basalt fibre designed using RSM and the proposed GPR-AWMOPSO algorithm, the permanent deformation resistance, crack resistance as well as moisture stability were conducted at high and low temperatures.For the two groups of designed mixture proportions, three replicate specimens were prepared for the rutting experiment at 60 • C, splitting experiment at −10 • C as well as an immersion Marshall experiment and freeze-thaw splitting experiment, and the corresponding results are illustrated in Figure 8.It can be observed that compared with the specimens designed by RSM, the dynamic stability of specimens designed by GPR-AWMOPSO increases by 3.59%, the stiffness modulus increases by 2.43%, and the residual Marshall stability and tensile strength ratio (TSR) increase by 3.48% and 2.16%, respectively.The comparison results verify that the proposed GPR-AWMOPSO algorithm is the best option for maximizing the performances of basalt fibre reinforced asphalt mixture.
thaw splitting experiment, and the corresponding results are illustrated in Figure 8.It can be observed that compared with the specimens designed by RSM, the dynamic stability of specimens designed by GPR-AWMOPSO increases by 3.59%, the stiffness modulus increases by 2.43%, and the residual Marshall stability and tensile strength ratio (TSR) increase by 3.48% and 2.16%, respectively.The comparison results verify that the proposed GPR-AWMOPSO algorithm is the best option for maximizing the performances of basalt fibre reinforced asphalt mixture.

Case Study II: Tri-Objective Optimization
Diatomite has been proved to dramatically improve the high-temperature performance of asphalt mixture [49][50][51][52].Hence, the incorporation of diatomite and basalt fibres would comprehensively enhance the resistance to high-temperature rutting, low-temperature cracking as well as moisture stability of asphalt mixture.There are four basic components of asphalt mixture, i.e., asphalt, mineral, diatomite, basalt fibre.With the selected gradation, it needs to take into account both volumetric properties and pavement performances at high and low temperatures when the proportions of asphalt mixture with diatomite and basalt fibre are studied.

Data Description
Raw materials for asphalt mixture with diatomite and basalt fibre included A-90# asphalt, coarse and fine aggregates, mineral filler, diatomite, basalt fibre.The corresponding technical properties as well as origins have been illustrated in the previous studies, which could satisfy the requirements of the Chinese specification JTG F40-2004.In this study, Marshall cylinder specimens were prepared following JTG E20-2011, in which the median gradation of AC-13 was selected.
In the previous study [53], the orthogonal experiment L16 (4 3 ) was designed as a three-factor and four-level experiment.The diatomite content (X1), basalt fibre content (X2), and asphalt-aggregate ratio (X3) were regarded as orthogonal input variables; VA (Y1), VFA (Y2), VMA (Y3), MS (Y4) and splitting strength at −10 °C (Sb, Y5) were used as output variables.These 16 designed groups of asphalt mixture specimens were tested and the orthogonal experimental results including volumetric properties as well as high-temperature and low-temperature performances are listed in Table 6.

Case Study II: Tri-Objective Optimization
Diatomite has been proved to dramatically improve the high-temperature performance of asphalt mixture [49][50][51][52].Hence, the incorporation of diatomite and basalt fibres would comprehensively enhance the resistance to high-temperature rutting, low-temperature cracking as well as moisture stability of asphalt mixture.There are four basic components of asphalt mixture, i.e., asphalt, mineral, diatomite, basalt fibre.With the selected gradation, it needs to take into account both volumetric properties and pavement performances at high and low temperatures when the proportions of asphalt mixture with diatomite and basalt fibre are studied.

Data Description
Raw materials for asphalt mixture with diatomite and basalt fibre included A-90# asphalt, coarse and fine aggregates, mineral filler, diatomite, basalt fibre.The corresponding technical properties as well as origins have been illustrated in the previous studies, which could satisfy the requirements of the Chinese specification JTG F40-2004.In this study, Marshall cylinder specimens were prepared following JTG E20-2011, in which the median gradation of AC-13 was selected.
In the previous study [53], the orthogonal experiment L 16 (4 3 ) was designed as a threefactor and four-level experiment.The diatomite content (X 1 ), basalt fibre content (X 2 ), and asphalt-aggregate ratio (X 3 ) were regarded as orthogonal input variables; VA (Y 1 ), VFA (Y 2 ), VMA (Y 3 ), MS (Y 4 ) and splitting strength at −10 • C (S b , Y 5 ) were used as output variables.These 16 designed groups of asphalt mixture specimens were tested and the orthogonal experimental results including volumetric properties as well as high-temperature and low-temperature performances are listed in Table 6.(1) Objective Function Considering the extremely variable climate in the seasonal frozen area in northeast China, high-and low-temperature performances are treated as optimization targets.The objective functions for asphalt mixture containing diatomite and basalt fibre are modelled using GPR based on ML methods.The objective functions of high-and low-temperature performances as well as costs are calculated by the following equations: where Y 4 (MS) represents the high-temperature performance and Y 5 (S b ) represents the low-temperature performance; F 3 (X) is the cost, the prices per kilogram of diatomite, basalt fibre and asphalt are 1.35, 12 and 2.2, respectively.
(2) Constraints Similarly, following JTG F40-2004 and the previous study [53], the range constraints and volumetric properties constraints for AC-13 specimens modified by diatomite and basalt fibre are summarized in Tables 7 and 8, respectively.The corresponding constraint expressions are written as: (1) System Models Based on GPR For the GPR system model, 12 groups are selected randomly from 16 orthogonal groups and also divided into training and testing dataset.Then we can predict the volumetric and physical properties of diatomite and basalt fibre composite-modified asphalt mixture from 25% of the selected dataset by using the training dataset from 75% of the selected dataset randomly.
The RMSE and R 2 can be calculated using Equations ( 21) and ( 22), and the prediction results for the volumetric and physical properties of asphalt mixture containing diatomite and basalt fibre are shown in Table 9. Figure 9 shows the predictive values by the GPR model versus the experimental values for training and testing sets.The regression analysis indicates that the optimal GPR model was successful in learning the relationship between the input and output variables.There are lower RMSE values and higher R 2 values for both training and testing datasets.It can be indicated that the GPR model can accurately model the relationship between input and output variables.Thus, it could be confirmed that the optimal GPR model has been well trained based on the comparison between the volumetric and physical properties of training and testing datasets.(2) Proportion Optimization Results for Modified Asphalt Mixture Based on GPR-AWMOPSO According to the above AWMOPSO algorithm, the optimal solution can be obtained for the proportion of asphalt mixture containing diatomite and basalt fibre.The tri-objective modified asphalt mixture optimization results are shown in Figure 10a.GPR (2) Proportion Optimization Results for Modified Asphalt Mixture Based on GPR-AWMOPSO According to the above AWMOPSO algorithm, the optimal solution can be obtained for the proportion of asphalt mixture containing diatomite and basalt fibre.The tri-objective modified asphalt mixture optimization results are shown in Figure 10a.GPR is also used as the objective function for modelling proportion design of asphalt mixture, The Pareto front based on mechanical performances (i.e., Marshall stability and splitting strength at −10 • C) and cost could be obtained.It is observed that these points are distributed within a reasonable range, showing a good effectiveness of GPR-AWMOPSO.Meanwhile, it is seen that the modified asphalt mixture cost would increase with increasing mechanical performance.In total, several non-dominated solutions (the optimal mixture proportions) are selected.In Figure 10a, the solution with higher mechanical performance on the Pareto front is more expensive.In general, the optimal solution depends on the mixture design consideration.Solution B on the Pareto front is more expensive but has higher mechanical performances, whereas solution C has lower mechanical performances with a correspondingly lower cost.(3) Performance Verification The high-temperature permanent deformation resistance, low-temperature crack resistance as well as moisture stability were also conducted to evaluate the comprehensive pavement performances of composite modified asphalt mixture designed by the orthogonal method and the proposed GPR-AWMOPSO algorithm.For the two groups of designed mixture proportions, three replicate specimens were prepared for the rutting experiment at 60 °C, splitting experiment at −10 °C as well as immersion Marshall and freeze-thaw splitting experiments, and the corresponding results are shown in Figure 11.It can be observed that compared with the specimens designed by the orthogonal method, the dynamic stability of specimens designed by GPR-AWMOPSO increases by 6.03%, the tensile strain increases by 9.61%, and the residual Marshall stability and tensile strength ratio (TSR) increase by 1.02% and 1.18%, respectively.The comparison results To verify the proposed GPR-AWMOPSO algorithm, the optimized proportion of asphalt mixture is compared with the results of the orthogonal design, traditional MOPSO and MADS.The three inputs and five outputs variables are summarized in Table 10.From Table 10, it can be found that the optimal results yield higher high-temperature and low-temperature physical performances without a noticeable increase in the mixture proportions.Meanwhile, these results are satisfied well with the constraints following JTG F40-2004.The optimal solution generally relies on the considerations of mixture design.Compared with the traditional orthogonal design, traditional MOPSO and MADS, the proposed GPR-AWMOPSO algorithm has higher physical performances.When maximizing performances is the main objective, the proposed GPR-AWMOPSO algorithm is the best option.Figure 10b displays the swarm maximum value of the objective function (at the swarm best position) versus iteration for traditional MOPSO, MADS and GPR-AWMOPSO.After several numbers of iteration, it can be found that a significant increase in the maximum objective function value occurs, indicating that both traditional MOPSO, MADS and GPR-AWMOPSO algorithms are efficient and possess good precision in the proportion optimization of composite modified asphalt mixture.In terms of calculation speed, the convergence rates of GPR-AWMOPSO and MADS are similar.Moreover, it is observed that the optimization rate based on the proposed GPR-AWMOPSO has a faster convergence rate compared with the traditional MOPSO.
(3) Performance Verification The high-temperature permanent deformation resistance, low-temperature crack resistance as well as moisture stability were also conducted to evaluate the comprehensive pavement performances of composite modified asphalt mixture designed by the orthogonal method and the proposed GPR-AWMOPSO algorithm.For the two groups of designed mixture proportions, three replicate specimens were prepared for the rutting experiment at 60 • C, splitting experiment at −10 • C as well as immersion Marshall and freeze-thaw splitting experiments, and the corresponding results are shown in Figure 11.It can be observed that compared with the specimens designed by the orthogonal method, the dynamic stability of specimens designed by GPR-AWMOPSO increases by 6.03%, the tensile strain increases by 9.61%, and the residual Marshall stability and tensile strength ratio (TSR) increase by 1.02% and 1.18%, respectively.The comparison results verify that the proposed GPR-AWMOPSO algorithm is the best option for maximizing the performances of composite modified asphalt mixture.

Case Study III: Optimization with Four Objectives
The existing literature shows that rubber aggregate and basalt fibre have been proved to help modify concrete [54].Waste rubber materials can be recycled to reduce environmental pollution.On the other hand, it can also improve the performance of concrete.

Case Study III: Optimization with Four Objectives
The existing literature shows that rubber aggregate and basalt fibre have been proved to help modify concrete [54].Waste rubber materials can be recycled to reduce environmental pollution.On the other hand, it can also improve the performance of concrete.

Data Description
The rubber and basalt fibre modified concrete was also optimized by using RSM.Based on the Box Behnken design, the experiment with three independent factors at three experimental levels was designed.The preparation parameters are regarded as the waterbinder ratio (X 1 ) and basalt fibre content (X 2 ) as well as rubber content (X 3 ); and slump (Y 1 ), flexural strength (Y 2 ) and compressive strength (Y 3 ) represent response variables.The experimental design and outputs with 17 groups of rubber and basalt fibre-modified concretes are listed in detail in the study [54].

MOO Problem Formulation for Rubber-Basalt Fibre Composite Modified Concrete
(1) Objective Function According to the requirements of rubber and basalt fibre-modified concrete, the objective functions are calculated as follows: (2) Constraints Following the previous study [54], the range constraints for concretes modified by rubber and basalt fibre are summarized and written as:

Results of Concrete Mix Proportion Optimization using GPR-AWMOPSO
(1) System Models Based on GPR For the GPR system model, 12 groups are selected randomly from 17 groups of RSM experimental design and also divided into training and testing dataset.Then we can predict the performance properties of rubber and basalt fibre composite modified concrete from 25% of the selected dataset by using the training dataset from 75% of the selected dataset randomly.
Figure 12 shows the predictive values by the GPR model versus the experimental values for training and testing sets.The regression analysis indicates that the optimal GPR model with lower RMSE values and higher R 2 values was successful in learning the relationship between the input and output variables.It can be indicated that the GPR model can accurately model the relationship between input and output variables, which could be used for the performance prediction of rubber and basalt fibre modified concrete.(2) Proportion Optimization Results Based on GPR-AWMOPSO According to the above AWMOPSO algorithm, the optimal solution can be obtained for the proportion of concretes modified by rubber and basalt fibre.The modified concrete optimization results are shown in Figure 13.GPR is used as the objective function for modelling proportion design of asphalt mixture.The Pareto front based on response (2) Proportion Optimization Results Based on GPR-AWMOPSO According to the above AWMOPSO algorithm, the optimal solution can be obtained for the proportion of concretes modified by rubber and basalt fibre.The modified concrete optimization results are shown in Figure 13.GPR is used as the objective function for modelling proportion design of asphalt mixture.The Pareto front based on response variables and cost is obtained.These points are distributed within a reasonable range, showing a good effectiveness of GPR-AWMOPSO.In total, several non-dominated solutions (the optimal mixture proportions) are selected.In Figure 13, the solution with higher mechanical performance on the Pareto front is more expensive.In general, the optimal solution depends on the mixture design consideration.According to the closeness coefficient of TOPSIS on the Pareto front, the point A has the highest TOPSIS score.Solution B on the Pareto front is more expensive and has higher compression strength and flexural strength but lower slump.However, solution C has lower compression strength and flexural strength with a correspondingly lower cost.It also shows that solution A by TOPSIS has an intermediate cost and response variables compared with solutions B and C. With the mechanical performance including compression strength and flexural strength as the primary goal, point B will be the best option to achieve the maximum mechanical performance.While when minimizing cost is the main objective, point C will be the best option.Moreover, Table 11 also summarizes the comparison results of three inputs and three outputs variables.From Table 10, compared with the results of the orthogonal de-sign, traditional MOPSO and MADS, it can be found that the optimal results yield higher mechanical performance including compression strength and flexural strength.

Conclusions
In this study, a novel AWMOPSO algorithm based on GPR is proposed for the proportion design and optimization of asphalt mixture, which could be adopted to deal with multi-variable, multi-level problems with multiple constraints.The objective function based on GPR was established and evaluated in terms of RMSE and R 2 .The results showed that the optimal GPR can accurately model the relationship between input and output variables.The well trained GPR model was then applied to the AWMOPSO design of asphalt mixture (SMA-13) containing basalt fibre as well as asphalt mixture (AC-13) containing diatomite and basalt fibre.The proportions of modified asphalt mixtures were considered as input variables and the volumetric and physical properties were output variables.After the optimization, the comprehensive pavement performances were enhanced in terms of the high-temperature permanent deformation resistance, low-temperature crack resistance as well as moisture stability.Therefore, the proposed GPR-AWMOPSO algorithm is regarded as the best option and efficient for maximizing the performances of modified asphalt mixture.
The main objective of this study is the verification of the GPR-AWMOPSO algorithm

Figure 1 .
Figure 1.The study structure of metaheuristic optimization based on the machine learning (ML) model.

Figure 1 .
Figure 1.The study structure of metaheuristic optimization based on the machine learning (ML) model.

Figure 3 .
Figure 3.The update of velocity and position at t iteration regarding a bi-dimensional problem.

Figure 4 .
Figure 4.The flowchart of the procedure of GPR-adaptive weight multi-objective particle swarm optimization (AWMOPSO): (a) Experiments; (b) System model based on GPR; and (c) Search information based on MOPSO.

Figure 4 .
Figure 4.The flowchart of the procedure of GPR-adaptive weight multi-objective particle swarm optimization (AWMOPSO): (a) Experiments; (b) System model based on GPR; and (c) Search information based on MOPSO.
i y denote the predictive and observed values, y  and y are the mean values of the observed values and predictive values.Figure 5 shows the reliable results of GPR model versus size of training dataset.It can be observed that the reliability of the GPR model is better with the size of training dataset increasing.Therefore, in the following study, the training dataset consisting of 8 groups would be selected randomly from the experimental groups for the GPR model, which possesses a good prediction effect.

Figure 5 .
Figure 5.The reliable results of the GPR model versus size of training dataset.

Figure 5 .
Figure 5.The reliable results of the GPR model versus size of training dataset.

Figure 7 .
Figure 7.The proportion optimization results for basalt fibre modified asphalt mixture using GPR-AWMOPSO: (a) Pareto front showing two objectives (F1 and F2) and (b) the primary objective function value (F1) versus iteration.

Figure 7 .
Figure 7.The proportion optimization results for basalt fibre modified asphalt mixture using GPR-AWMOPSO: (a) Pareto front showing two objectives (F 1 and F 2 ) and (b) the primary objective function value (F 1 ) versus iteration.

Figure 8 .
Figure 8. Performance verification between the RSM and GPR-AWMOPSO algorithm: (a) high-and low-temperature performances; and (b) moisture stability.

Figure 8 .
Figure 8. Performance verification between the RSM and GPR-AWMOPSO algorithm: (a) high-and low-temperature performances; and (b) moisture stability.

Figure 10 .
Figure 10.The proportion optimization results for diatomite and basalt fibre modified asphalt mixture using GPR-AWMOPSO: (a) Pareto front showing three objectives (F 1 , F 2 and F 3 ) and (b) the primary objective function value (F 1 + F 2 ) versus iteration.

Figure 11 .
Figure 11.Performance verification between the orthogonal experimental design and GPR-AWMOPSO algorithm: (a) high-and low-temperature performances; and (b) moisture stability.

Figure 11 .
Figure 11.Performance verification between the orthogonal experimental design and GPR-AWMOPSO algorithm: (a) highand low-temperature performances; and (b) moisture stability.

Figure 12 Figure 12 .
Figure12shows the predictive values by the GPR model versus the experimental values for training and testing sets.The regression analysis indicates that the optimal GPR model with lower RMSE values and higher R 2 values was successful in learning the relationship between the input and output variables.It can be indicated that the GPR model can accurately model the relationship between input and output variables, which could be used for the performance prediction of rubber and basalt fibre modified concrete.

Figure 12 .
Figure 12.Comparison of the predictive values by GPR versus the experimental values: (a) Y 1 ; (b) Y 2 ; and (c) Y 3 .

Figure 13 .
Figure 13.The proportion optimization results for rubber and basalt fibre modified concrete using GPR-AWMOPSO: (a) Pareto front showing three objectives (F 1 , F 2 , and F 3 ); (b) Pareto front showing three objectives (F 1 , F 2 , and F 4 ) and (c) Pareto front showing three objectives (F 2 , F 3 and F 4 ).

Table 2 .
Range constrains of input variables.

Table 3 .
Volumetric and physical property constrains of output variables.

Table 4 .
Root-mean-square error (RMSE) and R 2 values for volumetric and physical properties of basalt fibre reinforced asphalt mixture.

Table 7 .
Range constrains of input variables.

Table 8 .
Volumetric and physical property constrains of the output variables.

Table 9 .
RMSE and R 2 for volumetric and physical properties of asphalt mixture.